/[MITgcm]/MITgcm_contrib/gael/matlab_class/gcmfaces_calc/calc_budget_mass.m
ViewVC logotype

Annotation of /MITgcm_contrib/gael/matlab_class/gcmfaces_calc/calc_budget_mass.m

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


Revision 1.4 - (hide annotations) (download)
Sun Feb 7 22:19:52 2016 UTC (9 years, 5 months ago) by gforget
Branch: MAIN
Changes since 1.3: +15 -6 lines
- budgO.cont, budgO.hconv, etc. now in extensive form.

1 gforget 1.3 function [budgO,budgI,budgOI]=calc_budget_mass(kBudget);
2 gforget 1.1 % CALC_BUDGET_MASS(kBudget,doMoreBudgetOutput)
3    
4     gcmfaces_global;
5    
6     %get variables from caller routine:
7     %----------------------------------
8    
9     global myparms;
10    
11     list_variables={'ETAN','SIheff','SIhsnow','oceFWflx','SIatmFW','oceFWflx',...
12     'UVELMASS','VVELMASS',...
13     'ADVxHEFF','ADVxSNOW','DFxEHEFF','DFxESNOW',...
14     'ADVyHEFF','ADVySNOW','DFyEHEFF','DFyESNOW'};
15    
16     if kBudget>1; list_variables={list_variables{:},'WVELMASS'}; end;
17    
18     for vv=1:length(list_variables);
19     v = evalin('caller',list_variables{vv});
20     eval([list_variables{vv} '=v;']);
21     end;
22     clear v;
23    
24     test3d=length(size(UVELMASS{1}))>2;
25    
26     %compute mapped budget:
27     %----------------------
28    
29     %mass = myparms.rhoconst * sea level
30     contOCN=ETAN*myparms.rhoconst;
31     contICE=(SIheff*myparms.rhoi+SIhsnow*myparms.rhosn);
32     %for deep ocean layer :
33     if kBudget>1&myparms.useNLFS<2;
34     contOCN=0;
35     elseif kBudget>1;%rstar case
36     tmp1=mk3D(mygrid.DRF,mygrid.hFacC).*mygrid.hFacC;
37     tmp2=sum(tmp1(:,:,kBudget:length(mygrid.RC)),3)./mygrid.Depth;
38     contOCN=tmp2.*ETAN*myparms.rhoconst;
39     end;
40     %
41     contTOT=contOCN+contICE;
42     %
43     if test3d;
44     tmp1=mk3D(mygrid.DRF,mygrid.hFacC).*mygrid.hFacC;
45     tmp2=tmp1./mk3D(mygrid.Depth,tmp1);
46     tend=tmp2.*mk3D(ETAN,tmp2)*myparms.rhoconst;
47     else;
48     tend=contOCN;
49     end;
50 gforget 1.2 budgO.tend=mk3D(mygrid.RAC,tend).*tend;%kg/s
51     budgI.tend=mygrid.RAC.*contICE;%kg/s
52 gforget 1.1
53     %vertical divergence (air-sea fluxes or vertical advection)
54     zdivOCN=oceFWflx;
55     zdivICE=SIatmFW-oceFWflx;
56     %in virtual salt flux we omit :
57     if ~myparms.useRFWF; zdivOCN=0*zdivOCN; end;
58     %for deep ocean layer :
59     if kBudget>1; zdivOCN=-WVELMASS*myparms.rhoconst; end;
60     %
61     zdivTOT=zdivOCN+zdivICE;
62     %
63     if test3d;
64     trWtop=-WVELMASS*myparms.rhoconst;
65     %trWtop(:,:,1)=zdivOCN;
66     trWbot=trWtop(:,:,2:length(mygrid.RC));
67     trWbot(:,:,length(mygrid.RC))=0;
68     %
69 gforget 1.2 budgO.trWtop=mk3D(mygrid.RAC,trWtop).*trWtop;
70     budgO.trWbot=mk3D(mygrid.RAC,trWbot).*trWbot;%kg/s
71 gforget 1.1 else;
72 gforget 1.2 budgO.trWtop=mygrid.RAC.*zdivOCN; budgO.trWbot=mygrid.RAC*0;%kg/s
73 gforget 1.1 end;
74 gforget 1.2 budgI.trWtop=mygrid.RAC.*(zdivICE+zdivOCN); budgI.trWbot=mygrid.RAC.*zdivOCN;%kg/s
75 gforget 1.1
76     %horizontal divergence (advection and ice diffusion)
77     if test3d;
78     %3D UVELMASS,VVELMASS are multiplied by DRF
79     %(2D diagnostics are expectedly vertically integrated by MITgcm)
80     tmp1=mk3D(mygrid.DRF,UVELMASS);
81     UVELMASS=tmp1.*UVELMASS;
82     VVELMASS=tmp1.*VVELMASS;
83     end;
84     dxg=mk3D(mygrid.DXG,VVELMASS); dyg=mk3D(mygrid.DYG,UVELMASS);
85     tmpUo=myparms.rhoconst*dyg.*UVELMASS;
86     tmpVo=myparms.rhoconst*dxg.*VVELMASS;
87     hdivOCN=calc_UV_conv(nansum(tmpUo,3),nansum(tmpVo,3));
88 gforget 1.4 tmpUi=(myparms.rhoi*DFxEHEFF+myparms.rhosn*DFxESNOW+...
89     myparms.rhoi*ADVxHEFF+myparms.rhosn*ADVxSNOW);
90     tmpVi=(myparms.rhoi*DFyEHEFF+myparms.rhosn*DFyESNOW+...
91     myparms.rhoi*ADVyHEFF+myparms.rhosn*ADVySNOW);
92 gforget 1.1 hdivICE=calc_UV_conv(tmpUi,tmpVi); %dh needed is alerady in DFxEHEFF etc
93     hdivTOT=hdivOCN+hdivICE;
94     %
95 gforget 1.2 budgO.trU=tmpUo; budgO.trV=tmpVo;%kg/s
96     budgI.trU=tmpUi; budgI.trV=tmpVi;%kg/s
97    
98 gforget 1.3 %outputs:
99 gforget 1.4 if test3d;
100     tmp1=mk3D(mygrid.DRF,mygrid.hFacC).*mygrid.hFacC;
101     tmp1=mk3D(mygrid.RAC,mygrid.hFacC).*tmp1;
102     else;
103     tmp1=mygrid.RAC;
104     end;
105    
106     budgO.cont=tmp1.*contOCN;
107 gforget 1.3 budgO.hconv=hdivOCN;
108 gforget 1.4 budgO.zconv=tmp1.*zdivOCN;
109 gforget 1.3
110 gforget 1.4 budgI.cont=tmp1.*contICE;
111 gforget 1.3 budgI.hconv=hdivICE;
112 gforget 1.4 budgI.zconv=tmp1.*zdivICE;
113 gforget 1.3
114     budgOI.cont=budgO.cont+budgI.cont;
115     budgOI.hconv=budgO.hconv+budgI.hconv;
116     budgOI.zconv=budgO.zconv+budgI.zconv;
117    
118    
119    
120 gforget 1.1

  ViewVC Help
Powered by ViewVC 1.1.22