Diagnostics/DiagUtility/calc_PsiCube.m

function [psi,mskG,ylat]=calc_PsiCube(delM,uu,vv,dxg,dyg,hFacW,hFacS,nBas,dBug);
% [psi,mskG,ylat]=calc_PsiCube(delM,uu,vv,[hFacW,hFacS],[nBas],[dBug]);
%- IMPORTANT: must multiply (u,v) by hFacW,S BEFORE using this script !
%   (so that it can be used in r* coordinate with (h*u,hv)_timeAv in input)
% delM= -delP/g for atmos ; =delZ for ocean (delR)

%- Units: dx,dy /1e6 ; delR /1e3 [hPa] ; psi in 10^9 kg/s
%  Atmos in p : use g=9.81 ; ocean in z : use g=-1;

krd=1; kMsep=1; jprt=0;
nr=length(delM);
Tprt=0;

if (nargin < 9), dBug=0; end
if (nargin < 8), nBas=0; end
if (nargin < 6), kfac=0;
else kfac=1; end;

if Tprt, TimeT0=clock; end

if krd > 0,
%rac='/home/jmc/grid_cs32/' ;
 rac='/u/u2/jmc/grid_cs32/' ;
%-- broken lines file, 1rst version ; 2nd version (including latitude strip):

%- load: bkl_Ylat,bkl_Npts,bkl_Flg,bkl_Iuv,bkl_Juv,bkl_Xsg,bkl_Ysg
% bk_lineF=[rac,'isoLat_cube32_59'];
% load(bk_lineF);
% bkl_IJuv=bkl_Iuv+ncx*(bkl_Juv-1);

%- load: bkl_Ylat, bkl_Npts, bkl_Flg, bkl_IJuv, bkl_Xsg, bkl_Ysg, bkl_Zon
bk_lineF=[rac,'isoLat_cs32_59.mat'];
load('isoLat_cs32_59.mat');
if dBug, fprintf([' load bk_line definition from: ',bk_lineF]); end

%- load the grid dx,dy , convert to 10^6 m :
 %dxg=rdmds([rac,'DXG']);
 dxg=dxg*1.e-6;
 %dyg=rdmds([rac,'DYG']);
 dyg=dyg*1.e-6;
 ncx=size(dxg,1); nc=size(dxg,2);
 dxg=reshape(dxg,ncx*nc,1); dyg=reshape(dyg,ncx*nc,1);
 if dBug, fprintf(' AND dxg,dyg'); end
          
if nBas > 0,
%- load Ocean Basin mask (& separation line):
 mskBw=rdda([rac,'maskW_bas.bin'],[ncx*nc 3],1,'real*4','b');
 mskBs=rdda([rac,'maskS_bas.bin'],[ncx*nc 3],1,'real*4','b');
 if nBas==2, 
  mskBw(:,2)=mskBw(:,2)+mskBw(:,3);
  mskBs(:,2)=mskBs(:,2)+mskBs(:,3);
  mskBw=min(1,mskBw); mskBs=min(1,mskBs);
 end
%- load: np_Sep, ij_Sep, tp_Sep:
 sep_lineF=[rac,'sepBas_cs32_60'];
 load(sep_lineF);
 if dBug, fprintf([' + bassin mask & Sep.line:',sep_lineF]); end
end
 if dBug, fprintf('\n'); end
end

if Tprt, TimeT1=clock; end

%- compute the horizontal transport ut,vt :
if length(size(uu)) < 4, Nit=1; else Nit=size(uu,4); end;
uu=reshape(uu,ncx*nc,nr,Nit); vv=reshape(vv,ncx*nc,nr,Nit);
ydim=length(bkl_Ylat); ylat=bkl_Ylat;
psi=zeros(ydim+2,nr+1,1+nBas,Nit); 
mskZ=zeros(ydim+2,nr+1,1+nBas); % = mask of Psi 
mskV=zeros(ydim+2,nr,1+nBas);   % = mask of the Merid.Transport
mskG=zeros(ydim+1,nr,1+nBas);   % = mask of the Ground

