anchor.diff:        deallocate(tmpP)
         allocate(tmpP(grid%nC))
       end if
  +    tmpP(:)=p(1:grid%nC)
  +    forall(i=1:grid%nE)
  +      fluxRho(:,i)=rho(i)*u(:,i)
  +    end forall
  +    call findAdv(grid,rho,fluxRho,flowRho) ! vary only the Rhie-Chow part of flowRho in RHS
       call pressureEq%solve(tmpP)
       p(1:grid%nC)=tmpP(:)
     end subroutine
  
positive.code:     integer(C_INT)::pressureRHS !< error code
      double precision,pointer::x(:) !< fortran pointer associated with oldVector
      double precision,pointer::y(:) !< fortran pointer associated with newVector
      
      call associateVector(oldVector,x)
      call associateVector(newVector,y)
  
negative.code:   
    !> solve the pressure, during which the laP is updated
    subroutine solvePressure()
      use modAdvection
      double precision,allocatable,save::tmpP(:)
      
      if(.not.allocated(tmpP))then
  