Parent Directory
| 
Revision Log
| 
Revision Graph
| 
Patch
--- 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
| ViewVC Help | |
| Powered by ViewVC 1.1.22 |