anchor.diff:        allocate(flow(grid%nP))
       end if
       flow(:)=0d0
  +    if(present(sensP))then
  +      if(.not.allocated(sensP))then
  +        allocate(sensP(2,grid%nP))
  +      end if
  +      sensP(:,:)=0d0
  +    end if
       !$omp parallel do default(shared)&
       !$omp& private(m,n,up,dn,rhoUp,rhoDn,dRho,r,rhof,vMN,vMP,vPN,eps,flux)
       do i=1,grid%nP
  
positive.code:         &       -dt/(1d0/rho(m)+1d0/rho(n))*&
          &        (presF(:,m)/grid%v(m)/rho(m)+presF(:,n)/grid%v(n)/rho(n))&
          &       +dt*grid%normP(:,i)*(p(m)-p(n))/dot_product(vMN,grid%normP(:,i))
        else ! boundary pairs
          flux(:)=rhoUp*0.5d0*(u(:,m)+u(:,n))
          if(abs(dot_product(grid%normP(:,i),flux(:)))>tiny(1d0))then
            vMP(:)=grid%pP(:,i)-grid%p(:,m)
            flux(:)=flux(:)-0.5d0*dt*presF(:,m)/grid%v(m)&
            &       +0.5d0*dt*grid%normP(:,i)*(p(m)-p(n))/(2d0*dot_product(vMP,grid%normP(:,i)))
          end if ! Rhie-Chow does not initiate penetration through boundary pair
        end if
        flow(i)=grid%aP(i)*dot_product(grid%normP(:,i),flux(:))
  
negative.code:     double precision,intent(in)::dt !< time step size
      double precision,allocatable,intent(inout)::flow(:) !< mass flow rate output
      double precision,intent(in),optional::gradRho(:,:) !< gradient of rho
      double precision,allocatable,intent(inout),optional::sensP(:,:) !< flow's sensitivity on p
      integer::up,dn
      double precision::vMN(DIMS),vMP(DIMS),vPN(DIMS),flux(DIMS),eps,rhoUp,rhoDn,dRho,r,rhof
      logical::sameDirection
  