%- define ufac,vfac for each bassin:
ufac=zeros([size(bkl_Flg) 1+nBas]);
vfac=zeros([size(bkl_Flg) 1+nBas]);
ufac(:,:,1)=rem(bkl_Flg,2) ; vfac(:,:,1)=fix(bkl_Flg/2) ;
for jl=1:ydim,
  ie=bkl_Npts(jl);
  for b=1:nBas, 
   ufac(1:ie,jl,1+b)=mskBw(bkl_IJuv(1:ie,jl),b).*ufac(1:ie,jl,1);
   vfac(1:ie,jl,1+b)=mskBs(bkl_IJuv(1:ie,jl),b).*vfac(1:ie,jl,1);
 end;
end;

%- compute transport ; integrate folowing broken-lines
for nt=1:Nit, 
 for k=nr:-1:1,

  ut=dyg.*uu(:,k,nt); 
  vt=dxg.*vv(:,k,nt);
  for jl=1:ydim,
   if jl == jprt, fprintf(' jl= %2i , lat= %8.3f , Npts= %3i :\n', ...
                            jl,ylat(jl),bkl_Npts(jl)); end
   ie=bkl_Npts(jl);
   for b=1:1+nBas, 
    vz=sum( ufac(1:ie,jl,b).*ut(bkl_IJuv(1:ie,jl)) ...
           +vfac(1:ie,jl,b).*vt(bkl_IJuv(1:ie,jl)) );
    psi(jl+1,k,b,nt)=psi(jl+1,k+1,b,nt) - delM(k)*vz ;
   end
  end

 end
end
   
if Tprt, TimeT2=clock; end

 %- compute the mask :
if kfac == 1, 
  hFacW=reshape(hFacW,ncx*nc,nr);
  hFacS=reshape(hFacS,ncx*nc,nr);
  ufac=abs(ufac) ; vfac=abs(vfac);
 for jl=1:ydim, 
  ie=bkl_Npts(jl); 
  hw=zeros(ie,nr); hs=zeros(ie,nr);
  hw=hFacW(bkl_IJuv(1:ie,jl),:);
  hs=hFacS(bkl_IJuv(1:ie,jl),:);
  for b=1:1+nBas,
   for k=1:nr,
%   for ii=1:bkl_Npts(jl);
%    ij=bkl_IJuv(ii,jl);
%    mskV(jl+1,k,b)=mskV(jl+1,k,b)+ufac(ii,jl,b)*hFacW(ij,k)+vfac(ii,jl,b)*hFacS(ij,k);
%   end ;
     tmpv=ufac(1:ie,jl,b).*hw(:,k)+vfac(1:ie,jl,b).*hs(:,k);
     mskV(jl+1,k,b)=mskV(jl+1,k,b)+max(tmpv);
   end ;
  end ;
 end
 mskV=ceil(mskV); mskV=min(1,mskV);
%- build the real mask (=mskG, ground) used to draw the continent with "surf":
%   position=centered , dim= ydim+1 x nr
 mskG=mskV(1:ydim+1,:,:)+mskV(2:ydim+2,:,:); mskG=min(1,mskG);

 if Tprt, TimeT3=clock; end

if kMsep & nBas > 0,
 mskW=1+min(1,ceil(hFacW));
 mskS=1+min(1,ceil(hFacS));
 for b=1:nBas,
  bs=b; b1=1+bs; b2=2+rem(bs,nBas);
  if nBas == 2, bs=b+b-1; b1=2; b2=3 ; end
  for j=1:ydim+1,
   for i=1:np_Sep(bs,j),
    ij=ij_Sep(bs,j,i); typ=abs(tp_Sep(bs,j,i));
    if typ == 1, 
     mskG(j,:,b1)=mskG(j,:,b1).*mskW(ij,:);
     mskG(j,:,b2)=mskG(j,:,b2).*mskW(ij,:);
    elseif typ == 2,
     mskG(j,:,b1)=mskG(j,:,b1).*mskS(ij,:);
     mskG(j,:,b2)=mskG(j,:,b2).*mskS(ij,:);
    end
   end
  end
 end
 mskG=min(2,mskG);
else
 if Tprt, TimeT3=clock; end
end

if Tprt, TimeT4=clock; end

%- to keep psi=0 on top & bottom
 mskZ(:,[2:nr+1],:)=mskV; 
 mskZ(:,[1:nr],:)=mskZ(:,[1:nr],:)+mskV;
