MITgcm_contrib/jmc_script/grph_StD.m

 prefix='dynStD';
%prefix='iceStD';
 namA={'g02','b02'};
%namA={'r17a'};
 Nexp=size(namA,2);
%-

% $Header: /u/gcmpack/MITgcm_contrib/jmc_script/grph_StD.m,v 1.8 2015/09/06 18:15:09 jmc Exp $
% $Name:  $

nItMx=1e10*ones(1,Nexp); %nItMx(3)=11;
%nItMx=2400*ones(1,Nexp);
namLg=namA ; namLg=strrep(namLg,'_','\_');
undef=123456.7;
%-----------
%- test if the variable krd is define :
if size(who('krd'),1) > 0,
 fprintf('krd is defined and = %i \n',krd);
else
 fprintf('krd undefined ; set to 1 \n'); krd=1 ;
end
if krd > 0,
%- define list of fields to read in:
%listV={'Eta','U','V','W','T','S','DETADT2','RELHUM','Phi'};
%listV={'Eta','W','T','S','CONVADJ','DETADT2'};
%listV={'Eta','UE_VEL_C','VN_VEL_C','W','T','DETADT2','Phi'};
%- or take all them:
 clear listV ; listV='all_flds';
%-----------

%- start to read the longest record:
  n=1; rf=-1; if strcmp(char(listV),'all_flds'), rf=0; end
  [ntA(n),rList,tim,vv1,listV] = ...
    read_StD(prefix,namA(n),listV);
  nIt=ntA(n); nk=size(vv1,1); nRg=size(vv1,3);
%- set global dims: & load vvA --> vvB
  nbV=size(listV,2);
  nrec=nIt; n3d=nk; nReg=nRg;
  vvA=zeros(n3d,nrec,nReg,5,nbV,Nexp); tiA=zeros(nrec,2,Nexp);
  vvA(1:nk,1:nIt,1:nRg,:,:,n)=vv1; tiA(1:nIt,:,n)=tim;
%----
 for n=2:Nexp,
  [ntA(n),rList,tim,vv1,listV] = ...
    read_StD(prefix,namA(n),listV);
  nIt=ntA(n); nk=size(vv1,1); nRg=size(vv1,3);
  if (nrec < nIt),
    fprintf('\n');
    error([' Nb of records=',int2str(nIt),' exceeds nrec=',int2str(nrec)]);
  end
  if (n3d < nk),
    fprintf('\n');
    error([' Nb of Levels=',int2str(nk),' exceeds n3d=',int2str(n3d)]);
  end
  vvA(1:nk,1:nIt,1:nRg,:,:,n)=vv1; tiA(1:nIt,:,n)=tim;
 end;
 if krd == 2,
   fprintf('save to "sav_StD.mat" file ...');
   save('sav_StD.mat','vvA','tiA','ntA','rList','listV');
   fprintf(' done\n')
 end
elseif krd < 0,
 fprintf('load from "sav_StD.mat" file ...');
 load sav_StD
 fprintf(' done\n'); nbV=size(listV,2);
end
if krd ~= 0,
 ttA=squeeze(tiA(:,2,:));
 ttA=ttA/3600; titT='hrs';  ttA=ttA/24; titT='days';
 ttA=ttA/30 ; titT='month'; ttA=ttA/12 ; titT='year';
end
%=========================================================

ttax1=0 ; ttax2=0 ; ttay=zeros(nbV,2);
%-- fixed time axis bound :
% ttax1=15.; ttax2=20.;
%-- fixed Y axis bound :
% ttay(4,:)=[0 0.6];
%-----------
fprintf('Total length: ntA=');fprintf(' %i ,',ntA); fprintf(' \n');
for n=1:Nexp,
fprintf(' exp %i : time(d):%10.2f ->%10.2f \n', n,ttA(1,n),ttA(ntA(n),n) );
end;
%--

list_on=zeros(1,nbV);
nbG=18;
nbG=min(nbG,nbV); list_on(1:nbG)=1 ;
%if nbG < nbV, list_on=0; list_on(nbG:nbV)=1; end %- to get 2nd set of plots
%list_on(1:6)=[1 1 1 1 1 1];

