| 1 |
function [PsiB,ylat]=calcBolusPsiCube(varargin); |
| 2 |
|
| 3 |
% [PsiB,ylat]=calcBolusPsiCube(d,g,GMform,blkFile); |
| 4 |
% |
| 5 |
% Compute eddy-induced streamfunction from Gent and McWilliams scheme |
| 6 |
% |
| 7 |
% Input arguments: |
| 8 |
% The incoming field data (d) and grid data (g) must be in a structured |
| 9 |
% array format (which is the format that comes from rdmnc): |
| 10 |
% d [structure] Kwx, Kwy (Skew flux form) or GM_PsiX, GM_PsiY (advective form) |
| 11 |
% g [structure] drF,rA,dxC,dyC,dxG,dyG,HFacW,HFacS |
| 12 |
% GMform [string] GM form: 'Skew' or 'Advc' |
| 13 |
% blkFile [string] Broken line file |
| 14 |
% mask [structure] Optional: Mask field for computation per basin, it assumes that |
| 15 |
% maskW and maskS are provided in a structure |
| 16 |
% Output arguments: |
| 17 |
% PsiB : bolus streamfunction at interface level (in Sv) |
| 18 |
% ylat : meridional coordinate of PsiB |
| 19 |
% |
| 20 |
% Comments: |
| 21 |
% -For Skew-flux form: |
| 22 |
% PsiB computed from Kwx & Kwy divided by 2. |
| 23 |
% first average Kwx and Kwy at u- and v-points: |
| 24 |
% psiX=(rAc*Kwx)_i / (dXc*dYg) ; psiY=(rAc*Kwy)_j / dYc ; |
| 25 |
% and then "zonally" average along broken lines |
| 26 |
% -For Advective form: |
| 27 |
% PsiB computed from PsiX and PsiY |
| 28 |
% just need to "zonally" average along broken lines |
| 29 |
% |
| 30 |
%--------------------------------------------------------------------- |
| 31 |
|
| 32 |
masking=0; |
| 33 |
|
| 34 |
% Read input parameters. |
| 35 |
d = varargin{1}; |
| 36 |
g = varargin{2}; |
| 37 |
GMform = varargin{3}; |
| 38 |
blkFile = varargin{4}; |
| 39 |
if length(varargin) == 5 |
| 40 |
mask = varargin{5}; |
| 41 |
masking = 1; |
| 42 |
end |
| 43 |
|
| 44 |
|
| 45 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 46 |
% Prepare grid stuff % |
| 47 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 48 |
|
| 49 |
nc = size(g.XC,2); |
| 50 |
nr = length(g.drF); |
| 51 |
|
| 52 |
switch GMform |
| 53 |
case 'Skew' |
| 54 |
nt = size(d.GM_Kwx,4); |
| 55 |
case 'Advc' |
| 56 |
nt = size(d.GM_PsiX,4); |
| 57 |
end |
| 58 |
|
| 59 |
%--- areas : |
| 60 |
ra = g.rA; |
| 61 |
dxc = reshape(g.dxC(1:6*nc,1:nc),[6*nc*nc,1]); |
| 62 |
dyc = reshape(g.dyC(1:6*nc,1:nc),[6*nc*nc,1]); |
| 63 |
dxg = reshape(g.dxG(1:6*nc,1:nc),[6*nc*nc,1]); |
| 64 |
dyg = reshape(g.dyG(1:6*nc,1:nc),[6*nc*nc,1]); |
| 65 |
|
| 66 |
rAu=dxc.*dyg; |
| 67 |
rAv=dyc.*dxg; |
| 68 |
|
| 69 |
%--- masks : |
| 70 |
hw = reshape(g.HFacW(1:6*nc,1:nc,1:nr),[6*nc*nc,nr]); |
| 71 |
hs = reshape(g.HFacS(1:6*nc,1:nc,1:nr),[6*nc*nc,nr]); |
| 72 |
mskw=ceil(hw); mskw=min(1,mskw); |
| 73 |
msks=ceil(hs); msks=min(1,msks); |
| 74 |
|
| 75 |
mskWloc = ones(6*nc*nc,1); |
| 76 |
mskSloc = ones(6*nc*nc,1); |
| 77 |
if masking == 1 |
| 78 |
mskWloc=reshape(mask.maskW(:,:,1),6*nc*nc,1); |
| 79 |
mskSloc=reshape(mask.maskS(:,:,1),6*nc*nc,1); |
| 80 |
end |
| 81 |
|
| 82 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 83 |
% Read/Prepare GM fields % |
| 84 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 85 |
psiX_all = zeros(6*nc*nc,nr,nt); |
| 86 |
psiY_all = zeros(6*nc*nc,nr,nt); |
| 87 |
|
| 88 |
switch GMform |
| 89 |
case 'Skew' |
| 90 |
|
| 91 |
kwx_all = 0.5*d.GM_Kwx; |
| 92 |
kwy_all = 0.5*d.GM_Kwy; |
| 93 |
|
| 94 |
for it = 1:nt |
| 95 |
kwx = kwx_all(:,:,:,it); |
| 96 |
kwy = kwy_all(:,:,:,it); |
| 97 |
|
| 98 |
%-- K*ra + add 1 overlap : |
| 99 |
kwx = repmat(ra,[1 1 nr]).*kwx; |
| 100 |
kwy = repmat(ra,[1 1 nr]).*kwy; |
| 101 |
v6X = split_C_cub(kwx,1); |
| 102 |
v6Y = split_C_cub(kwy,1); |
| 103 |
k6x = v6X(:,[2:nc+1],:,:); |
| 104 |
k6y = v6Y([2:nc+1],:,:,:); |
| 105 |
|
| 106 |
%----------------- |
| 107 |
clear v6X v6Y |
| 108 |
v6X = (k6x([2:nc+1],:,:,:) + k6x([1:nc],:,:,:))/2; |
| 109 |
v6Y = (k6y(:,[2:nc+1],:,:) + k6y(:,[1:nc],:,:))/2; |
| 110 |
|
| 111 |
psiX = zeros(6*nc,nc,nr); |
| 112 |
psiY = zeros(6*nc,nc,nr); |
| 113 |
|
| 114 |
for n = 1:6 |
| 115 |
is = 1+nc*(n-1);ie=nc*n; |
| 116 |
psiX([is:ie],[1:nc],[1:nr]) = v6X(:,:,:,n); |
| 117 |
psiY([is:ie],[1:nc],[1:nr]) = v6Y(:,:,:,n); |
| 118 |
end |
| 119 |
|
| 120 |
psiX = reshape(psiX,6*nc*nc,nr); |
| 121 |
psiY = reshape(psiY,6*nc*nc,nr); |
| 122 |
|
| 123 |
psiX_all(:,:,it) = mskw .* psiX ./ repmat(rAu,[1,nr]); |
| 124 |
psiY_all(:,:,it) = msks .* psiY ./ repmat(rAv,[1,nr]); |
| 125 |
|
| 126 |
end |
| 127 |
|
| 128 |
case 'Advc' |
| 129 |
|
| 130 |
psiX_all = reshape(d.GM_PsiX(1:6*nc,:,:,1:nt),6*nc*nc,nr,nt); |
| 131 |
psiY_all = reshape(d.GM_PsiY(:,1:nc,:,1:nt) ,6*nc*nc,nr,nt); |
| 132 |
|
| 133 |
otherwise |
| 134 |
disp(['C est Portnawak: GMform should be Skew or Advc: ',GMform]) |
| 135 |
end |
| 136 |
|
| 137 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 138 |
% Zonally integrate along broken lines % |
| 139 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 140 |
|
| 141 |
load(blkFile); |
| 142 |
ydim = length(bkl_Ylat); |
| 143 |
ylat = [-90,bkl_Ylat,90]; |
| 144 |
ufac= rem(bkl_Flg,2); |
| 145 |
vfac= fix(bkl_Flg/2); |
| 146 |
|
| 147 |
PsiB= zeros(ydim+2,nr+1,nt); |
| 148 |
|
| 149 |
for it = 1:nt |
| 150 |
for k = 1:nr |
| 151 |
psixt=dyg.*psiX_all(:,k,it).*mskWloc; |
| 152 |
psiyt=dxg.*psiY_all(:,k,it).*mskSloc; |
| 153 |
for jl = 1:ydim |
| 154 |
ie = bkl_Npts(jl); |
| 155 |
PsiB(jl+1,k,it) = sum( ufac(1:ie,jl).*psixt(bkl_IJuv(1:ie,jl)) ... |
| 156 |
+ vfac(1:ie,jl).*psiyt(bkl_IJuv(1:ie,jl)) ); |
| 157 |
end |
| 158 |
end |
| 159 |
end |
| 160 |
PsiB = 1e-6*PsiB; |
Parent Directory
| 
Revision Log
| 
Revision Graph