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) FLD,X,Y,Z -- FLD is the weighted mean, X,Y,Z are locations
%               that can serve for displaying the result.
%       (3) FLD,W -- FLD is the weighted mean, W are weights computed 
%               for each element of FLD. W can be used to simply compute 
%               the global weighted average as nansum(W(:).*FLD(:))/nansum(W(:))
%       (4) FLD,X,Y,Z,W
%usage: 
%       input 2a and/or 2b is necessary
%       input 3a and 3b are optional and mutually exclusive
%          if neither is specified, then it is assumed that fld
%          is a 3D field (or a series of them) and weight are
%          the usual grid cell volumes
%       output 2,3,4,5 are optional and allocated depending on call
%       by assumption, mygrid has already been loaded to memory

gcmfaces_global;

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

%initialize output:
[n1,n2,n3,n4]=size(fld.f1);

%test inputs/outputs consistency:
tmp1=isempty(whos('LONS'))&isempty(whos('LATS'));
if tmp1;
  error('wrong input specification : specify LONS and/or LATS');
end;
tmp1=~isempty(whos('weights'))+~isempty(whos('level'));
if tmp1>1;
  error('wrong input specification : omit weights or level');
end;
if n3~=length(mygrid.RC)&isempty(whos('level'))&isempty(whos('weights'));
  error('wrong input specification : specify weights or level');
end;

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

%and corresponding vertical positions:
Z=NaN; if n3>1; Z=[]; Z(1,1,:)=[1:n3]; end;
if isempty(whos('weights'))&isempty(whos('level'));
  Z(1,1,:)=mygrid.RC;
elseif ~isempty(whos('level'));
  Z(1,1,:)=mygrid.RC(level);
end;
[t1,t2,t3,t4]=size(X);
Z=repmat(Z,[t1 t2 1 t4]);

%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;

%check for potential inconsistencies in specified weights :
[w1,w2,w3,w4]=size(weights.f1);
if w3==1&n3>1; weights=repmat(weights,[1 1 n3 n4]); end;
if w4==1&n4>1; weights=repmat(weights,[1 1 1 n4]); end;
test1=sum(weights>0&isnan(fld));
if test1>0; error('non-zero weights found for NaN point'); end;

