--- MITgcm_contrib/dfer/matlab_stuff/calcEulerPsiCube.m 2018/03/07 22:01:18 1.2 +++ MITgcm_contrib/dfer/matlab_stuff/calcEulerPsiCube.m 2018/04/19 00:03:59 1.3 @@ -1,41 +1,43 @@ function [psi,mskG,ylat] = calcEulerPsiCube(varargin); -% [psi,mskG,ylat] = calcEulerPsiCube(d,g,flu,rstar,blkFile,[optional]); +% [psi,mskG,ylat] = calcEulerPsiCube(d,g,flu,rstar,blkFile,[mask]); % -% Input arguements: -% d [Field data] Velocity field (Mass-weighted if rstar=1): -% UVELMASS,VVELMASS for rstar=1 -% UVEL,VVEL for rstar=0 -% g [Grid data ] drF,dxG,dyG,HFacW,HFacS -% flu (str) 'O' or 'A' for ocean or atmosphere -% rstar (int) 0 or 1 if you are using r* coordinates or not -% blkFile (str) Broken line file ('isoLat_cs32_59.mat') -% The incoming field data (d) and grid data (g) must be in a -% structured array format -% Optional parameters: -% mask (struct) mask (structured array including maskW and maskS) +% Input arguments: +% d [structure] Velocity field (Mass-weighted if rstar=1): +% UVELMASS,VVELMASS for rstar=1 +% UVEL,VVEL for rstar=0 +% g [structure] drF,dxG,dyG,HFacW,HFacS +% flu [string] 'O' or 'A' for ocean or atmosphere +% rstar [integer] 1 or 0 if you are using r* coordinates or not +% blkFile [string] Broken line file ('isoLat_cs32_59.mat') +% +% Optional inputs: +% mask [structure] Optional: Mask field for computation per basin, it assumes that +% maskW and maskS are provided in a structure % % Output fields: -% psi Overturning (eg [61,6,nt]) -% mskG Land mask (eg [60,5]) -% ylat Latitude coordinate of psi (eg [61,1]) +% psi Overturning +% mskG Land mask +% ylat Latitude coordinate of psi % % Description: -% Caculates overturning stream function (psi). For the atmosphere, data +% Calculates overturning streamfunction (psi). For the atmosphere, data % is must be in p-coordinates and the output is the mass transport psi -% [10^9 kg/s]. For the ocean, data should be in z-coordinates and the -% output is the volume transport psi [10^6 m^3/s = Sv]. If the rstar -% parameters is on, hu and hv are used, if off, the hfacw*.u and hfacs*.v -% are used (the multiplication being done inside the function). +% [10^9 kg/s]. For the ocean, data should be in z-coordinates and the +% output is the volume transport psi [10^6 m^3/s = Sv]. If rstar +% is on (=1), UVELMASS and VVELMASS must be used. If off (rstar=0), +% the hfacw*.UVEL and hfacs*.VVEL are used (the multiplication being +% done inside the function). % % 'psi' is tabulated on broken lines at the interface between cells in -% the vertical. 'mskG' is for the area between broken lines and between +% the vertical. 'mskG' is for the area between broken lines and between % the cell interfaces in the vertical. % % Defaults that can be overriden. grav = 9.81; masking=0; +nBas=0; % Read input parameters. d = varargin{1}; @@ -48,7 +50,6 @@ masking = 1; end - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Prepare / reform incoming data % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -86,19 +87,10 @@ end % Load broken information. -% I looked at calc_psiH_CS.m, and did not find it very clear. -% May be you can try to see what is in -% MITgcm/utils/matlab/cs_grid/bk_line/use_psiLine.m -% it s shorter, and slightly better. load(blkFile); ydim = length(bkl_Ylat); ylat = [-90,bkl_Ylat,90]; -% kMsep=1; -% if (nargin < 6), kfac=0; -% else kfac=1; end; -nBas=0; - % Prepare arrays. psi = zeros(ydim+2,nr+1,1+nBas,nt); mskZ = zeros(ydim+2,nr+1,1+nBas); % Mask of psi @@ -115,14 +107,6 @@ 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 mass/volume stream function % @@ -156,76 +140,3 @@ if isequal(flu,'O'), psi = 1e-6*squeeze(psi); end if isequal(flu,'A'), psi =-1e-9/grav*squeeze(psi); end - - -% % Compute the mask : -% if kfac == 1, -% ufac=abs(ufac) ; vfac=abs(vfac); -% for jl=1:ydim, -% ie=bkl_Npts(jl); -% hw=zeros(ie,nr); hs=zeros(ie,nr); -% hw=hw(bkl_IJuv(1:ie,jl),:); % Would need correction! -% hs=hs(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)*hw(ij,k)+vfac(ii,jl,b)*hs(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 kMsep & nBas > 0, -% mskW=1+min(1,ceil(hw)); -% mskS=1+min(1,ceil(hs)); -% 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); -% 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 nt == 1, -% psi=squeeze(psi); mskV=squeeze(mskV); mskZ=squeeze(mskZ); -% psi( find(mskZ==0) )=NaN ; -% else -% for nt=1:nt, -% 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 nBas < 1 | nt == 1, psi=squeeze(psi); end -% end