matlab_class/gcmfaces_calc/calc_mermean_T.m

function [fldOut,X,Y,weightOut]=calc_mermean_T(fldIn,method,fldType);
% CALC_MERMEAN_T(budgIn,method,fldType)
%    computes meridional average of fldIn (or its fields recursively).
%
%    If method is 1 (default) then mskCedge (from mygrid.LONS_MASKS) 
%    and volume elements used; if method is 2 then mskCedge and surface 
%    elements are used; if method is -1 or -2 then mskC is used (from mygrid) 
%    instead of mskCedge to define the averaging footpring.
%    
%    If fldType is 'intensive' (default) then fldIn is mutliplied by 
%    RAC (method=2 or -2) or RAC.*hFacC*DRF (method=1 or -1).

global mygrid;

if ~isfield(mygrid,'LONS_MASKS');
    error('please execute gcmfaces_lines_zonal to define mygrid.LONS_MASKS');
end;

if isempty(who('fldType')); fldType='intensive'; end;
if isempty(who('method')); method=1; end;

if isa(fldIn,'struct');
  list0=fieldnames(fldIn);
  fldOut=[];
  for vv=1:length(list0);
    tmp1=getfield(fldIn,list0{vv});
    if isa(tmp1,'gcmfaces');
      [tmp2,X,Y,weightOut]=calc_mermean_T(tmp1,method,fldType);
      fldOut=setfield(fldOut,list0{vv},tmp2);
    end;
  end;
  return;
end;

%initialize output:
n3=max(size(fldIn.f1,3),1); n4=max(size(fldIn.f1,4),1);
fldOut=NaN*squeeze(zeros(length(mygrid.LONS_MASKS),n3,n4));
weightOut=NaN*squeeze(zeros(length(mygrid.LONS_MASKS),n3,n4));

%use array format to speed up computation below:
fldIn=convert2gcmfaces(fldIn);
n1=size(fldIn,1); n2=size(fldIn,2);
fldIn=reshape(fldIn,n1*n2,n3*n4);