%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=convert2array(weights);
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 ~isempty(whos('LONS'))&~isempty(whos('LATS'));
     mm=find(latvec>=imag(valranges(ix,iy))&latvec<imag(valranges(ix,iy+1))...
            &lonvec>=real(valranges(ix,iy))&lonvec<real(valranges(ix+1,iy)));
   elseif ~isempty(whos('LATS'));
     mm=find(latvec>=imag(valranges(ix,iy))&latvec<imag(valranges(ix,iy+1)));
   else;
     mm=find(lonvec>=real(valranges(ix,iy))&lonvec<real(valranges(ix+1,iy)));
   end;

   %do the area weighed average along this band: 
   tmp1=nansum(fld(mm,:),1); 
   tmp2=nansum(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;

%remove singleton dimensions:
X=squeeze(X); Y=squeeze(Y); Z=squeeze(Z); 
W=squeeze(W); FLD=squeeze(FLD);

%prepare outout:
if nargout==5;
  varargout={FLD,X,Y,Z,W};
elseif nargout==4;
  varargout={FLD,X,Y,Z};
elseif 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	gforget	1.3	% (2) FLD,X,Y,Z -- FLD is the weighted mean, X,Y,Z are locations
18			% that can serve for displaying the result.
19			% (3) FLD,W -- FLD is the weighted mean, W are weights computed
20			% for each element of FLD. W can be used to simply compute
21			% the global weighted average as nansum(W(:).*FLD(:))/nansum(W(:))
22			% (4) FLD,X,Y,Z,W
23	gforget	1.1	%usage:
24			% input 2a and/or 2b is necessary
25	gforget	1.3	% input 3a and 3b are optional and mutually exclusive
26			% if neither is specified, then it is assumed that fld
27			% is a 3D field (or a series of them) and weight are
28			% the usual grid cell volumes
29			% output 2,3,4,5 are optional and allocated depending on call
30			% by assumption, mygrid has already been loaded to memory
31	gforget	1.1
32	gforget	1.3	gcmfaces_global;
33	gforget	1.1
34			%get arguments
35			for ii=1:(nargin-1)/2;
36			tmp1=varargin{2ii-1}; eval([tmp1 '=varargin{2ii};']);
37			end;
38
39			%initialize output:
40	gforget	1.3	[n1,n2,n3,n4]=size(fld.f1);
41	gforget	1.1
42			%test inputs/outputs consistency:
43	gforget	1.3	tmp1=isempty(whos('LONS'))&isempty(whos('LATS'));
44			if tmp1;
45			error('wrong input specification : specify LONS and/or LATS');
46			end;
47	gforget	1.1	tmp1=~isempty(whos('weights'))+~isempty(whos('level'));
48	gforget	1.3	if tmp1>1;
49			error('wrong input specification : omit weights or level');
50			end;
51			if n3~=length(mygrid.RC)&isempty(whos('level'))&isempty(whos('weights'));
52			error('wrong input specification : specify weights or level');
53			end;
54	gforget	1.1
55			%prepare output etc.:
56			if ~isempty(whos('LONS'))&~isempty(whos('LATS')); %use complex
57	gforget	1.3	valranges=LONS'ones(size(LATS))+i(ones(size(LONS))'*LATS);
58			X=0.5*(LONS(1:end-1)+LONS(2:end))';
59			Y=0.5*(LATS(1:end-1)+LATS(2:end));
60	gforget	1.1	elseif ~isempty(whos('LONS'));
61	gforget	1.3	valranges=LONS';
62			X=0.5*(LONS(1:end-1)+LONS(2:end))';
63			Y=NaN;
64	gforget	1.1	else;
65	gforget	1.3	valranges=i*LATS;
66			X=NaN;
67			Y=0.5*(LATS(1:end-1)+LATS(2:end));
68	gforget	1.1	end;
69			nnranges=size(valranges);
70			nnout=max(size(valranges)-1,[1 1]);
71			%
72			FLD=NaN*ones([nnout n3 n4]);
73			W=NaN*ones([nnout n3 n4]);
74	gforget	1.3	X=repmat(X*ones(1,size(Y,2)),[1 1 n3 n4]);
75			Y=repmat(ones(size(X,1),1)*Y,[1 1 n3 n4]);
76
77			%and corresponding vertical positions:
78			Z=NaN; if n3>1; Z=[]; Z(1,1,:)=[1:n3]; end;
79			if isempty(whos('weights'))&isempty(whos('level'));
80			Z(1,1,:)=mygrid.RC;
81			elseif ~isempty(whos('level'));
82			Z(1,1,:)=mygrid.RC(level);
83			end;
84			[t1,t2,t3,t4]=size(X);
85			Z=repmat(Z,[t1 t2 1 t4]);
86	gforget	1.1
87			%select weights for average:
88			if isempty(whos('weights'))&isempty(whos('level'));
89			weights=mygrid.hFacC.*mk3D(mygrid.RAC,mygrid.hFacC);
90			weights=weights.*mk3D(mygrid.DRF,mygrid.hFacC);
91			elseif ~isempty(whos('level'));
92			weights=mygrid.hFacC(:,:,level).mygrid.RACmygrid.DRF(level);
93			end;
94	gforget	1.2
95	gforget	1.3	%check for potential inconsistencies in specified weights :
96			[w1,w2,w3,w4]=size(weights.f1);
97			if w3==1&n3>1; weights=repmat(weights,[1 1 n3 n4]); end;
98			if w4==1&n4>1; weights=repmat(weights,[1 1 1 n4]); end;
99			test1=sum(weights>0&isnan(fld));
100			if test1>0; error('non-zero weights found for NaN point'); end;
101
102	gforget	1.2	%switch to 2D array to speed up computation:
103			fld=convert2array(fld);
104			n1=size(fld,1); n2=size(fld,2);
105			fld=reshape(fld,n1n2,n3n4);
106			%same for the weights:
107	gforget	1.3	weights=convert2array(weights);
108	gforget	1.1	weights=reshape(weights,n1n2,n3n4);
109	gforget	1.2	%multiply one with the other
110	gforget	1.1	fld=fld.*weights;
111	gforget	1.2	%remove data mask
112			fld(isnan(fld))=0;
113	gforget	1.1
114			lonvec=reshape(convert2array(mygrid.XC),n1*n2,1);
115			latvec=reshape(convert2array(mygrid.YC),n1*n2,1);
116
117			for ix=1:nnout(1);
118			for iy=1:nnout(2);
119
120			%get list ofpoints that form a zonal band:
121	gforget	1.3	if ~isempty(whos('LONS'))&~isempty(whos('LATS'));
122			mm=find(latvec>=imag(valranges(ix,iy))&latvec<imag(valranges(ix,iy+1))...
123			&lonvec>=real(valranges(ix,iy))&lonvec<real(valranges(ix+1,iy)));
124			elseif ~isempty(whos('LATS'));
125	gforget	1.1	mm=find(latvec>=imag(valranges(ix,iy))&latvec<imag(valranges(ix,iy+1)));
126	gforget	1.3	else;
127	gforget	1.1	mm=find(lonvec>=real(valranges(ix,iy))&lonvec<real(valranges(ix+1,iy)));
128			end;
129
130			%do the area weighed average along this band:
131	gforget	1.3	tmp1=nansum(fld(mm,:),1);
132			tmp2=nansum(weights(mm,:),1);
133	gforget	1.1	tmp2(tmp2==0)=NaN;
134			tmp1=tmp1./tmp2;
135
136			%store:
137			FLD(ix,iy,:,:)=reshape(tmp1,n3,n4);
138			W(ix,iy,:,:)=reshape(tmp2,n3,n4);
139
140			end;
141			end;
142
143	gforget	1.3	%remove singleton dimensions:
144			X=squeeze(X); Y=squeeze(Y); Z=squeeze(Z);
145			W=squeeze(W); FLD=squeeze(FLD);
146
147			%prepare outout:
148			if nargout==5;
149			varargout={FLD,X,Y,Z,W};
150			elseif nargout==4;
151			varargout={FLD,X,Y,Z};
152			elseif nargout==2;
153	gforget	1.1	varargout={FLD,W};
154			elseif nargout==1;
155			varargout={FLD};
156			else;
157			varargout={};
158			end;
159
160
161