/[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.11 - (hide annotations) (download)
Mon Feb 2 22:05:44 2015 UTC (10 years, 5 months ago) by gforget
Branch: MAIN
Changes since 1.10: +33 -19 lines 
- diags_diff_snapshots.m : use rdmds_meta (new), when 3d snapshots multiply
  THETA, SALT by DRF*hFac*ETAN (r* case) before computing tendencies
- diags_driver.m : add comment
- diags_list_times.m : update file list
- diags_pre_process.m : treat case of budg3d_snap_set1 when 3d snapshots
- diags_set_D.m : accomodate 3d output (snapshots and interior fluxes)
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     budgO.specs.top='free surface';
112     if kBudget>1; budgO.specs.top=['interface no. ' num2str(kBudget)]; end;
113     budgO.specs.bottom='sea floor';
114     budgI.specs.top='ocn-ice to atm interface';
115     budgI.specs.bottom='free surface';
116    
117     %here we output tendencies and fluxes in kg/s
118     budgMo=budgO; budgMi=budgI;
119     budgMo.specs.units='kg/s';%ocean only
120     budgMi.specs.units='kg/s';%ice only
121     %here we output tendencies and fluxes in Watts
122     budgHo=budgO; budgHi=budgI;
123     budgHo.specs.units='W';%ocean only
124     budgHi.specs.units='W';%ice only
125     %here we output tendencies and fluxes in g/s
126     budgSo=budgO; budgSi=budgI;
127     budgSo.specs.units='g/s';%ocean only
128     budgSi.specs.units='g/s';%ice only
129     end;
130    
131 gforget 1.1 %compute mapped budget:
132     %----------------------
133    
134     %mass = myparms.rhoconst * sea level
135     contOCN=ETAN*myparms.rhoconst;
136     contICE=(SIheff*myparms.rhoi+SIhsnow*myparms.rhosn);
137     %for deep ocean layer :
138     if kBudget>1&myparms.useNLFS<2;
139     contOCN=0;
140     elseif kBudget>1;%rstar case
141     tmp1=mk3D(mygrid.DRF,mygrid.hFacC).*mygrid.hFacC;
142     tmp2=sum(tmp1(:,:,kBudget:length(mygrid.RC)),3)./mygrid.Depth;
143     contOCN=tmp2.*ETAN*myparms.rhoconst;
144     end;
145     %
146     contTOT=contOCN+contICE;
147 gforget 1.9 %
148     if doMoreBudgetOutput;
149     budgMo.tend=mygrid.RAC.*contOCN;%kg/s
150     budgMi.tend=mygrid.RAC.*contICE;%kg/s
151     end;
152    
153 gforget 1.1 %vertical divergence (air-sea fluxes or vertical advection)
154     zdivOCN=oceFWflx;
155     zdivICE=SIatmFW-oceFWflx;
156     %in virtual salt flux we omit :
157     if ~myparms.useRFWF; zdivOCN=0*zdivOCN; end;
158     %for deep ocean layer :
159     if kBudget>1; zdivOCN=-WVELMASS*myparms.rhoconst; end;
160     %
161     zdivTOT=zdivOCN+zdivICE;
162 gforget 1.9 %
163     if doMoreBudgetOutput;
164 gforget 1.10 budgMo.trWtop=mygrid.RAC.*zdivOCN; budgMo.trWbot=mygrid.RAC*0;%kg/s
165     budgMi.trWtop=mygrid.RAC.*(zdivICE+zdivOCN); budgMi.trWbot=mygrid.RAC.*zdivOCN;%kg/s
166 gforget 1.9 end;
167    
168 gforget 1.1 %horizontal divergence (advection and ice diffusion)
169 gforget 1.11 if test3d;
170     %3D UVELMASS,VVELMASS are multiplied by DRF
171     %(2D diagnostics are expectedly vertically integrated by MITgcm)
172     tmp1=mk3D(mygrid.DRF,UVELMASS);
173     UVELMASS=tmp1.*UVELMASS;
174     VVELMASS=tmp1.*VVELMASS;
175     end;
176 gforget 1.9 dxg=mk3D(mygrid.DXG,VVELMASS); dyg=mk3D(mygrid.DYG,UVELMASS);
177     tmpUo=myparms.rhoconst*dyg.*UVELMASS;
178     tmpVo=myparms.rhoconst*dxg.*VVELMASS;
179 gforget 1.11 hdivOCN=calc_UV_conv(nansum(tmpUo,3),nansum(tmpVo,3));
180 gforget 1.9 tmpUi=(myparms.rhoi*DFxEHEFF+myparms.rhosn*DFxESNOW+myparms.rhoi*ADVxHEFF+myparms.rhosn*ADVxSNOW);
181     tmpVi=(myparms.rhoi*DFyEHEFF+myparms.rhosn*DFyESNOW+myparms.rhoi*ADVyHEFF+myparms.rhosn*ADVySNOW);
182     hdivICE=calc_UV_conv(tmpUi,tmpVi); %dh needed is alerady in DFxEHEFF etc
183 gforget 1.1 hdivTOT=hdivOCN+hdivICE;
184 gforget 1.9 if doMoreBudgetOutput;
185     budgMo.trU=tmpUo; budgMo.trV=tmpVo;%kg/s
186     budgMi.trU=tmpUi; budgMi.trV=tmpVi;%kg/s
187     end;
188 gforget 1.1
189     %bottom pressure for comparison:
190     bp=myparms.rhoconst/9.81*PHIBOT;
191    
192     %compute global integrals:
193     %-------------------------
194     msk=mygrid.mskC(:,:,kBudget);
195     glo_vol_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
196     glo_vol_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
197     glo_vol_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
198     glo_bp=nansum(bp.*msk.*mygrid.RAC)/nansum(msk.*mygrid.RAC);
199    
200     %compute northern hemisphere integrals:
201     msk=mygrid.mskC(:,:,kBudget).*(mygrid.YC>0);
202     north_vol_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
203     north_vol_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
204     north_vol_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
205     north_bp=nansum(bp.*msk.*mygrid.RAC)/nansum(msk.*mygrid.RAC);
206    
207     %and southern hemisphere integrals:
208     msk=mygrid.mskC(:,:,kBudget).*(mygrid.YC<=0);
209     south_vol_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
210     south_vol_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
211     south_vol_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
212     south_bp=nansum(bp.*msk.*mygrid.RAC)/nansum(msk.*mygrid.RAC);
213    
214     %=======HEAT=======
215 gforget 1.11
216 gforget 1.1 contOCN=myparms.rcp*THETA-myparms.rcp*AB_gT;
217     contICE=-myparms.flami*(SIheff*myparms.rhoi+SIhsnow*myparms.rhosn);
218 gforget 1.9 %
219     if doMoreBudgetOutput;
220 gforget 1.11 budgHo.tend=mk3D(mygrid.RAC,contOCN).*contOCN;%Watt
221 gforget 1.9 budgHi.tend=mygrid.RAC.*contICE;%Watt
222     end;
223 gforget 1.11 contOCN=nansum(contOCN,3);
224     contTOT=contOCN+contICE;
225 gforget 1.9
226 gforget 1.1 %vertical divergence (air-sea fluxes or vertical adv/dif)
227     zdivOCN=TFLUX;
228     zdivICE=-(SItflux+TFLUX-TRELAX);
229     %in linear surface we omit :
230     if ~myparms.useNLFS; zdivOCN=zdivOCN-myparms.rcp*WTHMASS; end;
231     %in virtual salt flux we omit :
232     if ~myparms.useRFWF|~myparms.useNLFS; zdivICE=zdivICE+SIaaflux; end;
233     %working approach for real fresh water (?) and virtual salt flux
234     if 0; zdivICE=-oceQnet-SIatmQnt-myparms.flami*(SIsnPrcp-SIacSubl); end;
235     %for deep ocean layer :
236     if kBudget>1;
237     zdivOCN=-(ADVr_TH+DFrE_TH+DFrI_TH)./mygrid.RAC*myparms.rcp;
238     dd=mygrid.RF(kBudget); msk=mygrid.mskC(:,:,kBudget);
239     swfrac=0.62*exp(dd/0.6)+(1-0.62)*exp(dd/20);
240     if dd<-200; swfrac=0; end;
241     zdivOCN=zdivOCN+swfrac*oceQsw;%.*msk;
242     end;
243     %
244     zdivTOT=zdivOCN+zdivICE;
245 gforget 1.9 %
246     if doMoreBudgetOutput;
247 gforget 1.10 budgHo.trWtop=mygrid.RAC.*zdivOCN; budgHo.trWbot=mygrid.RAC*0;%kg/s
248     budgHi.trWtop=mygrid.RAC.*(zdivICE+zdivOCN); budgHi.trWbot=mygrid.RAC.*zdivOCN;%kg/s
249 gforget 1.9 end;
250 gforget 1.11
251 gforget 1.1 %horizontal divergence (advection and diffusion)
252 gforget 1.9 tmpUo=myparms.rcp*(ADVx_TH+DFxE_TH); tmpVo=myparms.rcp*(ADVy_TH+DFyE_TH);
253 gforget 1.11 hdivOCN=calc_UV_conv(nansum(tmpUo,3),nansum(tmpVo,3));
254 gforget 1.9 tmpUi=-myparms.flami*(myparms.rhoi*DFxEHEFF+myparms.rhosn*DFxESNOW+myparms.rhoi*ADVxHEFF+myparms.rhosn*ADVxSNOW);
255     tmpVi=-myparms.flami*(myparms.rhoi*DFyEHEFF+myparms.rhosn*DFyESNOW+myparms.rhoi*ADVyHEFF+myparms.rhosn*ADVySNOW);
256     hdivICE=calc_UV_conv(tmpUi,tmpVi); %no dh needed here
257 gforget 1.1 hdivTOT=hdivOCN+hdivICE;
258 gforget 1.9 if doMoreBudgetOutput;
259     budgHo.trU=tmpUo; budgHo.trV=tmpVo;%Watt
260     budgHi.trU=tmpUi; budgHi.trV=tmpVi;%Watt
261     end;
262 gforget 1.1
263     %compute global integrals:
264     %-------------------------
265     msk=mygrid.mskC(:,:,kBudget);
266     glo_heat_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
267     glo_heat_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
268     glo_heat_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
269    
270     %compute northern hemisphere integrals:
271     msk=mygrid.mskC(:,:,kBudget).*(mygrid.YC>0);
272     north_heat_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
273     north_heat_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
274     north_heat_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
275    
276     %and southern hemisphere integrals:
277     msk=mygrid.mskC(:,:,kBudget).*(mygrid.YC<=0);
278     south_heat_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
279     south_heat_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
280     south_heat_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
281 gforget 1.11
282 gforget 1.1 %=======SALT=======
283 gforget 1.9
284 gforget 1.1 contOCN=myparms.rhoconst*SALT-myparms.rhoconst*AB_gS;
285     contICE=myparms.SIsal0*myparms.rhoi*SIheff;
286 gforget 1.9 %
287     if doMoreBudgetOutput;
288 gforget 1.11 budgSo.tend=mk3D(mygrid.RAC,3).*contOCN;%g/s
289 gforget 1.9 budgSi.tend=mygrid.RAC.*contICE;%g/s
290     end;
291 gforget 1.11 contOCN=nansum(contOCN,3);
292     contTOT=contOCN+contICE;
293 gforget 1.9
294 gforget 1.1 %vertical divergence (air-sea fluxes or vertical adv/dif)
295     zdivOCN=SFLUX+oceSPflx;
296     zdivICE=-zdivOCN+SRELAX;
297     %in linear surface we omit :
298     if ~myparms.useNLFS; zdivOCN=zdivOCN-myparms.rhoconst*WSLTMASS; end;
299     %working approach for real fresh water (?) and virtual salt flux
300     if ~myparms.useRFWF|~myparms.useNLFS; zdivICE=-oceSflux; end;
301     %for deep ocean layer :
302     if kBudget>1;
303     zdivOCN=-(ADVr_SLT+DFrE_SLT+DFrI_SLT)./mygrid.RAC*myparms.rhoconst;
304     zdivOCN=zdivOCN+oceSPtnd;%.*msk;
305     end;
306     zdivTOT=zdivOCN+zdivICE;
307 gforget 1.9 %
308     if doMoreBudgetOutput;
309 gforget 1.10 budgSo.trWtop=mygrid.RAC.*zdivOCN; budgSo.trWbot=mygrid.RAC*0;%kg/s
310     budgSi.trWtop=0*mygrid.RAC; budgSi.trWbot=budgSo.trWtop;%kg/s
311 gforget 1.9 end;
312    
313 gforget 1.1 %horizontal divergence (advection and diffusion)
314 gforget 1.9 tmpUo=myparms.rhoconst*(ADVx_SLT+DFxE_SLT);
315     tmpVo=myparms.rhoconst*(ADVy_SLT+DFyE_SLT);
316 gforget 1.11 hdivOCN=calc_UV_conv(nansum(tmpUo,3),nansum(tmpVo,3));
317 gforget 1.9 tmpUi=myparms.SIsal0*(myparms.rhoi*DFxEHEFF+myparms.rhoi*ADVxHEFF);
318     tmpVi=myparms.SIsal0*(myparms.rhoi*DFyEHEFF+myparms.rhoi*ADVyHEFF);
319     hdivICE=calc_UV_conv(tmpUi,tmpVi); %no dh needed here
320 gforget 1.1 hdivTOT=hdivOCN+hdivICE;
321 gforget 1.9 if doMoreBudgetOutput;
322     budgSo.trU=tmpUo; budgSo.trV=tmpVo;%g/s
323     budgSi.trU=tmpUi; budgSi.trV=tmpVi;%g/s
324     end;
325 gforget 1.1
326     %compute global integrals:
327     %-------------------------
328     msk=mygrid.mskC(:,:,kBudget);
329     glo_salt_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
330     glo_salt_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
331     glo_salt_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
332    
333     %compute northern hemisphere integrals:
334     msk=mygrid.mskC(:,:,kBudget).*(mygrid.YC>0);
335     north_salt_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
336     north_salt_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
337     north_salt_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
338    
339     %and southern hemisphere integrals:
340     msk=mygrid.mskC(:,:,kBudget).*(mygrid.YC<=0);
341     south_salt_tot=calc_budget_mean_mask(contTOT,zdivTOT,hdivTOT,msk);
342     south_salt_ocn=calc_budget_mean_mask(contOCN,zdivOCN,hdivOCN,msk);
343     south_salt_ice=calc_budget_mean_mask(contICE,zdivICE,hdivICE,msk);
344 gforget 1.2
345 gforget 1.9 if doMoreBudgetOutput;
346     %list of budgets to output
347     listbudg={'budgMo','budgHo','budgSo'};
348     if kBudget==1; listbudg={listbudg{:},'budgMi','budgHi','budgSi'}; end;
349     %the actual output
350     for iibudg=1:length(listbudg);
351     %set directory name
352     dirbudg=dirMat;
353     if ~isempty(strfind(dirMat,['diags_set_' setDiags '/']))
354     dirbudg=fullfile(dirMat,'..',filesep);
355     end;
356     sufbudg='';
357     if kBudget>1; sufbudg=num2str(kBudget); end;
358     dirbudg=fullfile(dirbudg,['diags_set_' listbudg{iibudg} sufbudg],filesep);
359     %
360     if ~isdir(dirbudg); mkdir(dirbudg); end;
361     %set file name
362     filebudg=[listbudg{iibudg} '_' num2str(tt) '.mat'];
363     %output to file
364     eval(['tmpbudg=' listbudg{iibudg} ';']);
365     save([dirbudg filebudg],'-struct','tmpbudg');
366     end;
367     end;
368    
369 gforget 1.2 %===================== COMPUTATIONAL SEQUENCE ENDS =========================%
370     %===================== PLOTTING SEQUENCE BEGINS =========================%
371    
372     elseif userStep==-1;%plotting
373    
374     if isempty(setDiagsParams);
375     choicePlot={'all'};
376     elseif isnumeric(setDiagsParams{1})&length(setDiagsParams)==1;
377     choicePlot={'all'};
378     elseif isnumeric(setDiagsParams{1});
379     choicePlot={setDiagsParams{2:end}};
380     else;
381     choicePlot=setDiagsParams;
382     end;
383    
384 gforget 1.3 tt=[1:length(alldiag.listTimes)];
385     TT=alldiag.listTimes(tt);
386     nt=length(TT);
387    
388 gforget 1.2 if (kBudget==1)&(sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'mass')));
389    
390     %1.1) ocean+seaice mass budgets
391     %------------------------------
392     figureL;
393     %global volume budget:
394 heimbach 1.6 subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_vol_tot,'kg/m^2','Global Mean Mass (incl. ice)');
395 gforget 1.2 %add bp:
396     dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.glo_bp);
397     plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
398     %northern hemisphere budget:
399 heimbach 1.6 subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_vol_tot,'kg/m^2','Northern Mean Mass (incl. ice)');
400 gforget 1.2 %add bp:
401     dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.north_bp);
402     plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
403     %southern hemisphere budget:
404 heimbach 1.6 subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_vol_tot,'kg/m^2','Southern Mean Mass (incl. ice)');
405 gforget 1.2 %add bp:
406     dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.south_bp);
407     plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
408     %add to tex file
409     myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
410 heimbach 1.6 myCaption={myCaption{:},'mass budget (ocean+ice) in kg/m$^2$.'};
411 gforget 1.2 if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
412    
413     %1.2) ice mass budgets
414     %---------------------
415     figureL;
416 heimbach 1.6 subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_vol_ice,'kg/m^2','Global Mean Mass (only ice)');
417 gforget 1.2 dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.glo_bp);
418     plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
419 heimbach 1.6 subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_vol_ice,'kg/m^2','Northern Mean Mass (only ice)');
420 gforget 1.2 dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.north_bp);
421     plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
422 heimbach 1.6 subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_vol_ice,'kg/m^2','Southern Mean Mass (only ice)');
423 gforget 1.2 dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.south_bp);
424     plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
425     %add to tex file
426     myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
427 heimbach 1.6 myCaption={myCaption{:},'mass budget (ice only) in kg/m$^2$.'};
428 gforget 1.2 if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
429    
430     end;
431    
432     if (sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'mass')));
433    
434     %1.3) ocean mass budgets
435     %-----------------------
436     figureL;
437     %global volume budget:
438 heimbach 1.6 subplot(3,1,1); disp_budget_mean_mask(TT,alldiag.glo_vol_ocn,'kg/m^2','Global Mean Mass (only ocean)');
439 gforget 1.2 dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.glo_bp);
440     plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
441 heimbach 1.6 subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_vol_ocn,'kg/m^2','Northern Mean Mass (only ocean)');
442 gforget 1.2 dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.north_bp);
443     plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
444 heimbach 1.6 subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_vol_ocn,'kg/m^2','Southern Mean Mass (only ocean)');
445 gforget 1.2 dt=median(diff(TT))*86400; bp=dt*cumsum(alldiag.south_bp);
446     plot(TT,bp,'k'); aa=legend; bb=get(aa,'String'); bb={bb{:},'bp'}; legend(bb,'Orientation','horizontal');
447     %add to tex file
448     myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
449 heimbach 1.6 myCaption={myCaption{:},'mass budget (ocean only) in kg/m$^2$.'};
450 gforget 1.2 if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
451    
452     end;
453    
454     if (kBudget==1)&(sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'heat')));
455    
456     %2.1) ocean+seaice heat budgets
457     %------------------------------
458     figureL;
459 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)');
460     subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_heat_tot,'J/m^2','Northern Mean Ocean Heat (incl. ice)');
461     subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_heat_tot,'J/m^2','Southern Mean Ocean Heat (incl. ice)');
462 gforget 1.2 %add to tex file
463     myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
464 heimbach 1.6 myCaption={myCaption{:},'heat budget (ocean+ice) in J/m$^2$.'};
465 gforget 1.2 if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
466    
467     %2.2) ice heat budgets
468     %---------------------
469     figureL;
470 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)');
471     subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_heat_ice,'J/m^2','Northern Mean Ocean Heat (only ice)');
472     subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_heat_ice,'J/m^2','Southern Mean Ocean Heat (only ice)');
473 gforget 1.2 %add to tex file
474     myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
475 heimbach 1.6 myCaption={myCaption{:},'heat budget (ice only) in J/m$^2$.'};
476 gforget 1.2 if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
477    
478     end;
479    
480     if (sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'heat')));
481    
482     %2.3) ocean heat budgets
483     %-----------------------
484     figureL;
485 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)');
486     subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_heat_ocn,'J/m^2','Northern Mean Ocean Heat (only ocean)');
487     subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_heat_ocn,'J/m^2','Southern Mean Ocean Heat (only ocean)');
488 gforget 1.2 %add to tex file
489     myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
490 heimbach 1.6 myCaption={myCaption{:},'heat budget (ocean only) in J/m$^2$.'};
491 gforget 1.2 if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
492    
493     end;
494    
495     if (kBudget==1)&(sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'salt')));
496    
497     %3.1) ocean+seaice salt budgets
498     %------------------------------
499     figureL;
500 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)');
501     subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_salt_tot,'g/m^2','Northern Mean Ocean Salt (incl. ice)');
502     subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_salt_tot,'g/m^2','Southern Mean Ocean Salt (incl. ice)');
503 gforget 1.2 %add to tex file
504     myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
505 heimbach 1.6 myCaption={myCaption{:},'salt budget (ocean+ice) in g/m$^2$.'};
506 gforget 1.2 if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
507    
508     %2.2) ice salt budgets
509     %---------------------
510     figureL;
511 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)');
512     subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_salt_ice,'g/m^2','Northern Mean Ocean Salt (only ice)');
513     subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_salt_ice,'g/m^2','Southern Mean Ocean Salt (only 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{:},'salt budget (ice only) in g/m$^2$.'};
517 gforget 1.2 if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
518    
519     end;
520    
521    
522     if (sum(strcmp(choicePlot,'all'))|sum(strcmp(choicePlot,'salt')));
523    
524     %3.3) ocean salt budgets
525     %-----------------------
526     figureL;
527 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)');
528     subplot(3,1,2); disp_budget_mean_mask(TT,alldiag.north_salt_ocn,'g/m^2','Northern Mean Ocean Salt (only ocean)');
529     subplot(3,1,3); disp_budget_mean_mask(TT,alldiag.south_salt_ocn,'g/m^2','Southern Mean Ocean Salt (only ocean)');
530 gforget 1.2 %add to tex file
531     myCaption={myYmeanTxt,' global (upper) north (mid) and south (lower), '};
532 heimbach 1.6 myCaption={myCaption{:},'salt budget (ocean only) in g/m$^2$.'};
533 gforget 1.2 if addToTex&multiTimes; write2tex(fileTex,2,myCaption,gcf); elseif ~multiTimes; close; end;
534    
535     end;
536    
537 gforget 1.1 end;
  ViewVC Help
Powered by ViewVC 1.1.22