matlab_class/gcmfaces_calc/calc_budget_salt.m

function [budgO,budgI,budgOI]=calc_budget_salt(kBudget);
% CALC_BUDGET_SALT(kBudget)

gcmfaces_global;

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

global myparms;

list_variables={'SALT','AB_gS','SRELAX','SIheff',...
'SFLUX','oceSPflx','oceSflux','WSLTMASS',...
'ADVx_SLT','DFxE_SLT','ADVy_SLT','DFyE_SLT',...
'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(ADVx_SLT{1}))>2;

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

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

budgO.tend=myparms.rhoconst*SALT-myparms.rhoconst*AB_gS;
budgI.tend=myparms.SIsal0*myparms.rhoi*SIheff;
%
tmptend=mk3D(mygrid.RAC,budgO.tend).*budgO.tend;%g/s
budgO.fluxes.tend=tmptend;
budgO.tend=nansum(tmptend,3);
budgI.tend=mygrid.RAC.*budgI.tend;%g/s
%
budgOI.tend=budgO.tend+budgI.tend;

%vertical divergence (air-sea fluxes or vertical adv/dif)
budgO.zconv=SFLUX+oceSPflx;
budgI.zconv=-budgO.zconv+SRELAX;
%in linear surface we omit :
if ~myparms.useNLFS; budgO.zconv=budgO.zconv-myparms.rhoconst*WSLTMASS; end;
%working approach for real fresh water (?) and virtual salt flux
if ~myparms.useRFWF|~myparms.useNLFS; budgI.zconv=-oceSflux; end;
%
budgO.zdia=budgO.zconv;
%for deep ocean layer :
if kBudget>1;
  budgO.zconv=-(ADVr_SLT+DFrE_SLT+DFrI_SLT)./mygrid.RAC*myparms.rhoconst;
  budgO.zconv=budgO.zconv+oceSPtnd;%.*msk;
  budgO.zdia=-(DFrE_SLT+DFrI_SLT)./mygrid.RAC*myparms.rhoconst;
  budgO.zdia=budgO.zdia+oceSPtnd;%.*msk;
end;
%
if test3d;
  nr=length(mygrid.RC);
  trWtop=-(ADVr_SLT+DFrE_SLT+DFrI_SLT)*myparms.rhoconst;
  tmp1=mk3D(oceSPflx,oceSPtnd)-cumsum(oceSPtnd,3);
  tmp1=tmp1.*mk3D(mygrid.RAC,tmp1);
  trWtop(:,:,2:nr)=trWtop(:,:,2:nr)+tmp1(:,:,1:nr-1);
  %
  trWtop(:,:,1)=budgO.zconv.*mygrid.RAC;
  trWbot=trWtop(:,:,2:length(mygrid.RC));
  trWbot(:,:,length(mygrid.RC))=0;
  %
  budgO.fluxes.trWtop=trWtop;%kg/s
  budgO.fluxes.trWbot=trWbot;%kg/s
else;
  budgO.fluxes.trWtop=-mygrid.RAC.*budgO.zconv;
  budgO.fluxes.trWbot=mygrid.RAC*0;%kg/s
  budgO.fluxes.diaWtop=-mygrid.RAC.*budgO.zdia;
  budgO.fluxes.diaWbot=mygrid.RAC*0;%kg/s
end;
budgI.fluxes.trWtop=0*mygrid.RAC;
budgI.fluxes.trWbot=budgO.fluxes.trWtop(:,:,1);%kg/s
%
budgO.zconv=mk3D(mygrid.RAC,budgO.zconv).*budgO.zconv;%Watt
budgI.zconv=mygrid.RAC.*budgI.zconv;%Watt
%
budgOI.zconv=budgO.zconv+budgI.zconv;

%horizontal divergence (advection and diffusion)
tmpUo=myparms.rhoconst*(ADVx_SLT+DFxE_SLT); 
tmpVo=myparms.rhoconst*(ADVy_SLT+DFyE_SLT);
budgO.hconv=calc_UV_conv(nansum(tmpUo,3),nansum(tmpVo,3));
%
tmpUoD=myparms.rhoconst*(DFxE_SLT);
tmpVoD=myparms.rhoconst*(DFyE_SLT);
budgO.hdia=calc_UV_conv(nansum(tmpUoD,3),nansum(tmpVoD,3));
%
tmpUi=myparms.SIsal0*(myparms.rhoi*DFxEHEFF+myparms.rhoi*ADVxHEFF);
tmpVi=myparms.SIsal0*(myparms.rhoi*DFyEHEFF+myparms.rhoi*ADVyHEFF);
budgI.hconv=calc_UV_conv(tmpUi,tmpVi); %no dh needed here
budgOI.hconv=budgO.hconv+budgI.hconv;
%
budgO.fluxes.trU=tmpUo; budgO.fluxes.trV=tmpVo;%g/s
budgO.fluxes.diaU=tmpUoD; budgO.fluxes.diaV=tmpVoD;%g/s
budgI.fluxes.trU=tmpUi; budgI.fluxes.trV=tmpVi;%g/s

1	gforget	1.3	function [budgO,budgI,budgOI]=calc_budget_salt(kBudget);
2	gforget	1.1	% CALC_BUDGET_SALT(kBudget)
3
4			gcmfaces_global;
5
6			%get variables from caller routine:
7			%----------------------------------
8
9			global myparms;
10
11			list_variables={'SALT','AB_gS','SRELAX','SIheff',...
12			'SFLUX','oceSPflx','oceSflux','WSLTMASS',...
13			'ADVx_SLT','DFxE_SLT','ADVy_SLT','DFyE_SLT',...
14			'ADVxHEFF','ADVxSNOW','DFxEHEFF','DFxESNOW',...
15			'ADVyHEFF','ADVySNOW','DFyEHEFF','DFyESNOW'};
16
17			for vv=1:length(list_variables);
18			v = evalin('caller',list_variables{vv});
19			eval([list_variables{vv} '=v;']);
20			end;
21			clear v;
22
23			test3d=length(size(ADVx_SLT{1}))>2;
24
25	gforget	1.8	if test3d\|kBudget>1;
26			list_variables={'oceSPtnd','ADVr_SLT','DFrE_SLT',...
27			'DFrI_SLT','ADVr_SLT','DFrE_SLT','DFrI_SLT'};
28			for vv=1:length(list_variables);
29			v = evalin('caller',list_variables{vv});
30			eval([list_variables{vv} '=v;']);
31			end;
32			clear v;
33			end;
34
35	gforget	1.1	%compute mapped budget:
36			%----------------------
37
38	gforget	1.5	budgO.tend=myparms.rhoconstSALT-myparms.rhoconstAB_gS;
39			budgI.tend=myparms.SIsal0myparms.rhoiSIheff;
40	gforget	1.1	%
41	gforget	1.8	tmptend=mk3D(mygrid.RAC,budgO.tend).*budgO.tend;%g/s
42			budgO.fluxes.tend=tmptend;
43			budgO.tend=nansum(tmptend,3);
44	gforget	1.5	budgI.tend=mygrid.RAC.*budgI.tend;%g/s
45	gforget	1.1	%
46	gforget	1.5	budgOI.tend=budgO.tend+budgI.tend;
47	gforget	1.1
48			%vertical divergence (air-sea fluxes or vertical adv/dif)
49	gforget	1.5	budgO.zconv=SFLUX+oceSPflx;
50			budgI.zconv=-budgO.zconv+SRELAX;
51	gforget	1.1	%in linear surface we omit :
52	gforget	1.5	if ~myparms.useNLFS; budgO.zconv=budgO.zconv-myparms.rhoconst*WSLTMASS; end;
53	gforget	1.1	%working approach for real fresh water (?) and virtual salt flux
54	gforget	1.5	if ~myparms.useRFWF\|~myparms.useNLFS; budgI.zconv=-oceSflux; end;
55	gforget	1.6	%
56			budgO.zdia=budgO.zconv;
57	gforget	1.1	%for deep ocean layer :
58			if kBudget>1;
59	gforget	1.5	budgO.zconv=-(ADVr_SLT+DFrE_SLT+DFrI_SLT)./mygrid.RAC*myparms.rhoconst;
60			budgO.zconv=budgO.zconv+oceSPtnd;%.*msk;
61	gforget	1.6	budgO.zdia=-(DFrE_SLT+DFrI_SLT)./mygrid.RAC*myparms.rhoconst;
62			budgO.zdia=budgO.zdia+oceSPtnd;%.*msk;
63	gforget	1.1	end;
64			%
65			if test3d;
66			nr=length(mygrid.RC);
67			trWtop=-(ADVr_SLT+DFrE_SLT+DFrI_SLT)*myparms.rhoconst;
68			tmp1=mk3D(oceSPflx,oceSPtnd)-cumsum(oceSPtnd,3);
69			tmp1=tmp1.*mk3D(mygrid.RAC,tmp1);
70			trWtop(:,:,2:nr)=trWtop(:,:,2:nr)+tmp1(:,:,1:nr-1);
71			%
72	gforget	1.5	trWtop(:,:,1)=budgO.zconv.*mygrid.RAC;
73	gforget	1.1	trWbot=trWtop(:,:,2:length(mygrid.RC));
74			trWbot(:,:,length(mygrid.RC))=0;
75			%
76	gforget	1.5	budgO.fluxes.trWtop=trWtop;%kg/s
77			budgO.fluxes.trWbot=trWbot;%kg/s
78	gforget	1.1	else;
79	gforget	1.7	budgO.fluxes.trWtop=-mygrid.RAC.*budgO.zconv;
80	gforget	1.5	budgO.fluxes.trWbot=mygrid.RAC*0;%kg/s
81	gforget	1.7	budgO.fluxes.diaWtop=-mygrid.RAC.*budgO.zdia;
82	gforget	1.6	budgO.fluxes.diaWbot=mygrid.RAC*0;%kg/s
83	gforget	1.1	end;
84	gforget	1.5	budgI.fluxes.trWtop=0*mygrid.RAC;
85			budgI.fluxes.trWbot=budgO.fluxes.trWtop(:,:,1);%kg/s
86			%
87			budgO.zconv=mk3D(mygrid.RAC,budgO.zconv).*budgO.zconv;%Watt
88			budgI.zconv=mygrid.RAC.*budgI.zconv;%Watt
89	gforget	1.8	%
90	gforget	1.5	budgOI.zconv=budgO.zconv+budgI.zconv;
91	gforget	1.1
92			%horizontal divergence (advection and diffusion)
93			tmpUo=myparms.rhoconst*(ADVx_SLT+DFxE_SLT);
94			tmpVo=myparms.rhoconst*(ADVy_SLT+DFyE_SLT);
95	gforget	1.5	budgO.hconv=calc_UV_conv(nansum(tmpUo,3),nansum(tmpVo,3));
96	gforget	1.6	%
97			tmpUoD=myparms.rhoconst*(DFxE_SLT);
98			tmpVoD=myparms.rhoconst*(DFyE_SLT);
99			budgO.hdia=calc_UV_conv(nansum(tmpUoD,3),nansum(tmpVoD,3));
100			%
101	gforget	1.1	tmpUi=myparms.SIsal0(myparms.rhoiDFxEHEFF+myparms.rhoi*ADVxHEFF);
102			tmpVi=myparms.SIsal0(myparms.rhoiDFyEHEFF+myparms.rhoi*ADVyHEFF);
103	gforget	1.5	budgI.hconv=calc_UV_conv(tmpUi,tmpVi); %no dh needed here
104			budgOI.hconv=budgO.hconv+budgI.hconv;
105	gforget	1.1	%
106	gforget	1.5	budgO.fluxes.trU=tmpUo; budgO.fluxes.trV=tmpVo;%g/s
107	gforget	1.6	budgO.fluxes.diaU=tmpUoD; budgO.fluxes.diaV=tmpVoD;%g/s
108	gforget	1.5	budgI.fluxes.trU=tmpUi; budgI.fluxes.trV=tmpVi;%g/s
109	gforget	1.1