| 1 |
function [psi,mskG,ylat] = calcEulerPsiCube(varargin); |
| 2 |
|
| 3 |
% [psi,mskG,ylat] = calcEulerPsiCube(d,g,flu,rstar,blkFile,[optional]); |
| 4 |
% |
| 5 |
% Input arguements: |
| 6 |
% d [Field data] Velocity field (Mass-weighted if rstar=1): |
| 7 |
% UVELMASS,VVELMASS for rstar=1 |
| 8 |
% UVEL,VVEL for rstar=0 |
| 9 |
% g [Grid data ] drF,dxG,dyG,HFacW,HFacS |
| 10 |
% flu (str) 'O' or 'A' for ocean or atmosphere |
| 11 |
% rstar (int) 0 or 1 if you are using r* coordinates or not |
| 12 |
% blkFile (str) 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 |
| 15 |
% Optional parameters: |
| 16 |
% mask (struct) mask (structured array including maskW and maskS) |
| 17 |
% |
| 18 |
% Output fields: |
| 19 |
% psi Overturning (eg [61,6,nt]) |
| 20 |
% mskG Land mask (eg [60,5]) |
| 21 |
% ylat Latitude coordinate of psi (eg [61,1]) |
| 22 |
% |
| 23 |
% Description: |
| 24 |
% Caculates overturning stream function (psi). For the atmosphere, data |
| 25 |
% 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 |
| 27 |
% output is the volume transport psi [10^6 m^3/s = Sv]. If the rstar |
| 28 |
% parameters is on, hu and hv are used, if off, the hfacw*.u and hfacs*.v |
| 29 |
% are used (the multiplication being done inside the function). |
| 30 |
% |
| 31 |
% 'psi' is tabulated on broken lines at the interface between cells in |
| 32 |
% the vertical. 'mskG' is for the area between broken lines and between |
| 33 |
% the cell interfaces in the vertical. |
| 34 |
% |
| 35 |
|
| 36 |
% Defaults that can be overriden. |
| 37 |
grav = 9.81; |
| 38 |
masking=0; |
| 39 |
|
| 40 |
% Read input parameters. |
| 41 |
d = varargin{1}; |
| 42 |
g = varargin{2}; |
| 43 |
flu = varargin{3}; |
| 44 |
rstar = varargin{4}; |
| 45 |
blkFile = varargin{5}; |
| 46 |
if length(varargin) == 6 |
| 47 |
mask = varargin{6}; |
| 48 |
masking = 1; |
| 49 |
end |
| 50 |
|
| 51 |
|
| 52 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 53 |
% Prepare / reform incoming data % |
| 54 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 55 |
|
| 56 |
nc = size(g.XC,2); |
| 57 |
nr = length(g.drF); |
| 58 |
|
| 59 |
delM = g.drF; |
| 60 |
dxg = reshape(g.dxG(1:6*nc,1:nc),[6*nc*nc,1]); |
| 61 |
dyg = reshape(g.dyG(1:6*nc,1:nc),[6*nc*nc,1]); |
| 62 |
if rstar |
| 63 |
nt = size(d.UVELMASS,4); |
| 64 |
hu = reshape(d.UVELMASS(1:6*nc,1:nc,1:nr,1:nt),[6*nc*nc,nr,nt]); |
| 65 |
hv = reshape(d.VVELMASS(1:6*nc,1:nc,1:nr,1:nt),[6*nc*nc,nr,nt]); |
| 66 |
else |
| 67 |
nt = size(d.UVEL,4); |
| 68 |
hw = reshape(g.HFacW(1:6*nc,1:nc,1:nr),[6*nc*nc,nr]); |
| 69 |
hs = reshape(g.HFacS(1:6*nc,1:nc,1:nr),[6*nc*nc,nr]); |
| 70 |
hu = reshape(d.UVEL(1:6*nc,1:nc,1:nr,1:nt),[6*nc*nc,nr,nt]); |
| 71 |
hv = reshape(d.VVEL(1:6*nc,1:nc,1:nr,1:nt),[6*nc*nc,nr,nt]); |
| 72 |
for it = 1:nt |
| 73 |
hu(:,:,it) = hw.*hu(:,:,it); |
| 74 |
hv(:,:,it) = hs.*hv(:,:,it); |
| 75 |
end |
| 76 |
end |
| 77 |
|
| 78 |
mskWloc = ones(6*nc*nc,1); |
| 79 |
mskSloc = ones(6*nc*nc,1); |
| 80 |
|
| 81 |
if masking == 1 |
| 82 |
mskWloc=reshape(mask.maskW(:,:,1),6*nc*nc,1); |
| 83 |
mskSloc=reshape(mask.maskS(:,:,1),6*nc*nc,1); |
| 84 |
%hu = repmat(reshape(mask.maskW,6*nc*nc,1),[1 nr nt]) .* hu; |
| 85 |
%hv = repmat(reshape(mask.maskS,6*nc*nc,1),[1 nr nt]) .* hv; |
| 86 |
end |
| 87 |
|
| 88 |
% Load broken information. |
| 89 |
% I looked at calc_psiH_CS.m, and did not find it very clear. |
| 90 |
% May be you can try to see what is in |
| 91 |
% MITgcm/utils/matlab/cs_grid/bk_line/use_psiLine.m |
| 92 |
% it s shorter, and slightly better. |
| 93 |
load(blkFile); |
| 94 |
ydim = length(bkl_Ylat); |
| 95 |
ylat = [-90,bkl_Ylat,90]; |
| 96 |
|
| 97 |
% kMsep=1; |
| 98 |
% if (nargin < 6), kfac=0; |
| 99 |
% else kfac=1; end; |
| 100 |
nBas=0; |
| 101 |
|
| 102 |
% Prepare arrays. |
| 103 |
psi = zeros(ydim+2,nr+1,1+nBas,nt); |
| 104 |
mskZ = zeros(ydim+2,nr+1,1+nBas); % Mask of psi |
| 105 |
mskV = zeros(ydim+2,nr,1+nBas); % Mask of the merid. transport |
| 106 |
mskG = zeros(ydim+1,nr,1+nBas); % Mask of the ground |
| 107 |
|
| 108 |
% The variable "bkl_Flg" is -1/1 if edge (on a given broken) has a u point |
| 109 |
% and -2/2 if it has a v point. Positive/negative values contribute |
| 110 |
% positively/negatively to northward heat transport (this depends on the |
| 111 |
% oreientation of the cell). A zero value indicates an end of edges that |
| 112 |
% contribute to a broken line. The u and v information is parced into two |
| 113 |
% seperate fields, ufac and vfac (-2/2 are reduced to -1/1 for vfac). |
| 114 |
ufac = zeros([size(bkl_Flg),1+nBas]); |
| 115 |
vfac = zeros([size(bkl_Flg),1+nBas]); |
| 116 |
ufac(:,:,1) = rem(bkl_Flg,2); |
| 117 |
vfac(:,:,1) = fix(bkl_Flg/2); |
| 118 |
% for jl=1:ydim, |
| 119 |
% ie=bkl_Npts(jl); |
| 120 |
% for b=1:nBas, |
| 121 |
% ufac(1:ie,jl,1+b)=mskBw(bkl_IJuv(1:ie,jl),b).*ufac(1:ie,jl,1); |
| 122 |
% vfac(1:ie,jl,1+b)=mskBs(bkl_IJuv(1:ie,jl),b).*vfac(1:ie,jl,1); |
| 123 |
% end; |
| 124 |
% end; |
| 125 |
|
| 126 |
|
| 127 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 128 |
% Compute mass/volume stream function % |
| 129 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 130 |
|
| 131 |
% Compute volume transport through broken lines a hence psi. ut/vt is the |
| 132 |
% velocity times the edge length it is passing through. The sum of this |
| 133 |
% quantity along a broken line (vz) times the cell height is the volume |
| 134 |
% transport through broken line at one layer (delM(k)*vz). psi is then |
| 135 |
% the value of the volume transport through the level above subtracted |
| 136 |
% from the value of psi above. |
| 137 |
for it = 1:nt |
| 138 |
for k = nr:-1:1 |
| 139 |
ut = dyg.*hu(:,k,it).*mskWloc; |
| 140 |
vt = dxg.*hv(:,k,it).*mskSloc; |
| 141 |
for jl = 1:ydim |
| 142 |
ie = bkl_Npts(jl); |
| 143 |
for b = 1:1+nBas |
| 144 |
vz = sum( ufac(1:ie,jl,b).*ut(bkl_IJuv(1:ie,jl)) ... |
| 145 |
+ vfac(1:ie,jl,b).*vt(bkl_IJuv(1:ie,jl)) ); |
| 146 |
psi(jl+1,k,b,it) = psi(jl+1,k+1,b,it) - delM(k)*vz; |
| 147 |
end |
| 148 |
end |
| 149 |
end |
| 150 |
end |
| 151 |
|
| 152 |
psi = squeeze(psi); |
| 153 |
|
| 154 |
%% For Ocean, result in Sv (10^6 m3/s) |
| 155 |
%% For Atmos, results in 10^9 kg/s |
| 156 |
if isequal(flu,'O'), psi = 1e-6*squeeze(psi); end |
| 157 |
if isequal(flu,'A'), psi =-1e-9/grav*squeeze(psi); end |
| 158 |
|
| 159 |
|
| 160 |
|
| 161 |
% % Compute the mask : |
| 162 |
% if kfac == 1, |
| 163 |
% ufac=abs(ufac) ; vfac=abs(vfac); |
| 164 |
% for jl=1:ydim, |
| 165 |
% ie=bkl_Npts(jl); |
| 166 |
% hw=zeros(ie,nr); hs=zeros(ie,nr); |
| 167 |
% hw=hw(bkl_IJuv(1:ie,jl),:); % Would need correction! |
| 168 |
% hs=hs(bkl_IJuv(1:ie,jl),:); |
| 169 |
% for b=1:1+nBas, |
| 170 |
% for k=1:nr, |
| 171 |
% % for ii=1:bkl_Npts(jl); |
| 172 |
% % ij=bkl_IJuv(ii,jl); |
| 173 |
% % mskV(jl+1,k,b)=mskV(jl+1,k,b)+ufac(ii,jl,b)*hw(ij,k)+vfac(ii,jl,b)*hs(ij,k); |
| 174 |
% % end ; |
| 175 |
% tmpv=ufac(1:ie,jl,b).*hw(:,k)+vfac(1:ie,jl,b).*hs(:,k); |
| 176 |
% mskV(jl+1,k,b)=mskV(jl+1,k,b)+max(tmpv); |
| 177 |
% end |
| 178 |
% end |
| 179 |
% end |
| 180 |
% mskV=ceil(mskV); mskV=min(1,mskV); |
| 181 |
% %- build the real mask (=mskG, ground) used to draw the continent with "surf": |
| 182 |
% % position=centered , dim= ydim+1 x nr |
| 183 |
% mskG=mskV(1:ydim+1,:,:)+mskV(2:ydim+2,:,:); mskG=min(1,mskG); |
| 184 |
% |
| 185 |
% if kMsep & nBas > 0, |
| 186 |
% mskW=1+min(1,ceil(hw)); |
| 187 |
% mskS=1+min(1,ceil(hs)); |
| 188 |
% for b=1:nBas, |
| 189 |
% bs=b; b1=1+bs; b2=2+rem(bs,nBas); |
| 190 |
% if nBas == 2, bs=b+b-1; b1=2; b2=3 ; end |
| 191 |
% for j=1:ydim+1, |
| 192 |
% for i=1:np_Sep(bs,j), |
| 193 |
% ij=ij_Sep(bs,j,i); typ=abs(tp_Sep(bs,j,i)); |
| 194 |
% if typ == 1, |
| 195 |
% mskG(j,:,b1)=mskG(j,:,b1).*mskW(ij,:); |
| 196 |
% mskG(j,:,b2)=mskG(j,:,b2).*mskW(ij,:); |
| 197 |
% elseif typ == 2, |
| 198 |
% mskG(j,:,b1)=mskG(j,:,b1).*mskS(ij,:); |
| 199 |
% mskG(j,:,b2)=mskG(j,:,b2).*mskS(ij,:); |
| 200 |
% end |
| 201 |
% end |
| 202 |
% end |
| 203 |
% end |
| 204 |
% mskG=min(2,mskG); |
| 205 |
% end |
| 206 |
% |
| 207 |
% %- to keep psi=0 on top & bottom |
| 208 |
% mskZ(:,[2:nr+1],:)=mskV; |
| 209 |
% mskZ(:,[1:nr],:)=mskZ(:,[1:nr],:)+mskV; |
| 210 |
% %- to keep psi=0 on lateral boundaries : |
| 211 |
% mskZ([1:ydim],:,:)=mskZ([1:ydim],:,:)+mskZ([2:ydim+1],:,:); |
| 212 |
% mskZ([2:ydim+1],:,:)=mskZ([2:ydim+1],:,:)+mskZ([3:ydim+2],:,:); |
| 213 |
% mskZ=ceil(mskZ); mskZ=min(1,mskZ); |
| 214 |
% if kMsep & nBas > 0, |
| 215 |
% mskM=zeros(ydim+2,nr,1+nBas); mskM(2:ydim+2,:,:)=min(2-mskG,1); |
| 216 |
% mskM(1:ydim+1,:,:)=mskM(1:ydim+1,:,:)+mskM(2:ydim+2,:,:); |
| 217 |
% mskZ(:,1:nr,:)=min(mskZ(:,1:nr,:),mskM); |
| 218 |
% end |
| 219 |
% %- apply the mask (and remove dim = 1) : |
| 220 |
% if nt == 1, |
| 221 |
% psi=squeeze(psi); mskV=squeeze(mskV); mskZ=squeeze(mskZ); |
| 222 |
% psi( find(mskZ==0) )=NaN ; |
| 223 |
% else |
| 224 |
% for nt=1:nt, |
| 225 |
% psi1=psi(:,:,:,nt); psi1( find(mskZ==0) )=NaN ; psi(:,:,:,nt)=psi1; |
| 226 |
% end |
| 227 |
% if nBas < 1, psi=squeeze(psi); mskV=squeeze(mskV); mskZ=squeeze(mskZ); end |
| 228 |
% end |
| 229 |
% else |
| 230 |
% if nBas < 1 | nt == 1, psi=squeeze(psi); end |
| 231 |
% end |
Parent Directory
| 
Revision Log
| 
Revision Graph