matlab_class/sample_analysis/example_budget.m

function []=example_budget();
%EXAMPLE_TRANSPORTS illustrates ECCO v4 budgets
%
%stand-alone call: addpath gcmfaces/sample_analysis/; example_budget;
%
%needed input files:
%  wget --recursive ftp://mit.ecco-group.org/gforget/nctiles_budget_2d
%  mv mit.ecco-group.org/gforget/nctiles_budget_2d sample_input/.
%  rm -rf mit.ecco-group.org

gcmfaces_global;

if myenv.verbose>0;
    gcmfaces_msg('===============================================');
    gcmfaces_msg(['*** entering example_budget ' ...
        'load budget terms from nctiles files ' ...
        'and compute hemispheric budgets'],'');
end;

%%%%%%%%%%%%%%%%%
%load grid:
%%%%%%%%%%%%%%%%%

%expected location:
myenv.nctilesdir=fullfile('sample_input',filesep,'nctiles_budget_2d',filesep);
%if nctiles_budget_2d is not found then try old location:
if ~isdir(myenv.nctilesdir);
    %if not found then try old location:
    tmpdir=fullfile(pwd,'/gcmfaces/sample_input/nctiles_budget_2d/');
    if isdir(tmpdir); myenv.nctilesdir=tmpdir; end;
end;
%if nctiles_budget_2d is still not found then issue warning and skip example_budget
if ~isdir(myenv.nctilesdir);
    diags=[];
    help example_budget;
    warning(['skipping example_budget (missing ' myenv.nctilesdir ')']);
    return;
end;

dirName=fullfile(myenv.nctilesdir,filesep,'budgMo',filesep);

if ~isdir(dirName);
    help gcmfaces_demo;
    warning(['skipping example_budget (missing ' dirName ')']);
    return;
end;

if isempty(which('ncload'));
    warning(['skipping example_budget (missing ncload that is part of MITprof)']);
    return;
end;

%%%%%%%%%%%%%%%%%%
%main computation:
%%%%%%%%%%%%%%%%%%

%load grid:
if isempty(mygrid);
   grid_load;
end;

%select budget of interest
nameBudg='budgMo';

%load budget terms
fileName=[myenv.nctilesdir nameBudg filesep nameBudg(1:6)];
tend=read_nctiles(fileName,'tend');
trU=read_nctiles(fileName,'trU');
trV=read_nctiles(fileName,'trV');
trWtop=read_nctiles(fileName,'trWtop');
trWbot=read_nctiles(fileName,'trWbot');

%load dt (time increments) vector (not a gcmfaces object)
ncload([fileName '.0001.nc'],'dt');

%get budget descitption and units
ncid = netcdf.open([fileName '.0001.nc'],'NC_NOWRITE');
descr = netcdf.getAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'description');
varid = netcdf.inqVarID(ncid,'tend');
units = netcdf.getAtt(ncid,varid,'units');
netcdf.close(ncid);

if myenv.verbose>0; gcmfaces_msg(['* display Northern Hemisphere budget for ' nameBudg ' (' descr ')']);end;

%define northern hemisphere as domain of integration
nameMask='Northern Hemisphere';
mask=mygrid.mskC(:,:,1).*(mygrid.YC>0); 
areaMask=mygrid.RAC.*mask;

%edit plot title accordingly
tmp1=strfind(descr,'-- ECCO v4');
descr=[descr(1:tmp1-1) 'for: ' nameMask];

%compute northern hemisphere integrals
budg.tend=NaN*dt;
budg.hconv=NaN*dt;
budg.zconv=NaN*dt;
for tt=1:length(dt);
    %compute flux convergence
    hconv=calc_UV_conv(trU(:,:,tt),trV(:,:,tt));
    zconv=trWtop(:,:,tt)-trWbot(:,:,tt);
    %compute sum over domain
    budg.tend(tt)=nansum(tend(:,:,tt).*mask)/nansum(areaMask);
    budg.hconv(tt)=nansum(hconv.*mask)/nansum(areaMask);
    budg.zconv(tt)=nansum(zconv.*mask)/nansum(areaMask);
end;

%display result
figureL;
subplot(3,1,1); set(gca,'FontSize',12);
plot(cumsum(dt.*budg.tend));
grid on; xlabel('month'); ylabel([units ' x s']);
legend('content anomaly');
subplot(3,1,2); set(gca,'FontSize',12);
plot(cumsum(dt.*budg.tend)); hold on;
plot(cumsum(dt.*budg.hconv),'r'); 
plot(cumsum(dt.*budg.zconv),'g');
grid on; xlabel('month'); ylabel([units ' x s']);
legend('content anomaly','horizontal convergence','vertical convergence');
subplot(3,1,3); set(gca,'FontSize',12);
plot(cumsum(dt.*(budg.tend-budg.hconv-budg.zconv))); 
grid on; xlabel('month'); ylabel([units ' x s']);
legend('budget residual');

if myenv.verbose>0;
    gcmfaces_msg('*** leaving example_budget ','');
    gcmfaces_msg('===============================================');
end;

