matlab_class/gcmfaces_calc/calc_barostream.m

function [fldBAR]=calc_barostream(fldU,fldV,noDiv,list_factors);
%object:    compute barotropic streamfunction
%inputs:    fldU and fldV are the fields of grid point transport
%optional:  noDiv (default is 1). If 1 then remove the divergent 
%               part of the flow field first. If 0 then dont.
%           list_factors (default is {'dh','dz'})
%output:    FLD is the streamfunction
%notes:     the result is converted to Sv

global mygrid;

if nargin<3; noDiv=1; end;
if nargin<4; list_factors={'dh','dz'}; end;

%0) prepare fldU/fldV (transport fields):
n3=max(size(fldU.f1,3),1); n4=max(size(fldV.f1,4),1);

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_barostream : non supported factor\n'); return;
  end;
end;

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

%apply mask:
fldU=sum(fldU,3).*mygrid.mskW(:,:,1); 
fldV=sum(fldV,3).*mygrid.mskS(:,:,1);
%take out thedivergent part of the flow:
if noDiv;
  [fldUdiv,fldVdiv,fldDivPot]=diffsmooth2D_div_inv(fldU,fldV);
  fldU=fldU-fldUdiv; fldV=fldV-fldVdiv;
end;

%1) compute streamfunction face by face:
[fldU,fldV]=exch_UV(fldU,fldV); fldU(isnan(fldU))=0; fldV(isnan(fldV))=0;
tmp1=cumsum(fldV,1); for iF=1:fldU.nFaces; tmp2=tmp1{iF}; tmp1{iF}=[zeros(1,size(tmp2,2));tmp2]; end;
tmp2=diff(tmp1,1,2)+fldU; tmp3=cumsum(mean(tmp2,1)); 
%to check divergence implied errors:  
%figure; for iF=1:fldU.nFaces; subplot(3,2,iF); plot(std(tmp2{iF},0,1)); end;
for iF=1:fldU.nFaces; tmp2=tmp1{iF}; tmp1{iF}=tmp2-ones(size(tmp2,1),1)*[0 tmp3{iF}]; end;
bf_step1=tmp1;

if fldU.nFaces==1;
  bf_step2=bf_step1;
else;
%2) match edges:
%... set face number field
TMP1=tmp1; for iF=1:TMP1.nFaces; TMP1{iF}(:)=iF; end; 
TMP2=exch_T_N(TMP1);%!!! this is a trick, since TMP1 is (n+1 X n+1) and loc. at vorticity points
tmp1=bf_step1;
for iF=1:fldU.nFaces-1;
  tmp2=exch_T_N(tmp1);%!!! same trick 
  tmp3=tmp2{iF+1}; tmp3(3:end-2,3:end-2)=NaN;%mask out interior points
  TMP3=TMP2{iF+1}; tmp3(find(TMP3>iF+1))=NaN;%mask out edges points coming from unadjusted faces
  tmp3(:,1)=tmp3(:,1)-tmp3(:,2); tmp3(:,end)=tmp3(:,end)-tmp3(:,end-1);%compare edge points
  tmp3(1,:)=tmp3(1,:)-tmp3(2,:); tmp3(end,:)=tmp3(end,:)-tmp3(end-1,:);%compare edge points
  tmp3(2:end-1,2:end-1)=NaN;%mask out remaining interior points
  tmp1{iF+1}=tmp1{iF+1}+nanmedian(tmp3(:));%adjust the face data
end;
bf_step2=tmp1;
end;
%3) put streamfunction at cell center
tmp1=bf_step2;
tmp2=tmp1; for iF=1:tmp1.nFaces; tmp3=tmp2{iF}; tmp3=(tmp3(:,1:end-1)+tmp3(:,2:end))/2;
tmp3=(tmp3(1:end-1,:)+tmp3(2:end,:))/2; tmp2{iF}=tmp3; end;
bf_step3=tmp2;
%4) set 0 on land on average:
tmp1=convert2array(bf_step3); 
tmp2=convert2array(mygrid.mskC(:,:,1)); tmp2=find(isnan(tmp2)&~isnan(tmp1)); 
tmp2=median(tmp1(tmp2)); if isnan(tmp2); tmp2=nanmedian(tmp1(:)); end; 
bf_step4=(bf_step3-tmp2).*mygrid.mskC(:,:,1);

%5) return the result:
fldBAR=bf_step4;

%convert to Sv and change sign:
fldBAR=1e-6*fldBAR;

