diff --git a/EOS/breakout/actual_eos.H b/EOS/breakout/actual_eos.H
index cae836ae0b..0674931872 100644
--- a/EOS/breakout/actual_eos.H
+++ b/EOS/breakout/actual_eos.H
@@ -50,7 +50,7 @@ void actual_eos (I input, T& state)
 {
     static_assert(std::is_same_v<I, eos_input_t>, "input must be an eos_input_t");
 
-    const Real R = C::k_B * C::n_A;
+    const amrex::Real R = C::k_B * C::n_A;
 
     // Calculate mu. This is the only difference between
     // this EOS and gamma_law.
@@ -62,8 +62,8 @@ void actual_eos (I input, T& state)
     {
 
         // dens, temp and xmass are inputs
-        Real cv = R / (state.mu * (gamma_const - 1.0_rt));
-        Real e = cv * state.T;
+        amrex::Real cv = R / (state.mu * (gamma_const - 1.0_rt));
+        amrex::Real e = cv * state.T;
         if constexpr (has_energy<T>::value) {
             state.cv = cv;
             state.e = e;
@@ -100,7 +100,7 @@ void actual_eos (I input, T& state)
         // dens, pres, and xmass are inputs
 
         if constexpr (has_pressure<T>::value) {
-            Real poverrho = state.p / state.rho;
+            amrex::Real poverrho = state.p / state.rho;
             state.T = poverrho * state.mu * (1.0_rt / R);
             if constexpr (has_energy<T>::value) {
                 state.e = poverrho * (1.0_rt / (gamma_const - 1.0_rt));
@@ -115,7 +115,7 @@ void actual_eos (I input, T& state)
         // dens, energy, and xmass are inputs
 
         if constexpr (has_energy<T>::value) {
-            Real poverrho = (gamma_const - 1.0_rt) * state.e;
+            amrex::Real poverrho = (gamma_const - 1.0_rt) * state.e;
             state.T = poverrho * state.mu * (1.0_rt / R);
 
             if constexpr (has_pressure<T>::value) {
diff --git a/EOS/eos_composition.H b/EOS/eos_composition.H
index be39fa0492..99b91df977 100644
--- a/EOS/eos_composition.H
+++ b/EOS/eos_composition.H
@@ -9,16 +9,16 @@
 #include <actual_network.H>
 #endif
 
-using namespace amrex;
+using namespace amrex::literals;
 
 #ifdef AUX_THERMO
 using namespace AuxZero;
 #endif
 
 struct eos_xderivs_t {
-  Real dedX[NumSpec];
-  Real dpdX[NumSpec];
-  Real dhdX[NumSpec];
+  amrex::Real dedX[NumSpec];
+  amrex::Real dpdX[NumSpec];
+  amrex::Real dhdX[NumSpec];
 };
 
 #ifdef AUX_THERMO
@@ -65,7 +65,7 @@ void composition (T& state)
 
 #else
 
-  Real sum = 0;
+  amrex::Real sum = 0;
   for (int n = 0; n < NumSpec; n++) {
     sum += state.xn[n] * zion[n] * aion_inv[n];
   }
diff --git a/EOS/gamma_law/actual_eos.H b/EOS/gamma_law/actual_eos.H
index dc075b67c9..f7a56ee622 100644
--- a/EOS/gamma_law/actual_eos.H
+++ b/EOS/gamma_law/actual_eos.H
@@ -53,13 +53,13 @@ void actual_eos (I input, T& state)
     static_assert(std::is_same_v<I, eos_input_t>, "input must be an eos_input_t");
 
     // Get the mass of a nucleon from m_u.
-    const Real m_nucleon = C::m_u;
+    const amrex::Real m_nucleon = C::m_u;
 
     if constexpr (has_xn<T>::value) {
         if (eos_assume_neutral) {
             state.mu = state.abar;
         } else {
-            Real sum = 0.0;
+            amrex::Real sum = 0.0;
             for (int n = 0; n < NumSpec; n++) {
                 sum += (zion[n] + 1.0) * state.xn[n] * aion_inv[n];
             }
@@ -193,13 +193,13 @@ void actual_eos (I input, T& state)
     // Now we have the density and temperature (and mass fractions /
     // mu), regardless of the inputs.
 
-    Real Tinv = 1.0 / state.T;
-    Real rhoinv = 1.0 / state.rho;
+    amrex::Real Tinv = 1.0 / state.T;
+    amrex::Real rhoinv = 1.0 / state.rho;
 
     // Compute the pressure simply from the ideal gas law, and the
     // specific internal energy using the gamma-law EOS relation.
-    Real pressure = state.rho * state.T * C::k_B / (state.mu * m_nucleon);
-    Real energy = pressure / (eos_gamma - 1.0) * rhoinv;
+    amrex::Real pressure = state.rho * state.T * C::k_B / (state.mu * m_nucleon);
+    amrex::Real energy = pressure / (eos_gamma - 1.0) * rhoinv;
     if constexpr (has_pressure<T>::value) {
         state.p = pressure;
     }
@@ -215,7 +215,7 @@ void actual_eos (I input, T& state)
     // entropy (per gram) of an ideal monoatomic gas (the Sackur-Tetrode equation)
     // NOTE: this expression is only valid for gamma = 5/3.
     if constexpr (has_entropy<T>::value) {
-        const Real fac = 1.0 / std::pow(2.0 * M_PI * C::hbar * C::hbar, 1.5);
+        const amrex::Real fac = 1.0 / std::pow(2.0 * M_PI * C::hbar * C::hbar, 1.5);
 
         state.s = (C::k_B / (state.mu * m_nucleon)) *
                   (2.5 + std::log((std::pow(state.mu * m_nucleon, 2.5) * rhoinv) *
diff --git a/EOS/helmholtz/actual_eos.H b/EOS/helmholtz/actual_eos.H
index 205e1a6097..dfbe9aa7e7 100644
--- a/EOS/helmholtz/actual_eos.H
+++ b/EOS/helmholtz/actual_eos.H
@@ -43,63 +43,63 @@ const std::string eos_name = "helmholtz";
 // quintic hermite polynomial functions
 // psi0 and its derivatives
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-Real psi0 (Real z)
+amrex::Real psi0 (amrex::Real z)
 {
     return z * z * z * (z * (-6.0e0_rt * z + 15.0e0_rt) -10.0e0_rt) + 1.0e0_rt;
 }
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-Real dpsi0 (Real z)
+amrex::Real dpsi0 (amrex::Real z)
 {
     return z * z * (z * (-30.0e0_rt*z + 60.0e0_rt) - 30.0e0_rt);
 }
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-Real ddpsi0 (Real z)
+amrex::Real ddpsi0 (amrex::Real z)
 {
     return z * (z * (-120.0e0_rt * z + 180.0e0_rt) - 60.0e0_rt);
 }
 
 // psi1 and its derivatives
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-Real psi1 (Real z)
+amrex::Real psi1 (amrex::Real z)
 {
     return z * ( z * z * (z * (-3.0e0_rt * z + 8.0e0_rt) - 6.0e0_rt) + 1.0e0_rt);
 }
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-Real dpsi1 (Real z)
+amrex::Real dpsi1 (amrex::Real z)
 {
     return z * z * ( z * (-15.0e0_rt * z + 32.0e0_rt) - 18.0e0_rt) + 1.0e0_rt;
 }
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-Real ddpsi1 (Real z)
+amrex::Real ddpsi1 (amrex::Real z)
 {
     return z * (z * (-60.0e0_rt * z + 96.0e0_rt) - 36.0e0_rt);
 }
 
 // psi2 and its derivatives
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-Real psi2 (Real z)
+amrex::Real psi2 (amrex::Real z)
 {
     return 0.5e0_rt * z * z * ( z * ( z * (-z + 3.0e0_rt) - 3.0e0_rt) + 1.0e0_rt);
 }
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-Real dpsi2 (Real z)
+amrex::Real dpsi2 (amrex::Real z)
 {
     return 0.5e0_rt * z * (z * (z * (-5.0e0_rt * z + 12.0e0_rt) - 9.0e0_rt) + 2.0e0_rt);
 }
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-Real ddpsi2 (Real z)
+amrex::Real ddpsi2 (amrex::Real z)
 {
     return 0.5e0_rt * (z * ( z * (-20.0e0_rt * z + 36.0e0_rt) - 18.0e0_rt) + 2.0e0_rt);
 }
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-void fwt (const Real* fi, const Real* wt, Real* fwtr)
+void fwt (const amrex::Real* fi, const amrex::Real* wt, amrex::Real* fwtr)
 {
     fwtr[0] = fi[ 0]*wt[0] + fi[ 1]*wt[1] + fi[ 2]*wt[2] + fi[18]*wt[3] + fi[19]*wt[4] + fi[20]*wt[5];
     fwtr[1] = fi[ 3]*wt[0] + fi[ 5]*wt[1] + fi[ 7]*wt[2] + fi[21]*wt[3] + fi[23]*wt[4] + fi[25]*wt[5];
@@ -112,26 +112,26 @@ void fwt (const Real* fi, const Real* wt, Real* fwtr)
 // cubic hermite polynomial functions
 // psi0 & derivatives
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-Real xpsi0(Real z)
+amrex::Real xpsi0(amrex::Real z)
 {
     return z * z * (2.0e0_rt * z - 3.0e0_rt) + 1.0_rt;
 }
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-Real xdpsi0(Real z)
+amrex::Real xdpsi0(amrex::Real z)
 {
     return z * (6.0e0_rt * z - 6.0e0_rt);
 }
 
 // psi1 & derivatives
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-Real xpsi1(Real z)
+amrex::Real xpsi1(amrex::Real z)
 {
     return z * (z * (z - 2.0e0_rt) + 1.0e0_rt);
 }
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-Real xdpsi1(Real z)
+amrex::Real xdpsi1(amrex::Real z)
 {
     return z * (3.0e0_rt * z - 4.0e0_rt) + 1.0e0_rt;
 }
@@ -145,13 +145,13 @@ void apply_electrons (T& state)
     using namespace helmholtz;
 
     // assume complete ionization
-    [[maybe_unused]] Real ytot1 = 1.0e0_rt / state.abar;
+    [[maybe_unused]] amrex::Real ytot1 = 1.0e0_rt / state.abar;
 
     // define mu -- the total mean molecular weight including both electrons and ions
     state.mu = 1.0e0_rt / (1.0e0_rt / state.abar + 1.0e0_rt / state.mu_e);
 
     // enter the table with ye*den
-    Real din = state.y_e * state.rho;
+    amrex::Real din = state.y_e * state.rho;
 
     // hash locate this temperature and density
     int jat = int((std::log10(state.T) - tlo) * tstpi) + 1;
@@ -159,7 +159,7 @@ void apply_electrons (T& state)
     int iat = int((std::log10(din) - dlo) * dstpi) + 1;
     iat = amrex::max(1, amrex::min(iat, imax-1)) - 1;
 
-    Real fi[36];
+    amrex::Real fi[36];
 
     // access the table locations only once
     for (int i = 0; i < 9; ++i) {
@@ -170,13 +170,13 @@ void apply_electrons (T& state)
     }
 
     // various differences
-    Real xt  = amrex::max((state.T - t[jat]) * dti_sav[jat], 0.0e0_rt);
-    Real xd  = amrex::max((din - d[iat]) * ddi_sav[iat], 0.0e0_rt);
-    Real mxt = 1.0e0_rt - xt;
-    Real mxd = 1.0e0_rt - xd;
+    amrex::Real xt  = amrex::max((state.T - t[jat]) * dti_sav[jat], 0.0e0_rt);
+    amrex::Real xd  = amrex::max((din - d[iat]) * ddi_sav[iat], 0.0e0_rt);
+    amrex::Real mxt = 1.0e0_rt - xt;
+    amrex::Real mxd = 1.0e0_rt - xd;
 
     // the six density and six temperature basis functions
-    Real sit[6];
+    amrex::Real sit[6];
 
     sit[0] = psi0(xt);
     sit[1] = psi1(xt) * dt_sav[jat];
@@ -186,7 +186,7 @@ void apply_electrons (T& state)
     sit[4] = -psi1(mxt) * dt_sav[jat];
     sit[5] =  psi2(mxt) * dt2_sav[jat];
 
-    Real sid[6];
+    amrex::Real sid[6];
 
     sid[0] =  psi0(xd);
     sid[1] =  psi1(xd) * dd_sav[iat];
@@ -197,7 +197,7 @@ void apply_electrons (T& state)
     sid[5] =  psi2(mxd) * dd2_sav[iat];
 
     // derivatives of the weight functions
-    Real dsit[6];
+    amrex::Real dsit[6];
 
     dsit[0] =  dpsi0(xt) * dti_sav[jat];
     dsit[1] =  dpsi1(xt);
@@ -207,7 +207,7 @@ void apply_electrons (T& state)
     dsit[4] =  dpsi1(mxt);
     dsit[5] = -dpsi2(mxt) * dt_sav[jat];
 
-    Real dsid[6];
+    amrex::Real dsid[6];
 
     dsid[0] =  dpsi0(xd) * ddi_sav[iat];
     dsid[1] =  dpsi1(xd);
@@ -218,7 +218,7 @@ void apply_electrons (T& state)
     dsid[5] = -dpsi2(mxd) * dd_sav[iat];
 
     // second derivatives of the weight functions
-    Real ddsit[6];
+    amrex::Real ddsit[6];
 
     ddsit[0] =  ddpsi0(xt) * dt2i_sav[jat];
     ddsit[1] =  ddpsi1(xt) * dti_sav[jat];
@@ -234,12 +234,12 @@ void apply_electrons (T& state)
     // and then compute the result as
     // (table data * temperature terms) * density terms.
 
-    Real fwtr[6];
+    amrex::Real fwtr[6];
 
     fwt(fi, sit, fwtr);
 
-    Real free = 0.e0_rt;
-    Real df_d = 0.e0_rt;
+    amrex::Real free = 0.e0_rt;
+    amrex::Real df_d = 0.e0_rt;
 
     for (int i = 0; i <= 5; ++i) {
         // the free energy
@@ -251,8 +251,8 @@ void apply_electrons (T& state)
 
     fwt(fi, dsit, fwtr);
 
-    Real df_t = 0.e0_rt;
-    Real df_dt = 0.e0_rt;
+    amrex::Real df_t = 0.e0_rt;
+    amrex::Real df_dt = 0.e0_rt;
 
     for (int i = 0; i <= 5; ++i) {
         // derivative with respect to temperature
@@ -264,7 +264,7 @@ void apply_electrons (T& state)
 
     fwt(fi, ddsit, fwtr);
 
-    Real df_tt = 0.e0_rt;
+    amrex::Real df_tt = 0.e0_rt;
     for (int i = 0; i <= 5; ++i) {
         // derivative with respect to temperature**2
         df_tt = df_tt + fwtr[i] * sid[i];
@@ -300,7 +300,7 @@ void apply_electrons (T& state)
 
     // Reuse subexpressions that would go into computing the
     // cubic interpolation.
-    Real wdt[16];
+    amrex::Real wdt[16];
 
     for (int i = 0; i <= 3; ++i) {
         wdt[i     ] = sid[0] * sit[i];
@@ -338,7 +338,7 @@ void apply_electrons (T& state)
     fi[15] = dpdf[jat+1][iat+1][3];
 
     // pressure derivative with density
-    Real dpepdd = 0.0e0_rt;
+    amrex::Real dpepdd = 0.0e0_rt;
     for (int i = 0; i <= 15; ++i) {
         dpepdd = dpepdd + fi[i] * wdt[i];
     }
@@ -366,7 +366,7 @@ void apply_electrons (T& state)
     fi[15] = ef[jat+1][iat+1][3];
 
     // electron chemical potential etaele
-    Real etaele = 0.0e0_rt;
+    amrex::Real etaele = 0.0e0_rt;
     for (int i = 0; i <= 15; ++i) {
         etaele = etaele + fi[i] * wdt[i];
     }
@@ -393,7 +393,7 @@ void apply_electrons (T& state)
     fi[15] = xf[jat+1][iat+1][3];
 
     // electron + positron number densities
-    Real xnefer = 0.0e0_rt;
+    amrex::Real xnefer = 0.0e0_rt;
     for (int i = 0; i <= 15; ++i) {
         xnefer = xnefer + fi[i] * wdt[i];
     }
@@ -404,25 +404,25 @@ void apply_electrons (T& state)
     // floating point limit of the subtraction of two large terms.
     // since state.dpdr doesn't enter the maxwell relations at all, use the
     // bicubic interpolation done above instead of this one
-    Real x       = din * din;
-    Real pele    = x * df_d;
-    Real dpepdt  = x * df_dt;
-    [[maybe_unused]] Real s       = dpepdd/state.y_e - 2.0e0_rt * din * df_d;
-    [[maybe_unused]] Real dpepda  = -ytot1 * (2.0e0_rt * pele + s * din);
-    [[maybe_unused]] Real dpepdz  = state.rho*ytot1*(2.0e0_rt * din * df_d  +  s);
+    amrex::Real x       = din * din;
+    amrex::Real pele    = x * df_d;
+    amrex::Real dpepdt  = x * df_dt;
+    [[maybe_unused]] amrex::Real s       = dpepdd/state.y_e - 2.0e0_rt * din * df_d;
+    [[maybe_unused]] amrex::Real dpepda  = -ytot1 * (2.0e0_rt * pele + s * din);
+    [[maybe_unused]] amrex::Real dpepdz  = state.rho*ytot1*(2.0e0_rt * din * df_d  +  s);
 
     x            = state.y_e * state.y_e;
-    Real sele    = -df_t * state.y_e;
-    Real dsepdt  = -df_tt * state.y_e;
-    Real dsepdd  = -df_dt * x;
-    [[maybe_unused]] Real dsepda  = ytot1 * (state.y_e * df_dt * din - sele);
-    [[maybe_unused]] Real dsepdz  = -ytot1 * (state.y_e * df_dt * state.rho  + df_t);
-
-    Real eele    = state.y_e*free + state.T * sele;
-    Real deepdt  = state.T * dsepdt;
-    Real deepdd  = x * df_d + state.T * dsepdd;
-    [[maybe_unused]] Real deepda  = -state.y_e * ytot1 * (free +  df_d * din) + state.T * dsepda;
-    [[maybe_unused]] Real deepdz  = ytot1* (free + state.y_e * df_d * state.rho) + state.T * dsepdz;
+    amrex::Real sele    = -df_t * state.y_e;
+    amrex::Real dsepdt  = -df_tt * state.y_e;
+    amrex::Real dsepdd  = -df_dt * x;
+    [[maybe_unused]] amrex::Real dsepda  = ytot1 * (state.y_e * df_dt * din - sele);
+    [[maybe_unused]] amrex::Real dsepdz  = -ytot1 * (state.y_e * df_dt * state.rho  + df_t);
+
+    amrex::Real eele    = state.y_e*free + state.T * sele;
+    amrex::Real deepdt  = state.T * dsepdt;
+    amrex::Real deepdd  = x * df_d + state.T * dsepdd;
+    [[maybe_unused]] amrex::Real deepda  = -state.y_e * ytot1 * (free +  df_d * din) + state.T * dsepda;
+    [[maybe_unused]] amrex::Real deepdz  = ytot1* (free + state.y_e * df_d * state.rho) + state.T * dsepdz;
 
     if constexpr (has_pressure<T>::value) {
         state.p    = state.p + pele;
@@ -477,40 +477,40 @@ void apply_ions (T& state)
 {
     using namespace helmholtz;
 
-    constexpr Real pi      = 3.1415926535897932384e0_rt;
-    constexpr Real sioncon = (2.0e0_rt * pi * amu * kerg)/(h*h);
-    constexpr Real kergavo = kerg * avo_eos;
-
-    Real deni = 1.0e0_rt / state.rho;
-    Real tempi = 1.0e0_rt / state.T;
-
-    Real ytot1   = 1.0e0_rt / state.abar;
-    Real xni     = avo_eos * ytot1 * state.rho;
-    Real dxnidd  = avo_eos * ytot1;
-    [[maybe_unused]] Real dxnida  = -xni * ytot1;
-
-    Real kt = kerg * state.T;
-
-    Real pion    = xni * kt;
-    Real dpiondd = dxnidd * kt;
-    Real dpiondt = xni * kerg;
-    [[maybe_unused]] Real dpionda = dxnida * kt;
-    [[maybe_unused]] Real dpiondz = 0.0e0_rt;
-
-    Real eion    = 1.5e0_rt * pion * deni;
-    Real deiondd = (1.5e0_rt * dpiondd - eion) * deni;
-    Real deiondt = 1.5e0_rt * dpiondt * deni;
-    [[maybe_unused]] Real deionda = 1.5e0_rt * dpionda * deni;
-    [[maybe_unused]] Real deiondz = 0.0e0_rt;
-
-    Real x       = state.abar * state.abar * std::sqrt(state.abar) * deni / avo_eos;
-    Real s       = sioncon * state.T;
-    Real z       = x * s * std::sqrt(s);
-    Real y       = std::log(z);
-    Real sion    = (pion * deni + eion) * tempi + kergavo * ytot1 * y;
-    Real dsiondd = (dpiondd * deni - pion * deni * deni + deiondd) * tempi -
+    constexpr amrex::Real pi      = 3.1415926535897932384e0_rt;
+    constexpr amrex::Real sioncon = (2.0e0_rt * pi * amu * kerg)/(h*h);
+    constexpr amrex::Real kergavo = kerg * avo_eos;
+
+    amrex::Real deni = 1.0e0_rt / state.rho;
+    amrex::Real tempi = 1.0e0_rt / state.T;
+
+    amrex::Real ytot1   = 1.0e0_rt / state.abar;
+    amrex::Real xni     = avo_eos * ytot1 * state.rho;
+    amrex::Real dxnidd  = avo_eos * ytot1;
+    [[maybe_unused]] amrex::Real dxnida  = -xni * ytot1;
+
+    amrex::Real kt = kerg * state.T;
+
+    amrex::Real pion    = xni * kt;
+    amrex::Real dpiondd = dxnidd * kt;
+    amrex::Real dpiondt = xni * kerg;
+    [[maybe_unused]] amrex::Real dpionda = dxnida * kt;
+    [[maybe_unused]] amrex::Real dpiondz = 0.0e0_rt;
+
+    amrex::Real eion    = 1.5e0_rt * pion * deni;
+    amrex::Real deiondd = (1.5e0_rt * dpiondd - eion) * deni;
+    amrex::Real deiondt = 1.5e0_rt * dpiondt * deni;
+    [[maybe_unused]] amrex::Real deionda = 1.5e0_rt * dpionda * deni;
+    [[maybe_unused]] amrex::Real deiondz = 0.0e0_rt;
+
+    amrex::Real x       = state.abar * state.abar * std::sqrt(state.abar) * deni / avo_eos;
+    amrex::Real s       = sioncon * state.T;
+    amrex::Real z       = x * s * std::sqrt(s);
+    amrex::Real y       = std::log(z);
+    amrex::Real sion    = (pion * deni + eion) * tempi + kergavo * ytot1 * y;
+    amrex::Real dsiondd = (dpiondd * deni - pion * deni * deni + deiondd) * tempi -
                    kergavo * deni * ytot1;
-    Real dsiondt = (dpiondt * deni + deiondt) * tempi -
+    amrex::Real dsiondt = (dpiondt * deni + deiondt) * tempi -
                    (pion * deni + eion) * tempi * tempi +
                    1.5e0_rt * kergavo * tempi * ytot1;
 
@@ -553,19 +553,19 @@ void apply_radiation (T& state)
 {
     using namespace helmholtz;
 
-    constexpr Real clight  = 2.99792458e10_rt;
+    constexpr amrex::Real clight  = 2.99792458e10_rt;
 #ifdef RADIATION
-    constexpr Real ssol    = 0.0e0_rt;
+    constexpr amrex::Real ssol    = 0.0e0_rt;
 #else
-    constexpr Real ssol    = 5.67051e-5_rt;
+    constexpr amrex::Real ssol    = 5.67051e-5_rt;
 #endif
-    constexpr Real asol    = 4.0e0_rt * ssol / clight;
-    constexpr Real asoli3  = asol/3.0e0_rt;
+    constexpr amrex::Real asol    = 4.0e0_rt * ssol / clight;
+    constexpr amrex::Real asoli3  = asol/3.0e0_rt;
 
-    Real deni = 1.0e0_rt / state.rho;
-    Real tempi = 1.0e0_rt / state.T;
+    amrex::Real deni = 1.0e0_rt / state.rho;
+    amrex::Real tempi = 1.0e0_rt / state.T;
 
-    Real prad = asoli3 * state.T * state.T * state.T * state.T;
+    amrex::Real prad = asoli3 * state.T * state.T * state.T * state.T;
 
     // In low density material, this radiation pressure becomes unphysically high.
     // For rho ~ 1 g/cc and T ~ 1e9, then this radiation pressure will lead to a
@@ -579,20 +579,20 @@ void apply_radiation (T& state)
         prad = prad * 0.5e0_rt * (1.0e0_rt + std::tanh((state.rho - prad_limiter_rho_c) / prad_limiter_delta_rho));
     }
 
-    Real dpraddd = 0.0e0_rt;
-    Real dpraddt = 4.0e0_rt * prad * tempi;
-    [[maybe_unused]] Real dpradda = 0.0e0_rt;
-    [[maybe_unused]] Real dpraddz = 0.0e0_rt;
+    amrex::Real dpraddd = 0.0e0_rt;
+    amrex::Real dpraddt = 4.0e0_rt * prad * tempi;
+    [[maybe_unused]] amrex::Real dpradda = 0.0e0_rt;
+    [[maybe_unused]] amrex::Real dpraddz = 0.0e0_rt;
 
-    Real erad    = 3.0e0_rt * prad*deni;
-    Real deraddd = -erad * deni;
-    Real deraddt = 3.0e0_rt * dpraddt * deni;
-    [[maybe_unused]] Real deradda = 0.0e0_rt;
-    [[maybe_unused]] Real deraddz = 0.0e0_rt;
+    amrex::Real erad    = 3.0e0_rt * prad*deni;
+    amrex::Real deraddd = -erad * deni;
+    amrex::Real deraddt = 3.0e0_rt * dpraddt * deni;
+    [[maybe_unused]] amrex::Real deradda = 0.0e0_rt;
+    [[maybe_unused]] amrex::Real deraddz = 0.0e0_rt;
 
-    Real srad    = (prad * deni + erad) * tempi;
-    Real dsraddd = (dpraddd * deni - prad * deni * deni + deraddd) * tempi;
-    Real dsraddt = (dpraddt * deni + deraddt - srad) * tempi;
+    amrex::Real srad    = (prad * deni + erad) * tempi;
+    amrex::Real dsraddd = (dpraddd * deni - prad * deni * deni + deraddd) * tempi;
+    amrex::Real dsraddt = (dpraddt * deni + deraddt - srad) * tempi;
 
     // Note that unlike the other terms, radiation
     // sets these terms instead of adding to them,
@@ -638,67 +638,67 @@ void apply_coulomb_corrections (T& state)
     using namespace helmholtz;
 
     // Constants used for the Coulomb corrections
-    constexpr Real a1 = -0.898004e0_rt;
-    constexpr Real b1 =  0.96786e0_rt;
-    constexpr Real c1 =  0.220703e0_rt;
-    constexpr Real d1 = -0.86097e0_rt;
-    constexpr Real e1 =  2.5269e0_rt;
-    constexpr Real a2 =  0.29561e0_rt;
-    constexpr Real b2 =  1.9885e0_rt;
-    constexpr Real c2 =  0.288675e0_rt;
-    constexpr Real qe   = 4.8032042712e-10_rt;
-    constexpr Real esqu = qe * qe;
-    constexpr Real onethird = 1.0e0_rt/3.0e0_rt;
-    constexpr Real forth = 4.0e0_rt/3.0e0_rt;
-    constexpr Real pi    = 3.1415926535897932384e0_rt;
-
-    [[maybe_unused]] Real pcoul    = 0.e0_rt;
-    [[maybe_unused]] Real dpcouldd = 0.e0_rt;
-    [[maybe_unused]] Real dpcouldt = 0.e0_rt;
-    [[maybe_unused]] Real ecoul    = 0.e0_rt;
-    [[maybe_unused]] Real decouldd = 0.e0_rt;
-    [[maybe_unused]] Real decouldt = 0.e0_rt;
-    [[maybe_unused]] Real scoul    = 0.e0_rt;
-    [[maybe_unused]] Real dscouldd = 0.e0_rt;
-    [[maybe_unused]] Real dscouldt = 0.e0_rt;
-
-    [[maybe_unused]] Real dpcoulda = 0.e0_rt;
-    [[maybe_unused]] Real dpcouldz = 0.e0_rt;
-    [[maybe_unused]] Real decoulda = 0.e0_rt;
-    [[maybe_unused]] Real decouldz = 0.e0_rt;
-
-    Real s, x, y, z; // temporary variables
+    constexpr amrex::Real a1 = -0.898004e0_rt;
+    constexpr amrex::Real b1 =  0.96786e0_rt;
+    constexpr amrex::Real c1 =  0.220703e0_rt;
+    constexpr amrex::Real d1 = -0.86097e0_rt;
+    constexpr amrex::Real e1 =  2.5269e0_rt;
+    constexpr amrex::Real a2 =  0.29561e0_rt;
+    constexpr amrex::Real b2 =  1.9885e0_rt;
+    constexpr amrex::Real c2 =  0.288675e0_rt;
+    constexpr amrex::Real qe   = 4.8032042712e-10_rt;
+    constexpr amrex::Real esqu = qe * qe;
+    constexpr amrex::Real onethird = 1.0e0_rt/3.0e0_rt;
+    constexpr amrex::Real forth = 4.0e0_rt/3.0e0_rt;
+    constexpr amrex::Real pi    = 3.1415926535897932384e0_rt;
+
+    [[maybe_unused]] amrex::Real pcoul    = 0.e0_rt;
+    [[maybe_unused]] amrex::Real dpcouldd = 0.e0_rt;
+    [[maybe_unused]] amrex::Real dpcouldt = 0.e0_rt;
+    [[maybe_unused]] amrex::Real ecoul    = 0.e0_rt;
+    [[maybe_unused]] amrex::Real decouldd = 0.e0_rt;
+    [[maybe_unused]] amrex::Real decouldt = 0.e0_rt;
+    [[maybe_unused]] amrex::Real scoul    = 0.e0_rt;
+    [[maybe_unused]] amrex::Real dscouldd = 0.e0_rt;
+    [[maybe_unused]] amrex::Real dscouldt = 0.e0_rt;
+
+    [[maybe_unused]] amrex::Real dpcoulda = 0.e0_rt;
+    [[maybe_unused]] amrex::Real dpcouldz = 0.e0_rt;
+    [[maybe_unused]] amrex::Real decoulda = 0.e0_rt;
+    [[maybe_unused]] amrex::Real decouldz = 0.e0_rt;
+
+    amrex::Real s, x, y, z; // temporary variables
 
     // uniform background corrections only
     // from yakovlev & shalybkov 1989
     // lami is the average ion separation
     // plasg is the plasma coupling parameter
 
-    Real ytot1 = 1.0e0_rt / state.abar;
-    Real xni     = avo_eos * ytot1 * state.rho;
-    Real dxnidd  = avo_eos * ytot1;
-    [[maybe_unused]] Real dxnida  = -xni * ytot1;
+    amrex::Real ytot1 = 1.0e0_rt / state.abar;
+    amrex::Real xni     = avo_eos * ytot1 * state.rho;
+    amrex::Real dxnidd  = avo_eos * ytot1;
+    [[maybe_unused]] amrex::Real dxnida  = -xni * ytot1;
 
-    Real kt      = kerg * state.T;
-    Real ktinv   = 1.0e0_rt / kt;
+    amrex::Real kt      = kerg * state.T;
+    amrex::Real ktinv   = 1.0e0_rt / kt;
 
     z             = forth * pi;
     s             = z * xni;
-    Real dsdd     = z * dxnidd;
-    [[maybe_unused]] Real dsda     = z * dxnida;
+    amrex::Real dsdd     = z * dxnidd;
+    [[maybe_unused]] amrex::Real dsda     = z * dxnida;
 
-    Real lami     = 1.0e0_rt / std::pow(s, onethird);
-    Real inv_lami = 1.0e0_rt / lami;
+    amrex::Real lami     = 1.0e0_rt / std::pow(s, onethird);
+    amrex::Real inv_lami = 1.0e0_rt / lami;
     z             = -onethird * lami;
-    Real lamidd   = z * dsdd / s;
-    [[maybe_unused]] Real lamida   = z * dsda / s;
+    amrex::Real lamidd   = z * dsdd / s;
+    [[maybe_unused]] amrex::Real lamida   = z * dsda / s;
 
-    Real plasg    = state.zbar * state.zbar * esqu * ktinv * inv_lami;
+    amrex::Real plasg    = state.zbar * state.zbar * esqu * ktinv * inv_lami;
     z             = -plasg * inv_lami;
-    Real plasgdd  = z * lamidd;
-    Real plasgdt  = -plasg*ktinv * kerg;
-    [[maybe_unused]] Real plasgda  = z * lamida;
-    [[maybe_unused]] Real plasgdz  = 2.0e0_rt * plasg/state.zbar;
+    amrex::Real plasgdd  = z * lamidd;
+    amrex::Real plasgdt  = -plasg*ktinv * kerg;
+    [[maybe_unused]] amrex::Real plasgda  = z * lamida;
+    [[maybe_unused]] amrex::Real plasgdz  = 2.0e0_rt * plasg/state.zbar;
 
     // .yakovlev & shalybkov 1989 equations 82, 85, 86, 87
     if (plasg >= 1.0e0_rt)
@@ -734,11 +734,11 @@ void apply_coulomb_corrections (T& state)
     else if (plasg < 1.0e0_rt)
     {
 
-        Real pion    = xni * kt;
-        Real dpiondd = dxnidd * kt;
-        Real dpiondt = xni * kerg;
-        Real dpionda = dxnida * kt;
-        Real dpiondz = 0.0e0_rt;
+        amrex::Real pion    = xni * kt;
+        amrex::Real dpiondd = dxnidd * kt;
+        amrex::Real dpiondt = xni * kerg;
+        amrex::Real dpionda = dxnida * kt;
+        amrex::Real dpiondz = 0.0e0_rt;
 
         x        = plasg * std::sqrt(plasg);
         y        = std::pow(plasg, b2);
@@ -768,8 +768,8 @@ void apply_coulomb_corrections (T& state)
     // Disable Coulomb corrections if they cause
     // the energy or pressure to go negative.
 
-    Real p_temp = std::numeric_limits<Real>::max();
-    Real e_temp = std::numeric_limits<Real>::max();
+    amrex::Real p_temp = std::numeric_limits<amrex::Real>::max();
+    amrex::Real e_temp = std::numeric_limits<amrex::Real>::max();
 
     if constexpr (has_pressure<T>::value) {
         p_temp = state.p + pcoul;
@@ -832,7 +832,7 @@ void apply_coulomb_corrections (T& state)
 template <typename I, typename T>
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
 void prepare_for_iterations (I input, T& state,
-                             bool& single_iter, Real& v_want, Real& v1_want, Real& v2_want,
+                             bool& single_iter, amrex::Real& v_want, amrex::Real& v1_want, amrex::Real& v2_want,
                              int& var, int& dvar, int& var1, int& var2)
 {
     using namespace helmholtz;
@@ -926,13 +926,13 @@ void prepare_for_iterations (I input, T& state,
 template <typename T>
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
 void single_iter_update (T& state, int var, int dvar,
-                         Real v_want, bool& converged)
+                         amrex::Real v_want, bool& converged)
 {
 
     using namespace helmholtz;
     using namespace EOS;
 
-    Real x, v = 0.0_rt, dvdx = 0.0_rt, xtol, smallx;
+    amrex::Real x, v = 0.0_rt, dvdx = 0.0_rt, xtol, smallx;
 
     if (dvar == itemp) {
 
@@ -1001,7 +1001,7 @@ void single_iter_update (T& state, int var, int dvar,
     }
 
     // Now do the calculation for the next guess for T/rho
-    Real xnew = x - (v - v_want) / dvdx;
+    amrex::Real xnew = x - (v - v_want) / dvdx;
 
     // Don't let the temperature/density change by more than a factor of two
     xnew = amrex::max(0.5_rt * x, amrex::min(xnew, 2.0_rt * x));
@@ -1021,7 +1021,7 @@ void single_iter_update (T& state, int var, int dvar,
 
     // Compute the error from the last iteration
 
-    Real error = std::abs((xnew - x) / x);
+    amrex::Real error = std::abs((xnew - x) / x);
 
     if (error < xtol) converged = true;
 }
@@ -1031,17 +1031,17 @@ void single_iter_update (T& state, int var, int dvar,
 template <typename T>
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
 void double_iter_update (T& state, int var1, int var2,
-                         Real v1_want, Real v2_want, bool& converged)
+                         amrex::Real v1_want, amrex::Real v2_want, bool& converged)
 {
     using namespace helmholtz;
     using namespace EOS;
 
-    Real v1 = 0.0_rt, dv1dt = 0.0_rt, dv1dr = 0.0_rt, v2 = 0.0_rt, dv2dt = 0.0_rt, dv2dr = 0.0_rt;
+    amrex::Real v1 = 0.0_rt, dv1dt = 0.0_rt, dv1dr = 0.0_rt, v2 = 0.0_rt, dv2dt = 0.0_rt, dv2dr = 0.0_rt;
 
     // Figure out which variables we're using
 
-    Real told = state.T;
-    Real rold = state.rho;
+    amrex::Real told = state.T;
+    amrex::Real rold = state.rho;
 
     if (var1 == ipres) {
         if constexpr (has_pressure<T>::value) {
@@ -1102,8 +1102,8 @@ void double_iter_update (T& state, int var1, int var2,
     }
 
     // Two functions, f and g, to iterate over
-    Real v1i = v1_want - v1;
-    Real v2i = v2_want - v2;
+    amrex::Real v1i = v1_want - v1;
+    amrex::Real v2i = v2_want - v2;
 
     //
     // 0 = f + dfdr * delr + dfdt * delt
@@ -1111,11 +1111,11 @@ void double_iter_update (T& state, int var1, int var2,
     //
 
     // note that dfi/dT = - df/dT
-    Real delr = (-v1i * dv2dt + v2i * dv1dt) / (dv2dr * dv1dt - dv2dt * dv1dr);
+    amrex::Real delr = (-v1i * dv2dt + v2i * dv1dt) / (dv2dr * dv1dt - dv2dt * dv1dr);
 
-    Real rnew = rold + delr;
+    amrex::Real rnew = rold + delr;
 
-    Real tnew = told + (v1i - dv1dr * delr) / dv1dt;
+    amrex::Real tnew = told + (v1i - dv1dr * delr) / dv1dt;
 
     // Don't let the temperature or density change by more
     // than a factor of two
@@ -1131,8 +1131,8 @@ void double_iter_update (T& state, int var1, int var2,
     state.T = tnew;
 
     // Compute the errors
-    Real error1 = std::abs((rnew - rold) / rold);
-    Real error2 = std::abs((tnew - told) / told);
+    amrex::Real error1 = std::abs((rnew - rold) / rold);
+    amrex::Real error2 = std::abs((tnew - told) / told);
 
     if (error1 < dtol && error2 < ttol) converged = true;
 }
@@ -1142,7 +1142,7 @@ void double_iter_update (T& state, int var1, int var2,
 template <typename I, typename T>
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
 void finalize_state (I input, T& state,
-                     Real v_want, Real v1_want, Real v2_want)
+                     amrex::Real v_want, amrex::Real v1_want, amrex::Real v2_want)
 {
     using namespace helmholtz;
 
@@ -1157,8 +1157,8 @@ void finalize_state (I input, T& state,
         state.cv = state.dedT;
 
         if constexpr (has_pressure<T>::value) {
-            Real chit = state.T / state.p * state.dpdT;
-            Real chid = state.dpdr * state.rho / state.p;
+            amrex::Real chit = state.T / state.p * state.dpdT;
+            amrex::Real chid = state.dpdr * state.rho / state.p;
 
             state.gam1 = (chit * (state.p / state.rho)) * (chit / (state.T * state.cv)) + chid;
             state.cp = state.cv * state.gam1 / chid;
@@ -1246,7 +1246,7 @@ void actual_eos (I input, T& state)
 
     bool single_iter, converged;
     int var{}, dvar{}, var1, var2;
-    Real v_want{}, v1_want{}, v2_want{};
+    amrex::Real v_want{}, v1_want{}, v2_want{};
 
     prepare_for_iterations(input, state, single_iter, v_want, v1_want, v2_want, var, dvar, var1, var2);
 
@@ -1304,8 +1304,8 @@ void actual_eos_init ()
 {
     using namespace helmholtz;
 
-    Real dth, dt2, dti, dt2i;
-    Real dd, dd2, ddi, dd2i;
+    amrex::Real dth, dt2, dti, dt2i;
+    amrex::Real dd, dd2, ddi, dd2i;
 
     // Read in the runtime parameters
 
@@ -1317,10 +1317,10 @@ void actual_eos_init ()
     //    read the helmholtz free energy table
 
     for (int j = 0; j < jmax; ++j) {
-        Real tsav = tlo + j * tstp;
+        amrex::Real tsav = tlo + j * tstp;
         t[j] = std::pow(10.0e0_rt, tsav);
         for (int i = 0; i < imax; ++i) {
-            Real dsav = dlo + i * dstp;
+            amrex::Real dsav = dlo + i * dstp;
             d[i] = std::pow(10.0e0_rt, dsav);
         }
     }
@@ -1330,10 +1330,10 @@ void actual_eos_init ()
     // buffer, broadcast that, and then copy that into the managed
     // memory.
 
-    amrex::Vector<Real> f_local(static_cast<size_t>(9) * imax * jmax);
-    amrex::Vector<Real> dpdf_local(static_cast<size_t>(4) * imax * jmax);
-    amrex::Vector<Real> ef_local(static_cast<size_t>(4) * imax * jmax);
-    amrex::Vector<Real> xf_local(static_cast<size_t>(4) * imax * jmax);
+    amrex::Vector<amrex::Real> f_local(static_cast<size_t>(9) * imax * jmax);
+    amrex::Vector<amrex::Real> dpdf_local(static_cast<size_t>(4) * imax * jmax);
+    amrex::Vector<amrex::Real> ef_local(static_cast<size_t>(4) * imax * jmax);
+    amrex::Vector<amrex::Real> xf_local(static_cast<size_t>(4) * imax * jmax);
 
     if (amrex::ParallelDescriptor::IOProcessor()) {
 
diff --git a/EOS/multigamma/actual_eos.H b/EOS/multigamma/actual_eos.H
index 8508b284d5..e673474f7c 100644
--- a/EOS/multigamma/actual_eos.H
+++ b/EOS/multigamma/actual_eos.H
@@ -73,11 +73,11 @@ void actual_eos (I input, T& state)
     static_assert(std::is_same_v<I, eos_input_t>, "input must be an eos_input_t");
 
     // Get the mass of a nucleon from Avogadro's number.
-    const Real m_nucleon = 1.0_rt / C::n_A;
+    const amrex::Real m_nucleon = 1.0_rt / C::n_A;
 
     // Special gamma factors
-    Real sumY_gm1 = 0.0_rt;
-    Real sumYg_gm1 = 0.0_rt;
+    amrex::Real sumY_gm1 = 0.0_rt;
+    amrex::Real sumYg_gm1 = 0.0_rt;
     for (int n = 0; n < NumSpec; ++n) {
         sumY_gm1 = sumY_gm1 + state.xn[n] * aion_inv[n] / (gammas[n] - 1.0_rt);
         sumYg_gm1 = sumYg_gm1 + state.xn[n] * gammas[n] * aion_inv[n] / (gammas[n] - 1.0_rt);
@@ -88,7 +88,7 @@ void actual_eos (I input, T& state)
     // and temp as needed from the inputs.
     //-------------------------------------------------------------------------
 
-    Real temp, dens;
+    amrex::Real temp, dens;
 
     switch (input) {
 
@@ -201,9 +201,9 @@ void actual_eos (I input, T& state)
 
     state.T   = temp;
     state.rho = dens;
-    Real pres = dens * (C::k_B / m_nucleon) * temp / state.abar;
-    Real ener = (C::k_B / m_nucleon) * temp * sumY_gm1;
-    Real enth = ener + pres / dens;
+    amrex::Real pres = dens * (C::k_B / m_nucleon) * temp / state.abar;
+    amrex::Real ener = (C::k_B / m_nucleon) * temp * sumY_gm1;
+    amrex::Real enth = ener + pres / dens;
 
     // Compute the pressure simply from the ideal gas law, and the
     // specific internal energy using the gamma-law EOS relation.
@@ -229,29 +229,29 @@ void actual_eos (I input, T& state)
     }
 
     // Compute the thermodynamic derivatives and specific heats
-    Real dpdT = pres / temp;
-    Real dpdr = pres / dens;
+    amrex::Real dpdT = pres / temp;
+    amrex::Real dpdr = pres / dens;
     if constexpr (has_pressure<T>::value) {
         state.dpdT = dpdT;
         state.dpdr = dpdr;
     }
 
-    Real dedT = ener / temp;
-    Real dedr = 0.0_rt;
+    amrex::Real dedT = ener / temp;
+    amrex::Real dedr = 0.0_rt;
     if constexpr (has_energy<T>::value) {
         state.dedT = dedT;
         state.dedr = dedr;
     }
 
-    Real dsdT = 0.0_rt;
-    Real dsdr = 0.0_rt;
+    amrex::Real dsdT = 0.0_rt;
+    amrex::Real dsdr = 0.0_rt;
     if constexpr (has_entropy<T>::value) {
         state.dsdT = dsdT;
         state.dsdr = dsdr;
     }
 
-    Real dhdT = dedT + dpdT / dens;
-    Real dhdr = 0.0_rt;
+    amrex::Real dhdT = dedT + dpdT / dens;
+    amrex::Real dhdr = 0.0_rt;
     if constexpr (has_enthalpy<T>::value) {
         state.dhdT = dhdT;
         state.dhdr = dhdr;
diff --git a/EOS/polytrope/actual_eos.H b/EOS/polytrope/actual_eos.H
index f517b9dc7a..3f0108a8e5 100644
--- a/EOS/polytrope/actual_eos.H
+++ b/EOS/polytrope/actual_eos.H
@@ -92,7 +92,7 @@ bool is_input_valid (I input)
 //---------------------------------------------------------------------------
 
 inline
-void eos_get_polytrope_parameters (int& polytrope_out, Real& gamma_out, Real& K_out, Real& mu_e_out)
+void eos_get_polytrope_parameters (int& polytrope_out, amrex::Real& gamma_out, amrex::Real& K_out, amrex::Real& mu_e_out)
 {
 
     polytrope_out = polytrope;
@@ -124,25 +124,25 @@ void actual_eos (I input, T& state)
 {
     static_assert(std::is_same_v<I, eos_input_t>, "input must be an eos_input_t");
 
-    Real dens = state.rho;
-    Real temp = state.T;
+    amrex::Real dens = state.rho;
+    amrex::Real temp = state.T;
 
-    Real pres = 1.0_rt;
+    amrex::Real pres = 1.0_rt;
     if constexpr (has_pressure<T>::value) {
         pres = state.p;
     }
 
-    Real enth = 1.0_rt;
+    amrex::Real enth = 1.0_rt;
     if constexpr (has_enthalpy<T>::value) {
         enth = state.h;
     }
 
-    Real eint = 1.0_rt;
+    amrex::Real eint = 1.0_rt;
     if constexpr (has_energy<T>::value) {
         eint = state.e;
     }
 
-    Real entr = 1.0_rt;
+    amrex::Real entr = 1.0_rt;
     if constexpr (has_entropy<T>::value) {
         entr = state.s;
     }
diff --git a/EOS/primordial_chem/actual_eos.H b/EOS/primordial_chem/actual_eos.H
index e03733076d..3d50fa76b9 100644
--- a/EOS/primordial_chem/actual_eos.H
+++ b/EOS/primordial_chem/actual_eos.H
@@ -139,12 +139,12 @@ void actual_eos (I input, T& state)
 {
     static_assert(std::is_same_v<I, eos_input_t>, "input must be either eos_input_rt or eos_input_re");
 
-    const Real gasconstant = C::n_A * C::k_B;
-    const Real protonmass = C::m_p;
+    const amrex::Real gasconstant = C::n_A * C::k_B;
+    const amrex::Real protonmass = C::m_p;
     // Special gamma factors
-    Real sum_Abarinv = 0.0_rt;
-    Real sum_gammasinv = 0.0_rt;
-    Real rhotot = 0.0_rt;
+    amrex::Real sum_Abarinv = 0.0_rt;
+    amrex::Real sum_gammasinv = 0.0_rt;
+    amrex::Real rhotot = 0.0_rt;
 
     for (int n = 0; n < NumSpec; ++n) {
         rhotot += state.xn[n]*spmasses[n];
@@ -164,9 +164,9 @@ void actual_eos (I input, T& state)
     // and temp as needed from the inputs.
     //-------------------------------------------------------------------------
 
-    Real temp = NAN;
-    Real dens = NAN;
-    Real eint = NAN;
+    amrex::Real temp = NAN;
+    amrex::Real dens = NAN;
+    amrex::Real eint = NAN;
 
     switch (input) {
 
@@ -272,7 +272,7 @@ void actual_eos (I input, T& state)
     }
 
     if constexpr (has_pressure<T>::value) {
-        Real pressure = state.rho * eint / sum_gammasinv;
+        amrex::Real pressure = state.rho * eint / sum_gammasinv;
         state.p = pressure;
 
         state.dpdT = pressure / temp;
@@ -283,8 +283,8 @@ void actual_eos (I input, T& state)
         }
     }
 
-    Real dedT = sum_gammasinv * sum_Abarinv * gasconstant;
-    Real dedr = 0.0_rt;
+    amrex::Real dedT = sum_gammasinv * sum_Abarinv * gasconstant;
+    amrex::Real dedr = 0.0_rt;
     if constexpr (has_energy<T>::value) {
         state.dedT = dedT;
         state.dedr = dedr;
diff --git a/EOS/tillotson/actual_eos.H b/EOS/tillotson/actual_eos.H
index 653416d4f7..c9df21832a 100644
--- a/EOS/tillotson/actual_eos.H
+++ b/EOS/tillotson/actual_eos.H
@@ -62,25 +62,25 @@ void actual_eos (I input, T& state)
         state.e = amrex::max(eos_e_0 + eos_c_v * state.T, EOSData::mine);
     }
 
-    Real eta = state.rho / eos_rho_0;
-    Real mu = eta - 1.0_rt;
-    Real omega = state.e / (eos_e_0 * eta * eta) + 1.0_rt;
-    Real z = (1.0_rt / eta - 1.0_rt);
+    amrex::Real eta = state.rho / eos_rho_0;
+    amrex::Real mu = eta - 1.0_rt;
+    amrex::Real omega = state.e / (eos_e_0 * eta * eta) + 1.0_rt;
+    amrex::Real z = (1.0_rt / eta - 1.0_rt);
 
-    Real g_c = eos_la + eos_lb / omega;
-    Real g_e = eos_la + eos_lb / omega * std::exp(-eos_beta * z * z);
+    amrex::Real g_c = eos_la + eos_lb / omega;
+    amrex::Real g_e = eos_la + eos_lb / omega * std::exp(-eos_beta * z * z);
 
-    Real P_c = g_c * state.e * state.rho + eos_A * mu + eos_B * mu * mu;
-    Real P_e = g_e * state.e * state.rho + eos_A * mu * std::exp(-eos_alpha * z - eos_beta * z * z);
+    amrex::Real P_c = g_c * state.e * state.rho + eos_A * mu + eos_B * mu * mu;
+    amrex::Real P_e = g_e * state.e * state.rho + eos_A * mu * std::exp(-eos_alpha * z - eos_beta * z * z);
 
     // Floor the pressure since the above expressions can be negative
 
     P_c = amrex::max(P_c, EOSData::minp);
     P_e = amrex::max(P_e, EOSData::minp);
 
-    Real cs_squared_c = (g_c + 2.0_rt * eos_lb * (omega - 1.0_rt) / (omega * omega)) * state.e +
+    amrex::Real cs_squared_c = (g_c + 2.0_rt * eos_lb * (omega - 1.0_rt) / (omega * omega)) * state.e +
                         (eos_A + 2.0_rt * eos_B * mu) / eos_rho_0;
-    Real cs_squared_e = (g_e + (2.0_rt * eos_lb / omega) * std::exp(-eos_beta * z * z) *
+    amrex::Real cs_squared_e = (g_e + (2.0_rt * eos_lb / omega) * std::exp(-eos_beta * z * z) *
                          ((omega - 1.0_rt) / omega + eos_beta * z / eta)) * state.e +
                         (eos_A / eos_rho_0) * std::exp(-eos_alpha * z - eos_beta * z * z) *
                         (1.0_rt + mu / (eta * eta) * (eos_alpha + 2.0_rt * eos_beta * z));
@@ -90,12 +90,12 @@ void actual_eos (I input, T& state)
     cs_squared_c = amrex::max(cs_squared_c, EOSData::minp / state.rho);
     cs_squared_e = amrex::max(cs_squared_e, EOSData::minp / state.rho);
 
-    Real dPde_c = eos_la * state.rho + eos_lb * state.rho / (omega * omega);
-    Real dPde_e = eos_la * state.rho + eos_lb * state.rho / (omega * omega) * std::exp(-eos_beta * z * z);
+    amrex::Real dPde_c = eos_la * state.rho + eos_lb * state.rho / (omega * omega);
+    amrex::Real dPde_e = eos_la * state.rho + eos_lb * state.rho / (omega * omega) * std::exp(-eos_beta * z * z);
 
     // At this point we should have valid (rho, e); compute P = P(rho, e).
 
-    Real pres, dpdr_e, dpde;
+    amrex::Real pres, dpdr_e, dpde;
 
     if (state.rho >= eos_rho_0 || (state.rho < eos_rho_0 && state.e < eos_e_s)) {
         // Region I, II (compressed state)
@@ -114,9 +114,9 @@ void actual_eos (I input, T& state)
     else {
         // Region III (interpolate between compressed and expanded state)
 
-        Real denom = (eos_e_s_prime - eos_e_s);
-        Real xi_c = (state.e - eos_e_s) / denom;
-        Real xi_e = (eos_e_s_prime - state.e) / denom;
+        amrex::Real denom = (eos_e_s_prime - eos_e_s);
+        amrex::Real xi_c = (state.e - eos_e_s) / denom;
+        amrex::Real xi_e = (eos_e_s_prime - state.e) / denom;
 
         pres = P_c * xi_c + P_e * xi_e;
         dpdr_e = cs_squared_c * xi_c + cs_squared_e * xi_e;
diff --git a/interfaces/burn_type.H b/interfaces/burn_type.H
index 5618c52e7f..05eef0def1 100644
--- a/interfaces/burn_type.H
+++ b/interfaces/burn_type.H
@@ -12,7 +12,7 @@
 
 #include <ArrayUtilities.H>
 
-using namespace amrex;
+using namespace amrex::literals;
 using namespace network_rp;
 
 // A generic structure holding data necessary to do a nuclear burn.
@@ -32,7 +32,7 @@ const int net_ienuc = NumSpec + 1;
 
 // this is the data type that is used to get the ydots from the actual
 // RHS of the network, regardless of Strang or SDC
-typedef amrex::Array1D<Real, 1, neqs> YdotNetArray1D;
+typedef amrex::Array1D<amrex::Real, 1, neqs> YdotNetArray1D;
 
 // these indices represent the order that the conserved state comes
 // into the ODE integration from the hydro code.
@@ -65,26 +65,26 @@ struct burn_t
 
   // on exit, this will be set to the integration time we achieved
 
-  Real time{};
+  amrex::Real time{};
 
   // this first group are the quantities the network RHS uses
 
-  Real rho{};
-  Real T{};
-  Real e{};
-  Real xn[NumSpec]{};
+  amrex::Real rho{};
+  amrex::Real T{};
+  amrex::Real e{};
+  amrex::Real xn[NumSpec]{};
 #if NAUX_NET > 0
-  Real aux[NumAux]{};
+  amrex::Real aux[NumAux]{};
 #endif
 
   // scaling used to make e we integrate dimensionless
-  Real e_scale{};
+  amrex::Real e_scale{};
 
   // derivatives needed for calling the EOS
-  Real dedr{};
-  Real dedT{};
-  Real dedA{};
-  Real dedZ{};
+  amrex::Real dedr{};
+  amrex::Real dedT{};
+  amrex::Real dedA{};
+  amrex::Real dedZ{};
 
   // for diagnostics / error reporting
 
@@ -93,7 +93,7 @@ struct burn_t
   int k{};
 #ifdef NONAKA_PLOT
   int level{};
-  Real reference_time{};
+  amrex::Real reference_time{};
 #ifdef STRANG
   int strang_half{};
 #endif
@@ -104,17 +104,17 @@ struct burn_t
   // y is the input conserved state.  We will keep this state updated
   // in time as we integrate, such that upon output it will be the
   // final conserved state.
-  Real y[SVAR]{};
+  amrex::Real y[SVAR]{};
 
   // we need to store a copy of the original state as well so we can
   // handle the non-evolved state evolution
-  Real rho_orig{};
+  amrex::Real rho_orig{};
 
-  Real rhoe_orig{};
+  amrex::Real rhoe_orig{};
 
   // ydot_a are the advective terms that will modify the state y due
   // to hydrodynamics over the timestep.
-  Real ydot_a[SVAR]{};
+  amrex::Real ydot_a[SVAR]{};
 
   int sdc_iter{};
   int num_sdc_iters{};
@@ -122,25 +122,25 @@ struct burn_t
   // for drive_initial_convection, we will fix T during the
   // integration to a passed in value.  We will interpret a positive
   // T_fixed as setting this feature.
-  Real T_fixed{-1.0};
+  amrex::Real T_fixed{-1.0};
 
   // all coupling types need the specific heats to transform the
   // reaction Jacobian elements from T to e
-  Real cv{};
+  amrex::Real cv{};
 
   // dx is useful for estimating timescales for equilibriation
-  Real dx{};
+  amrex::Real dx{};
 
-  Real mu{};
-  Real mu_e{};
-  Real y_e{};
-  Real eta{};
-  Real abar{};
-  Real zbar{};
+  amrex::Real mu{};
+  amrex::Real mu_e{};
+  amrex::Real y_e{};
+  amrex::Real eta{};
+  amrex::Real abar{};
+  amrex::Real zbar{};
 
 #ifdef NSE_NET
-  Real mu_p{};
-  Real mu_n{};
+  amrex::Real mu_p{};
+  amrex::Real mu_n{};
 #endif
 
 #ifdef NSE
@@ -298,7 +298,7 @@ void normalize_abundances_sdc_burn (BurnT& state)
 
     // we have state.y[SRHO] to constrain to
 
-    Real rhoX_sum = 0.0_rt;
+    amrex::Real rhoX_sum = 0.0_rt;
     for (int n = 0; n < NumSpec; n++) {
         state.y[SFS+n] = amrex::max(amrex::min(state.y[SFS+n], state.y[SRHO]), 0.0_rt);
         rhoX_sum += state.y[SFS+n];
@@ -333,9 +333,9 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 void normalize_abundances_burn (BurnT& state)
 {
 
-    Real sum = 0.0_rt;
+    amrex::Real sum = 0.0_rt;
     for (double &x : state.xn) {
-        x = amrex::max(small_x, amrex::min(1.0_rt, x));
+        x = amrex::max(network_rp::small_x, amrex::min(1.0_rt, x));
         sum += x;
     }
     for (double &x : state.xn) {
diff --git a/interfaces/eos.H b/interfaces/eos.H
index 850f55b394..c5ecfd7cf4 100644
--- a/interfaces/eos.H
+++ b/interfaces/eos.H
@@ -9,13 +9,12 @@
 #include <actual_eos.H>
 #include <AMReX_Algorithm.H>
 
-using namespace amrex;
 
 // EOS initialization routine: read in general EOS parameters, then
 // call any specific initialization used by the EOS.
 
 AMREX_INLINE
-void eos_init (Real& small_temp_in, Real& small_dens_in)
+void eos_init (amrex::Real& small_temp_in, amrex::Real& small_dens_in)
 {
   // Allocate and set default values
 
@@ -54,8 +53,8 @@ void eos_init (Real& small_temp_in, Real& small_dens_in)
 AMREX_INLINE
 void eos_init ()
 {
-    Real small_temp_tmp = -1.e200;
-    Real small_dens_tmp = -1.e200;
+    amrex::Real small_temp_tmp = -1.e200;
+    amrex::Real small_dens_tmp = -1.e200;
     eos_init(small_temp_tmp, small_dens_tmp);
 }
 
diff --git a/interfaces/rhs_type.H b/interfaces/rhs_type.H
index 887b194b79..d51854400f 100644
--- a/interfaces/rhs_type.H
+++ b/interfaces/rhs_type.H
@@ -21,12 +21,12 @@ struct rhs_t {
     int species_E{-1};
     int species_F{-1};
 
-    Real number_A{0.0_rt};
-    Real number_B{0.0_rt};
-    Real number_C{0.0_rt};
-    Real number_D{0.0_rt};
-    Real number_E{0.0_rt};
-    Real number_F{0.0_rt};
+    amrex::Real number_A{0.0_rt};
+    amrex::Real number_B{0.0_rt};
+    amrex::Real number_C{0.0_rt};
+    amrex::Real number_D{0.0_rt};
+    amrex::Real number_E{0.0_rt};
+    amrex::Real number_F{0.0_rt};
 
     int exponent_A{0};
     int exponent_B{0};
@@ -55,9 +55,9 @@ constexpr amrex::Real tab_thi = 10.0e0_rt;
 constexpr int tab_per_decade = 2000;
 constexpr int nrattab = static_cast<int>(tab_thi - tab_tlo) * tab_per_decade + 1;
 constexpr int tab_imax = static_cast<int>(tab_thi - tab_tlo) * tab_per_decade + 1;
-constexpr amrex::Real tab_tstp = (tab_thi - tab_tlo) / static_cast<Real>(tab_imax - 1);
+constexpr amrex::Real tab_tstp = (tab_thi - tab_tlo) / static_cast<amrex::Real>(tab_imax - 1);
 
-extern AMREX_GPU_MANAGED Array1D<amrex::Real, 1, nrattab> ttab;
+extern AMREX_GPU_MANAGED amrex::Array1D<amrex::Real, 1, nrattab> ttab;
 
 // A struct that contains the terms needed to evaluate a tabulated rate.
 struct rate_tab_t
@@ -105,7 +105,7 @@ struct rhs_state_t
 #endif
     amrex::Real y_e;
     amrex::Real eta;
-    Array1D<amrex::Real, 1, NumSpec> y;
+    amrex::Array1D<amrex::Real, 1, NumSpec> y;
 };
 
 struct rate_t
diff --git a/interfaces/tfactors.H b/interfaces/tfactors.H
index 983b2d43c1..ba54d9eef2 100644
--- a/interfaces/tfactors.H
+++ b/interfaces/tfactors.H
@@ -4,7 +4,7 @@
 #include <AMReX.H>
 #include <cmath>
 
-using namespace amrex;
+using namespace amrex::literals;
 
 struct tf_t {
     amrex::Real temp;
diff --git a/networks/rhs.H b/networks/rhs.H
index 7c3acfb74a..01b721bbc3 100644
--- a/networks/rhs.H
+++ b/networks/rhs.H
@@ -27,8 +27,8 @@ namespace RHS
 {
 
 // Rate tabulation data.
-extern AMREX_GPU_MANAGED Array3D<Real, 1, Rates::NumRates, 1, 2, 1, nrattab> rattab;
-extern AMREX_GPU_MANAGED Array3D<Real, 1, Rates::NumRates, 1, 2, 1, nrattab> drattabdt;
+extern AMREX_GPU_MANAGED amrex::Array3D<amrex::Real, 1, Rates::NumRates, 1, 2, 1, nrattab> rattab;
+extern AMREX_GPU_MANAGED amrex::Array3D<amrex::Real, 1, Rates::NumRates, 1, 2, 1, nrattab> drattabdt;
 
 // Calculate an integer factorial.
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
@@ -273,14 +273,14 @@ constexpr int is_jacobian_term_used ()
 }
 
 typedef amrex::Array1D<int, 1, INT_NEQS> IArray1D;
-typedef amrex::Array1D<Real, 1, INT_NEQS> RArray1D;
+typedef amrex::Array1D<amrex::Real, 1, INT_NEQS> RArray1D;
 typedef ArrayUtil::MathArray2D<1, INT_NEQS, 1, INT_NEQS> RArray2D;
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
 void dgesl (RArray2D& a1, RArray1D& b1)
 {
-    auto const& a = reinterpret_cast<Array2D<Real, 0, INT_NEQS-1, 0, INT_NEQS-1>&>(a1);
-    auto& b = reinterpret_cast<Array1D<Real, 0, INT_NEQS-1>&>(b1);
+    auto const& a = reinterpret_cast<amrex::Array2D<amrex::Real, 0, INT_NEQS-1, 0, INT_NEQS-1>&>(a1);
+    auto& b = reinterpret_cast<amrex::Array1D<amrex::Real, 0, INT_NEQS-1>&>(b1);
 
     // solve a * x = b
     // first solve l * y = b
@@ -288,7 +288,7 @@ void dgesl (RArray2D& a1, RArray1D& b1)
     {
         constexpr int k = n1;
 
-        Real t = b(k);
+        amrex::Real t = b(k);
         constexpr_for<k+1, INT_NEQS>([&] (auto n2)
         {
             constexpr int j = n2;
@@ -305,7 +305,7 @@ void dgesl (RArray2D& a1, RArray1D& b1)
         constexpr int k = INT_NEQS - kb - 1;
 
         b(k) = b(k) / a(k,k);
-        Real t = -b(k);
+        amrex::Real t = -b(k);
 
         constexpr_for<0, k>([&] (auto j)
         {
@@ -319,7 +319,7 @@ void dgesl (RArray2D& a1, RArray1D& b1)
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
 void dgefa (RArray2D& a1)
 {
-    auto& a = reinterpret_cast<Array2D<Real, 0, INT_NEQS-1, 0, INT_NEQS-1>&>(a1);
+    auto& a = reinterpret_cast<amrex::Array2D<amrex::Real, 0, INT_NEQS-1, 0, INT_NEQS-1>&>(a1);
 
     // LU factorization in-place without pivoting.
 
@@ -328,7 +328,7 @@ void dgefa (RArray2D& a1)
         [[maybe_unused]] constexpr int k = n1;
 
         // compute multipliers
-        Real t = -1.0_rt / a(k,k);
+        amrex::Real t = -1.0_rt / a(k,k);
         constexpr_for<k+1, INT_NEQS>([&] (auto n2)
         {
             [[maybe_unused]] constexpr int j = n2;
@@ -422,8 +422,8 @@ void apply_density_scaling (const rhs_state_t& state, rate_t& rates)
     static_assert(reverse_exponent <= 3);
 
     // We know the exponent is an integer so we can construct the density term explicitly.
-    Real density_term_forward = 1.0;
-    Real density_term_reverse = 1.0;
+    amrex::Real density_term_forward = 1.0;
+    amrex::Real density_term_reverse = 1.0;
 
     if constexpr (forward_exponent == 1) {
         density_term_forward = state.rho;
@@ -469,11 +469,11 @@ void apply_screening (const rhs_state_t& state, rate_t& rates)
     if constexpr (data.exponent_A == 1 && data.exponent_B == 1 && data.exponent_C == 0) {
         // Forward reaction is A + B, screen using these two species
 
-        constexpr Real Z1 = NetworkProperties::zion(data.species_A);
-        constexpr Real A1 = NetworkProperties::aion(data.species_A);
+        constexpr amrex::Real Z1 = NetworkProperties::zion(data.species_A);
+        constexpr amrex::Real A1 = NetworkProperties::aion(data.species_A);
 
-        constexpr Real Z2 = NetworkProperties::zion(data.species_B);
-        constexpr Real A2 = NetworkProperties::aion(data.species_B);
+        constexpr amrex::Real Z2 = NetworkProperties::zion(data.species_B);
+        constexpr amrex::Real A2 = NetworkProperties::aion(data.species_B);
 
         constexpr auto scn_fac = scrn::calculate_screen_factor(Z1, A1, Z2, A2);
 
@@ -481,7 +481,7 @@ void apply_screening (const rhs_state_t& state, rate_t& rates)
         // compiler to evaluate the screen factor at compile time.
         static_assert(scn_fac.z1 == Z1);
 
-        Real sc, scdt;
+        amrex::Real sc, scdt;
         actual_screen(state.pstate, scn_fac, sc, scdt);
 
         if constexpr (data.screen_forward_reaction == 1) {
@@ -498,14 +498,14 @@ void apply_screening (const rhs_state_t& state, rate_t& rates)
     if constexpr (data.exponent_A == 2 && data.exponent_B == 0 && data.exponent_C == 0) {
         // Forward reaction is A + A, screen using just this species
 
-        constexpr Real Z1 = NetworkProperties::zion(data.species_A);
-        constexpr Real A1 = NetworkProperties::aion(data.species_A);
+        constexpr amrex::Real Z1 = NetworkProperties::zion(data.species_A);
+        constexpr amrex::Real A1 = NetworkProperties::aion(data.species_A);
 
         constexpr auto scn_fac = scrn::calculate_screen_factor(Z1, A1, Z1, A1);
 
         static_assert(scn_fac.z1 == Z1);
 
-        Real sc, scdt;
+        amrex::Real sc, scdt;
         actual_screen(state.pstate, scn_fac, sc, scdt);
 
         if constexpr (data.screen_forward_reaction == 1) {
@@ -523,30 +523,30 @@ void apply_screening (const rhs_state_t& state, rate_t& rates)
         // Forward reaction is triple alpha or an equivalent, screen using A + A
         // and then A + X where X has twice the number of protons and neutrons.
 
-        constexpr Real Z1 = NetworkProperties::zion(data.species_A);
-        constexpr Real A1 = NetworkProperties::aion(data.species_A);
+        constexpr amrex::Real Z1 = NetworkProperties::zion(data.species_A);
+        constexpr amrex::Real A1 = NetworkProperties::aion(data.species_A);
 
         constexpr auto scn_fac1 = scrn::calculate_screen_factor(Z1, A1, Z1, A1);
 
         static_assert(scn_fac1.z1 == Z1);
 
-        Real sc1, sc1dt;
+        amrex::Real sc1, sc1dt;
         actual_screen(state.pstate, scn_fac1, sc1, sc1dt);
 
-        constexpr Real Z2 = 2.0_rt * Z1;
-        constexpr Real A2 = 2.0_rt * A1;
+        constexpr amrex::Real Z2 = 2.0_rt * Z1;
+        constexpr amrex::Real A2 = 2.0_rt * A1;
 
         constexpr auto scn_fac2 = scrn::calculate_screen_factor(Z1, A1, Z2, A2);
 
         static_assert(scn_fac2.z1 == Z1);
 
-        Real sc2, sc2dt;
+        amrex::Real sc2, sc2dt;
         actual_screen(state.pstate, scn_fac2, sc2, sc2dt);
 
         // Compute combined screening factor
 
-        Real sc = sc1 * sc2;
-        Real scdt = sc1dt * sc2 + sc1 * sc2dt;
+        amrex::Real sc = sc1 * sc2;
+        amrex::Real scdt = sc1dt * sc2 + sc1 * sc2dt;
 
         if constexpr (data.screen_forward_reaction == 1) {
             rates.frdt = rates.frdt * sc + rates.fr * scdt;
@@ -587,7 +587,7 @@ void tabulate_rates ()
     using namespace Rates;
 
     for (int i = 1; i <= tab_imax; ++i) {
-        Real temp = tab_tlo + static_cast<Real>(i-1) * tab_tstp;
+        amrex::Real temp = tab_tlo + static_cast<amrex::Real>(i-1) * tab_tstp;
         temp = std::pow(10.0e0_rt, temp);
 
         ttab(i) = temp;
@@ -692,7 +692,7 @@ void evaluate_tabulated_rate (const rhs_state_t& state, rate_t& rates)
 // in the multiplication.
 template<int species, int rate, int use_T_derivatives>
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-constexpr std::pair<Real, Real> rhs_term (const burn_t& state, const rate_t& rates)
+constexpr std::pair<amrex::Real, amrex::Real> rhs_term (const burn_t& state, const rate_t& rates)
 {
     constexpr rhs_t data = rhs_data(rate);
 
@@ -700,14 +700,14 @@ constexpr std::pair<Real, Real> rhs_term (const burn_t& state, const rate_t& rat
     // reverse reactions, which is the same regardless of which species
     // we're producing or consuming.
 
-    Real forward_term = use_T_derivatives ? rates.frdt : rates.fr;
+    amrex::Real forward_term = use_T_derivatives ? rates.frdt : rates.fr;
 
     if constexpr (data.species_A >= 0) {
-        Real Y_A = state.xn[data.species_A-1] * aion_inv[data.species_A-1];
+        amrex::Real Y_A = state.xn[data.species_A-1] * aion_inv[data.species_A-1];
 
         static_assert(data.exponent_A <= 3);
 
-        Real df = 1.0;
+        amrex::Real df = 1.0;
 
         if constexpr (data.exponent_A == 1) {
             df = Y_A;
@@ -721,18 +721,18 @@ constexpr std::pair<Real, Real> rhs_term (const burn_t& state, const rate_t& rat
 
         if constexpr (data.apply_identical_particle_factor) {
             constexpr int identical_particle_factor = factorial(data.exponent_A);
-            df *= 1.0_rt / static_cast<Real>(identical_particle_factor);
+            df *= 1.0_rt / static_cast<amrex::Real>(identical_particle_factor);
         }
 
         forward_term *= df;
     }
 
     if constexpr (data.species_B >= 0) {
-        Real Y_B = state.xn[data.species_B-1] * aion_inv[data.species_B-1];
+        amrex::Real Y_B = state.xn[data.species_B-1] * aion_inv[data.species_B-1];
 
         static_assert(data.exponent_B <= 3);
 
-        Real df = 1.0;
+        amrex::Real df = 1.0;
 
         if constexpr (data.exponent_B == 1) {
             df = Y_B;
@@ -746,18 +746,18 @@ constexpr std::pair<Real, Real> rhs_term (const burn_t& state, const rate_t& rat
 
         if constexpr (data.apply_identical_particle_factor) {
             constexpr int identical_particle_factor = factorial(data.exponent_B);
-            df *= 1.0_rt / static_cast<Real>(identical_particle_factor);
+            df *= 1.0_rt / static_cast<amrex::Real>(identical_particle_factor);
         }
 
         forward_term *= df;
     }
 
     if constexpr (data.species_C >= 0) {
-        Real Y_C = state.xn[data.species_C-1] * aion_inv[data.species_C-1];
+        amrex::Real Y_C = state.xn[data.species_C-1] * aion_inv[data.species_C-1];
 
         static_assert(data.exponent_C <= 3);
 
-        Real df = 1.0;
+        amrex::Real df = 1.0;
 
         if constexpr (data.exponent_C == 1) {
             df = Y_C;
@@ -771,20 +771,20 @@ constexpr std::pair<Real, Real> rhs_term (const burn_t& state, const rate_t& rat
 
         if constexpr (data.apply_identical_particle_factor) {
             constexpr int identical_particle_factor = factorial(data.exponent_C);
-            df *= 1.0_rt / static_cast<Real>(identical_particle_factor);
+            df *= 1.0_rt / static_cast<amrex::Real>(identical_particle_factor);
         }
 
         forward_term *= df;
     }
 
-    Real reverse_term = use_T_derivatives ? rates.rrdt : rates.rr;
+    amrex::Real reverse_term = use_T_derivatives ? rates.rrdt : rates.rr;
 
     if constexpr (data.species_D >= 0) {
-        Real Y_D = state.xn[data.species_D-1] * aion_inv[data.species_D-1];
+        amrex::Real Y_D = state.xn[data.species_D-1] * aion_inv[data.species_D-1];
 
         static_assert(data.exponent_D <= 3);
 
-        Real dr = 1.0;
+        amrex::Real dr = 1.0;
 
         if constexpr (data.exponent_D == 1) {
             dr = Y_D;
@@ -798,16 +798,16 @@ constexpr std::pair<Real, Real> rhs_term (const burn_t& state, const rate_t& rat
 
         if constexpr (data.apply_identical_particle_factor) {
             constexpr int identical_particle_factor = factorial(data.exponent_D);
-            dr *= 1.0_rt / static_cast<Real>(identical_particle_factor);
+            dr *= 1.0_rt / static_cast<amrex::Real>(identical_particle_factor);
         }
 
         reverse_term *= dr;
     }
 
     if constexpr (data.species_E >= 0) {
-        Real Y_E = state.xn[data.species_E-1] * aion_inv[data.species_E-1];
+        amrex::Real Y_E = state.xn[data.species_E-1] * aion_inv[data.species_E-1];
 
-        Real dr = 1.0;
+        amrex::Real dr = 1.0;
 
         if constexpr (data.exponent_E == 1) {
             dr = Y_E;
@@ -821,16 +821,16 @@ constexpr std::pair<Real, Real> rhs_term (const burn_t& state, const rate_t& rat
 
         if constexpr (data.apply_identical_particle_factor) {
             constexpr int identical_particle_factor = factorial(data.exponent_E);
-            dr *= 1.0_rt / static_cast<Real>(identical_particle_factor);
+            dr *= 1.0_rt / static_cast<amrex::Real>(identical_particle_factor);
         }
 
         reverse_term *= dr;
     }
 
     if constexpr (data.species_F >= 0) {
-        Real Y_F = state.xn[data.species_F-1] * aion_inv[data.species_F-1];
+        amrex::Real Y_F = state.xn[data.species_F-1] * aion_inv[data.species_F-1];
 
-        Real dr = 1.0;
+        amrex::Real dr = 1.0;
 
         if constexpr (data.exponent_F == 1) {
             dr = Y_F;
@@ -844,7 +844,7 @@ constexpr std::pair<Real, Real> rhs_term (const burn_t& state, const rate_t& rat
 
         if constexpr (data.apply_identical_particle_factor) {
             constexpr int identical_particle_factor = factorial(data.exponent_F);
-            dr *= 1.0_rt / static_cast<Real>(identical_particle_factor);
+            dr *= 1.0_rt / static_cast<amrex::Real>(identical_particle_factor);
         }
 
         reverse_term *= dr;
@@ -895,11 +895,11 @@ constexpr std::pair<Real, Real> rhs_term (const burn_t& state, const rate_t& rat
 // contribution if spec2 is one of (D, E, F).
 template<int spec1, int spec2, int rate>
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-constexpr Real jac_term (const burn_t& state, const rate_t& rates)
+constexpr amrex::Real jac_term (const burn_t& state, const rate_t& rates)
 {
     constexpr rhs_t data = rhs_data(rate);
 
-    Real forward_term = 0.0_rt;
+    amrex::Real forward_term = 0.0_rt;
 
     if constexpr (is_rate_used<spec1, rate>() &&
                   (spec2 == data.species_A || spec2 == data.species_B || spec2 == data.species_C)) {
@@ -907,7 +907,7 @@ constexpr Real jac_term (const burn_t& state, const rate_t& rates)
         forward_term = rates.fr;
 
         if constexpr (data.species_A >= 0) {
-            Real Y_A = state.xn[data.species_A-1] * aion_inv[data.species_A-1];
+            amrex::Real Y_A = state.xn[data.species_A-1] * aion_inv[data.species_A-1];
 
             constexpr int exponent = data.exponent_A;
 
@@ -921,7 +921,7 @@ constexpr Real jac_term (const burn_t& state, const rate_t& rates)
 
             static_assert(exp <= 3);
 
-            Real df = 1.0;
+            amrex::Real df = 1.0;
 
             if constexpr (exp == 1) {
                 df = Y_A;
@@ -939,14 +939,14 @@ constexpr Real jac_term (const burn_t& state, const rate_t& rates)
 
             if constexpr (data.apply_identical_particle_factor) {
                 constexpr int identical_particle_factor = factorial(exponent);
-                df *= 1.0_rt / static_cast<Real>(identical_particle_factor);
+                df *= 1.0_rt / static_cast<amrex::Real>(identical_particle_factor);
             }
 
             forward_term *= df;
         }
 
         if constexpr (data.species_B >= 0) {
-            Real Y_B = state.xn[data.species_B-1] * aion_inv[data.species_B-1];
+            amrex::Real Y_B = state.xn[data.species_B-1] * aion_inv[data.species_B-1];
 
             constexpr int exponent = data.exponent_B;
 
@@ -954,7 +954,7 @@ constexpr Real jac_term (const burn_t& state, const rate_t& rates)
 
             static_assert(exp <= 3);
 
-            Real df = 1.0;
+            amrex::Real df = 1.0;
 
             if constexpr (exp == 1) {
                 df = Y_B;
@@ -972,14 +972,14 @@ constexpr Real jac_term (const burn_t& state, const rate_t& rates)
 
             if constexpr (data.apply_identical_particle_factor) {
                 constexpr int identical_particle_factor = factorial(exponent);
-                df *= 1.0_rt / static_cast<Real>(identical_particle_factor);
+                df *= 1.0_rt / static_cast<amrex::Real>(identical_particle_factor);
             }
 
             forward_term *= df;
         }
 
         if constexpr (data.species_C >= 0) {
-            Real Y_C = state.xn[data.species_C-1] * aion_inv[data.species_C-1];
+            amrex::Real Y_C = state.xn[data.species_C-1] * aion_inv[data.species_C-1];
 
             constexpr int exponent = data.exponent_C;
 
@@ -987,7 +987,7 @@ constexpr Real jac_term (const burn_t& state, const rate_t& rates)
 
             static_assert(exp <= 3);
 
-            Real df = 1.0;
+            amrex::Real df = 1.0;
 
             if constexpr (exp == 1) {
                 df = Y_C;
@@ -1005,7 +1005,7 @@ constexpr Real jac_term (const burn_t& state, const rate_t& rates)
 
             if constexpr (data.apply_identical_particle_factor) {
                 constexpr int identical_particle_factor = factorial(exponent);
-                df *= 1.0_rt / static_cast<Real>(identical_particle_factor);
+                df *= 1.0_rt / static_cast<amrex::Real>(identical_particle_factor);
             }
 
             forward_term *= df;
@@ -1013,7 +1013,7 @@ constexpr Real jac_term (const burn_t& state, const rate_t& rates)
 
     }
 
-    Real reverse_term = 0.0_rt;
+    amrex::Real reverse_term = 0.0_rt;
 
     if constexpr (is_rate_used<spec1, rate>() &&
                   (spec2 == data.species_D || spec2 == data.species_E || spec2 == data.species_F)) {
@@ -1021,7 +1021,7 @@ constexpr Real jac_term (const burn_t& state, const rate_t& rates)
         reverse_term = rates.rr;
 
         if constexpr (data.species_D >= 0) {
-            Real Y_D = state.xn[data.species_D-1] * aion_inv[data.species_D-1];
+            amrex::Real Y_D = state.xn[data.species_D-1] * aion_inv[data.species_D-1];
 
             constexpr int exponent = data.exponent_D;
 
@@ -1029,7 +1029,7 @@ constexpr Real jac_term (const burn_t& state, const rate_t& rates)
 
             static_assert(exp <= 3);
 
-            Real dr = 1.0;
+            amrex::Real dr = 1.0;
 
             if constexpr (exp == 1) {
                 dr = Y_D;
@@ -1047,14 +1047,14 @@ constexpr Real jac_term (const burn_t& state, const rate_t& rates)
 
             if constexpr (data.apply_identical_particle_factor) {
                 constexpr int identical_particle_factor = factorial(exponent);
-                dr *= 1.0_rt / static_cast<Real>(identical_particle_factor);
+                dr *= 1.0_rt / static_cast<amrex::Real>(identical_particle_factor);
             }
 
             reverse_term *= dr;
         }
 
         if constexpr (data.species_E >= 0) {
-            Real Y_E = state.xn[data.species_E-1] * aion_inv[data.species_E-1];
+            amrex::Real Y_E = state.xn[data.species_E-1] * aion_inv[data.species_E-1];
 
             constexpr int exponent = data.exponent_E;
 
@@ -1062,7 +1062,7 @@ constexpr Real jac_term (const burn_t& state, const rate_t& rates)
 
             static_assert(exp <= 3);
 
-            Real dr = 1.0;
+            amrex::Real dr = 1.0;
 
             if constexpr (exp == 1) {
                 dr = Y_E;
@@ -1080,14 +1080,14 @@ constexpr Real jac_term (const burn_t& state, const rate_t& rates)
 
             if constexpr (data.apply_identical_particle_factor) {
                 constexpr int identical_particle_factor = factorial(exponent);
-                dr *= 1.0_rt / static_cast<Real>(identical_particle_factor);
+                dr *= 1.0_rt / static_cast<amrex::Real>(identical_particle_factor);
             }
 
             reverse_term *= dr;
         }
 
         if constexpr (data.species_F >= 0) {
-            Real Y_F = state.xn[data.species_F-1] * aion_inv[data.species_F-1];
+            amrex::Real Y_F = state.xn[data.species_F-1] * aion_inv[data.species_F-1];
 
             constexpr int exponent = data.exponent_F;
 
@@ -1095,7 +1095,7 @@ constexpr Real jac_term (const burn_t& state, const rate_t& rates)
 
             static_assert(exp <= 3);
 
-            Real dr = 1.0;
+            amrex::Real dr = 1.0;
 
             if constexpr (exp == 1) {
                 dr = Y_F;
@@ -1113,7 +1113,7 @@ constexpr Real jac_term (const burn_t& state, const rate_t& rates)
 
             if constexpr (data.apply_identical_particle_factor) {
                 constexpr int identical_particle_factor = factorial(exponent);
-                dr *= 1.0_rt / static_cast<Real>(identical_particle_factor);
+                dr *= 1.0_rt / static_cast<amrex::Real>(identical_particle_factor);
             }
 
             reverse_term *= dr;
@@ -1123,7 +1123,7 @@ constexpr Real jac_term (const burn_t& state, const rate_t& rates)
 
     // Now compute the total contribution to this species.
 
-    Real term = 0.0_rt;
+    amrex::Real term = 0.0_rt;
 
     if constexpr (data.species_A == spec1) {
         term = data.number_A * (reverse_term - forward_term);
@@ -1238,7 +1238,7 @@ void rhs_init ()
 // for each term in ydot).
 template<int nrhs>
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-void rhs (burn_t& burn_state, Array1D<Real, 1, nrhs>& ydot)
+void rhs (burn_t& burn_state, amrex::Array1D<amrex::Real, 1, nrhs>& ydot)
 {
     static_assert(nrhs == neqs || nrhs == 2 * neqs);
 
@@ -1276,7 +1276,7 @@ void rhs (burn_t& burn_state, Array1D<Real, 1, nrhs>& ydot)
     constexpr int intermediate_array_size = num_intermediate > 0 ? num_intermediate : 1;
 
     // Define forward and reverse (and d/dT) rate arrays.
-    Array1D<rate_t, 1, intermediate_array_size> intermediate_rates;
+    amrex::Array1D<rate_t, 1, intermediate_array_size> intermediate_rates;
 
     // Fill all intermediate rates first.
     constexpr_for<1, Rates::NumRates+1>([&] (auto n)
@@ -1350,11 +1350,11 @@ void rhs (burn_t& burn_state, Array1D<Real, 1, nrhs>& ydot)
     // Evaluate the neutrino cooling.
 #ifdef NEUTRINOS
     constexpr int do_derivatives{0};
-    Real sneut, dsneutdt, dsneutdd, dsnuda, dsnudz;
+    amrex::Real sneut, dsneutdt, dsneutdd, dsnuda, dsnudz;
     sneut5<do_derivatives>(burn_state.T, burn_state.rho, burn_state.abar, burn_state.zbar,
                            sneut, dsneutdt, dsneutdd, dsnuda, dsnudz);
 #else
-    Real sneut = 0.0;
+    amrex::Real sneut = 0.0;
 #endif
 
     // Compute the energy RHS term.
@@ -1420,7 +1420,7 @@ void jac (burn_t& burn_state, ArrayUtil::MathArray2D<1, neqs, 1, neqs>& jac)
     constexpr int intermediate_array_size = num_intermediate > 0 ? num_intermediate : 1;
 
     // Define forward and reverse (and d/dT) rate arrays.
-    Array1D<rate_t, 1, intermediate_array_size> intermediate_rates;
+    amrex::Array1D<rate_t, 1, intermediate_array_size> intermediate_rates;
 
     rate_t rates1, rates2, rates3;
 
@@ -1496,12 +1496,12 @@ void jac (burn_t& burn_state, ArrayUtil::MathArray2D<1, neqs, 1, neqs>& jac)
 
     // Evaluate the neutrino cooling.
 #ifdef NEUTRINOS
-    Real sneut, dsneutdt, dsneutdd, dsnuda, dsnudz;
+    amrex::Real sneut, dsneutdt, dsneutdd, dsnuda, dsnudz;
     constexpr int do_derivatives{1};
     sneut5<do_derivatives>(burn_state.T, burn_state.rho, burn_state.abar, burn_state.zbar,
                            sneut, dsneutdt, dsneutdd, dsnuda, dsnudz);
 #else
-    Real sneut = 0.0, dsneutdt = 0.0, dsneutdd = 0.0, dsnuda = 0.0, dsnudz = 0.0;
+    amrex::Real sneut = 0.0, dsneutdt = 0.0, dsneutdd = 0.0, dsnuda = 0.0, dsnudz = 0.0;
     amrex::ignore_unused(sneut, dsneutdd);
 #endif
 
@@ -1512,7 +1512,7 @@ void jac (burn_t& burn_state, ArrayUtil::MathArray2D<1, neqs, 1, neqs>& jac)
         constexpr int species = j;
 
         // Energy generation rate Jacobian elements with respect to species.
-        Real b1 = (-burn_state.abar * burn_state.abar * dsnuda + (NetworkProperties::zion(species) - burn_state.zbar) * burn_state.abar * dsnudz);
+        amrex::Real b1 = (-burn_state.abar * burn_state.abar * dsnuda + (NetworkProperties::zion(species) - burn_state.zbar) * burn_state.abar * dsnudz);
         jac(net_ienuc, species) = -b1;
 
         constexpr_for<1, NumSpec+1>([&] (auto i)
@@ -1536,7 +1536,7 @@ void jac (burn_t& burn_state, ArrayUtil::MathArray2D<1, neqs, 1, neqs>& jac)
 // namespace. This should be retired later when nothing still depends on those names.
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-void actual_rhs (burn_t& state, Array1D<Real, 1, neqs>& ydot)
+void actual_rhs (burn_t& state, amrex::Array1D<amrex::Real, 1, neqs>& ydot)
 {
     RHS::rhs(state, ydot);
 }
diff --git a/networks/rhs.cpp b/networks/rhs.cpp
index 269e2486f4..b02f935d53 100644
--- a/networks/rhs.cpp
+++ b/networks/rhs.cpp
@@ -2,8 +2,8 @@
 
 #ifdef NEW_NETWORK_IMPLEMENTATION
 
-AMREX_GPU_MANAGED Array3D<Real, 1, Rates::NumRates, 1, 2, 1, RHS::nrattab> RHS::rattab;
-AMREX_GPU_MANAGED Array3D<Real, 1, Rates::NumRates, 1, 2, 1, RHS::nrattab> RHS::drattabdt;
-AMREX_GPU_MANAGED Array1D<Real, 1, RHS::nrattab> RHS::ttab;
+AMREX_GPU_MANAGED Array3D<amrex::Real, 1, Rates::NumRates, 1, 2, 1, RHS::nrattab> RHS::rattab;
+AMREX_GPU_MANAGED Array3D<amrex::Real, 1, Rates::NumRates, 1, 2, 1, RHS::nrattab> RHS::drattabdt;
+AMREX_GPU_MANAGED Array1D<amrex::Real, 1, RHS::nrattab> RHS::ttab;
 
 #endif
diff --git a/networks/rprox/rprox_rates.H b/networks/rprox/rprox_rates.H
index 1d8dbe5f9e..7d9394ac5f 100644
--- a/networks/rprox/rprox_rates.H
+++ b/networks/rprox/rprox_rates.H
@@ -5,25 +5,25 @@
 #include <tfactors.H>
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-void rate_p_c12_to_n13 (const tf_t& tfactors, Real& rate, Real& dratedt)
+void rate_p_c12_to_n13 (const tf_t& tfactors, amrex::Real& rate, amrex::Real& dratedt)
 {
     rate = 0.0_rt;
     dratedt = 0.0_rt;
 
-    Real ct9i = tfactors.t9i * (-3.77849e0_rt);
-    Real ct9i13 = tfactors.t9i13 * (-5.10735e0_rt);
-    Real ct913 = tfactors.t913 * (-2.24111e0_rt);
-    Real ct9 = tfactors.t9 * (0.148883e0_rt);
-    Real clnt9 = tfactors.lnt9 * (-1.5e0_rt);
+    amrex::Real ct9i = tfactors.t9i * (-3.77849e0_rt);
+    amrex::Real ct9i13 = tfactors.t9i13 * (-5.10735e0_rt);
+    amrex::Real ct913 = tfactors.t913 * (-2.24111e0_rt);
+    amrex::Real ct9 = tfactors.t9 * (0.148883e0_rt);
+    amrex::Real clnt9 = tfactors.lnt9 * (-1.5e0_rt);
 
-    Real r0 = std::exp(17.5428e0_rt +
+    amrex::Real r0 = std::exp(17.5428e0_rt +
                        ct9i +
                        ct9i13 +
                        ct913 +
                        ct9 +
                        clnt9);
 
-    Real dr0dt = r0 * tfactors.t9i * (-ct9i -
+    amrex::Real dr0dt = r0 * tfactors.t9i * (-ct9i -
                                       (1.0_rt / 3.0_rt) * ct9i13 +
                                       (1.0_rt / 3.0_rt) * ct913 +
                                       ct9 +
@@ -32,17 +32,17 @@ void rate_p_c12_to_n13 (const tf_t& tfactors, Real& rate, Real& dratedt)
     ct9i13 = tfactors.t9i13 * (-13.692e0_rt);
     ct913 = tfactors.t913 * (-0.230881e0_rt);
     ct9 = tfactors.t9 * (4.44362e0_rt);
-    Real ct953 = tfactors.t953 * (-3.15898e0_rt);
+    amrex::Real ct953 = tfactors.t953 * (-3.15898e0_rt);
     clnt9 = tfactors.lnt9 * (-0.666667e0_rt);
 
-    Real r1 = std::exp(17.1482e0_rt +
+    amrex::Real r1 = std::exp(17.1482e0_rt +
                        ct9i13 +
                        ct913 +
                        ct9 +
                        ct953 +
                        clnt9);
 
-    Real dr1dt = r1 * tfactors.t9i * (-(1.0_rt / 3.0_rt) * ct9i13 +
+    amrex::Real dr1dt = r1 * tfactors.t9i * (-(1.0_rt / 3.0_rt) * ct9i13 +
                                       (1.0_rt / 3.0_rt) * ct913 +
                                       ct9 +
                                       (5.0_rt / 3.0_rt) * ct953 +
@@ -54,19 +54,19 @@ void rate_p_c12_to_n13 (const tf_t& tfactors, Real& rate, Real& dratedt)
 
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-void rate_f17_to_p_o16 (const tf_t& tfactors, Real& rate, Real& dratedt)
+void rate_f17_to_p_o16 (const tf_t& tfactors, amrex::Real& rate, amrex::Real& dratedt)
 {
     rate = 0.0_rt;
     dratedt = 0.0_rt;
 
-    Real ct9i = tfactors.t9i * (-6.96583e0_rt);
-    Real ct9i13 = tfactors.t9i13 * (-16.696e0_rt);
-    Real ct913 = tfactors.t913 * (-1.16252e0_rt);
-    Real ct9 = tfactors.t9 * (0.267703e0_rt);
-    Real ct953 = tfactors.t953 * (-0.0338411e0_rt);
-    Real clnt9 = tfactors.lnt9 * (0.833333e0_rt);
+    amrex::Real ct9i = tfactors.t9i * (-6.96583e0_rt);
+    amrex::Real ct9i13 = tfactors.t9i13 * (-16.696e0_rt);
+    amrex::Real ct913 = tfactors.t913 * (-1.16252e0_rt);
+    amrex::Real ct9 = tfactors.t9 * (0.267703e0_rt);
+    amrex::Real ct953 = tfactors.t953 * (-0.0338411e0_rt);
+    amrex::Real clnt9 = tfactors.lnt9 * (0.833333e0_rt);
 
-    Real r0 = std::exp(40.9135e0_rt +
+    amrex::Real r0 = std::exp(40.9135e0_rt +
                        ct9i +
                        ct9i13 +
                        ct913 +
@@ -74,7 +74,7 @@ void rate_f17_to_p_o16 (const tf_t& tfactors, Real& rate, Real& dratedt)
                        ct953 +
                        clnt9);
 
-    Real dr0dt = r0 * tfactors.t9i * (-ct9i -
+    amrex::Real dr0dt = r0 * tfactors.t9i * (-ct9i -
                                       (1.0_rt / 3.0_rt) * ct9i13 +
                                       (1.0_rt / 3.0_rt) * ct913 +
                                       ct9 +
@@ -87,7 +87,7 @@ void rate_f17_to_p_o16 (const tf_t& tfactors, Real& rate, Real& dratedt)
 
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-void rate_f17_to_o17 ([[maybe_unused]] const tf_t& tfactors, Real& rate, Real& dratedt)
+void rate_f17_to_o17 ([[maybe_unused]] const tf_t& tfactors, amrex::Real& rate, amrex::Real& dratedt)
 {
     rate = std::exp(-4.53302e0_rt);
     dratedt = 0.0_rt;
@@ -95,19 +95,19 @@ void rate_f17_to_o17 ([[maybe_unused]] const tf_t& tfactors, Real& rate, Real& d
 
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-void rate_p_f17_to_ne18 (const tf_t& tfactors, Real& rate, Real& dratedt)
+void rate_p_f17_to_ne18 (const tf_t& tfactors, amrex::Real& rate, amrex::Real& dratedt)
 {
     rate = 0.0_rt;
     dratedt = 0.0_rt;
 
-    Real ct9i = tfactors.t9i * (-0.0323504e0_rt);
-    Real ct9i13 = tfactors.t9i13 * (-14.2191e0_rt);
-    Real ct913 = tfactors.t913 * (34.0647e0_rt);
-    Real ct9 = tfactors.t9 * (-16.5698e0_rt);
-    Real ct953 = tfactors.t953 * (2.48116e0_rt);
-    Real clnt9 = tfactors.lnt9 * (-2.13376e0_rt);
+    amrex::Real ct9i = tfactors.t9i * (-0.0323504e0_rt);
+    amrex::Real ct9i13 = tfactors.t9i13 * (-14.2191e0_rt);
+    amrex::Real ct913 = tfactors.t913 * (34.0647e0_rt);
+    amrex::Real ct9 = tfactors.t9 * (-16.5698e0_rt);
+    amrex::Real ct953 = tfactors.t953 * (2.48116e0_rt);
+    amrex::Real clnt9 = tfactors.lnt9 * (-2.13376e0_rt);
 
-    Real r0 = std::exp(-7.84708e0_rt +
+    amrex::Real r0 = std::exp(-7.84708e0_rt +
                        ct9i +
                        ct9i13 +
                        ct913 +
@@ -115,7 +115,7 @@ void rate_p_f17_to_ne18 (const tf_t& tfactors, Real& rate, Real& dratedt)
                        ct953 +
                        clnt9);
 
-    Real dr0dt = r0 * tfactors.t9i * (-ct9i -
+    amrex::Real dr0dt = r0 * tfactors.t9i * (-ct9i -
                                       (1.0_rt / 3.0_rt) * ct9i13 +
                                       (1.0_rt / 3.0_rt) * ct913 +
                                       ct9 +
@@ -129,14 +129,14 @@ void rate_p_f17_to_ne18 (const tf_t& tfactors, Real& rate, Real& dratedt)
     ct953 = tfactors.t953 * (-0.0440377e0_rt);
     clnt9 = tfactors.lnt9 * (-7.36014e0_rt);
 
-    Real r1 = std::exp(27.5778e0_rt +
+    amrex::Real r1 = std::exp(27.5778e0_rt +
                        ct9i +
                        ct9i13 +
                        ct913 +
                        ct9 +
                        ct953 +
                        clnt9);
-    Real dr1dt = r1 * tfactors.t9i * (-ct9i -
+    amrex::Real dr1dt = r1 * tfactors.t9i * (-ct9i -
                                       (1.0_rt / 3.0_rt) * ct9i13 +
                                       (1.0_rt / 3.0_rt) * ct913 +
                                       ct9 +
@@ -149,19 +149,19 @@ void rate_p_f17_to_ne18 (const tf_t& tfactors, Real& rate, Real& dratedt)
 
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-void rate_he4_he4_he4_to_c12 (const tf_t& tfactors, Real& rate, Real& dratedt)
+void rate_he4_he4_he4_to_c12 (const tf_t& tfactors, amrex::Real& rate, amrex::Real& dratedt)
 {
     rate = 0.0_rt;
     dratedt = 0.0_rt;
 
-    Real ct9i = tfactors.t9i * (-1.02446e0_rt);
-    Real ct9i13 = tfactors.t9i13 * (-23.57e0_rt);
-    Real ct913 = tfactors.t913 * (20.4886e0_rt);
-    Real ct9 = tfactors.t9 * (-12.9882e0_rt);
-    Real ct953 = tfactors.t953 * (-20.0e0_rt);
-    Real clnt9 = tfactors.lnt9 * (-2.16667e0_rt);
+    amrex::Real ct9i = tfactors.t9i * (-1.02446e0_rt);
+    amrex::Real ct9i13 = tfactors.t9i13 * (-23.57e0_rt);
+    amrex::Real ct913 = tfactors.t913 * (20.4886e0_rt);
+    amrex::Real ct9 = tfactors.t9 * (-12.9882e0_rt);
+    amrex::Real ct953 = tfactors.t953 * (-20.0e0_rt);
+    amrex::Real clnt9 = tfactors.lnt9 * (-2.16667e0_rt);
 
-    Real r0 = std::exp(-11.7884e0_rt +
+    amrex::Real r0 = std::exp(-11.7884e0_rt +
                        ct9i +
                        ct9i13 +
                        ct913 +
@@ -169,7 +169,7 @@ void rate_he4_he4_he4_to_c12 (const tf_t& tfactors, Real& rate, Real& dratedt)
                        ct953 +
                        clnt9);
 
-    Real dr0dt = r0 * tfactors.t9i * (-ct9i -
+    amrex::Real dr0dt = r0 * tfactors.t9i * (-ct9i -
                                       (1.0_rt / 3.0_rt) * ct9i13 +
                                       (1.0_rt / 3.0_rt) * ct913 +
                                       ct9 +
@@ -182,14 +182,14 @@ void rate_he4_he4_he4_to_c12 (const tf_t& tfactors, Real& rate, Real& dratedt)
     ct953 = tfactors.t953 * (-10.0e0_rt);
     clnt9 = tfactors.lnt9 * (-1.33333e0_rt);
 
-    Real r1 = std::exp(-0.971052e0_rt +
+    amrex::Real r1 = std::exp(-0.971052e0_rt +
                        ct9i13 +
                        ct913 +
                        ct9 +
                        ct953 +
                        clnt9);
 
-    Real dr1dt = r1 * tfactors.t9i * (-(1.0_rt / 3.0_rt) * ct9i13 +
+    amrex::Real dr1dt = r1 * tfactors.t9i * (-(1.0_rt / 3.0_rt) * ct9i13 +
                                       (1.0_rt / 3.0_rt) * ct913 +
                                       ct9 +
                                       (5.0_rt / 3.0_rt) * ct953 +
@@ -202,7 +202,7 @@ void rate_he4_he4_he4_to_c12 (const tf_t& tfactors, Real& rate, Real& dratedt)
     ct953 = tfactors.t953 * (0.0879816e0_rt);
     clnt9 = tfactors.lnt9 * (-13.1653e0_rt);
 
-    Real r2 = std::exp(-24.3505e0_rt +
+    amrex::Real r2 = std::exp(-24.3505e0_rt +
                        ct9i +
                        ct9i13 +
                        ct913 +
@@ -210,7 +210,7 @@ void rate_he4_he4_he4_to_c12 (const tf_t& tfactors, Real& rate, Real& dratedt)
                        ct953 +
                        clnt9);
 
-    Real dr2dt = r2 * tfactors.t9i * (-ct9i -
+    amrex::Real dr2dt = r2 * tfactors.t9i * (-ct9i -
                                       (1.0_rt / 3.0_rt) * ct9i13 +
                                       (1.0_rt / 3.0_rt) * ct913 +
                                       ct9 +
@@ -223,25 +223,25 @@ void rate_he4_he4_he4_to_c12 (const tf_t& tfactors, Real& rate, Real& dratedt)
 
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-void rate_p_n14_to_o15 (const tf_t& tfactors, Real& rate, Real& dratedt)
+void rate_p_n14_to_o15 (const tf_t& tfactors, amrex::Real& rate, amrex::Real& dratedt)
 {
     rate = 0.0_rt;
     dratedt = 0.0_rt;
 
-    Real ct9i13 = tfactors.t9i13 * (-15.193e0_rt);
-    Real ct913 = tfactors.t913 * (-4.63975e0_rt);
-    Real ct9 = tfactors.t9 * (9.73458e0_rt);
-    Real ct953 = tfactors.t953 * (-9.55051e0_rt);
-    Real clnt9 = tfactors.lnt9 * (0.333333e0_rt);
+    amrex::Real ct9i13 = tfactors.t9i13 * (-15.193e0_rt);
+    amrex::Real ct913 = tfactors.t913 * (-4.63975e0_rt);
+    amrex::Real ct9 = tfactors.t9 * (9.73458e0_rt);
+    amrex::Real ct953 = tfactors.t953 * (-9.55051e0_rt);
+    amrex::Real clnt9 = tfactors.lnt9 * (0.333333e0_rt);
 
-    Real r0 = std::exp(20.1169e0_rt +
+    amrex::Real r0 = std::exp(20.1169e0_rt +
                        ct9i13 +
                        ct913 +
                        ct9 +
                        ct953 +
                        clnt9);
 
-    Real dr0dt = r0 * tfactors.t9i * (-(1.0_rt / 3.0_rt) * ct9i13 +
+    amrex::Real dr0dt = r0 * tfactors.t9i * (-(1.0_rt / 3.0_rt) * ct9i13 +
                                       (1.0_rt / 3.0_rt) * ct913 +
                                       ct9 +
                                       (5.0_rt / 3.0_rt) * ct953 +
@@ -253,37 +253,37 @@ void rate_p_n14_to_o15 (const tf_t& tfactors, Real& rate, Real& dratedt)
     ct953 = tfactors.t953 * (-0.987565e0_rt);
     clnt9 = tfactors.lnt9 * (-0.666667e0_rt);
 
-    Real r1 = std::exp(17.01e0_rt +
+    amrex::Real r1 = std::exp(17.01e0_rt +
                        ct9i13 +
                        ct913 +
                        ct9 +
                        ct953 +
                        clnt9);
 
-    Real dr1dt = r1 * tfactors.t9i * (-(1.0_rt / 3.0_rt) * ct9i13 +
+    amrex::Real dr1dt = r1 * tfactors.t9i * (-(1.0_rt / 3.0_rt) * ct9i13 +
                                       (1.0_rt / 3.0_rt) * ct913 +
                                       ct9 +
                                       (5.0_rt / 3.0_rt) * ct953 +
                                       (-0.666667e0_rt));
 
-    Real ct9i = tfactors.t9i * (-4.891e0_rt);
+    amrex::Real ct9i = tfactors.t9i * (-4.891e0_rt);
     clnt9 = tfactors.lnt9 * (0.0682e0_rt);
 
-    Real r2 = std::exp(6.73578e0_rt +
+    amrex::Real r2 = std::exp(6.73578e0_rt +
                        ct9i +
                        clnt9);
 
-    Real dr2dt = r2 * tfactors.t9i * (-ct9i +
+    amrex::Real dr2dt = r2 * tfactors.t9i * (-ct9i +
                                       (0.0682e0_rt));
 
     ct9i = tfactors.t9i * (-2.998e0_rt);
     clnt9 = tfactors.lnt9 * (-1.5e0_rt);
 
-    Real r3 = std::exp(7.65444e0_rt +
+    amrex::Real r3 = std::exp(7.65444e0_rt +
                        ct9i +
                        clnt9);
 
-    Real dr3dt = r3 * tfactors.t9i * (-ct9i +
+    amrex::Real dr3dt = r3 * tfactors.t9i * (-ct9i +
                                       (-1.5e0_rt));
 
     rate = r0 + r1 + r2 + r3;
@@ -292,69 +292,69 @@ void rate_p_n14_to_o15 (const tf_t& tfactors, Real& rate, Real& dratedt)
 
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-void rate_he4_ne18_to_p_na21 (const tf_t& tfactors, Real& rate, Real& dratedt)
+void rate_he4_ne18_to_p_na21 (const tf_t& tfactors, amrex::Real& rate, amrex::Real& dratedt)
 {
     rate = 0.0_rt;
     dratedt = 0.0_rt;
 
-    Real ct9i = tfactors.t9i * (-24.7176e0_rt);
-    Real clnt9 = tfactors.lnt9 * (-1.5e0_rt);
+    amrex::Real ct9i = tfactors.t9i * (-24.7176e0_rt);
+    amrex::Real clnt9 = tfactors.lnt9 * (-1.5e0_rt);
 
-    Real r0 = std::exp(19.4058e0_rt +
+    amrex::Real r0 = std::exp(19.4058e0_rt +
                        ct9i +
                        clnt9);
 
-    Real dr0dt = r0 * tfactors.t9i * (-ct9i +
+    amrex::Real dr0dt = r0 * tfactors.t9i * (-ct9i +
                                       (-1.5e0_rt));
 
     ct9i = tfactors.t9i * (-0.452576e0_rt);
     clnt9 = tfactors.lnt9 * (-1.5e0_rt);
 
-    Real r1 = std::exp(-137.358e0_rt +
+    amrex::Real r1 = std::exp(-137.358e0_rt +
                        ct9i +
                        clnt9);
 
-    Real dr1dt = r1 * tfactors.t9i * (-ct9i +
+    amrex::Real dr1dt = r1 * tfactors.t9i * (-ct9i +
                                       (-1.5e0_rt));
 
     ct9i = tfactors.t9i * (-10.885e0_rt);
     clnt9 = tfactors.lnt9 * (-1.5e0_rt);
 
-    Real r2 = std::exp(6.39797e0_rt +
+    amrex::Real r2 = std::exp(6.39797e0_rt +
                        ct9i +
                        clnt9);
 
-    Real dr2dt = r2 * tfactors.t9i * (-ct9i +
+    amrex::Real dr2dt = r2 * tfactors.t9i * (-ct9i +
                                       (-1.5e0_rt));
 
     ct9i = tfactors.t9i * (-7.45009e0_rt);
     clnt9 = tfactors.lnt9 * (-1.5e0_rt);
 
-    Real r3 = std::exp(-1.15641e0_rt +
+    amrex::Real r3 = std::exp(-1.15641e0_rt +
                        ct9i +
                        +clnt9);
 
-    Real dr3dt = r3 * tfactors.t9i * (-ct9i +
+    amrex::Real dr3dt = r3 * tfactors.t9i * (-ct9i +
                                       (-1.5e0_rt));
 
     ct9i = tfactors.t9i * (-5.97632e0_rt);
     clnt9 = tfactors.lnt9 * (-1.5e0_rt);
 
-    Real r4 = std::exp(-6.65137e0_rt +
+    amrex::Real r4 = std::exp(-6.65137e0_rt +
                        ct9i +
                        clnt9);
 
-    Real dr4dt = r4 * tfactors.t9i * (-ct9i +
+    amrex::Real dr4dt = r4 * tfactors.t9i * (-ct9i +
                                       (-1.5e0_rt));
 
     ct9i = tfactors.t9i * (-3.14078e0_rt);
-    Real ct9i13 = tfactors.t9i13 * (-17.7759e0_rt);
-    Real ct913 = tfactors.t913 * (36.0724e0_rt);
-    Real ct9 = tfactors.t9 * (-5.34039e0_rt);
-    Real ct953 = tfactors.t953 * (0.382679e0_rt);
+    amrex::Real ct9i13 = tfactors.t9i13 * (-17.7759e0_rt);
+    amrex::Real ct913 = tfactors.t913 * (36.0724e0_rt);
+    amrex::Real ct9 = tfactors.t9 * (-5.34039e0_rt);
+    amrex::Real ct953 = tfactors.t953 * (0.382679e0_rt);
     clnt9 = tfactors.lnt9 * (-1.5e0_rt);
 
-    Real r5 = std::exp(-13.3392e0_rt +
+    amrex::Real r5 = std::exp(-13.3392e0_rt +
                        ct9i +
                        ct9i13 +
                        ct913 +
@@ -362,7 +362,7 @@ void rate_he4_ne18_to_p_na21 (const tf_t& tfactors, Real& rate, Real& dratedt)
                        ct953 +
                        clnt9);
 
-    Real dr5dt = r5 * tfactors.t9i * (-ct9i -
+    amrex::Real dr5dt = r5 * tfactors.t9i * (-ct9i -
                                       (1.0_rt / 3.0_rt) * ct9i13 +
                                       (1.0_rt / 3.0_rt) * ct913 +
                                       ct9 +
@@ -372,11 +372,11 @@ void rate_he4_ne18_to_p_na21 (const tf_t& tfactors, Real& rate, Real& dratedt)
     ct9i = tfactors.t9i * (-2.81989e0_rt);
     clnt9 = tfactors.lnt9 * (-1.5e0_rt);
 
-    Real r6 = std::exp(-28.6929e0_rt +
+    amrex::Real r6 = std::exp(-28.6929e0_rt +
                        ct9i +
                        clnt9);
 
-    Real dr6dt = r6 * tfactors.t9i * (-ct9i +
+    amrex::Real dr6dt = r6 * tfactors.t9i * (-ct9i +
                                       (-1.5e0_rt));
 
     rate = r0 + r1 + r2 + r3 + r4 + r5 + r6;
@@ -385,7 +385,7 @@ void rate_he4_ne18_to_p_na21 (const tf_t& tfactors, Real& rate, Real& dratedt)
 
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-void rate_ne18_to_f18 ([[maybe_unused]] const tf_t& tfactors, Real& rate, Real& dratedt)
+void rate_ne18_to_f18 ([[maybe_unused]] const tf_t& tfactors, amrex::Real& rate, amrex::Real& dratedt)
 {
     rate = std::exp(-0.880534e0_rt);
     dratedt = 0.0_rt;
@@ -393,7 +393,7 @@ void rate_ne18_to_f18 ([[maybe_unused]] const tf_t& tfactors, Real& rate, Real&
 
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-void rate_ne19_to_f19 ([[maybe_unused]] const tf_t& tfactors, Real& rate, Real& dratedt)
+void rate_ne19_to_f19 ([[maybe_unused]] const tf_t& tfactors, amrex::Real& rate, amrex::Real& dratedt)
 {
     rate = std::exp(-3.21258e0_rt);
     dratedt = 0.0_rt;
@@ -401,44 +401,44 @@ void rate_ne19_to_f19 ([[maybe_unused]] const tf_t& tfactors, Real& rate, Real&
 
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-void rate_p_ne19_to_na20 (const tf_t& tfactors, Real& rate, Real& dratedt)
+void rate_p_ne19_to_na20 (const tf_t& tfactors, amrex::Real& rate, amrex::Real& dratedt)
 {
     rate = 0.0_rt;
     dratedt = 0.0_rt;
 
-    Real ct9i = tfactors.t9i * (-5.07623e0_rt);
-    Real ct913 = tfactors.t913 * (1.23704e0_rt);
-    Real ct9 = tfactors.t9 * (0.337618e0_rt);
-    Real ct953 = tfactors.t953 * (-0.0562825e0_rt);
-    Real clnt9 = tfactors.lnt9 * (-1.5e0_rt);
+    amrex::Real ct9i = tfactors.t9i * (-5.07623e0_rt);
+    amrex::Real ct913 = tfactors.t913 * (1.23704e0_rt);
+    amrex::Real ct9 = tfactors.t9 * (0.337618e0_rt);
+    amrex::Real ct953 = tfactors.t953 * (-0.0562825e0_rt);
+    amrex::Real clnt9 = tfactors.lnt9 * (-1.5e0_rt);
 
-    Real r0 = std::exp(5.63289e0_rt +
+    amrex::Real r0 = std::exp(5.63289e0_rt +
                        ct9i +
                        ct913 +
                        ct9 +
                        ct953 +
                        clnt9);
 
-    Real dr0dt = r0 * tfactors.t9i * (-ct9i +
+    amrex::Real dr0dt = r0 * tfactors.t9i * (-ct9i +
                                       (1.0_rt / 3.0_rt) * ct913 +
                                       ct9 +
                                       (5.0_rt / 3.0_rt) * ct953 +
                                       (-1.5e0_rt));
 
-    Real ct9i13 = tfactors.t9i13 * (-19.5908e0_rt);
+    amrex::Real ct9i13 = tfactors.t9i13 * (-19.5908e0_rt);
     ct913 = tfactors.t913 * (-2.37696e0_rt);
     ct9 = tfactors.t9 * (3.26815e0_rt);
     ct953 = tfactors.t953 * (-1.06524e0_rt);
     clnt9 = tfactors.lnt9 * (-0.666667e0_rt);
 
-    Real r1 = std::exp(17.822e0_rt +
+    amrex::Real r1 = std::exp(17.822e0_rt +
                        ct9i13 +
                        ct913 +
                        ct9 +
                        ct953 +
                        clnt9);
 
-    Real dr1dt = r1 * tfactors.t9i * (-(1.0_rt / 3.0_rt) * ct9i13 +
+    amrex::Real dr1dt = r1 * tfactors.t9i * (-(1.0_rt / 3.0_rt) * ct9i13 +
                                       (1.0_rt / 3.0_rt) * ct913 +
                                       ct9 +
                                       (5.0_rt / 3.0_rt) * ct953 +
@@ -450,45 +450,45 @@ void rate_p_ne19_to_na20 (const tf_t& tfactors, Real& rate, Real& dratedt)
 
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-void rate_he4_o14_to_p_f17 (const tf_t& tfactors, Real& rate, Real& dratedt)
+void rate_he4_o14_to_p_f17 (const tf_t& tfactors, amrex::Real& rate, amrex::Real& dratedt)
 {
     rate = 0.0_rt;
     dratedt = 0.0_rt;
 
-    Real ct9i = tfactors.t9i * (-12.0223e0_rt);
-    Real clnt9 = tfactors.lnt9 * (-1.5e0_rt);
+    amrex::Real ct9i = tfactors.t9i * (-12.0223e0_rt);
+    amrex::Real clnt9 = tfactors.lnt9 * (-1.5e0_rt);
 
-    Real r0 = std::exp(12.1289e0_rt +
+    amrex::Real r0 = std::exp(12.1289e0_rt +
                        ct9i +
                        clnt9);
 
-    Real dr0dt = r0 * tfactors.t9i * (-ct9i +
+    amrex::Real dr0dt = r0 * tfactors.t9i * (-ct9i +
                                       (-1.5e0_rt));
 
     ct9i = tfactors.t9i * (-26.0e0_rt);
     clnt9 = tfactors.lnt9 * (-1.5e0_rt);
 
-    Real r1 = std::exp(18.6518e0_rt +
+    amrex::Real r1 = std::exp(18.6518e0_rt +
                        ct9i +
                        clnt9);
 
-    Real dr1dt = r1 * tfactors.t9i * (-ct9i +
+    amrex::Real dr1dt = r1 * tfactors.t9i * (-ct9i +
                                       (-1.5e0_rt));
 
-    Real ct9i13 = tfactors.t9i13 * (-39.388e0_rt);
-    Real ct913 = tfactors.t913 * (-17.4673e0_rt);
-    Real ct9 = tfactors.t9 * (35.3029e0_rt);
-    Real ct953 = tfactors.t953 * (-24.8162e0_rt);
+    amrex::Real ct9i13 = tfactors.t9i13 * (-39.388e0_rt);
+    amrex::Real ct913 = tfactors.t913 * (-17.4673e0_rt);
+    amrex::Real ct9 = tfactors.t9 * (35.3029e0_rt);
+    amrex::Real ct953 = tfactors.t953 * (-24.8162e0_rt);
     clnt9 = tfactors.lnt9 * (-0.666667e0_rt);
 
-    Real r2 = std::exp(40.8358e0_rt +
+    amrex::Real r2 = std::exp(40.8358e0_rt +
                        ct9i13 +
                        ct913 +
                        ct9 +
                        ct953 +
                        clnt9);
 
-    Real dr2dt = r2 * tfactors.t9i * (-(1.0_rt / 3.0_rt) * ct9i13 +
+    amrex::Real dr2dt = r2 * tfactors.t9i * (-(1.0_rt / 3.0_rt) * ct9i13 +
                                       (1.0_rt / 3.0_rt) * ct913 +
                                       ct9 +
                                       (5.0_rt / 3.0_rt) * ct953 +
@@ -497,31 +497,31 @@ void rate_he4_o14_to_p_f17 (const tf_t& tfactors, Real& rate, Real& dratedt)
     ct9i = tfactors.t9i * (-22.51e0_rt);
     clnt9 = tfactors.lnt9 * (-1.5e0_rt);
 
-    Real r3 = std::exp(16.3087e0_rt +
+    amrex::Real r3 = std::exp(16.3087e0_rt +
                        ct9i +
                        clnt9);
 
-    Real dr3dt = r3 * tfactors.t9i * (-ct9i +
+    amrex::Real dr3dt = r3 * tfactors.t9i * (-ct9i +
                                       (-1.5e0_rt));
 
     ct9i = tfactors.t9i * (-13.6e0_rt);
     clnt9 = tfactors.lnt9 * (-1.5e0_rt);
 
-    Real r4 = std::exp(11.1184e0_rt +
+    amrex::Real r4 = std::exp(11.1184e0_rt +
                        ct9i +
                        clnt9);
 
-    Real dr4dt = r4 * tfactors.t9i * (-ct9i +
+    amrex::Real dr4dt = r4 * tfactors.t9i * (-ct9i +
                                       (-1.5e0_rt));
 
     ct9i = tfactors.t9i * (-0.453036e0_rt);
     clnt9 = tfactors.lnt9 * (-1.5e0_rt);
 
-    Real r5 = std::exp(-106.091e0_rt +
+    amrex::Real r5 = std::exp(-106.091e0_rt +
                        ct9i +
                        clnt9);
 
-    Real dr5dt = r5 * tfactors.t9i * (-ct9i +
+    amrex::Real dr5dt = r5 * tfactors.t9i * (-ct9i +
                                       (-1.5e0_rt));
 
     rate = r0 + r1 + r2 + r3 + r4 + r5;
@@ -530,7 +530,7 @@ void rate_he4_o14_to_p_f17 (const tf_t& tfactors, Real& rate, Real& dratedt)
 
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-void rate_o14_to_n14 ([[maybe_unused]] const tf_t& tfactors, Real& rate, Real& dratedt)
+void rate_o14_to_n14 ([[maybe_unused]] const tf_t& tfactors, amrex::Real& rate, amrex::Real& dratedt)
 {
     rate = std::exp(-4.62412e0_rt);
     dratedt = 0.0_rt;
@@ -538,42 +538,42 @@ void rate_o14_to_n14 ([[maybe_unused]] const tf_t& tfactors, Real& rate, Real& d
 
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-void rate_he4_o15_to_ne19 (const tf_t& tfactors, Real& rate, Real& dratedt)
+void rate_he4_o15_to_ne19 (const tf_t& tfactors, amrex::Real& rate, amrex::Real& dratedt)
 {
     rate = 0.0_rt;
     dratedt = 0.0_rt;
 
-    Real ct9i = tfactors.t9i * (-5.88439e0_rt);
-    Real clnt9 = tfactors.lnt9 * (-1.5e0_rt);
+    amrex::Real ct9i = tfactors.t9i * (-5.88439e0_rt);
+    amrex::Real clnt9 = tfactors.lnt9 * (-1.5e0_rt);
 
-    Real r0 = std::exp(-0.0452465e0_rt +
+    amrex::Real r0 = std::exp(-0.0452465e0_rt +
                        ct9i +
                        clnt9);
 
-    Real dr0dt = r0 * tfactors.t9i * (-ct9i +
+    amrex::Real dr0dt = r0 * tfactors.t9i * (-ct9i +
                                       (-1.5e0_rt));
 
-    Real ct9i13 = tfactors.t9i13 * (-39.578e0_rt);
-    Real ct953 = tfactors.t953 * (-3.0e0_rt);
+    amrex::Real ct9i13 = tfactors.t9i13 * (-39.578e0_rt);
+    amrex::Real ct953 = tfactors.t953 * (-3.0e0_rt);
     clnt9 = tfactors.lnt9 * (-0.666667e0_rt);
 
-    Real r1 = std::exp(26.2914e0_rt +
+    amrex::Real r1 = std::exp(26.2914e0_rt +
                        ct9i13 +
                        ct953 +
                        clnt9);
 
-    Real dr1dt = r1 * tfactors.t9i * (-(1.0_rt / 3.0_rt) * ct9i13 +
+    amrex::Real dr1dt = r1 * tfactors.t9i * (-(1.0_rt / 3.0_rt) * ct9i13 +
                                       (5.0_rt / 3.0_rt) * ct953 +
                                       (-0.666667e0_rt));
 
     ct9i = tfactors.t9i * (-4.20439e0_rt);
     ct9i13 = tfactors.t9i13 * (-3.24609e0_rt);
-    Real ct913 = tfactors.t913 * (44.4647e0_rt);
-    Real ct9 = tfactors.t9 * (-9.79962e0_rt);
+    amrex::Real ct913 = tfactors.t913 * (44.4647e0_rt);
+    amrex::Real ct9 = tfactors.t9 * (-9.79962e0_rt);
     ct953 = tfactors.t953 * (0.841782e0_rt);
     clnt9 = tfactors.lnt9 * (-1.5e0_rt);
 
-    Real r2 = std::exp(-32.2496e0_rt +
+    amrex::Real r2 = std::exp(-32.2496e0_rt +
                        ct9i +
                        ct9i13 +
                        ct913 +
@@ -581,7 +581,7 @@ void rate_he4_o15_to_ne19 (const tf_t& tfactors, Real& rate, Real& dratedt)
                        ct953 +
                        clnt9);
 
-    Real dr2dt = r2 * tfactors.t9i * (-ct9i -
+    amrex::Real dr2dt = r2 * tfactors.t9i * (-ct9i -
                                       (1.0_rt / 3.0_rt) * ct9i13 +
                                       (1.0_rt / 3.0_rt) * ct913 +
                                       ct9 +
@@ -594,7 +594,7 @@ void rate_he4_o15_to_ne19 (const tf_t& tfactors, Real& rate, Real& dratedt)
 
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-void rate_o15_to_n15 ([[maybe_unused]] const tf_t& tfactors, Real& rate, Real& dratedt)
+void rate_o15_to_n15 ([[maybe_unused]] const tf_t& tfactors, amrex::Real& rate, amrex::Real& dratedt)
 {
     rate = std::exp(-5.1725e0_rt);
     dratedt = 0.0_rt;
@@ -602,35 +602,35 @@ void rate_o15_to_n15 ([[maybe_unused]] const tf_t& tfactors, Real& rate, Real& d
 
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-void rate_he4_o16_to_ne20 (const tf_t& tfactors, Real& rate, Real& dratedt)
+void rate_he4_o16_to_ne20 (const tf_t& tfactors, amrex::Real& rate, amrex::Real& dratedt)
 {
     rate = 0.0_rt;
     dratedt = 0.0_rt;
 
-    Real ct9i = tfactors.t9i * (-10.3585e0_rt);
-    Real clnt9 = tfactors.lnt9 * (-1.5e0_rt);
+    amrex::Real ct9i = tfactors.t9i * (-10.3585e0_rt);
+    amrex::Real clnt9 = tfactors.lnt9 * (-1.5e0_rt);
 
-    Real r0 = std::exp(3.88571e0_rt +
+    amrex::Real r0 = std::exp(3.88571e0_rt +
                        ct9i +
                        clnt9);
 
-    Real dr0dt = r0 * tfactors.t9i * (-ct9i +
+    amrex::Real dr0dt = r0 * tfactors.t9i * (-ct9i +
                                       (-1.5e0_rt));
 
-    Real ct9i13 = tfactors.t9i13 * (-39.7262e0_rt);
-    Real ct913 = tfactors.t913 * (-0.210799e0_rt);
-    Real ct9 = tfactors.t9 * (0.442879e0_rt);
-    Real ct953 = tfactors.t953 * (-0.0797753e0_rt);
+    amrex::Real ct9i13 = tfactors.t9i13 * (-39.7262e0_rt);
+    amrex::Real ct913 = tfactors.t913 * (-0.210799e0_rt);
+    amrex::Real ct9 = tfactors.t9 * (0.442879e0_rt);
+    amrex::Real ct953 = tfactors.t953 * (-0.0797753e0_rt);
     clnt9 = tfactors.lnt9 * (-0.666667e0_rt);
 
-    Real r1 = std::exp(23.903e0_rt +
+    amrex::Real r1 = std::exp(23.903e0_rt +
                        ct9i13 +
                        ct913 +
                        ct9 +
                        ct953 +
                        clnt9);
 
-    Real dr1dt = r1 * tfactors.t9i * (-(1.0_rt / 3.0_rt) * ct9i13 +
+    amrex::Real dr1dt = r1 * tfactors.t9i * (-(1.0_rt / 3.0_rt) * ct9i13 +
                                       (1.0_rt / 3.0_rt) * ct913 +
                                       ct9 +
                                       (5.0_rt / 3.0_rt) * ct953 +
@@ -642,14 +642,14 @@ void rate_he4_o16_to_ne20 (const tf_t& tfactors, Real& rate, Real& dratedt)
     ct953 = tfactors.t953 * (-0.00107508e0_rt);
     clnt9 = tfactors.lnt9 * (-1.5e0_rt);
 
-    Real r2 = std::exp(9.50848e0_rt +
+    amrex::Real r2 = std::exp(9.50848e0_rt +
                        ct9i +
                        ct913 +
                        ct9 +
                        ct953 +
                        clnt9);
 
-    Real dr2dt = r2 * tfactors.t9i * (-ct9i +
+    amrex::Real dr2dt = r2 * tfactors.t9i * (-ct9i +
                                       (1.0_rt / 3.0_rt) * ct913 +
                                       ct9 +
                                       (5.0_rt / 3.0_rt) * ct953 +
@@ -661,25 +661,25 @@ void rate_he4_o16_to_ne20 (const tf_t& tfactors, Real& rate, Real& dratedt)
 
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-void rate_p_o16_to_f17 (const tf_t& tfactors, Real& rate, Real& dratedt)
+void rate_p_o16_to_f17 (const tf_t& tfactors, amrex::Real& rate, amrex::Real& dratedt)
 {
     rate = 0.0_rt;
     dratedt = 0.0_rt;
 
-    Real ct9i13 = tfactors.t9i13 * (-16.696e0_rt);
-    Real ct913 = tfactors.t913 * (-1.16252e0_rt);
-    Real ct9 = tfactors.t9 * (0.267703e0_rt);
-    Real ct953 = tfactors.t953 * (-0.0338411e0_rt);
-    Real clnt9 = tfactors.lnt9 * (-0.666667e0_rt);
+    amrex::Real ct9i13 = tfactors.t9i13 * (-16.696e0_rt);
+    amrex::Real ct913 = tfactors.t913 * (-1.16252e0_rt);
+    amrex::Real ct9 = tfactors.t9 * (0.267703e0_rt);
+    amrex::Real ct953 = tfactors.t953 * (-0.0338411e0_rt);
+    amrex::Real clnt9 = tfactors.lnt9 * (-0.666667e0_rt);
 
-    Real r0 = std::exp(19.0904e0_rt +
+    amrex::Real r0 = std::exp(19.0904e0_rt +
                        ct9i13 +
                        ct913 +
                        ct9 +
                        ct953 +
                        clnt9);
 
-    Real dr0dt = r0 * tfactors.t9i * (-(1.0_rt / 3.0_rt) * ct9i13 +
+    amrex::Real dr0dt = r0 * tfactors.t9i * (-(1.0_rt / 3.0_rt) * ct9i13 +
                                       (1.0_rt / 3.0_rt) * ct913 +
                                       ct9 +
                                       (5.0_rt / 3.0_rt) * ct953 +
@@ -691,25 +691,25 @@ void rate_p_o16_to_f17 (const tf_t& tfactors, Real& rate, Real& dratedt)
 
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-void rate_he4_si26_to_p_p29 (const tf_t& tfactors, Real& rate, Real& dratedt)
+void rate_he4_si26_to_p_p29 (const tf_t& tfactors, amrex::Real& rate, amrex::Real& dratedt)
 {
     rate = 0.0_rt;
     dratedt = 0.0_rt;
 
-    Real ct9i13 = tfactors.t9i13 * (-59.3013e0_rt);
-    Real ct913 = tfactors.t913 * (0.480742e0_rt);
-    Real ct9 = tfactors.t9 * (-0.834505e0_rt);
-    Real ct953 = tfactors.t953 * (0.0621841e0_rt);
-    Real clnt9 = tfactors.lnt9 * (-0.666667e0_rt);
+    amrex::Real ct9i13 = tfactors.t9i13 * (-59.3013e0_rt);
+    amrex::Real ct913 = tfactors.t913 * (0.480742e0_rt);
+    amrex::Real ct9 = tfactors.t9 * (-0.834505e0_rt);
+    amrex::Real ct953 = tfactors.t953 * (0.0621841e0_rt);
+    amrex::Real clnt9 = tfactors.lnt9 * (-0.666667e0_rt);
 
-    Real r0 = std::exp(48.8732e0_rt +
+    amrex::Real r0 = std::exp(48.8732e0_rt +
                        ct9i13 +
                        ct913 +
                        ct9 +
                        ct953 +
                        clnt9);
 
-    Real dr0dt = r0 * tfactors.t9i * (-(1.0_rt / 3.0_rt) * ct9i13 +
+    amrex::Real dr0dt = r0 * tfactors.t9i * (-(1.0_rt / 3.0_rt) * ct9i13 +
                                       (1.0_rt / 3.0_rt) * ct913 +
                                       ct9 +
                                       (5.0_rt / 3.0_rt) * ct953 +
@@ -721,19 +721,19 @@ void rate_he4_si26_to_p_p29 (const tf_t& tfactors, Real& rate, Real& dratedt)
 
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-void rate_he4_ti44_to_p_v47 (const tf_t& tfactors, Real& rate, Real& dratedt)
+void rate_he4_ti44_to_p_v47 (const tf_t& tfactors, amrex::Real& rate, amrex::Real& dratedt)
 {
     rate = 0.0_rt;
     dratedt = 0.0_rt;
 
-    Real ct9i = tfactors.t9i * (-9.07869e0_rt);
-    Real ct9i13 = tfactors.t9i13 * (5.56533e0_rt);
-    Real ct913 = tfactors.t913 * (18.4415e0_rt);
-    Real ct9 = tfactors.t9 * (-4.10095e0_rt);
-    Real ct953 = tfactors.t953 * (0.24244e0_rt);
-    Real clnt9 = tfactors.lnt9 * (16.0516e0_rt);
+    amrex::Real ct9i = tfactors.t9i * (-9.07869e0_rt);
+    amrex::Real ct9i13 = tfactors.t9i13 * (5.56533e0_rt);
+    amrex::Real ct913 = tfactors.t913 * (18.4415e0_rt);
+    amrex::Real ct9 = tfactors.t9 * (-4.10095e0_rt);
+    amrex::Real ct953 = tfactors.t953 * (0.24244e0_rt);
+    amrex::Real clnt9 = tfactors.lnt9 * (16.0516e0_rt);
 
-    Real r0 = std::exp(-34.2468e0_rt +
+    amrex::Real r0 = std::exp(-34.2468e0_rt +
                        ct9i +
                        ct9i13 +
                        ct913 +
@@ -741,7 +741,7 @@ void rate_he4_ti44_to_p_v47 (const tf_t& tfactors, Real& rate, Real& dratedt)
                        ct953 +
                        clnt9);
 
-    Real dr0dt = r0 * tfactors.t9i * (-ct9i -
+    amrex::Real dr0dt = r0 * tfactors.t9i * (-ct9i -
                                       (1.0_rt / 3.0_rt) * ct9i13 +
                                       (1.0_rt / 3.0_rt) * ct913 +
                                       ct9 +
diff --git a/screening/screen.H b/screening/screen.H
index 7f79f6fb14..fae39862d0 100644
--- a/screening/screen.H
+++ b/screening/screen.H
@@ -20,9 +20,7 @@
 #include <screen_data.H>
 #include <extern_parameters.H>
 
-using namespace amrex;
-using namespace integrator_rp;
-using namespace screening_rp;
+using namespace amrex::literals;
 
 #if SCREEN_METHOD == SCREEN_METHOD_null
 const std::string screen_name = "null";
@@ -38,22 +36,22 @@ const std::string screen_name = "chabrier1998";
 
 struct plasma_state_t {
 
-    Real qlam0z;
-    Real qlam0zdt;
-    //Real qlam0zdd;
+    amrex::Real qlam0z;
+    amrex::Real qlam0zdt;
+    //amrex::Real qlam0zdd;
 
-    Real taufac;
-    Real taufacdt;
+    amrex::Real taufac;
+    amrex::Real taufacdt;
 
-    Real aa;
-    Real daadt;
-    //Real daadd;
+    amrex::Real aa;
+    amrex::Real daadt;
+    //amrex::Real daadd;
 
-    Real temp;
-    Real zbar;
-    Real z2bar;
-    Real n_e;
-    Real gamma_e_fac;
+    amrex::Real temp;
+    amrex::Real zbar;
+    amrex::Real z2bar;
+    amrex::Real n_e;
+    amrex::Real gamma_e_fac;
 };
 
 inline
@@ -89,18 +87,18 @@ screening_finalize() {
 template <int do_T_derivatives>
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 void
-fill_plasma_state(plasma_state_t& state, const Real temp,
-                  const Real dens, Array1D<Real, 1, NumSpec> const& y) {
+fill_plasma_state(plasma_state_t& state, const amrex::Real temp,
+                  const amrex::Real dens, amrex::Array1D<amrex::Real, 1, NumSpec> const& y) {
 
-    Real sum = 0.0_rt;
+    amrex::Real sum = 0.0_rt;
     for (int n = 1; n <= NumSpec; n++) {
         sum += y(n);
     }
-    Real abar = 1.0_rt / sum;
-    Real ytot = sum;
+    amrex::Real abar = 1.0_rt / sum;
+    amrex::Real ytot = sum;
 
     sum = 0.0_rt;
-    Real sum2 = 0.0_rt;
+    amrex::Real sum2 = 0.0_rt;
     for (int n = 1; n <= NumSpec; n++) {
         sum += zion[n-1]*y(n);
         sum2 += zion[n-1]*zion[n-1]*y(n);
@@ -109,25 +107,25 @@ fill_plasma_state(plasma_state_t& state, const Real temp,
     // Part of Eq.6 in Itoh:1979
     // 4.248719e3 = (27*pi^2*e^4*m_u/(2*k_B*hbar^2))^(1/3)
     // the extra (1/3) to make tau -> tau/3
-    const Real co2 = (1.0_rt/3.0_rt) * 4.248719e3_rt;
+    const amrex::Real co2 = (1.0_rt/3.0_rt) * 4.248719e3_rt;
 
-    Real zbar = sum * abar;
-    Real z2bar = sum2 * abar;
+    amrex::Real zbar = sum * abar;
+    amrex::Real z2bar = sum2 * abar;
 
     // ntot
-    Real rr = dens * ytot;
-    Real tempi = 1.0_rt / temp;
-    [[maybe_unused]] Real dtempi;
+    amrex::Real rr = dens * ytot;
+    amrex::Real tempi = 1.0_rt / temp;
+    [[maybe_unused]] amrex::Real dtempi;
     if constexpr (do_T_derivatives) {
         dtempi = -tempi * tempi;
     }
 
     // Part of Eq. 19 in Graboske:1973
     // pp = sqrt( \tilde{z}*(rho/u_I/T) )
-    Real pp = std::sqrt(rr*tempi*(z2bar + zbar));
-    [[maybe_unused]] Real dppdt;
+    amrex::Real pp = std::sqrt(rr*tempi*(z2bar + zbar));
+    [[maybe_unused]] amrex::Real dppdt;
     if constexpr (do_T_derivatives) {
-        Real qq = 0.5_rt/pp *(z2bar + zbar);
+        amrex::Real qq = 0.5_rt/pp *(z2bar + zbar);
         dppdt = qq*rr*dtempi;
     }
 
@@ -143,7 +141,7 @@ fill_plasma_state(plasma_state_t& state, const Real temp,
         state.taufacdt = -(1.0_rt/3.0_rt) * state.taufac * tempi;
     }
 
-    Real xni = std::cbrt(rr * zbar);
+    amrex::Real xni = std::cbrt(rr * zbar);
 
     // Part of Eq.4 in Itoh:1979
     // 2.27493e5 = e^2 / ( (3*m_u/(4pi))^(1/3) *k_B )
@@ -161,16 +159,16 @@ fill_plasma_state(plasma_state_t& state, const Real temp,
     state.n_e = zbar * rr * C::n_A;
 
     // precomputed part of Gamma_e, from Chugunov 2009 eq. 6
-    constexpr Real gamma_e_constants =
+    constexpr amrex::Real gamma_e_constants =
         C::q_e*C::q_e/C::k_B * gcem::pow(4.0_rt/3.0_rt*M_PI, 1.0_rt/3.0_rt);
     state.gamma_e_fac = gamma_e_constants * std::cbrt(state.n_e);
 }
 
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 void
-fill_plasma_state(plasma_state_t& state, const Real temp,
-                  const Real dens, Array1D<Real, 1, NumSpec> const& y) {
-    if (jacobian == 1) {
+fill_plasma_state(plasma_state_t& state, const amrex::Real temp,
+                  const amrex::Real dens, amrex::Array1D<amrex::Real, 1, NumSpec> const& y) {
+    if (integrator_rp::jacobian == 1) {
         constexpr int do_T_derivatives = 1;
         fill_plasma_state<do_T_derivatives>(state, temp, dens, y);
     } else {
@@ -184,7 +182,7 @@ template <int do_T_derivatives>
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
 void actual_screen5 (const plasma_state_t& state,
                      const scrn::screen_factors_t& scn_fac,
-                     Real& scor, Real& scordt)
+                     amrex::Real& scor, amrex::Real& scordt)
 {
     // this subroutine calculates screening factors and their derivatives
     // for nuclear reaction rates in the weak, intermediate and strong regimes.
@@ -203,19 +201,19 @@ void actual_screen5 (const plasma_state_t& state,
 
 
     // fact = 2^(1/3)
-    const Real fact    = 1.25992104989487e0_rt;
-    const Real gamefx  = 0.3e0_rt;          // lower gamma limit for intermediate screening
-    const Real gamefs  = 0.8e0_rt;          // upper gamma limit for intermediate screening
-    const Real h12_max = 300.e0_rt;
+    const amrex::Real fact    = 1.25992104989487e0_rt;
+    const amrex::Real gamefx  = 0.3e0_rt;          // lower gamma limit for intermediate screening
+    const amrex::Real gamefs  = 0.8e0_rt;          // upper gamma limit for intermediate screening
+    const amrex::Real h12_max = 300.e0_rt;
 
     // Get the ion data based on the input index
-    Real z1 = scn_fac.z1;
-    Real z2 = scn_fac.z2;
+    amrex::Real z1 = scn_fac.z1;
+    amrex::Real z2 = scn_fac.z2;
 
     // calculate individual screening factors
-    Real bb = z1 * z2;
-    Real gamp = state.aa;
-    [[maybe_unused]] Real gampdt;
+    amrex::Real bb = z1 * z2;
+    amrex::Real gamp = state.aa;
+    [[maybe_unused]] amrex::Real gampdt;
     if constexpr (do_T_derivatives) {
         gampdt = state.daadt;
     }
@@ -223,12 +221,12 @@ void actual_screen5 (const plasma_state_t& state,
     // In Eq.4 in Itoh:1979, this term is 2*Z_1*Z_2/(Z_1^(1/3) + Z_2^(1/3))
     // However here we follow Wallace:1982 Eq. A13, which is Z_1*Z_2*(2/(Z_1+Z_2))^(1/3)
 
-    Real qq = fact * bb * scn_fac.zs13inv;
+    amrex::Real qq = fact * bb * scn_fac.zs13inv;
 
     // Full Equation of Wallace:1982 Eq. A13
 
-    Real gamef = qq * gamp;
-    [[maybe_unused]] Real gamefdt;
+    amrex::Real gamef = qq * gamp;
+    [[maybe_unused]] amrex::Real gamefdt;
     if constexpr (do_T_derivatives) {
         gamefdt = qq * gampdt;
     }
@@ -237,8 +235,8 @@ void actual_screen5 (const plasma_state_t& state,
     // the extra 1/3 factor is there for convenience.
     // tau12 = Eq.6 / 3
 
-    Real tau12 = state.taufac * scn_fac.aznut;
-    [[maybe_unused]] Real tau12dt;
+    amrex::Real tau12 = state.taufac * scn_fac.aznut;
+    [[maybe_unused]] amrex::Real tau12dt;
     if constexpr (do_T_derivatives) {
         tau12dt = state.taufacdt * scn_fac.aznut;
     }
@@ -247,8 +245,8 @@ void actual_screen5 (const plasma_state_t& state,
 
     // alph12 = 3*gamma_ij/tau_ij
 
-    Real alph12 = gamef * qq;
-    [[maybe_unused]] Real alph12dt;
+    amrex::Real alph12 = gamef * qq;
+    [[maybe_unused]] amrex::Real alph12dt;
     if constexpr (do_T_derivatives) {
         alph12dt = (gamefdt - alph12*tau12dt) * qq;
     }
@@ -282,14 +280,14 @@ void actual_screen5 (const plasma_state_t& state,
     // Full version of Eq. 19 in Graboske:1973 by considering weak regime
     // and Wallace:1982 Eq. A14. Here the degeneracy factor is assumed to be 1.
 
-    Real h12w = bb * state.qlam0z;
-    [[maybe_unused]] Real dh12wdt;
+    amrex::Real h12w = bb * state.qlam0z;
+    [[maybe_unused]] amrex::Real dh12wdt;
     if constexpr (do_T_derivatives) {
         dh12wdt = bb * state.qlam0zdt;
     }
 
-    Real h12 = h12w;
-    [[maybe_unused]] Real dh12dt;
+    amrex::Real h12 = h12w;
+    [[maybe_unused]] amrex::Real dh12dt;
     if constexpr (do_T_derivatives) {
         dh12dt = dh12wdt;
     }
@@ -300,60 +298,60 @@ void actual_screen5 (const plasma_state_t& state,
 
       // gamma_ij^(1/4)
 
-        Real gamp14 = std::pow(gamp, 0.25_rt);
-        Real rr = 1.0_rt/gamp;
+        amrex::Real gamp14 = std::pow(gamp, 0.25_rt);
+        amrex::Real rr = 1.0_rt/gamp;
 
         // Here we follow Eq. A9 in Wallace:1982
         // See Eq. 25 Alastuey:1978, Eq. 16 and 17 in Jancovici:1977 for reference
-        Real cc = 0.896434e0_rt * gamp * scn_fac.zhat
+        amrex::Real cc = 0.896434e0_rt * gamp * scn_fac.zhat
             - 3.44740e0_rt * gamp14 * scn_fac.zhat2
             - 0.5551e0_rt * (std::log(gamp) + scn_fac.lzav)
             - 2.996e0_rt;
 
-        [[maybe_unused]] Real dccdt;
+        [[maybe_unused]] amrex::Real dccdt;
         if constexpr (do_T_derivatives) {
             qq = 0.25_rt * gamp14 * rr;
-            Real gamp14dt = qq * gampdt;
+            amrex::Real gamp14dt = qq * gampdt;
             dccdt = 0.896434e0_rt * gampdt * scn_fac.zhat
                 - 3.44740e0_rt * gamp14dt * scn_fac.zhat2
                 - 0.5551e0_rt *rr * gampdt;
         }
 
         // (3gamma_ij/tau_ij)^3
-        Real a3 = alph12 * alph12 * alph12;
-        Real da3 = 3.0e0_rt * alph12 * alph12;
+        amrex::Real a3 = alph12 * alph12 * alph12;
+        amrex::Real da3 = 3.0e0_rt * alph12 * alph12;
 
         // Part of Eq. 28 in Alastuey:1978
         qq = 0.014e0_rt + 0.0128e0_rt*alph12;
 
         // Part of Eq. 28 in Alastuey:1978
         rr = (5.0_rt/32.0_rt) - alph12*qq;
-        [[maybe_unused]] Real drrdt;
+        [[maybe_unused]] amrex::Real drrdt;
         if constexpr (do_T_derivatives) {
-            Real dqqdt  = 0.0128e0_rt*alph12dt;
+            amrex::Real dqqdt  = 0.0128e0_rt*alph12dt;
             drrdt  = -(alph12dt*qq + alph12*dqqdt);
         }
 
         // Part of Eq. 28 in Alastuey:1978
-        Real ss = tau12*rr;
+        amrex::Real ss = tau12*rr;
 
         // Part of Eq. 31 in Alastuey:1978
-        Real tt = -0.0098e0_rt + 0.0048e0_rt*alph12;
+        amrex::Real tt = -0.0098e0_rt + 0.0048e0_rt*alph12;
 
         // Part of Eq. 31 in Alastuey:1978
-        Real uu = 0.0055e0_rt + alph12*tt;
+        amrex::Real uu = 0.0055e0_rt + alph12*tt;
 
         // Part of Eq. 31 in Alastuey:1978
-        Real vv = gamef * alph12 * uu;
+        amrex::Real vv = gamef * alph12 * uu;
 
         // Exponent of Eq. 32 in Alastuey:1978, which uses Eq.28 and Eq.31
         // Strong screening factor
         h12 = cc - a3 * (ss + vv);
         if constexpr (do_T_derivatives) {
-            Real dssdt  = tau12dt*rr + tau12*drrdt;
-            Real dttdt  = 0.0048e0_rt*alph12dt;
-            Real duudt  = alph12dt*tt + alph12*dttdt;
-            Real dvvdt = gamefdt*alph12*uu + gamef*alph12dt*uu + gamef*alph12*duudt;
+            amrex::Real dssdt  = tau12dt*rr + tau12*drrdt;
+            amrex::Real dttdt  = 0.0048e0_rt*alph12dt;
+            amrex::Real duudt  = alph12dt*tt + alph12*dttdt;
+            amrex::Real dvvdt = gamefdt*alph12*uu + gamef*alph12dt*uu + gamef*alph12*duudt;
             rr = da3 * (ss + vv);
             dh12dt  = dccdt - rr*alph12dt - a3*(dssdt + dvvdt);
         }
@@ -366,8 +364,8 @@ void actual_screen5 (const plasma_state_t& state,
             drrdt = ss*alph12dt;
         }
 
-        Real xlgfac;
-        [[maybe_unused]] Real dxlgfacdt;
+        amrex::Real xlgfac;
+        [[maybe_unused]] amrex::Real dxlgfacdt;
 
         // In extreme case, rr is 0.77, see conclusion in Alastuey:1978
         if (rr >= 0.77e0_rt) {
@@ -392,7 +390,7 @@ void actual_screen5 (const plasma_state_t& state,
         // If gamma_ij < upper limit of intermediate regime
         // then it is in the intermediate regime, else strong screening.
         if (gamef <= gamefs) {
-            Real dgamma = 1.0e0_rt/(gamefs - gamefx);
+            amrex::Real dgamma = 1.0e0_rt/(gamefs - gamefx);
 
             rr =  dgamma*(gamefs - gamef);
 
@@ -405,7 +403,7 @@ void actual_screen5 (const plasma_state_t& state,
             h12 = h12w*rr + vv*ss;
             if constexpr (do_T_derivatives) {
                 drrdt  = -dgamma*gamefdt;
-                Real dssdt = dgamma*gamefdt;
+                amrex::Real dssdt = dgamma*gamefdt;
                 dh12dt = dh12wdt*rr + h12w*drrdt + dh12dt*ss + vv*dssdt;
             }
         }
@@ -432,7 +430,7 @@ template <int do_T_derivatives>
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
 void chugunov2007 (const plasma_state_t& state,
                    const scrn::screen_factors_t& scn_fac,
-                   Real& scor, Real& scordt)
+                   amrex::Real& scor, amrex::Real& scordt)
 {
     // Calculates screening factors based on Chugunov et al. 2007, following the
     // the approach in Yakovlev 2006 to extend to a multi-component plasma.
@@ -450,7 +448,7 @@ void chugunov2007 (const plasma_state_t& state,
     // scor    = screening correction
     // scordt  = derivative of screening correction with temperature
 
-    Real tmp;
+    amrex::Real tmp;
     // Plasma temperature T_p
     // This formula comes from working backwards from zeta_ij (Chugunov 2009 eq. 12)
     // through Chugunov 2007 eq. 3 to Chugunov 2007 eq. 2.
@@ -470,25 +468,25 @@ void chugunov2007 (const plasma_state_t& state,
     //   Z^2 -> zbar^2
     //   n_i -> ntot
     //   m_i -> m_u * abar
-    Real mu12 = scn_fac.a1 * scn_fac.a2 / (scn_fac.a1 + scn_fac.a2);
-    Real z_factor = scn_fac.z1 * scn_fac.z2;
-    Real n_i = state.n_e / scn_fac.ztilde3;
-    Real m_i = 2.0_rt * mu12 / C::n_A;
+    amrex::Real mu12 = scn_fac.a1 * scn_fac.a2 / (scn_fac.a1 + scn_fac.a2);
+    amrex::Real z_factor = scn_fac.z1 * scn_fac.z2;
+    amrex::Real n_i = state.n_e / scn_fac.ztilde3;
+    amrex::Real m_i = 2.0_rt * mu12 / C::n_A;
 
-    constexpr Real T_p_factor = C::hbar/C::k_B*C::q_e*gcem::sqrt(4.0_rt*GCEM_PI);
-    Real T_p = T_p_factor * std::sqrt(z_factor * n_i / m_i);
+    constexpr amrex::Real T_p_factor = C::hbar/C::k_B*C::q_e*gcem::sqrt(4.0_rt*GCEM_PI);
+    amrex::Real T_p = T_p_factor * std::sqrt(z_factor * n_i / m_i);
 
     // Normalized temperature
-    Real inv_T_p = 1.0_rt / T_p;
-    Real T_norm = state.temp * inv_T_p;
-    Real dT_norm_dT = inv_T_p;
+    amrex::Real inv_T_p = 1.0_rt / T_p;
+    amrex::Real T_norm = state.temp * inv_T_p;
+    amrex::Real dT_norm_dT = inv_T_p;
 
     // The fit has only been verified down to T ~ 0.1 T_p, below which the rate
     // should be nearly temperature-independent (in the pycnonuclear regime),
     // and we clip the temperature to 0.1 T_p at small T.
     // start the transition here
-    constexpr Real T_norm_fade = 0.2_rt;
-    constexpr Real T_norm_min = 0.1_rt;
+    constexpr amrex::Real T_norm_fade = 0.2_rt;
+    constexpr amrex::Real T_norm_min = 0.1_rt;
 
     if (T_norm < T_norm_min) {
         // clip temperature to the minimum value
@@ -496,26 +494,26 @@ void chugunov2007 (const plasma_state_t& state,
         dT_norm_dT = 0.0_rt;
     } else if (T_norm <= T_norm_fade) {
         // blend using a cosine, after MESA
-        constexpr Real delta_T = T_norm_fade - T_norm_min;
+        constexpr amrex::Real delta_T = T_norm_fade - T_norm_min;
         tmp = M_PI * (T_norm - T_norm_min) / delta_T;
-        Real f = 0.5_rt * (1.0_rt - std::cos(tmp));
+        amrex::Real f = 0.5_rt * (1.0_rt - std::cos(tmp));
         if constexpr (do_T_derivatives) {
-            Real df_dT = 0.5_rt * M_PI / delta_T * std::sin(tmp) * dT_norm_dT;
+            amrex::Real df_dT = 0.5_rt * M_PI / delta_T * std::sin(tmp) * dT_norm_dT;
             dT_norm_dT = -df_dT * T_norm_min + df_dT * T_norm + f * dT_norm_dT;
         }
         T_norm = (1.0_rt - f) * T_norm_min + f * T_norm;
     }
 
     // Coulomb coupling parameter from Yakovlev 2006 eq. 10
-    Real Gamma = state.gamma_e_fac*scn_fac.z1*scn_fac.z2 / (scn_fac.ztilde*T_norm*T_p);
-    [[maybe_unused]] Real dGamma_dT;
+    amrex::Real Gamma = state.gamma_e_fac*scn_fac.z1*scn_fac.z2 / (scn_fac.ztilde*T_norm*T_p);
+    [[maybe_unused]] amrex::Real dGamma_dT;
     if constexpr (do_T_derivatives) {
         dGamma_dT = -Gamma / T_norm * dT_norm_dT;
     }
 
     // The fit for Gamma is only applicable up to ~600, so smoothly cap its value
-    constexpr Real Gamma_fade = 590;
-    constexpr Real Gamma_max = 600;
+    constexpr amrex::Real Gamma_fade = 590;
+    constexpr amrex::Real Gamma_max = 600;
 
     if (Gamma > Gamma_max) {
         // clip Gamma to the max value
@@ -523,57 +521,57 @@ void chugunov2007 (const plasma_state_t& state,
         dGamma_dT = 0.0_rt;
     } else if (Gamma >= Gamma_fade) {
         // blend using a cosine, after MESA
-        constexpr Real delta_gamma = Gamma_max - Gamma_fade;
+        constexpr amrex::Real delta_gamma = Gamma_max - Gamma_fade;
         tmp = M_PI * (Gamma - Gamma_fade) / delta_gamma;
-        Real f = 0.5_rt * (1.0_rt - std::cos(tmp));
+        amrex::Real f = 0.5_rt * (1.0_rt - std::cos(tmp));
         if constexpr (do_T_derivatives) {
-            Real df_dT = 0.5_rt * M_PI / delta_gamma * std::sin(tmp) * dGamma_dT;
+            amrex::Real df_dT = 0.5_rt * M_PI / delta_gamma * std::sin(tmp) * dGamma_dT;
             dGamma_dT = dGamma_dT - (df_dT * Gamma + f * dGamma_dT) + df_dT * Gamma_max;
         }
         Gamma = (1.0_rt - f) * Gamma + f * Gamma_max;
     }
 
     // Chugunov 2007 eq. 3
-    Real zeta = std::cbrt(4.0_rt / (3.0_rt * GCEM_PI*GCEM_PI * T_norm*T_norm));
-    [[maybe_unused]] Real dzeta_dT;
+    amrex::Real zeta = std::cbrt(4.0_rt / (3.0_rt * GCEM_PI*GCEM_PI * T_norm*T_norm));
+    [[maybe_unused]] amrex::Real dzeta_dT;
     if constexpr (do_T_derivatives) {
         dzeta_dT = -2.0_rt / (3.0_rt * T_norm) * zeta * dT_norm_dT;
     }
 
     // Gamma tilde from Chugunov 2007 eq. 21
-    constexpr Real fit_alpha = 0.022_rt;
-    Real fit_beta = 0.41_rt - 0.6_rt / Gamma;
-    Real fit_gamma = 0.06_rt + 2.2_rt / Gamma;
+    constexpr amrex::Real fit_alpha = 0.022_rt;
+    amrex::Real fit_beta = 0.41_rt - 0.6_rt / Gamma;
+    amrex::Real fit_gamma = 0.06_rt + 2.2_rt / Gamma;
     // Polynomial term in Gamma tilde
-    Real poly = 1.0_rt + zeta*(fit_alpha + zeta*(fit_beta + fit_gamma*zeta));
-    [[maybe_unused]] Real dpoly_dT;
+    amrex::Real poly = 1.0_rt + zeta*(fit_alpha + zeta*(fit_beta + fit_gamma*zeta));
+    [[maybe_unused]] amrex::Real dpoly_dT;
     if constexpr (do_T_derivatives) {
         tmp = dGamma_dT / (Gamma * Gamma);
-        Real dfit_beta_dT = 0.6_rt * tmp;
-        Real dfit_gamma_dT = -2.2_rt * tmp;
+        amrex::Real dfit_beta_dT = 0.6_rt * tmp;
+        amrex::Real dfit_gamma_dT = -2.2_rt * tmp;
         dpoly_dT = (fit_alpha + 2.0_rt*zeta*fit_beta + 3.0_rt*fit_gamma*zeta*zeta) * dzeta_dT
                    + zeta*zeta*(dfit_beta_dT + dfit_gamma_dT*zeta);
     }
 
-    Real gamtilde = Gamma / std::cbrt(poly);
+    amrex::Real gamtilde = Gamma / std::cbrt(poly);
     // this is gamtilde * dlog(gamtilde)/dT
-    [[maybe_unused]] Real dgamtilde_dT;
+    [[maybe_unused]] amrex::Real dgamtilde_dT;
     if constexpr (do_T_derivatives) {
         dgamtilde_dT = gamtilde * (dGamma_dT / Gamma - dpoly_dT / poly / 3.0_rt);
     }
 
     // fit parameters just after Chugunov 2007 eq. 19
-    constexpr Real A1 = 2.7822_rt;
-    constexpr Real A2 = 98.34_rt;
-    constexpr Real A3 = gcem::sqrt(3.0_rt) - A1 / gcem::sqrt(A2);
-    const Real B1 = -1.7476_rt;
-    const Real B2 = 66.07_rt;
-    const Real B3 = 1.12_rt;
-    const Real B4 = 65_rt;
-    Real gamtilde2 = gamtilde * gamtilde;
-
-    Real term1, term2, term3, term4;
-    [[maybe_unused]] Real dterm1_dT, dterm2_dT, dterm3_dT, dterm4_dT;
+    constexpr amrex::Real A1 = 2.7822_rt;
+    constexpr amrex::Real A2 = 98.34_rt;
+    constexpr amrex::Real A3 = gcem::sqrt(3.0_rt) - A1 / gcem::sqrt(A2);
+    const amrex::Real B1 = -1.7476_rt;
+    const amrex::Real B2 = 66.07_rt;
+    const amrex::Real B3 = 1.12_rt;
+    const amrex::Real B4 = 65_rt;
+    amrex::Real gamtilde2 = gamtilde * gamtilde;
+
+    amrex::Real term1, term2, term3, term4;
+    [[maybe_unused]] amrex::Real dterm1_dT, dterm2_dT, dterm3_dT, dterm4_dT;
     // Chugunov 2007 eq. 19
     term1 = 1.0_rt / std::sqrt(A2 + gamtilde);
     if constexpr (do_T_derivatives) {
@@ -597,20 +595,20 @@ void chugunov2007 (const plasma_state_t& state,
         dterm4_dT = B4 / (tmp * tmp) * 2.0_rt * gamtilde * dgamtilde_dT;
     }
 
-    Real inner = A1 * term1 + A3 * term2;
+    amrex::Real inner = A1 * term1 + A3 * term2;
 
-    Real gamtilde32 = std::pow(gamtilde, 1.5_rt);
-    Real h = gamtilde32 * inner + B1 * term3 + B3 * term4;
-    [[maybe_unused]] Real dh_dT;
+    amrex::Real gamtilde32 = std::pow(gamtilde, 1.5_rt);
+    amrex::Real h = gamtilde32 * inner + B1 * term3 + B3 * term4;
+    [[maybe_unused]] amrex::Real dh_dT;
     if constexpr (do_T_derivatives) {
-        Real dinner_dT = A1 * dterm1_dT + A3 * dterm2_dT;
-        Real dgamtilde32_dT = 1.5_rt * std::sqrt(gamtilde) * dgamtilde_dT;
+        amrex::Real dinner_dT = A1 * dterm1_dT + A3 * dterm2_dT;
+        amrex::Real dgamtilde32_dT = 1.5_rt * std::sqrt(gamtilde) * dgamtilde_dT;
         dh_dT = dgamtilde32_dT * inner + gamtilde32 * dinner_dT
                 + B1 * dterm3_dT + B3 * dterm4_dT;
     }
 
     // machine limit the output
-    constexpr Real h_max = 300.e0_rt;
+    constexpr amrex::Real h_max = 300.e0_rt;
     h = amrex::min(h, h_max);
     scor = std::exp(h);
 
@@ -626,7 +624,7 @@ void chugunov2007 (const plasma_state_t& state,
 #elif SCREEN_METHOD == SCREEN_METHOD_chugunov2009
 template <int do_T_derivatives>
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-void chugunov2009_f0 (const Real gamma, const Real dlog_dT, Real& f, Real& df_dT)
+void chugunov2009_f0 (const amrex::Real gamma, const amrex::Real dlog_dT, amrex::Real& f, amrex::Real& df_dT)
 {
     // Calculates the free energy per ion in a OCP, from Chugunov and DeWitt 2009
     // equation 24.
@@ -640,17 +638,17 @@ void chugunov2009_f0 (const Real gamma, const Real dlog_dT, Real& f, Real& df_dT
     // df_dT   = derivative of free energy with temperature
 
     // fit parameters
-    constexpr Real A1 = -0.907_rt;
-    constexpr Real A2 = 0.62954_rt;
-    constexpr Real A3 = -gcem::sqrt(3.0_rt) / 2.0_rt - A1 / gcem::sqrt(A2);
-    constexpr Real B1 = 0.00456_rt;
-    constexpr Real B2 = 211.6_rt;
-    constexpr Real B3 = -1e-4_rt;
-    constexpr Real B4 = 0.00462_rt;
-    Real gamma_12 = std::sqrt(gamma);
-
-    Real term1, term2, term3, term4, term5;
-    Real dterm1_dgamma, dterm2_dgamma, dterm3_dgamma, dterm4_dgamma, dterm5_dgamma;
+    constexpr amrex::Real A1 = -0.907_rt;
+    constexpr amrex::Real A2 = 0.62954_rt;
+    constexpr amrex::Real A3 = -gcem::sqrt(3.0_rt) / 2.0_rt - A1 / gcem::sqrt(A2);
+    constexpr amrex::Real B1 = 0.00456_rt;
+    constexpr amrex::Real B2 = 211.6_rt;
+    constexpr amrex::Real B3 = -1e-4_rt;
+    constexpr amrex::Real B4 = 0.00462_rt;
+    amrex::Real gamma_12 = std::sqrt(gamma);
+
+    amrex::Real term1, term2, term3, term4, term5;
+    amrex::Real dterm1_dgamma, dterm2_dgamma, dterm3_dgamma, dterm4_dgamma, dterm5_dgamma;
 
     term1 = gamma_12 * std::sqrt(A2 + gamma);
     if constexpr (do_T_derivatives) {
@@ -690,7 +688,7 @@ template <int do_T_derivatives>
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
 void chugunov2009 (const plasma_state_t& state,
                    const scrn::screen_factors_t& scn_fac,
-                   Real& scor, Real& scordt)
+                   amrex::Real& scor, amrex::Real& scordt)
 {
     // Calculates screening factors based on Chugunov and DeWitt 2009, PhRvC, 80, 014611
 
@@ -702,88 +700,88 @@ void chugunov2009 (const plasma_state_t& state,
     // scor    = screening correction
     // scordt  = derivative of screening correction with temperature
 
-    Real z1z2 = scn_fac.z1 * scn_fac.z2;
-    Real zcomp = scn_fac.z1 + scn_fac.z2;
+    amrex::Real z1z2 = scn_fac.z1 * scn_fac.z2;
+    amrex::Real zcomp = scn_fac.z1 + scn_fac.z2;
 
     // Gamma_e from eq. 6
-    Real Gamma_e = state.gamma_e_fac / state.temp;
+    amrex::Real Gamma_e = state.gamma_e_fac / state.temp;
     // work in terms of log derivatives, since it's much simpler
-    Real dlog_Gamma_dT = -1.0_rt / state.temp;
+    amrex::Real dlog_Gamma_dT = -1.0_rt / state.temp;
 
     // Coulomb coupling parameters for ions and compound nucleus, eqs. 7 & 9
-    Real Gamma_1 = Gamma_e * scn_fac.z1_53;
-    Real Gamma_2 = Gamma_e * scn_fac.z2_53;
-    Real Gamma_comp = Gamma_e * scn_fac.zs53;
+    amrex::Real Gamma_1 = Gamma_e * scn_fac.z1_53;
+    amrex::Real Gamma_2 = Gamma_e * scn_fac.z2_53;
+    amrex::Real Gamma_comp = Gamma_e * scn_fac.zs53;
 
-    Real Gamma_12 = Gamma_e * z1z2 / scn_fac.ztilde;
+    amrex::Real Gamma_12 = Gamma_e * z1z2 / scn_fac.ztilde;
 
     // Coulomb barrier penetrability, eq. 10
 
-    constexpr Real tau_factor = gcem::pow(
+    constexpr amrex::Real tau_factor = gcem::pow(
         27.0_rt/2.0_rt * amrex::Math::powi<2>(M_PI*C::q_e*C::q_e/C::hbar)
         / (C::n_A*C::k_B), 1.0_rt/3.0_rt);
 
-    Real tau_12 = tau_factor * scn_fac.aznut / std::cbrt(state.temp);
-    Real dlog_tau_12_dT;
+    amrex::Real tau_12 = tau_factor * scn_fac.aznut / std::cbrt(state.temp);
+    amrex::Real dlog_tau_12_dT;
     if constexpr (do_T_derivatives) {
         dlog_tau_12_dT = -1.0_rt / state.temp / 3.0_rt;
     }
 
     // eq. 12
-    Real zeta = 3.0_rt * Gamma_12 / tau_12;
-    Real dzeta_dT;
+    amrex::Real zeta = 3.0_rt * Gamma_12 / tau_12;
+    amrex::Real dzeta_dT;
     if constexpr (do_T_derivatives) {
         dzeta_dT = zeta * (dlog_Gamma_dT - dlog_tau_12_dT);
     }
 
     // additional fit parameters, eq. 25
-    Real y_12 = 4.0_rt * z1z2 / (zcomp * zcomp);
-    Real c1 = 0.013_rt * y_12 * y_12;
-    Real c2 = 0.406_rt * std::pow(y_12, 0.14_rt);
-    Real c3 = 0.062_rt * std::pow(y_12, 0.19_rt) + 1.8_rt / Gamma_12;
+    amrex::Real y_12 = 4.0_rt * z1z2 / (zcomp * zcomp);
+    amrex::Real c1 = 0.013_rt * y_12 * y_12;
+    amrex::Real c2 = 0.406_rt * std::pow(y_12, 0.14_rt);
+    amrex::Real c3 = 0.062_rt * std::pow(y_12, 0.19_rt) + 1.8_rt / Gamma_12;
 
-    Real poly = 1.0_rt + zeta*(c1 + zeta*(c2 + c3*zeta));
-    Real t_12 = std::cbrt(poly);
-    Real dlog_dT = 0.0_rt;
+    amrex::Real poly = 1.0_rt + zeta*(c1 + zeta*(c2 + c3*zeta));
+    amrex::Real t_12 = std::cbrt(poly);
+    amrex::Real dlog_dT = 0.0_rt;
     if constexpr (do_T_derivatives) {
-        Real dc3_dT = -1.8_rt / Gamma_12 * dlog_Gamma_dT;
-        Real dpoly_dT = (c1 + zeta*(2.0_rt*c2 + 3.0_rt*c3*zeta))*dzeta_dT +
+        amrex::Real dc3_dT = -1.8_rt / Gamma_12 * dlog_Gamma_dT;
+        amrex::Real dpoly_dT = (c1 + zeta*(2.0_rt*c2 + 3.0_rt*c3*zeta))*dzeta_dT +
                         dc3_dT*zeta*zeta*zeta;
-        Real dlog_t_12_dT = dpoly_dT / (3.0_rt * poly);
+        amrex::Real dlog_t_12_dT = dpoly_dT / (3.0_rt * poly);
         dlog_dT = dlog_Gamma_dT - dlog_t_12_dT;
     }
 
     // strong screening enhancement factor, eq. 23, replacing tau_ij with t_ij
     // Using Gamma/tau_ij gives extremely low values, while Gamma/t_ij gives
     // values similar to those from Chugunov 2007.
-    Real term1, term2, term3;
-    Real dterm1_dT = 0.0_rt, dterm2_dT = 0.0_rt, dterm3_dT = 0.0_rt;
+    amrex::Real term1, term2, term3;
+    amrex::Real dterm1_dT = 0.0_rt, dterm2_dT = 0.0_rt, dterm3_dT = 0.0_rt;
     chugunov2009_f0<do_T_derivatives>(Gamma_1 / t_12, dlog_dT, term1, dterm1_dT);
     chugunov2009_f0<do_T_derivatives>(Gamma_2 / t_12, dlog_dT, term2, dterm2_dT);
     chugunov2009_f0<do_T_derivatives>(Gamma_comp / t_12, dlog_dT, term3, dterm3_dT);
-    Real h_fit = term1 + term2 - term3;
-    Real dh_fit_dT;
+    amrex::Real h_fit = term1 + term2 - term3;
+    amrex::Real dh_fit_dT;
     if constexpr (do_T_derivatives) {
         dh_fit_dT = dterm1_dT + dterm2_dT - dterm3_dT;
     }
 
     // weak screening correction term, eq. A3
-    Real corr_C = 3.0_rt*z1z2 * std::sqrt(state.z2bar/state.zbar) /
+    amrex::Real corr_C = 3.0_rt*z1z2 * std::sqrt(state.z2bar/state.zbar) /
                   (scn_fac.zs52 - scn_fac.z1_52 - scn_fac.z2_52);
 
     // corrected enhancement factor, eq. A4
-    Real Gamma_12_2 = Gamma_12 * Gamma_12;
-    Real numer = corr_C + Gamma_12_2;
-    Real denom = 1.0_rt + Gamma_12_2;
-    Real h12 = numer / denom * h_fit;
-    Real dh12_dT;
+    amrex::Real Gamma_12_2 = Gamma_12 * Gamma_12;
+    amrex::Real numer = corr_C + Gamma_12_2;
+    amrex::Real denom = 1.0_rt + Gamma_12_2;
+    amrex::Real h12 = numer / denom * h_fit;
+    amrex::Real dh12_dT;
     if constexpr (do_T_derivatives) {
-        Real dGamma_12_2_dT = 2 * Gamma_12_2 * dlog_Gamma_dT;
+        amrex::Real dGamma_12_2_dT = 2 * Gamma_12_2 * dlog_Gamma_dT;
         dh12_dT = h12 * (dGamma_12_2_dT/numer - dGamma_12_2_dT/denom + dh_fit_dT/h_fit);
     }
 
     // machine limit the output
-    const Real h12_max = 300.e0_rt;
+    const amrex::Real h12_max = 300.e0_rt;
     h12 = amrex::min(h12, h12_max);
     scor = std::exp(h12);
 
@@ -799,20 +797,20 @@ void chugunov2009 (const plasma_state_t& state,
 #elif SCREEN_METHOD == SCREEN_METHOD_chabrier1998
 template <int do_T_derivatives>
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-void chabrier1998_helmholtz_F(const Real gamma, const Real dgamma_dT,
-                              Real& f, Real& df_dT) {
+void chabrier1998_helmholtz_F(const amrex::Real gamma, const amrex::Real dgamma_dT,
+                              amrex::Real& f, amrex::Real& df_dT) {
     // Helmholtz free energy, See Chabrier & Potekhin 1998 Eq. 28
 
     // Fitted parameters, see Chabrier & Potekhin 1998 Sec.IV
-    constexpr Real A_1 = -0.9052_rt;
-    constexpr Real A_2 = 0.6322_rt;
-    constexpr Real A_3 = -0.5_rt * gcem::sqrt(3.0_rt) - A_1 / gcem::sqrt(A_2);
+    constexpr amrex::Real A_1 = -0.9052_rt;
+    constexpr amrex::Real A_2 = 0.6322_rt;
+    constexpr amrex::Real A_3 = -0.5_rt * gcem::sqrt(3.0_rt) - A_1 / gcem::sqrt(A_2);
 
     // Precompute some expressions that are reused in the derivative
-    const Real sqrt_gamma = std::sqrt(gamma);
-    const Real sqrt_1_gamma_A2 = std::sqrt(1.0_rt + gamma/A_2);
-    const Real sqrt_gamma_A2_gamma = std::sqrt(gamma * (A_2 + gamma));
-    const Real sqrt_gamma_A2 = std::sqrt(gamma/A_2);
+    const amrex::Real sqrt_gamma = std::sqrt(gamma);
+    const amrex::Real sqrt_1_gamma_A2 = std::sqrt(1.0_rt + gamma/A_2);
+    const amrex::Real sqrt_gamma_A2_gamma = std::sqrt(gamma * (A_2 + gamma));
+    const amrex::Real sqrt_gamma_A2 = std::sqrt(gamma/A_2);
 
     f = A_1 * (sqrt_gamma_A2_gamma -
                A_2 * std::log(sqrt_gamma_A2 + sqrt_1_gamma_A2)) +
@@ -832,7 +830,7 @@ template <int do_T_derivatives>
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
 void chabrier1998 (const plasma_state_t& state,
                    const scrn::screen_factors_t& scn_fac,
-                   Real& scor, Real& scordt)
+                   amrex::Real& scor, amrex::Real& scordt)
 {
     // Calculates screening factors based on Chabrier & Potekhin 1998,
     // Calder2007 and partly screen5 routine mentioned in Alastuey 1978.
@@ -847,15 +845,15 @@ void chabrier1998 (const plasma_state_t& state,
 
     // Eq. 2 in Chabrier & Potekhin 1998
 
-    Real Gamma_e = state.gamma_e_fac / state.temp;
+    amrex::Real Gamma_e = state.gamma_e_fac / state.temp;
 
     // See Calder2007 appendix Eq. A9
 
-    Real Gamma1 = Gamma_e * scn_fac.z1_53;
-    Real Gamma2 = Gamma_e * scn_fac.z2_53;
-    Real Gamma12 = Gamma_e * scn_fac.zs53;
+    amrex::Real Gamma1 = Gamma_e * scn_fac.z1_53;
+    amrex::Real Gamma2 = Gamma_e * scn_fac.z2_53;
+    amrex::Real Gamma12 = Gamma_e * scn_fac.zs53;
 
-    Real Gamma1dT{}, Gamma2dT{}, Gamma12dT{};
+    amrex::Real Gamma1dT{}, Gamma2dT{}, Gamma12dT{};
 
     if constexpr (do_T_derivatives) {
         Gamma1dT = -Gamma1 / state.temp;
@@ -865,8 +863,8 @@ void chabrier1998 (const plasma_state_t& state,
 
     // Helmholtz free energy
 
-    Real f1, f2, f12;
-    Real f1dT, f2dT, f12dT;
+    amrex::Real f1, f2, f12;
+    amrex::Real f1dT, f2dT, f12dT;
 
     chabrier1998_helmholtz_F<do_T_derivatives>(Gamma1, Gamma1dT, f1, f1dT);
     chabrier1998_helmholtz_F<do_T_derivatives>(Gamma2, Gamma2dT, f2, f2dT);
@@ -875,19 +873,19 @@ void chabrier1998 (const plasma_state_t& state,
     // Now we add quantum correction terms discussed in Alastuey 1978.
     // Notice in Alastuey 1978, they have a different classical term,
     // which is implemented in the strong screening limit of our screen5 routine.
-    Real quantum_corr_1 = 0.0_rt;
+    amrex::Real quantum_corr_1 = 0.0_rt;
 
-    Real quantum_corr_2 = 0.0_rt;
+    amrex::Real quantum_corr_2 = 0.0_rt;
 
-    [[maybe_unused]] Real quantum_corr_1_dT = 0.0_rt;
-    [[maybe_unused]] Real quantum_corr_2_dT = 0.0_rt;
+    [[maybe_unused]] amrex::Real quantum_corr_1_dT = 0.0_rt;
+    [[maybe_unused]] amrex::Real quantum_corr_2_dT = 0.0_rt;
 
-    if (enable_chabrier1998_quantum_corr) {
+    if (screening_rp::enable_chabrier1998_quantum_corr) {
         // See Wallace1982, Eq. A13
 
-        Real Gamma_eff = std::cbrt(2.0_rt) * scn_fac.z1 * scn_fac.z2 *
+        amrex::Real Gamma_eff = std::cbrt(2.0_rt) * scn_fac.z1 * scn_fac.z2 *
                          scn_fac.zs13inv * Gamma_e;
-        [[maybe_unused]] Real Gamma_eff_dT;
+        [[maybe_unused]] amrex::Real Gamma_eff_dT;
 
         if constexpr (do_T_derivatives) {
             Gamma_eff_dT = -Gamma_eff / state.temp;
@@ -895,18 +893,18 @@ void chabrier1998 (const plasma_state_t& state,
 
         // TAU/3, see Wallace1982, Eq. A2
 
-        Real tau12 = state.taufac * scn_fac.aznut;
+        amrex::Real tau12 = state.taufac * scn_fac.aznut;
 
-        [[maybe_unused]] Real tau12dT;
+        [[maybe_unused]] amrex::Real tau12dT;
         if constexpr (do_T_derivatives) {
             tau12dT = state.taufacdt * scn_fac.aznut;
         }
 
         // see Calder 2007 Eq. A8
 
-        Real b_fac = Gamma_eff / tau12;
+        amrex::Real b_fac = Gamma_eff / tau12;
 
-        [[maybe_unused]] Real b_fac_dT;
+        [[maybe_unused]] amrex::Real b_fac_dT;
         if constexpr (do_T_derivatives) {
             b_fac_dT = (Gamma_eff_dT - b_fac * tau12dT) / tau12;
         }
@@ -940,14 +938,14 @@ void chabrier1998 (const plasma_state_t& state,
     // is that we replaced the classical term which is f1 + f2 - f12
     // using results from Chabrier&Potekhin1998.
 
-    Real h12 = f1 + f2 - f12 + quantum_corr_1 + quantum_corr_2;
+    amrex::Real h12 = f1 + f2 - f12 + quantum_corr_1 + quantum_corr_2;
 
-    [[maybe_unused]] Real dh12dT;
+    [[maybe_unused]] amrex::Real dh12dT;
     if constexpr (do_T_derivatives) {
         dh12dT = f1dT + f2dT - f12dT + quantum_corr_1_dT + quantum_corr_2_dT;
     }
 
-    Real h12_max = 300.0_rt;
+    amrex::Real h12_max = 300.0_rt;
     h12 = amrex::min(h12_max, h12);
 
     scor = std::exp(h12);
@@ -966,7 +964,7 @@ template <int do_T_derivatives>
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 void actual_screen(const plasma_state_t& state,
                    const scrn::screen_factors_t& scn_fac,
-                   Real& scor, Real& scordt)
+                   amrex::Real& scor, amrex::Real& scordt)
 {
 #if SCREEN_METHOD == SCREEN_METHOD_null
     // null screening
@@ -987,9 +985,9 @@ void actual_screen(const plasma_state_t& state,
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 void actual_screen(const plasma_state_t& state,
                    const scrn::screen_factors_t& scn_fac,
-                   Real& scor, Real& scordt)
+                   amrex::Real& scor, amrex::Real& scordt)
 {
-    if (jacobian == 1) {
+    if (integrator_rp::jacobian == 1) {
         constexpr int do_T_derivatives = 1;
         actual_screen<do_T_derivatives>(state, scn_fac, scor, scordt);
     } else {
diff --git a/screening/screen_data.H b/screening/screen_data.H
index 571d730860..d81c8bc961 100644
--- a/screening/screen_data.H
+++ b/screening/screen_data.H
@@ -7,7 +7,7 @@
 #include <network.H>
 #include <microphysics_math.H>
 
-using namespace amrex;
+using namespace amrex::literals;
 
 namespace scrn {
     class screen_factors_t {