/[MITgcm]/MITgcm_contrib/gael/matlab_class/gcmfaces_diags/diags_set_E.m
ViewVC logotype

Annotation of /MITgcm_contrib/gael/matlab_class/gcmfaces_diags/diags_set_E.m

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


Revision 1.4 - (hide annotations) (download)
Wed Feb 12 05:10:26 2014 UTC (11 years, 5 months ago) by gforget
Branch: MAIN
Changes since 1.3: +18 -13 lines
(fix issues reported by H. Zhang)

- diags_set_E.m : fix cumulative sum in case of multiple records.
- diags_set_F.m : if not llc90, then bypass parts that require basin mask.
- diags_set_F.m : fix typos.
- diags_set_MLD.m : replace 240 with nt.
- diags_set_SEAICE.m : replace 240 with nt.

1 gforget 1.1
2     %select kBudget:
3     if ~isempty(setDiagsParams);
4     kBudget=setDiagsParams{1};
5     else;
6     kBudget=1;
7     end;
8    
9     %override default file name:
10     %---------------------------
11     tmp1=setDiags;
12     if kBudget>1;
13 gforget 1.2 tmp1=sprintf('E%02i',kBudget);
14 gforget 1.1 end;
15     fileMat=['diags_set_' tmp1];
16    
17     if userStep==1;%diags to be computed
18     listDiags=['zm_area zm_vol_ocn zm_vol_tot zm_vol_ice'];
19     listDiags=[listDiags ' zm_heat_ocn zm_heat_tot zm_heat_ice zm_heat_ocn_diff'];
20     listDiags=[listDiags ' zm_salt_ocn zm_salt_tot zm_salt_ice zm_salt_ocn_diff'];
21     elseif userStep==2;%input files and variables
22     listFlds={ 'ETAN','SIheff','SIhsnow','THETA ','SALT ','PHIBOT'};
23     listFlds={listFlds{:},'SIatmFW ','oceFWflx','SItflux','TFLUX','SFLUX','oceSPflx','SRELAX'};
24     listFlds={listFlds{:},'oceQnet ','SIatmQnt','SIaaflux','SIsnPrcp','SIacSubl'};
25     listFlds={listFlds{:},'TRELAX','WTHMASS','WSLTMASS','oceSflux','oceQsw','oceSPtnd'};
26     if kBudget>1;
27     listFlds={listFlds{:},'ADVr_TH','DFrE_TH','DFrI_TH','ADVr_SLT','DFrE_SLT','DFrI_SLT','WVELMASS'};
28     end;
29     listFlds={listFlds{:},'UVELMASS','VVELMASS','AB_gT','AB_gS'};
30     listFlds={listFlds{:},'ADVx_TH ','ADVy_TH ','DFxE_TH ','DFyE_TH '};
31     listFlds={listFlds{:},'ADVx_SLT','ADVy_SLT','DFxE_SLT','DFyE_SLT'};
32     listFlds={listFlds{:},'ADVxHEFF','ADVyHEFF','DFxEHEFF','DFyEHEFF'};
33     listFlds={listFlds{:},'ADVxSNOW','ADVySNOW','DFxESNOW','DFyESNOW'};
34     listFldsNames=deblank(listFlds);
35     %
36     listFiles={'rate_budg2d_snap_set1','budg2d_hflux_set1','budg2d_zflux_set1','budg2d_zflux_set2'};
37     if kBudget==1;
38     listFiles={listFiles{:},'rate_budg2d_snap_set2','budg2d_hflux_set2'};
39     else;
40     tmp1=sprintf('rate_budg2d_snap_set3_%02i',kBudget);
41     tmp2=sprintf('budg2d_zflux_set3_%02i',kBudget);
42     tmp3=sprintf('budg2d_hflux_set3_%02i',kBudget);
43     listFiles={listFiles{:},tmp1,tmp2,tmp3};
44     end;
45 gforget 1.2 listSubdirs={[dirMat 'BUDG/' ],[dirMat '../BUDG/' ],[dirModel 'diags/BUDG/']};
46 gforget 1.1 elseif userStep==3;%computational part;
47    
48     %preliminary tests
49     test1=isempty(dir([dirModel 'diags/BUDG/budg2d_snap_set1*']));
50 gforget 1.2 test2=isempty(dir([dirMat 'BUDG/rate_budg2d_snap_set1*']))&...
51     isempty(dir([dirMat '../BUDG/rate_budg2d_snap_set1*']));
52 gforget 1.1
53 gforget 1.2 if (strcmp(setDiags,'E')&test1&test2);
54 gforget 1.1 fprintf('\n abort : global and regional budgets, due to missing \n');
55     fprintf(['\n ' dirModel 'diags/BUDG/budg2d_snap_set1* \n']);
56     return;
57     end;
58    
59 gforget 1.2 if (strcmp(setDiags,'E')&test2);
60 gforget 1.1 fprintf('\n abort : global and regional budgets, due to missing \n');
61     fprintf(['\n ' dirModel 'diags/BUDG/rate_budg2d_snap_set1* \n']);
62     return;
63     end;
64    
65     %override default file name:
66     %---------------------------
67     tmp1=setDiags;
68     if kBudget>1;
69     tmp1=sprintf('E%02i',kBudget);
70     end;
71     fileMat=['diags_set_' tmp1 '_' num2str(tt) '.mat'];
72    
73     %fill in optional fields:
74     %------------------------
75     if isempty(who('TRELAX')); TRELAX=0*mygrid.XC; end;
76     if isempty(who('SRELAX')); SRELAX=0*mygrid.XC; end;
77     if isempty(who('AB_gT')); AB_gT=0*mygrid.XC; end;
78     if isempty(who('AB_gS')); AB_gS=0*mygrid.XC; end;
79     if isempty(who('oceSPtnd')); oceSPtnd=0*mygrid.XC; end;
80     if isempty(who('oceSPflx')); oceSPflx=0*mygrid.XC; end;
81     if isempty(who('PHIBOT')); PHIBOT=0*mygrid.XC; end;
82    
83     %=======MASS=========
84    
85     %compute mapped budget:
86     %----------------------
87    
88     %mass = myparms.rhoconst * sea level
89     contOCN=ETAN*myparms.rhoconst;
90     contICE=(SIheff*myparms.rhoi+SIhsnow*myparms.rhosn);
91     %for deep ocean layer :
92     if kBudget>1&myparms.useNLFS<2;
93     contOCN=0;
94     elseif kBudget>1;%rstar case
95     tmp1=mk3D(mygrid.DRF,mygrid.hFacC).*mygrid.hFacC;
96     tmp2=sum(tmp1(:,:,kBudget:length(mygrid.RC)),3)./mygrid.Depth;
97     contOCN=tmp2.*ETAN*myparms.rhoconst;
98     end;
99     %
100     contTOT=contOCN+contICE;
101     %vertical divergence (air-sea fluxes or vertical advection)
102     zdivOCN=oceFWflx;
103     zdivICE=SIatmFW-oceFWflx;
104     %in virtual salt flux we omit :
105     if ~myparms.useRFWF; zdivOCN=0*zdivOCN; end;
106    
107     %for deep ocean layer :
108     if kBudget>1; zdivOCN=-WVELMASS*myparms.rhoconst; end;
109     %
110     zdivTOT=zdivOCN+zdivICE;
111     %horizontal divergence (advection and ice diffusion)
112 gforget 1.3 hdivOCN=myparms.rhoconst*calc_UV_conv(UVELMASS,VVELMASS,{'dh'}); %for 2D, already vertically integrated, fields
113 gforget 1.1 tmpU=(myparms.rhoi*DFxEHEFF+myparms.rhosn*DFxESNOW+myparms.rhoi*ADVxHEFF+myparms.rhosn*ADVxSNOW);
114     tmpV=(myparms.rhoi*DFyEHEFF+myparms.rhosn*DFyESNOW+myparms.rhoi*ADVyHEFF+myparms.rhosn*ADVySNOW);
115 gforget 1.3 hdivICE=calc_UV_conv(tmpU,tmpV); %no dh needed here
116 gforget 1.1 hdivTOT=hdivOCN+hdivICE;
117     %bottom pressure for comparison:
118     bp=myparms.rhoconst/9.81*PHIBOT;
119    
120     [zm_vol_tot,zm_area]=calc_budget_mean_zonal(contTOT,zdivTOT,hdivTOT);
121     zm_vol_ocn=calc_budget_mean_zonal(contOCN,zdivOCN,hdivOCN);
122     zm_vol_ice=calc_budget_mean_zonal(contICE,zdivICE,hdivICE);
123    
124     %=======HEAT=======
125    
126     contOCN=myparms.rcp*THETA-myparms.rcp*AB_gT;
127     contICE=-myparms.flami*(SIheff*myparms.rhoi+SIhsnow*myparms.rhosn);
128     contTOT=contOCN+contICE;
129     %vertical divergence (air-sea fluxes or vertical adv/dif)
130     zdivOCN=TFLUX;
131     zdivICE=-(SItflux+TFLUX-TRELAX);
132     %in linear surface we omit :
133     if ~myparms.useNLFS; zdivOCN=zdivOCN-myparms.rcp*WTHMASS; end;
134     %in virtual salt flux we omit :
135     if ~myparms.useRFWF|~myparms.useNLFS; zdivICE=zdivICE+SIaaflux; end;
136     %working approach for real fresh water (?) and virtual salt flux
137     if 0; zdivICE=-oceQnet-SIatmQnt-myparms.flami*(SIsnPrcp-SIacSubl); end;
138     %for deep ocean layer :
139     if kBudget>1;
140     zdivOCN=-(ADVr_TH+DFrE_TH+DFrI_TH)./mygrid.RAC*myparms.rcp;
141     dd=mygrid.RF(kBudget); msk=mygrid.mskC(:,:,kBudget);
142     swfrac=0.62*exp(dd/0.6)+(1-0.62)*exp(dd/20);
143     if dd<-200; swfrac=0; end;
144     zdivOCN=zdivOCN+swfrac*oceQsw;%.*msk;
145     end;
146     %
147     zdivTOT=zdivOCN+zdivICE;
148     %horizontal divergence (advection and diffusion)
149     tmpU=myparms.rcp*(ADVx_TH+DFxE_TH); tmpV=myparms.rcp*(ADVy_TH+DFyE_TH);
150 gforget 1.3 hdivOCN=calc_UV_conv(tmpU,tmpV);
151 gforget 1.1 tmpU=-myparms.flami*(myparms.rhoi*DFxEHEFF+myparms.rhosn*DFxESNOW+myparms.rhoi*ADVxHEFF+myparms.rhosn*ADVxSNOW);
152     tmpV=-myparms.flami*(myparms.rhoi*DFyEHEFF+myparms.rhosn*DFyESNOW+myparms.rhoi*ADVyHEFF+myparms.rhosn*ADVySNOW);
153 gforget 1.3 hdivICE=calc_UV_conv(tmpU,tmpV); %no dh needed here
154 gforget 1.1 hdivTOT=hdivOCN+hdivICE;
155    
156     zm_heat_tot=calc_budget_mean_zonal(contTOT,zdivTOT,hdivTOT);
157     zm_heat_ocn=calc_budget_mean_zonal(contOCN,zdivOCN,hdivOCN);
158     zm_heat_ice=calc_budget_mean_zonal(contICE,zdivICE,hdivICE);
159    
160     %ocean diffusion alone
161     tmpU=myparms.rcp*(DFxE_TH); tmpV=myparms.rcp*(DFyE_TH);
162 gforget 1.3 hdivOCN=calc_UV_conv(tmpU,tmpV);
163 gforget 1.1 zm_heat_ocn_diff=calc_budget_mean_zonal(0*contOCN,0*zdivOCN,hdivOCN);
164    
165     %=======SALT=======
166    
167     contOCN=myparms.rhoconst*SALT-myparms.rhoconst*AB_gS;
168     contICE=myparms.SIsal0*myparms.rhoi*SIheff;
169     contTOT=contOCN+contICE;
170     %vertical divergence (air-sea fluxes or vertical adv/dif)
171     zdivOCN=SFLUX+oceSPflx;
172     zdivICE=-zdivOCN+SRELAX;
173     %in linear surface we omit :
174     if ~myparms.useNLFS; zdivOCN=zdivOCN-myparms.rhoconst*WSLTMASS; end;
175     %working approach for real fresh water (?) and virtual salt flux
176     if ~myparms.useRFWF|~myparms.useNLFS; zdivICE=-oceSflux; end;
177     %for deep ocean layer :
178     if kBudget>1;
179     zdivOCN=-(ADVr_SLT+DFrE_SLT+DFrI_SLT)./mygrid.RAC*myparms.rhoconst;
180     zdivOCN=zdivOCN+oceSPtnd;%.*msk;
181     end;
182     zdivTOT=zdivOCN+zdivICE;
183     %horizontal divergence (advection and diffusion)
184     tmpU=myparms.rhoconst*(ADVx_SLT+DFxE_SLT); tmpV=myparms.rhoconst*(ADVy_SLT+DFyE_SLT);
185 gforget 1.3 hdivOCN=calc_UV_conv(tmpU,tmpV);
186 gforget 1.1 tmpU=myparms.SIsal0*(myparms.rhoi*DFxEHEFF+myparms.rhoi*ADVxHEFF);
187     tmpV=myparms.SIsal0*(myparms.rhoi*DFyEHEFF+myparms.rhoi*ADVyHEFF);
188 gforget 1.3 hdivICE=calc_UV_conv(tmpU,tmpV); %no dh needed here
189 gforget 1.1 hdivTOT=hdivOCN+hdivICE;
190    
191     zm_salt_tot=calc_budget_mean_zonal(contTOT,zdivTOT,hdivTOT);
192     zm_salt_ocn=calc_budget_mean_zonal(contOCN,zdivOCN,hdivOCN);
193     zm_salt_ice=calc_budget_mean_zonal(contICE,zdivICE,hdivICE);
194    
195     %ocean diffusion alone
196     tmpU=myparms.rhoconst*(DFxE_SLT); tmpV=myparms.rhoconst*(DFyE_SLT);
197 gforget 1.3 hdivOCN=calc_UV_conv(tmpU,tmpV);
198 gforget 1.1 zm_salt_ocn_diff=calc_budget_mean_zonal(0*contOCN,0*zdivOCN,hdivOCN);
199    
200     %===================== COMPUTATIONAL SEQUENCE ENDS =========================%
201     %===================== PLOTTING SEQUENCE BEGINS =========================%
202    
203     elseif userStep==-1;%plotting
204    
205     if isempty(setDiagsParams);
206     choicePlot={'all'};
207     elseif isnumeric(setDiagsParams{1})&length(setDiagsParams)==1;
208     choicePlot={'all'};
209     elseif isnumeric(setDiagsParams{1});
210     choicePlot={setDiagsParams{2:end}};
211     else;
212     choicePlot=setDiagsParams;
213     end;
214    
215     tt=[1:length(alldiag.listTimes)];
216     TT=alldiag.listTimes(tt);
217     nt=length(TT);
218    
219     if (kBudget==1)&(sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'mass')));
220    
221     %1.1) ocean+seaice mass budgets
222     %------------------------------
223     figureL;
224     %volume budget:
225     subplot(2,1,1); disp_budget_mean_zonal(mygrid.LATS,alldiag.zm_vol_tot,'kg/m2','Mass (incl. ice)');
226     %cumulative integral:
227 gforget 1.4 tmp1=repmat(alldiag.zm_area,[3 1 1]);
228     tmp1=tmp1.*alldiag.zm_vol_tot; cumbudg=cumsum(tmp1,2);
229 gforget 1.1 subplot(2,1,2); disp_budget_mean_zonal(mygrid.LATS,cumbudg,'kg','Mass (incl. ice)');
230     %add to tex file
231     myCaption={myYmeanTxt,'mass budget (ocean+ice) at each latitude in kg/m2 (upper) and integrated from South (lower).'};
232     if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); end;
233    
234     %1.2) ice mass budgets
235     %---------------------
236     figureL;
237     %volume budget:
238     subplot(2,1,1); disp_budget_mean_zonal(mygrid.LATS,alldiag.zm_vol_ice,'kg/m2','Mass (only ice)');
239     %cumulative integral:
240 gforget 1.4 tmp1=repmat(alldiag.zm_area,[3 1 1]);
241     tmp1=tmp1.*alldiag.zm_vol_ice; cumbudg=cumsum(tmp1,2);
242 gforget 1.1 subplot(2,1,2); disp_budget_mean_zonal(mygrid.LATS,cumbudg,'kg','Mass (only ice)');
243     %add to tex file
244     myCaption={myYmeanTxt,'mass budget (only ice) at each latitude in kg/m2 (upper) and integrated from South (lower).'};
245     if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); end;
246     end;
247    
248     if (sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'mass')));
249    
250     %1.3) ocean mass budgets
251     %-----------------------
252     figureL;
253     %volume budget:
254     subplot(2,1,1); disp_budget_mean_zonal(mygrid.LATS,alldiag.zm_vol_ocn,'kg/m2','Mass (ocean only)');
255     %cumulative integral:
256 gforget 1.4 tmp1=repmat(alldiag.zm_area,[3 1 1]);
257     tmp1=tmp1.*alldiag.zm_vol_ocn; cumbudg=cumsum(tmp1,2);
258 gforget 1.1 subplot(2,1,2); disp_budget_mean_zonal(mygrid.LATS,cumbudg,'kg','Mass (ocean only)');
259     %add to tex file
260     myCaption={myYmeanTxt,'mass budget (ocean only) at each latitude in kg/m2 (upper) and integrated from South (lower).'};
261     if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); end;
262    
263     end;
264    
265     if (kBudget==1)&(sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'heat')));
266    
267     %1.1) ocean+seaice heat budgets
268     %------------------------------
269     figureL;
270     %heat budget:
271     subplot(2,1,1); disp_budget_mean_zonal(mygrid.LATS,alldiag.zm_heat_tot,'J/m2','Heat (incl. ice)');
272     %cumulative integral:
273 gforget 1.4 tmp1=repmat(alldiag.zm_area,[3 1 1]);
274     tmp1=tmp1.*alldiag.zm_heat_tot; cumbudg=cumsum(tmp1,2);
275 gforget 1.1 subplot(2,1,2); disp_budget_mean_zonal(mygrid.LATS,cumbudg,'J','Heat (incl. ice)');
276     %add to tex file
277     myCaption={myYmeanTxt,'heat budget (ocean+ice) at each latitude in J/m2 (upper) and integrated from South (lower).'};
278     if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); end;
279    
280     %1.2) ice heat budgets
281     %---------------------
282     figureL;
283     %heat budget:
284     subplot(2,1,1); disp_budget_mean_zonal(mygrid.LATS,alldiag.zm_heat_ice,'J/m2','Heat (only ice)');
285     %cumulative integral:
286 gforget 1.4 tmp1=repmat(alldiag.zm_area,[3 1 1]);
287     tmp1=tmp1.*alldiag.zm_heat_ice; cumbudg=cumsum(tmp1,2);
288 gforget 1.1 subplot(2,1,2); disp_budget_mean_zonal(mygrid.LATS,cumbudg,'J','Heat (only ice)');
289     %add to tex file
290     myCaption={myYmeanTxt,'heat budget (only ice) at each latitude in J/m2 (upper) and integrated from South (lower).'};
291     if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); end;
292     end;
293    
294     if (sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'heat')));
295    
296     %1.3) ocean heat budgets
297     %-----------------------
298     figureL;
299     %heat budget:
300     subplot(2,1,1); disp_budget_mean_zonal(mygrid.LATS,alldiag.zm_heat_ocn,'J/m2','Heat (ocean only)');
301     %cumulative integral:
302 gforget 1.4 tmp1=repmat(alldiag.zm_area,[3 1 1]);
303     tmp1=tmp1.*alldiag.zm_heat_ocn; cumbudg=cumsum(tmp1,2);
304 gforget 1.1 subplot(2,1,2); disp_budget_mean_zonal(mygrid.LATS,cumbudg,'J','Heat (ocean only)');
305     %add to tex file
306     myCaption={myYmeanTxt,'heat budget (ocean only) at each latitude in J/m2 (upper) and integrated from South (lower).'};
307     if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); end;
308    
309     end;
310    
311     if (kBudget==1)&(sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'salt')));
312    
313     %1.1) ocean+seaice salt budgets
314     %------------------------------
315     figureL;
316     %salt budget:
317     subplot(2,1,1); disp_budget_mean_zonal(mygrid.LATS,alldiag.zm_salt_tot,'g/m2','Salt (incl. ice)');
318     %cumulative integral:
319 gforget 1.4 tmp1=repmat(alldiag.zm_area,[3 1 1]);
320     tmp1=tmp1.*alldiag.zm_salt_tot; cumbudg=cumsum(tmp1,2);
321 gforget 1.1 subplot(2,1,2); disp_budget_mean_zonal(mygrid.LATS,cumbudg,'g','Salt (incl. ice)');
322     %add to tex file
323     myCaption={myYmeanTxt,'salt budget (ocean+ice) at each latitude in g/m2 (upper) and integrated from South (lower).'};
324     if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); end;
325    
326     %1.2) ice salt budgets
327     %---------------------
328     figureL;
329     %salt budget:
330     subplot(2,1,1); disp_budget_mean_zonal(mygrid.LATS,alldiag.zm_salt_ice,'g/m2','Salt (only ice)');
331     %cumulative integral:
332 gforget 1.4 tmp1=repmat(alldiag.zm_area,[3 1 1]);
333     tmp1=tmp1.*alldiag.zm_salt_ice; cumbudg=cumsum(tmp1,2);
334 gforget 1.1 subplot(2,1,2); disp_budget_mean_zonal(mygrid.LATS,cumbudg,'g','Salt (only ice)');
335     %add to tex file
336     myCaption={myYmeanTxt,'salt budget (only ice) at each latitude in g/m2 (upper) and integrated from South (lower).'};
337     if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); end;
338     end;
339    
340     if (sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'salt')));
341    
342     %1.3) ocean salt budgets
343     %-----------------------
344     figureL;
345     %salt budget:
346     subplot(2,1,1); disp_budget_mean_zonal(mygrid.LATS,alldiag.zm_salt_ocn,'g/m2','Salt (ocean only)');
347     %cumulative integral:
348 gforget 1.4 tmp1=repmat(alldiag.zm_area,[3 1 1]);
349     tmp1=tmp1.*alldiag.zm_salt_ocn; cumbudg=cumsum(tmp1,2);
350 gforget 1.1 subplot(2,1,2); disp_budget_mean_zonal(mygrid.LATS,cumbudg,'g','Salt (ocean only)');
351     %add to tex file
352     myCaption={myYmeanTxt,'salt budget (ocean only) at each latitude in g/m2 (upper) and integrated from South (lower).'};
353     if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); end;
354    
355     end;
356    
357     end;
358    

  ViewVC Help
Powered by ViewVC 1.1.22