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;