anchor.diff:      double precision,parameter::CFL_FLOW=0.5d0
       double precision,parameter::CFL_DIFFUSION=0.5d0
       
  -    ! TODO calculate dt, cantera sound speed
       dt=min(dt,minval(CFL_ACCOUSTIC*grid%v(:)**(1d0/3d0)&
  -    &                     /sqrt(1.4d0*p(1:grid%nC)/rho(1:grid%nC))))
  +    &                     /sqrt(gamm*p(1:grid%nC)/rho(1:grid%nC))))
       dt=min(dt,minval(CFL_FLOW*grid%v(:)**(1d0/3d0)&
       &                /norm2(u(:,1:grid%nC),1)))
       dt=min(dt,minval(CFL_DIFFUSION*grid%v(:)**(2d0/3d0)&
  
positive.code:     end do
      call recoverState(p,u,temp,Y,rho,rhou,rhoE)
      t=0d0
      dt=1d-7 ! TODO replace the good old 1d-7 with user input
      tNext=tInt
      iOut=0
      nRetry=0
  
negative.code:   public
    
    integer,parameter::DIMS=3 !< three dimensions
    
    integer,parameter::MAXIT_MOMENTUM=40 !< max number of momentum prediction equation iteration
    integer,parameter::MAXIT_PRESSURE=20 !< max number of pressure correction equation iteration
    integer,parameter::MAXIT_DENSITY=20 !< max number of density equation iteration
  