!:physfluxRequestFluxes:!
        CASE (RequestFluxes)
          IF (MINVAL(q1D(1-mbc:mx+mbc,1))<=zero) THEN
            PRINT *,'Error: Negative or zero density in physflux.'
            STOP
          END IF
          ! Compute cell centered quantities required in Riemann solve & entropy fix
          rho(iCell) = 1./q1D(iCell,1)
          u(iCell) = q1D(iCell,mu)*rho(iCell)
          v(iCell) = q1D(iCell,mv)*rho(iCell)
          rho(iCell) = q1D(iCell,1)
          sqrtrho(iCell) = SQRT(rho(iCell))
          pres(iCell) = gamma1*( q1D(iCell,4) - &
                                 half*(q1D(iCell,mu)*u(iCell) + q1D(iCell,mv)*v(iCell) ) )
          c = gamma*pres/rho
          IF (MINVAL(c)<=zero) THEN
            PRINT *,'Error: Non-physical centered sound velocity in physflux, RequestFluxes.'
            STOP
          END IF
          c = SQRT(c)
          h(iCell) = (q1D(iCell,4) + pres(iCell)) / rho(iCell)
          ! Compute the physical fluxes at cell centers
          f(iCell,1)  = q1D(iCell,mu)
          f(iCell,mu) = f(iCell,1)*u(iCell)+pres(iCell)
          f(iCell,mv) = f(iCell,1)*v(iCell)
          f(iCell,4)  = f(iCell,1)*h(iCell)