1 |
gforget |
1.1 |
function []=example_budget(choiceV3orV4); |
2 |
|
|
%object: illustrates budget analysis |
3 |
|
|
%inputs: choiceV3orV4 ('v3' or 'v4') selects the sample GRID |
4 |
|
|
|
5 |
|
|
gcmfaces_global; |
6 |
|
|
if myenv.verbose>0; |
7 |
|
|
gcmfaces_msg('==============================================='); |
8 |
|
|
gcmfaces_msg(['*** entering example_budget ' ... |
9 |
|
|
'that will load a set of budget terms from file ' ... |
10 |
|
|
'and compute a series of derived diagnostics '],''); |
11 |
|
|
end; |
12 |
|
|
|
13 |
|
|
%%%%%%%%%%%%%%%%% |
14 |
|
|
%load parameters: |
15 |
|
|
%%%%%%%%%%%%%%%%% |
16 |
|
|
|
17 |
|
|
if myenv.verbose>0; |
18 |
|
|
gcmfaces_msg('* set grid files path and format, and number of faces for this grid'); |
19 |
|
|
end; |
20 |
|
|
myenv.nctiles=1; |
21 |
|
|
if strcmp(choiceV3orV4,'v4'); |
22 |
|
|
nF=5; |
23 |
|
|
fileFormat='compact'; |
24 |
|
|
dir0=fullfile([myenv.gcmfaces_dir '../']); |
25 |
|
|
dirGrid=fullfile(dir0,'GRID/'); |
26 |
|
|
myenv.nctilesdir=fullfile(dir0,'release1/nctiles_budget/'); |
27 |
|
|
else; |
28 |
|
|
error('files are missing'); |
29 |
|
|
nF=1; |
30 |
|
|
dir0=fullfile(myenv.gcmfaces_dir,'sample_input/'); |
31 |
|
|
fileFormat='straight'; |
32 |
|
|
dirGrid=fullfile(dir0,'grid_occa/'); |
33 |
|
|
myenv.nctilesdir=fullfile(dir0,'nctiles_occa/'); |
34 |
|
|
end; |
35 |
|
|
if myenv.verbose>0; |
36 |
|
|
gcmfaces_msg('* call grid_load : load grid to memory (mygrid) according to'); |
37 |
|
|
gcmfaces_msg(['dirGrid = ' dirGrid],' '); |
38 |
|
|
gcmfaces_msg(['nFaces = ' num2str(nF)],' '); |
39 |
|
|
gcmfaces_msg(['fileFormat = ' fileFormat],' '); |
40 |
|
|
% fprintf([' > dirGrid = ' dirGrid '\n']); |
41 |
|
|
end; |
42 |
|
|
|
43 |
|
|
%%%%%%%%%%%%%%%%%% |
44 |
|
|
%main computation: |
45 |
|
|
%%%%%%%%%%%%%%%%%% |
46 |
|
|
|
47 |
|
|
%select budget of interest |
48 |
|
|
nameBudg='budgMo'; |
49 |
|
|
|
50 |
|
|
%load budget terms |
51 |
|
|
fileName=[myenv.nctilesdir nameBudg filesep nameBudg(1:6)]; |
52 |
|
|
tend=read_nctiles(fileName,'tend'); |
53 |
|
|
trU=read_nctiles(fileName,'trU'); |
54 |
|
|
trV=read_nctiles(fileName,'trV'); |
55 |
|
|
trWtop=read_nctiles(fileName,'trWtop'); |
56 |
|
|
trWbot=read_nctiles(fileName,'trWbot'); |
57 |
|
|
|
58 |
|
|
%load dt (time increments) vector (not a gcmfaces object) |
59 |
|
|
ncload([fileName '.0001.nc'],'dt'); |
60 |
|
|
|
61 |
|
|
%get budget descitption and units |
62 |
|
|
ncid = netcdf.open([fileName '.0001.nc'],'NC_NOWRITE'); |
63 |
|
|
descr = netcdf.getAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'description'); |
64 |
|
|
varid = netcdf.inqVarID(ncid,'tend'); |
65 |
|
|
units = netcdf.getAtt(ncid,varid,'units'); |
66 |
|
|
netcdf.close(ncid); |
67 |
|
|
|
68 |
|
|
%define northern hemisphere as domain of integration |
69 |
|
|
nameMask='Northern Hemisphere'; |
70 |
|
|
mask=mygrid.mskC(:,:,1).*(mygrid.YC>0); |
71 |
|
|
areaMask=mygrid.RAC.*mask; |
72 |
|
|
|
73 |
|
|
%edit plot title accordingly |
74 |
|
|
tmp1=strfind(descr,'-- ECCO v4'); |
75 |
|
|
descr=[descr(1:tmp1-1) 'for: ' nameMask]; |
76 |
|
|
|
77 |
|
|
%compute northern hemisphere integrals |
78 |
|
|
budg.tend=NaN*dt; |
79 |
|
|
budg.hconv=NaN*dt; |
80 |
|
|
budg.zconv=NaN*dt; |
81 |
|
|
for tt=1:length(dt); |
82 |
|
|
%compute flux convergence |
83 |
|
|
hconv=calc_UV_conv(trU(:,:,tt),trV(:,:,tt)); |
84 |
|
|
zconv=trWtop(:,:,tt)-trWbot(:,:,tt); |
85 |
|
|
%compute sum over domain |
86 |
|
|
budg.tend(tt)=nansum(tend(:,:,tt).*mask)/nansum(areaMask); |
87 |
|
|
budg.hconv(tt)=nansum(hconv.*mask)/nansum(areaMask); |
88 |
|
|
budg.zconv(tt)=nansum(zconv.*mask)/nansum(areaMask); |
89 |
|
|
end; |
90 |
|
|
|
91 |
|
|
%display result |
92 |
|
|
figureL; |
93 |
|
|
subplot(3,1,1); set(gca,'FontSize',12); |
94 |
|
|
plot(cumsum(dt.*budg.tend)); |
95 |
|
|
grid on; xlabel('month'); ylabel([units '.s']); |
96 |
|
|
legend('content anomaly'); title(descr); |
97 |
|
|
subplot(3,1,2); set(gca,'FontSize',12); |
98 |
|
|
plot(cumsum(dt.*budg.tend)); hold on; |
99 |
|
|
plot(cumsum(dt.*budg.hconv),'r'); |
100 |
|
|
plot(cumsum(dt.*budg.zconv),'g'); |
101 |
|
|
grid on; xlabel('month'); ylabel([units '.s']); |
102 |
|
|
legend('content anomaly','horizontal convergence','vertical convergence'); |
103 |
|
|
subplot(3,1,3); set(gca,'FontSize',12); |
104 |
|
|
plot(cumsum(dt.*(budg.tend-budg.hconv-budg.zconv))); |
105 |
|
|
grid on; xlabel('month'); ylabel([units '.s']); |
106 |
|
|
legend('budget residual'); |
107 |
|
|
|