anchor.diff:      allocate(presF1(DIMS,grid%nC))
       allocate(condQ(grid%nC))
       allocate(flowRho(grid%nP))
  -    allocate(flowRho1(grid%nP))
       allocate(advRho(grid%nC))
       allocate(advRhou(DIMS,grid%nC))
       allocate(flowRhoH(grid%nP))
  
positive.code:       end forall
      end if
      call setBC()
      call findMassFlow(grid,rhou,flowRho)
      call findGrad(grid,p,gradP)
      call addRhieChow(grid,p,gradP,dt,flowRho)
      call findAdv(grid,flowRho,advRho)
    end subroutine
    
  
negative.code:     end if
      tmpP(:)=p(1:grid%nC)
      call setBC()
      call pressureEq%solve(tmpP,info=info)
      needRetry=info<0
      call pressureEq%getNIt(n)
  