%- to keep psi=0 on lateral boundaries :
 mskZ([1:ydim],:,:)=mskZ([1:ydim],:,:)+mskZ([2:ydim+1],:,:);
 mskZ([2:ydim+1],:,:)=mskZ([2:ydim+1],:,:)+mskZ([3:ydim+2],:,:);
 mskZ=ceil(mskZ); mskZ=min(1,mskZ);
 if kMsep & nBas > 0,
  mskM=zeros(ydim+2,nr,1+nBas); mskM(2:ydim+2,:,:)=min(2-mskG,1); 
  mskM(1:ydim+1,:,:)=mskM(1:ydim+1,:,:)+mskM(2:ydim+2,:,:);
  mskZ(:,1:nr,:)=min(mskZ(:,1:nr,:),mskM); 
 end
%- apply the mask (and remove dim = 1) :
 if Nit == 1,
  psi=squeeze(psi); mskV=squeeze(mskV); mskZ=squeeze(mskZ);
  psi( find(mskZ==0) )=NaN ;
 else
  for nt=1:Nit,
    psi1=psi(:,:,:,nt); psi1( find(mskZ==0) )=NaN ; psi(:,:,:,nt)=psi1;
  end;
  if nBas < 1, psi=squeeze(psi); mskV=squeeze(mskV); mskZ=squeeze(mskZ); end
 end
else
 if Tprt, TimeT3=TimeT2; TimeT4=TimeT2; end
 if nBas < 1 | Nit == 1, psi=squeeze(psi); end
end;
%----------------- 

if Tprt, TimeT5=clock;
 fprintf([' ---- Load, Comp.1,2,3,4 Total time Psi:', ...
          ' %6.3f %6.3f %6.3f %6.3f %6.3f %9.6f \n'],...
         etime(TimeT1,TimeT0), etime(TimeT2,TimeT1), ...
         etime(TimeT3,TimeT2), etime(TimeT4,TimeT3), ...
         etime(TimeT5,TimeT4), etime(TimeT5,TimeT0) );
end