isA=ones(1,Nexp); ieA=ntA;
%- limit the length : for search of isA <->1500y: find(ttA(:,2) == 1500)
%isA=isA*31 ; % drop the 1rst mnth (1 Monitor/d)
isA=isA*2 ; % drop the 1rst mnth (1 Monitor/30d)
%isA(1)=31 ; isA(2)=4 ; % drop the 1rst mnth
%ieA(:)=240; %isA(:)=1;

linA(1,:)='k-'; % ieA(1)=60 ; % ieA(1)=1152 ;
linA(2,:)='b-';
linA(3,:)='r-';
linA(4,:)='g-';
linA(5,:)='m-';
linA(6,:)='c-';

ieA=min(ieA,nItMx);
%titall='AIM , Cubic-G (32x32) , cpl-FM Forcing' ;
titall='Global Ocean, Cubic-G (32x32) , CORE Forc (2)' ;
%titall='Dyncore test-case 5 (cs-32)' ;

%=========================================================
ng=0; fxb=100; fyb=60; fdx=100; fdy=40; fsc=1.;
%fyb=-360; fxb=-2600; fdy=60; fsc=1.5;

for jv=1:nbV,
%-------------------
 flag=list_on(jv); kl=0;
 vv1=vvA(:,:,:,:,jv,:); namV=char(listV(jv));
 titv=strrep(namV,'_','\_');
%if strcmp(namV,'Eta'), vv1=vv1/100; titv='Eta [mb]'; end
%if strcmp(namV,'T'), kl=1; end         % <-- to get surf.Temp
 if kl > 0, titv=[titv,'\_',int2str(kl)];
   fprintf([' var= ',namV,' at level k= %i \n'],kl);
 end

 if flag == 1
%--
  ng=ng+1;
%- reset "ng" to jv for fix fig number (independent of which one is in list_on):
  %ng=jv;
  figure(ng); set(ng,'position',[fxb+fdx*ng fyb+fdy*ng [500 700]*fsc]);clf;
  var=squeeze(vv1(1+kl,:,1,:,:)); dd=zeros(5,Nexp); av=zeros(5,Nexp);
  for n=1:Nexp,
   tmp=var(isA(n):ieA(n),:,n); [I]=find(tmp(:,1)==undef); tmp(I,:)=0;
   av(:,n)=sum(tmp);
   if length(I) < size(tmp,1), av(:,n)=av(:,n)/(size(tmp,1)-length(I)); end
   tmp(I,:)=NaN; dd(:,n)=max(tmp)-min(tmp);
  end
  var(find(var==undef))=NaN;
  for nv=1:4,
    subplot(410+nv); ttmn=' Mx-mn:'; ttav=' Av:';
    for n=1:Nexp,
      plot(ttA(isA(n):ieA(n),n),var(isA(n):ieA(n),nv,n),linA(n,:));
      if n == 1, hold on ; end ;
      ttmn=sprintf([ttmn,' %2.1e ;'],dd(nv,n));
      ttav=sprintf([ttav,' %3.2e ;'],av(nv,n));
    end ; hold off ;
    AA=axis ; dAA=AA(4)-AA(3);
    if AA(3)*AA(4) <= 0, AA(3)=min(AA(3),-dAA/10); AA(4)=max(AA(4),dAA/10); end
    if ttax1 < ttax2, AA(1)=ttax1; AA(2)=ttax2; end;
    axis(AA); grid ;
    if nv == 1, title(['Avr ',titv,'  ',ttmn]); end
    if nv == 2, title(['Std-Dev ',titv,'  ',ttav]); end
    if nv == 3, title(['min ',titv,'  ',ttav]);
                legend(namLg,'Location','best'); end
    if nv == 4, title(['Max ',titv,'  ',ttav]); end
  end ; xlabel(titT);
%--
  axes('position',[.01,.01,.99,.99],'Visible','off');
  T=text(0.5,0.97,titall);
  set(T,'HorizontalAlignment','center','FontSize',12);
  Td=text(0.99,0.01,date);
  set(Td,'HorizontalAlignment','right','FontSize',6);
%---
 end

%-------------------
end

