anchor.diff:      end if
       tmpP(:)=p(1:grid%nC)
       call setBC()
  -    call findMassFlow(grid,rhou,flowRho)
  +    call findMassFlow(grid,rhou,flowRho1)
       call pressureEq%solve(tmpP,info=info)
       needRetry=info<0
       call pressureEq%getNIt(n)
  
positive.code:     allocate(presF1(DIMS,grid%nC))
      allocate(condQ(grid%nC))
      allocate(flowRho(grid%nP))
      allocate(advRho(grid%nC))
      allocate(advRho1(grid%nC))
      allocate(advRhou(DIMS,grid%nC))
      allocate(advRhoH(grid%nC))
      ! load initial condition
  
negative.code:     p(1:grid%nC)=x(1:grid%nC)
      call setBC()
      call findGrad(grid,p,gradP)
      flowRho(:)=flowRho1(:)
      call addRhieChow(grid,p,gradP,dt,flowRho)
      call findAdv(grid,flowRho,advRho)
      call findDiff(grid,p,[(1d0,i=1,grid%nC)],laP)
      forall(i=1:grid%nC)
        y(i)=p(i)-Rgas*temp(i)*(rho0(i)+dt/grid%v(i)*(advRho(i)+dt*(laP(i)-laP1(i))))
      end forall
      pressureRHS=merge(1,0,any(ieee_is_nan(y).or.(.not.ieee_is_finite(y))))
      if(c_associated(dat))then
  