MITgcm_contrib/jmc_script/plot_StD.m

 prefix='dynStD';
%prefix='oceStD';
 pCoords=0;
 namA='c06';
 Nexp=1; Nc=size(namA,2);
 nAvr=1;
%- to plot annual mean instead of 10.d aver:
%nAvr=36;
%--

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

nItMx=1e10*ones(1,Nexp); %nItMx(3)=11;
%nItMx=60*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','U','V','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,kList] = ...
    read_StD(prefix,namA(n,:),listV);
  nIt=ntA(n); nk=size(vv1,1); nRg=size(vv1,3);
  vv1(find(vv1==undef))=NaN;
%- 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);
  vv1(find(vv1==undef))=NaN;
  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';
%- change to plot annual mean:
 if nAvr > 1 & rem(nrec,nAvr)==0,
  nrec=nrec/nAvr; ntA=ntA/nAvr;
  vvA=reshape(   vvA,[n3d nAvr nrec nReg 5 nbV Nexp]);
  vvA=reshape(mean(vvA,2),[n3d nrec nReg 5 nbV Nexp]);
  ttA=reshape(   ttA,[nAvr nrec Nexp]);
  ttA=reshape(mean(ttA,1),[nrec Nexp]);
 end
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=10;
nbG=min(nbG,nbV); list_on(1:nbG)=1 ;
%list_on(1:6)=[1 1 1 1 1 1];
%list_on(5:7)=0;

isA=ones(1,Nexp); ieA=ntA;
%- limit the length : for search of isA <->1500y: find(ttA(:,2) == 1500)
%isA=isA*3 ; % drop the 1rst mnth (1 Monitor/10.d)
%ieA(:)=360; 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=['Exp: ',namLg(1,:)];

%=========================================================

%-default: dxRed=0; dyRed=0; dxB=0.1; dyB=0.9;
dxRed=0; dyRed=0.03; dxB=0.02; dyB=0.9;
[xyP,xyB]=def_subP(-4,dxRed,dyRed,dxB,dyB);
xyP(:,2)=xyP(:,2)+0.010;
xyB(:,2)=xyB(:,2)+0.010;

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

for ng=1:nbV,
%-------------------
 yax=[1:nk-1]; if pCoords == 0, yax=-[1:nk-1]; end
 flag=list_on(ng);
 vv1=vvA(:,:,:,:,ng,:); namV=char(listV(ng));
 titv=strrep(namV,'_','\_');
 if strcmp(namV,'Eta') & pCoords == 1, vv1=vv1/100; titv='Eta [mb]'; end
 if strcmp(namV,'T') & pCoords == 1,
   namfil=['../res_',namA(1:end),'/RC']; D=dir([namfil,'.data']);
   if size(D,1) == 1,
     rC=rdmds(namfil);
     fprintf(' convert Pot.Temp to Temp.:');
%    fprintf(' convert Pot.Temp to Temp.:'); fprintf(' %i',size(vv1));
     kappa=2/7; facP=squeeze(rC)/1.e+5; facP=facP.^kappa;
     var=facP*ones(1,nrec*4*Nexp); var=reshape(var,[nk-1 nrec 1 4 Nexp]);
     vv1([2:nk],:,1,[1:4],:)=vv1([2:nk],:,1,[1:4],:).*var;
%    for k=2:nk, vv1(k,:,1,[1:4],:)=vv1(k,:,1,[1:4],:)*facP(k-1); end
     fprintf(' OK\n');
   else
     fprintf(' no file: %s\n',namfil);
   end
 end

 if flag == 1
