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

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

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


Revision 1.6 - (hide annotations) (download)
Tue Feb 9 01:34:47 2016 UTC (9 years, 5 months ago) by gforget
Branch: MAIN
Changes since 1.5: +18 -6 lines
- add diabatic flux and convergence computations

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 gforget 1.5 budgO.tend=myparms.rcp*THETA-myparms.rcp*AB_gT;
35     budgI.tend=-myparms.flami*(SIheff*myparms.rhoi+SIhsnow*myparms.rhosn);
36 gforget 1.1 %
37 gforget 1.5 budgO.tend=mk3D(mygrid.RAC,budgO.tend).*budgO.tend;%Watt
38     budgI.tend=mygrid.RAC.*budgI.tend;%Watt
39 gforget 1.1 %
40 gforget 1.5 budgO.tend=nansum(budgO.tend,3);
41     budgOI.tend=budgO.tend+budgI.tend;
42 gforget 1.1
43     %vertical divergence (air-sea fluxes or vertical adv/dif)
44 gforget 1.5 budgO.zconv=TFLUX+geothFlux;
45     budgI.zconv=-(SItflux+TFLUX-TRELAX);
46 gforget 1.1 %in linear surface we omit :
47 gforget 1.5 if ~myparms.useNLFS; budgO.zconv=budgO.zconv-myparms.rcp*WTHMASS; end;
48 gforget 1.1 %in virtual salt flux we omit :
49 gforget 1.5 if ~myparms.useRFWF|~myparms.useNLFS; budgI.zconv=budgI.zconv+SIaaflux; end;
50 gforget 1.1 %working approach for real fresh water (?) and virtual salt flux
51 gforget 1.5 if 0; budgI.zconv=-oceQnet-SIatmQnt-myparms.flami*(SIsnPrcp-SIacSubl); end;
52 gforget 1.6 %
53     budgO.zdia=budgO.zconv;
54 gforget 1.1 %for deep ocean layer :
55     if kBudget>1;
56 gforget 1.6 budgO.zconv=-(ADVr_TH+DFrE_TH+DFrI_TH)./mygrid.RAC*myparms.rcp;
57     budgO.zdia=-(DFrE_TH+DFrI_TH)./mygrid.RAC*myparms.rcp;
58     dd=mygrid.RF(kBudget); msk=mygrid.mskC(:,:,kBudget);
59     swfrac=0.62*exp(dd/0.6)+(1-0.62)*exp(dd/20);
60     if dd<-200; swfrac=0; end;
61     budgO.zconv=budgO.zconv+swfrac*oceQsw+geothFlux;%.*msk;
62     budgO.zdia=budgO.zdia+swfrac*oceQsw+geothFlux;%.*msk;
63 gforget 1.1 end;
64     %
65 gforget 1.6 %notes: - geothFlux remains to be accounted for in 3D case
66     % - diaWtop, diaWbot remain to be implemented in 3D case
67 gforget 1.1 if test3d;
68     trWtop=-(ADVr_TH+DFrE_TH+DFrI_TH)*myparms.rcp;
69     %
70     dd=mygrid.RF(1:end-1);
71     swfrac=0.62*exp(dd/0.6)+(1-0.62)*exp(dd/20);
72     swfrac(dd<-200)=0;
73     swtop=mk3D(swfrac,trWtop).*mk3D(mygrid.RAC.*oceQsw,trWtop);
74     swtop(isnan(mygrid.mskC))=0;
75     trWtop=trWtop+swtop;
76     %
77 gforget 1.5 trWtop(:,:,1)=budgO.zconv.*mygrid.RAC;
78 gforget 1.1 trWbot=trWtop(:,:,2:length(mygrid.RC));
79     trWbot(:,:,length(mygrid.RC))=0;
80     %
81 gforget 1.5 budgO.fluxes.trWtop=trWtop;%Watt
82     budgO.fluxes.trWbot=trWbot;%Watt
83 gforget 1.1 else;
84 gforget 1.5 budgO.fluxes.trWtop=mygrid.RAC.*(budgO.zconv-geothFlux);
85     budgO.fluxes.trWbot=mygrid.RAC.*geothFlux;%Watt
86 gforget 1.6 budgO.fluxes.diaWtop=mygrid.RAC.*(budgO.zdia-geothFlux);
87     budgO.fluxes.diaWbot=mygrid.RAC.*geothFlux;%Watt
88 gforget 1.1 end;
89 gforget 1.5 budgI.fluxes.trWtop=mygrid.RAC.*(budgI.zconv+budgO.zconv);
90     budgI.fluxes.trWbot=mygrid.RAC.*budgO.zconv;%Watt
91     %
92     budgO.zconv=mk3D(mygrid.RAC,budgO.zconv).*budgO.zconv;%Watt
93     budgI.zconv=mygrid.RAC.*budgI.zconv;%Watt
94     budgOI.zconv=budgO.zconv+budgI.zconv;
95 gforget 1.1
96     %horizontal divergence (advection and diffusion)
97     tmpUo=myparms.rcp*(ADVx_TH+DFxE_TH); tmpVo=myparms.rcp*(ADVy_TH+DFyE_TH);
98 gforget 1.5 budgO.hconv=calc_UV_conv(nansum(tmpUo,3),nansum(tmpVo,3));
99 gforget 1.6 %
100     tmpUoD=myparms.rcp*DFxE_TH; tmpVoD=myparms.rcp*DFyE_TH;
101     budgO.hdia=calc_UV_conv(nansum(tmpUoD,3),nansum(tmpVoD,3));
102     %
103 gforget 1.1 tmpUi=-myparms.flami*(myparms.rhoi*DFxEHEFF+myparms.rhosn*DFxESNOW+myparms.rhoi*ADVxHEFF+myparms.rhosn*ADVxSNOW);
104     tmpVi=-myparms.flami*(myparms.rhoi*DFyEHEFF+myparms.rhosn*DFyESNOW+myparms.rhoi*ADVyHEFF+myparms.rhosn*ADVySNOW);
105 gforget 1.5 budgI.hconv=calc_UV_conv(tmpUi,tmpVi); %no dh needed here
106     budgOI.hconv=budgO.hconv+budgI.hconv;
107 gforget 1.1 %
108 gforget 1.5 budgO.fluxes.trU=tmpUo; budgO.fluxes.trV=tmpVo;%Watt
109 gforget 1.6 budgO.fluxes.diaU=tmpUoD; budgO.fluxes.diaV=tmpVoD;%Watt
110 gforget 1.5 budgI.fluxes.trU=tmpUi; budgI.fluxes.trV=tmpVi;%Watt
111 gforget 1.1

  ViewVC Help
Powered by ViewVC 1.1.22