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	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	% If specifying 'weights' you likely want to set it to 0 on land.
12	% (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) 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	%usage:
24	% input 2a and/or 2b is necessary
25	% 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
32	gcmfaces_global;
33
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	[n1,n2,n3,n4]=size(fld.f1);
41
42	%test inputs/outputs consistency:
43	tmp1=isempty(whos('LONS'))&isempty(whos('LATS'));
44	if tmp1;
45	error('wrong input specification : specify LONS and/or LATS');
46	end;
47	tmp1=~isempty(whos('weights'))+~isempty(whos('level'));
48	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
55	%prepare output etc.:
56	if ~isempty(whos('LONS'))&~isempty(whos('LATS')); %use complex
57	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	elseif ~isempty(whos('LONS'));
61	valranges=LONS';
62	X=0.5*(LONS(1:end-1)+LONS(2:end))';
63	Y=NaN;
64	else;
65	valranges=i*LATS;
66	X=NaN;
67	Y=0.5*(LATS(1:end-1)+LATS(2:end));
68	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	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
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
95	%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	%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	weights=convert2array(weights);
108	weights=reshape(weights,n1n2,n3n4);
109	%multiply one with the other
110	fld=fld.*weights;
111	%remove data mask
112	fld(isnan(fld))=0;
113
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	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	mm=find(latvec>=imag(valranges(ix,iy))&latvec<imag(valranges(ix,iy+1)));
126	else;
127	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	tmp1=nansum(fld(mm,:),1);
132	tmp2=nansum(weights(mm,:),1);
133	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	%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	varargout={FLD,W};
154	elseif nargout==1;
155	varargout={FLD};
156	else;
157	varargout={};
158	end;
159
160
161