dfer/matlab_stuff/calcBolusVelCube.m

function [ub,vb,wb]=calcBolusVelCube(d,g,GMform);

% [ub,vb] = calcBolusVelCube(d,g,GMform);
%
% Input arguments:  
%   The incoming field data (d) and grid data (g) must be in a structured
%   array format (which is the format that comes from rdmnc):
%       d  [Field data]   Kwx,Kwy
%       g  [Grid data ]   drF,rA,dxC,dyC,dxG,dyG,HFacW,HFacS
%       GMform  [string]  GM form, 'Skew' or 'Advc'
%
% Output arguments:
%       ub, vb: GM-Bolus mass-weigthed velocity (i.e include 
%               implicitly hFac factor)  
%
% Comments:
%       For Skew-Flux form: uses Kwx & Kwy divided by 2 
%       compute Volume Stream function psiX,psiY above uVel.vVel 
%       (at interface between 2 levels), units=m^3/s :
%       psiX=(rAc*kwx)_i / dXc ; psiY=(rAc*kwy)_j / dYc ;
%       and then the bolus velocity (m/s):
%       ub = d_k(psiX)/rAw/drF ; vb = d_k(psiY)/rAs/drF ;
%
%---------------------------------------------------------------------


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                     Prepare / reform  incoming data                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nc = size(g.XC,2);
nr = length(g.drF);

switch GMform
case 'Skew'
    nt = size(d.GM_Kwx,4);
case 'Advc'
    nt = size(d.GM_PsiX,4);
end

dr  = g.drF;
hw = reshape(g.HFacW(1:6*nc,1:nc,1:nr),[6*nc*nc,nr]);
hs = reshape(g.HFacS(1:6*nc,1:nc,1:nr),[6*nc*nc,nr]);
ra  = reshape(g.rA(1:6*nc,1:nc) ,[6*nc*nc,1]);
rAs = reshape(g.rAs(1:6*nc,1:nc) ,[6*nc*nc,1]);
rAw = reshape(g.rAw(1:6*nc,1:nc) ,[6*nc*nc,1]);
dxc = reshape(g.dxC(1:6*nc,1:nc),[6*nc*nc,1]);
dyc = reshape(g.dyC(1:6*nc,1:nc),[6*nc*nc,1]);
dxg = reshape(g.dxG(1:6*nc,1:nc),[6*nc*nc,1]);
dyg = reshape(g.dyG(1:6*nc,1:nc),[6*nc*nc,1]);

%--- recip_hFac & mask :
mw=ceil(hw); mw=min(1,mw);
ms=ceil(hs); ms=min(1,ms);

hw(find(hw==0))=Inf;
hs(find(hs==0))=Inf;
hw_recip=1./hw; %hw_recip(find(hw==0))=0;
hs_recip=1./hs; %hs_recip(find(hs==0))=0;

ub_all = zeros(6*nc,nc,nr,nt);
vb_all = zeros(6*nc,nc,nr,nt);
wb_all = zeros(6*nc,nc,nr,nt);

switch GMform

%%%%%%% Skew-flux form case
case 'Skew'
kwx_all = reshape(d.GM_Kwx,[6*nc*nc,nr,nt]);
kwy_all = reshape(d.GM_Kwy,[6*nc*nc,nr,nt]);

kwx_all = 0.5*kwx_all;
kwy_all = 0.5*kwy_all;

