matlab_class/gcmfaces_calc/calc_bolus.m

function [bolusU,bolusV,bolusW]=calc_bolus(GM_PsiX,GM_PsiY);
%object:    compute bolus velocty field (bolusU,bolusV,bolusW)
%           from gm streamfunction (GM_PsiX,GM_PsiY).
%input:     GM_PsiX,GM_PsiY
%output:    bolusU,bolusV,bolusW

gcmfaces_global;
nr=length(mygrid.RC);

%compute bolus velocity:
GM_PsiX(:,:,nr+1)=0*GM_PsiX(:,:,end);
GM_PsiY(:,:,nr+1)=0*GM_PsiY(:,:,end);
bolusU=0*mygrid.hFacW;
bolusV=0*mygrid.hFacS;
for k=1:nr;
    bolusU(:,:,k)=(GM_PsiX(:,:,k+1)-GM_PsiX(:,:,k))/mygrid.DRF(k);
    bolusV(:,:,k)=(GM_PsiY(:,:,k+1)-GM_PsiY(:,:,k))/mygrid.DRF(k);
end;
bolusU=bolusU.*(~isnan(mygrid.mskW));
bolusV=bolusV.*(~isnan(mygrid.mskS));

%and its vertical part
%   (seems correct, leading to 0 divergence)
tmp_x=GM_PsiX(:,:,1:nr).*repmat(mygrid.DYG,[1 1 nr]);
tmp_y=GM_PsiY(:,:,1:nr).*repmat(mygrid.DXG,[1 1 nr]);
[tmp_x,tmp_y]=exch_UV(tmp_x,tmp_y);
tmp_a=repmat(mygrid.RAC,[1 1 nr]);
tmp_w=0*tmp_x;
for iFace=1:mygrid.nFaces;
    tmp_w{iFace}=( tmp_x{iFace}(2:end,:,:)-tmp_x{iFace}(1:end-1,:,:) )+...
        ( tmp_y{iFace}(:,2:end,:)-tmp_y{iFace}(:,1:end-1,:) );
    tmp_w{iFace}=tmp_w{iFace}./tmp_a{iFace};
end;
bolusW=tmp_w.*(~isnan(mygrid.mskC));
1	gforget	1.1	function [bolusU,bolusV,bolusW]=calc_bolus(GM_PsiX,GM_PsiY);
2			%object: compute bolus velocty field (bolusU,bolusV,bolusW)
3			% from gm streamfunction (GM_PsiX,GM_PsiY).
4			%input: GM_PsiX,GM_PsiY
5			%output: bolusU,bolusV,bolusW
6
7			gcmfaces_global;
8			nr=length(mygrid.RC);
9
10			%compute bolus velocity:
11			GM_PsiX(:,:,nr+1)=0*GM_PsiX(:,:,end);
12			GM_PsiY(:,:,nr+1)=0*GM_PsiY(:,:,end);
13			bolusU=0*mygrid.hFacW;
14			bolusV=0*mygrid.hFacS;
15			for k=1:nr;
16			bolusU(:,:,k)=(GM_PsiX(:,:,k+1)-GM_PsiX(:,:,k))/mygrid.DRF(k);
17			bolusV(:,:,k)=(GM_PsiY(:,:,k+1)-GM_PsiY(:,:,k))/mygrid.DRF(k);
18			end;
19			bolusU=bolusU.*(~isnan(mygrid.mskW));
20			bolusV=bolusV.*(~isnan(mygrid.mskS));
21
22			%and its vertical part
23			% (seems correct, leading to 0 divergence)
24			tmp_x=GM_PsiX(:,:,1:nr).*repmat(mygrid.DYG,[1 1 nr]);
25			tmp_y=GM_PsiY(:,:,1:nr).*repmat(mygrid.DXG,[1 1 nr]);
26			[tmp_x,tmp_y]=exch_UV(tmp_x,tmp_y);
27			tmp_a=repmat(mygrid.RAC,[1 1 nr]);
28			tmp_w=0*tmp_x;
29			for iFace=1:mygrid.nFaces;
30			tmp_w{iFace}=( tmp_x{iFace}(2:end,:,:)-tmp_x{iFace}(1:end-1,:,:) )+...
31			( tmp_y{iFace}(:,2:end,:)-tmp_y{iFace}(:,1:end-1,:) );
32			tmp_w{iFace}=tmp_w{iFace}./tmp_a{iFace};
33			end;
34			bolusW=tmp_w.*(~isnan(mygrid.mskC));