/[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.16 - (hide annotations) (download)
Sat Nov 7 14:17:25 2015 UTC (9 years, 8 months ago) by gforget
Branch: MAIN
CVS Tags: checkpoint65r, checkpoint65q
Changes since 1.15: +9 -5 lines
- add geothFlux in heat budget computation.

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

  ViewVC Help
Powered by ViewVC 1.1.22