for it = 1:nt
    kwx = kwx_all(:,:,it);
    kwy = kwy_all(:,:,it);

    %-- K*ra + add 1 overlap :
    kwx = (ra*ones(1,nr)).*kwx;
    kwy = (ra*ones(1,nr)).*kwy;
    kwx = reshape(kwx,[6*nc,nc,nr]);
    kwy = reshape(kwy,[6*nc,nc,nr]);
    v6X = split_C_cub(kwx,1);
    v6Y = split_C_cub(kwy,1);
    k6x = v6X(:,[2:nc+1],:,:);
    k6y = v6Y([2:nc+1],:,:,:);

    %----------------- 
    v6X = zeros(nc,nc,nr,6);
    v6Y = zeros(nc,nc,nr,6);

    v6X([1:nc],:,:,:) = k6x([2:nc+1],:,:,:) + k6x([1:nc],:,:,:);
    v6Y(:,[1:nc],:,:) = k6y(:,[2:nc+1],:,:) + k6y(:,[1:nc],:,:);

    v6X = v6X/2;
    v6Y = v6Y/2;

    psiX = zeros(6*nc,nc,nr+1);
    psiY = zeros(6*nc,nc,nr+1);

    for n = 1:6
        is = 1+nc*(n-1);ie=nc*n;
        psiX([is:ie],[1:nc],[1:nr]) = v6X([1:nc],[1:nc],[1:nr],n);
        psiY([is:ie],[1:nc],[1:nr]) = v6Y([1:nc],[1:nc],[1:nr],n);
    end

    psiX = reshape(psiX,6*nc*nc,nr+1);
    psiY = reshape(psiY,6*nc*nc,nr+1);

    psiX(:,[1:nr]) = mw.*psiX(:,[1:nr]);
    psiY(:,[1:nr]) = ms.*psiY(:,[1:nr]);
    ub = psiX(:,[2:nr+1]) - psiX(:,[1:nr]);
    vb = psiY(:,[2:nr+1]) - psiY(:,[1:nr]);

    dr = reshape(dr,[1,length(dr)]);
%    ub = reshape(hw_recip.*ub./(rAw*dr),[6*nc,nc,nr]);
    ub = reshape(ub./(rAw*dr),[6*nc,nc,nr]);
%    vb = reshape(hs_recip.*vb./(rAs*dr),[6*nc,nc,nr]);
    vb = reshape(vb./(rAs*dr),[6*nc,nc,nr]);

    ub_all(:,:,:,it) = ub;
    vb_all(:,:,:,it) = vb;

%%%%%%%%%%%%%
    [u6t,v6t]  =split_UV_cub(ub,vb,0,1);
    [hw6t,hs6t]=split_UV_cub(reshape(hw,6*nc,nc,nr),reshape(hs,6*nc,nc,nr),0,1);
    [dy6t,dx6t]=split_UV_cub(reshape(dyg,6*nc,nc),reshape(dxg,6*nc,nc),0,1);

    %F6tX = u6t.*hw6t.*permute(repmat(dy6t,[1 1 1 nr]),[1 2 4 3]);
    %F6tY = v6t.*hs6t.*permute(repmat(dx6t,[1 1 1 nr]),[1 2 4 3]);
    F6tX = u6t.*permute(repmat(dy6t,[1 1 1 nr]),[1 2 4 3]);
    F6tY = v6t.*permute(repmat(dx6t,[1 1 1 nr]),[1 2 4 3]);

    Hdiv = zeros(nc,nc,nr,6);
    Hdiv = F6tX([2:nc+1],:,:,:) - F6tX([1:nc],:,:,:) ...
         + F6tY(:,[2:nc+1],:,:) - F6tY(:,[1:nc],:,:);
    for k=1:nr
      Hdiv(:,:,k,:) = -Hdiv(:,:,k,:)*dr(k);
    end

    %psiX = zeros(6*nc,nc,nr);
    %for n = 1:6
    %  is = 1+nc*(n-1);ie=nc*n;
    %  psiX([is:ie],[1:nc],[1:nr]) = Hdiv([1:nc],[1:nc],[1:nr],n);
    %end
    psiX = reshape(permute(Hdiv,[1 4 2 3]),6*nc,nc,nr);
    %wb = psiX ./ repmat(reshape(ra,6*nc,nc),[1 1 nr]);

    wb = zeros(6*nc,nc,nr+1);
    for k=nr:-1:1
     wb(:,:,k)= psiX(:,:,k) ./ reshape(ra,6*nc,nc) + wb(:,:,k+1);
    end

    wb_all(:,:,:,it) = wb(:,:,1:30);
%%%%%%%%%%%%
end

%%%%%%% Advective form case
case 'Advc'

PsiX_all = reshape(d.GM_PsiX(1:6*nc,1:nc,1:nr,:),[6*nc*nc,nr,nt]);
PsiY_all = reshape(d.GM_PsiY(1:6*nc,1:nc,1:nr,:),[6*nc*nc,nr,nt]);

