matlab_class/gcmfaces_calc/calc_zonmean_T.m

function [FLD,X,Y]=calc_zonmean_T(fld);
%object:    compute zonal mean
%inputs:    fld is the field of interest
%output:    FLD is the zonal mean field
%
%notes:     mygrid.LATS_MASKS is the set of quasi longitudinal lines along which
%               means will computed, as computed in gcmfaces_lines_zonal

global mygrid;

%initialize output:
n3=max(size(fld.f1,3),1); n4=max(size(fld.f1,4),1);
FLD=NaN*squeeze(zeros(length(mygrid.LATS_MASKS),n3,n4));

%use array format to speed up computation below:
fld=convert2array(fld);
n1=size(fld,1); n2=size(fld,2);
fld=reshape(fld,n1*n2,n3*n4);
rac=reshape(convert2array(mygrid.RAC),n1*n2,1)*ones(1,n3*n4);
if n3==length(mygrid.RC);
    hFacC=reshape(convert2array(mygrid.hFacC),n1*n2,n3);
    hFacC=repmat(hFacC,[1 n4]);
else;
    hFacC=reshape(convert2array(mygrid.mskC(:,:,1)),n1*n2,1)*ones(1,n3*n4);
    hFacC(isnan(hFacC))=0;
end;
%masked area only:
rac(isnan(fld))=0;
hFacC(isnan(fld))=0;
fld(isnan(fld))=0;

rac=rac.*hFacC;
fld=fld.*rac;

for iy=1:length(mygrid.LATS_MASKS);
    
    %get list ofpoints that form a zonal band:
    mm=convert2array(mygrid.LATS_MASKS(iy).mskCedge);
    mm=find(~isnan(mm)&mm~=0);
    
    %do the area weighed average along this band:
    tmp1=sum(fld(mm,:),1);
    tmp2=sum(rac(mm,:),1);
    tmp2(tmp2==0)=NaN;
    tmp1=tmp1./tmp2;
    
    %store:
    FLD(iy,:,:)=reshape(tmp1,n3,n4);
    
end;

X=[]; Y=[];
if size(FLD,2)==length(mygrid.RC);
    X=mygrid.LATS*ones(1,length(mygrid.RC));
    Y=ones(length(mygrid.LATS),1)*(mygrid.RC');
elseif size(FLD,2)==1;
    X=mygrid.LATS;
    Y=ones(length(mygrid.LATS),1);
end;

1	gforget	1.7	function [FLD,X,Y]=calc_zonmean_T(fld);
2	gforget	1.3	%object: compute zonal mean
3			%inputs: fld is the field of interest
4			%output: FLD is the zonal mean field
5	gforget	1.4	%
6			%notes: mygrid.LATS_MASKS is the set of quasi longitudinal lines along which
7			% means will computed, as computed in gcmfaces_lines_zonal
8	gforget	1.1
9			global mygrid;
10
11			%initialize output:
12			n3=max(size(fld.f1,3),1); n4=max(size(fld.f1,4),1);
13	gforget	1.4	FLD=NaN*squeeze(zeros(length(mygrid.LATS_MASKS),n3,n4));
14	gforget	1.1
15			%use array format to speed up computation below:
16	gforget	1.7	fld=convert2array(fld);
17			n1=size(fld,1); n2=size(fld,2);
18			fld=reshape(fld,n1n2,n3n4);
19	gforget	1.1	rac=reshape(convert2array(mygrid.RAC),n1n2,1)ones(1,n3*n4);
20	gforget	1.7	if n3==length(mygrid.RC);
21	gforget	1.8	hFacC=reshape(convert2array(mygrid.hFacC),n1*n2,n3);
22			hFacC=repmat(hFacC,[1 n4]);
23	gforget	1.7	else;
24			hFacC=reshape(convert2array(mygrid.mskC(:,:,1)),n1n2,1)ones(1,n3*n4);
25			hFacC(isnan(hFacC))=0;
26	gforget	1.1	end;
27	gforget	1.6	%masked area only:
28			rac(isnan(fld))=0;
29			hFacC(isnan(fld))=0;
30			fld(isnan(fld))=0;
31
32	gforget	1.1	rac=rac.*hFacC;
33			fld=fld.*rac;
34
35	gforget	1.7	for iy=1:length(mygrid.LATS_MASKS);
36
37			%get list ofpoints that form a zonal band:
38			mm=convert2array(mygrid.LATS_MASKS(iy).mskCedge);
39			mm=find(~isnan(mm)&mm~=0);
40
41			%do the area weighed average along this band:
42			tmp1=sum(fld(mm,:),1);
43			tmp2=sum(rac(mm,:),1);
44			tmp2(tmp2==0)=NaN;
45			tmp1=tmp1./tmp2;
46
47			%store:
48			FLD(iy,:,:)=reshape(tmp1,n3,n4);
49
50			end;
51	gforget	1.1
52	gforget	1.7	X=[]; Y=[];
53			if size(FLD,2)==length(mygrid.RC);
54			X=mygrid.LATS*ones(1,length(mygrid.RC));
55			Y=ones(length(mygrid.LATS),1)*(mygrid.RC');
56			elseif size(FLD,2)==1;
57			X=mygrid.LATS;
58			Y=ones(length(mygrid.LATS),1);
59			end;
60	gforget	1.1