/[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.14 - (hide annotations) (download)
Tue Jun 9 21:18:24 2015 UTC (10 years, 1 month ago) by gforget
Branch: MAIN
Changes since 1.13: +1 -1 lines 
- diags_set*.m : add [dirModel 'diags/'] to listSubdirs for case when diags is un-sorted.
- diags_pre_process.m : remove addition of pwd ahead of dirModel, link MITprof files
  from dirModel in case they were not found in hard-coded directory.
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.14 listSubdirs={[dirMat 'BUDG/' ],[dirMat '../BUDG/' ],[dirModel 'diags/BUDG/'],[dirModel 'diags/']};
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 gforget 1.13 swtop=mk3D(swfrac,trWtop).*mk3D(mygrid.RAC.*oceQsw,trWtop);
276 gforget 1.12 swtop(isnan(mygrid.mskC))=0;
277     trWtop=trWtop+swtop;
278     %
279 gforget 1.13 trWtop(:,:,1)=zdivOCN.*mygrid.RAC;
280 gforget 1.12 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 gforget 1.13 nr=length(mygrid.RC);
351 gforget 1.12 trWtop=-(ADVr_SLT+DFrE_SLT+DFrI_SLT)*myparms.rhoconst;
352 gforget 1.13 tmp1=mk3D(oceSPflx,oceSPtnd)-cumsum(oceSPtnd,3);
353     tmp1=tmp1.*mk3D(mygrid.RAC,tmp1);
354     trWtop(:,:,2:nr)=trWtop(:,:,2:nr)+tmp1(:,:,1:nr-1);
355 gforget 1.12 %
356     trWtop(:,:,1)=zdivOCN.*mygrid.RAC;
357     trWbot=trWtop(:,:,2:length(mygrid.RC));
358     trWbot(:,:,length(mygrid.RC))=0;
359     %
360     budgSo.trWtop=trWtop;%kg/s
361     budgSo.trWbot=trWbot;%kg/s
362     else;
363 gforget 1.10 budgSo.trWtop=mygrid.RAC.*zdivOCN; budgSo.trWbot=mygrid.RAC*0;%kg/s
364 gforget 1.12 end;
365     budgSi.trWtop=0*mygrid.RAC; budgSi.trWbot=budgSo.trWtop(:,:,1);%kg/s
366 gforget 1.9 end;
367 gforget 1.12
368 gforget 1.1 %horizontal divergence (advection and diffusion)
369 gforget 1.9 tmpUo=myparms.rhoconst*(ADVx_SLT+DFxE_SLT);
370     tmpVo=myparms.rhoconst*(ADVy_SLT+DFyE_SLT);
371 gforget 1.11 hdivOCN=calc_UV_conv(nansum(tmpUo,3),nansum(tmpVo,3));
372 gforget 1.9 tmpUi=myparms.SIsal0*(myparms.rhoi*DFxEHEFF+myparms.rhoi*ADVxHEFF);
373     tmpVi=myparms.SIsal0*(myparms.rhoi*DFyEHEFF+myparms.rhoi*ADVyHEFF);
374     hdivICE=calc_UV_conv(tmpUi,tmpVi); %no dh needed here
375 gforget 1.1 hdivTOT=hdivOCN+hdivICE;
376 gforget 1.9 if doMoreBudgetOutput;
377     budgSo.trU=tmpUo; budgSo.trV=tmpVo;%g/s
378     budgSi.trU=tmpUi; budgSi.trV=tmpVi;%g/s
379     end;
380 gforget 1.1
381     %compute global integrals:
382     %-------------------------
383     msk=mygrid.mskC(:,:,kBudget);
384     glo_salt_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
385     glo_salt_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
386     glo_salt_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
387    
388     %compute northern hemisphere integrals:
389     msk=mygrid.mskC(:,:,kBudget).*(mygrid.YC>0);
390     north_salt_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
391     north_salt_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
392     north_salt_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
393    
394     %and southern hemisphere integrals:
395     msk=mygrid.mskC(:,:,kBudget).*(mygrid.YC<=0);
396     south_salt_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
397     south_salt_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
398     south_salt_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
399 gforget 1.2
400 gforget 1.9 if doMoreBudgetOutput;
401     %list of budgets to output
402     listbudg={'budgMo','budgHo','budgSo'};
403     if kBudget==1; listbudg={listbudg{:},'budgMi','budgHi','budgSi'}; end;
404     %the actual output
405     for iibudg=1:length(listbudg);
406     %set directory name
407     dirbudg=dirMat;
408     if ~isempty(strfind(dirMat,['diags_set_' setDiags '/']))
409     dirbudg=fullfile(dirMat,'..',filesep);
410     end;
411     sufbudg='';
412     if kBudget>1; sufbudg=num2str(kBudget); end;
413     dirbudg=fullfile(dirbudg,['diags_set_' listbudg{iibudg} sufbudg],filesep);
414     %
415     if ~isdir(dirbudg); mkdir(dirbudg); end;
416     %set file name
417     filebudg=[listbudg{iibudg} '_' num2str(tt) '.mat'];
418     %output to file
419     eval(['tmpbudg=' listbudg{iibudg} ';']);
420     save([dirbudg filebudg],'-struct','tmpbudg');
421     end;
422     end;
423 gforget 1.13
424 gforget 1.2 %===================== COMPUTATIONAL SEQUENCE ENDS =========================%
425     %===================== PLOTTING SEQUENCE BEGINS =========================%
426    
427     elseif userStep==-1;%plotting
428    
429     if isempty(setDiagsParams);
430     choicePlot={'all'};
431     elseif isnumeric(setDiagsParams{1})&length(setDiagsParams)==1;
432     choicePlot={'all'};
433     elseif isnumeric(setDiagsParams{1});
434     choicePlot={setDiagsParams{2:end}};
435     else;
436     choicePlot=setDiagsParams;
437     end;
438    
439 gforget 1.3 tt=[1:length(alldiag.listTimes)];
440     TT=alldiag.listTimes(tt);
441     nt=length(TT);
442    
443 gforget 1.2 if (kBudget==1)&(sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'mass')));
444    
445     %1.1) ocean+seaice mass budgets
446     %------------------------------
447     figureL;
448     %global volume budget:
449 heimbach 1.6 subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_vol_tot,'kg/m^2','Global Mean Mass (incl. ice)');
450 gforget 1.2 %add bp:
451     dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.glo_bp);
452     plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
453     %northern hemisphere budget:
454 heimbach 1.6 subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_vol_tot,'kg/m^2','Northern Mean Mass (incl. ice)');
455 gforget 1.2 %add bp:
456     dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.north_bp);
457     plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
458     %southern hemisphere budget:
459 heimbach 1.6 subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_vol_tot,'kg/m^2','Southern Mean Mass (incl. ice)');
460 gforget 1.2 %add bp:
461     dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.south_bp);
462     plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
463     %add to tex file
464     myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
465 heimbach 1.6 myCaption={myCaption{:},'mass budget (ocean+ice) in kg/m$^2$.'};
466 gforget 1.2 if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
467    
468     %1.2) ice mass budgets
469     %---------------------
470     figureL;
471 heimbach 1.6 subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_vol_ice,'kg/m^2','Global Mean Mass (only ice)');
472 gforget 1.2 dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.glo_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,2); disp_budget_mean_mask(TT,alldiag.north_vol_ice,'kg/m^2','Northern Mean Mass (only ice)');
475 gforget 1.2 dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.north_bp);
476     plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
477 heimbach 1.6 subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_vol_ice,'kg/m^2','Southern Mean Mass (only ice)');
478 gforget 1.2 dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.south_bp);
479     plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
480     %add to tex file
481     myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
482 heimbach 1.6 myCaption={myCaption{:},'mass budget (ice only) in kg/m$^2$.'};
483 gforget 1.2 if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
484    
485     end;
486    
487     if (sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'mass')));
488    
489     %1.3) ocean mass budgets
490     %-----------------------
491     figureL;
492     %global volume budget:
493 heimbach 1.6 subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_vol_ocn,'kg/m^2','Global Mean Mass (only ocean)');
494 gforget 1.2 dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.glo_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,2); disp_budget_mean_mask(TT,alldiag.north_vol_ocn,'kg/m^2','Northern Mean Mass (only ocean)');
497 gforget 1.2 dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.north_bp);
498     plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
499 heimbach 1.6 subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_vol_ocn,'kg/m^2','Southern Mean Mass (only ocean)');
500 gforget 1.2 dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.south_bp);
501     plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
502     %add to tex file
503     myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
504 heimbach 1.6 myCaption={myCaption{:},'mass budget (ocean only) in kg/m$^2$.'};
505 gforget 1.2 if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
506    
507     end;
508    
509     if (kBudget==1)&(sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'heat')));
510    
511     %2.1) ocean+seaice heat budgets
512     %------------------------------
513     figureL;
514 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)');
515     subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_heat_tot,'J/m^2','Northern Mean Ocean Heat (incl. ice)');
516     subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_heat_tot,'J/m^2','Southern Mean Ocean Heat (incl. ice)');
517 gforget 1.2 %add to tex file
518     myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
519 heimbach 1.6 myCaption={myCaption{:},'heat budget (ocean+ice) in J/m$^2$.'};
520 gforget 1.2 if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
521    
522     %2.2) ice heat budgets
523     %---------------------
524     figureL;
525 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)');
526     subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_heat_ice,'J/m^2','Northern Mean Ocean Heat (only ice)');
527     subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_heat_ice,'J/m^2','Southern Mean Ocean Heat (only ice)');
528 gforget 1.2 %add to tex file
529     myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
530 heimbach 1.6 myCaption={myCaption{:},'heat budget (ice only) in J/m$^2$.'};
531 gforget 1.2 if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
532    
533     end;
534    
535     if (sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'heat')));
536    
537     %2.3) ocean heat budgets
538     %-----------------------
539     figureL;
540 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)');
541     subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_heat_ocn,'J/m^2','Northern Mean Ocean Heat (only ocean)');
542     subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_heat_ocn,'J/m^2','Southern Mean Ocean Heat (only ocean)');
543 gforget 1.2 %add to tex file
544     myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
545 heimbach 1.6 myCaption={myCaption{:},'heat budget (ocean only) in J/m$^2$.'};
546 gforget 1.2 if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
547    
548     end;
549    
550     if (kBudget==1)&(sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'salt')));
551    
552     %3.1) ocean+seaice salt budgets
553     %------------------------------
554     figureL;
555 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)');
556     subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_salt_tot,'g/m^2','Northern Mean Ocean Salt (incl. ice)');
557     subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_salt_tot,'g/m^2','Southern Mean Ocean Salt (incl. ice)');
558 gforget 1.2 %add to tex file
559     myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
560 heimbach 1.6 myCaption={myCaption{:},'salt budget (ocean+ice) in g/m$^2$.'};
561 gforget 1.2 if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
562    
563     %2.2) ice salt budgets
564     %---------------------
565     figureL;
566 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)');
567     subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_salt_ice,'g/m^2','Northern Mean Ocean Salt (only ice)');
568     subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_salt_ice,'g/m^2','Southern Mean Ocean Salt (only ice)');
569 gforget 1.2 %add to tex file
570     myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
571 heimbach 1.6 myCaption={myCaption{:},'salt budget (ice only) in g/m$^2$.'};
572 gforget 1.2 if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
573    
574     end;
575    
576    
577     if (sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'salt')));
578    
579     %3.3) ocean salt budgets
580     %-----------------------
581     figureL;
582 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)');
583     subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_salt_ocn,'g/m^2','Northern Mean Ocean Salt (only ocean)');
584     subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_salt_ocn,'g/m^2','Southern Mean Ocean Salt (only ocean)');
585 gforget 1.2 %add to tex file
586     myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
587 heimbach 1.6 myCaption={myCaption{:},'salt budget (ocean only) in g/m$^2$.'};
588 gforget 1.2 if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
589    
590     end;
591    
592 gforget 1.1 end;
  ViewVC Help
Powered by ViewVC 1.1.22