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);

%0.1) compute streamfunction mask:
if noDiv;
 mskU=1*~isnan(fldU);
 mskV=1*~isnan(fldV);
 [mskU,mskV]=exch_UV(mskU,mskV);
 mskU=abs(mskU); mskV=abs(mskV);
 mskBF=0*exch_T_N(mygrid.RAC);
 for iF=1:fldU.nFaces;
  tmp2=mskV{iF};
  tmp3a=[ones(1,size(tmp2,2));tmp2];
  tmp3b=[tmp2;ones(1,size(tmp2,2))];
  mskBF{iF}=tmp3a.*tmp3b;
 end;
end;

%take out the divergent 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;

%1.1) reset one land value per face to 0:
if noDiv;
 for iF=1:fldU.nFaces;
  tmpA=tmp1{iF};
  tmpB=mskBF{iF};
  [ii,jj]=find(tmpB==0);
  ii=ii(1); jj=jj(1);
  tmpA(:,jj)=tmpA(:,jj)-tmpB(ii,jj);
  for kk=1:jj-1; 
    tmpC=tmpA(:,kk+1)-tmpA(:,kk);
    tmpD=fldU{iF}(:,kk);
    tmpE=nanmedian(tmpC+tmpD);
    tmpA(:,kk)=tmpA(:,kk)-tmpE;
  end;
  for kk=jj+1:size(tmpA,2);
    tmpC=tmpA(:,kk)-tmpA(:,kk-1);
    tmpD=fldU{iF}(:,kk-1);
    tmpE=nanmedian(tmpC+tmpD);
    tmpA(:,kk)=tmpA(:,kk)-tmpE;
  end;
  tmp1{iF}=tmpA;
 end;
end;

