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

p = genpath('gcmfaces/'); addpath(p);
p = genpath('MITprof/'); addpath(p);

gcmfaces_global;

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