| 1 |
dfer |
1.2 |
function [ub,vb,wb]=calcBolusVelCube(d,g,GMform); |
| 2 |
dfer |
1.1 |
|
| 3 |
|
|
% [ub,vb] = calcBolusVelCube(d,g,GMform); |
| 4 |
|
|
% |
| 5 |
|
|
% Input arguments: |
| 6 |
|
|
% The incoming field data (d) and grid data (g) must be in a structured |
| 7 |
|
|
% array format (which is the format that comes from rdmnc): |
| 8 |
|
|
% d [Field data] Kwx,Kwy |
| 9 |
|
|
% g [Grid data ] drF,rA,dxC,dyC,dxG,dyG,HFacW,HFacS |
| 10 |
|
|
% GMform [string] GM form, 'Skew' or 'Advc' |
| 11 |
|
|
% |
| 12 |
|
|
% Output arguments: |
| 13 |
|
|
% ub, vb: GM-Bolus mass-weigthed velocity (i.e include |
| 14 |
|
|
% implicitly hFac factor) |
| 15 |
|
|
% |
| 16 |
|
|
% Comments: |
| 17 |
|
|
% For Skew-Flux form: uses Kwx & Kwy divided by 2 |
| 18 |
|
|
% compute Volume Stream function psiX,psiY above uVel.vVel |
| 19 |
|
|
% (at interface between 2 levels), units=m^3/s : |
| 20 |
|
|
% psiX=(rAc*kwx)_i / dXc ; psiY=(rAc*kwy)_j / dYc ; |
| 21 |
|
|
% and then the bolus velocity (m/s): |
| 22 |
|
|
% ub = d_k(psiX)/rAw/drF ; vb = d_k(psiY)/rAs/drF ; |
| 23 |
|
|
% |
| 24 |
|
|
%--------------------------------------------------------------------- |
| 25 |
|
|
|
| 26 |
|
|
|
| 27 |
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 28 |
|
|
% Prepare / reform incoming data % |
| 29 |
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 30 |
|
|
|
| 31 |
|
|
nc = size(g.XC,2); |
| 32 |
|
|
nr = length(g.drF); |
| 33 |
|
|
|
| 34 |
|
|
switch GMform |
| 35 |
|
|
case 'Skew' |
| 36 |
|
|
nt = size(d.GM_Kwx,4); |
| 37 |
|
|
case 'Advc' |
| 38 |
|
|
nt = size(d.GM_PsiX,4); |
| 39 |
|
|
end |
| 40 |
|
|
|
| 41 |
|
|
dr = g.drF; |
| 42 |
|
|
hw = reshape(g.HFacW(1:6*nc,1:nc,1:nr),[6*nc*nc,nr]); |
| 43 |
|
|
hs = reshape(g.HFacS(1:6*nc,1:nc,1:nr),[6*nc*nc,nr]); |
| 44 |
|
|
ra = reshape(g.rA(1:6*nc,1:nc) ,[6*nc*nc,1]); |
| 45 |
|
|
rAs = reshape(g.rAs(1:6*nc,1:nc) ,[6*nc*nc,1]); |
| 46 |
|
|
rAw = reshape(g.rAw(1:6*nc,1:nc) ,[6*nc*nc,1]); |
| 47 |
dfer |
1.2 |
dxc = reshape(g.dxC(1:6*nc,1:nc),[6*nc*nc,1]); |
| 48 |
|
|
dyc = reshape(g.dyC(1:6*nc,1:nc),[6*nc*nc,1]); |
| 49 |
|
|
dxg = reshape(g.dxG(1:6*nc,1:nc),[6*nc*nc,1]); |
| 50 |
|
|
dyg = reshape(g.dyG(1:6*nc,1:nc),[6*nc*nc,1]); |
| 51 |
dfer |
1.1 |
|
| 52 |
|
|
%--- recip_hFac & mask : |
| 53 |
|
|
mw=ceil(hw); mw=min(1,mw); |
| 54 |
|
|
ms=ceil(hs); ms=min(1,ms); |
| 55 |
|
|
|
| 56 |
dfer |
1.2 |
hw(find(hw==0))=Inf; |
| 57 |
|
|
hs(find(hs==0))=Inf; |
| 58 |
|
|
hw_recip=1./hw; %hw_recip(find(hw==0))=0; |
| 59 |
|
|
hs_recip=1./hs; %hs_recip(find(hs==0))=0; |
| 60 |
|
|
|
| 61 |
|
|
ub_all = zeros(6*nc,nc,nr,nt); |
| 62 |
|
|
vb_all = zeros(6*nc,nc,nr,nt); |
| 63 |
|
|
wb_all = zeros(6*nc,nc,nr,nt); |
| 64 |
dfer |
1.1 |
|
| 65 |
|
|
switch GMform |
| 66 |
|
|
|
| 67 |
|
|
%%%%%%% Skew-flux form case |
| 68 |
|
|
case 'Skew' |
| 69 |
|
|
kwx_all = reshape(d.GM_Kwx,[6*nc*nc,nr,nt]); |
| 70 |
|
|
kwy_all = reshape(d.GM_Kwy,[6*nc*nc,nr,nt]); |
| 71 |
|
|
|
| 72 |
|
|
kwx_all = 0.5*kwx_all; |
| 73 |
|
|
kwy_all = 0.5*kwy_all; |
| 74 |
|
|
|
| 75 |
|
|
for it = 1:nt |
| 76 |
|
|
kwx = kwx_all(:,:,it); |
| 77 |
|
|
kwy = kwy_all(:,:,it); |
| 78 |
|
|
|
| 79 |
|
|
%-- K*ra + add 1 overlap : |
| 80 |
|
|
kwx = (ra*ones(1,nr)).*kwx; |
| 81 |
|
|
kwy = (ra*ones(1,nr)).*kwy; |
| 82 |
|
|
kwx = reshape(kwx,[6*nc,nc,nr]); |
| 83 |
|
|
kwy = reshape(kwy,[6*nc,nc,nr]); |
| 84 |
|
|
v6X = split_C_cub(kwx,1); |
| 85 |
|
|
v6Y = split_C_cub(kwy,1); |
| 86 |
|
|
k6x = v6X(:,[2:nc+1],:,:); |
| 87 |
|
|
k6y = v6Y([2:nc+1],:,:,:); |
| 88 |
|
|
|
| 89 |
|
|
%----------------- |
| 90 |
|
|
v6X = zeros(nc,nc,nr,6); |
| 91 |
|
|
v6Y = zeros(nc,nc,nr,6); |
| 92 |
|
|
|
| 93 |
|
|
v6X([1:nc],:,:,:) = k6x([2:nc+1],:,:,:) + k6x([1:nc],:,:,:); |
| 94 |
|
|
v6Y(:,[1:nc],:,:) = k6y(:,[2:nc+1],:,:) + k6y(:,[1:nc],:,:); |
| 95 |
|
|
|
| 96 |
|
|
v6X = v6X/2; |
| 97 |
|
|
v6Y = v6Y/2; |
| 98 |
|
|
|
| 99 |
|
|
psiX = zeros(6*nc,nc,nr+1); |
| 100 |
|
|
psiY = zeros(6*nc,nc,nr+1); |
| 101 |
|
|
|
| 102 |
|
|
for n = 1:6 |
| 103 |
|
|
is = 1+nc*(n-1);ie=nc*n; |
| 104 |
|
|
psiX([is:ie],[1:nc],[1:nr]) = v6X([1:nc],[1:nc],[1:nr],n); |
| 105 |
|
|
psiY([is:ie],[1:nc],[1:nr]) = v6Y([1:nc],[1:nc],[1:nr],n); |
| 106 |
|
|
end |
| 107 |
|
|
|
| 108 |
|
|
psiX = reshape(psiX,6*nc*nc,nr+1); |
| 109 |
|
|
psiY = reshape(psiY,6*nc*nc,nr+1); |
| 110 |
|
|
|
| 111 |
|
|
psiX(:,[1:nr]) = mw.*psiX(:,[1:nr]); |
| 112 |
|
|
psiY(:,[1:nr]) = ms.*psiY(:,[1:nr]); |
| 113 |
|
|
ub = psiX(:,[2:nr+1]) - psiX(:,[1:nr]); |
| 114 |
|
|
vb = psiY(:,[2:nr+1]) - psiY(:,[1:nr]); |
| 115 |
|
|
|
| 116 |
|
|
dr = reshape(dr,[1,length(dr)]); |
| 117 |
|
|
% ub = reshape(hw_recip.*ub./(rAw*dr),[6*nc,nc,nr]); |
| 118 |
|
|
ub = reshape(ub./(rAw*dr),[6*nc,nc,nr]); |
| 119 |
|
|
% vb = reshape(hs_recip.*vb./(rAs*dr),[6*nc,nc,nr]); |
| 120 |
|
|
vb = reshape(vb./(rAs*dr),[6*nc,nc,nr]); |
| 121 |
dfer |
1.2 |
|
| 122 |
dfer |
1.1 |
ub_all(:,:,:,it) = ub; |
| 123 |
|
|
vb_all(:,:,:,it) = vb; |
| 124 |
dfer |
1.2 |
|
| 125 |
|
|
%%%%%%%%%%%%% |
| 126 |
|
|
[u6t,v6t] =split_UV_cub(ub,vb,0,1); |
| 127 |
|
|
[hw6t,hs6t]=split_UV_cub(reshape(hw,6*nc,nc,nr),reshape(hs,6*nc,nc,nr),0,1); |
| 128 |
|
|
[dy6t,dx6t]=split_UV_cub(reshape(dyg,6*nc,nc),reshape(dxg,6*nc,nc),0,1); |
| 129 |
|
|
|
| 130 |
|
|
%F6tX = u6t.*hw6t.*permute(repmat(dy6t,[1 1 1 nr]),[1 2 4 3]); |
| 131 |
|
|
%F6tY = v6t.*hs6t.*permute(repmat(dx6t,[1 1 1 nr]),[1 2 4 3]); |
| 132 |
|
|
F6tX = u6t.*permute(repmat(dy6t,[1 1 1 nr]),[1 2 4 3]); |
| 133 |
|
|
F6tY = v6t.*permute(repmat(dx6t,[1 1 1 nr]),[1 2 4 3]); |
| 134 |
|
|
|
| 135 |
|
|
Hdiv = zeros(nc,nc,nr,6); |
| 136 |
|
|
Hdiv = F6tX([2:nc+1],:,:,:) - F6tX([1:nc],:,:,:) ... |
| 137 |
|
|
+ F6tY(:,[2:nc+1],:,:) - F6tY(:,[1:nc],:,:); |
| 138 |
|
|
for k=1:nr |
| 139 |
|
|
Hdiv(:,:,k,:) = -Hdiv(:,:,k,:)*dr(k); |
| 140 |
|
|
end |
| 141 |
|
|
|
| 142 |
|
|
%psiX = zeros(6*nc,nc,nr); |
| 143 |
|
|
%for n = 1:6 |
| 144 |
|
|
% is = 1+nc*(n-1);ie=nc*n; |
| 145 |
|
|
% psiX([is:ie],[1:nc],[1:nr]) = Hdiv([1:nc],[1:nc],[1:nr],n); |
| 146 |
|
|
%end |
| 147 |
|
|
psiX = reshape(permute(Hdiv,[1 4 2 3]),6*nc,nc,nr); |
| 148 |
|
|
%wb = psiX ./ repmat(reshape(ra,6*nc,nc),[1 1 nr]); |
| 149 |
|
|
|
| 150 |
|
|
wb = zeros(6*nc,nc,nr+1); |
| 151 |
|
|
for k=nr:-1:1 |
| 152 |
|
|
wb(:,:,k)= psiX(:,:,k) ./ reshape(ra,6*nc,nc) + wb(:,:,k+1); |
| 153 |
|
|
end |
| 154 |
|
|
|
| 155 |
|
|
wb_all(:,:,:,it) = wb(:,:,1:30); |
| 156 |
|
|
%%%%%%%%%%%% |
| 157 |
dfer |
1.1 |
end |
| 158 |
|
|
|
| 159 |
|
|
%%%%%%% Advective form case |
| 160 |
|
|
case 'Advc' |
| 161 |
|
|
|
| 162 |
|
|
PsiX_all = reshape(d.GM_PsiX(1:6*nc,1:nc,1:nr,:),[6*nc*nc,nr,nt]); |
| 163 |
|
|
PsiY_all = reshape(d.GM_PsiY(1:6*nc,1:nc,1:nr,:),[6*nc*nc,nr,nt]); |
| 164 |
|
|
|
| 165 |
|
|
dr3d = ones(6*nc*nc,1)*reshape(dr,[1,length(dr)]); |
| 166 |
|
|
|
| 167 |
|
|
for it = 1:nt |
| 168 |
|
|
|
| 169 |
|
|
psiX = zeros(6*nc*nc,nr+1); |
| 170 |
|
|
psiY = zeros(6*nc*nc,nr+1); |
| 171 |
|
|
|
| 172 |
|
|
psiX(:,1:nr) = mw.*PsiX_all(:,:,it); |
| 173 |
|
|
psiY(:,1:nr) = ms.*PsiY_all(:,:,it); |
| 174 |
|
|
|
| 175 |
|
|
% psiX(:,[1:nr]) = mw.*psiX(:,[1:nr]); |
| 176 |
|
|
% psiY(:,[1:nr]) = ms.*psiY(:,[1:nr]); |
| 177 |
|
|
ub = psiX(:,[2:nr+1]) - psiX(:,[1:nr]); |
| 178 |
|
|
vb = psiY(:,[2:nr+1]) - psiY(:,[1:nr]); |
| 179 |
|
|
|
| 180 |
|
|
% ub = reshape(hw_recip.*ub./(rAw*dr),[6*nc,nc,nr]); |
| 181 |
|
|
% ub = reshape(ub./(rAw*dr),[6*nc,nc,nr]); |
| 182 |
|
|
ub = reshape(ub./dr3d,[6*nc,nc,nr]); |
| 183 |
|
|
% vb = reshape(hs_recip.*vb./(rAs*dr),[6*nc,nc,nr]); |
| 184 |
|
|
% vb = reshape(vb./(rAs*dr),[6*nc,nc,nr]); |
| 185 |
|
|
vb = reshape(vb./dr3d,[6*nc,nc,nr]); |
| 186 |
|
|
|
| 187 |
|
|
ub_all(:,:,:,it) = ub; |
| 188 |
|
|
vb_all(:,:,:,it) = vb; |
| 189 |
|
|
end |
| 190 |
|
|
|
| 191 |
|
|
otherwise |
| 192 |
|
|
disp('Ce portnawak') |
| 193 |
|
|
end |
| 194 |
|
|
|
| 195 |
|
|
ub = ub_all; |
| 196 |
|
|
vb = vb_all; |
| 197 |
dfer |
1.2 |
wb = wb_all; |
Parent Directory
| 
Revision Log
| 
Revision Graph