1	gforget	1.7	function [fldBAR]=calc_barostream(fldU,fldV,noDiv,list_factors);
2	gforget	1.5	%object: compute barotropic streamfunction
3			%inputs: fldU and fldV are the fields of grid point transport
4	gforget	1.7	%optional: noDiv (default is 1). If 1 then remove the divergent
5			% part of the flow field first. If 0 then dont.
6			% list_factors (default is {'dh','dz'})
7	gforget	1.5	%output: FLD is the streamfunction
8	gforget	1.7	%notes: the result is converted to Sv
9	gforget	1.1
10			global mygrid;
11
12	gforget	1.7	if nargin<3; noDiv=1; end;
13			if nargin<4; list_factors={'dh','dz'}; end;
14
15	gforget	1.2	%0) prepare fldU/fldV (transport fields):
16	gforget	1.1	n3=max(size(fldU.f1,3),1); n4=max(size(fldV.f1,4),1);
17
18			fldU(isnan(fldU))=0; fldV(isnan(fldV))=0;
19
20	gforget	1.7
21			dxg=mk3D(mygrid.DXG,fldU); dyg=mk3D(mygrid.DYG,fldU);
22			if size(fldU.f1,3)==length(mygrid.DRF); drf=mk3D(mygrid.DRF,fldU); else; drf=fldU; drf(:)=1; end;
23			facW=drf; facW(:)=1; facS=facW;
24			for ii=1:length(list_factors);
25			tmp1=list_factors{ii};
26			if strcmp(tmp1,'dh'); facW=facW.dyg; facS=facS.dxg;
27			elseif strcmp(tmp1,'dz'); facW=facW.drf; facS=facS.drf;
28			elseif strcmp(tmp1,'hfac'); facW=facW.mygrid.hFacW; facS=facS.mygrid.hFacS;
29			elseif isempty(tmp1); 1;
30	gforget	1.8	else; fprintf('error in calc_barostream : non supported factor\n'); return;
31	gforget	1.7	end;
32			end;
33
34	gforget	1.1	for k4=1:n4;
35	gforget	1.7	fldU(:,:,:,k4)=fldU(:,:,:,k4).*facW;
36			fldV(:,:,:,k4)=fldV(:,:,:,k4).*facS;
37	gforget	1.1	end;
38
39	gforget	1.3	%apply mask:
40			fldU=sum(fldU,3).*mygrid.mskW(:,:,1);
41			fldV=sum(fldV,3).*mygrid.mskS(:,:,1);
42			%take out thedivergent part of the flow:
43	gforget	1.7	if noDiv;
44	gforget	1.3	[fldUdiv,fldVdiv,fldDivPot]=diffsmooth2D_div_inv(fldU,fldV);
45			fldU=fldU-fldUdiv; fldV=fldV-fldVdiv;
46			end;
47	gforget	1.1
48	gforget	1.2	%1) compute streamfunction face by face:
49			[fldU,fldV]=exch_UV(fldU,fldV); fldU(isnan(fldU))=0; fldV(isnan(fldV))=0;
50			tmp1=cumsum(fldV,1); for iF=1:fldU.nFaces; tmp2=tmp1{iF}; tmp1{iF}=[zeros(1,size(tmp2,2));tmp2]; end;
51			tmp2=diff(tmp1,1,2)+fldU; tmp3=cumsum(mean(tmp2,1));
52	gforget	1.3	%to check divergence implied errors:
53			%figure; for iF=1:fldU.nFaces; subplot(3,2,iF); plot(std(tmp2{iF},0,1)); end;
54	gforget	1.2	for iF=1:fldU.nFaces; tmp2=tmp1{iF}; tmp1{iF}=tmp2-ones(size(tmp2,1),1)*[0 tmp3{iF}]; end;
55			bf_step1=tmp1;
56
57	gforget	1.3	if fldU.nFaces==1;
58	gforget	1.4	bf_step2=bf_step1;
59	gforget	1.3	else;
60	gforget	1.2	%2) match edges:
61			%... set face number field
62			TMP1=tmp1; for iF=1:TMP1.nFaces; TMP1{iF}(:)=iF; end;
63			TMP2=exch_T_N(TMP1);%!!! this is a trick, since TMP1 is (n+1 X n+1) and loc. at vorticity points
64			tmp1=bf_step1;
65			for iF=1:fldU.nFaces-1;
66			tmp2=exch_T_N(tmp1);%!!! same trick
67			tmp3=tmp2{iF+1}; tmp3(3:end-2,3:end-2)=NaN;%mask out interior points
68			TMP3=TMP2{iF+1}; tmp3(find(TMP3>iF+1))=NaN;%mask out edges points coming from unadjusted faces
69			tmp3(:,1)=tmp3(:,1)-tmp3(:,2); tmp3(:,end)=tmp3(:,end)-tmp3(:,end-1);%compare edge points
70			tmp3(1,:)=tmp3(1,:)-tmp3(2,:); tmp3(end,:)=tmp3(end,:)-tmp3(end-1,:);%compare edge points
71			tmp3(2:end-1,2:end-1)=NaN;%mask out remaining interior points
72			tmp1{iF+1}=tmp1{iF+1}+nanmedian(tmp3(:));%adjust the face data
73	gforget	1.1	end;
74	gforget	1.2	bf_step2=tmp1;
75	gforget	1.4	end;
76	gforget	1.2	%3) put streamfunction at cell center
77			tmp1=bf_step2;
78			tmp2=tmp1; for iF=1:tmp1.nFaces; tmp3=tmp2{iF}; tmp3=(tmp3(:,1:end-1)+tmp3(:,2:end))/2;
79			tmp3=(tmp3(1:end-1,:)+tmp3(2:end,:))/2; tmp2{iF}=tmp3; end;
80			bf_step3=tmp2;
81			%4) set 0 on land on average:
82			tmp1=convert2array(bf_step3);
83			tmp2=convert2array(mygrid.mskC(:,:,1)); tmp2=find(isnan(tmp2)&~isnan(tmp1));
84			tmp2=median(tmp1(tmp2)); if isnan(tmp2); tmp2=nanmedian(tmp1(:)); end;
85			bf_step4=(bf_step3-tmp2).*mygrid.mskC(:,:,1);
86	gforget	1.1
87	gforget	1.2	%5) return the result:
88			fldBAR=bf_step4;
89	gforget	1.1
90	gforget	1.7	%convert to Sv and change sign:
91			fldBAR=1e-6*fldBAR;
92