/[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.1 - (hide annotations) (download)
Sat Feb 6 15:31:21 2016 UTC (9 years, 5 months ago) by gforget
Branch: MAIN
- move budget computations from diags_set_D.m to separate functions.

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

  ViewVC Help
Powered by ViewVC 1.1.22