anchor.diff:      p(1:grid%nC)=x(1:grid%nC)
       call setBC()
       call findGrad(grid,p,gradP)
  -    forall(i=1:grid%nE)
  -      fluxRho(:,i)=rho(i)*u(:,i)
  -    end forall
  -    call findAdv(grid,rho,fluxRho,flowRho)
  -    call addRhieChow(grid,rho,p,gradP,rho,dt,flowRho)
  +    if(.not.allocated(tmpFlowRho))then
  +      allocate(tmpFlowRho(grid%nC))
  +    else if(size(tmpFlowRho)/=grid%nC)then
  +      deallocate(tmpFlowRho)
  +      allocate(tmpFlowRho(grid%nC))
  +    end if
  +    tmpFlowRho(:)=flowRho(:)
  +    call addRhieChow(grid,rho,p,gradP,rho,dt,tmpFlowRho)
       call findDiff(grid,p,[(1d0,i=1,grid%nC)],laP)
       ! blah blah
       pressureRHS=0
  
positive.code:   
    !> solve the pressure, during which the laP is updated
    subroutine solvePressure()
      double precision,allocatable,save::tmpP(:)
      
      if(.not.allocated(tmpP))then
  
negative.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
      double precision,allocatable,save::tmpFlowRho(:)
      
      call associateVector(oldVector,x)
      call associateVector(newVector,y)
  