diff --git a/src/hydro/rsolvers/hydro_hlle.hpp b/src/hydro/rsolvers/hydro_hlle.hpp
index 814acc63..922a4d80 100644
--- a/src/hydro/rsolvers/hydro_hlle.hpp
+++ b/src/hydro/rsolvers/hydro_hlle.hpp
@@ -50,7 +50,7 @@ struct Riemann<Fluid::euler, RiemannSolver::hlle> {
     Real gm1 = gamma - 1.0;
     Real igm1 = 1.0 / gm1;
     parthenon::par_for_inner(member, il, iu, [&](const int i) {
-      Real wli[(NHYDRO)], wri[(NHYDRO)];
+      Real wli[(NHYDRO)], wri[(NHYDRO)], wroe[(NHYDRO)];
       Real fl[(NHYDRO)], fr[(NHYDRO)], flxi[(NHYDRO)];
       //--- Step 1.  Load L/R states into local variables
       wli[IDN] = wl(IDN, i);
@@ -65,34 +65,37 @@ struct Riemann<Fluid::euler, RiemannSolver::hlle> {
       wri[IV3] = wr(ivz, i);
       wri[IPR] = wr(IPR, i);
 
-      //--- Step 2.  Compute middle state estimates with PVRS (Toro 10.5.2)
-      Real al, ar, el, er;
-      Real cl = eos.SoundSpeed(wli);
-      Real cr = eos.SoundSpeed(wri);
-      el = wli[IPR] * igm1 +
-           0.5 * wli[IDN] * (SQR(wli[IV1]) + SQR(wli[IV2]) + SQR(wli[IV3]));
-      er = wri[IPR] * igm1 +
-           0.5 * wri[IDN] * (SQR(wri[IV1]) + SQR(wri[IV2]) + SQR(wri[IV3]));
-      Real rhoa = .5 * (wli[IDN] + wri[IDN]); // average density
-      Real ca = .5 * (cl + cr);               // average sound speed
-      Real pmid = .5 * (wli[IPR] + wri[IPR] + (wli[IV1] - wri[IV1]) * rhoa * ca);
-
-      //--- Step 3.  Compute sound speed in L,R
-      Real ql, qr;
-      ql = (pmid <= wli[IPR])
-               ? 1.0
-               : (1.0 + (gamma + 1) / std::sqrt(2 * gamma) * (pmid / wli[IPR] - 1.0));
-      qr = (pmid <= wri[IPR])
-               ? 1.0
-               : (1.0 + (gamma + 1) / std::sqrt(2 * gamma) * (pmid / wri[IPR] - 1.0));
-
-      //--- Step 4. Compute the max/min wave speeds based on L/R states
-
-      al = wli[IV1] - cl * ql;
-      ar = wri[IV1] + cr * qr;
-
-      Real bp = ar > 0.0 ? ar : 0.0;
-      Real bm = al < 0.0 ? al : 0.0;
+      //--- Step 2.  Compute Roe-averaged state
+      Real sqrtdl = std::sqrt(wli[IDN]);
+      Real sqrtdr = std::sqrt(wri[IDN]);
+      Real isdlpdr = 1.0 / (sqrtdl + sqrtdr);
+
+      wroe[IDN] = sqrtdl * sqrtdr;
+      wroe[IV1] = (sqrtdl * wli[IV1] + sqrtdr * wri[IV1]) * isdlpdr;
+      wroe[IV2] = (sqrtdl * wli[IV2] + sqrtdr * wri[IV2]) * isdlpdr;
+      wroe[IV3] = (sqrtdl * wli[IV3] + sqrtdr * wri[IV3]) * isdlpdr;
+
+      // Following Roe(1981), the enthalpy H=(E+P)/d is averaged for adiabatic flows,
+      // rather than E or P directly.  sqrtdl*hl = sqrtdl*(el+pl)/dl = (el+pl)/sqrtdl
+      const Real el = wli[IPR] * igm1 +
+                      0.5 * wli[IDN] * (SQR(wli[IV1]) + SQR(wli[IV2]) + SQR(wli[IV3]));
+      const Real er = wri[IPR] * igm1 +
+                      0.5 * wri[IDN] * (SQR(wri[IV1]) + SQR(wri[IV2]) + SQR(wri[IV3]));
+      const Real hroe = ((el + wli[IPR]) / sqrtdl + (er + wri[IPR]) / sqrtdr) * isdlpdr;
+
+      //--- Step 3.  Compute sound speed in L,R, and Roe-averaged states
+
+      const Real cl = eos.SoundSpeed(wli);
+      const Real cr = eos.SoundSpeed(wri);
+      Real q = hroe - 0.5 * (SQR(wroe[IV1]) + SQR(wroe[IV2]) + SQR(wroe[IV3]));
+      const Real a = (q < 0.0) ? 0.0 : std::sqrt(gm1 * q);
+
+      //--- Step 4. Compute the max/min wave speeds based on L/R and Roe-averaged values
+      const Real al = std::min((wroe[IV1] - a), (wli[IV1] - cl));
+      const Real ar = std::max((wroe[IV1] + a), (wri[IV1] + cr));
+
+      Real bp = ar > 0.0 ? ar : TINY_NUMBER;
+      Real bm = al < 0.0 ? al : TINY_NUMBER;
 
       //-- Step 5. Compute L/R fluxes along lines bm/bp: F_L - (S_L)U_L; F_R - (S_R)U_R
       Real vxl = wli[IV1] - bm;