diff --git a/Util/exact_riemann/exact_riemann_rarefaction.H b/Source/hydro/exact_riemann_rarefaction.H
similarity index 100%
rename from Util/exact_riemann/exact_riemann_rarefaction.H
rename to Source/hydro/exact_riemann_rarefaction.H
diff --git a/Util/exact_riemann/exact_riemann_sample.H b/Source/hydro/exact_riemann_sample.H
similarity index 95%
rename from Util/exact_riemann/exact_riemann_sample.H
rename to Source/hydro/exact_riemann_sample.H
index c255ea1cce..a147763f6c 100644
--- a/Util/exact_riemann/exact_riemann_sample.H
+++ b/Source/hydro/exact_riemann_sample.H
@@ -14,7 +14,7 @@ riemann_sample(const amrex::Real rho_l, const amrex::Real u_l, const amrex::Real
                const amrex::Real rho_r, const amrex::Real u_r, const amrex::Real p_r, const amrex::Real* xn_r,
                const amrex::Real ustar, const amrex::Real pstar,
                const amrex::Real W_l, const amrex::Real W_r,
-               const amrex::Real x, const amrex::Real xjump, const amrex::Real time,
+               const amrex::Real xi,
                amrex::Real& rho, amrex::Real& u, amrex::Real& p, amrex::Real* xn) {
 
     // get the initial sound speeds
@@ -47,11 +47,6 @@ riemann_sample(const amrex::Real rho_l, const amrex::Real u_l, const amrex::Real
 
     // This follows from the discussion around C&G Eq. 15
 
-    // compute xi = x/t -- this is the similarity variable for the
-    // solution
-
-    amrex::Real xi = (x - xjump) / time;
-
     // check which side of the contact we need to worry about
 
     amrex::Real chi = std::copysign(1.0_rt, xi - ustar);
diff --git a/Util/exact_riemann/exact_riemann_shock.H b/Source/hydro/exact_riemann_shock.H
similarity index 100%
rename from Util/exact_riemann/exact_riemann_shock.H
rename to Source/hydro/exact_riemann_shock.H
diff --git a/Util/exact_riemann/exact_riemann_star_state.H b/Source/hydro/exact_riemann_star_state.H
similarity index 100%
rename from Util/exact_riemann/exact_riemann_star_state.H
rename to Source/hydro/exact_riemann_star_state.H
diff --git a/Source/hydro/riemann.cpp b/Source/hydro/riemann.cpp
index dc6eb0df98..2c2f284b80 100644
--- a/Source/hydro/riemann.cpp
+++ b/Source/hydro/riemann.cpp
@@ -71,7 +71,9 @@ Castro::cmpflx_plus_godunov(const Box& bx,
     {
 
 
-        if (riemann_solver == 0 || riemann_solver == 1) {
+        if (riemann_solver == riemann_constants::TWO_SHOCK_CGF ||
+            riemann_solver == riemann_constants::TWO_SHOCK_CG ||
+            riemann_solver == riemann_constants::EXACT) {
             // approximate state Riemann solvers
 
             // first find the interface state on the current interface
@@ -120,7 +122,7 @@ Castro::cmpflx_plus_godunov(const Box& bx,
                 }
             }
 
-        } else if (riemann_solver == 2) {
+        } else if (riemann_solver == riemann_constants::HLLC) {
             // HLLC
             HLL::HLLC(i, j, k, idir,
                       qm, qp,
diff --git a/Source/hydro/riemann_constants.H b/Source/hydro/riemann_constants.H
index b29b92e2a5..eec8a6e501 100644
--- a/Source/hydro/riemann_constants.H
+++ b/Source/hydro/riemann_constants.H
@@ -6,6 +6,14 @@
 using namespace amrex::literals;
 
 namespace riemann_constants {
+
+    enum RiemannSolver : std::uint8_t {
+        TWO_SHOCK_CGF=0,
+        TWO_SHOCK_CG,
+        HLLC,
+        EXACT
+    };
+
     constexpr amrex::Real smlp1 = 1.e-10_rt;
     constexpr amrex::Real small = 1.e-8_rt;
     constexpr amrex::Real smallu = 1.e-12_rt;
diff --git a/Source/hydro/riemann_solvers.H b/Source/hydro/riemann_solvers.H
index 19b93287a7..b547dc11ef 100644
--- a/Source/hydro/riemann_solvers.H
+++ b/Source/hydro/riemann_solvers.H
@@ -4,6 +4,8 @@
 #include <Castro_util.H>
 #include <riemann_type.H>
 #include <riemann_2shock_solvers.H>
+#include <exact_riemann_star_state.H>
+#include <exact_riemann_sample.H>
 #ifdef HYBRID_MOMENTUM
 #include <hybrid.H>
 #endif
@@ -318,18 +320,87 @@ riemann_state(const int i, const int j, const int k, const int idir,
 
 
   // Solve Riemann problem
-  if (riemann_solver == 0) {
+  if (riemann_solver == riemann_constants::TWO_SHOCK_CGF) {
       // Colella, Glaz, & Ferguson solver
 
       TwoShock::riemannus(ql, qr, raux, qint);
 
-  } else if (riemann_solver == 1) {
+  } else if (riemann_solver == riemann_constants::TWO_SHOCK_CG) {
       // Colella & Glaz solver
 
 #ifndef RADIATION
       TwoShock::riemanncg(ql, qr, raux, qint);
 #endif
 
+  } else if (riemann_solver == riemann_constants::EXACT) {
+      // Exact Riemann solver
+
+      // solve for the star state
+
+      amrex::Real ustar{};
+      amrex::Real pstar{};
+      amrex::Real W_l{};
+      amrex::Real W_r{};
+
+      amrex::Real xn_l[NumSpec];
+      for (int n = 0; n < NumSpec; ++n) {
+          xn_l[n] = qm(i,j,k,QFS+n);
+      }
+      amrex::Real xn_r[NumSpec];
+      for (int n = 0; n < NumSpec; ++n) {
+          xn_r[n] = qp(i,j,k,QFS+n);
+      }
+
+      riemann_star_state(ql.rho, ql.un, ql.p, xn_l,
+                         qr.rho, qr.un, qr.p, xn_r,
+                         ustar, pstar, W_l, W_r);
+
+      // sample the solution at xi = 0
+
+      amrex::Real xi{};
+
+      amrex::Real rho{};
+      amrex::Real u{};
+      amrex::Real p{};
+      amrex::Real xn_s[NumSpec]{};
+
+      riemann_sample(ql.rho, ql.un, ql.p, xn_l,
+                     qr.rho, qr.un, qr.p, xn_r,
+                     ustar, pstar, W_l, W_r,
+                     xi,
+                     rho, u, p, xn_s);
+
+      // pack the RiemannState interface state
+
+      qint.rho = rho;
+      qint.un = u;
+      qint.p = p;
+
+      // get the energy
+      eos_rep_t eos_state;
+      eos_state.rho = qint.rho;
+      eos_state.T = castro::T_guess;
+      eos_state.p = qint.p;
+      for (int n = 0; n < NumSpec; ++n) {
+          eos_state.xn[n] = xn_s[n];
+      }
+
+      eos(eos_input_rp, eos_state);
+
+      qint.rhoe = qint.rho * eos_state.e;
+      qint.gamc = eos_state.gam1;
+
+      // we need to manually do the transverse sampling based on the
+      // contact
+
+      if (ustar > 0) {
+          qint.ut = ql.ut;
+          qint.utt = ql.utt;
+      } else {
+          qint.ut = qr.ut;
+          qint.utt = qr.utt;
+      }
+
 #ifndef AMREX_USE_GPU
   } else {
       amrex::Error("ERROR: invalid value of riemann_solver");
diff --git a/Util/exact_riemann/README.md b/Util/exact_riemann/README.md
index 013a81693a..16fca8674b 100644
--- a/Util/exact_riemann/README.md
+++ b/Util/exact_riemann/README.md
@@ -1,9 +1,12 @@
-This is an exact Riemann solver for a general equation of state.  It
-follows the outline for Colella & Glaz 1985, section 1.  This is too
-slow to use in an actual hydro run, but instead is intended to
-generate exact solutions to the Riemann problem for comparison with
-Castro shocktube output.  Several inputs files for Helmholtz EOS-based
-shocktubes are provided.
+# Exact Riemann solver
+
+This is a driver for the exact Riemann solver for a general equation
+of state.  The main implementation is in Source/hydro.
+
+The exact Riemann solver follows the outline for Colella & Glaz 1985,
+section 1 and this driver is intended to generate exact solutions to
+the Riemann problem for comparison with Castro shocktube output.
+Several inputs files for Helmholtz EOS-based shocktubes are provided.
 
 This solver is used in Zingale & Katz (2015):
 
diff --git a/Util/exact_riemann/exact_riemann.H b/Util/exact_riemann/exact_riemann.H
index b1cf14e832..93e9140e52 100644
--- a/Util/exact_riemann/exact_riemann.H
+++ b/Util/exact_riemann/exact_riemann.H
@@ -78,10 +78,15 @@ exact_riemann() {
 
         amrex::Real rho, u, p, xn_s[NumSpec];
 
+        // compute xi = x/t -- this is the similarity variable for the
+        // solution
+
+        amrex::Real xi = (x - problem::xjump) / problem::t;
+
         riemann_sample(problem::rho_l, problem::u_l, problem::p_l, xn,
                        problem::rho_r, problem::u_r, problem::p_r, xn,
                        ustar, pstar, W_l, W_r,
-                       x, problem::xjump, problem::t,
+                       xi,
                        rho, u, p, xn_s);
 
         // get the thermodynamics for this state for output