/[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.3 - (show annotations) (download)
Tue Feb 8 17:10:21 2022 UTC (3 years, 5 months ago) by jmc
Branch: MAIN
Changes since 1.2: +2 -2 lines
fix previous modif

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.2 2022/02/08 17:06:00 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 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 titv1='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 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 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 %title(['mean ',titv1,titv2]);
246 title(['mean ',titv1,titv2,' ; ',ttav]);
247 %--
248 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 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