1	gforget	1.4	function []=example_budget();
2	gforget	1.6	%EXAMPLE_TRANSPORTS illustrates ECCO v4 budgets
3			%
4			%stand-alone call: addpath gcmfaces/sample_analysis/; example_budget;
5	gforget	1.7	%
6			%needed input files:
7			% wget --recursive ftp://mit.ecco-group.org/gforget/nctiles_budget_2d
8			% mv mit.ecco-group.org/gforget/nctiles_budget_2d sample_input/.
9			% rm -rf mit.ecco-group.org
10	gforget	1.6
11	gforget	1.1	gcmfaces_global;
12	gforget	1.6
13	gforget	1.1	if myenv.verbose>0;
14			gcmfaces_msg('===============================================');
15			gcmfaces_msg(['*** entering example_budget ' ...
16	gforget	1.9	'load budget terms from nctiles files ' ...
17			'and compute hemispheric budgets'],'');
18	gforget	1.1	end;
19
20			%%%%%%%%%%%%%%%%%
21	gforget	1.6	%load grid:
22	gforget	1.1	%%%%%%%%%%%%%%%%%
23
24	gforget	1.7	%expected location:
25			myenv.nctilesdir=fullfile('sample_input',filesep,'nctiles_budget_2d',filesep);
26			%if nctiles_budget_2d is not found then try old location:
27			if ~isdir(myenv.nctilesdir);
28			%if not found then try old location:
29			tmpdir=fullfile(pwd,'/gcmfaces/sample_input/nctiles_budget_2d/');
30			if isdir(tmpdir); myenv.nctilesdir=tmpdir; end;
31	gforget	1.1	end;
32	gforget	1.7	%if nctiles_budget_2d is still not found then issue warning and skip example_budget
33	gforget	1.6	if ~isdir(myenv.nctilesdir);
34			diags=[];
35	gforget	1.7	help example_budget;
36	gforget	1.6	warning(['skipping example_budget (missing ' myenv.nctilesdir ')']);
37			return;
38	gforget	1.1	end;
39
40	gforget	1.2	dirName=fullfile(myenv.nctilesdir,filesep,'budgMo',filesep);
41	gforget	1.6
42	gforget	1.2	if ~isdir(dirName);
43	gforget	1.5	help gcmfaces_demo;
44	gforget	1.2	warning(['skipping example_budget (missing ' dirName ')']);
45			return;
46			end;
47
48	gforget	1.3	if isempty(which('ncload'));
49	gforget	1.2	warning(['skipping example_budget (missing ncload that is part of MITprof)']);
50			return;
51			end;
52
53	gforget	1.1	%%%%%%%%%%%%%%%%%%
54			%main computation:
55			%%%%%%%%%%%%%%%%%%
56
57	gforget	1.7	%load grid:
58			if isempty(mygrid);
59			grid_load;
60			end;
61
62	gforget	1.1	%select budget of interest
63			nameBudg='budgMo';
64
65			%load budget terms
66			fileName=[myenv.nctilesdir nameBudg filesep nameBudg(1:6)];
67			tend=read_nctiles(fileName,'tend');
68			trU=read_nctiles(fileName,'trU');
69			trV=read_nctiles(fileName,'trV');
70			trWtop=read_nctiles(fileName,'trWtop');
71			trWbot=read_nctiles(fileName,'trWbot');
72
73			%load dt (time increments) vector (not a gcmfaces object)
74			ncload([fileName '.0001.nc'],'dt');
75
76			%get budget descitption and units
77			ncid = netcdf.open([fileName '.0001.nc'],'NC_NOWRITE');
78			descr = netcdf.getAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'description');
79			varid = netcdf.inqVarID(ncid,'tend');
80			units = netcdf.getAtt(ncid,varid,'units');
81			netcdf.close(ncid);
82
83	gforget	1.6	if myenv.verbose>0; gcmfaces_msg(['* display Northern Hemisphere budget for ' nameBudg ' (' descr ')']);end;
84
85	gforget	1.1	%define northern hemisphere as domain of integration
86			nameMask='Northern Hemisphere';
87			mask=mygrid.mskC(:,:,1).*(mygrid.YC>0);
88			areaMask=mygrid.RAC.*mask;
89
90			%edit plot title accordingly
91			tmp1=strfind(descr,'-- ECCO v4');
92			descr=[descr(1:tmp1-1) 'for: ' nameMask];
93
94			%compute northern hemisphere integrals
95			budg.tend=NaN*dt;
96			budg.hconv=NaN*dt;
97			budg.zconv=NaN*dt;
98			for tt=1:length(dt);
99			%compute flux convergence
100			hconv=calc_UV_conv(trU(:,:,tt),trV(:,:,tt));
101			zconv=trWtop(:,:,tt)-trWbot(:,:,tt);
102			%compute sum over domain
103			budg.tend(tt)=nansum(tend(:,:,tt).*mask)/nansum(areaMask);
104			budg.hconv(tt)=nansum(hconv.*mask)/nansum(areaMask);
105			budg.zconv(tt)=nansum(zconv.*mask)/nansum(areaMask);
106			end;
107
108			%display result
109			figureL;
110			subplot(3,1,1); set(gca,'FontSize',12);
111			plot(cumsum(dt.*budg.tend));
112	gforget	1.6	grid on; xlabel('month'); ylabel([units ' x s']);
113			legend('content anomaly');
114	gforget	1.1	subplot(3,1,2); set(gca,'FontSize',12);
115			plot(cumsum(dt.*budg.tend)); hold on;
116			plot(cumsum(dt.*budg.hconv),'r');
117			plot(cumsum(dt.*budg.zconv),'g');
118	gforget	1.6	grid on; xlabel('month'); ylabel([units ' x s']);
119	gforget	1.1	legend('content anomaly','horizontal convergence','vertical convergence');
120			subplot(3,1,3); set(gca,'FontSize',12);
121			plot(cumsum(dt.*(budg.tend-budg.hconv-budg.zconv)));
122	gforget	1.6	grid on; xlabel('month'); ylabel([units ' x s']);
123	gforget	1.1	legend('budget residual');
124
125	gforget	1.6	if myenv.verbose>0;
126			gcmfaces_msg('*** leaving example_budget ','');
127			gcmfaces_msg('===============================================');
128			end;
129