function [bolusU,bolusV,bolusW]=calc_bolus_mod(GM_PsiX,GM_PsiY,nfx,nfy);
%object:    compute bolus velocty field (bolusU,bolusV,bolusW)
%           from gm streamfunction (GM_PsiX,GM_PsiY).
%input:     GM_PsiX,GM_PsiY
%output:    bolusU,bolusV,bolusW

global mygrid
nr=length(mygrid.RC);
sz=size(mygrid.RAC);
nx=sz(1);ny=sz(2);

%replace NaNs with 0s:
GM_PsiX(isnan(GM_PsiX))=0;
GM_PsiY(isnan(GM_PsiY))=0;

%compute bolus velocity:
GM_PsiX=cat(3,GM_PsiX,zeros(nx,ny));	%nr+1
GM_PsiY=cat(3,GM_PsiY,zeros(nx,ny));	%nr+1
%GM_PsiX(:,:,nr+1)=0*GM_PsiX(:,:,end);
%GM_PsiY(:,:,nr+1)=0*GM_PsiY(:,:,end);
bolusU=0*mygrid.hFacW;
bolusV=0*mygrid.hFacS;
for k=1:nr;
    bolusU(:,:,k)=(GM_PsiX(:,:,k+1)-GM_PsiX(:,:,k))/mygrid.DRF(k);
    bolusV(:,:,k)=(GM_PsiY(:,:,k+1)-GM_PsiY(:,:,k))/mygrid.DRF(k);
end;
bolusU=bolusU.*(~isnan(mygrid.mskW));
bolusV=bolusV.*(~isnan(mygrid.mskS));

%and its vertical part
%   (seems correct, leading to 0 divergence)
tmp_x=GM_PsiX(:,:,1:nr).*repmat(mygrid.DYG,[1 1 nr]);
tmp_y=GM_PsiY(:,:,1:nr).*repmat(mygrid.DXG,[1 1 nr]);
%[tmp_x,tmp_y]=exch_UV(tmp_x,tmp_y);
[tmp_x,tmp_y]=llc_uv_addpadding(tmp_x,tmp_y,nfx,nfy);
tmp_a=repmat(mygrid.RAC,[1 1 nr]);tmp_a=get_aste_faces(mygrid.RAC,nfx,nfy);
for iFace=1:mygrid.nFaces;
    tmp_w{iFace}=0*tmp_x{iFace};
    tmp_w{iFace}=( tmp_x{iFace}(2:end,:,:)-tmp_x{iFace}(1:end-1,:,:) )+...
        ( tmp_y{iFace}(:,2:end,:)-tmp_y{iFace}(:,1:end-1,:) );
    tmp_w{iFace}=tmp_w{iFace}./tmp_a{iFace};
end;

%put bolusW back to compact:
if(nfx(2)>0);
  if(nfx(4)>0);
    if(nfx(3)>0);
      bolusW=cat(2,tmp_w{1},tmp_w{2},tmp_w{3},reshape(tmp_w{4},nfy(4),nfx(4),nr),reshape(tmp_w{5},nfy(5),nfx(5),nr));
    else;
      bolusW=cat(2,tmp_w{1},tmp_w{2},reshape(tmp_w{4},nfy(4),nfx(4),nr),reshape(tmp_w{5},nfy(5),nfx(5),nr));
    end;
  else;
    if(nfx(3)>0);
      bolusW=cat(2,tmp_w{1},tmp_w{2},tmp_w{3},reshape(tmp_w{5},nfy(5),nfx(5),nr));
    else;
      bolusW=cat(2,tmp_w{1},tmp_w{2},reshape(tmp_w{5},nfy(5),nfx(5),nr));
    end;
  end;
elseif(nfx(2)==0);
  if(nfx(4)>0);
    if(nfx(3)>0);
      bolusW=cat(2,tmp_w{1},tmp_w{3},reshape(tmp_w{4},nfy(4),nfx(4),nr),reshape(tmp_w{5},nfy(5),nfx(5),nr));
    else;
      bolusW=cat(2,tmp_w{1},reshape(tmp_w{4},nfy(4),nfx(4),nr),reshape(tmp_w{5},nfy(5),nfx(5),nr));
    end;
  else;
    if(nfx(3)>0);
      bolusW=cat(2,tmp_w{1},tmp_w{3},reshape(tmp_w{5},nfy(5),nfx(5),nr));
    else;
      bolusW=cat(2,tmp_w{1},reshape(tmp_w{5},nfy(5),nfx(5),nr));
    end;
  end;
end;

bolusW=bolusW.*(~isnan(mygrid.mskC));
