/[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.1 - (hide annotations) (download)
Tue Jan 8 17:04:15 2013 UTC (12 years, 6 months ago) by gforget
Branch: MAIN
- diags_set_LAYERS.m : treat basins as done in z-coordinate case, and add display section.
- diags_set_E.m : zonal mean budgets.

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

  ViewVC Help
Powered by ViewVC 1.1.22