diff --git a/Diagnostics/DustCollapse/GNUmakefile b/Diagnostics/DustCollapse/GNUmakefile
index 8facc39491..5c543008bf 100644
--- a/Diagnostics/DustCollapse/GNUmakefile
+++ b/Diagnostics/DustCollapse/GNUmakefile
@@ -10,8 +10,6 @@ USE_OMP = FALSE
 
 USE_REACT = FALSE
 
-USE_ACC = FALSE
-
 # programs to be compiled
 ALL: dustcollapse_$(DIM)d.ex
 
diff --git a/Diagnostics/Radiation/GNUmakefile b/Diagnostics/Radiation/GNUmakefile
index 2fd07f63a6..04449ffb4b 100644
--- a/Diagnostics/Radiation/GNUmakefile
+++ b/Diagnostics/Radiation/GNUmakefile
@@ -10,8 +10,6 @@ USE_OMP = FALSE
 
 USE_REACT = FALSE
 
-USE_ACC = FALSE
-
 ALL: radhelp
 
 radhelp:
diff --git a/Diagnostics/Sedov/GNUmakefile b/Diagnostics/Sedov/GNUmakefile
index d3a5bf5095..1773875f5e 100644
--- a/Diagnostics/Sedov/GNUmakefile
+++ b/Diagnostics/Sedov/GNUmakefile
@@ -10,8 +10,6 @@ USE_OMP = FALSE
 
 USE_REACT = FALSE
 
-USE_ACC = FALSE
-
 # programs to be compiled
 ALL: sedov_$(DIM)d.ex
 
diff --git a/Exec/Make.Castro b/Exec/Make.Castro
index 342f5fe97d..910a577563 100644
--- a/Exec/Make.Castro
+++ b/Exec/Make.Castro
@@ -150,12 +150,6 @@ build_status:
 # The default is to include the sponge functionality
 DEFINES += -DSPONGE
 
-# OpenACC support
-ifeq ($(USE_ACC), TRUE)
-  DEFINES += -DACC
-endif
-
-
 
 #------------------------------------------------------------------------------
 # Castro directories
