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:       deallocate(tmpP)
        allocate(tmpP(grid%nC))
      end if
      call pressureEq%solve(tmpP)
      p(1:grid%nC)=tmpP(:)
    end subroutine
  
negative.code:     use iso_c_binding
      use modNumerics
      use modGradient
      use modDiffusion
      use modRhieChow
      type(C_PTR),value::oldVector !< old N_Vector
  