anchor.diff:      end if
       tmpP(:)=p(1:grid%nC)
       call setBC()
  -    call findAdv(grid,rho,rhou,advRho)
  +    call findMassFlow(grid,rhou,flowRho)
       call pressureEq%solve(tmpP,info=info)
       needRetry=info<0
       call pressureEq%getNIt(n)
  
positive.code:       deallocate(tmpAdvRho)
        allocate(tmpAdvRho(grid%nC))
      end if
      tmpAdvRho(:)=advRho(:)
      call addRhieChow(grid,rho,p,gradP,rho,dt,tmpAdvRho)
      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)*(tmpAdvRho(i)+dt*(laP(i)-laP1(i))))
  
negative.code:   double precision,allocatable::condQ(:) !< heat conduction into cell [W]
    double precision,allocatable::carrier(:,:) !< {rho,rhou,rhoH} of each cell [*]
    double precision,allocatable::advRho(:) !< advection rate of rho [kg/s]
    double precision,allocatable::flowRho(:) !< flow rate of rho through pairs [kg/s]
    double precision,allocatable::fluxRhou(:,:,:) !< flux of rhou [kg/m/s^2]
    double precision,allocatable::advRhou(:,:) !< advection rate of rhou into cell [kg*m/s^2]
    double precision,allocatable::fluxRhoH(:,:) !< flux of rhoH [J/m^2/s]
  