matlab_class/gcmfaces_calc/calc_budget_heat.m

function [budgO,budgI,budgOI]=calc_budget_heat(kBudget);
% CALC_BUDGET_HEAT(kBudget)

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'};

if kBudget>1; 
  list_variables={list_variables{:},'oceQsw','ADVr_TH','DFrE_TH',...
                  'DFrI_TH','ADVr_TH','DFrE_TH','DFrI_TH'};
end;

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;

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

contOCN=myparms.rcp*THETA-myparms.rcp*AB_gT;
contICE=-myparms.flami*(SIheff*myparms.rhoi+SIhsnow*myparms.rhosn);
%
budgO.tend=mk3D(mygrid.RAC,contOCN).*contOCN;%Watt
budgI.tend=mygrid.RAC.*contICE;%Watt
%
contOCN=nansum(contOCN,3);
contTOT=contOCN+contICE;

%vertical divergence (air-sea fluxes or vertical adv/dif)
zdivOCN=TFLUX+geothFlux;
zdivICE=-(SItflux+TFLUX-TRELAX);
%in linear surface we omit :
if ~myparms.useNLFS; zdivOCN=zdivOCN-myparms.rcp*WTHMASS; end;
%in virtual salt flux we omit :
if ~myparms.useRFWF|~myparms.useNLFS; zdivICE=zdivICE+SIaaflux; end;
%working approach for real fresh water (?) and virtual salt flux
if 0; zdivICE=-oceQnet-SIatmQnt-myparms.flami*(SIsnPrcp-SIacSubl); end;
%for deep ocean layer :
if kBudget>1;
zdivOCN=-(ADVr_TH+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;
zdivOCN=zdivOCN+swfrac*oceQsw+geothFlux;%.*msk;
end;
%
zdivTOT=zdivOCN+zdivICE;
%
%note: geothFlux remains to be accounted for in the following
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)=zdivOCN.*mygrid.RAC;
  trWbot=trWtop(:,:,2:length(mygrid.RC));
  trWbot(:,:,length(mygrid.RC))=0;
  %
  budgO.trWtop=trWtop;%Watt
  budgO.trWbot=trWbot;%Watt
else;
  budgO.trWtop=mygrid.RAC.*zdivOCN; budgO.trWbot=mygrid.RAC*0;%Watt
end;
budgI.trWtop=mygrid.RAC.*(zdivICE+zdivOCN); budgI.trWbot=mygrid.RAC.*zdivOCN;%Watt

%horizontal divergence (advection and diffusion)
tmpUo=myparms.rcp*(ADVx_TH+DFxE_TH); tmpVo=myparms.rcp*(ADVy_TH+DFyE_TH);
hdivOCN=calc_UV_conv(nansum(tmpUo,3),nansum(tmpVo,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);
hdivICE=calc_UV_conv(tmpUi,tmpVi); %no dh needed here
hdivTOT=hdivOCN+hdivICE;
%
budgO.trU=tmpUo; budgO.trV=tmpVo;%Watt
budgI.trU=tmpUi; budgI.trV=tmpVi;%Watt

%outputs:
if test3d;
  tmp1=mk3D(mygrid.DRF,mygrid.hFacC).*mygrid.hFacC;
  tmp1=mk3D(mygrid.RAC,mygrid.hFacC).*tmp1;
else;
  tmp1=mygrid.RAC;
end;

budgO.cont=tmp1.*contOCN;
budgO.hconv=hdivOCN;
budgO.zconv=tmp1.*zdivOCN;

budgI.cont=tmp1.*contICE;
budgI.hconv=hdivICE;
budgI.zconv=tmp1.*zdivICE;

budgOI.cont=budgO.cont+budgI.cont;
budgOI.hconv=budgO.hconv+budgI.hconv;
budgOI.zconv=budgO.zconv+budgI.zconv;

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