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.
%       (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')); 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]);

%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);
msk=1*~isnan(fld); fld(isnan(fld))=0;
%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;
weights=reshape(convert2array(weights),n1*n2*n3,1)*ones(1,n4);
weights=reshape(weights,n1*n2,n3*n4);
%
fld=fld.*weights;

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