diff --git a/Exec/science/Detonation/problem_initialize_state_data.H b/Exec/science/Detonation/problem_initialize_state_data.H
index f594e8f15e..8460f583d6 100644
--- a/Exec/science/Detonation/problem_initialize_state_data.H
+++ b/Exec/science/Detonation/problem_initialize_state_data.H
@@ -22,7 +22,15 @@ void problem_initialize_state_data (int i, int j, int k,
 
     state(i,j,k,URHO) = problem::dens;
 
-    Real sigma = 1.0_rt / (1.0_rt + std::exp(-(c_T - xcen)/ width));
+    Real sigma_arg = -(c_T - xcen) / width;
+    Real sigma;
+    // need to avoid FP overflow for sigma_arg >= 709
+    // 1/(1 + exp(100)) ~= 3e-44, which is much smaller than machine epsilon
+    if (sigma_arg < 100) {
+        sigma = 1.0_rt / (1.0_rt + std::exp(sigma_arg));
+    } else {
+        sigma = 0.0_rt;
+    }
 
     state(i,j,k,UTEMP) = problem::T_l + (problem::T_r - problem::T_l) * (1.0_rt - sigma);
 
diff --git a/Exec/unit_tests/model_burner/GNUmakefile b/Exec/unit_tests/model_burner/GNUmakefile
index d159960fe0..94f3300a5b 100644
--- a/Exec/unit_tests/model_burner/GNUmakefile
+++ b/Exec/unit_tests/model_burner/GNUmakefile
@@ -7,7 +7,6 @@ COMP	         = gnu
 
 USE_MPI          = TRUE
 USE_OMP          = FALSE
-USE_ACC          = FALSE
 
 USE_DIFFUSION    = FALSE
 USE_GRAV         = FALSE
diff --git a/Source/driver/Castro.H b/Source/driver/Castro.H
index 30976c248d..8c7b9c13f8 100644
--- a/Source/driver/Castro.H
+++ b/Source/driver/Castro.H
@@ -16,10 +16,6 @@ using namespace castro;
 
 #include <fundamental_constants.H>
 
-constexpr int HISTORY_SIZE=40;
-constexpr int PSTAR_BISECT_FACTOR = 5;
-
-
 #include <AMReX_BC_TYPES.H>
 #include <AMReX_AmrLevel.H>
 #include <AMReX_iMultiFab.H>
diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp
index f478f8402d..e79cf0574e 100644
--- a/Source/driver/Castro.cpp
+++ b/Source/driver/Castro.cpp
@@ -58,6 +58,8 @@
 #include <ambient.H>
 #include <castro_limits.H>
 
+#include <riemann_constants.H>
+
 using namespace amrex;
 
 bool         Castro::signalStopJob = false;
@@ -392,7 +394,7 @@ Castro::read_params ()
 #endif
 
     if (riemann_solver == 1) {
-        if (riemann_shock_maxiter > HISTORY_SIZE) {
+        if (riemann_shock_maxiter > riemann_constants::HISTORY_SIZE) {
             amrex::Error("riemann_shock_maxiter > HISTORY_SIZE");
         }
 
diff --git a/Source/hydro/Make.package b/Source/hydro/Make.package
index 6b6ba6c7e8..c5813e27c3 100644
--- a/Source/hydro/Make.package
+++ b/Source/hydro/Make.package
@@ -30,6 +30,7 @@ CEXE_sources += riemann.cpp
 CEXE_headers += HLL_solvers.H
 CEXE_headers += riemann_solvers.H
 CEXE_headers += riemann_2shock_solvers.H
+CEXE_headers += riemann_constants.H
 CEXE_headers += riemann_type.H
 CEXE_headers += slope.H
 CEXE_headers += reconstruction.H
diff --git a/Source/hydro/riemann_2shock_solvers.H b/Source/hydro/riemann_2shock_solvers.H
index f81b495f29..6a70d0e117 100644
--- a/Source/hydro/riemann_2shock_solvers.H
+++ b/Source/hydro/riemann_2shock_solvers.H
@@ -5,6 +5,7 @@
 
 using namespace amrex::literals;
 
+#include <riemann_constants.H>
 #include <riemann_type.H>
 #ifdef RADIATION
 #include <Radiation.H>
@@ -60,7 +61,7 @@ namespace TwoShock {
                     const amrex::Real gdot, const amrex::Real gmin, const amrex::Real gmax,
                     const int lriemann_shock_maxiter, const amrex::Real lriemann_pstar_tol,
                     amrex::Real& pstar, amrex::Real& gamstar,
-                    bool& converged, amrex::GpuArray<amrex::Real, PSTAR_BISECT_FACTOR*HISTORY_SIZE>& pstar_hist_extra) {
+                    bool& converged, amrex::GpuArray<amrex::Real, riemann_constants::PSTAR_BISECT_FACTOR * riemann_constants::HISTORY_SIZE>& pstar_hist_extra) {
 
       // we want to zero
       // f(p*) = u*_l(p*) - u*_r(p*)
@@ -106,7 +107,7 @@ namespace TwoShock {
       converged = false;
       amrex::Real pstar_c = 0.0;
 
-      for (int iter = 0; iter < PSTAR_BISECT_FACTOR * lriemann_shock_maxiter; iter++) {
+      for (int iter = 0; iter < riemann_constants::PSTAR_BISECT_FACTOR * lriemann_shock_maxiter; iter++) {
 
         pstar_c = 0.5_rt * (pstar_lo + pstar_hi);
         pstar_hist_extra[iter] = pstar_c;
@@ -167,7 +168,7 @@ namespace TwoShock {
       constexpr amrex::Real weakwv = 1.e-3_rt;
 
     #ifndef AMREX_USE_GPU
-      amrex::GpuArray<amrex::Real, HISTORY_SIZE> pstar_hist;
+      amrex::GpuArray<amrex::Real, riemann_constants::HISTORY_SIZE> pstar_hist;
     #endif
 
 
@@ -336,7 +337,7 @@ namespace TwoShock {
               pstarl = amrex::max(pstarl, small_pres);
               pstaru = amrex::max(pstaru, small_pres);
 
-              amrex::GpuArray<amrex::Real, PSTAR_BISECT_FACTOR*HISTORY_SIZE> pstar_hist_extra;
+              amrex::GpuArray<amrex::Real, riemann_constants::PSTAR_BISECT_FACTOR * riemann_constants::HISTORY_SIZE> pstar_hist_extra;
 
               pstar_bisection(pstarl, pstaru,
                               ql.un, ql.p, taul, gamel, clsql,
@@ -352,7 +353,7 @@ namespace TwoShock {
                       std::cout << iter_l << " " << pstar_hist[iter_l] << std::endl;
                   }
                   std::cout << "pstar extra history: " << std::endl;
-                  for (int iter_l = 0; iter_l < PSTAR_BISECT_FACTOR * riemann_shock_maxiter; iter_l++) {
+                  for (int iter_l = 0; iter_l < riemann_constants::PSTAR_BISECT_FACTOR * riemann_shock_maxiter; iter_l++) {
                       std::cout << iter_l << " " << pstar_hist_extra[iter_l] << std::endl;
                   }
 
diff --git a/Source/hydro/riemann_constants.H b/Source/hydro/riemann_constants.H
new file mode 100644
index 0000000000..f06776238f
--- /dev/null
+++ b/Source/hydro/riemann_constants.H
@@ -0,0 +1,16 @@
+#ifndef RIEMANN_CONSTANTS_H
+#define RIEMANN_CONSTANTS_H
+
+#include <AMReX_REAL.H>
+
+using namespace amrex::literals;
+
+namespace riemann_constants {
+    constexpr amrex::Real smlp1 = 1.e-10_rt;
+    constexpr amrex::Real small = 1.e-8_rt;
+    constexpr amrex::Real smallu = 1.e-12_rt;
+    constexpr int HISTORY_SIZE=40;
+    constexpr int PSTAR_BISECT_FACTOR = 5;
+}
+
+#endif
diff --git a/Source/hydro/riemann_type.H b/Source/hydro/riemann_type.H
index c8be48f47a..658befe151 100644
--- a/Source/hydro/riemann_type.H
+++ b/Source/hydro/riemann_type.H
@@ -5,13 +5,6 @@
 
 using namespace amrex::literals;
 
-namespace riemann_constants {
-    const amrex::Real smlp1 = 1.e-10_rt;
-    const amrex::Real small = 1.e-8_rt;
-    const amrex::Real smallu = 1.e-12_rt;
-}
-
-
 struct RiemannState
 {
     amrex::Real rho;
diff --git a/Source/sdc/sdc_newton_solve.H b/Source/sdc/sdc_newton_solve.H
index 580b3e6c43..7a796e270c 100644
--- a/Source/sdc/sdc_newton_solve.H
+++ b/Source/sdc/sdc_newton_solve.H
@@ -159,8 +159,7 @@ sdc_newton_solve(const Real dt_m,
 
         // solve the linear system: Jac dU = -f
 #ifdef NEW_NETWORK_IMPLEMENTATION
-        RHS::dgefa(Jac);
-        info = 0;
+        info = RHS::dgefa(Jac);
 #else
         IArray1D ipvt;
         dgefa<NumSpec+1>(Jac, ipvt, info);
diff --git a/Util/exact_riemann/GNUmakefile b/Util/exact_riemann/GNUmakefile
index 8be2709e1a..fbf934926a 100644
--- a/Util/exact_riemann/GNUmakefile
+++ b/Util/exact_riemann/GNUmakefile
@@ -33,6 +33,6 @@ Blocs   := .
 
 # we explicitly want runtime_parameters.H so we have access to
 # Castro's runtime parameter system
-Blocs += $(CASTRO_HOME)/Source/driver
+Blocs += $(CASTRO_HOME)/Source/driver $(CASTRO_HOME)/Source/hydro
 
 include $(CASTRO_HOME)/Exec/Make.Castro
diff --git a/Util/exact_riemann/README.md b/Util/exact_riemann/README.md
index ada8e281ec..013a81693a 100644
--- a/Util/exact_riemann/README.md
+++ b/Util/exact_riemann/README.md
@@ -12,3 +12,14 @@ https://ui.adsabs.harvard.edu/abs/2015ApJS..216...31Z/abstract
 and more details are given there.
 
 To build the solver, simply type 'make' in this directory.
+
+
+Notes:
+
+* ``test2`` is a vacuum
+
+  The velocity in the star region should be 0, but it seems that
+  slight differences in the rarefaction integration from the left and
+  from the right give different roundoff, resulting in a small, finite
+  velocity.  This can be controlled a bit by tightening the tolerance
+  in the rarefaction integration.
diff --git a/Util/exact_riemann/main.cpp b/Util/exact_riemann/main.cpp
index 7f75f6fc36..4707027763 100644
--- a/Util/exact_riemann/main.cpp
+++ b/Util/exact_riemann/main.cpp
@@ -41,9 +41,7 @@ int main(int argc, char *argv[]) {
     network_init();
 #endif
 
-    amrex::Real small_temp = 1.e-200;
-    amrex::Real small_dens = 1.e-200;
-    eos_init(small_temp, small_dens);
+    eos_init(castro::small_temp, castro::small_dens);
 
     exact_riemann();
 
diff --git a/Util/exact_riemann/riemann_shock.H b/Util/exact_riemann/riemann_shock.H
index 3a74b6e3d3..a6bb04eb9e 100644
--- a/Util/exact_riemann/riemann_shock.H
+++ b/Util/exact_riemann/riemann_shock.H
@@ -6,6 +6,7 @@
 #include <eos_type.H>
 #include <eos.H>
 
+#include <riemann_constants.H>
 
 AMREX_INLINE
 void
@@ -62,7 +63,7 @@ newton_shock(amrex::Real& W_s, const amrex::Real pstar,
     converged = false;
 
     int iter = 1;
-    while (! converged && iter < riemann_solver::max_iters) {
+    while (! converged && iter < castro::riemann_shock_maxiter) {
 
         amrex::Real rhostar_s, f, fprime;
 
@@ -71,7 +72,7 @@ newton_shock(amrex::Real& W_s, const amrex::Real pstar,
 
         amrex::Real dW = -f / fprime;
 
-        if (std::abs(dW) < riemann_solver::tol * W_s) {
+        if (std::abs(dW) < castro::riemann_pstar_tol * W_s) {
             converged = true;
         }
 
@@ -96,7 +97,7 @@ shock(const amrex::Real pstar,
 
     amrex::ignore_unused(u_s);
 
-    amrex::Real rhostar_hist[riemann_solver::max_iters], Ws_hist[riemann_solver::max_iters];
+    amrex::Real rhostar_hist[riemann_constants::HISTORY_SIZE], Ws_hist[riemann_constants::HISTORY_SIZE];
 
     // compute the Z_s function for a shock following C&G Eq. 20 and
     // 23.  Here the "_s" variables are the state (L or R) that we are
@@ -147,7 +148,7 @@ shock(const amrex::Real pstar,
     amrex::Real rhostar_s;
 
     if (taustar_s < 0.0_rt) {
-        rhostar_s = riemann_solver::smallrho;
+        rhostar_s = riemann_constants::small;
         W_s = std::sqrt((pstar - p_s) / (1.0_rt / rho_s - 1.0_rt / rhostar_s));
     }
 
@@ -162,7 +163,7 @@ shock(const amrex::Real pstar,
     // now did we converge?
 
     if (! converged) {
-        for (int i = 0; i < riemann_solver::max_iters; ++i) {
+        for (int i = 0; i < castro::riemann_shock_maxiter; ++i) {
             std::cout << i << " " << rhostar_hist[i] << " " << Ws_hist[i] << std::endl;
         }
 
diff --git a/Util/exact_riemann/riemann_star_state.H b/Util/exact_riemann/riemann_star_state.H
index 2448247129..31e1952722 100644
--- a/Util/exact_riemann/riemann_star_state.H
+++ b/Util/exact_riemann/riemann_star_state.H
@@ -142,12 +142,11 @@ riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex::
 
         amrex::Real pstar_new = pstar - Z_l * Z_r * (ustar_r - ustar_l) / (Z_l + Z_r);
 
-        // estimate the error in the current star solution
-        amrex::Real err1 = std::abs(ustar_r - ustar_l);
-        amrex::Real err2 = pstar_new - pstar;
+        // estimate the error in the current pstar solution
 
-        if (err1 < tol * amrex::max(std::abs(ustar_l), std::abs(ustar_r)) &&
-            err2 < tol * pstar) {
+        amrex::Real error = std::abs(pstar_new - pstar);
+
+        if (error < tol * pstar) {
             converged = true;
         }