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

Annotation of /MITgcm_contrib/jmc_script/grph_MON.m

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


Revision 1.4 - (hide 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 jmc 1.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 jmc 1.4 % $Header: /u/gcmpack/MITgcm_contrib/jmc_script/grph_MON.m,v 1.3 2022/02/08 17:10:21 jmc Exp $
14 jmc 1.2 % $Name: $
15 jmc 1.1
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 jmc 1.4 %- in case PE is missing from Monitor, use Mean-Eta and this "gdH" ratio to get PE
22 jmc 1.1 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 jmc 1.4 % 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 jmc 1.1 %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 jmc 1.4 if ddKe > 0, keA(:,3,n)=keA(:,4,n); %- take PE from Monitor (since it's there)
72 jmc 1.1 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 jmc 1.4 %-- 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 jmc 1.1 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 jmc 1.4 if ng == 6, flag=1*list_on(7) ; titv='W-vel [m/s]' ; vvA=wmA ; end
126 jmc 1.1 %if ng == 5, flag=1*list_on(2) ; titv='Eta [mb]'; vvA=etA/100 ; end
127 jmc 1.4 %if ng == 6, flag=1*list_on(7) ; titv='W-vel [Pa/s]'; vvA=wmA ; end
128 jmc 1.1 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 jmc 1.4 % LL(n)=plot(ttA(isA(n):ieA(n),n),var(isA(n):ieA(n),n),char(linA(n)));
207 jmc 1.1 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 jmc 1.4 % set(LL(1),'LineWidth',2);
212 jmc 1.1 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 jmc 1.3 titv1='surf-Cor : ' ; vvA=scA ;
233 jmc 1.1 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 jmc 1.2 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 jmc 1.1 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 jmc 1.2 %title(['mean ',titv1,titv2]);
254     title(['mean ',titv1,titv2,' ; ',ttav]);
255 jmc 1.1 %--
256 jmc 1.2 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 jmc 1.1 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