%--
  figure(ng); set(ng,'position',[fxb+fdx*ng fyb+fdy*ng [500 700]*fsc]);clf;
  colormap jet
  if kList(ng) == 1,
   var=squeeze(vv1(1,:,1,:,:));
   dd=squeeze(max(var)-min(var)); av=squeeze(mean(var));
   if Nexp == 1, av=av'; dd=dd'; end ;
   for nv=1:4,
    axes('position',xyP(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]); xlabel(titT); end
    if nv == 2, title(['Std-Dev ',titv,'  ',ttav]); end
    if nv == 3, title(['min ',titv,'  ',ttav]);
                legend(namLg(1:Nexp,:),'Location','best'); end
    if nv == 4, title(['Max ',titv,'  ',ttav]); end
   %if nv == 2, title(['Del-2 ',titv,'  ',ttav]); end
   end ; %xlabel(titT);
  else
   n=1; k1=2;
   if strcmp(namV,'CONVADJ') || strcmp(namV,'DRHODR'),
     k1=3; yax=yax(2:nk-1);
   end
   for nv=1:4,
    axes('position',xyP(nv,:));
    var=squeeze(vv1(k1:nk,:,1,nv,n))'; mnV=min(var(:)); MxV=max(var(:));
    ccv=c_levs(mnV,MxV,-12); %ccv=c_levs(mnV,MxV,-20);
    if MxV > mnV,
     [cs,h]=contour(ttA(isA(n):ieA(n),n),yax,var(isA(n):ieA(n),:)',ccv);
    %clabel(cs);isoline0(h);
     BB=colorbar; set(BB,'Position',xyB(nv,:));
    end
    if nv == 1, title(['Avr ',titv]); xlabel(titT); end
    if nv == 2, title(['Std-Dev ',titv]); end
    if nv == 3, title(['min ',titv]); ; end
    if nv == 4, title(['Max ',titv]); end
    AA=axis; dAA=AA(4)-AA(3);
    ttmn=sprintf('mn= %4.3g , Mx= %4.3g',mnV,MxV);
    text(AA(1)*.4+AA(2)*.6,AA(3)-0.20*dAA,ttmn);
   end ; %xlabel(titT);
  end
%--
  axes('position',[.01,.01,.99,.99],'Visible','off');
  T=text(0.2,0.98,titall);
  set(T,'HorizontalAlignment','center','FontSize',12);
  Td=text(0.01,0.01,date);
  set(Td,'HorizontalAlignment','left','FontSize',6);
%---
 end

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

%=========================================================
1	jmc	1.2	prefix='dynStD';
2	jmc	1.5	%prefix='oceStD';
3	jmc	1.3	pCoords=0;
4			namA='c06';
5	jmc	1.1	Nexp=1; Nc=size(namA,2);
6	jmc	1.5	nAvr=1;
7			%- to plot annual mean instead of 10.d aver:
8			%nAvr=36;
9	jmc	1.1	%--
10
11	jmc	1.6	% $Header: /u/gcmpack/MITgcm_contrib/jmc_script/plot_StD.m,v 1.5 2015/09/06 18:15:09 jmc Exp $
12	jmc	1.1	% $Name: $
13
14			nItMx=1e10*ones(1,Nexp); %nItMx(3)=11;
15			%nItMx=60*ones(1,Nexp);
16			namLg=namA ; namLg=strrep(namLg,'_','\_');
17	jmc	1.4	undef=123456.7;
18	jmc	1.1	%-----------
19			%- test if the variable krd is define :
20			if size(who('krd'),1) > 0,
21			fprintf('krd is defined and = %i \n',krd);
22			else
23			fprintf('krd undefined ; set to 1 \n'); krd=1 ;
24			end
25			if krd > 0,
26			%- define list of fields to read in:
27			%listV={'Eta','U','V','W','T','S','DETADT2','RELHUM','Phi'};
28			%listV={'Eta','U','V','W','T','S','CONVADJ','DETADT2'};
29			%listV={'Eta','UE_VEL_C','VN_VEL_C','W','T','DETADT2','Phi'};
30			%- or take all them:
31			clear listV ; listV='all_flds';
32			%-----------
33
34			%- start to read the longest record:
35			n=1; rf=-1; if strcmp(char(listV),'all_flds'), rf=0; end
36			[ntA(n),rList,tim,vv1,listV,kList] = ...
37			read_StD(prefix,namA(n,:),listV);
38			nIt=ntA(n); nk=size(vv1,1); nRg=size(vv1,3);
39	jmc	1.4	vv1(find(vv1==undef))=NaN;
40	jmc	1.1	%- set global dims: & load vvA --> vvB
41			nbV=size(listV,2);
42			nrec=nIt; n3d=nk; nReg=nRg;
43			vvA=zeros(n3d,nrec,nReg,5,nbV,Nexp); tiA=zeros(nrec,2,Nexp);
44			vvA(1:nk,1:nIt,1:nRg,:,:,n)=vv1; tiA(1:nIt,:,n)=tim;
45			%----
46			for n=2:Nexp,
47			[ntA(n),rList,tim,vv1,listV] = ...
48			read_StD(prefix,namA(n,:),listV);
49			nIt=ntA(n); nk=size(vv1,1); nRg=size(vv1,3);
50	jmc	1.4	vv1(find(vv1==undef))=NaN;
51	jmc	1.1	if (nrec < nIt),
52			fprintf('\n');
53			error([' Nb of records=',int2str(nIt),' exceeds nrec=',int2str(nrec)]);
54			end
55			if (n3d < nk),
56			fprintf('\n');
57			error([' Nb of Levels=',int2str(nk),' exceeds n3d=',int2str(n3d)]);
58			end
59			vvA(1:nk,1:nIt,1:nRg,:,:,n)=vv1; tiA(1:nIt,:,n)=tim;
60			end;
61			if krd == 2,
62			fprintf('save to "sav_StD.mat" file ...');
63			save('sav_StD.mat','vvA','tiA','ntA','rList','listV');
64			fprintf(' done\n')
65	jmc	1.3	end
66	jmc	1.1	elseif krd < 0,
67			fprintf('load from "sav_StD.mat" file ...');
68			load sav_StD
69			fprintf(' done\n'); nbV=size(listV,2);
70			end
71			if krd ~= 0,
72			ttA=squeeze(tiA(:,2,:));
73	jmc	1.3	ttA=ttA/3600; titT='hrs'; ttA=ttA/24 ; titT='days';
74	jmc	1.2	ttA=ttA/30 ; titT='month'; ttA=ttA/12 ; titT='year';
75	jmc	1.5	%- change to plot annual mean:
76			if nAvr > 1 & rem(nrec,nAvr)==0,
77			nrec=nrec/nAvr; ntA=ntA/nAvr;
78			vvA=reshape( vvA,[n3d nAvr nrec nReg 5 nbV Nexp]);
79			vvA=reshape(mean(vvA,2),[n3d nrec nReg 5 nbV Nexp]);
80			ttA=reshape( ttA,[nAvr nrec Nexp]);
81			ttA=reshape(mean(ttA,1),[nrec Nexp]);
82			end
83	jmc	1.1	end
84			%=========================================================
85
86			ttax1=0 ; ttax2=0 ; ttay=zeros(nbV,2);
87			%-- fixed time axis bound :
88			% ttax1=15.; ttax2=20.;
89			%-- fixed Y axis bound :
90			% ttay(4,:)=[0 0.6];
91			%-----------
92			fprintf('Total length: ntA=');fprintf(' %i ,',ntA); fprintf(' \n');
93			for n=1:Nexp,
94			fprintf(' exp %i : time(d):%10.2f ->%10.2f \n', n,ttA(1,n),ttA(ntA(n),n) );
95			end;
96	jmc	1.5	%-
97	jmc	1.1
98			list_on=zeros(1,nbV);
99	jmc	1.5	nbG=10;
100	jmc	1.1	nbG=min(nbG,nbV); list_on(1:nbG)=1 ;
101			%list_on(1:6)=[1 1 1 1 1 1];
102	jmc	1.3	%list_on(5:7)=0;
103	jmc	1.1
104			isA=ones(1,Nexp); ieA=ntA;
105	jmc	1.3	%- limit the length : for search of isA <->1500y: find(ttA(:,2) == 1500)
106			%isA=isA*3 ; % drop the 1rst mnth (1 Monitor/10.d)
107	jmc	1.1	%ieA(:)=360; isA(:)=1;
108
109			linA(1,:)='k-'; % ieA(1)=60 ; % ieA(1)=1152 ;
110			linA(2,:)='b-';
111			linA(3,:)='r-';
112			linA(4,:)='g-';
113			linA(5,:)='m-';
114			linA(6,:)='c-';
115
116			ieA=min(ieA,nItMx);
117			titall=['Exp: ',namLg(1,:)];
118
119			%=========================================================
120	jmc	1.3
121			%-default: dxRed=0; dyRed=0; dxB=0.1; dyB=0.9;
122			dxRed=0; dyRed=0.03; dxB=0.02; dyB=0.9;
123			[xyP,xyB]=def_subP(-4,dxRed,dyRed,dxB,dyB);
124			xyP(:,2)=xyP(:,2)+0.010;
125			xyB(:,2)=xyB(:,2)+0.010;
126
127	jmc	1.6	fxb=100; fyb=60; fdx=100; fdy=40; fsc=1.;
128			%fyb=-360; fxb=-2600; fdy=60; fsc=1.5;
129	jmc	1.5
130	jmc	1.1	for ng=1:nbV,
131			%-------------------
132	jmc	1.3	yax=[1:nk-1]; if pCoords == 0, yax=-[1:nk-1]; end
133			flag=list_on(ng);
134	jmc	1.1	vv1=vvA(:,:,:,:,ng,:); namV=char(listV(ng));
135			titv=strrep(namV,'_','\_');
136	jmc	1.3	if strcmp(namV,'Eta') & pCoords == 1, vv1=vv1/100; titv='Eta [mb]'; end
137			if strcmp(namV,'T') & pCoords == 1,
138			namfil=['../res_',namA(1:end),'/RC']; D=dir([namfil,'.data']);
139	jmc	1.1	if size(D,1) == 1,
140			rC=rdmds(namfil);
141			fprintf(' convert Pot.Temp to Temp.:');
142			% fprintf(' convert Pot.Temp to Temp.:'); fprintf(' %i',size(vv1));
143			kappa=2/7; facP=squeeze(rC)/1.e+5; facP=facP.^kappa;
144			var=facPones(1,nrec4*Nexp); var=reshape(var,[nk-1 nrec 1 4 Nexp]);
145			vv1([2:nk],:,1,[1:4],:)=vv1([2:nk],:,1,[1:4],:).*var;
146			% for k=2:nk, vv1(k,:,1,[1:4],:)=vv1(k,:,1,[1:4],:)*facP(k-1); end
147			fprintf(' OK\n');
148			else
149			fprintf(' no file: %s\n',namfil);
150			end
151			end
152
153			if flag == 1
154			%--
155	jmc	1.6	figure(ng); set(ng,'position',[fxb+fdxng fyb+fdyng [500 700]*fsc]);clf;
156	jmc	1.5	colormap jet
157	jmc	1.1	if kList(ng) == 1,
158	jmc	1.3	var=squeeze(vv1(1,:,1,:,:));
159	jmc	1.1	dd=squeeze(max(var)-min(var)); av=squeeze(mean(var));
160			if Nexp == 1, av=av'; dd=dd'; end ;
161			for nv=1:4,
162	jmc	1.3	axes('position',xyP(nv,:)); ttmn=' Mx-mn:'; ttav=' Av:';
163	jmc	1.1	for n=1:Nexp,
164			plot(ttA(isA(n):ieA(n),n),var(isA(n):ieA(n),nv,n),linA(n,:));
165			if n == 1, hold on ; end ;
166			ttmn=sprintf([ttmn,' %2.1e ;'],dd(nv,n));
167			ttav=sprintf([ttav,' %3.2e ;'],av(nv,n));
168			end ; hold off ;
169			AA=axis ; dAA=AA(4)-AA(3);
170			if AA(3)*AA(4) <= 0, AA(3)=min(AA(3),-dAA/10); AA(4)=max(AA(4),dAA/10); end
171	jmc	1.3	if ttax1 < ttax2, AA(1)=ttax1; AA(2)=ttax2; end;
172	jmc	1.1	axis(AA); grid ;
173	jmc	1.3	if nv == 1, title(['Avr ',titv,' ',ttmn]); xlabel(titT); end
174	jmc	1.1	if nv == 2, title(['Std-Dev ',titv,' ',ttav]); end
175	jmc	1.6	if nv == 3, title(['min ',titv,' ',ttav]);
176			legend(namLg(1:Nexp,:),'Location','best'); end
177	jmc	1.1	if nv == 4, title(['Max ',titv,' ',ttav]); end
178			%if nv == 2, title(['Del-2 ',titv,' ',ttav]); end
179	jmc	1.3	end ; %xlabel(titT);
180	jmc	1.1	else
181	jmc	1.3	n=1; k1=2;
182			if strcmp(namV,'CONVADJ') \|\| strcmp(namV,'DRHODR'),
183			k1=3; yax=yax(2:nk-1);
184			end
185	jmc	1.1	for nv=1:4,
186	jmc	1.3	axes('position',xyP(nv,:));
187			var=squeeze(vv1(k1:nk,:,1,nv,n))'; mnV=min(var(:)); MxV=max(var(:));
188			ccv=c_levs(mnV,MxV,-12); %ccv=c_levs(mnV,MxV,-20);
189			if MxV > mnV,
190			[cs,h]=contour(ttA(isA(n):ieA(n),n),yax,var(isA(n):ieA(n),:)',ccv);
191			%clabel(cs);isoline0(h);
192			BB=colorbar; set(BB,'Position',xyB(nv,:));
193			end
194			if nv == 1, title(['Avr ',titv]); xlabel(titT); end
195	jmc	1.1	if nv == 2, title(['Std-Dev ',titv]); end
196			if nv == 3, title(['min ',titv]); ; end
197			if nv == 4, title(['Max ',titv]); end
198	jmc	1.3	AA=axis; dAA=AA(4)-AA(3);
199	jmc	1.1	ttmn=sprintf('mn= %4.3g , Mx= %4.3g',mnV,MxV);
200	jmc	1.3	text(AA(1).4+AA(2).6,AA(3)-0.20*dAA,ttmn);
201			end ; %xlabel(titT);
202	jmc	1.1	end
203			%--
204			axes('position',[.01,.01,.99,.99],'Visible','off');
205	jmc	1.3	T=text(0.2,0.98,titall);
206	jmc	1.1	set(T,'HorizontalAlignment','center','FontSize',12);
207	jmc	1.3	Td=text(0.01,0.01,date);
208			set(Td,'HorizontalAlignment','left','FontSize',6);
209			%---
210	jmc	1.1	end
211
212	jmc	1.3	%-------------------
213	jmc	1.1	end
214
215			%=========================================================