matlab_class/gcmfaces_calc/calc_budget_mass.m

function [budgO,budgI,budgOI]=calc_budget_mass(kBudget);
% CALC_BUDGET_MASS(kBudget,doMoreBudgetOutput)

gcmfaces_global;

%get variables from caller routine:
%----------------------------------

global myparms;

list_variables={'ETAN','SIheff','SIhsnow','oceFWflx','SIatmFW','oceFWflx',...
                'UVELMASS','VVELMASS',...
                'ADVxHEFF','ADVxSNOW','DFxEHEFF','DFxESNOW',...
                'ADVyHEFF','ADVySNOW','DFyEHEFF','DFyESNOW'};

for vv=1:length(list_variables);
  v = evalin('caller',list_variables{vv});
  eval([list_variables{vv} '=v;']);
end;
clear v;

test3d=length(size(UVELMASS{1}))>2;

if test3d|kBudget>1;
  list_variables={'WVELMASS'};
  for vv=1:length(list_variables);
    v = evalin('caller',list_variables{vv});
    eval([list_variables{vv} '=v;']);
  end;
  clear v;
end;

%compute mapped budget:
%----------------------

%mass = myparms.rhoconst * sea level
budgO.tend=ETAN*myparms.rhoconst;
budgI.tend=(SIheff*myparms.rhoi+SIhsnow*myparms.rhosn);
%for deep ocean layer :
if kBudget>1&myparms.useNLFS<2;
  budgO.tend=0;
elseif kBudget>1;%rstar case
  tmp1=mk3D(mygrid.DRF,mygrid.hFacC).*mygrid.hFacC;
  tmp2=sum(tmp1(:,:,kBudget:length(mygrid.RC)),3)./mygrid.Depth;
  budgO.tend=tmp2.*ETAN*myparms.rhoconst;
end;
%
if test3d;
  tmp1=mk3D(mygrid.DRF,mygrid.hFacC).*mygrid.hFacC;
  tmp2=tmp1./mk3D(mygrid.Depth,tmp1);
  budgO.tend=tmp2.*mk3D(ETAN,tmp2)*myparms.rhoconst;
end;
%
tmptend=mk3D(mygrid.RAC,budgO.tend).*budgO.tend;%kg/s
budgO.fluxes.tend=tmptend;
budgO.tend=nansum(tmptend,3);
budgI.tend=mygrid.RAC.*budgI.tend;%kg/s
budgOI.tend=budgO.tend+budgI.tend;

%vertical divergence (air-sea fluxes or vertical advection)
budgO.zconv=oceFWflx;
budgI.zconv=SIatmFW-oceFWflx;
%in virtual salt flux we omit :
if ~myparms.useRFWF; budgO.zconv=0*budgO.zconv; end;
%for deep ocean layer :
if kBudget>1; budgO.zconv=-WVELMASS*myparms.rhoconst; end;
%
if test3d;
  trWtop=-WVELMASS*myparms.rhoconst;
  %trWtop(:,:,1)=budgO.zconv;
  trWbot=trWtop(:,:,2:length(mygrid.RC));
  trWbot(:,:,length(mygrid.RC))=0;
  %
  budgO.fluxes.trWtop=mk3D(mygrid.RAC,trWtop).*trWtop;
  budgO.fluxes.trWbot=mk3D(mygrid.RAC,trWbot).*trWbot;%kg/s
else;
  budgO.fluxes.trWtop=-mygrid.RAC.*budgO.zconv; 
  budgO.fluxes.trWbot=mygrid.RAC*0;%kg/s
end;
budgI.fluxes.trWtop=-mygrid.RAC.*(budgI.zconv+budgO.zconv); 
budgI.fluxes.trWbot=-mygrid.RAC.*budgO.zconv;%kg/s
%
budgO.zconv=mk3D(mygrid.RAC,budgO.zconv).*budgO.zconv;%kg/s
budgI.zconv=mygrid.RAC.*budgI.zconv;%kg/s
budgOI.zconv=budgO.zconv+budgI.zconv;