dr3d = ones(6*nc*nc,1)*reshape(dr,[1,length(dr)]);

for it = 1:nt

    psiX = zeros(6*nc*nc,nr+1);
    psiY = zeros(6*nc*nc,nr+1);

    psiX(:,1:nr) = mw.*PsiX_all(:,:,it);
    psiY(:,1:nr) = ms.*PsiY_all(:,:,it);

%    psiX(:,[1:nr]) = mw.*psiX(:,[1:nr]);
%    psiY(:,[1:nr]) = ms.*psiY(:,[1:nr]);
    ub = psiX(:,[2:nr+1]) - psiX(:,[1:nr]);
    vb = psiY(:,[2:nr+1]) - psiY(:,[1:nr]);

%    ub = reshape(hw_recip.*ub./(rAw*dr),[6*nc,nc,nr]);
%    ub = reshape(ub./(rAw*dr),[6*nc,nc,nr]);
    ub = reshape(ub./dr3d,[6*nc,nc,nr]);
%    vb = reshape(hs_recip.*vb./(rAs*dr),[6*nc,nc,nr]);
%    vb = reshape(vb./(rAs*dr),[6*nc,nc,nr]);
    vb = reshape(vb./dr3d,[6*nc,nc,nr]);
 
    ub_all(:,:,:,it) = ub;
    vb_all(:,:,:,it) = vb;
end

otherwise 
   disp('Ce portnawak')
end

ub = ub_all;
vb = vb_all;
wb = wb_all;
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=(rAckwx)_i / dXc ; psiY=(rAckwy)_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:6nc,1:nc,1:nr),[6nc*nc,nr]);
43			hs = reshape(g.HFacS(1:6nc,1:nc,1:nr),[6nc*nc,nr]);
44			ra = reshape(g.rA(1:6nc,1:nc) ,[6nc*nc,1]);
45			rAs = reshape(g.rAs(1:6nc,1:nc) ,[6nc*nc,1]);
46			rAw = reshape(g.rAw(1:6nc,1:nc) ,[6nc*nc,1]);
47	dfer	1.2	dxc = reshape(g.dxC(1:6nc,1:nc),[6nc*nc,1]);
48			dyc = reshape(g.dyC(1:6nc,1:nc),[6nc*nc,1]);
49			dxg = reshape(g.dxG(1:6nc,1:nc),[6nc*nc,1]);
50			dyg = reshape(g.dyG(1:6nc,1:nc),[6nc*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,[6ncnc,nr,nt]);
70			kwy_all = reshape(d.GM_Kwy,[6ncnc,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 = (raones(1,nr)).kwx;
81			kwy = (raones(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=ncn;
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,6ncnc,nr+1);
109			psiY = reshape(psiY,6ncnc,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./(rAwdr),[6*nc,nc,nr]);
118			ub = reshape(ub./(rAwdr),[6nc,nc,nr]);
119			% vb = reshape(hs_recip.vb./(rAsdr),[6*nc,nc,nr]);
120			vb = reshape(vb./(rAsdr),[6nc,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,6nc,nc,nr),reshape(hs,6nc,nc,nr),0,1);
128			[dy6t,dx6t]=split_UV_cub(reshape(dyg,6nc,nc),reshape(dxg,6nc,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=ncn;
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:6nc,1:nc,1:nr,:),[6nc*nc,nr,nt]);
163			PsiY_all = reshape(d.GM_PsiY(1:6nc,1:nc,1:nr,:),[6nc*nc,nr,nt]);
164
165			dr3d = ones(6ncnc,1)*reshape(dr,[1,length(dr)]);
166
167			for it = 1:nt
168
169			psiX = zeros(6ncnc,nr+1);
170			psiY = zeros(6ncnc,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./(rAwdr),[6*nc,nc,nr]);
181			% ub = reshape(ub./(rAwdr),[6nc,nc,nr]);
182			ub = reshape(ub./dr3d,[6*nc,nc,nr]);
183			% vb = reshape(hs_recip.vb./(rAsdr),[6*nc,nc,nr]);
184			% vb = reshape(vb./(rAsdr),[6nc,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;