matlab_class/gcmfaces_calc/calc_boxmean_T.m

function [varargout]=calc_boxmean_T(fld,varargin);
%purpose: compute a weighted average of a field
%
%inputs:
%       (1) fld is a 3D (or 2D) field in gcmfaces format at the grid 
%               center where missing values are assigned a NaN  
%       (2a) 'LATS',LATS is a line vector of latitutde interval edges
%       (2b) 'LONS',LONS is a line vector of longitutde interval edges
%       (3a) 'weights',weights is a weights field. If it is not specified
%               then we use the grid cell volume in the weighted average.
%               If specifying 'weights' you likely want to set it to 0 on land.
%       (3b) 'level',kk is a specification of the vertical level. This
%               is only used when fld is a 2D field to define weights
%outputs:
%       (1) FLD is the field/vector/value of the weighted mean. It's 
%               size depends on the input size
%       (2) W is the weights for each element of FLD. It can be used to simply
%               compute the global weighted average as nansum(W(:).*FLD(:))/nansum(W(:))
%usage: 
%       input 2a and/or 2b is necessary
%       input 3a or 3b is optional, and are mutually exclusive, except if 
%               fld is 2D when input 3b is needed and 3a is excluded 
%       output 2 is optional

%by assumption: grid_load is done
global mygrid;

%get arguments
for ii=1:(nargin-1)/2;
        tmp1=varargin{2*ii-1}; eval([tmp1 '=varargin{2*ii};']);
end;

%initialize output:
n3=max(size(fld.f1,3),1); n4=max(size(fld.f1,4),1);

%test inputs/outputs consistency:
tmp1=~isempty(whos('LONS'))+~isempty(whos('LATS'));
if tmp1~=1&tmp1~=2; error('wrong input specification'); end;
tmp1=~isempty(whos('weights'))+~isempty(whos('level'));
if tmp1>1; error('wrong input specification'); end;
if nargout>2; error('wrong output specification'); end
if n3==1&isempty(whos('level'))&isempty(whos('weights')); error('wrong input specification'); end;

%prepare output etc.:
if ~isempty(whos('LONS'))&~isempty(whos('LATS')); %use complex
valranges=LONS'*ones(size(LATS))+i*(ones(size(LONS))'*LATS); 
elseif ~isempty(whos('LONS'));
valranges=LONS'; 
else;
valranges=i*LATS;
end;
nnranges=size(valranges);
nnout=max(size(valranges)-1,[1 1]);
%
FLD=NaN*ones([nnout n3 n4]);
W=NaN*ones([nnout n3 n4]);

%select weights for average:
if isempty(whos('weights'))&isempty(whos('level'));
  weights=mygrid.hFacC.*mk3D(mygrid.RAC,mygrid.hFacC);
  weights=weights.*mk3D(mygrid.DRF,mygrid.hFacC);
elseif ~isempty(whos('level'));
  weights=mygrid.hFacC(:,:,level).*mygrid.RAC*mygrid.DRF(level);
end;

%multiply with data mask:
msk=repmat(1*~isnan(fld),[1 1 n3 n4]); 
weights=weights.*msk;
%switch to 2D array to speed up computation:
fld=convert2array(fld);
n1=size(fld,1); n2=size(fld,2);
fld=reshape(fld,n1*n2,n3*n4);
%same for the weights:
weights=reshape(convert2array(weights),n1*n2*n3,1)*ones(1,n4);
weights=reshape(weights,n1*n2,n3*n4);
%multiply one with the other
fld=fld.*weights;
%remove data mask
fld(isnan(fld))=0;

lonvec=reshape(convert2array(mygrid.XC),n1*n2,1);
latvec=reshape(convert2array(mygrid.YC),n1*n2,1);

for ix=1:nnout(1);
for iy=1:nnout(2);

   %get list ofpoints that form a zonal band:
   if nnout(1)==1;
     mm=find(latvec>=imag(valranges(ix,iy))&latvec<imag(valranges(ix,iy+1)));
   elseif nnout(2)==1;
     mm=find(lonvec>=real(valranges(ix,iy))&lonvec<real(valranges(ix+1,iy)));
   else;
     mm=find(latvec>=imag(valranges(ix,iy))&latvec<imag(valranges(ix,iy+1))...
            &lonvec>=real(valranges(ix,iy))&lonvec<real(valranges(ix+1,iy)));
   end;

   %do the area weighed average along this band: 
   tmp1=sum(fld(mm,:),1); 
   tmp2=sum(weights(mm,:),1); 
   tmp2(tmp2==0)=NaN;
   tmp1=tmp1./tmp2;

   %store:
   FLD(ix,iy,:,:)=reshape(tmp1,n3,n4);
   W(ix,iy,:,:)=reshape(tmp2,n3,n4);

end; 
end;

if nargout==2;
  varargout={FLD,W};
elseif nargout==1;
  varargout={FLD};
else;
  varargout={};
end;