%set rac and hFacC according to method
if abs(method)==1;
  rac=reshape(convert2gcmfaces(mygrid.RAC),n1*n2,1)*ones(1,n3*n4);
  if n3==length(mygrid.RC);
      hFacC=reshape(convert2gcmfaces(mygrid.hFacC),n1*n2,n3);
      hFacC=repmat(hFacC,[1 n4]);
      DRF=repmat(mygrid.DRF',[n1*n2 n4]);
  else;
      hFacC=reshape(convert2gcmfaces(mygrid.mskC(:,:,1)),n1*n2,1)*ones(1,n3*n4);
      hFacC(isnan(hFacC))=0;
      DRF=repmat(mygrid.DRF(1),[n1*n2 n3*n4]);
  end;
  weight=rac.*hFacC.*DRF;
else;
  weight=mygrid.mskC(:,:,1).*mygrid.RAC;
  weight=reshape(convert2gcmfaces(weight),n1*n2,1)*ones(1,n3*n4);
end;

%masked area only:
weight(isnan(fldIn))=0;
weight(isnan(weight))=0;
mask=weight; mask(weight~=0)=1;
fldIn(isnan(fldIn))=0;

ny=length(mygrid.LONS_MASKS);
for iy=1:ny;

  if method>0;    
    %get list of points that form a zonal band:
    mm=convert2gcmfaces(mygrid.LONS_MASKS(iy).mskCedge);
    mm=find(~isnan(mm)&mm~=0);
  else;
      error('method < 0 remains to be implemented');
  end;

  if strcmp(fldType,'intensive');
    tmp1=nansum(fldIn(mm,:).*weight(mm,:),1)./nansum(weight(mm,:),1);
    tmp2=nansum(weight(mm,:),1);
  else;
    tmp1=nansum(fldIn(mm,:).*mask(mm,:),1)./nansum(weight(mm,:),1);
    tmp2=nansum(weight(mm,:),1);
  end;

  %store:
  fldOut(iy,:,:)=reshape(tmp1,n3,n4);
  weightOut(iy,:,:)=reshape(tmp2,n3,n4);
    
end;

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

1	gforget	1.1	function [fldOut,X,Y,weightOut]=calc_mermean_T(fldIn,method,fldType);
2			% CALC_MERMEAN_T(budgIn,method,fldType)
3			% computes meridional average of fldIn (or its fields recursively).
4			%
5			% If method is 1 (default) then mskCedge (from mygrid.LONS_MASKS)
6			% and volume elements used; if method is 2 then mskCedge and surface
7			% elements are used; if method is -1 or -2 then mskC is used (from mygrid)
8			% instead of mskCedge to define the averaging footpring.
9			%
10			% If fldType is 'intensive' (default) then fldIn is mutliplied by
11			% RAC (method=2 or -2) or RAC.hFacCDRF (method=1 or -1).
12
13			global mygrid;
14
15			if ~isfield(mygrid,'LONS_MASKS');
16			error('please execute gcmfaces_lines_zonal to define mygrid.LONS_MASKS');
17			end;
18
19			if isempty(who('fldType')); fldType='intensive'; end;
20			if isempty(who('method')); method=1; end;
21
22			if isa(fldIn,'struct');
23			list0=fieldnames(fldIn);
24			fldOut=[];
25			for vv=1:length(list0);
26			tmp1=getfield(fldIn,list0{vv});
27			if isa(tmp1,'gcmfaces');
28			[tmp2,X,Y,weightOut]=calc_mermean_T(tmp1,method,fldType);
29			fldOut=setfield(fldOut,list0{vv},tmp2);
30			end;
31			end;
32			return;
33			end;
34
35			%initialize output:
36			n3=max(size(fldIn.f1,3),1); n4=max(size(fldIn.f1,4),1);
37			fldOut=NaN*squeeze(zeros(length(mygrid.LONS_MASKS),n3,n4));
38			weightOut=NaN*squeeze(zeros(length(mygrid.LONS_MASKS),n3,n4));
39
40			%use array format to speed up computation below:
41			fldIn=convert2gcmfaces(fldIn);
42			n1=size(fldIn,1); n2=size(fldIn,2);
43			fldIn=reshape(fldIn,n1n2,n3n4);
44
45			%set rac and hFacC according to method
46			if abs(method)==1;
47			rac=reshape(convert2gcmfaces(mygrid.RAC),n1n2,1)ones(1,n3*n4);
48			if n3==length(mygrid.RC);
49			hFacC=reshape(convert2gcmfaces(mygrid.hFacC),n1*n2,n3);
50			hFacC=repmat(hFacC,[1 n4]);
51			DRF=repmat(mygrid.DRF',[n1*n2 n4]);
52			else;
53			hFacC=reshape(convert2gcmfaces(mygrid.mskC(:,:,1)),n1n2,1)ones(1,n3*n4);
54			hFacC(isnan(hFacC))=0;
55			DRF=repmat(mygrid.DRF(1),[n1n2 n3n4]);
56			end;
57			weight=rac.hFacC.DRF;
58			else;
59			weight=mygrid.mskC(:,:,1).*mygrid.RAC;
60			weight=reshape(convert2gcmfaces(weight),n1n2,1)ones(1,n3*n4);
61			end;
62
63			%masked area only:
64			weight(isnan(fldIn))=0;
65			weight(isnan(weight))=0;
66			mask=weight; mask(weight~=0)=1;
67			fldIn(isnan(fldIn))=0;
68
69			ny=length(mygrid.LONS_MASKS);
70			for iy=1:ny;
71
72			if method>0;
73			%get list of points that form a zonal band:
74			mm=convert2gcmfaces(mygrid.LONS_MASKS(iy).mskCedge);
75			mm=find(~isnan(mm)&mm~=0);
76			else;
77			error('method < 0 remains to be implemented');
78			end;
79
80			if strcmp(fldType,'intensive');
81			tmp1=nansum(fldIn(mm,:).*weight(mm,:),1)./nansum(weight(mm,:),1);
82			tmp2=nansum(weight(mm,:),1);
83			else;
84			tmp1=nansum(fldIn(mm,:).*mask(mm,:),1)./nansum(weight(mm,:),1);
85			tmp2=nansum(weight(mm,:),1);
86			end;
87
88			%store:
89			fldOut(iy,:,:)=reshape(tmp1,n3,n4);
90			weightOut(iy,:,:)=reshape(tmp2,n3,n4);
91
92			end;
93
94			X=[]; Y=[];
95			if size(fldOut,2)==length(mygrid.RC);
96			X=mygrid.LONS*ones(1,length(mygrid.RC));
97			Y=ones(length(mygrid.LONS),1)*(mygrid.RC');
98			elseif size(fldOut,2)==1;
99			X=mygrid.LONS;
100			Y=ones(length(mygrid.LONS),1);
101			end;
102