%horizontal divergence (advection and ice diffusion)
if test3d; 
  %3D UVELMASS,VVELMASS are multiplied by DRF
  %(2D diagnostics are expectedly vertically integrated by MITgcm)
  tmp1=mk3D(mygrid.DRF,UVELMASS);
  UVELMASS=tmp1.*UVELMASS;
  VVELMASS=tmp1.*VVELMASS;
end;
dxg=mk3D(mygrid.DXG,VVELMASS); dyg=mk3D(mygrid.DYG,UVELMASS);
tmpUo=myparms.rhoconst*dyg.*UVELMASS;
tmpVo=myparms.rhoconst*dxg.*VVELMASS;
budgO.hconv=calc_UV_conv(nansum(tmpUo,3),nansum(tmpVo,3));
tmpUi=(myparms.rhoi*DFxEHEFF+myparms.rhosn*DFxESNOW+...
       myparms.rhoi*ADVxHEFF+myparms.rhosn*ADVxSNOW);
tmpVi=(myparms.rhoi*DFyEHEFF+myparms.rhosn*DFyESNOW+...
       myparms.rhoi*ADVyHEFF+myparms.rhosn*ADVySNOW);
budgI.hconv=calc_UV_conv(tmpUi,tmpVi); %dh needed is alerady in DFxEHEFF etc
%
budgOI.hconv=budgO.hconv;
budgOI.hconv(:,:,1)=budgOI.hconv(:,:,1)+budgI.hconv;
%
budgO.fluxes.trU=tmpUo; budgO.fluxes.trV=tmpVo;%kg/s
budgI.fluxes.trU=tmpUi; budgI.fluxes.trV=tmpVi;%kg/s

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			for vv=1:length(list_variables);
17			v = evalin('caller',list_variables{vv});
18			eval([list_variables{vv} '=v;']);
19			end;
20			clear v;
21
22			test3d=length(size(UVELMASS{1}))>2;
23
24	gforget	1.7	if test3d\|kBudget>1;
25			list_variables={'WVELMASS'};
26			for vv=1:length(list_variables);
27			v = evalin('caller',list_variables{vv});
28			eval([list_variables{vv} '=v;']);
29			end;
30			clear v;
31			end;
32
33	gforget	1.1	%compute mapped budget:
34			%----------------------
35
36			%mass = myparms.rhoconst * sea level
37	gforget	1.5	budgO.tend=ETAN*myparms.rhoconst;
38			budgI.tend=(SIheffmyparms.rhoi+SIhsnowmyparms.rhosn);
39	gforget	1.1	%for deep ocean layer :
40			if kBudget>1&myparms.useNLFS<2;
41	gforget	1.5	budgO.tend=0;
42	gforget	1.1	elseif kBudget>1;%rstar case
43	gforget	1.5	tmp1=mk3D(mygrid.DRF,mygrid.hFacC).*mygrid.hFacC;
44			tmp2=sum(tmp1(:,:,kBudget:length(mygrid.RC)),3)./mygrid.Depth;
45			budgO.tend=tmp2.ETANmyparms.rhoconst;
46	gforget	1.1	end;
47			%
48			if test3d;
49			tmp1=mk3D(mygrid.DRF,mygrid.hFacC).*mygrid.hFacC;
50			tmp2=tmp1./mk3D(mygrid.Depth,tmp1);
51	gforget	1.5	budgO.tend=tmp2.mk3D(ETAN,tmp2)myparms.rhoconst;
52	gforget	1.1	end;
53	gforget	1.5	%
54	gforget	1.7	tmptend=mk3D(mygrid.RAC,budgO.tend).*budgO.tend;%kg/s
55			budgO.fluxes.tend=tmptend;
56			budgO.tend=nansum(tmptend,3);
57	gforget	1.5	budgI.tend=mygrid.RAC.*budgI.tend;%kg/s
58			budgOI.tend=budgO.tend+budgI.tend;
59	gforget	1.1
60			%vertical divergence (air-sea fluxes or vertical advection)
61	gforget	1.5	budgO.zconv=oceFWflx;
62			budgI.zconv=SIatmFW-oceFWflx;
63	gforget	1.1	%in virtual salt flux we omit :
64	gforget	1.5	if ~myparms.useRFWF; budgO.zconv=0*budgO.zconv; end;
65	gforget	1.1	%for deep ocean layer :
66	gforget	1.5	if kBudget>1; budgO.zconv=-WVELMASS*myparms.rhoconst; end;
67	gforget	1.1	%
68			if test3d;
69			trWtop=-WVELMASS*myparms.rhoconst;
70	gforget	1.5	%trWtop(:,:,1)=budgO.zconv;
71	gforget	1.1	trWbot=trWtop(:,:,2:length(mygrid.RC));
72			trWbot(:,:,length(mygrid.RC))=0;
73			%
74	gforget	1.5	budgO.fluxes.trWtop=mk3D(mygrid.RAC,trWtop).*trWtop;
75			budgO.fluxes.trWbot=mk3D(mygrid.RAC,trWbot).*trWbot;%kg/s
76	gforget	1.1	else;
77	gforget	1.6	budgO.fluxes.trWtop=-mygrid.RAC.*budgO.zconv;
78	gforget	1.5	budgO.fluxes.trWbot=mygrid.RAC*0;%kg/s
79	gforget	1.1	end;
80	gforget	1.6	budgI.fluxes.trWtop=-mygrid.RAC.*(budgI.zconv+budgO.zconv);
81			budgI.fluxes.trWbot=-mygrid.RAC.*budgO.zconv;%kg/s
82	gforget	1.5	%
83			budgO.zconv=mk3D(mygrid.RAC,budgO.zconv).*budgO.zconv;%kg/s
84			budgI.zconv=mygrid.RAC.*budgI.zconv;%kg/s
85			budgOI.zconv=budgO.zconv+budgI.zconv;
86	gforget	1.1
87			%horizontal divergence (advection and ice diffusion)
88			if test3d;
89			%3D UVELMASS,VVELMASS are multiplied by DRF
90			%(2D diagnostics are expectedly vertically integrated by MITgcm)
91			tmp1=mk3D(mygrid.DRF,UVELMASS);
92			UVELMASS=tmp1.*UVELMASS;
93			VVELMASS=tmp1.*VVELMASS;
94			end;
95			dxg=mk3D(mygrid.DXG,VVELMASS); dyg=mk3D(mygrid.DYG,UVELMASS);
96			tmpUo=myparms.rhoconstdyg.UVELMASS;
97			tmpVo=myparms.rhoconstdxg.VVELMASS;
98	gforget	1.5	budgO.hconv=calc_UV_conv(nansum(tmpUo,3),nansum(tmpVo,3));
99	gforget	1.4	tmpUi=(myparms.rhoiDFxEHEFF+myparms.rhosnDFxESNOW+...
100			myparms.rhoiADVxHEFF+myparms.rhosnADVxSNOW);
101			tmpVi=(myparms.rhoiDFyEHEFF+myparms.rhosnDFyESNOW+...
102			myparms.rhoiADVyHEFF+myparms.rhosnADVySNOW);
103	gforget	1.5	budgI.hconv=calc_UV_conv(tmpUi,tmpVi); %dh needed is alerady in DFxEHEFF etc
104	gforget	1.7	%
105			budgOI.hconv=budgO.hconv;
106			budgOI.hconv(:,:,1)=budgOI.hconv(:,:,1)+budgI.hconv;
107	gforget	1.1	%
108	gforget	1.5	budgO.fluxes.trU=tmpUo; budgO.fluxes.trV=tmpVo;%kg/s
109			budgI.fluxes.trU=tmpUi; budgI.fluxes.trV=tmpVi;%kg/s
110	gforget	1.1