anchor.diff:        rf(2)=-grid%IntfNorm(1,i)*tf(3)+grid%IntfNorm(3,i)*tf(1)
         rf(3)=grid%IntfNorm(1,i)*tf(2)-grid%IntfNorm(2,i)*tf(1)
         Afs=grid%IntfArea(i)/dot_product(sf,grid%IntfNorm(:,i))
  -      flowRate(:)=Afs*((v(:,n)-v(:,m))/dPF&
  -      &           -matmul(transpose(gradIntf(:,:,i)),dot_product(sf,tf)*tf+dot_product(sf,rf)*rf))
  +      flowRate(:)=gammIntf(i)*(Afs*((v(:,n)-v(:,m))/dPF&
  +      &                        -matmul(transpose(gradIntf(:,:,i)),&
  +      &                                dot_product(sf,tf)*tf+dot_product(sf,rf)*rf)))
         findDiffusSDVect(:,m)=findDiffusSDVect(:,m)+flowRate(:)
         findDiffusSDVect(:,n)=findDiffusSDVect(:,n)-flowRate(:)
       end do
       deallocate(gradIntf)
  +    deallocate(gammIntf)
     end function
     
     !> find diffusion driven by scalar v through interface using surface decomposition scheme
  -  !> \f[ \int_{intf} \nabla v \cdot \hat{n} dA \f]
  -  function findDiffusSDScal(v,grid,grad)
  +  !> \f[ \int_{intf} \Gamma \nabla v \cdot \hat{n} dA \f]
  +  function findDiffusSDScal(gamm,gammbind,v,grid,grad)
       use moduleGrid
  +    double precision,intent(in)::gamm(:) !< the diffusivity
  +    integer,intent(in)::gammbind !< bind diffusivity with node/block/interface
       double precision,intent(in)::v(:) !< the scalar which drives the diffusion
       type(typeGrid),intent(inout)::grid !< the grid on which v is defined
       double precision,intent(in)::grad(DIMS,size(v)) !< the gradient of v
  
positive.code: contains
    
    !> find diffusion driven by vector v through interface using orthogonal scheme
    !> \f[ \int_{intf} \nabla \mathbf{v} \cdot \hat{n} dA \f]
    function findDiffusORTHVect(v,grid)
      use moduleGrid
      double precision,intent(in)::v(:,:) !< the vector which drives the diffusion
      type(typeGrid),intent(inout)::grid !< the grid on which v is defined
      double precision findDiffusORTHVect(size(v,1),grid%nBlock) !< increment due to diffusion
      double precision flowRate(size(v,1))
      
      call grid%updateIntfPos()
      call grid%updateIntfArea()
      call grid%updateBlockPos()
      findDiffusORTHVect(:,:)=0d0
      do i=1,grid%nIntf
        m=grid%IntfNeibBlock(1,i)
        n=grid%IntfNeibBlock(2,i)
        flowRate(:)=(v(:,n)-v(:,m))*grid%IntfArea(i)/norm2(grid%BlockPos(:,n)-grid%BlockPos(:,m))
        findDiffusORTHVect(:,m)=findDiffusORTHVect(:,m)+flowRate(:)
        findDiffusORTHVect(:,n)=findDiffusORTHVect(:,n)-flowRate(:)
      end do
    end function
    
    !> find diffusion driven by scalar v through interface using orthogonal scheme
    !> \f[ \int_{intf} \nabla v \cdot \hat{n} dA \f]
    function findDiffusORTHScal(v,grid)
      use moduleGrid
      double precision,intent(in)::v(:) !< the scalar which drives the diffusion
      type(typeGrid),intent(inout)::grid !< the grid on which v is defined
      double precision findDiffusORTHScal(grid%nBlock) !< increment due to diffusion
      double precision vv(1,size(v)),vrst(1,grid%nBlock)
      
      vv(1,:)=v(:)
      vrst=findDiffusORTHVect(vv,grid)
      findDiffusORTHScal(:)=vrst(1,:)
    end function
    
    !> find diffusion driven by vector v through interface using surface decomposition scheme
    !> \f[ \int_{intf} \nabla \mathbf{v} \cdot \hat{n} dA \f]
    function findDiffusSDVect(v,grid,grad)
      use moduleGrid
      use moduleInterpolation
      double precision,intent(in)::v(:,:) !< the vector which drives the diffusion
      type(typeGrid),intent(inout)::grid !< the grid on which v is defined
      double precision,intent(in)::grad(DIMS,size(v,1),size(v,2)) !< the gradient of v
      double precision findDiffusSDVect(size(v,1),grid%nBlock) !< increment due to diffusion
      double precision flowRate(size(v,1)),dPF,sf(DIMS),tf(DIMS),rf(DIMS),Afs
      double precision,allocatable::gradIntf(:,:,:)
      
      call grid%updateIntfPos()
      call grid%updateIntfArea()
  
negative.code:     
      vv(1,:)=v(:)
      vgrad(:,1,:)=grad(:,:)
      vrst=findDiffusSDVect(gamm,gammbind,vv,grid,grad)
      findDiffusSDScal(:)=vrst(1,:)
    end function
    
  