matlab_class/gcmfaces_calc/calc_zonmean_T.m

function [fldOut,X,Y,weightOut]=calc_zonmean_T(fldIn,method,fldType);
% CALC_ZONMEAN_T(budgIn,method,fldType)
%    computes zonal average of fldIn (or its fields recursively).
%
%    If method is 1 (default) then mskCedge (from mygrid.LATS_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).

gcmfaces_global;

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

%check that LATS_MASKS has already been defined:
if ~isfield(mygrid,'LATS_MASKS');
    fprintf('one-time initialization of gcmfaces_lines_zonal: begin\n');
    gcmfaces_lines_zonal;
    fprintf('one-time initialization of gcmfaces_lines_zonal: end\n');
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_zonmean_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.LATS_MASKS),n3,n4));
weightOut=NaN*squeeze(zeros(length(mygrid.LATS_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.LATS_MASKS);
for iy=1:ny;

  if method>0;    
    %get list of points that form a zonal band:
    mm=convert2gcmfaces(mygrid.LATS_MASKS(iy).mskCedge);
    mm=find(~isnan(mm)&mm~=0);
  else;
    if iy>1&iy<ny;
      tmpMin=0.5*(mygrid.LATS(iy-1)+mygrid.LATS(iy));
      tmpMax=0.5*(mygrid.LATS(iy)+mygrid.LATS(iy+1));
    elseif iy==1;
      tmpMin=-Inf;
      tmpMax=0.5*(mygrid.LATS(iy)+mygrid.LATS(iy+1));
    elseif iy==ny;
      tmpMin=0.5*(mygrid.LATS(iy-1)+mygrid.LATS(iy));
      tmpMax=+Inf;
    end;
    mm=convert2gcmfaces(mygrid.YC>=tmpMin&mygrid.YC<tmpMax);
    mm=find(~isnan(mm)&mm~=0);
  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:
  if ~isempty(mm);
    fldOut(iy,:,:)=reshape(tmp1,n3,n4);
    weightOut(iy,:,:)=reshape(tmp2,n3,n4);
  end;
   
end;

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

1	function [fldOut,X,Y,weightOut]=calc_zonmean_T(fldIn,method,fldType);
2	% CALC_ZONMEAN_T(budgIn,method,fldType)
3	% computes zonal average of fldIn (or its fields recursively).
4	%
5	% If method is 1 (default) then mskCedge (from mygrid.LATS_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	gcmfaces_global;
14
15	if isempty(who('fldType')); fldType='intensive'; end;
16	if isempty(who('method')); method=1; end;
17
18	%check that LATS_MASKS has already been defined:
19	if ~isfield(mygrid,'LATS_MASKS');
20	fprintf('one-time initialization of gcmfaces_lines_zonal: begin\n');
21	gcmfaces_lines_zonal;
22	fprintf('one-time initialization of gcmfaces_lines_zonal: end\n');
23	end;
24
25	if isa(fldIn,'struct');
26	list0=fieldnames(fldIn);
27	fldOut=[];
28	for vv=1:length(list0);
29	tmp1=getfield(fldIn,list0{vv});
30	if isa(tmp1,'gcmfaces');
31	[tmp2,X,Y,weightOut]=calc_zonmean_T(tmp1,method,fldType);
32	fldOut=setfield(fldOut,list0{vv},tmp2);
33	end;
34	end;
35	return;
36	end;
37
38	%initialize output:
39	n3=max(size(fldIn.f1,3),1); n4=max(size(fldIn.f1,4),1);
40	fldOut=NaN*squeeze(zeros(length(mygrid.LATS_MASKS),n3,n4));
41	weightOut=NaN*squeeze(zeros(length(mygrid.LATS_MASKS),n3,n4));
42
43	%use array format to speed up computation below:
44	fldIn=convert2gcmfaces(fldIn);
45	n1=size(fldIn,1); n2=size(fldIn,2);
46	fldIn=reshape(fldIn,n1n2,n3n4);
47
48	%set rac and hFacC according to method
49	if abs(method)==1;
50	rac=reshape(convert2gcmfaces(mygrid.RAC),n1n2,1)ones(1,n3*n4);
51	if n3==length(mygrid.RC);
52	hFacC=reshape(convert2gcmfaces(mygrid.hFacC),n1*n2,n3);
53	hFacC=repmat(hFacC,[1 n4]);
54	DRF=repmat(mygrid.DRF',[n1*n2 n4]);
55	else;
56	hFacC=reshape(convert2gcmfaces(mygrid.mskC(:,:,1)),n1n2,1)ones(1,n3*n4);
57	hFacC(isnan(hFacC))=0;
58	DRF=repmat(mygrid.DRF(1),[n1n2 n3n4]);
59	end;
60	weight=rac.hFacC.DRF;
61	else;
62	weight=mygrid.mskC(:,:,1).*mygrid.RAC;
63	weight=reshape(convert2gcmfaces(weight),n1n2,1)ones(1,n3*n4);
64	end;
65
66	%masked area only:
67	weight(isnan(fldIn))=0;
68	weight(isnan(weight))=0;
69	mask=weight; mask(weight~=0)=1;
70	fldIn(isnan(fldIn))=0;
71
72	ny=length(mygrid.LATS_MASKS);
73	for iy=1:ny;
74
75	if method>0;
76	%get list of points that form a zonal band:
77	mm=convert2gcmfaces(mygrid.LATS_MASKS(iy).mskCedge);
78	mm=find(~isnan(mm)&mm~=0);
79	else;
80	if iy>1&iy<ny;
81	tmpMin=0.5*(mygrid.LATS(iy-1)+mygrid.LATS(iy));
82	tmpMax=0.5*(mygrid.LATS(iy)+mygrid.LATS(iy+1));
83	elseif iy==1;
84	tmpMin=-Inf;
85	tmpMax=0.5*(mygrid.LATS(iy)+mygrid.LATS(iy+1));
86	elseif iy==ny;
87	tmpMin=0.5*(mygrid.LATS(iy-1)+mygrid.LATS(iy));
88	tmpMax=+Inf;
89	end;
90	mm=convert2gcmfaces(mygrid.YC>=tmpMin&mygrid.YC<tmpMax);
91	mm=find(~isnan(mm)&mm~=0);
92	end;
93
94	if strcmp(fldType,'intensive');
95	tmp1=nansum(fldIn(mm,:).*weight(mm,:),1)./nansum(weight(mm,:),1);
96	tmp2=nansum(weight(mm,:),1);
97	else;
98	tmp1=nansum(fldIn(mm,:).*mask(mm,:),1)./nansum(weight(mm,:),1);
99	tmp2=nansum(weight(mm,:),1);
100	end;
101
102	%store:
103	if ~isempty(mm);
104	fldOut(iy,:,:)=reshape(tmp1,n3,n4);
105	weightOut(iy,:,:)=reshape(tmp2,n3,n4);
106	end;
107
108	end;
109
110	X=[]; Y=[];
111	if size(fldOut,2)==length(mygrid.RC);
112	X=mygrid.LATS*ones(1,length(mygrid.RC));
113	Y=ones(length(mygrid.LATS),1)*(mygrid.RC');
114	elseif size(fldOut,2)==1;
115	X=mygrid.LATS;
116	Y=ones(length(mygrid.LATS),1);
117	end;
118