matlab_class/gcmfaces_calc/calc_budget_heat.m

function [budgO,budgI,budgOI]=calc_budget_heat(kBudget);
% CALC_BUDGET_HEAT(kBudget)
%
% note: within this routine `THETA', `SIheff', and `SIhsnow' denote 
%  the corresponding tendencies as computed by diags_diff_snapshots.m
%  rather than the state variables themselves.

gcmfaces_global;

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

global myparms;

list_variables={'THETA','AB_gT','TRELAX','SIheff','SIhsnow',...
                'TFLUX','geothFlux','SItflux','SIaaflux','oceQnet',...
                'SIatmQnt','SIsnPrcp','SIacSubl','WTHMASS',...
                'ADVx_TH','DFxE_TH','ADVy_TH','DFyE_TH',...
                '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_TH{1}))>2;

if test3d|kBudget>1;
  list_variables={'oceQsw','ADVr_TH','DFrE_TH',...
                  'DFrI_TH','ADVr_TH','DFrE_TH','DFrI_TH'};
  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.rcp*THETA-myparms.rcp*AB_gT;
budgI.tend=-myparms.flami*(SIheff*myparms.rhoi+SIhsnow*myparms.rhosn);
%
tmptend=mk3D(mygrid.RAC,budgO.tend).*budgO.tend;%Watt
budgO.fluxes.tend=tmptend;
budgO.tend=nansum(tmptend,3);
budgI.tend=mygrid.RAC.*budgI.tend;%Watt
%
budgOI.tend=budgO.tend+budgI.tend;

%vertical divergence (air-sea fluxes or vertical adv/dif)
budgO.zconv=TFLUX+geothFlux;
budgI.zconv=-(SItflux+TFLUX-TRELAX);
%in linear surface we omit :
if ~myparms.useNLFS; budgO.zconv=budgO.zconv-myparms.rcp*WTHMASS; end;
%in virtual salt flux we omit :
if ~myparms.useRFWF|~myparms.useNLFS; budgI.zconv=budgI.zconv+SIaaflux; end;
%working approach for real fresh water (?) and virtual salt flux
if 0; budgI.zconv=-oceQnet-SIatmQnt-myparms.flami*(SIsnPrcp-SIacSubl); end;
%
budgO.zdia=budgO.zconv;
%for deep ocean layer :
if kBudget>1;
  budgO.zconv=-(ADVr_TH+DFrE_TH+DFrI_TH)./mygrid.RAC*myparms.rcp;
  budgO.zdia=-(DFrE_TH+DFrI_TH)./mygrid.RAC*myparms.rcp;
  dd=mygrid.RF(kBudget); msk=mygrid.mskC(:,:,kBudget);
  swfrac=0.62*exp(dd/0.6)+(1-0.62)*exp(dd/20);
  if dd<-200; swfrac=0; end;
  budgO.zconv=budgO.zconv+swfrac*oceQsw+geothFlux;%.*msk;
  budgO.zdia=budgO.zdia+swfrac*oceQsw+geothFlux;%.*msk;
end;
%
%notes: - geothFlux remains to be accounted for in 3D case
%       - diaWtop, diaWbot remain to be implemented in 3D case
if test3d;
  trWtop=-(ADVr_TH+DFrE_TH+DFrI_TH)*myparms.rcp;
  %
  dd=mygrid.RF(1:end-1);
  swfrac=0.62*exp(dd/0.6)+(1-0.62)*exp(dd/20);
  swfrac(dd<-200)=0;
  swtop=mk3D(swfrac,trWtop).*mk3D(mygrid.RAC.*oceQsw,trWtop);
  swtop(isnan(mygrid.mskC))=0;
  trWtop=trWtop+swtop;
  %
  trWtop(:,:,1)=budgO.zconv.*mygrid.RAC;
  trWbot=trWtop(:,:,2:length(mygrid.RC));
  trWbot(:,:,length(mygrid.RC))=0;
  %
  budgO.fluxes.trWtop=trWtop;%Watt
  budgO.fluxes.trWbot=trWbot;%Watt
else;
  budgO.fluxes.trWtop=-mygrid.RAC.*(budgO.zconv-geothFlux);
  budgO.fluxes.trWbot=mygrid.RAC.*geothFlux;%Watt
  budgO.fluxes.diaWtop=-mygrid.RAC.*(budgO.zdia-geothFlux);
  budgO.fluxes.diaWbot=mygrid.RAC.*geothFlux;%Watt
end;
budgI.fluxes.trWtop=-mygrid.RAC.*(budgI.zconv+budgO.zconv);
budgI.fluxes.trWbot=-mygrid.RAC.*budgO.zconv;%Watt
%
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.rcp*(ADVx_TH+DFxE_TH); tmpVo=myparms.rcp*(ADVy_TH+DFyE_TH);
budgO.hconv=calc_UV_conv(nansum(tmpUo,3),nansum(tmpVo,3));
%
tmpUoD=myparms.rcp*DFxE_TH; tmpVoD=myparms.rcp*DFyE_TH;
budgO.hdia=calc_UV_conv(nansum(tmpUoD,3),nansum(tmpVoD,3));
%
tmpUi=-myparms.flami*(myparms.rhoi*DFxEHEFF+myparms.rhosn*DFxESNOW+myparms.rhoi*ADVxHEFF+myparms.rhosn*ADVxSNOW);
tmpVi=-myparms.flami*(myparms.rhoi*DFyEHEFF+myparms.rhosn*DFyESNOW+myparms.rhoi*ADVyHEFF+myparms.rhosn*ADVySNOW);
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;%Watt
budgO.fluxes.diaU=tmpUoD; budgO.fluxes.diaV=tmpVoD;%Watt
budgI.fluxes.trU=tmpUi; budgI.fluxes.trV=tmpVi;%Watt

1	gforget	1.3	function [budgO,budgI,budgOI]=calc_budget_heat(kBudget);
2	gforget	1.1	% CALC_BUDGET_HEAT(kBudget)
3	gforget	1.9	%
4			% note: within this routine `THETA', `SIheff', and `SIhsnow' denote
5			% the corresponding tendencies as computed by diags_diff_snapshots.m
6			% rather than the state variables themselves.
7	gforget	1.1
8			gcmfaces_global;
9
10			%get variables from caller routine:
11			%----------------------------------
12
13			global myparms;
14
15			list_variables={'THETA','AB_gT','TRELAX','SIheff','SIhsnow',...
16			'TFLUX','geothFlux','SItflux','SIaaflux','oceQnet',...
17			'SIatmQnt','SIsnPrcp','SIacSubl','WTHMASS',...
18			'ADVx_TH','DFxE_TH','ADVy_TH','DFyE_TH',...
19			'ADVxHEFF','ADVxSNOW','DFxEHEFF','DFxESNOW',...
20			'ADVyHEFF','ADVySNOW','DFyEHEFF','DFyESNOW'};
21
22			for vv=1:length(list_variables);
23			v = evalin('caller',list_variables{vv});
24			eval([list_variables{vv} '=v;']);
25			end;
26			clear v;
27
28			test3d=length(size(ADVx_TH{1}))>2;
29
30	gforget	1.8	if test3d\|kBudget>1;
31			list_variables={'oceQsw','ADVr_TH','DFrE_TH',...
32			'DFrI_TH','ADVr_TH','DFrE_TH','DFrI_TH'};
33			for vv=1:length(list_variables);
34			v = evalin('caller',list_variables{vv});
35			eval([list_variables{vv} '=v;']);
36			end;
37			clear v;
38			end;
39
40	gforget	1.1	%compute mapped budget:
41			%----------------------
42
43	gforget	1.5	budgO.tend=myparms.rcpTHETA-myparms.rcpAB_gT;
44			budgI.tend=-myparms.flami(SIheffmyparms.rhoi+SIhsnow*myparms.rhosn);
45	gforget	1.1	%
46	gforget	1.8	tmptend=mk3D(mygrid.RAC,budgO.tend).*budgO.tend;%Watt
47			budgO.fluxes.tend=tmptend;
48			budgO.tend=nansum(tmptend,3);
49	gforget	1.5	budgI.tend=mygrid.RAC.*budgI.tend;%Watt
50	gforget	1.1	%
51	gforget	1.5	budgOI.tend=budgO.tend+budgI.tend;
52	gforget	1.1
53			%vertical divergence (air-sea fluxes or vertical adv/dif)
54	gforget	1.5	budgO.zconv=TFLUX+geothFlux;
55			budgI.zconv=-(SItflux+TFLUX-TRELAX);
56	gforget	1.1	%in linear surface we omit :
57	gforget	1.5	if ~myparms.useNLFS; budgO.zconv=budgO.zconv-myparms.rcp*WTHMASS; end;
58	gforget	1.1	%in virtual salt flux we omit :
59	gforget	1.5	if ~myparms.useRFWF\|~myparms.useNLFS; budgI.zconv=budgI.zconv+SIaaflux; end;
60	gforget	1.1	%working approach for real fresh water (?) and virtual salt flux
61	gforget	1.5	if 0; budgI.zconv=-oceQnet-SIatmQnt-myparms.flami*(SIsnPrcp-SIacSubl); end;
62	gforget	1.6	%
63			budgO.zdia=budgO.zconv;
64	gforget	1.1	%for deep ocean layer :
65			if kBudget>1;
66	gforget	1.6	budgO.zconv=-(ADVr_TH+DFrE_TH+DFrI_TH)./mygrid.RAC*myparms.rcp;
67			budgO.zdia=-(DFrE_TH+DFrI_TH)./mygrid.RAC*myparms.rcp;
68			dd=mygrid.RF(kBudget); msk=mygrid.mskC(:,:,kBudget);
69			swfrac=0.62exp(dd/0.6)+(1-0.62)exp(dd/20);
70			if dd<-200; swfrac=0; end;
71			budgO.zconv=budgO.zconv+swfracoceQsw+geothFlux;%.msk;
72			budgO.zdia=budgO.zdia+swfracoceQsw+geothFlux;%.msk;
73	gforget	1.1	end;
74			%
75	gforget	1.6	%notes: - geothFlux remains to be accounted for in 3D case
76			% - diaWtop, diaWbot remain to be implemented in 3D case
77	gforget	1.1	if test3d;
78			trWtop=-(ADVr_TH+DFrE_TH+DFrI_TH)*myparms.rcp;
79			%
80			dd=mygrid.RF(1:end-1);
81			swfrac=0.62exp(dd/0.6)+(1-0.62)exp(dd/20);
82			swfrac(dd<-200)=0;
83			swtop=mk3D(swfrac,trWtop).mk3D(mygrid.RAC.oceQsw,trWtop);
84			swtop(isnan(mygrid.mskC))=0;
85			trWtop=trWtop+swtop;
86			%
87	gforget	1.5	trWtop(:,:,1)=budgO.zconv.*mygrid.RAC;
88	gforget	1.1	trWbot=trWtop(:,:,2:length(mygrid.RC));
89			trWbot(:,:,length(mygrid.RC))=0;
90			%
91	gforget	1.5	budgO.fluxes.trWtop=trWtop;%Watt
92			budgO.fluxes.trWbot=trWbot;%Watt
93	gforget	1.1	else;
94	gforget	1.7	budgO.fluxes.trWtop=-mygrid.RAC.*(budgO.zconv-geothFlux);
95	gforget	1.5	budgO.fluxes.trWbot=mygrid.RAC.*geothFlux;%Watt
96	gforget	1.7	budgO.fluxes.diaWtop=-mygrid.RAC.*(budgO.zdia-geothFlux);
97	gforget	1.6	budgO.fluxes.diaWbot=mygrid.RAC.*geothFlux;%Watt
98	gforget	1.1	end;
99	gforget	1.7	budgI.fluxes.trWtop=-mygrid.RAC.*(budgI.zconv+budgO.zconv);
100			budgI.fluxes.trWbot=-mygrid.RAC.*budgO.zconv;%Watt
101	gforget	1.5	%
102			budgO.zconv=mk3D(mygrid.RAC,budgO.zconv).*budgO.zconv;%Watt
103			budgI.zconv=mygrid.RAC.*budgI.zconv;%Watt
104			budgOI.zconv=budgO.zconv+budgI.zconv;
105	gforget	1.1
106			%horizontal divergence (advection and diffusion)
107			tmpUo=myparms.rcp(ADVx_TH+DFxE_TH); tmpVo=myparms.rcp(ADVy_TH+DFyE_TH);
108	gforget	1.5	budgO.hconv=calc_UV_conv(nansum(tmpUo,3),nansum(tmpVo,3));
109	gforget	1.6	%
110			tmpUoD=myparms.rcpDFxE_TH; tmpVoD=myparms.rcpDFyE_TH;
111			budgO.hdia=calc_UV_conv(nansum(tmpUoD,3),nansum(tmpVoD,3));
112			%
113	gforget	1.1	tmpUi=-myparms.flami(myparms.rhoiDFxEHEFF+myparms.rhosnDFxESNOW+myparms.rhoiADVxHEFF+myparms.rhosn*ADVxSNOW);
114			tmpVi=-myparms.flami(myparms.rhoiDFyEHEFF+myparms.rhosnDFyESNOW+myparms.rhoiADVyHEFF+myparms.rhosn*ADVySNOW);
115	gforget	1.5	budgI.hconv=calc_UV_conv(tmpUi,tmpVi); %no dh needed here
116			budgOI.hconv=budgO.hconv+budgI.hconv;
117	gforget	1.1	%
118	gforget	1.5	budgO.fluxes.trU=tmpUo; budgO.fluxes.trV=tmpVo;%Watt
119	gforget	1.6	budgO.fluxes.diaU=tmpUoD; budgO.fluxes.diaV=tmpVoD;%Watt
120	gforget	1.5	budgI.fluxes.trU=tmpUi; budgI.fluxes.trV=tmpVi;%Watt
121	gforget	1.1