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
%
%note: in the nctiles_budget_2d files from 03-Feb-2015 the sign convention 
%   for trWtop and trWbot was positive downward; this sign convention was
%   reversed in later versions for consistency with MITgcm convention.

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