function [ub,vb]=calc_Bolus_CS(rAc,dXc,dYc,dXg,dYg,delR,kwx,kwy,hFacW,hFacS); % [ub,vb]=calc_Bolus_CS(rAc,dXc,dYc,dXg,dYg,delR,kwx,kwy,hFacW,hFacS); %-- GM-Bolus transp. : % from Kwx & Kwy (Skew-Flx => /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)/dYg/drF ; vb = d_k(psiY)/dXg/drF ; %--------------------------------------------------------------------- % set_axis [ncx,nc]=size(rAc); nr=length(delR); nPg=6*nc*nc; ncp=nc+1; rAu=dXc.*dYg; rAu=reshape(rAu,nPg,1); rAv=dYc.*dXg; rAv=reshape(rAv,nPg,1); %-- K*rAc + add 1 overlap : tmp=reshape(rAc,nPg,1)*ones(1,nr); kwx=reshape(kwx,nPg,nr); kwy=reshape(kwy,nPg,nr); tmp=0.5*tmp; kwx=tmp.*kwx; kwy=tmp.*kwy; kwx=reshape(kwx,ncx,nc,nr); kwy=reshape(kwy,ncx,nc,nr); v6X=split_C_cub(kwx,1); v6Y=split_C_cub(kwy,1); k6x=v6X(:,[2:ncp],:,:); k6y=v6Y([2:ncp],:,:,:); %--- recip_hFac & mask : hFacW=reshape(hFacW,nPg,nr); hFacS=reshape(hFacS,nPg,nr); rhFcW=1./hFacW; rhFcW(find(hFacW==0))=0; rhFcS=1./hFacS; rhFcS(find(hFacS==0))=0; hFacW=ceil(hFacW); hFacW=min(1,hFacW); hFacS=ceil(hFacS); hFacS=min(1,hFacS); %----------------- v6X=zeros(nc,nc,nr,6); v6X([1:nc],:,:,:)=k6x([2:ncp],:,:,:)+k6x([1:nc],:,:,:); v6X=v6X/2; psiX=zeros(ncx,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); end psiX=reshape(psiX,nPg,nr+1); psiX(:,[1:nr])=hFacW.*psiX(:,[1:nr]); ub=zeros(nPg,nr); ub=psiX(:,[2:nr+1])-psiX(:,[1:nr]); delR = reshape(delR,[1,length(delR)]); tmp=rAu*delR; ub=ub./tmp; ub=rhFcW.*ub; ub=reshape(ub,ncx,nc,nr); %----------------- v6Y=zeros(nc,nc,nr,6); v6Y(:,[1:nc],:,:)=k6y(:,[2:ncp],:,:)+k6y(:,[1:nc],:,:); v6Y=v6Y/2; psiY=zeros(ncx,nc,nr+1); for n=1:6 is=1+nc*(n-1);ie=nc*n; psiY([is:ie],[1:nc],[1:nr])=v6Y([1:nc],[1:nc],[1:nr],n); end psiY=reshape(psiY,nPg,nr+1); psiY(:,[1:nr])=hFacS.*psiY(:,[1:nr]); vb=zeros(nPg,nr); vb=psiY(:,[2:nr+1])-psiY(:,[1:nr]); tmp=rAv*delR; vb=vb./tmp; vb=rhFcS.*vb; vb=reshape(vb,ncx,nc,nr); %----------------- return