/[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.15 - (hide annotations) (download)
Mon Oct 5 16:51:45 2015 UTC (9 years, 9 months ago) by gforget
Branch: MAIN
CVS Tags: checkpoint65p
Changes since 1.14: +4 -2 lines
- improved treatment of directory names, switches, and driver.

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

  ViewVC Help
Powered by ViewVC 1.1.22