return
1	enderton	1.1	function [psi,mskG,ylat]=calc_PsiCube(delM,uu,vv,dxg,dyg,hFacW,hFacS,nBas,dBug);
2			% [psi,mskG,ylat]=calc_PsiCube(delM,uu,vv,[hFacW,hFacS],[nBas],[dBug]);
3			%- IMPORTANT: must multiply (u,v) by hFacW,S BEFORE using this script !
4			% (so that it can be used in r* coordinate with (h*u,hv)_timeAv in input)
5			% delM= -delP/g for atmos ; =delZ for ocean (delR)
6
7			%- Units: dx,dy /1e6 ; delR /1e3 [hPa] ; psi in 10^9 kg/s
8			% Atmos in p : use g=9.81 ; ocean in z : use g=-1;
9
10			krd=1; kMsep=1; jprt=0;
11			nr=length(delM);
12			Tprt=0;
13
14			if (nargin < 9), dBug=0; end
15			if (nargin < 8), nBas=0; end
16			if (nargin < 6), kfac=0;
17			else kfac=1; end;
18
19			if Tprt, TimeT0=clock; end
20
21			if krd > 0,
22			%rac='/home/jmc/grid_cs32/' ;
23			rac='/u/u2/jmc/grid_cs32/' ;
24			%-- broken lines file, 1rst version ; 2nd version (including latitude strip):
25
26			%- load: bkl_Ylat,bkl_Npts,bkl_Flg,bkl_Iuv,bkl_Juv,bkl_Xsg,bkl_Ysg
27			% bk_lineF=[rac,'isoLat_cube32_59'];
28			% load(bk_lineF);
29			% bkl_IJuv=bkl_Iuv+ncx*(bkl_Juv-1);
30
31			%- load: bkl_Ylat, bkl_Npts, bkl_Flg, bkl_IJuv, bkl_Xsg, bkl_Ysg, bkl_Zon
32			bk_lineF=[rac,'isoLat_cs32_59.mat'];
33			load('isoLat_cs32_59.mat');
34			if dBug, fprintf([' load bk_line definition from: ',bk_lineF]); end
35
36			%- load the grid dx,dy , convert to 10^6 m :
37			%dxg=rdmds([rac,'DXG']);
38			dxg=dxg*1.e-6;
39			%dyg=rdmds([rac,'DYG']);
40			dyg=dyg*1.e-6;
41			ncx=size(dxg,1); nc=size(dxg,2);
42			dxg=reshape(dxg,ncxnc,1); dyg=reshape(dyg,ncxnc,1);
43			if dBug, fprintf(' AND dxg,dyg'); end
44
45			if nBas > 0,
46			%- load Ocean Basin mask (& separation line):
47			mskBw=rdda([rac,'maskW_bas.bin'],[ncxnc 3],1,'real4','b');
48			mskBs=rdda([rac,'maskS_bas.bin'],[ncxnc 3],1,'real4','b');
49			if nBas==2,
50			mskBw(:,2)=mskBw(:,2)+mskBw(:,3);
51			mskBs(:,2)=mskBs(:,2)+mskBs(:,3);
52			mskBw=min(1,mskBw); mskBs=min(1,mskBs);
53			end
54			%- load: np_Sep, ij_Sep, tp_Sep:
55			sep_lineF=[rac,'sepBas_cs32_60'];
56			load(sep_lineF);
57			if dBug, fprintf([' + bassin mask & Sep.line:',sep_lineF]); end
58			end
59			if dBug, fprintf('\n'); end
60			end
61
62			if Tprt, TimeT1=clock; end
63
64			%- compute the horizontal transport ut,vt :
65			if length(size(uu)) < 4, Nit=1; else Nit=size(uu,4); end;
66			uu=reshape(uu,ncxnc,nr,Nit); vv=reshape(vv,ncxnc,nr,Nit);
67			ydim=length(bkl_Ylat); ylat=bkl_Ylat;
68			psi=zeros(ydim+2,nr+1,1+nBas,Nit);
69			mskZ=zeros(ydim+2,nr+1,1+nBas); % = mask of Psi
70			mskV=zeros(ydim+2,nr,1+nBas); % = mask of the Merid.Transport
71			mskG=zeros(ydim+1,nr,1+nBas); % = mask of the Ground
72
73			%- define ufac,vfac for each bassin:
74			ufac=zeros([size(bkl_Flg) 1+nBas]);
75			vfac=zeros([size(bkl_Flg) 1+nBas]);
76			ufac(:,:,1)=rem(bkl_Flg,2) ; vfac(:,:,1)=fix(bkl_Flg/2) ;
77			for jl=1:ydim,
78			ie=bkl_Npts(jl);
79			for b=1:nBas,
80			ufac(1:ie,jl,1+b)=mskBw(bkl_IJuv(1:ie,jl),b).*ufac(1:ie,jl,1);
81			vfac(1:ie,jl,1+b)=mskBs(bkl_IJuv(1:ie,jl),b).*vfac(1:ie,jl,1);
82			end;
83			end;
84
85			%- compute transport ; integrate folowing broken-lines
86			for nt=1:Nit,
87			for k=nr:-1:1,
88
89			ut=dyg.*uu(:,k,nt);
90			vt=dxg.*vv(:,k,nt);
91			for jl=1:ydim,
92			if jl == jprt, fprintf(' jl= %2i , lat= %8.3f , Npts= %3i :\n', ...
93			jl,ylat(jl),bkl_Npts(jl)); end
94			ie=bkl_Npts(jl);
95			for b=1:1+nBas,
96			vz=sum( ufac(1:ie,jl,b).*ut(bkl_IJuv(1:ie,jl)) ...
97			+vfac(1:ie,jl,b).*vt(bkl_IJuv(1:ie,jl)) );
98			psi(jl+1,k,b,nt)=psi(jl+1,k+1,b,nt) - delM(k)*vz ;
99			end
100			end
101
102			end
103			end
104
105			if Tprt, TimeT2=clock; end
106
107			%- compute the mask :
108			if kfac == 1,
109			hFacW=reshape(hFacW,ncx*nc,nr);
110			hFacS=reshape(hFacS,ncx*nc,nr);
111			ufac=abs(ufac) ; vfac=abs(vfac);
112			for jl=1:ydim,
113			ie=bkl_Npts(jl);
114			hw=zeros(ie,nr); hs=zeros(ie,nr);
115			hw=hFacW(bkl_IJuv(1:ie,jl),:);
116			hs=hFacS(bkl_IJuv(1:ie,jl),:);
117			for b=1:1+nBas,
118			for k=1:nr,
119			% for ii=1:bkl_Npts(jl);
120			% ij=bkl_IJuv(ii,jl);
121			% mskV(jl+1,k,b)=mskV(jl+1,k,b)+ufac(ii,jl,b)hFacW(ij,k)+vfac(ii,jl,b)hFacS(ij,k);
122			% end ;
123			tmpv=ufac(1:ie,jl,b).hw(:,k)+vfac(1:ie,jl,b).hs(:,k);
124			mskV(jl+1,k,b)=mskV(jl+1,k,b)+max(tmpv);
125			end ;
126			end ;
127			end
128			mskV=ceil(mskV); mskV=min(1,mskV);
129			%- build the real mask (=mskG, ground) used to draw the continent with "surf":
130			% position=centered , dim= ydim+1 x nr
131			mskG=mskV(1:ydim+1,:,:)+mskV(2:ydim+2,:,:); mskG=min(1,mskG);
132
133			if Tprt, TimeT3=clock; end
134
135			if kMsep & nBas > 0,
136			mskW=1+min(1,ceil(hFacW));
137			mskS=1+min(1,ceil(hFacS));
138			for b=1:nBas,
139			bs=b; b1=1+bs; b2=2+rem(bs,nBas);
140			if nBas == 2, bs=b+b-1; b1=2; b2=3 ; end
141			for j=1:ydim+1,
142			for i=1:np_Sep(bs,j),
143			ij=ij_Sep(bs,j,i); typ=abs(tp_Sep(bs,j,i));
144			if typ == 1,
145			mskG(j,:,b1)=mskG(j,:,b1).*mskW(ij,:);
146			mskG(j,:,b2)=mskG(j,:,b2).*mskW(ij,:);
147			elseif typ == 2,
148			mskG(j,:,b1)=mskG(j,:,b1).*mskS(ij,:);
149			mskG(j,:,b2)=mskG(j,:,b2).*mskS(ij,:);
150			end
151			end
152			end
153			end
154			mskG=min(2,mskG);
155			else
156			if Tprt, TimeT3=clock; end
157			end
158
159			if Tprt, TimeT4=clock; end
160
161			%- to keep psi=0 on top & bottom
162			mskZ(:,[2:nr+1],:)=mskV;
163			mskZ(:,[1:nr],:)=mskZ(:,[1:nr],:)+mskV;
164			%- to keep psi=0 on lateral boundaries :
165			mskZ([1:ydim],:,:)=mskZ([1:ydim],:,:)+mskZ([2:ydim+1],:,:);
166			mskZ([2:ydim+1],:,:)=mskZ([2:ydim+1],:,:)+mskZ([3:ydim+2],:,:);
167			mskZ=ceil(mskZ); mskZ=min(1,mskZ);
168			if kMsep & nBas > 0,
169			mskM=zeros(ydim+2,nr,1+nBas); mskM(2:ydim+2,:,:)=min(2-mskG,1);
170			mskM(1:ydim+1,:,:)=mskM(1:ydim+1,:,:)+mskM(2:ydim+2,:,:);
171			mskZ(:,1:nr,:)=min(mskZ(:,1:nr,:),mskM);
172			end
173			%- apply the mask (and remove dim = 1) :
174			if Nit == 1,
175			psi=squeeze(psi); mskV=squeeze(mskV); mskZ=squeeze(mskZ);
176			psi( find(mskZ==0) )=NaN ;
177			else
178			for nt=1:Nit,
179			psi1=psi(:,:,:,nt); psi1( find(mskZ==0) )=NaN ; psi(:,:,:,nt)=psi1;
180			end;
181			if nBas < 1, psi=squeeze(psi); mskV=squeeze(mskV); mskZ=squeeze(mskZ); end
182			end
183			else
184			if Tprt, TimeT3=TimeT2; TimeT4=TimeT2; end
185			if nBas < 1 \| Nit == 1, psi=squeeze(psi); end
186			end;
187			%-----------------
188
189			if Tprt, TimeT5=clock;
190			fprintf([' ---- Load, Comp.1,2,3,4 Total time Psi:', ...
191			' %6.3f %6.3f %6.3f %6.3f %6.3f %9.6f \n'],...
192			etime(TimeT1,TimeT0), etime(TimeT2,TimeT1), ...
193			etime(TimeT3,TimeT2), etime(TimeT4,TimeT3), ...
194			etime(TimeT5,TimeT4), etime(TimeT5,TimeT0) );
195			end
196
197			return