Warning, /utils/matlab/Graphix/GraphixUtility/calc_PsiCube.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit 497302ac on 2005-11-25 06:36:32 UTC
892fe9e2f3 Dani*0001 function [psi,mskG,ylat]=calc_PsiCube(delM,uu,vv,dxg,dyg,hFacW,hFacS,nBas,dBug);
497302ac9d Dani*0002
892fe9e2f3 Dani*0003 % [psi,mskG,ylat]=calc_PsiCube(delM,uu,vv,[hFacW,hFacS],[nBas],[dBug]);
0004 %- IMPORTANT: must multiply (u,v) by hFacW,S BEFORE using this script !
0005 % (so that it can be used in r* coordinate with (h*u,hv)_timeAv in input)
0006 % delM= -delP/g for atmos ; =delZ for ocean (delR)
0007
0008 %- Units: dx,dy /1e6 ; delR /1e3 [hPa] ; psi in 10^9 kg/s
0009 % Atmos in p : use g=9.81 ; ocean in z : use g=-1;
0010
0011 krd=1; kMsep=1; jprt=0;
0012 nr=length(delM);
0013 Tprt=0;
0014
0015 if (nargin < 9), dBug=0; end
0016 if (nargin < 8), nBas=0; end
0017 if (nargin < 6), kfac=0;
0018 else kfac=1; end;
0019
0020 if Tprt, TimeT0=clock; end
0021
0022 if krd > 0,
0023 %- load: bkl_Ylat, bkl_Npts, bkl_Flg, bkl_IJuv, bkl_Xsg, bkl_Ysg, bkl_Zon
0024 load('isoLat_cs32_59.mat');
0025
0026 %- load the grid dx,dy , convert to 10^6 m :
497302ac9d Dani*0027 dxg=dxg*1.e-6;
0028 dyg=dyg*1.e-6;
0029 ncx=size(dxg,1);
0030 nc=size(dxg,2);
0031 dxg=reshape(dxg,ncx*nc,1); dyg=reshape(dyg,ncx*nc,1);
0032 if dBug, fprintf(' AND dxg,dyg'); end
892fe9e2f3 Dani*0033
0034 if nBas > 0,
0035 %- load Ocean Basin mask (& separation line):
0036 mskBw=rdda([rac,'maskW_bas.bin'],[ncx*nc 3],1,'real*4','b');
0037 mskBs=rdda([rac,'maskS_bas.bin'],[ncx*nc 3],1,'real*4','b');
0038 if nBas==2,
0039 mskBw(:,2)=mskBw(:,2)+mskBw(:,3);
0040 mskBs(:,2)=mskBs(:,2)+mskBs(:,3);
0041 mskBw=min(1,mskBw); mskBs=min(1,mskBs);
0042 end
0043 %- load: np_Sep, ij_Sep, tp_Sep:
0044 sep_lineF=[rac,'sepBas_cs32_60'];
0045 load(sep_lineF);
0046 if dBug, fprintf([' + bassin mask & Sep.line:',sep_lineF]); end
0047 end
0048 if dBug, fprintf('\n'); end
0049 end
0050
0051 if Tprt, TimeT1=clock; end
0052
0053 %- compute the horizontal transport ut,vt :
0054 if length(size(uu)) < 4, Nit=1; else Nit=size(uu,4); end;
0055 uu=reshape(uu,ncx*nc,nr,Nit); vv=reshape(vv,ncx*nc,nr,Nit);
0056 ydim=length(bkl_Ylat); ylat=bkl_Ylat;
0057 psi=zeros(ydim+2,nr+1,1+nBas,Nit);
0058 mskZ=zeros(ydim+2,nr+1,1+nBas); % = mask of Psi
0059 mskV=zeros(ydim+2,nr,1+nBas); % = mask of the Merid.Transport
0060 mskG=zeros(ydim+1,nr,1+nBas); % = mask of the Ground
0061
0062 %- define ufac,vfac for each bassin:
0063 ufac=zeros([size(bkl_Flg) 1+nBas]);
0064 vfac=zeros([size(bkl_Flg) 1+nBas]);
0065 ufac(:,:,1)=rem(bkl_Flg,2) ; vfac(:,:,1)=fix(bkl_Flg/2) ;
0066 for jl=1:ydim,
0067 ie=bkl_Npts(jl);
0068 for b=1:nBas,
0069 ufac(1:ie,jl,1+b)=mskBw(bkl_IJuv(1:ie,jl),b).*ufac(1:ie,jl,1);
0070 vfac(1:ie,jl,1+b)=mskBs(bkl_IJuv(1:ie,jl),b).*vfac(1:ie,jl,1);
0071 end;
0072 end;
0073
0074 %- compute transport ; integrate folowing broken-lines
0075 for nt=1:Nit,
0076 for k=nr:-1:1,
0077
0078 ut=dyg.*uu(:,k,nt);
0079 vt=dxg.*vv(:,k,nt);
0080 for jl=1:ydim,
0081 if jl == jprt, fprintf(' jl= %2i , lat= %8.3f , Npts= %3i :\n', ...
0082 jl,ylat(jl),bkl_Npts(jl)); end
0083 ie=bkl_Npts(jl);
0084 for b=1:1+nBas,
0085 vz=sum( ufac(1:ie,jl,b).*ut(bkl_IJuv(1:ie,jl)) ...
0086 +vfac(1:ie,jl,b).*vt(bkl_IJuv(1:ie,jl)) );
0087 psi(jl+1,k,b,nt)=psi(jl+1,k+1,b,nt) - delM(k)*vz ;
0088 end
0089 end
0090
0091 end
0092 end
0093
0094 if Tprt, TimeT2=clock; end
0095
0096 %- compute the mask :
0097 if kfac == 1,
0098 hFacW=reshape(hFacW,ncx*nc,nr);
0099 hFacS=reshape(hFacS,ncx*nc,nr);
0100 ufac=abs(ufac) ; vfac=abs(vfac);
0101 for jl=1:ydim,
0102 ie=bkl_Npts(jl);
0103 hw=zeros(ie,nr); hs=zeros(ie,nr);
0104 hw=hFacW(bkl_IJuv(1:ie,jl),:);
0105 hs=hFacS(bkl_IJuv(1:ie,jl),:);
0106 for b=1:1+nBas,
0107 for k=1:nr,
0108 % for ii=1:bkl_Npts(jl);
0109 % ij=bkl_IJuv(ii,jl);
0110 % mskV(jl+1,k,b)=mskV(jl+1,k,b)+ufac(ii,jl,b)*hFacW(ij,k)+vfac(ii,jl,b)*hFacS(ij,k);
0111 % end ;
0112 tmpv=ufac(1:ie,jl,b).*hw(:,k)+vfac(1:ie,jl,b).*hs(:,k);
0113 mskV(jl+1,k,b)=mskV(jl+1,k,b)+max(tmpv);
0114 end ;
0115 end ;
0116 end
0117 mskV=ceil(mskV); mskV=min(1,mskV);
0118 %- build the real mask (=mskG, ground) used to draw the continent with "surf":
0119 % position=centered , dim= ydim+1 x nr
0120 mskG=mskV(1:ydim+1,:,:)+mskV(2:ydim+2,:,:); mskG=min(1,mskG);
0121
0122 if Tprt, TimeT3=clock; end
0123
0124 if kMsep & nBas > 0,
0125 mskW=1+min(1,ceil(hFacW));
0126 mskS=1+min(1,ceil(hFacS));
0127 for b=1:nBas,
0128 bs=b; b1=1+bs; b2=2+rem(bs,nBas);
0129 if nBas == 2, bs=b+b-1; b1=2; b2=3 ; end
0130 for j=1:ydim+1,
0131 for i=1:np_Sep(bs,j),
0132 ij=ij_Sep(bs,j,i); typ=abs(tp_Sep(bs,j,i));
0133 if typ == 1,
0134 mskG(j,:,b1)=mskG(j,:,b1).*mskW(ij,:);
0135 mskG(j,:,b2)=mskG(j,:,b2).*mskW(ij,:);
0136 elseif typ == 2,
0137 mskG(j,:,b1)=mskG(j,:,b1).*mskS(ij,:);
0138 mskG(j,:,b2)=mskG(j,:,b2).*mskS(ij,:);
0139 end
0140 end
0141 end
0142 end
0143 mskG=min(2,mskG);
0144 else
0145 if Tprt, TimeT3=clock; end
0146 end
0147
0148 if Tprt, TimeT4=clock; end
0149
0150 %- to keep psi=0 on top & bottom
0151 mskZ(:,[2:nr+1],:)=mskV;
0152 mskZ(:,[1:nr],:)=mskZ(:,[1:nr],:)+mskV;
0153 %- to keep psi=0 on lateral boundaries :
0154 mskZ([1:ydim],:,:)=mskZ([1:ydim],:,:)+mskZ([2:ydim+1],:,:);
0155 mskZ([2:ydim+1],:,:)=mskZ([2:ydim+1],:,:)+mskZ([3:ydim+2],:,:);
0156 mskZ=ceil(mskZ); mskZ=min(1,mskZ);
0157 if kMsep & nBas > 0,
0158 mskM=zeros(ydim+2,nr,1+nBas); mskM(2:ydim+2,:,:)=min(2-mskG,1);
0159 mskM(1:ydim+1,:,:)=mskM(1:ydim+1,:,:)+mskM(2:ydim+2,:,:);
0160 mskZ(:,1:nr,:)=min(mskZ(:,1:nr,:),mskM);
0161 end
0162 %- apply the mask (and remove dim = 1) :
0163 if Nit == 1,
0164 psi=squeeze(psi); mskV=squeeze(mskV); mskZ=squeeze(mskZ);
0165 psi( find(mskZ==0) )=NaN ;
0166 else
0167 for nt=1:Nit,
0168 psi1=psi(:,:,:,nt); psi1( find(mskZ==0) )=NaN ; psi(:,:,:,nt)=psi1;
0169 end;
0170 if nBas < 1, psi=squeeze(psi); mskV=squeeze(mskV); mskZ=squeeze(mskZ); end
0171 end
0172 else
0173 if Tprt, TimeT3=TimeT2; TimeT4=TimeT2; end
0174 if nBas < 1 | Nit == 1, psi=squeeze(psi); end
0175 end;
0176 %-----------------
0177
0178 if Tprt, TimeT5=clock;
0179 fprintf([' ---- Load, Comp.1,2,3,4 Total time Psi:', ...
0180 ' %6.3f %6.3f %6.3f %6.3f %6.3f %9.6f \n'],...
0181 etime(TimeT1,TimeT0), etime(TimeT2,TimeT1), ...
0182 etime(TimeT3,TimeT2), etime(TimeT4,TimeT3), ...
0183 etime(TimeT5,TimeT4), etime(TimeT5,TimeT0) );
0184 end
0185
0186 return
