get_cmax Function

pure real(dp) function get_cmax(u, x, flux_dim) result(wC)
  
  real(dp), intent(in)  :: u(nw_phys)
  real(dp), intent(in)  :: x(1:ndim)
  integer, intent(in)   :: flux_dim

  real(dp) :: tmp1,tmp2,v2,csound2,vidim,cmax,cmin
  real(dp) :: rho,rhoh,pth,E

    ! input u is in primitive form
    ! auxiliaries are filled here
    rho=u(iw_rho)
    rhoh=u(xi_)/u(lfac_)**2
    pth=u(iw_e)
    !!tmp1=u(iw_rho)
    !!tmp2=u(xi_)/u(lfac_)**2.0d0
    !!v2=u(iw_e)
    !! begin: call srhd_get_csound2_prim_eos(tmp1,tmp2,v2,csound2)
    if(srhd_eos) then
       E = (rhoh+dsqrt(rhoh**2+(srhd_gamma**2-1.0_dp)&
            *rho**2))/(srhd_gamma+1.0_dp)
       csound2=(pth*((srhd_gamma+1.0_dp)&
                          +gamma_1*(rho/E)**2))&
                      /(2.0d0*rhoh)
    else
       csound2=srhd_gamma*pth/rhoh
    end if
    !! end call srhd_get_csound2_prim_eos(tmp1,tmp2,v2,csound2)

    v2    = 1.0d0-1.0d0/u(lfac_)**2
    vidim = u(iw_mom(flux_dim))/u(lfac_)
    tmp2  = vidim**2
    tmp1  = 1.0d0-v2*csound2-tmp2*(1.0d0-csound2)
    ! the radicand is positive analytically, as in estimate_speeds_minmax; clip
    ! it so an inadmissible state (c_s^2 > 1 from a collapsed rhoh) yields a
    ! degenerate speed rather than a NaN. get_cmax feeds setdt, where a single
    ! NaN survives the global MPI_ALLREDUCE and takes dt down for every rank.
    tmp2  = dsqrt(max(csound2*(1.0_dp-v2)*tmp1, 0.0_dp))
    tmp1  = vidim*(1.0_dp-csound2)
    cmax  = (tmp1+tmp2)/(1.0_dp-v2*csound2)
    cmin  = (tmp1-tmp2)/(1.0_dp-v2*csound2)
    ! Limit by speed of light
    cmin = max(cmin, - 1.0d0)
    cmin = min(cmin,   1.0d0)
    cmax = max(cmax, - 1.0d0)
    cmax = min(cmax,   1.0d0)
    ! now take extremal value only for dt limit
    wC = max(dabs(cmax),dabs(cmin))

end function get_cmax