| 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,ncx*nc,1); dyg=reshape(dyg,ncx*nc,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'],[ncx*nc 3],1,'real*4','b'); |
| 48 |
mskBs=rdda([rac,'maskS_bas.bin'],[ncx*nc 3],1,'real*4','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,ncx*nc,nr,Nit); vv=reshape(vv,ncx*nc,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 |
Parent Directory
| 
Revision Log
| 
Revision Graph