%=========================================================
1	prefix='dynStD';
2	%prefix='iceStD';
3	namA={'g02','b02'};
4	%namA={'r17a'};
5	Nexp=size(namA,2);
6	%-
7
8	% $Header: /u/gcmpack/MITgcm_contrib/jmc_script/grph_StD.m,v 1.8 2015/09/06 18:15:09 jmc Exp $
9	% $Name: $
10
11	nItMx=1e10*ones(1,Nexp); %nItMx(3)=11;
12	%nItMx=2400*ones(1,Nexp);
13	namLg=namA ; namLg=strrep(namLg,'_','\_');
14	undef=123456.7;
15	%-----------
16	%- test if the variable krd is define :
17	if size(who('krd'),1) > 0,
18	fprintf('krd is defined and = %i \n',krd);
19	else
20	fprintf('krd undefined ; set to 1 \n'); krd=1 ;
21	end
22	if krd > 0,
23	%- define list of fields to read in:
24	%listV={'Eta','U','V','W','T','S','DETADT2','RELHUM','Phi'};
25	%listV={'Eta','W','T','S','CONVADJ','DETADT2'};
26	%listV={'Eta','UE_VEL_C','VN_VEL_C','W','T','DETADT2','Phi'};
27	%- or take all them:
28	clear listV ; listV='all_flds';
29	%-----------
30
31	%- start to read the longest record:
32	n=1; rf=-1; if strcmp(char(listV),'all_flds'), rf=0; end
33	[ntA(n),rList,tim,vv1,listV] = ...
34	read_StD(prefix,namA(n),listV);
35	nIt=ntA(n); nk=size(vv1,1); nRg=size(vv1,3);
36	%- set global dims: & load vvA --> vvB
37	nbV=size(listV,2);
38	nrec=nIt; n3d=nk; nReg=nRg;
39	vvA=zeros(n3d,nrec,nReg,5,nbV,Nexp); tiA=zeros(nrec,2,Nexp);
40	vvA(1:nk,1:nIt,1:nRg,:,:,n)=vv1; tiA(1:nIt,:,n)=tim;
41	%----
42	for n=2:Nexp,
43	[ntA(n),rList,tim,vv1,listV] = ...
44	read_StD(prefix,namA(n),listV);
45	nIt=ntA(n); nk=size(vv1,1); nRg=size(vv1,3);
46	if (nrec < nIt),
47	fprintf('\n');
48	error([' Nb of records=',int2str(nIt),' exceeds nrec=',int2str(nrec)]);
49	end
50	if (n3d < nk),
51	fprintf('\n');
52	error([' Nb of Levels=',int2str(nk),' exceeds n3d=',int2str(n3d)]);
53	end
54	vvA(1:nk,1:nIt,1:nRg,:,:,n)=vv1; tiA(1:nIt,:,n)=tim;
55	end;
56	if krd == 2,
57	fprintf('save to "sav_StD.mat" file ...');
58	save('sav_StD.mat','vvA','tiA','ntA','rList','listV');
59	fprintf(' done\n')
60	end
61	elseif krd < 0,
62	fprintf('load from "sav_StD.mat" file ...');
63	load sav_StD
64	fprintf(' done\n'); nbV=size(listV,2);
65	end
66	if krd ~= 0,
67	ttA=squeeze(tiA(:,2,:));
68	ttA=ttA/3600; titT='hrs'; ttA=ttA/24; titT='days';
69	ttA=ttA/30 ; titT='month'; ttA=ttA/12 ; titT='year';
70	end
71	%=========================================================
72
73	ttax1=0 ; ttax2=0 ; ttay=zeros(nbV,2);
74	%-- fixed time axis bound :
75	% ttax1=15.; ttax2=20.;
76	%-- fixed Y axis bound :
77	% ttay(4,:)=[0 0.6];
78	%-----------
79	fprintf('Total length: ntA=');fprintf(' %i ,',ntA); fprintf(' \n');
80	for n=1:Nexp,
81	fprintf(' exp %i : time(d):%10.2f ->%10.2f \n', n,ttA(1,n),ttA(ntA(n),n) );
82	end;
83	%--
84
85	list_on=zeros(1,nbV);
86	nbG=18;
87	nbG=min(nbG,nbV); list_on(1:nbG)=1 ;
88	%if nbG < nbV, list_on=0; list_on(nbG:nbV)=1; end %- to get 2nd set of plots
89	%list_on(1:6)=[1 1 1 1 1 1];
90
91	isA=ones(1,Nexp); ieA=ntA;
92	%- limit the length : for search of isA <->1500y: find(ttA(:,2) == 1500)
93	%isA=isA*31 ; % drop the 1rst mnth (1 Monitor/d)
94	isA=isA*2 ; % drop the 1rst mnth (1 Monitor/30d)
95	%isA(1)=31 ; isA(2)=4 ; % drop the 1rst mnth
96	%ieA(:)=240; %isA(:)=1;
97
98	linA(1,:)='k-'; % ieA(1)=60 ; % ieA(1)=1152 ;
99	linA(2,:)='b-';
100	linA(3,:)='r-';
101	linA(4,:)='g-';
102	linA(5,:)='m-';
103	linA(6,:)='c-';
104
105	ieA=min(ieA,nItMx);
106	%titall='AIM , Cubic-G (32x32) , cpl-FM Forcing' ;
107	titall='Global Ocean, Cubic-G (32x32) , CORE Forc (2)' ;
108	%titall='Dyncore test-case 5 (cs-32)' ;
109
110	%=========================================================
111	ng=0; fxb=100; fyb=60; fdx=100; fdy=40; fsc=1.;
112	%fyb=-360; fxb=-2600; fdy=60; fsc=1.5;
113
114	for jv=1:nbV,
115	%-------------------
116	flag=list_on(jv); kl=0;
117	vv1=vvA(:,:,:,:,jv,:); namV=char(listV(jv));
118	titv=strrep(namV,'_','\_');
119	%if strcmp(namV,'Eta'), vv1=vv1/100; titv='Eta [mb]'; end
120	%if strcmp(namV,'T'), kl=1; end % <-- to get surf.Temp
121	if kl > 0, titv=[titv,'\_',int2str(kl)];
122	fprintf([' var= ',namV,' at level k= %i \n'],kl);
123	end
124
125	if flag == 1
126	%--
127	ng=ng+1;
128	%- reset "ng" to jv for fix fig number (independent of which one is in list_on):
129	%ng=jv;
130	figure(ng); set(ng,'position',[fxb+fdxng fyb+fdyng [500 700]*fsc]);clf;
131	var=squeeze(vv1(1+kl,:,1,:,:)); dd=zeros(5,Nexp); av=zeros(5,Nexp);
132	for n=1:Nexp,
133	tmp=var(isA(n):ieA(n),:,n); [I]=find(tmp(:,1)==undef); tmp(I,:)=0;
134	av(:,n)=sum(tmp);
135	if length(I) < size(tmp,1), av(:,n)=av(:,n)/(size(tmp,1)-length(I)); end
136	tmp(I,:)=NaN; dd(:,n)=max(tmp)-min(tmp);
137	end
138	var(find(var==undef))=NaN;
139	for nv=1:4,
140	subplot(410+nv); ttmn=' Mx-mn:'; ttav=' Av:';
141	for n=1:Nexp,
142	plot(ttA(isA(n):ieA(n),n),var(isA(n):ieA(n),nv,n),linA(n,:));
143	if n == 1, hold on ; end ;
144	ttmn=sprintf([ttmn,' %2.1e ;'],dd(nv,n));
145	ttav=sprintf([ttav,' %3.2e ;'],av(nv,n));
146	end ; hold off ;
147	AA=axis ; dAA=AA(4)-AA(3);
148	if AA(3)*AA(4) <= 0, AA(3)=min(AA(3),-dAA/10); AA(4)=max(AA(4),dAA/10); end
149	if ttax1 < ttax2, AA(1)=ttax1; AA(2)=ttax2; end;
150	axis(AA); grid ;
151	if nv == 1, title(['Avr ',titv,' ',ttmn]); end
152	if nv == 2, title(['Std-Dev ',titv,' ',ttav]); end
153	if nv == 3, title(['min ',titv,' ',ttav]);
154	legend(namLg,'Location','best'); end
155	if nv == 4, title(['Max ',titv,' ',ttav]); end
156	end ; xlabel(titT);
157	%--
158	axes('position',[.01,.01,.99,.99],'Visible','off');
159	T=text(0.5,0.97,titall);
160	set(T,'HorizontalAlignment','center','FontSize',12);
161	Td=text(0.99,0.01,date);
162	set(Td,'HorizontalAlignment','right','FontSize',6);
163	%---
164	end
165
166	%-------------------
167	end
168
169	%=========================================================