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

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

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


Revision 1.12 - (hide annotations) (download)
Tue Feb 3 00:53:23 2015 UTC (10 years, 5 months ago) by gforget
Branch: MAIN
Changes since 1.11: +62 -10 lines
- update optional output of budget in extensive form (kg/s,J/s,g/s)
  so that it works in three-dimensional case (test3d).

1 gforget 1.1
2     %select kBudget:
3     if ~isempty(setDiagsParams);
4     kBudget=setDiagsParams{1};
5     else;
6     kBudget=1;
7     end;
8    
9 gforget 1.9 doMoreBudgetOutput=0;
10    
11 gforget 1.2 %override default file name:
12     %---------------------------
13     tmp1=setDiags;
14     if kBudget>1;
15     tmp1=sprintf('D%02i',kBudget);
16 gforget 1.1 end;
17 gforget 1.2 fileMat=['diags_set_' tmp1];
18 gforget 1.1
19     if userStep==1;%diags to be computed
20     listDiags=['glo_vol_ocn glo_vol_tot glo_vol_ice glo_bp'];
21     listDiags=[listDiags ' north_vol_ocn north_vol_tot north_vol_ice north_bp'];
22     listDiags=[listDiags ' south_vol_ocn south_vol_tot south_vol_ice south_bp'];
23     listDiags=[listDiags ' glo_heat_ocn glo_heat_tot glo_heat_ice'];
24     listDiags=[listDiags ' north_heat_ocn north_heat_tot north_heat_ice'];
25     listDiags=[listDiags ' south_heat_ocn south_heat_tot south_heat_ice'];
26     listDiags=[listDiags ' glo_salt_ocn glo_salt_tot glo_salt_ice'];
27     listDiags=[listDiags ' north_salt_ocn north_salt_tot north_salt_ice'];
28     listDiags=[listDiags ' south_salt_ocn south_salt_tot south_salt_ice'];
29 gforget 1.2
30 gforget 1.1 elseif userStep==2;%input files and variables
31 gforget 1.11 tmp1=fullfile(dirModel,'diags',filesep,'BUDG',filesep,'budg2d_snap_set2*meta');
32     test3d=isempty(dir(tmp1));
33     %
34 gforget 1.1 listFlds={ 'ETAN','SIheff','SIhsnow','THETA ','SALT ','PHIBOT'};
35     listFlds={listFlds{:},'SIatmFW ','oceFWflx','SItflux','TFLUX','SFLUX','oceSPflx','SRELAX'};
36     listFlds={listFlds{:},'oceQnet ','SIatmQnt','SIaaflux','SIsnPrcp','SIacSubl'};
37     listFlds={listFlds{:},'TRELAX','WTHMASS','WSLTMASS','oceSflux','oceQsw','oceSPtnd'};
38 gforget 1.11 if kBudget>1|test3d;
39 gforget 1.1 listFlds={listFlds{:},'ADVr_TH','DFrE_TH','DFrI_TH','ADVr_SLT','DFrE_SLT','DFrI_SLT','WVELMASS'};
40     end;
41 gforget 1.4 listFlds={listFlds{:},'SDIAG1','SDIAG2','SDIAG3'};
42 gforget 1.1 listFlds={listFlds{:},'UVELMASS','VVELMASS','AB_gT','AB_gS'};
43     listFlds={listFlds{:},'ADVx_TH ','ADVy_TH ','DFxE_TH ','DFyE_TH '};
44     listFlds={listFlds{:},'ADVx_SLT','ADVy_SLT','DFxE_SLT','DFyE_SLT'};
45     listFlds={listFlds{:},'ADVxHEFF','ADVyHEFF','DFxEHEFF','DFyEHEFF'};
46     listFlds={listFlds{:},'ADVxSNOW','ADVySNOW','DFxESNOW','DFyESNOW'};
47     listFldsNames=deblank(listFlds);
48     %
49     listFiles={'rate_budg2d_snap_set1','budg2d_hflux_set1','budg2d_zflux_set1','budg2d_zflux_set2'};
50 gforget 1.11 if test3d;
51     listFiles={listFiles{:},'rate_budg3d_snap_set1','budg3d_hflux_set1','budg3d_zflux_set1'};
52     elseif kBudget==1;
53 gforget 1.1 listFiles={listFiles{:},'rate_budg2d_snap_set2','budg2d_hflux_set2'};
54     else;
55     tmp1=sprintf('rate_budg2d_snap_set3_%02i',kBudget);
56     tmp2=sprintf('budg2d_zflux_set3_%02i',kBudget);
57     tmp3=sprintf('budg2d_hflux_set3_%02i',kBudget);
58     listFiles={listFiles{:},tmp1,tmp2,tmp3};
59     end;
60 gforget 1.5 listSubdirs={[dirMat 'BUDG/' ],[dirMat '../BUDG/' ],[dirModel 'diags/BUDG/']};
61 gforget 1.2
62 gforget 1.1 elseif userStep==3;%computational part;
63 gforget 1.2
64     %preliminary tests
65     test1=isempty(dir([dirModel 'diags/BUDG/budg2d_snap_set1*']));
66 gforget 1.5 test2=isempty(dir([dirMat 'BUDG/rate_budg2d_snap_set1*']))&...
67     isempty(dir([dirMat '../BUDG/rate_budg2d_snap_set1*']));
68 gforget 1.2 if (strcmp(setDiags,'D')&test1&test2);
69     fprintf('\n abort : global and regional budgets, due to missing \n');
70     fprintf(['\n ' dirModel 'diags/BUDG/budg2d_snap_set1* \n']);
71     return;
72     end;
73    
74     if (strcmp(setDiags,'D')&test2);
75     fprintf('\n abort : global and regional budgets, due to missing \n');
76     fprintf(['\n ' dirModel 'diags/BUDG/rate_budg2d_snap_set1* \n']);
77     return;
78     end;
79 gforget 1.1
80     %override default file name:
81     %---------------------------
82     tmp1=setDiags;
83     if kBudget>1;
84     tmp1=sprintf('D%02i',kBudget);
85     end;
86     fileMat=['diags_set_' tmp1 '_' num2str(tt) '.mat'];
87    
88     %fill in optional fields:
89     %------------------------
90 gforget 1.11 if isempty(who('TRELAX')); TRELAX=0; end;
91     if isempty(who('SRELAX')); SRELAX=0; end;
92     if isempty(who('AB_gT')); AB_gT=0; end;
93     if isempty(who('AB_gS')); AB_gS=0; end;
94     if isempty(who('oceSPtnd')); oceSPtnd=0; end;
95     if isempty(who('oceSPflx')); oceSPflx=0; end;
96     if isempty(who('PHIBOT')); PHIBOT=0; end;
97 gforget 1.4
98     %aliases from development phase (applies to 2012 core runs)
99     %---------------------------------------------------------
100     if ~isempty(who('SDIAG1')); SRELAX=SDIAG1; end;
101     if ~isempty(who('SDIAG2')); SIatmFW=SDIAG2; end;
102     if ~isempty(who('SDIAG3')); SItflux=SDIAG3; end;
103    
104 gforget 1.1
105     %=======MASS=========
106    
107 gforget 1.9 if doMoreBudgetOutput;
108     %indexing and sign convention:
109     %- MITgcm: fluxes are >0 downward, k=1 start at free surface
110     %- here, similarly: >0 downward, k=1 free surface k=2 sea floor
111 gforget 1.12 if ~test3d;
112     budgO.specs.top='free surface';
113     if kBudget>1; budgO.specs.top=['interface no. ' num2str(kBudget)]; end;
114     budgO.specs.bottom='sea floor';
115     else;
116     budgO.specs.top='interface k';
117     budgO.specs.bottom='interface k+1';
118     end;
119 gforget 1.9 budgI.specs.top='ocn-ice to atm interface';
120     budgI.specs.bottom='free surface';
121    
122     %here we output tendencies and fluxes in kg/s
123     budgMo=budgO; budgMi=budgI;
124     budgMo.specs.units='kg/s';%ocean only
125     budgMi.specs.units='kg/s';%ice only
126     %here we output tendencies and fluxes in Watts
127     budgHo=budgO; budgHi=budgI;
128     budgHo.specs.units='W';%ocean only
129     budgHi.specs.units='W';%ice only
130     %here we output tendencies and fluxes in g/s
131     budgSo=budgO; budgSi=budgI;
132     budgSo.specs.units='g/s';%ocean only
133     budgSi.specs.units='g/s';%ice only
134     end;
135    
136 gforget 1.1 %compute mapped budget:
137     %----------------------
138    
139     %mass = myparms.rhoconst * sea level
140     contOCN=ETAN*myparms.rhoconst;
141     contICE=(SIheff*myparms.rhoi+SIhsnow*myparms.rhosn);
142     %for deep ocean layer :
143     if kBudget>1&myparms.useNLFS<2;
144     contOCN=0;
145     elseif kBudget>1;%rstar case
146     tmp1=mk3D(mygrid.DRF,mygrid.hFacC).*mygrid.hFacC;
147     tmp2=sum(tmp1(:,:,kBudget:length(mygrid.RC)),3)./mygrid.Depth;
148     contOCN=tmp2.*ETAN*myparms.rhoconst;
149     end;
150     %
151     contTOT=contOCN+contICE;
152 gforget 1.9 %
153     if doMoreBudgetOutput;
154 gforget 1.12 if test3d;
155     tmp1=mk3D(mygrid.DRF,mygrid.hFacC).*mygrid.hFacC;
156     tmp2=tmp1./mk3D(mygrid.Depth,tmp1);
157     tend=tmp2.*mk3D(ETAN,tmp2)*myparms.rhoconst;
158     else;
159     tend=contOCN;
160     end;
161     budgMo.tend=mk3D(mygrid.RAC,tend).*tend;%kg/s
162 gforget 1.9 budgMi.tend=mygrid.RAC.*contICE;%kg/s
163     end;
164    
165 gforget 1.1 %vertical divergence (air-sea fluxes or vertical advection)
166     zdivOCN=oceFWflx;
167     zdivICE=SIatmFW-oceFWflx;
168     %in virtual salt flux we omit :
169     if ~myparms.useRFWF; zdivOCN=0*zdivOCN; end;
170     %for deep ocean layer :
171     if kBudget>1; zdivOCN=-WVELMASS*myparms.rhoconst; end;
172     %
173     zdivTOT=zdivOCN+zdivICE;
174 gforget 1.9 %
175     if doMoreBudgetOutput;
176 gforget 1.12 if test3d;
177     trWtop=-WVELMASS*myparms.rhoconst;
178     %trWtop(:,:,1)=zdivOCN;
179     trWbot=trWtop(:,:,2:length(mygrid.RC));
180     trWbot(:,:,length(mygrid.RC))=0;
181     %
182     budgMo.trWtop=mk3D(mygrid.RAC,trWtop).*trWtop;
183     budgMo.trWbot=mk3D(mygrid.RAC,trWbot).*trWbot;%kg/s
184     else;
185 gforget 1.10 budgMo.trWtop=mygrid.RAC.*zdivOCN; budgMo.trWbot=mygrid.RAC*0;%kg/s
186 gforget 1.12 end;
187     budgMi.trWtop=mygrid.RAC.*(zdivICE+zdivOCN); budgMi.trWbot=mygrid.RAC.*zdivOCN;%kg/s
188 gforget 1.9 end;
189    
190 gforget 1.1 %horizontal divergence (advection and ice diffusion)
191 gforget 1.11 if test3d;
192     %3D UVELMASS,VVELMASS are multiplied by DRF
193     %(2D diagnostics are expectedly vertically integrated by MITgcm)
194     tmp1=mk3D(mygrid.DRF,UVELMASS);
195     UVELMASS=tmp1.*UVELMASS;
196     VVELMASS=tmp1.*VVELMASS;
197     end;
198 gforget 1.9 dxg=mk3D(mygrid.DXG,VVELMASS); dyg=mk3D(mygrid.DYG,UVELMASS);
199     tmpUo=myparms.rhoconst*dyg.*UVELMASS;
200     tmpVo=myparms.rhoconst*dxg.*VVELMASS;
201 gforget 1.11 hdivOCN=calc_UV_conv(nansum(tmpUo,3),nansum(tmpVo,3));
202 gforget 1.9 tmpUi=(myparms.rhoi*DFxEHEFF+myparms.rhosn*DFxESNOW+myparms.rhoi*ADVxHEFF+myparms.rhosn*ADVxSNOW);
203     tmpVi=(myparms.rhoi*DFyEHEFF+myparms.rhosn*DFyESNOW+myparms.rhoi*ADVyHEFF+myparms.rhosn*ADVySNOW);
204     hdivICE=calc_UV_conv(tmpUi,tmpVi); %dh needed is alerady in DFxEHEFF etc
205 gforget 1.1 hdivTOT=hdivOCN+hdivICE;
206 gforget 1.9 if doMoreBudgetOutput;
207     budgMo.trU=tmpUo; budgMo.trV=tmpVo;%kg/s
208     budgMi.trU=tmpUi; budgMi.trV=tmpVi;%kg/s
209     end;
210 gforget 1.1
211     %bottom pressure for comparison:
212     bp=myparms.rhoconst/9.81*PHIBOT;
213    
214     %compute global integrals:
215     %-------------------------
216     msk=mygrid.mskC(:,:,kBudget);
217     glo_vol_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
218     glo_vol_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
219     glo_vol_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
220     glo_bp=nansum(bp.*msk.*mygrid.RAC)/nansum(msk.*mygrid.RAC);
221    
222     %compute northern hemisphere integrals:
223     msk=mygrid.mskC(:,:,kBudget).*(mygrid.YC>0);
224     north_vol_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
225     north_vol_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
226     north_vol_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
227     north_bp=nansum(bp.*msk.*mygrid.RAC)/nansum(msk.*mygrid.RAC);
228    
229     %and southern hemisphere integrals:
230     msk=mygrid.mskC(:,:,kBudget).*(mygrid.YC<=0);
231     south_vol_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
232     south_vol_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
233     south_vol_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
234     south_bp=nansum(bp.*msk.*mygrid.RAC)/nansum(msk.*mygrid.RAC);
235    
236     %=======HEAT=======
237 gforget 1.11
238 gforget 1.1 contOCN=myparms.rcp*THETA-myparms.rcp*AB_gT;
239     contICE=-myparms.flami*(SIheff*myparms.rhoi+SIhsnow*myparms.rhosn);
240 gforget 1.9 %
241     if doMoreBudgetOutput;
242 gforget 1.11 budgHo.tend=mk3D(mygrid.RAC,contOCN).*contOCN;%Watt
243 gforget 1.9 budgHi.tend=mygrid.RAC.*contICE;%Watt
244     end;
245 gforget 1.11 contOCN=nansum(contOCN,3);
246     contTOT=contOCN+contICE;
247 gforget 1.9
248 gforget 1.1 %vertical divergence (air-sea fluxes or vertical adv/dif)
249     zdivOCN=TFLUX;
250     zdivICE=-(SItflux+TFLUX-TRELAX);
251     %in linear surface we omit :
252     if ~myparms.useNLFS; zdivOCN=zdivOCN-myparms.rcp*WTHMASS; end;
253     %in virtual salt flux we omit :
254     if ~myparms.useRFWF|~myparms.useNLFS; zdivICE=zdivICE+SIaaflux; end;
255     %working approach for real fresh water (?) and virtual salt flux
256     if 0; zdivICE=-oceQnet-SIatmQnt-myparms.flami*(SIsnPrcp-SIacSubl); end;
257     %for deep ocean layer :
258     if kBudget>1;
259     zdivOCN=-(ADVr_TH+DFrE_TH+DFrI_TH)./mygrid.RAC*myparms.rcp;
260     dd=mygrid.RF(kBudget); msk=mygrid.mskC(:,:,kBudget);
261     swfrac=0.62*exp(dd/0.6)+(1-0.62)*exp(dd/20);
262     if dd<-200; swfrac=0; end;
263     zdivOCN=zdivOCN+swfrac*oceQsw;%.*msk;
264     end;
265     %
266     zdivTOT=zdivOCN+zdivICE;
267 gforget 1.9 %
268     if doMoreBudgetOutput;
269 gforget 1.12 if test3d;
270     trWtop=-(ADVr_TH+DFrE_TH+DFrI_TH)*myparms.rcp;
271     %
272     dd=mygrid.RF(1:end-1);
273     swfrac=0.62*exp(dd/0.6)+(1-0.62)*exp(dd/20);
274     swfrac(dd<-200)=0;
275     swtop=mk3D(swfrac,trWtop).*mk3D(oceQsw,trWtop);
276     swtop(isnan(mygrid.mskC))=0;
277     trWtop=trWtop+swtop;
278     %
279     trWtop(:,:,1)=zdivOCN;
280     trWbot=trWtop(:,:,2:length(mygrid.RC));
281     trWbot(:,:,length(mygrid.RC))=0;
282     %
283     budgHo.trWtop=trWtop;%Watt
284     budgHo.trWbot=trWbot;%Watt
285     else;
286     budgHo.trWtop=mygrid.RAC.*zdivOCN; budgHo.trWbot=mygrid.RAC*0;%Watt
287     end;
288     budgHi.trWtop=mygrid.RAC.*(zdivICE+zdivOCN); budgHi.trWbot=mygrid.RAC.*zdivOCN;%Watt
289 gforget 1.9 end;
290 gforget 1.11
291 gforget 1.1 %horizontal divergence (advection and diffusion)
292 gforget 1.9 tmpUo=myparms.rcp*(ADVx_TH+DFxE_TH); tmpVo=myparms.rcp*(ADVy_TH+DFyE_TH);
293 gforget 1.11 hdivOCN=calc_UV_conv(nansum(tmpUo,3),nansum(tmpVo,3));
294 gforget 1.9 tmpUi=-myparms.flami*(myparms.rhoi*DFxEHEFF+myparms.rhosn*DFxESNOW+myparms.rhoi*ADVxHEFF+myparms.rhosn*ADVxSNOW);
295     tmpVi=-myparms.flami*(myparms.rhoi*DFyEHEFF+myparms.rhosn*DFyESNOW+myparms.rhoi*ADVyHEFF+myparms.rhosn*ADVySNOW);
296     hdivICE=calc_UV_conv(tmpUi,tmpVi); %no dh needed here
297 gforget 1.1 hdivTOT=hdivOCN+hdivICE;
298 gforget 1.9 if doMoreBudgetOutput;
299     budgHo.trU=tmpUo; budgHo.trV=tmpVo;%Watt
300     budgHi.trU=tmpUi; budgHi.trV=tmpVi;%Watt
301     end;
302 gforget 1.1
303     %compute global integrals:
304     %-------------------------
305     msk=mygrid.mskC(:,:,kBudget);
306     glo_heat_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
307     glo_heat_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
308     glo_heat_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
309    
310     %compute northern hemisphere integrals:
311     msk=mygrid.mskC(:,:,kBudget).*(mygrid.YC>0);
312     north_heat_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
313     north_heat_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
314     north_heat_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
315    
316     %and southern hemisphere integrals:
317     msk=mygrid.mskC(:,:,kBudget).*(mygrid.YC<=0);
318     south_heat_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
319     south_heat_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
320     south_heat_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
321 gforget 1.11
322 gforget 1.1 %=======SALT=======
323 gforget 1.9
324 gforget 1.1 contOCN=myparms.rhoconst*SALT-myparms.rhoconst*AB_gS;
325     contICE=myparms.SIsal0*myparms.rhoi*SIheff;
326 gforget 1.9 %
327     if doMoreBudgetOutput;
328 gforget 1.12 budgSo.tend=mk3D(mygrid.RAC,contOCN).*contOCN;%g/s
329 gforget 1.9 budgSi.tend=mygrid.RAC.*contICE;%g/s
330     end;
331 gforget 1.11 contOCN=nansum(contOCN,3);
332     contTOT=contOCN+contICE;
333 gforget 1.9
334 gforget 1.1 %vertical divergence (air-sea fluxes or vertical adv/dif)
335     zdivOCN=SFLUX+oceSPflx;
336     zdivICE=-zdivOCN+SRELAX;
337     %in linear surface we omit :
338     if ~myparms.useNLFS; zdivOCN=zdivOCN-myparms.rhoconst*WSLTMASS; end;
339     %working approach for real fresh water (?) and virtual salt flux
340     if ~myparms.useRFWF|~myparms.useNLFS; zdivICE=-oceSflux; end;
341     %for deep ocean layer :
342     if kBudget>1;
343     zdivOCN=-(ADVr_SLT+DFrE_SLT+DFrI_SLT)./mygrid.RAC*myparms.rhoconst;
344     zdivOCN=zdivOCN+oceSPtnd;%.*msk;
345     end;
346     zdivTOT=zdivOCN+zdivICE;
347 gforget 1.9 %
348     if doMoreBudgetOutput;
349 gforget 1.12 if test3d;
350     trWtop=-(ADVr_SLT+DFrE_SLT+DFrI_SLT)*myparms.rhoconst;
351     trWtop=trWtop+oceSPtnd.*mk3D(mygrid.RAC,trWtop);
352     %
353     trWtop(:,:,1)=zdivOCN.*mygrid.RAC;
354     trWbot=trWtop(:,:,2:length(mygrid.RC));
355     trWbot(:,:,length(mygrid.RC))=0;
356     %
357     budgSo.trWtop=trWtop;%kg/s
358     budgSo.trWbot=trWbot;%kg/s
359     else;
360 gforget 1.10 budgSo.trWtop=mygrid.RAC.*zdivOCN; budgSo.trWbot=mygrid.RAC*0;%kg/s
361 gforget 1.12 end;
362     budgSi.trWtop=0*mygrid.RAC; budgSi.trWbot=budgSo.trWtop(:,:,1);%kg/s
363 gforget 1.9 end;
364 gforget 1.12
365 gforget 1.1 %horizontal divergence (advection and diffusion)
366 gforget 1.9 tmpUo=myparms.rhoconst*(ADVx_SLT+DFxE_SLT);
367     tmpVo=myparms.rhoconst*(ADVy_SLT+DFyE_SLT);
368 gforget 1.11 hdivOCN=calc_UV_conv(nansum(tmpUo,3),nansum(tmpVo,3));
369 gforget 1.9 tmpUi=myparms.SIsal0*(myparms.rhoi*DFxEHEFF+myparms.rhoi*ADVxHEFF);
370     tmpVi=myparms.SIsal0*(myparms.rhoi*DFyEHEFF+myparms.rhoi*ADVyHEFF);
371     hdivICE=calc_UV_conv(tmpUi,tmpVi); %no dh needed here
372 gforget 1.1 hdivTOT=hdivOCN+hdivICE;
373 gforget 1.9 if doMoreBudgetOutput;
374     budgSo.trU=tmpUo; budgSo.trV=tmpVo;%g/s
375     budgSi.trU=tmpUi; budgSi.trV=tmpVi;%g/s
376     end;
377 gforget 1.1
378     %compute global integrals:
379     %-------------------------
380     msk=mygrid.mskC(:,:,kBudget);
381     glo_salt_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
382     glo_salt_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
383     glo_salt_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
384    
385     %compute northern hemisphere integrals:
386     msk=mygrid.mskC(:,:,kBudget).*(mygrid.YC>0);
387     north_salt_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
388     north_salt_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
389     north_salt_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
390    
391     %and southern hemisphere integrals:
392     msk=mygrid.mskC(:,:,kBudget).*(mygrid.YC<=0);
393     south_salt_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
394     south_salt_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
395     south_salt_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
396 gforget 1.2
397 gforget 1.9 if doMoreBudgetOutput;
398     %list of budgets to output
399     listbudg={'budgMo','budgHo','budgSo'};
400     if kBudget==1; listbudg={listbudg{:},'budgMi','budgHi','budgSi'}; end;
401     %the actual output
402     for iibudg=1:length(listbudg);
403     %set directory name
404     dirbudg=dirMat;
405     if ~isempty(strfind(dirMat,['diags_set_' setDiags '/']))
406     dirbudg=fullfile(dirMat,'..',filesep);
407     end;
408     sufbudg='';
409     if kBudget>1; sufbudg=num2str(kBudget); end;
410     dirbudg=fullfile(dirbudg,['diags_set_' listbudg{iibudg} sufbudg],filesep);
411     %
412     if ~isdir(dirbudg); mkdir(dirbudg); end;
413     %set file name
414     filebudg=[listbudg{iibudg} '_' num2str(tt) '.mat'];
415     %output to file
416     eval(['tmpbudg=' listbudg{iibudg} ';']);
417     save([dirbudg filebudg],'-struct','tmpbudg');
418     end;
419     end;
420    
421 gforget 1.2 %===================== COMPUTATIONAL SEQUENCE ENDS =========================%
422     %===================== PLOTTING SEQUENCE BEGINS =========================%
423    
424     elseif userStep==-1;%plotting
425    
426     if isempty(setDiagsParams);
427     choicePlot={'all'};
428     elseif isnumeric(setDiagsParams{1})&length(setDiagsParams)==1;
429     choicePlot={'all'};
430     elseif isnumeric(setDiagsParams{1});
431     choicePlot={setDiagsParams{2:end}};
432     else;
433     choicePlot=setDiagsParams;
434     end;
435    
436 gforget 1.3 tt=[1:length(alldiag.listTimes)];
437     TT=alldiag.listTimes(tt);
438     nt=length(TT);
439    
440 gforget 1.2 if (kBudget==1)&(sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'mass')));
441    
442     %1.1) ocean+seaice mass budgets
443     %------------------------------
444     figureL;
445     %global volume budget:
446 heimbach 1.6 subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_vol_tot,'kg/m^2','Global Mean Mass (incl. ice)');
447 gforget 1.2 %add bp:
448     dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.glo_bp);
449     plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
450     %northern hemisphere budget:
451 heimbach 1.6 subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_vol_tot,'kg/m^2','Northern Mean Mass (incl. ice)');
452 gforget 1.2 %add bp:
453     dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.north_bp);
454     plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
455     %southern hemisphere budget:
456 heimbach 1.6 subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_vol_tot,'kg/m^2','Southern Mean Mass (incl. ice)');
457 gforget 1.2 %add bp:
458     dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.south_bp);
459     plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
460     %add to tex file
461     myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
462 heimbach 1.6 myCaption={myCaption{:},'mass budget (ocean+ice) in kg/m$^2$.'};
463 gforget 1.2 if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
464    
465     %1.2) ice mass budgets
466     %---------------------
467     figureL;
468 heimbach 1.6 subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_vol_ice,'kg/m^2','Global Mean Mass (only ice)');
469 gforget 1.2 dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.glo_bp);
470     plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
471 heimbach 1.6 subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_vol_ice,'kg/m^2','Northern Mean Mass (only ice)');
472 gforget 1.2 dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.north_bp);
473     plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
474 heimbach 1.6 subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_vol_ice,'kg/m^2','Southern Mean Mass (only ice)');
475 gforget 1.2 dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.south_bp);
476     plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
477     %add to tex file
478     myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
479 heimbach 1.6 myCaption={myCaption{:},'mass budget (ice only) in kg/m$^2$.'};
480 gforget 1.2 if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
481    
482     end;
483    
484     if (sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'mass')));
485    
486     %1.3) ocean mass budgets
487     %-----------------------
488     figureL;
489     %global volume budget:
490 heimbach 1.6 subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_vol_ocn,'kg/m^2','Global Mean Mass (only ocean)');
491 gforget 1.2 dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.glo_bp);
492     plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
493 heimbach 1.6 subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_vol_ocn,'kg/m^2','Northern Mean Mass (only ocean)');
494 gforget 1.2 dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.north_bp);
495     plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
496 heimbach 1.6 subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_vol_ocn,'kg/m^2','Southern Mean Mass (only ocean)');
497 gforget 1.2 dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.south_bp);
498     plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
499     %add to tex file
500     myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
501 heimbach 1.6 myCaption={myCaption{:},'mass budget (ocean only) in kg/m$^2$.'};
502 gforget 1.2 if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
503    
504     end;
505    
506     if (kBudget==1)&(sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'heat')));
507    
508     %2.1) ocean+seaice heat budgets
509     %------------------------------
510     figureL;
511 heimbach 1.6 subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_heat_tot,'J/m^2','Global Mean Ocean Heat (incl. ice)');
512     subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_heat_tot,'J/m^2','Northern Mean Ocean Heat (incl. ice)');
513     subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_heat_tot,'J/m^2','Southern Mean Ocean Heat (incl. ice)');
514 gforget 1.2 %add to tex file
515     myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
516 heimbach 1.6 myCaption={myCaption{:},'heat budget (ocean+ice) in J/m$^2$.'};
517 gforget 1.2 if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
518    
519     %2.2) ice heat budgets
520     %---------------------
521     figureL;
522 heimbach 1.6 subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_heat_ice,'J/m^2','Global Mean Ocean Heat (only ice)');
523     subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_heat_ice,'J/m^2','Northern Mean Ocean Heat (only ice)');
524     subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_heat_ice,'J/m^2','Southern Mean Ocean Heat (only ice)');
525 gforget 1.2 %add to tex file
526     myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
527 heimbach 1.6 myCaption={myCaption{:},'heat budget (ice only) in J/m$^2$.'};
528 gforget 1.2 if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
529    
530     end;
531    
532     if (sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'heat')));
533    
534     %2.3) ocean heat budgets
535     %-----------------------
536     figureL;
537 heimbach 1.6 subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_heat_ocn,'J/m^2','Global Mean Ocean Heat (only ocean)');
538     subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_heat_ocn,'J/m^2','Northern Mean Ocean Heat (only ocean)');
539     subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_heat_ocn,'J/m^2','Southern Mean Ocean Heat (only ocean)');
540 gforget 1.2 %add to tex file
541     myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
542 heimbach 1.6 myCaption={myCaption{:},'heat budget (ocean only) in J/m$^2$.'};
543 gforget 1.2 if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
544    
545     end;
546    
547     if (kBudget==1)&(sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'salt')));
548    
549     %3.1) ocean+seaice salt budgets
550     %------------------------------
551     figureL;
552 heimbach 1.6 subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_salt_tot,'g/m^2','Global Mean Ocean Salt (incl. ice)');
553     subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_salt_tot,'g/m^2','Northern Mean Ocean Salt (incl. ice)');
554     subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_salt_tot,'g/m^2','Southern Mean Ocean Salt (incl. ice)');
555 gforget 1.2 %add to tex file
556     myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
557 heimbach 1.6 myCaption={myCaption{:},'salt budget (ocean+ice) in g/m$^2$.'};
558 gforget 1.2 if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
559    
560     %2.2) ice salt budgets
561     %---------------------
562     figureL;
563 heimbach 1.6 subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_salt_ice,'g/m^2','Global Mean Ocean Salt (only ice)');
564     subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_salt_ice,'g/m^2','Northern Mean Ocean Salt (only ice)');
565     subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_salt_ice,'g/m^2','Southern Mean Ocean Salt (only ice)');
566 gforget 1.2 %add to tex file
567     myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
568 heimbach 1.6 myCaption={myCaption{:},'salt budget (ice only) in g/m$^2$.'};
569 gforget 1.2 if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
570    
571     end;
572    
573    
574     if (sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'salt')));
575    
576     %3.3) ocean salt budgets
577     %-----------------------
578     figureL;
579 gforget 1.8 subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_salt_ocn,'g/m^2','Global Mean Ocean Salt (only ocean)');
580     subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_salt_ocn,'g/m^2','Northern Mean Ocean Salt (only ocean)');
581     subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_salt_ocn,'g/m^2','Southern Mean Ocean Salt (only ocean)');
582 gforget 1.2 %add to tex file
583     myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
584 heimbach 1.6 myCaption={myCaption{:},'salt budget (ocean only) in g/m$^2$.'};
585 gforget 1.2 if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
586    
587     end;
588    
589 gforget 1.1 end;

  ViewVC Help
Powered by ViewVC 1.1.22