Issuing rollback() due to DESTROY without explicit disconnect() of DBD::mysql::db handle dbname=MITgcm at /usr/local/share/lxr/lib/LXR/Common.pm line 1224.
Last-Modified: Fri, 28 Nov 2024 06:11:57 GMT
Content-Type: text/html; charset=utf-8
MITgcm/MITgcm/utils/matlab/cs_grid/calcEulerPsiCube.m
7ab08d6865 Dani*0001 function [psi,mskG,ylat] = calcEulerPsiCube(varargin);
00020003 % [psi,mskG,ylat] = calcEulerPsiCube(d,g,flu,rstar,blkFile,[optional]);
0004 % [psi,mskG,ylat] = calcEulerPsiCube(d,g,flu,rstar,blkFile,'grav',9.81);
0005 %
0006 % Input arguements:
0007 % The incoming field data (d) and grid data (g) must be in a structured
0008 % array format (which is the format that comes from rdmnc):
0009 % d [Field data] hUtave,hVtave OR uVeltave,vVeltave
0010 % g [Grid data ] drF,dxG,dyG,HFacW,HFacS
0011 % Other input parameters:
0012 % flu (str) 'O' or 'A' for ocean or atmosphere
0013 % rstar (str) 0 or 1 if you are using r* coordinates or not
0014 % blkFile (str) Broken line file ('isoLat_cs32_59.mat')
0015 % Optional parameters:
0016 % 'grav' (num) Acceleration due to gravity (default 9.81);
0017 %
0018 % Output fields:
0019 % psi Overturning (eg [61,6,nt])
0020 % mskG Land mask (eg [60,5])
0021 % ylat Latitude coordinate of psi (eg [61,1])
0022 %
0023 % Description:
0024 % Caculates overturning stream function (psi). For the atmosphere, data
0025 % is must be in p-coordinates and the output is the mass transport psi
0026 % [10^9 kg/s]. For the ocean, data should be in z-coordinates and the
0027 % output is the volume transport psi [10^6 m^3/s = Sv]. If the rstar
0028 % parameters is on, hu and hv are used, if off, the hfacw*.u and hfacs*.v
0029 % are used (the multiplication being done inside the function).
0030 %
0031 % 'psi' is tabulated on broken lines at the interface between cells in
0032 % the vertical. 'mskG' is for the area between broken lines and between
0033 % the cell interfaces in the vertical.
0034 %
0035 % Original Author: Jean-Michel Campin
0036 % Modifications: Daniel Enderton
00370038 % Defaults that can be overriden.
0039 grav = 9.81;
00400041 % Read input parameters.
0042 d = varargin{1};
0043 g = varargin{2};
0044 flu = varargin{3};
0045 rstar = varargin{4};
0046 blkFile = varargin{5};
0047 for ii = 6:2:length(varargin)
0048 temp = varargin{ii+1}; eval([varargin{ii},' = temp;']);
0049 end
005000510052 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0053 % Prepare / reform incoming data %
0054 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00550056 nc = size(g.XC,2);
0057 nr = length(g.drF);
00580059 delM = g.drF;
0060 hw = reshape(g.HFacW(1:6*nc,1:nc,1:nr),[6*nc*nc,nr]);
0061 hs = reshape(g.HFacS(1:6*nc,1:nc,1:nr),[6*nc*nc,nr]);
0062 dxg = reshape(g.dxG(1:6*nc,1:nc),[6*nc*nc,1]).*1e-6;
0063 dyg = reshape(g.dyG(1:6*nc,1:nc),[6*nc*nc,1]).*1e-6;
0064 if rstar
0065 nt = size(d.hUtave,4);
0066 hu = reshape(d.hUtave(1:6*nc,1:nc,1:nr,1:nt),[6*nc*nc,nr,nt]);
0067 hv = reshape(d.hVtave(1:6*nc,1:nc,1:nr,1:nt),[6*nc*nc,nr,nt]);
0068 else
0069 nt = size(d.uVeltave,4);
0070 hu = reshape(d.uVeltave(1:6*nc,1:nc,1:nr,1:nt),[6*nc*nc,nr,nt]);
0071 hv = reshape(d.vVeltave(1:6*nc,1:nc,1:nr,1:nt),[6*nc*nc,nr,nt]);
0072 for it = 1:nt
0073 hu(:,:,it) = hw.*hu(:,:,it);
0074 hv(:,:,it) = hs.*hv(:,:,it);
0075 end
0076 end
00770078 % Load broken information.
0079 % I looked at calc_psiH_CS.m, and did not find it very clear.
0080 % May be you can try to see what is in
0081 % MITgcm/utils/matlab/cs_grid/bk_line/use_psiLine.m
0082 % it's shorter, and slightly better.
0083 load(blkFile);
0084 ydim = length(bkl_Ylat);
0085 ylat = [-90,bkl_Ylat,90];
00860087 % For the volume transport through a cell face, the area of the face times
0088 % the velocity through the face is calculated. For the ocean, the area is
0089 % calculated simply by the vertical cell length times the horizontal cell
0090 % length. In the atmosphere, pressure coordinates are used. Using the
0091 % hydrostatic relationship, delM = rho*dz = -dp/g. by leaving the rho on
0092 % the side of the dz and using this for the vertical cell length the volume
0093 % transport is converted to a mass transport as desired. Note the the
0094 % additional factor of 1e-3 is the difference between the 10^6 and 10^9
0095 % used for the ocean and atmosphere, respectively. The 1e-6 is applied at
0096 % the end for the final conversion into the desired units.
0097 if isequal(flu,'A'), delM = -(1e-3).*delM./grav; end
00980099 % kMsep=1;
0100 % if (nargin < 6), kfac=0;
0101 % else kfac=1; end;
0102 nBas=0;
01030104 % Load basin information.
0105 % if nBas > 0,
0106 % % Load ocean basin mask (& separation line):
0107 % mskBw=rdda([rac,'maskW_bas.bin'],[6*nc*nc 3],1,'real*4','b');
0108 % mskBs=rdda([rac,'maskS_bas.bin'],[6*nc*nc 3],1,'real*4','b');
0109 % if nBas==2,
0110 % mskBw(:,2)=mskBw(:,2)+mskBw(:,3);
0111 % mskBs(:,2)=mskBs(:,2)+mskBs(:,3);
0112 % mskBw=min(1,mskBw); mskBs=min(1,mskBs);
0113 % end
0114 % %- load: np_Sep, ij_Sep, tp_Sep:
0115 % sep_lineF=[rac,'sepBas_cs32_60'];
0116 % load(sep_lineF);
0117 % end
01180119 % Prepare arrays.
0120 psi = zeros(ydim+2,nr+1,1+nBas,nt);
0121 mskZ = zeros(ydim+2,nr+1,1+nBas); % Mask of psi
0122 mskV = zeros(ydim+2,nr,1+nBas); % Mask of the merid. transport
0123 mskG = zeros(ydim+1,nr,1+nBas); % Mask of the ground
01240125 % The variable "bkl_Flg" is -1/1 if edge (on a given broken) has a u point
0126 % and -2/2 if it has a v point. Positive/negative values contribute
0127 % positively/negatively to northward heat transport (this depends on the
0128 % oreientation of the cell). A zero value indicates an end of edges that
0129 % contribute to a broken line. The u and v information is parced into two
0130 % seperate fields, ufac and vfac (-2/2 are reduced to -1/1 for vfac).
0131 ufac = zeros([size(bkl_Flg),1+nBas]);
0132 vfac = zeros([size(bkl_Flg),1+nBas]);
0133 ufac(:,:,1) = rem(bkl_Flg,2);
0134 vfac(:,:,1) = fix(bkl_Flg/2);
0135 % for jl=1:ydim,
0136 % ie=bkl_Npts(jl);
0137 % for b=1:nBas,
0138 % ufac(1:ie,jl,1+b)=mskBw(bkl_IJuv(1:ie,jl),b).*ufac(1:ie,jl,1);
0139 % vfac(1:ie,jl,1+b)=mskBs(bkl_IJuv(1:ie,jl),b).*vfac(1:ie,jl,1);
0140 % end;
0141 % end;
014201430144 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0145 % Compute mass/volume stream function %
0146 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01470148 % Compute volume transport through broken lines a hence psi. ut/vt is the
0149 % velocity times the edge length it is passing through. The sum of this
0150 % quantity along a broken line (vz) times the cell height is the volume
0151 % transport through broken line at one layer (delM(k)*vz). psi is then
0152 % the value of the volume transport through the level above subtracted
0153 % from the value of psi above.
0154 for it = 1:nt
0155 for k = nr:-1:1
0156 ut = dyg.*hu(:,k,it);
0157 vt = dxg.*hv(:,k,it);
0158 for jl = 1:ydim
0159 ie = bkl_Npts(jl);
0160 for b = 1:1+nBas
0161 vz = sum( ufac(1:ie,jl,b).*ut(bkl_IJuv(1:ie,jl)) ...
0162 + vfac(1:ie,jl,b).*vt(bkl_IJuv(1:ie,jl)) );
0163 psi(jl+1,k,b,it) = psi(jl+1,k+1,b,it) - delM(k)*vz;
0164 end
0165 end
0166 end
0167 end
01680169 % % Compute the mask :
0170 % if kfac == 1,
0171 % ufac=abs(ufac) ; vfac=abs(vfac);
0172 % for jl=1:ydim,
0173 % ie=bkl_Npts(jl);
0174 % hw=zeros(ie,nr); hs=zeros(ie,nr);
0175 % hw=hw(bkl_IJuv(1:ie,jl),:); % Would need correction!
0176 % hs=hs(bkl_IJuv(1:ie,jl),:);
0177 % for b=1:1+nBas,
0178 % for k=1:nr,
0179 % % for ii=1:bkl_Npts(jl);
0180 % % ij=bkl_IJuv(ii,jl);
0181 % % mskV(jl+1,k,b)=mskV(jl+1,k,b)+ufac(ii,jl,b)*hw(ij,k)+vfac(ii,jl,b)*hs(ij,k);
0182 % % end ;
0183 % tmpv=ufac(1:ie,jl,b).*hw(:,k)+vfac(1:ie,jl,b).*hs(:,k);
0184 % mskV(jl+1,k,b)=mskV(jl+1,k,b)+max(tmpv);
0185 % end
0186 % end
0187 % end
0188 % mskV=ceil(mskV); mskV=min(1,mskV);
0189 % %- build the real mask (=mskG, ground) used to draw the continent with "surf":
0190 % % position=centered , dim= ydim+1 x nr
0191 % mskG=mskV(1:ydim+1,:,:)+mskV(2:ydim+2,:,:); mskG=min(1,mskG);
0192 %
0193 % if kMsep & nBas > 0,
0194 % mskW=1+min(1,ceil(hw));
0195 % mskS=1+min(1,ceil(hs));
0196 % for b=1:nBas,
0197 % bs=b; b1=1+bs; b2=2+rem(bs,nBas);
0198 % if nBas == 2, bs=b+b-1; b1=2; b2=3 ; end
0199 % for j=1:ydim+1,
0200 % for i=1:np_Sep(bs,j),
0201 % ij=ij_Sep(bs,j,i); typ=abs(tp_Sep(bs,j,i));
0202 % if typ == 1,
0203 % mskG(j,:,b1)=mskG(j,:,b1).*mskW(ij,:);
0204 % mskG(j,:,b2)=mskG(j,:,b2).*mskW(ij,:);
0205 % elseif typ == 2,
0206 % mskG(j,:,b1)=mskG(j,:,b1).*mskS(ij,:);
0207 % mskG(j,:,b2)=mskG(j,:,b2).*mskS(ij,:);
0208 % end
0209 % end
0210 % end
0211 % end
0212 % mskG=min(2,mskG);
0213 % end
0214 %
0215 % %- to keep psi=0 on top & bottom
0216 % mskZ(:,[2:nr+1],:)=mskV;
0217 % mskZ(:,[1:nr],:)=mskZ(:,[1:nr],:)+mskV;
0218 % %- to keep psi=0 on lateral boundaries :
0219 % mskZ([1:ydim],:,:)=mskZ([1:ydim],:,:)+mskZ([2:ydim+1],:,:);
0220 % mskZ([2:ydim+1],:,:)=mskZ([2:ydim+1],:,:)+mskZ([3:ydim+2],:,:);
0221 % mskZ=ceil(mskZ); mskZ=min(1,mskZ);
0222 % if kMsep & nBas > 0,
0223 % mskM=zeros(ydim+2,nr,1+nBas); mskM(2:ydim+2,:,:)=min(2-mskG,1);
0224 % mskM(1:ydim+1,:,:)=mskM(1:ydim+1,:,:)+mskM(2:ydim+2,:,:);
0225 % mskZ(:,1:nr,:)=min(mskZ(:,1:nr,:),mskM);
0226 % end
0227 % %- apply the mask (and remove dim = 1) :
0228 % if nt == 1,
0229 % psi=squeeze(psi); mskV=squeeze(mskV); mskZ=squeeze(mskZ);
0230 % psi( find(mskZ==0) )=NaN ;
0231 % else
0232 % for nt=1:nt,
0233 % psi1=psi(:,:,:,nt); psi1( find(mskZ==0) )=NaN ; psi(:,:,:,nt)=psi1;
0234 % end
0235 % if nBas < 1, psi=squeeze(psi); mskV=squeeze(mskV); mskZ=squeeze(mskZ); end
0236 % end
0237 % else
0238 % if nBas < 1 | nt == 1, psi=squeeze(psi); end
0239 % end
02400241 psi = squeeze(psi);
