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;

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:
%%%%%%%%%%%%%%%%%


if isempty(mygrid);
   grid_load;
end;

myenv.nctilesdir=fullfile('sample_input',filesep,'nctiles_budget_2d',filesep);

if ~isdir(myenv.nctilesdir);
    diags=[];
    help gcmfaces_demo;
    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:
%%%%%%%%%%%%%%%%%%

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