/[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.2 - (hide annotations) (download)
Tue Feb 8 17:06:00 2022 UTC (3 years, 5 months ago) by jmc
Branch: MAIN
Changes since 1.1: +13 -6 lines
fix "mean surf-Cor" plots

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

  ViewVC Help
Powered by ViewVC 1.1.22