%1.2) compute divergent part of the flow on average line by line:
if ~noDiv;
  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;
  %subtract from streamfunction:
  for iF=1:fldU.nFaces; tmp2=tmp1{iF}; tmp1{iF}=tmp2-ones(size(tmp2,1),1)*[0 tmp3{iF}]; end;
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 fist land point:
tmp1=convert2vector(bf_step3,'new'); 
tmp2=convert2vector(mygrid.mskC(:,:,1),'new');
tmp2=find(isnan(tmp2)&~isnan(tmp1)); 
%
%tmp2=median(tmp1(tmp2));%the original method
%tmp2=tmp1(tmp2(1));%a point in Antarctica (in LLC90 at least)
%closest land point closest to Boston:
tmp_lon=convert2vector(mygrid.XC,'new'); tmp_lon=tmp_lon(tmp2);
tmp_lat=convert2vector(mygrid.YC,'new'); tmp_lat=tmp_lat(tmp2);
tmp_dis=(tmp_lat-42.3601).^2+(tmp_lon-71.0589).^2;
tmp2=tmp1(tmp2(find(tmp_dis==min(tmp_dis))));
%
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	dxg=mk3D(mygrid.DXG,fldU); dyg=mk3D(mygrid.DYG,fldU);
21			if size(fldU.f1,3)==length(mygrid.DRF); drf=mk3D(mygrid.DRF,fldU); else; drf=fldU; drf(:)=1; end;
22			facW=drf; facW(:)=1; facS=facW;
23			for ii=1:length(list_factors);
24			tmp1=list_factors{ii};
25			if strcmp(tmp1,'dh'); facW=facW.dyg; facS=facS.dxg;
26			elseif strcmp(tmp1,'dz'); facW=facW.drf; facS=facS.drf;
27			elseif strcmp(tmp1,'hfac'); facW=facW.mygrid.hFacW; facS=facS.mygrid.hFacS;
28			elseif isempty(tmp1); 1;
29	gforget	1.8	else; fprintf('error in calc_barostream : non supported factor\n'); return;
30	gforget	1.7	end;
31			end;
32
33	gforget	1.1	for k4=1:n4;
34	gforget	1.7	fldU(:,:,:,k4)=fldU(:,:,:,k4).*facW;
35			fldV(:,:,:,k4)=fldV(:,:,:,k4).*facS;
36	gforget	1.1	end;
37
38	gforget	1.3	%apply mask:
39			fldU=sum(fldU,3).*mygrid.mskW(:,:,1);
40			fldV=sum(fldV,3).*mygrid.mskS(:,:,1);
41	gforget	1.9
42			%0.1) compute streamfunction mask:
43			if noDiv;
44			mskU=1*~isnan(fldU);
45			mskV=1*~isnan(fldV);
46			[mskU,mskV]=exch_UV(mskU,mskV);
47			mskU=abs(mskU); mskV=abs(mskV);
48			mskBF=0*exch_T_N(mygrid.RAC);
49			for iF=1:fldU.nFaces;
50			tmp2=mskV{iF};
51			tmp3a=[ones(1,size(tmp2,2));tmp2];
52			tmp3b=[tmp2;ones(1,size(tmp2,2))];
53			mskBF{iF}=tmp3a.*tmp3b;
54			end;
55			end;
56
57			%take out the divergent part of the flow:
58	gforget	1.7	if noDiv;
59	gforget	1.3	[fldUdiv,fldVdiv,fldDivPot]=diffsmooth2D_div_inv(fldU,fldV);
60			fldU=fldU-fldUdiv; fldV=fldV-fldVdiv;
61			end;
62	gforget	1.1
63	gforget	1.2	%1) compute streamfunction face by face:
64			[fldU,fldV]=exch_UV(fldU,fldV); fldU(isnan(fldU))=0; fldV(isnan(fldV))=0;
65			tmp1=cumsum(fldV,1); for iF=1:fldU.nFaces; tmp2=tmp1{iF}; tmp1{iF}=[zeros(1,size(tmp2,2));tmp2]; end;
66	gforget	1.9
67			%1.1) reset one land value per face to 0:
68			if noDiv;
69			for iF=1:fldU.nFaces;
70			tmpA=tmp1{iF};
71			tmpB=mskBF{iF};
72			[ii,jj]=find(tmpB==0);
73			ii=ii(1); jj=jj(1);
74			tmpA(:,jj)=tmpA(:,jj)-tmpB(ii,jj);
75			for kk=1:jj-1;
76			tmpC=tmpA(:,kk+1)-tmpA(:,kk);
77			tmpD=fldU{iF}(:,kk);
78			tmpE=nanmedian(tmpC+tmpD);
79			tmpA(:,kk)=tmpA(:,kk)-tmpE;
80			end;
81			for kk=jj+1:size(tmpA,2);
82			tmpC=tmpA(:,kk)-tmpA(:,kk-1);
83			tmpD=fldU{iF}(:,kk-1);
84			tmpE=nanmedian(tmpC+tmpD);
85			tmpA(:,kk)=tmpA(:,kk)-tmpE;
86			end;
87			tmp1{iF}=tmpA;
88			end;
89			end;
90
91			%1.2) compute divergent part of the flow on average line by line:
92			if ~noDiv;
93			tmp2=diff(tmp1,1,2)+fldU; tmp3=cumsum(mean(tmp2,1));
94			% to check divergence implied errors:
95			% figure; for iF=1:fldU.nFaces; subplot(3,2,iF); plot(std(tmp2{iF},0,1)); end;
96			%subtract from streamfunction:
97			for iF=1:fldU.nFaces; tmp2=tmp1{iF}; tmp1{iF}=tmp2-ones(size(tmp2,1),1)*[0 tmp3{iF}]; end;
98			end;
99
100	gforget	1.2	bf_step1=tmp1;
101
102	gforget	1.3	if fldU.nFaces==1;
103	gforget	1.4	bf_step2=bf_step1;
104	gforget	1.3	else;
105	gforget	1.2	%2) match edges:
106			%... set face number field
107			TMP1=tmp1; for iF=1:TMP1.nFaces; TMP1{iF}(:)=iF; end;
108			TMP2=exch_T_N(TMP1);%!!! this is a trick, since TMP1 is (n+1 X n+1) and loc. at vorticity points
109			tmp1=bf_step1;
110			for iF=1:fldU.nFaces-1;
111			tmp2=exch_T_N(tmp1);%!!! same trick
112			tmp3=tmp2{iF+1}; tmp3(3:end-2,3:end-2)=NaN;%mask out interior points
113			TMP3=TMP2{iF+1}; tmp3(find(TMP3>iF+1))=NaN;%mask out edges points coming from unadjusted faces
114			tmp3(:,1)=tmp3(:,1)-tmp3(:,2); tmp3(:,end)=tmp3(:,end)-tmp3(:,end-1);%compare edge points
115			tmp3(1,:)=tmp3(1,:)-tmp3(2,:); tmp3(end,:)=tmp3(end,:)-tmp3(end-1,:);%compare edge points
116			tmp3(2:end-1,2:end-1)=NaN;%mask out remaining interior points
117			tmp1{iF+1}=tmp1{iF+1}+nanmedian(tmp3(:));%adjust the face data
118	gforget	1.1	end;
119	gforget	1.2	bf_step2=tmp1;
120	gforget	1.4	end;
121	gforget	1.9
122	gforget	1.2	%3) put streamfunction at cell center
123			tmp1=bf_step2;
124			tmp2=tmp1; for iF=1:tmp1.nFaces; tmp3=tmp2{iF}; tmp3=(tmp3(:,1:end-1)+tmp3(:,2:end))/2;
125			tmp3=(tmp3(1:end-1,:)+tmp3(2:end,:))/2; tmp2{iF}=tmp3; end;
126			bf_step3=tmp2;
127	gforget	1.10	%4) set 0 on fist land point:
128			tmp1=convert2vector(bf_step3,'new');
129			tmp2=convert2vector(mygrid.mskC(:,:,1),'new');
130			tmp2=find(isnan(tmp2)&~isnan(tmp1));
131	gforget	1.11	%
132			%tmp2=median(tmp1(tmp2));%the original method
133			%tmp2=tmp1(tmp2(1));%a point in Antarctica (in LLC90 at least)
134			%closest land point closest to Boston:
135			tmp_lon=convert2vector(mygrid.XC,'new'); tmp_lon=tmp_lon(tmp2);
136			tmp_lat=convert2vector(mygrid.YC,'new'); tmp_lat=tmp_lat(tmp2);
137			tmp_dis=(tmp_lat-42.3601).^2+(tmp_lon-71.0589).^2;
138			tmp2=tmp1(tmp2(find(tmp_dis==min(tmp_dis))));
139			%
140	gforget	1.2	bf_step4=(bf_step3-tmp2).*mygrid.mskC(:,:,1);
141	gforget	1.1
142	gforget	1.2	%5) return the result:
143			fldBAR=bf_step4;
144	gforget	1.1
145	gforget	1.7	%convert to Sv and change sign:
146			fldBAR=1e-6*fldBAR;
147
148	gforget	1.9