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	function [ub,vb,wb]=calcBolusVelCube(d,g,GMform);
2
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	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
52	%--- recip_hFac & mask :
53	mw=ceil(hw); mw=min(1,mw);
54	ms=ceil(hs); ms=min(1,ms);
55
56	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
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
122	ub_all(:,:,:,it) = ub;
123	vb_all(:,:,:,it) = vb;
124
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	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	wb = wb_all;