1	gforget	1.1	function [varargout]=calc_boxmean_T(fld,varargin);
2			%purpose: compute a weighted average of a field
3			%
4			%inputs:
5			% (1) fld is a 3D (or 2D) field in gcmfaces format at the grid
6			% center where missing values are assigned a NaN
7			% (2a) 'LATS',LATS is a line vector of latitutde interval edges
8			% (2b) 'LONS',LONS is a line vector of longitutde interval edges
9			% (3a) 'weights',weights is a weights field. If it is not specified
10			% then we use the grid cell volume in the weighted average.
11	gforget	1.2	% If specifying 'weights' you likely want to set it to 0 on land.
12	gforget	1.1	% (3b) 'level',kk is a specification of the vertical level. This
13			% is only used when fld is a 2D field to define weights
14			%outputs:
15			% (1) FLD is the field/vector/value of the weighted mean. It's
16			% size depends on the input size
17			% (2) W is the weights for each element of FLD. It can be used to simply
18			% compute the global weighted average as nansum(W(:).*FLD(:))/nansum(W(:))
19			%usage:
20			% input 2a and/or 2b is necessary
21			% input 3a or 3b is optional, and are mutually exclusive, except if
22			% fld is 2D when input 3b is needed and 3a is excluded
23			% output 2 is optional
24
25			%by assumption: grid_load is done
26			global mygrid;
27
28			%get arguments
29			for ii=1:(nargin-1)/2;
30			tmp1=varargin{2ii-1}; eval([tmp1 '=varargin{2ii};']);
31			end;
32
33			%initialize output:
34			n3=max(size(fld.f1,3),1); n4=max(size(fld.f1,4),1);
35
36			%test inputs/outputs consistency:
37			tmp1=~isempty(whos('LONS'))+~isempty(whos('LATS'));
38			if tmp1~=1&tmp1~=2; error('wrong input specification'); end;
39			tmp1=~isempty(whos('weights'))+~isempty(whos('level'));
40			if tmp1>1; error('wrong input specification'); end;
41			if nargout>2; error('wrong output specification'); end
42	gforget	1.2	if n3==1&isempty(whos('level'))&isempty(whos('weights')); error('wrong input specification'); end;
43	gforget	1.1
44			%prepare output etc.:
45			if ~isempty(whos('LONS'))&~isempty(whos('LATS')); %use complex
46			valranges=LONS'ones(size(LATS))+i(ones(size(LONS))'*LATS);
47			elseif ~isempty(whos('LONS'));
48			valranges=LONS';
49			else;
50			valranges=i*LATS;
51			end;
52			nnranges=size(valranges);
53			nnout=max(size(valranges)-1,[1 1]);
54			%
55			FLD=NaN*ones([nnout n3 n4]);
56			W=NaN*ones([nnout n3 n4]);
57
58			%select weights for average:
59			if isempty(whos('weights'))&isempty(whos('level'));
60			weights=mygrid.hFacC.*mk3D(mygrid.RAC,mygrid.hFacC);
61			weights=weights.*mk3D(mygrid.DRF,mygrid.hFacC);
62			elseif ~isempty(whos('level'));
63			weights=mygrid.hFacC(:,:,level).mygrid.RACmygrid.DRF(level);
64			end;
65	gforget	1.2
66			%multiply with data mask:
67			msk=repmat(1*~isnan(fld),[1 1 n3 n4]);
68			weights=weights.*msk;
69			%switch to 2D array to speed up computation:
70			fld=convert2array(fld);
71			n1=size(fld,1); n2=size(fld,2);
72			fld=reshape(fld,n1n2,n3n4);
73			%same for the weights:
74	gforget	1.1	weights=reshape(convert2array(weights),n1n2n3,1)*ones(1,n4);
75			weights=reshape(weights,n1n2,n3n4);
76	gforget	1.2	%multiply one with the other
77	gforget	1.1	fld=fld.*weights;
78	gforget	1.2	%remove data mask
79			fld(isnan(fld))=0;
80	gforget	1.1
81			lonvec=reshape(convert2array(mygrid.XC),n1*n2,1);
82			latvec=reshape(convert2array(mygrid.YC),n1*n2,1);
83
84			for ix=1:nnout(1);
85			for iy=1:nnout(2);
86
87			%get list ofpoints that form a zonal band:
88			if nnout(1)==1;
89			mm=find(latvec>=imag(valranges(ix,iy))&latvec<imag(valranges(ix,iy+1)));
90			elseif nnout(2)==1;
91			mm=find(lonvec>=real(valranges(ix,iy))&lonvec<real(valranges(ix+1,iy)));
92			else;
93			mm=find(latvec>=imag(valranges(ix,iy))&latvec<imag(valranges(ix,iy+1))...
94			&lonvec>=real(valranges(ix,iy))&lonvec<real(valranges(ix+1,iy)));
95			end;
96
97			%do the area weighed average along this band:
98			tmp1=sum(fld(mm,:),1);
99			tmp2=sum(weights(mm,:),1);
100			tmp2(tmp2==0)=NaN;
101			tmp1=tmp1./tmp2;
102
103			%store:
104			FLD(ix,iy,:,:)=reshape(tmp1,n3,n4);
105			W(ix,iy,:,:)=reshape(tmp2,n3,n4);
106
107			end;
108			end;
109
110			if nargout==2;
111			varargout={FLD,W};
112			elseif nargout==1;
113			varargout={FLD};
114			else;
115			varargout={};
116			end;
117
118
119