/[MITgcm]/MITgcm_contrib/dfer/matlab_stuff/calcEulerPsiCube.m
ViewVC logotype

Diff of /MITgcm_contrib/dfer/matlab_stuff/calcEulerPsiCube.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.2 by dfer, Wed Mar 7 22:01:18 2018 UTC revision 1.3 by dfer, Thu Apr 19 00:03:59 2018 UTC
# Line 1  Line 1 
1  function [psi,mskG,ylat] = calcEulerPsiCube(varargin);  function [psi,mskG,ylat] = calcEulerPsiCube(varargin);
2    
3  % [psi,mskG,ylat] = calcEulerPsiCube(d,g,flu,rstar,blkFile,[optional]);  % [psi,mskG,ylat] = calcEulerPsiCube(d,g,flu,rstar,blkFile,[mask]);
4  %  %
5  % Input arguements:    % Input arguments:
6  %       d  [Field data]  Velocity field (Mass-weighted if rstar=1):  %       d  [structure]    Velocity field (Mass-weighted if rstar=1):
7  %                        UVELMASS,VVELMASS for rstar=1  %                         UVELMASS,VVELMASS for rstar=1
8  %                        UVEL,VVEL         for rstar=0  %                         UVEL,VVEL         for rstar=0
9  %       g  [Grid data ]  drF,dxG,dyG,HFacW,HFacS  %       g  [structure]    drF,dxG,dyG,HFacW,HFacS
10  %       flu      (str)   'O' or 'A' for ocean or atmosphere  %       flu [string]      'O' or 'A' for ocean or atmosphere
11  %       rstar    (int)   0 or 1 if you are using r* coordinates or not  %       rstar [integer]   1 or 0 if you are using r* coordinates or not
12  %       blkFile  (str)   Broken line file ('isoLat_cs32_59.mat')  %       blkFile [string]  Broken line file ('isoLat_cs32_59.mat')
13  %       The incoming field data (d) and grid data (g) must be in a  %
14  %       structured array format  % Optional inputs:
15  %     Optional parameters:  %       mask [structure]  Optional: Mask field for computation per basin, it assumes that
16  %       mask   (struct)  mask (structured array including maskW and maskS)  %                         maskW and maskS are provided in a structure
17  %  %
18  % Output fields:  % Output fields:
19  %       psi     Overturning (eg [61,6,nt])  %       psi     Overturning
20  %       mskG    Land mask (eg [60,5])  %       mskG    Land mask
21  %       ylat    Latitude coordinate of psi (eg [61,1])  %       ylat    Latitude coordinate of psi
22  %  %
23  % Description:  % Description:
24  %   Caculates overturning stream function (psi).  For the atmosphere, data  %   Calculates overturning streamfunction (psi). For the atmosphere, data
25  %   is must be in p-coordinates and the output is the mass transport psi  %   is must be in p-coordinates and the output is the mass transport psi
26  %   [10^9 kg/s].  For the ocean, data should be in z-coordinates and the  %   [10^9 kg/s]. For the ocean, data should be in z-coordinates and the
27  %   output is the volume transport psi [10^6 m^3/s = Sv].  If the rstar  %   output is the volume transport psi [10^6 m^3/s = Sv]. If rstar
28  %   parameters is on, hu and hv are used, if off, the hfacw*.u and hfacs*.v  %   is on (=1), UVELMASS and VVELMASS must be used. If off (rstar=0),
29  %   are used (the multiplication being done inside the function).  %   the hfacw*.UVEL and hfacs*.VVEL are used (the multiplication being
30    %   done inside the function).
31  %  %
32  %   'psi' is tabulated on broken lines at the interface between cells in  %   'psi' is tabulated on broken lines at the interface between cells in
33  %   the vertical.  'mskG' is for the area between broken lines and between  %   the vertical. 'mskG' is for the area between broken lines and between
34  %   the cell interfaces in the vertical.  %   the cell interfaces in the vertical.
35  %  %
36    
37  % Defaults that can be overriden.  % Defaults that can be overriden.
38  grav = 9.81;  grav = 9.81;
39  masking=0;  masking=0;
40    nBas=0;
41    
42  % Read input parameters.  % Read input parameters.
43  d = varargin{1};  d = varargin{1};
# Line 48  if length(varargin) == 6 Line 50  if length(varargin) == 6
50     masking = 1;     masking = 1;
51  end  end
52    
   
53  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
54  %                     Prepare / reform  incoming data                     %  %                     Prepare / reform  incoming data                     %
55  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Line 86  if masking == 1 Line 87  if masking == 1
87  end  end
88    
89  % Load broken information.  % 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.  
90  load(blkFile);  load(blkFile);
91  ydim = length(bkl_Ylat);  ydim = length(bkl_Ylat);
92  ylat = [-90,bkl_Ylat,90];  ylat = [-90,bkl_Ylat,90];
93    
 % kMsep=1;  
 % if (nargin < 6), kfac=0;  
 % else kfac=1; end;  
 nBas=0;  
   
94  % Prepare arrays.  % Prepare arrays.
95  psi = zeros(ydim+2,nr+1,1+nBas,nt);  psi = zeros(ydim+2,nr+1,1+nBas,nt);
96  mskZ = zeros(ydim+2,nr+1,1+nBas);  % Mask of psi  mskZ = zeros(ydim+2,nr+1,1+nBas);  % Mask of psi
# Line 115  ufac = zeros([size(bkl_Flg),1+nBas]); Line 107  ufac = zeros([size(bkl_Flg),1+nBas]);
107  vfac = zeros([size(bkl_Flg),1+nBas]);  vfac = zeros([size(bkl_Flg),1+nBas]);
108  ufac(:,:,1) = rem(bkl_Flg,2);  ufac(:,:,1) = rem(bkl_Flg,2);
109  vfac(:,:,1) = fix(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;  
   
110    
111  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
112  %                   Compute mass/volume stream function                   %  %                   Compute mass/volume stream function                   %
# Line 156  psi = squeeze(psi); Line 140  psi = squeeze(psi);
140  if isequal(flu,'O'), psi = 1e-6*squeeze(psi); end  if isequal(flu,'O'), psi = 1e-6*squeeze(psi); end
141  if isequal(flu,'A'), psi =-1e-9/grav*squeeze(psi); end  if isequal(flu,'A'), psi =-1e-9/grav*squeeze(psi); end
142    
   
   
 % % 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  
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.22