/[MITgcm]/MITgcm_contrib/jmc_script/plot_StD.m
ViewVC logotype

Contents of /MITgcm_contrib/jmc_script/plot_StD.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.4 - (show annotations) (download)
Thu Mar 5 20:57:26 2015 UTC (9 years, 1 month ago) by jmc
Branch: MAIN
Changes since 1.3: +4 -1 lines
ignore missing values

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

  ViewVC Help
Powered by ViewVC 1.1.22