matlab_class/gcmfaces_calc/calc_overturn.m

function [FLD]=calc_overturn(fldU,fldV,doFlip,list_factors);
%object:    compute meridional overturning streamfunction
%inputs:    fldU and fldV are the fields of grid point transport
%optional:  doFlip (default is 1). If 1 then flip the vertical 
%               axis back and forth, hence intergrating from the 
%               'bottom'. If 0 then dont.
%           list_factors (default is {'dh','dz'})
%output:    FLD is the streamfunction
%
%notes:     mygrid.LATS_MASKS is the set of quasi longitudinal lines along which
%               transports will integrated, as computed in gcmfaces_lines_zonal
%           the result is converted to Sv, and sign is changed.

gcmfaces_global;

if nargin<3; doFlip=1; end;
if nargin<4; list_factors={'dh','dz'}; 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;

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

%prepare fldU/fldV:
fldU(isnan(fldU))=0; fldV(isnan(fldV))=0;

dxg=mk3D(mygrid.DXG,fldU); dyg=mk3D(mygrid.DYG,fldU);
if size(fldU.f1,3)==length(mygrid.DRF); drf=mk3D(mygrid.DRF,fldU); else; drf=fldU; drf(:)=1; end;
facW=drf; facW(:)=1; facS=facW;
for ii=1:length(list_factors);
  tmp1=list_factors{ii};
  if strcmp(tmp1,'dh'); facW=facW.*dyg; facS=facS.*dxg;
  elseif strcmp(tmp1,'dz'); facW=facW.*drf; facS=facS.*drf;
  elseif strcmp(tmp1,'hfac'); facW=facW.*mygrid.hFacW; facS=facS.*mygrid.hFacS;
  elseif isempty(tmp1); 1;
  else; fprintf('error in calc_overturn : non supported factor\n'); return;
  end;
end;

for k4=1:n4;
fldU(:,:,:,k4)=fldU(:,:,:,k4).*facW;
fldV(:,:,:,k4)=fldV(:,:,:,k4).*facS;
end;

%use array format to speed up computation below:
fldU=convert2array(fldU); fldV=convert2array(fldV);

for iy=1:length(mygrid.LATS_MASKS); 

   %get list ofpoints that form a zonal band:
   mskW=mygrid.LATS_MASKS(iy).mskWedge;
   vecW=gcmfaces_subset(mskW,fldU);
   mskS=mygrid.LATS_MASKS(iy).mskSedge;
   vecS=gcmfaces_subset(mskS,fldV);
   trsp=nansum(vecW,1)+nansum(vecS,1);

   %store:
   if doFlip;
     FLD(iy,1:n3,:)=flipdim(cumsum(flipdim(trsp,2),2),2);
   else;
     FLD(iy,2:n3+1,:)=cumsum(trsp,2);
   end;

   %convert to Sv and change sign:
   FLD(iy,:,:)=-1e-6*FLD(iy,:,:);

end; 

if doFlip;
  FLD(:,end,:)=0;
else;
  FLD(:,1,:)=0;
end;

1	gforget	1.6	function [FLD]=calc_overturn(fldU,fldV,doFlip,list_factors);
2	gforget	1.2	%object: compute meridional overturning streamfunction
3			%inputs: fldU and fldV are the fields of grid point transport
4	gforget	1.6	%optional: doFlip (default is 1). If 1 then flip the vertical
5			% axis back and forth, hence intergrating from the
6			% 'bottom'. If 0 then dont.
7			% list_factors (default is {'dh','dz'})
8	gforget	1.2	%output: FLD is the streamfunction
9	gforget	1.3	%
10			%notes: mygrid.LATS_MASKS is the set of quasi longitudinal lines along which
11			% transports will integrated, as computed in gcmfaces_lines_zonal
12	gforget	1.6	% the result is converted to Sv, and sign is changed.
13	gforget	1.1
14	gforget	1.10	gcmfaces_global;
15	gforget	1.1
16	gforget	1.6	if nargin<3; doFlip=1; end;
17			if nargin<4; list_factors={'dh','dz'}; end;
18
19	gforget	1.10	%check that LATS_MASKS has already been defined:
20			if ~isfield(mygrid,'LATS_MASKS');
21			fprintf('one-time initialization of gcmfaces_lines_zonal: begin\n');
22			gcmfaces_lines_zonal;
23			fprintf('one-time initialization of gcmfaces_lines_zonal: end\n');
24			end;
25
26	gforget	1.1	%initialize output:
27			n3=max(size(fldU.f1,3),1); n4=max(size(fldV.f1,4),1);
28	gforget	1.3	FLD=NaN*squeeze(zeros(length(mygrid.LATS_MASKS),n3+1,n4));
29	gforget	1.1
30			%prepare fldU/fldV:
31			fldU(isnan(fldU))=0; fldV(isnan(fldV))=0;
32
33	gforget	1.6	dxg=mk3D(mygrid.DXG,fldU); dyg=mk3D(mygrid.DYG,fldU);
34			if size(fldU.f1,3)==length(mygrid.DRF); drf=mk3D(mygrid.DRF,fldU); else; drf=fldU; drf(:)=1; end;
35			facW=drf; facW(:)=1; facS=facW;
36			for ii=1:length(list_factors);
37			tmp1=list_factors{ii};
38			if strcmp(tmp1,'dh'); facW=facW.dyg; facS=facS.dxg;
39			elseif strcmp(tmp1,'dz'); facW=facW.drf; facS=facS.drf;
40			elseif strcmp(tmp1,'hfac'); facW=facW.mygrid.hFacW; facS=facS.mygrid.hFacS;
41			elseif isempty(tmp1); 1;
42	gforget	1.9	else; fprintf('error in calc_overturn : non supported factor\n'); return;
43	gforget	1.6	end;
44			end;
45
46	gforget	1.1	for k4=1:n4;
47	gforget	1.6	fldU(:,:,:,k4)=fldU(:,:,:,k4).*facW;
48			fldV(:,:,:,k4)=fldV(:,:,:,k4).*facS;
49	gforget	1.1	end;
50
51			%use array format to speed up computation below:
52			fldU=convert2array(fldU); fldV=convert2array(fldV);
53
54	gforget	1.3	for iy=1:length(mygrid.LATS_MASKS);
55	gforget	1.1
56			%get list ofpoints that form a zonal band:
57	gforget	1.5	mskW=mygrid.LATS_MASKS(iy).mskWedge;
58			vecW=gcmfaces_subset(mskW,fldU);
59			mskS=mygrid.LATS_MASKS(iy).mskSedge;
60			vecS=gcmfaces_subset(mskS,fldV);
61	gforget	1.7	trsp=nansum(vecW,1)+nansum(vecS,1);
62	gforget	1.1
63			%store:
64	gforget	1.8	if doFlip;
65			FLD(iy,1:n3,:)=flipdim(cumsum(flipdim(trsp,2),2),2);
66			else;
67			FLD(iy,2:n3+1,:)=cumsum(trsp,2);
68			end;
69	gforget	1.6
70			%convert to Sv and change sign:
71	gforget	1.8	FLD(iy,:,:)=-1e-6*FLD(iy,:,:);
72	gforget	1.1
73			end;
74
75	gforget	1.8	if doFlip;
76			FLD(:,end,:)=0;
77			else;
78			FLD(:,1,:)=0;
79			end;
80	gforget	1.1