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

Annotation of /MITgcm_contrib/gael/matlab_class/gcmfaces_calc/calc_budget_mass.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,budgMo,...
2     contICE,hdivICE,zdivICE,budgMi]=...
3     calc_budget_mass(kBudget);
4     % CALC_BUDGET_MASS(kBudget,doMoreBudgetOutput)
5    
6     gcmfaces_global;
7    
8     %get variables from caller routine:
9     %----------------------------------
10    
11     global myparms;
12    
13     list_variables={'ETAN','SIheff','SIhsnow','oceFWflx','SIatmFW','oceFWflx',...
14     'UVELMASS','VVELMASS',...
15     'ADVxHEFF','ADVxSNOW','DFxEHEFF','DFxESNOW',...
16     'ADVyHEFF','ADVySNOW','DFyEHEFF','DFyESNOW'};
17    
18     if kBudget>1; list_variables={list_variables{:},'WVELMASS'}; end;
19    
20     for vv=1:length(list_variables);
21     v = evalin('caller',list_variables{vv});
22     eval([list_variables{vv} '=v;']);
23     end;
24     clear v;
25    
26     test3d=length(size(UVELMASS{1}))>2;
27    
28     %compute mapped budget:
29     %----------------------
30    
31     %mass = myparms.rhoconst * sea level
32     contOCN=ETAN*myparms.rhoconst;
33     contICE=(SIheff*myparms.rhoi+SIhsnow*myparms.rhosn);
34     %for deep ocean layer :
35     if kBudget>1&myparms.useNLFS<2;
36     contOCN=0;
37     elseif kBudget>1;%rstar case
38     tmp1=mk3D(mygrid.DRF,mygrid.hFacC).*mygrid.hFacC;
39     tmp2=sum(tmp1(:,:,kBudget:length(mygrid.RC)),3)./mygrid.Depth;
40     contOCN=tmp2.*ETAN*myparms.rhoconst;
41     end;
42     %
43     contTOT=contOCN+contICE;
44     %
45     if test3d;
46     tmp1=mk3D(mygrid.DRF,mygrid.hFacC).*mygrid.hFacC;
47     tmp2=tmp1./mk3D(mygrid.Depth,tmp1);
48     tend=tmp2.*mk3D(ETAN,tmp2)*myparms.rhoconst;
49     else;
50     tend=contOCN;
51     end;
52     budgMo.tend=mk3D(mygrid.RAC,tend).*tend;%kg/s
53     budgMi.tend=mygrid.RAC.*contICE;%kg/s
54    
55     %vertical divergence (air-sea fluxes or vertical advection)
56     zdivOCN=oceFWflx;
57     zdivICE=SIatmFW-oceFWflx;
58     %in virtual salt flux we omit :
59     if ~myparms.useRFWF; zdivOCN=0*zdivOCN; end;
60     %for deep ocean layer :
61     if kBudget>1; zdivOCN=-WVELMASS*myparms.rhoconst; end;
62     %
63     zdivTOT=zdivOCN+zdivICE;
64     %
65     if test3d;
66     trWtop=-WVELMASS*myparms.rhoconst;
67     %trWtop(:,:,1)=zdivOCN;
68     trWbot=trWtop(:,:,2:length(mygrid.RC));
69     trWbot(:,:,length(mygrid.RC))=0;
70     %
71     budgMo.trWtop=mk3D(mygrid.RAC,trWtop).*trWtop;
72     budgMo.trWbot=mk3D(mygrid.RAC,trWbot).*trWbot;%kg/s
73     else;
74     budgMo.trWtop=mygrid.RAC.*zdivOCN; budgMo.trWbot=mygrid.RAC*0;%kg/s
75     end;
76     budgMi.trWtop=mygrid.RAC.*(zdivICE+zdivOCN); budgMi.trWbot=mygrid.RAC.*zdivOCN;%kg/s
77    
78     %horizontal divergence (advection and ice diffusion)
79     if test3d;
80     %3D UVELMASS,VVELMASS are multiplied by DRF
81     %(2D diagnostics are expectedly vertically integrated by MITgcm)
82     tmp1=mk3D(mygrid.DRF,UVELMASS);
83     UVELMASS=tmp1.*UVELMASS;
84     VVELMASS=tmp1.*VVELMASS;
85     end;
86     dxg=mk3D(mygrid.DXG,VVELMASS); dyg=mk3D(mygrid.DYG,UVELMASS);
87     tmpUo=myparms.rhoconst*dyg.*UVELMASS;
88     tmpVo=myparms.rhoconst*dxg.*VVELMASS;
89     hdivOCN=calc_UV_conv(nansum(tmpUo,3),nansum(tmpVo,3));
90     tmpUi=(myparms.rhoi*DFxEHEFF+myparms.rhosn*DFxESNOW+myparms.rhoi*ADVxHEFF+myparms.rhosn*ADVxSNOW);
91     tmpVi=(myparms.rhoi*DFyEHEFF+myparms.rhosn*DFyESNOW+myparms.rhoi*ADVyHEFF+myparms.rhosn*ADVySNOW);
92     hdivICE=calc_UV_conv(tmpUi,tmpVi); %dh needed is alerady in DFxEHEFF etc
93     hdivTOT=hdivOCN+hdivICE;
94     %
95     budgMo.trU=tmpUo; budgMo.trV=tmpVo;%kg/s
96     budgMi.trU=tmpUi; budgMi.trV=tmpVi;%kg/s
97    

  ViewVC Help
Powered by ViewVC 1.1.22