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=(trWbot-trWtop);

%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	function [prec]=check_budg(budgIn,kk);
2	%function [prec]=check_budg(budgIn,kk);
3	%
4	%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
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	units='1';
19	%dt=budgIn.t1-budgIn.t0;
20	%units = budgIn.specs.units;
21
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	%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	%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	zconv=(trWbot-trWtop);
48
49	%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	%compute sum over domain
59	budg.tend=nansum(tend.*mask);
60	budg.hconv=nansum(hconv.*mask);
61	budg.zconv=nansum(zconv.*mask);
62	%
63	budg.res=abs(budg.tend-budg.hconv-budg.zconv);
64	budg.norm=sqrt(budg.tend.^2+budg.hconv.^2+budg.zconv.^2);
65
66	%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