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

Contents of /MITgcm_contrib/jmc_script/grph_MON.m

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


Revision 1.4 - (show annotations) (download)
Wed Oct 19 15:29:58 2022 UTC (2 years, 8 months ago) by jmc
Branch: MAIN
CVS Tags: HEAD
Changes since 1.3: +13 -5 lines
add more comments + fix W-vel units

1
2 %namA={'s05'};
3 namA={'s19c','s19'};
4 %namA={'s19c','s19','s19a'};
5 %namA={'s18c','s18b','s18a','s01'};
6 %namA={'s19a','s19','s18c','s18b','s18a','s04'};
7 %namA={'s19c','s19','s19a','s18b','s18a','s01','s04'};
8 Nexp=size(namA,2); list_log=zeros(1,10);
9 %list_log=[1 0 0 1 1 1 1 1 0 0] ;
10 %-- set type of monitor output files: ncF=0 : ASCII output file ; ncF=1 : NetCDF file
11 ncF=zeros(1,Nexp);
12
13 % $Header: /u/gcmpack/MITgcm_contrib/jmc_script/grph_MON.m,v 1.3 2022/02/08 17:10:21 jmc Exp $
14 % $Name: $
15
16 nItMx=1e10*ones(1,Nexp); %nItMx(3)=11;
17 namLg=namA ; namLg=strrep(namLg,'_','\_');
18 %-----------
19 %- ngEn = Nb of Energy plot: = 2 (ke: Mx+Av) ; = 3 (+AM) or = 4 (add Pe ?) ;
20 ngEn=2;
21 %- in case PE is missing from Monitor, use Mean-Eta and this "gdH" ratio to get PE
22 gdH=9.81/1000 ; % ratio gravity / mean_H (ocean) ; mean_Bo / mean_Pground (atmos)
23 %- test if the variable krd is define :
24 if size(who('krd'),1) > 0,
25 fprintf('krd is defined and = %i \n',krd);
26 else
27 fprintf('krd undefined ; set to 1 \n'); krd=1 ;
28 end
29 if krd == 1,
30 % list_on : controls which field to read-in (and plot):
31 % 1 : KE ; 2 : Eta ; 3,4 : T,S ; 5,6 : U,V ; 7 : W ; 8 : CFL ; 9 : Vort ; 10 : "sc"
32 %list_on=[1 1 1 1 0 0 1 1 0 0] ;
33 list_on=[1 1 1 1 1 1 1 1 0 0] ;
34 %list_on=[1 1 0 0 0 0 0 0 1 1] ;
35 %- Warning: list_log works with fig # and do not match list_on !
36 %list_log=zeros(1,length(list_on));
37 %-----------
38
39 %- start to read the longest record:
40 n=1;
41 if ncF(n) == 1,
42 [ntA(n),ttA(:,n),keA(:,:,n),etA(:,:,n),tmA(:,:,n),smA(:,:,n), ...
43 umA(:,:,n),vmA(:,:,n),wmA(:,:,n),cfA(:,:,n),zmA(:,:,n),scA(:,:,n)] = ...
44 readnc_MON(char(namA(n)),list_on);
45 else
46 [ntA(n),ttA(:,n),keA(:,:,n),etA(:,:,n),tmA(:,:,n),smA(:,:,n), ...
47 umA(:,:,n),vmA(:,:,n),wmA(:,:,n),cfA(:,:,n),zmA(:,:,n),scA(:,:,n)] = ...
48 read_MON(char(namA(n)),list_on);
49 end
50 nrec=ntA(1);
51 for n=2:Nexp,
52 if ncF(n) == 1,
53 [ntA(n),ttA(:,n),keA(:,:,n),etA(:,:,n),tmA(:,:,n),smA(:,:,n), ...
54 umA(:,:,n),vmA(:,:,n),wmA(:,:,n),cfA(:,:,n),zmA(:,:,n),scA(:,:,n)] = ...
55 readnc_MON(char(namA(n)),list_on,nrec);
56 else
57 [ntA(n),ttA(:,n),keA(:,:,n),etA(:,:,n),tmA(:,:,n),smA(:,:,n), ...
58 umA(:,:,n),vmA(:,:,n),wmA(:,:,n),cfA(:,:,n),zmA(:,:,n),scA(:,:,n)] = ...
59 read_MON(char(namA(n)),list_on,nrec);
60 end
61 end;
62 vol=squeeze(keA(:,3,:)); if vol(1,1) ~= 0,
63 for n=1:Nexp, vol(:,n)=vol(:,n)/vol(1,n); end; end; % <- volume ratio = V/Vo
64 %vol=squeeze(etA(:,3,:)); vol=(vol+1000)/1000 ;
65 %------------------
66 if ngEn == 4,
67 tt2e=' (+) '; tt3e=' (+) '; tt4e=' (+) ';
68 for n=1:Nexp,
69 keA(:,1,n)=etA(:,2,n)-etA(:,1,n);
70 ddKe=max(keA(:,4,n));
71 if ddKe > 0, keA(:,3,n)=keA(:,4,n); %- take PE from Monitor (since it's there)
72 else keA(:,3,n)=etA(:,4,n).*etA(:,4,n);
73 keA(:,3,n)=keA(:,3,n)*gdH/2; end
74 keA(:,4,n)=keA(:,3,n)+keA(:,2,n);
75 tt2e=sprintf([tt2e,' %1.1e ;'],keA(1,2,n)); keA(:,2,n)=keA(:,2,n)-keA(1,2,n);
76 tt3e=sprintf([tt3e,' %1.1e ;'],keA(1,3,n)); keA(:,3,n)=keA(:,3,n)-keA(1,3,n);
77 tt4e=sprintf([tt4e,' %1.1e ;'],keA(1,4,n)); keA(:,4,n)=keA(:,4,n)-keA(1,4,n);
78 end
79 else tt2e=' ' ; tt3e=' ' ; tt4e=' ' ; end;
80 %-- set time units:
81 titT='s' ; ttA=ttA/60; titT='mn' ;
82 ttA=ttA/3600; titT='hrs' ; ttA=ttA/24; titT='days';
83 ttA=ttA/30 ; titT='month'; ttA=ttA/12 ; titT='year';
84 end
85 %-----------
86
87 ttax1=0 ; ttax2=0 ; ttay=zeros(size(list_on,2),2);
88 %-- fixed time axis bound :
89 % ttax1=3; ttax2=3.3;
90 %-- fixed Y axis bound :
91 % ttay(4,:)=[0 0.6];
92 %=========================================================
93 fprintf('Total length: ntA=');fprintf(' %i ,',ntA); fprintf(' \n');
94 for n=1:Nexp,
95 fprintf(' exp %i : time(d):%10.2f ->%10.2f \n', n,ttA(1,n),ttA(ntA(n),n) );
96 end;
97 %--
98
99 isA=ones(1,Nexp); ieA=ntA ;
100 %- limit the length : for search of isA <->1500y: find(ttA(:,2) == 1500)
101 %isA=isA*31 ; % drop the 1rst mnth
102 %ieA(1)=1681; ieA(2)=281;
103
104 linA={'k-','b-','r-','g-','m-','c-'};
105 %linA={'k-','k--','b-','r-','g-','m-','c-'};
106
107 ieA=min(ieA,nItMx);
108 %titall='AIM , Cubic-G (32x32) , NCEP Forc (2)' ;
109 titall='S.Ocean Section (320x50), CORE Forcing' ;
110 %titall='Enceladus with Ice, Lat-Lon half-sphere' ;
111
112 %=========================================================
113 ng=0; fxb=100; fyb=60; fsc=1.5;
114 %fxb=-2600; %fyb=160;
115
116 for ng=1:size(list_on,2)
117 %-------------------
118 flag=0;
119 if ng == 1, flag=2*list_on(1) ; titv='Ke [m2/s2]' ; vvA=keA ; end
120 if ng == 2, flag=1*list_on(3) ; titv='Theta [K]' ; vvA=tmA ; end
121 if ng == 3, flag=1*list_on(4) ; titv='Salin [psu]' ; vvA=smA ; end
122 %if ng == 3, flag=1*list_on(4) ; titv='Spe.Hum [g/kg]';vvA=smA; end
123 if ng == 4, flag=3*list_on(8) ; titv='CFL-max' ; vvA=cfA ; end
124 if ng == 5, flag=1*list_on(2) ; titv='Eta [m]' ; vvA=etA ; end
125 if ng == 6, flag=1*list_on(7) ; titv='W-vel [m/s]' ; vvA=wmA ; end
126 %if ng == 5, flag=1*list_on(2) ; titv='Eta [mb]'; vvA=etA/100 ; end
127 %if ng == 6, flag=1*list_on(7) ; titv='W-vel [Pa/s]'; vvA=wmA ; end
128 if ng == 7, flag=1*list_on(5) ; titv='U-vel [m/s]' ; vvA=umA ; end
129 if ng == 8, flag=1*list_on(6) ; titv='V-vel [m/s]' ; vvA=vmA ; end
130 if ng == 9, flag=4*list_on(9) ; titv= 'Vort [s-1]' ; vvA=zmA ; end
131
132 if flag == 1
133 %--
134 figure(ng); set(ng,'position',[fxb+100*ng fyb+40*ng [500 700]*fsc]);clf;
135 dd=squeeze(max(vvA)-min(vvA)); av=squeeze(mean(vvA));
136 if Nexp == 1, av=av'; dd=dd'; end ;
137 for nv=1:4,
138 subplot(410+nv); ttmn=' Mx-mn:'; ttav=' Av:'; ttyax='';
139 if list_log(ng) > 0,
140 var=zeros(nrec,Nexp); for n=1:Nexp, var(isA(n):ieA(n),n)=vvA(isA(n):ieA(n),nv,n); end
141 if min(var(:))*max(var(:)) >= 0, msk=var; ttyax='log';
142 var(find(msk==0))=1; var=log10(abs(var)); var(find(msk==0))=NaN;
143 end
144 else var=squeeze(vvA(:,nv,:)); end
145 for n=1:Nexp,
146 plot(ttA(isA(n):ieA(n),n),var(isA(n):ieA(n),n),char(linA(n)));
147 if n == 1, hold on ; end ;
148 ttmn=sprintf([ttmn,' %2.1e ;'],dd(nv,n));
149 ttav=sprintf([ttav,' %3.2e ;'],av(nv,n));
150 end ; hold off ;
151 if ttax1 < ttax2, AA=axis; axis([ttax1 ttax2 AA(3:4)]); end;
152 AA=axis ; dAA=AA(4)-AA(3);
153 if AA(3)*AA(4) <= 0, AA(3)=min(AA(3),-dAA/10); AA(4)=max(AA(4),dAA/10); end
154 axis(AA); grid ;
155 %AA=axis; text(AA(1)*.9+AA(2)*.1,AA(3)*0.1+AA(4)*0.9,ttmn);
156 if nv == 1, title(['min ',titv,' ',ttav]);
157 legend(namLg,'Location','best'); end
158 if nv == 2, title(['Max ',titv,' ',ttav]); end
159 if nv == 3, title(['Avr ',titv,' ',ttmn]); end
160 if nv == 4, title(['Std-Dev ',titv,' ',ttav]); end
161 %if nv == 4, title(['Del-2 ',titv,' ',ttav]); end
162 if length(ttyax) > 0, ylabel(ttyax); end
163 end ; xlabel(titT);
164 %--
165 axes('position',[.01,.01,.99,.99],'Visible','off');
166 T=text(0.5,0.97,titall);
167 set(T,'HorizontalAlignment','center','FontSize',12);
168 put_date;
169 %---
170 end
171
172 if flag == 2
173 %--
174 ngEx = ngEn ; if ngEn == 3, ngEx = 4; end ; ngEk = ngEx;
175 if list_on(10) == 1, ngEx = 4; ngEk = 2; end
176 if ngEx == 4,
177 figure(ng); set(ng,'position',[fxb+100*ng fyb+40*ng [500 700]*fsc]);clf;
178 else
179 figure(ng); set(ng,'position',[fxb+100*ng fyb+140*ng [500 500]*fsc]);clf;
180 end;
181 dd=squeeze(max(vvA)-min(vvA)); av=squeeze(mean(vvA));
182 if Nexp == 1, av=av'; dd=dd'; end ;
183 % for n=1:Nexp,
184 % for n=1:0,
185 % nt1=find(ttA(:,n)==4);
186 % vv1=min(vvA([nt1:ntA(n)],2,n));[I1]=find(vvA(:,2,n)==vv1);
187 % vv2=max(vvA([nt1:ntA(n)],3,n));[I2]=find(vvA(:,3,n)==vv2);
188 % if length(I1) ~= 1, I1, return ; end
189 % if length(I2) ~= 1, I2, return ; end
190 % fprintf('n= %i min(KE): %i %7.6f ; max(PE): %i %7.6f \n',n,I1,vv1,I2,vv2);
191 % fprintf(' %i KE: %7.6f ; PE: %7.6f ', I1-1,vvA(I1-1,2,n),vvA(I1-1,3,n));
192 % fprintf(' ; %i KE: %7.6f ; PE: %7.6f \n',I1+1,vvA(I1+1,2,n),vvA(I1+1,3,n));
193 % if n == 1, ii1=I1 ; end
194 % end
195 for np=1:ngEk,
196 nv=np; if ngEn == 3, nv=np+1; end
197 subplot(100*ngEx+10+np); ttmn=' Mx-mn:'; ttav=' Av:'; ttyax='';
198 if list_log(ng) > 0,
199 var=zeros(nrec,Nexp); for n=1:Nexp, var(isA(n):ieA(n),n)=vvA(isA(n):ieA(n),nv,n); end
200 if min(var(:))*max(var(:)) >= 0, msk=var; ttyax='log';
201 var(find(msk==0))=1; var=log10(abs(var)); var(find(msk==0))=NaN;
202 end
203 else var=squeeze(vvA(:,nv,:)); end
204 for n=1:Nexp,
205 plot(ttA(isA(n):ieA(n),n),var(isA(n):ieA(n),n),char(linA(n)));
206 % LL(n)=plot(ttA(isA(n):ieA(n),n),var(isA(n):ieA(n),n),char(linA(n)));
207 if n == 1, hold on ; end ;
208 ttmn=sprintf([ttmn,' %2.1e ;'],dd(nv,n));
209 ttav=sprintf([ttav,' %3.2e ;'],av(nv,n));
210 end ; hold off ;
211 % set(LL(1),'LineWidth',2);
212 if ttax1 < ttax2, AA=axis; axis([ttax1 ttax2 AA(3:4)]); end;
213 grid ;
214 if np == 1, legend(namLg,'Location','best'); end
215 if nv == 1, title(['Max ',titv,' ',ttav]); end
216 if nv == 2, title(['Avr ',titv,' ',ttmn]); end
217 if ngEn == 3,
218 if nv == 3, title(['AM-eta [kg/s] ',ttmn]); end
219 if nv == 4, title(['AM-Uzo [kg/s] ',ttmn]); end
220 if nv == 5, title(['AM-Tot [kg/s] ',ttmn]); end
221 else
222 if nv == 3, title('PE-eta [m2/s2]'); end
223 if nv == 4, title('Tot-En [m2/s2]'); end
224 end
225 %if nv == 1, title('eta Max-min'); end % legend(namLg,0); end
226 %if nv == 2, title(['Avr ',titv,' ',tt2e]);legend(namLg,0); end
227 %if nv == 3, title(['PE-eta [m2/s2]',' ',tt3e]); end
228 %if nv == 4, title(['Tot-En [m2/s2]',' ',tt4e]); end
229 if length(ttyax) > 0, ylabel(ttyax); end
230 end ;
231 if ngEn < 4 & list_on(10) == 1,
232 titv1='surf-Cor : ' ; vvA=scA ;
233 dd=squeeze(max(vvA)-min(vvA)); av=squeeze(mean(vvA));
234 if Nexp == 1, av=av'; dd=dd'; end ;
235 for nv=1:ngEx-ngEk,
236 vvM=mean(vvA(max(2,isA(n)):ieA(n),nv,:),1);
237 subplot(100*ngEx+10+nv+2); ttmn=' Mx-mn:'; ttav=' Av:';
238 if nv == 1, titv2=' \theta'; else titv2=' S '; end
239 if list_log(ng) > 0,
240 var=zeros(nrec,Nexp); for n=1:Nexp, var(isA(n):ieA(n),n)=vvA(isA(n):ieA(n),nv,n); end
241 if min(var(:))*max(var(:)) >= 0, msk=var; ttyax='log';
242 var(find(msk==0))=1; var=log10(abs(var)); var(find(msk==0))=NaN;
243 end
244 else var=squeeze(vvA(:,nv,:)); end
245 for n=1:Nexp,
246 plot(ttA(isA(n):ieA(n),n),var(isA(n):ieA(n),n),char(linA(n)));
247 if n == 1, hold on ; end ;
248 ttmn=sprintf([ttmn,' %2.1e ;'],dd(nv,n));
249 ttav=sprintf([ttav,' %3.2e ;'],av(nv,n));
250 end ; hold off ;
251 if ttax1 < ttax2, AA=axis; axis([ttax1 ttax2 AA(3:4)]); end;
252 grid ;
253 %title(['mean ',titv1,titv2]);
254 title(['mean ',titv1,titv2,' ; ',ttav]);
255 %--
256 titX=' '; for n=1:Nexp, titX=[titX,sprintf(' %4.3e,',vvM(n))]; end
257 AA=axis; TM=text(AA(2),AA(3)+0.1*(AA(4)-AA(3)),titX(1:end-1));
258 set(TM,'HorizontalAlignment','right');
259 %--
260 end ;
261 end ; xlabel(titT);
262 %--
263 axes('position',[.01,.01,.99,.99],'Visible','off');
264 T=text(0.5,0.97,titall);
265 set(T,'HorizontalAlignment','center','FontSize',12);
266 put_date;
267 %---
268 end
269
270 if flag == 3
271 %--
272 figure(ng); set(ng,'position',[fxb+100*ng fyb+40*ng [500 700]*fsc]);clf;
273 dd=squeeze(max(vvA)-min(vvA)); av=squeeze(mean(vvA));
274 if Nexp == 1, av=av'; dd=dd'; end
275 for nv=1:4,
276 subplot(410+nv); ttmn=' Mx-mn:'; ttav=' Av:'; ttyax='';
277 if list_log(ng) > 0,
278 var=zeros(nrec,Nexp); for n=1:Nexp, var(isA(n):ieA(n),n)=vvA(isA(n):ieA(n),nv,n); end
279 if min(var(:))*max(var(:)) >= 0, msk=var; ttyax='log';
280 var(find(msk==0))=1; var=log10(abs(var)); var(find(msk==0))=NaN;
281 end
282 else var=squeeze(vvA(:,nv,:)); end
283 for n=1:Nexp,
284 plot(ttA(isA(n):ieA(n),n),var(isA(n):ieA(n),n),char(linA(n)));
285 if n == 1, hold on ; end ;
286 ttmn=sprintf([ttmn,' %2.1e ;'],dd(nv,n));
287 ttav=sprintf([ttav,' %3.2e ;'],av(nv,n));
288 end ; hold off ;
289 if ttax1 < ttax2, AA=axis; axis([ttax1 ttax2 AA(3:4)]); end;
290 if ttay(ng,1) < ttay(ng,2), AA=axis; AA(3)=max(AA(3),ttay(ng,1));
291 AA(4)=min(AA(4),ttay(ng,2)); axis(AA); end;
292 grid ;
293 if nv == 1, title(['U ',titv,' ',ttav]);
294 legend(namLg,'Location','best'); end
295 if nv == 2, title(['V ',titv,' ',ttav]); end
296 if nv == 3, title(['W ',titv,' ',ttav]); end
297 if nv == 4, title(['W+h ',titv,' ',ttav]); end
298 if length(ttyax) > 0, ylabel(ttyax); end
299 end ; xlabel(titT);
300 %--
301 axes('position',[.01,.01,.99,.99],'Visible','off');
302 T=text(0.5,0.97,titall);
303 set(T,'HorizontalAlignment','center','FontSize',12);
304 put_date;
305 %---
306 end
307 %--
308
309 if flag == 4
310 %--
311 % tt3i=' (+) '; for n=1:Nexp, tt3i=sprintf([tt3i,' %1.1e ;'],vvA(1,3,n));
312 % vvA(:,3,n)=vvA(:,3,n)-vvA(1,3,n); end
313 % tt4i=' (+) '; for n=1:Nexp, tt4i=sprintf([tt4i,' %1.1e ;'],vvA(1,4,n));
314 % vvA(:,4,n)=vvA(:,4,n)-vvA(1,4,n); end
315 figure(ng); set(ng,'position',[fxb+100*ng fyb+40*ng [500 700]*fsc]);clf;
316 dd=squeeze(max(vvA)-min(vvA)); av=squeeze(mean(vvA));
317 if Nexp == 1, av=av'; dd=dd'; end
318 % for nv=1:6, subplot(610+nv);
319 for nv=1:4, subplot(410+nv);
320 % for nV=1:4, nv=nV; subplot(410+nv); nv=max(1,2*(nv-1));
321 ttmn=' Mx-mn:'; ttav=' Av:'; ttyax='';
322 if list_log(ng) > 0,
323 var=zeros(nrec,Nexp); for n=1:Nexp, var(isA(n):ieA(n),n)=vvA(isA(n):ieA(n),nv,n); end
324 if min(var(:))*max(var(:)) >= 0, msk=var; ttyax='log';
325 var(find(msk==0))=1; var=log10(abs(var)); var(find(msk==0))=NaN;
326 end
327 else var=squeeze(vvA(:,nv,:)); end
328 for n=1:Nexp,
329 plot(ttA(isA(n):ieA(n),n),var(isA(n):ieA(n),n),char(linA(n)));
330 if n == 1, hold on ; end ;
331 ttmn=sprintf([ttmn,' %2.1e ;'],dd(nv,n));
332 ttav=sprintf([ttav,' %3.2e ;'],av(nv,n));
333 end ; hold off ;
334 if ttax1 < ttax2, AA=axis; axis([ttax1 ttax2 AA(3:4)]); end;
335 grid ;
336 if nv == 1, title(['min R.',titv,' ',ttav]); end
337 if nv == 2, title(['Max R.',titv,' ',ttav]);
338 legend(namLg,'Location','best'); end
339 % if nv == 3, title(['av A.',titv,' ',ttmn]); end
340 % if nv == 4, title(['sd A.',titv,' ',ttav]); end
341 % if nv == 5, title(['av P.',titv,' ',ttmn]); end
342 % if nv == 6, title(['sd P.',titv,' ',ttav]); end
343 if nv == 3, title(['Std-Dev A.',titv,' ',ttav]); end
344 if nv == 4, title(['Std-Dev P.',titv,' ',ttav]); end
345 end ; xlabel(titT);
346 %--
347 axes('position',[.01,.01,.99,.99],'Visible','off');
348 T=text(0.5,0.97,titall);
349 set(T,'HorizontalAlignment','center','FontSize',12);
350 put_date;
351 %---
352 end
353 %--
354
355 %-------------------
356 end
357
358 %=========================================================

  ViewVC Help
Powered by ViewVC 1.1.22