Warning, /utils/matlab/cs_grid/calcHorizPsiCube.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit 4ec37fd8 on 2023-10-03 22:00:32 UTC
c97281acbe Dani*0001 function [psiH,xv,yv,psiH_cubeC] = calcHorizPsiCube(d,g,rstar,psiLineF);
ad5563f4b2 Dani*0002
1d9f99be15 Dani*0003 % [psiH,xv,yv,psiH_cubeC] = calcHorizPsiCube(d,g,rstar,psiLineF);
0004 %
4ec37fd829 Jean*0005 % Input arguements:
1d9f99be15 Dani*0006 % The incoming field data (d) and grid data (g) must be in a structured
0007 % array format (which is the format that comes from rdmnc):
0008 % d [Field data] hUtave,hVtave (rstar=1) OR uVeltave,vVeltave (=0)
0009 % g [Grid data ] drF,dxG,dyG,HFacW,HFacS
0010 % Other input parameters:
0011 % rstar (str) 0 or 1 if you are using r* coordinates or not
0012 % psiLineF (str) File ('psiLine_N2S_cs32')
0013 %
0014 % Output fields:
0015 % psiH Barotropic streamfunction, vector points (eg [6146,nt])
0016 % xv Longitude of vector points (eg [6146,nt])
0017 % yv Latitude of vector points (eg [6146,nt])
0018 % psiH_cubeC 'psiH' interpolated to cell center points (eg [192,32,nt])
0019 %
0020 % Description:
0021 % Caculates barotropic stream function (psi). Data should be in
0022 % z-coordinates and the output is the volume transport psi [10^6 m^3/s =
0023 % Sv]. If the rstar parameters is on, hu and hv are used, if off, the
0024 % hfacw*.u and hfacs*.v are used (the multiplication being done inside
0025 % the function).
0026 %
0027 % Original Author: Jean-Michel Campin
0028 % Modifications: Daniel Enderton
ad5563f4b2 Dani*0029
0030 nc = size(g.XC,2);
0031 nr = length(g.drF);
0032 nt = size(d.uVeltave,4);
0033
0034 xv = reshape(g.XG(1:6*nc,1:nc),[6*nc*nc,1]);
0035 yv = reshape(g.YG(1:6*nc,1:nc),[6*nc*nc,1]);
0036 xv(end+1)=xv(1); yv(end+1)=yv(1+2*nc);
0037 xv(end+1)=xv(1+3*nc); yv(end+1)=yv(1);
0038
0039 dxg = reshape(g.dxG(1:6*nc,1:nc),[6*nc*nc,1])*1.e-6;
0040 dyg = reshape(g.dyG(1:6*nc,1:nc),[6*nc*nc,1])*1.e-6;
0041 delR = reshape(g.drF,[1,nr]);
0042 hw = reshape(g.HFacW(1:6*nc,1:nc,1:nr),[6*nc*nc,nr]);
0043 hs = reshape(g.HFacS(1:6*nc,1:nc,1:nr),[6*nc*nc,nr]);
0044 if rstar
0045 hu = reshape(d.hUtave(1:6*nc,1:nc,1:nr,1:nt),[6*nc*nc,nr,nt]);
0046 hv = reshape(d.hVtave(1:6*nc,1:nc,1:nr,1:nt),[6*nc*nc,nr,nt]);
0047 else
0048 hu = reshape(d.uVeltave(1:6*nc,1:nc,1:nr,1:nt),[6*nc*nc,nr,nt]);
0049 hv = reshape(d.vVeltave(1:6*nc,1:nc,1:nr,1:nt),[6*nc*nc,nr,nt]);
0050 for it = 1:nt
0051 hu(:,:,it) = hw.*hu(:,:,it);
0052 hv(:,:,it) = hs.*hv(:,:,it);
0053 end
0054 end
0055
1d9f99be15 Dani*0056 % Load: 'psi_N','psi_C','psi_P','psiUV','psi_T'
ad5563f4b2 Dani*0057 load(psiLineF);
0058
0059
0060 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0061 % Calculate horizontal stream function %
0062 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0063
0064 % Compute depth integrated volume transport through faces.
c97281acbe Dani*0065 utrs = NaN.*zeros([6*nc*nc,nt]);
0066 vtrs = NaN.*zeros([6*nc*nc,nt]);
0067 for it = 1:nt
0068 utrs(:,it) = sum((dyg*delR).*hu(:,:,it),2);
0069 vtrs(:,it) = sum((dxg*delR).*hv(:,:,it),2);
0070 end
4ec37fd829 Jean*0071 %-- To use psiLine indices, put long-vector utrs,vtrs in "compact-format":
0072 utrs=reshape(permute(reshape(utrs,[nc 6 nc nt]),[1 3 2 4]),[nc*nc*6 nt]);
0073 vtrs=reshape(permute(reshape(vtrs,[nc 6 nc nt]),[1 3 2 4]),[nc*nc*6 nt]);
ad5563f4b2 Dani*0074
0075 % Compute barotropic stream function. A little description of what the
0076 % variables are and how the computation is done would be nice.
c97281acbe Dani*0077 psiH = zeros(6*nc*nc+2,nt);
0078 for it = 1:nt
0079 psiNx = size(psi_C,1);
0080 psiNy = size(psi_C,2);
0081 %nPt2 = ;
0082 ufac = rem(psi_T,2); % Mask for u transport.
0083 vfac = fix(psi_T/2); % Mask for v transport.
0084 for j = 2:psiNy
0085 for i = 1:psi_N(j)
0086 i0 = psi_P(i,j);
0087 i1 = psi_C(i,j);
0088 i2 = psiUV(i,j);
c599e640f4 Jean*0089 psiH(i1,it) = psiH(i0,it) ...
c97281acbe Dani*0090 + ufac(i,j)*utrs(i2,it) ...
0091 + vfac(i,j)*vtrs(i2,it);
0092 end
ad5563f4b2 Dani*0093 end
0094 end
0095
0096
0097 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0098 % Interpolate to grid tracer (cell center) points. %
0099 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0100
c97281acbe Dani*0101 psiH_cubeC = NaN.*zeros([6*nc,nc,nt]);
0102 for it = 1:nt
0103 psiH_extra = psiH(end-1:end,it);
4ec37fd829 Jean*0104 %-- below is the original version, now commented out:
0105 % temp = reshape(psiH(1:end-2,it),[6*nc,nc]);
c97281acbe Dani*0106
4ec37fd829 Jean*0107 % [nx ny nt]=size(temp); %- resetting "nt" to 1 inside loop it=1:nt is not great
0108 % temp=permute( reshape(temp,[nx/6 6 ny]),[1 3 2]);
0109 %-- now "long-vector" psiH is in "compact-format", no need to permute dims
0110 temp = reshape(psiH(1:end-2,it),[nc nc 6]);
0111 %--- update ends here.
c97281acbe Dani*0112 temp(end+1,:,:)=NaN;
0113 temp(:,end+1,:)=NaN;
0114 temp(end,:,[1 3 5])=temp(1,:,[2 4 6]);
0115 temp(:,end,[2 4 6])=temp(:,1,[3 5 1]);
0116 temp(:,end,[1 3 5])=squeeze(temp(1,end:-1:1,[3 5 1]));
0117 temp(end,:,[2 4 6])=squeeze(temp(end:-1:1,1,[4 6 2]));
0118 temp(1,end,[1 3 5]) = psiH_extra(1);
0119 temp(end,1,[2 4 6]) = psiH_extra(2);
0120
0121 temp = 0.25 .* (temp(1:nc ,1:nc ,:) + ...
0122 temp(2:nc+1,1:nc ,:) + ...
0123 temp(1:nc ,2:nc+1,:) + ...
0124 temp(2:nc+1,2:nc+1,:));
0125
0126 psiH_cubeC(:,:,it) = reshape(permute(temp,[1 3 2]),[6*nc,nc]);
c599e640f4 Jean*0127 end
