matlab_class/gcmfaces_misc/check_budg.m

function [prec]=check_budg(budgIn,kk);
%function [prec]=check_budg(budgIn,kk);
%
%inputs:  budgIn is a structure containting tend, trU, trV, trWtop, 
%               trWbot (calling sequence reported in check_loop.m)
%         kk is a vertical level (0 for all levels)
%outputs: prec contains 3 measures of residual errors

gcmfaces_global;

tend=budgIn.tend;
trU=budgIn.trU;
trV=budgIn.trV;
trWtop=budgIn.trWtop;
trWbot=budgIn.trWbot;

dt=1;
units='1';
%dt=budgIn.t1-budgIn.t0;
%units = budgIn.specs.units;

%define northern hemisphere as domain of integration
nameMask='Northern Hemisphere';
mask=mygrid.mskC.*mk3D(mygrid.YC>0,mygrid.mskC);
if length(size(tend{1}))<3; mask=mask(:,:,1); end;

%focus on one level
if kk>0;
  mask=mask(:,:,kk);
  tend=tend(:,:,kk);
  trU=trU(:,:,kk);
  trV=trV(:,:,kk);
  trWtop=trWtop(:,:,kk);
  trWbot=trWbot(:,:,kk);
end;

%edit plot title accordingly
descr=nameMask;

%compute northern hemisphere integrals
budg.tend=NaN*dt;
budg.hconv=NaN*dt;
budg.zconv=NaN*dt;

%compute flux convergence
hconv=calc_UV_conv(trU,trV,{});
zconv=-(trWtop-trWbot);%note the confusing upward sign convention ...

%compute local residuals
hconv(tend==0)=NaN;
zconv(tend==0)=NaN;
tend(tend==0)=NaN;

%compute local residuals
norm=sqrt(tend.^2+hconv.^2+zconv.^2);
res=abs(tend-hconv-zconv);

%compute sum over domain
budg.tend=nansum(tend.*mask);
budg.hconv=nansum(hconv.*mask);
budg.zconv=nansum(zconv.*mask);
%
budg.res=abs(budg.tend-budg.hconv-budg.zconv);
budg.norm=sqrt(budg.tend.^2+budg.hconv.^2+budg.zconv.^2);

%output result
prec(1,1)=nanmedian(res)/nanmedian(norm);
prec(1,2)=nanmax(res./norm);
prec(1,3)=budg.res/budg.norm;

1	gforget	1.3	function [prec]=check_budg(budgIn,kk);
2			%function [prec]=check_budg(budgIn,kk);
3	gforget	1.2	%
4	gforget	1.3	%inputs: budgIn is a structure containting tend, trU, trV, trWtop,
5			% trWbot (calling sequence reported in check_loop.m)
6			% kk is a vertical level (0 for all levels)
7			%outputs: prec contains 3 measures of residual errors
8	gforget	1.1
9			gcmfaces_global;
10
11			tend=budgIn.tend;
12			trU=budgIn.trU;
13			trV=budgIn.trV;
14			trWtop=budgIn.trWtop;
15			trWbot=budgIn.trWbot;
16
17			dt=1;
18	gforget	1.2	units='1';
19			%dt=budgIn.t1-budgIn.t0;
20			%units = budgIn.specs.units;
21	gforget	1.1
22			%define northern hemisphere as domain of integration
23			nameMask='Northern Hemisphere';
24			mask=mygrid.mskC.*mk3D(mygrid.YC>0,mygrid.mskC);
25			if length(size(tend{1}))<3; mask=mask(:,:,1); end;
26
27	gforget	1.2	%focus on one level
28			if kk>0;
29			mask=mask(:,:,kk);
30			tend=tend(:,:,kk);
31			trU=trU(:,:,kk);
32			trV=trV(:,:,kk);
33			trWtop=trWtop(:,:,kk);
34			trWbot=trWbot(:,:,kk);
35			end;
36
37	gforget	1.1	%edit plot title accordingly
38			descr=nameMask;
39
40			%compute northern hemisphere integrals
41			budg.tend=NaN*dt;
42			budg.hconv=NaN*dt;
43			budg.zconv=NaN*dt;
44
45			%compute flux convergence
46			hconv=calc_UV_conv(trU,trV,{});
47	gforget	1.2	zconv=-(trWtop-trWbot);%note the confusing upward sign convention ...
48
49	gforget	1.3	%compute local residuals
50			hconv(tend==0)=NaN;
51			zconv(tend==0)=NaN;
52			tend(tend==0)=NaN;
53
54			%compute local residuals
55			norm=sqrt(tend.^2+hconv.^2+zconv.^2);
56			res=abs(tend-hconv-zconv);
57
58	gforget	1.1	%compute sum over domain
59			budg.tend=nansum(tend.*mask);
60			budg.hconv=nansum(hconv.*mask);
61			budg.zconv=nansum(zconv.*mask);
62	gforget	1.3	%
63			budg.res=abs(budg.tend-budg.hconv-budg.zconv);
64			budg.norm=sqrt(budg.tend.^2+budg.hconv.^2+budg.zconv.^2);
65	gforget	1.1
66	gforget	1.3	%output result
67			prec(1,1)=nanmedian(res)/nanmedian(norm);
68			prec(1,2)=nanmax(res./norm);
69			prec(1,3)=budg.res/budg.norm;
70	gforget	1.1