dfer/matlab_stuff/calcBolusVelCube.m

function [ub,vb]=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]);

%rAw=dxc.*dyg;
%rAs=dyc.*dxg;

%--- 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;

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;
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;
1	dfer	1.1	function [ub,vb]=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			%rAw=dxc.*dyg;
53			%rAs=dyc.*dxg;
54
55			%--- recip_hFac & mask :
56			mw=ceil(hw); mw=min(1,mw);
57			ms=ceil(hs); ms=min(1,ms);
58
59			%hw(find(hw==0))=Inf;
60			%hs(find(hs==0))=Inf;
61			%hw_recip=1./hw; %hw_recip(find(hw==0))=0;
62			%hs_recip=1./hs; %hs_recip(find(hs==0))=0;
63
64			switch GMform
65
66			%%%%%%% Skew-flux form case
67			case 'Skew'
68			kwx_all = reshape(d.GM_Kwx,[6ncnc,nr,nt]);
69			kwy_all = reshape(d.GM_Kwy,[6ncnc,nr,nt]);
70
71			kwx_all = 0.5*kwx_all;
72			kwy_all = 0.5*kwy_all;
73
74			for it = 1:nt
75			kwx = kwx_all(:,:,it);
76			kwy = kwy_all(:,:,it);
77
78			%-- K*ra + add 1 overlap :
79			kwx = (raones(1,nr)).kwx;
80			kwy = (raones(1,nr)).kwy;
81			kwx = reshape(kwx,[6*nc,nc,nr]);
82			kwy = reshape(kwy,[6*nc,nc,nr]);
83			v6X = split_C_cub(kwx,1);
84			v6Y = split_C_cub(kwy,1);
85			k6x = v6X(:,[2:nc+1],:,:);
86			k6y = v6Y([2:nc+1],:,:,:);
87
88			%-----------------
89			v6X = zeros(nc,nc,nr,6);
90			v6Y = zeros(nc,nc,nr,6);
91
92			v6X([1:nc],:,:,:) = k6x([2:nc+1],:,:,:) + k6x([1:nc],:,:,:);
93			v6Y(:,[1:nc],:,:) = k6y(:,[2:nc+1],:,:) + k6y(:,[1:nc],:,:);
94
95			v6X = v6X/2;
96			v6Y = v6Y/2;
97
98			psiX = zeros(6*nc,nc,nr+1);
99			psiY = zeros(6*nc,nc,nr+1);
100
101			for n = 1:6
102			is = 1+nc(n-1);ie=ncn;
103			psiX([is:ie],[1:nc],[1:nr]) = v6X([1:nc],[1:nc],[1:nr],n);
104			psiY([is:ie],[1:nc],[1:nr]) = v6Y([1:nc],[1:nc],[1:nr],n);
105			end
106
107			psiX = reshape(psiX,6ncnc,nr+1);
108			psiY = reshape(psiY,6ncnc,nr+1);
109
110			psiX(:,[1:nr]) = mw.*psiX(:,[1:nr]);
111			psiY(:,[1:nr]) = ms.*psiY(:,[1:nr]);
112			ub = psiX(:,[2:nr+1]) - psiX(:,[1:nr]);
113			vb = psiY(:,[2:nr+1]) - psiY(:,[1:nr]);
114
115			dr = reshape(dr,[1,length(dr)]);
116			% ub = reshape(hw_recip.ub./(rAwdr),[6*nc,nc,nr]);
117			ub = reshape(ub./(rAwdr),[6nc,nc,nr]);
118			% vb = reshape(hs_recip.vb./(rAsdr),[6*nc,nc,nr]);
119			vb = reshape(vb./(rAsdr),[6nc,nc,nr]);
120
121			ub_all(:,:,:,it) = ub;
122			vb_all(:,:,:,it) = vb;
123			end
124
125			%%%%%%% Advective form case
126			case 'Advc'
127
128			PsiX_all = reshape(d.GM_PsiX(1:6nc,1:nc,1:nr,:),[6nc*nc,nr,nt]);
129			PsiY_all = reshape(d.GM_PsiY(1:6nc,1:nc,1:nr,:),[6nc*nc,nr,nt]);
130
131			dr3d = ones(6ncnc,1)*reshape(dr,[1,length(dr)]);
132
133			for it = 1:nt
134
135			psiX = zeros(6ncnc,nr+1);
136			psiY = zeros(6ncnc,nr+1);
137
138			psiX(:,1:nr) = mw.*PsiX_all(:,:,it);
139			psiY(:,1:nr) = ms.*PsiY_all(:,:,it);
140
141			% psiX(:,[1:nr]) = mw.*psiX(:,[1:nr]);
142			% psiY(:,[1:nr]) = ms.*psiY(:,[1:nr]);
143			ub = psiX(:,[2:nr+1]) - psiX(:,[1:nr]);
144			vb = psiY(:,[2:nr+1]) - psiY(:,[1:nr]);
145
146			% ub = reshape(hw_recip.ub./(rAwdr),[6*nc,nc,nr]);
147			% ub = reshape(ub./(rAwdr),[6nc,nc,nr]);
148			ub = reshape(ub./dr3d,[6*nc,nc,nr]);
149			% vb = reshape(hs_recip.vb./(rAsdr),[6*nc,nc,nr]);
150			% vb = reshape(vb./(rAsdr),[6nc,nc,nr]);
151			vb = reshape(vb./dr3d,[6*nc,nc,nr]);
152
153			ub_all(:,:,:,it) = ub;
154			vb_all(:,:,:,it) = vb;
155			end
156
157			otherwise
158			disp('Ce portnawak')
159			end
160
161			ub = ub_all;
162			vb = vb_all;