diff --git a/Source/CAMR.cpp b/Source/CAMR.cpp
index 54a7a94..6385318 100644
--- a/Source/CAMR.cpp
+++ b/Source/CAMR.cpp
@@ -184,8 +184,9 @@ CAMR::read_params()
   }
 
 #ifdef AMREX_USE_EB
-  if (do_mol == 0) {
-      amrex::Warning("EBGodunov is still a WIP");
+  // We only support PLM (not PPM) with using EB
+  if (do_mol == 0 && ppm_type != 0) {
+      amrex::Abort("Must use ppm_type = 0 when running with EBGodunov");
   }
 
   if ( (redistribution_type != "FluxRedist" ) &&
diff --git a/Source/Hydro/CAMR_construct_hydro_source.cpp b/Source/Hydro/CAMR_construct_hydro_source.cpp
index c749735..1d4c8c1 100644
--- a/Source/Hydro/CAMR_construct_hydro_source.cpp
+++ b/Source/Hydro/CAMR_construct_hydro_source.cpp
@@ -130,7 +130,15 @@ CAMR::construct_hydro_source (const MultiFab& S,
             const auto& src_in = sources_for_hydro.array(mfi);
             ParallelFor(
               qbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
-                hydro_srctoprim(i, j, k, qarr, qauxar, src_in, srcqarr, *lpmap);
+#ifdef AMREX_USE_EB
+                if (!flag_arr(i,j,k).isCovered()) {
+#endif
+                    hydro_srctoprim(i, j, k, qarr, qauxar, src_in, srcqarr, *lpmap);
+#ifdef AMREX_USE_EB
+                } else {
+                   for (int n=0; n<QVAR; n++) srcqarr(i,j,k,n) = 0.;
+                }
+#endif
               });
         }
 
diff --git a/Source/Hydro/Godunov/Godunov_2D_eb.cpp b/Source/Hydro/Godunov/Godunov_2D_eb.cpp
index 3ac6d90..eee3ac6 100644
--- a/Source/Hydro/Godunov/Godunov_2D_eb.cpp
+++ b/Source/Hydro/Godunov/Godunov_2D_eb.cpp
@@ -4,6 +4,7 @@
 #include "Hydro_utils_K.H"
 #include "flatten.H"
 #include "PLM.H"
+#include "PLM_eb.H"
 #include "PPM.H"
 
 #if (AMREX_SPACEDIM == 2)
@@ -107,8 +108,8 @@ Godunov_umeth_eb (
                 slope[n] = plm_slope_eb(i, j, k, n, 0, flag_arr, q, flat, iorder);
             }
           }
-          hydro_plm_d(i, j, k, 0, qxmarr, qxparr, slope, q, qaux(i, j, k, QC), dx, dt,
-                     small_dens, small_pres, *lpmap, apx);
+          hydro_plm_d_eb(i, j, k, 0, qxmarr, qxparr, slope, q, qaux(i, j, k, QC), dx, dt,
+                         small_dens, small_pres, *lpmap, apx);
 
           //
           // Y slopes and interp
@@ -121,29 +122,11 @@ Godunov_umeth_eb (
                   slope[n] = plm_slope_eb(i, j, k, n, 1, flag_arr, q, flat, iorder);
               }
           }
-          hydro_plm_d(i, j, k, 1, qymarr, qyparr, slope, q, qaux(i, j, k, QC), dy, dt,
-                     small_dens, small_pres, *lpmap, apy);
+          hydro_plm_d_eb(i, j, k, 1, qymarr, qyparr, slope, q, qaux(i, j, k, QC), dy, dt,
+                         small_dens, small_pres, *lpmap, apy);
         });
-
-  } else if (ppm_type == 1) {
-
-      // Compute the normal interface states by reconstructing
-      // the primitive variables using the piecewise parabolic method
-      // and doing characteristic tracing.  We do not apply the
-      // transverse terms here.
-
-      int idir = 0;
-      trace_ppm(
-        bxg2, idir, q, qaux, srcQ, qxmarr, qxparr, bxg2, dt, del, use_flattening,
-        small_dens, small_pres, lpmap);
-
-      idir = 1;
-      trace_ppm(
-        bxg2, idir, q, qaux, srcQ, qymarr, qyparr, bxg2, dt, del, use_flattening,
-        small_dens, small_pres, lpmap);
-
     } else {
-        amrex::Error("ppm_type must be 0 (PLM) or 1 (PPM)");
+        amrex::Error("ppm_type must be 0 (PLM) when using EB");
     }
 
   // *******************************************************************************
diff --git a/Source/Hydro/Godunov/Godunov_3D_eb.cpp b/Source/Hydro/Godunov/Godunov_3D_eb.cpp
index 869dd8a..2598615 100644
--- a/Source/Hydro/Godunov/Godunov_3D_eb.cpp
+++ b/Source/Hydro/Godunov/Godunov_3D_eb.cpp
@@ -5,6 +5,7 @@
 #include "Hydro_utils_K.H"
 #include "flatten.H"
 #include "PLM.H"
+#include "PLM_eb.H"
 #include "PPM.H"
 
 #if (AMREX_SPACEDIM == 3)
@@ -116,8 +117,8 @@ Godunov_umeth_eb (
               slope[n] = plm_slope_eb(i, j, k, n, 0, flag_arr, q, flat, iorder);
         }
         // cells -5,-5,-5
-        hydro_plm_d(i, j, k, idir, qxmarr, qxparr, slope, q, qaux(i, j, k, QC), dx, dt,
-                   small_dens, small_pres, *lpmap, apx);
+        hydro_plm_d_eb(i, j, k, idir, qxmarr, qxparr, slope, q, qaux(i, j, k, QC), dx, dt,
+                       small_dens, small_pres, *lpmap, apx);
 
         // Y slopes and interp
         idir = 1;
@@ -128,8 +129,8 @@ Godunov_umeth_eb (
               slope[n] = plm_slope_eb(i, j, k, n, 1, flag_arr, q, flat, iorder);
         }
         // cells -5,-5,-5
-        hydro_plm_d(i, j, k, idir, qymarr, qyparr, slope, q, qaux(i, j, k, QC), dy, dt,
-                   small_dens, small_pres, *lpmap, apy);
+        hydro_plm_d_eb(i, j, k, idir, qymarr, qyparr, slope, q, qaux(i, j, k, QC), dy, dt,
+                       small_dens, small_pres, *lpmap, apy);
 
         // Z slopes and interp
         idir = 2;
@@ -140,27 +141,11 @@ Godunov_umeth_eb (
               slope[n] = plm_slope_eb(i, j, k, n, 2, flag_arr, q, flat, iorder);
         }
         // cells -5,-5,-5
-        hydro_plm_d(i, j, k, idir, qzmarr, qzparr, slope, q, qaux(i, j, k, QC), dz, dt,
-                   small_dens, small_pres, *lpmap, apz);
+        hydro_plm_d_eb(i, j, k, idir, qzmarr, qzparr, slope, q, qaux(i, j, k, QC), dz, dt,
+                       small_dens, small_pres, *lpmap, apz);
       });
-  } else if (ppm_type == 1) {
-
-      // Compute the normal interface states by reconstructing
-      // the primitive variables using the piecewise parabolic method
-      // and doing characteristic tracing.  We do not apply the
-      // transverse terms here.
-
-      trace_ppm(bxg2, 0, q, qaux, srcQ, qxmarr, qxparr, bxg2, dt, del, use_flattening,
-                small_dens, small_pres, lpmap);
-
-      trace_ppm(bxg2, 1, q, qaux, srcQ, qymarr, qyparr, bxg2, dt, del, use_flattening,
-                small_dens, small_pres, lpmap);
-
-      trace_ppm(bxg2, 2, q, qaux, srcQ, qzmarr, qzparr, bxg2, dt, del, use_flattening,
-                small_dens, small_pres, lpmap);
-
   } else {
-      amrex::Error("ppm_type must be 0 (PLM) or 1 (PPM)");
+      amrex::Error("ppm_type must be 0 (PLM) when using EB");
   }
 
 
diff --git a/Source/Hydro/Godunov/Make.package b/Source/Hydro/Godunov/Make.package
index df1a20e..9f46eb2 100644
--- a/Source/Hydro/Godunov/Make.package
+++ b/Source/Hydro/Godunov/Make.package
@@ -9,4 +9,5 @@ CEXE_headers += PPM.H
 
 ifeq ($(USE_EB),TRUE)
 CEXE_sources += Godunov_$(DIM)D_eb.cpp
+CEXE_headers += PLM_eb.H
 endif
diff --git a/Source/Hydro/Godunov/PLM.H b/Source/Hydro/Godunov/PLM.H
index 0566872..7c8fccf 100644
--- a/Source/Hydro/Godunov/PLM.H
+++ b/Source/Hydro/Godunov/PLM.H
@@ -164,101 +164,6 @@ plm_pslope (
   }
 }
 
-
-#ifdef AMREX_USE_EB
-AMREX_GPU_DEVICE
-AMREX_FORCE_INLINE
-amrex::Real
-plm_slope_eb (
-  const int i,
-  const int j,
-  const int k,
-  const int n,
-  const int dir,
-  amrex::Array4<amrex::EBCellFlag const> const& flags,
-  amrex::Array4<const amrex::Real> const& q,
-  amrex::Real flat,
-  const int order)
-{
-  const amrex::IntVect iv{AMREX_D_DECL(i, j, k)};
-  const amrex::IntVect dvec(amrex::IntVect::TheDimensionVector(dir));
-  const amrex::IntVect ivm2(iv - 2 * dvec);
-  const amrex::IntVect ivm(iv - dvec);
-  const amrex::IntVect ivp(iv + dvec);
-  const amrex::IntVect ivp2(iv + 2 * dvec);
-
-  // We have enough cells in the x-direction to do 4th order slopes
-  // centered on (i,j,k)
-  if ( (order == 4) &&
-       !flags(iv).isCovered() &&
-       !flags(ivm).isCovered() &&
-       !flags(ivm2).isCovered() &&
-       !flags(ivp).isCovered() &&
-       !flags(ivp2).isCovered())
-  {
-    return plm_slope(i, j, k, n, dir, q, flat, order);
-  }
-  // We have enough cells in the x-direction to do 2nd order slopes
-  // centered on (i,j,k)
-  else if ( (order > 1) &&
-            !flags(iv).isCovered() &&
-            !flags(ivm).isCovered() &&
-            !flags(ivp).isCovered())
-  {
-    return plm_slope(i, j, k, n, dir, q, flat, 2);
-  } else {
-    return 0.0;
-  }
-}
-
-AMREX_GPU_DEVICE
-AMREX_FORCE_INLINE
-amrex::Real
-plm_pslope_eb (
-  const int i,
-  const int j,
-  const int k,
-  const int n,
-  const int dir,
-  amrex::Array4<amrex::EBCellFlag const> const& flags,
-  amrex::Array4<const amrex::Real> const& q,
-  const amrex::Real dx,
-  amrex::Array4<const amrex::Real> const& srcQ,
-  amrex::Real flat,
-  const int order)
-{
-  const amrex::IntVect iv{AMREX_D_DECL(i, j, k)};
-  const amrex::IntVect dvec(amrex::IntVect::TheDimensionVector(dir));
-  const amrex::IntVect ivm2(iv - 2 * dvec);
-  const amrex::IntVect ivm(iv - dvec);
-  const amrex::IntVect ivp(iv + dvec);
-  const amrex::IntVect ivp2(iv + 2 * dvec);
-
-  // We have enough cells in the x-direction to do 4th order slopes
-  // centered on (i,j,k) with all values at cell centers
-  if ((order == 4) &&
-      !flags(iv).isCovered() &&
-      !flags(ivm).isCovered() &&
-      !flags(ivm2).isCovered() &&
-      !flags(ivp).isCovered() &&
-      !flags(ivp2).isCovered())
-  {
-    return plm_pslope(i, j, k, n, dir, q, dx, srcQ, flat, order);
-  }
-  // We have enough cells in the x-direction to do 2nd order slopes
-  // centered on (i,j,k) with all values at cell centers
-  else if ((order > 1) &&
-           !flags(iv).isCovered() &&
-           !flags(ivm).isCovered() &&
-           !flags(ivp).isCovered())
-  {
-    return plm_pslope(i, j, k, n, dir, q, dx, srcQ, flat, 2);
-  } else {
-    return 0.0;
-  }
-}
-#endif
-
 AMREX_GPU_DEVICE
 AMREX_FORCE_INLINE
 void
@@ -276,11 +181,7 @@ hydro_plm_d (
   amrex::Real const dt,
   amrex::Real const small_dens,
   amrex::Real const small_pres,
-  PassMap const& pmap
-#ifdef AMREX_USE_EB
-  , amrex::Array4<const amrex::Real> const& area = {}
-#endif
-)
+  PassMap const& pmap)
 {
   amrex::ignore_unused(k);
   const amrex::IntVect iv{AMREX_D_DECL(i, j, k)};
@@ -314,9 +215,6 @@ hydro_plm_d (
     vel[dir] - cc, vel[dir], vel[dir] + cc};
 
   //  Construct the right state on the i-1/2 interface
-#ifdef AMREX_USE_EB
-  if(area(iv) > 0.0)
-#endif
   {
   const amrex::Real rho_ref =
     rho - 0.5 * (1.0 + dtdx * std::min(wv[0], 0.0)) * drho;
@@ -361,18 +259,8 @@ hydro_plm_d (
   qp(iv, QPRES) = std::max( qp(iv, QPRES), small_pres);
   qp(iv, QREINT) = rhoe_ref + (apright + amright) * enth * cs + azeright;
   }
-#ifdef AMREX_USE_EB
-  else {
-    for(int n = 0; n < QVAR; n++){
-      qp(iv, n) = 0.0;
-    }
-  }
-#endif
 
   //  Construct the left state on the i+1/2 interface
-#ifdef AMREX_USE_EB
-  if(area(ivp) > 0.0)
-#endif
   {
   const amrex::Real rho_ref =
     rho + 0.5 * (1.0 - dtdx * std::max(wv[2], 0.0)) * drho;
@@ -415,47 +303,22 @@ hydro_plm_d (
   qm(ivp, QPRES) = std::max( qm(ivp, QPRES), small_pres);
   qm(ivp, QREINT) = rhoe_ref + (apleft + amleft) * enth * cs + azeleft;
   }
-#ifdef AMREX_USE_EB
-  else {
-    for(int n = 0; n < QVAR; n++){
-      qm(ivp, n) = 0.0;
-    }
-  }
-#endif
 
   // Upwind the passive variables
   for (int ipassive = 0; ipassive < NPASSIVE; ++ipassive) {
     const int n = pmap.qpassMap[ipassive];
     const amrex::Real vel_ad = q(iv, QU + l_idx[0]);
     // Right state
-#ifdef AMREX_USE_EB
-    if (area(iv) > 0.0)
-#endif
     {
       const amrex::Real spzero = (vel_ad > 0) ? -1.0 : vel_ad * dtdx;
       const amrex::Real acmprght = 0.5 * (-1.0 - spzero) * slope[n];
       qp(iv, n) = q(iv, n) + acmprght;
     }
-#ifdef AMREX_USE_EB
-    else {
-      qp(iv, n) = 0.0;
-    }
-#endif
-
-#ifdef AMREX_USE_EB
-    if (area(iv) > 0.0)
-#endif
     {
       const amrex::Real spzero = vel_ad >= 0 ? vel_ad * dtdx : 1.0;
       const amrex::Real acmpleft = 0.5 * (1.0 - spzero) * slope[n];
       qm(ivp, n) = q(iv, n) + acmpleft;
     }
-#ifdef AMREX_USE_EB
-    else {
-      qm(ivp, n) = 0.0;
-    }
-#endif
   }
 }
-
 #endif
diff --git a/Source/Hydro/Godunov/PLM_eb.H b/Source/Hydro/Godunov/PLM_eb.H
new file mode 100644
index 0000000..225389d
--- /dev/null
+++ b/Source/Hydro/Godunov/PLM_eb.H
@@ -0,0 +1,288 @@
+#ifndef PLM_EB_H
+#define PLM_EB_H
+
+#ifdef AMREX_USE_EB
+
+#include <cmath>
+
+#include <AMReX_FArrayBox.H>
+#include <AMReX_IArrayBox.H>
+#include <AMReX_Box.H>
+#include <AMReX_EBCellFlag.H>
+
+#include "IndexDefines.H"
+
+AMREX_GPU_DEVICE
+AMREX_FORCE_INLINE
+amrex::Real
+plm_slope_eb (
+  const int i,
+  const int j,
+  const int k,
+  const int n,
+  const int dir,
+  amrex::Array4<amrex::EBCellFlag const> const& flags,
+  amrex::Array4<const amrex::Real> const& q,
+  amrex::Real flat,
+  const int order)
+{
+  const amrex::IntVect iv{AMREX_D_DECL(i, j, k)};
+  const amrex::IntVect dvec(amrex::IntVect::TheDimensionVector(dir));
+  const amrex::IntVect ivm2(iv - 2 * dvec);
+  const amrex::IntVect ivm(iv - dvec);
+  const amrex::IntVect ivp(iv + dvec);
+  const amrex::IntVect ivp2(iv + 2 * dvec);
+
+  // We have enough cells in the x-direction to do 4th order slopes
+  // centered on (i,j,k)
+  if ( (order == 4) &&
+       !flags(iv).isCovered() &&
+       !flags(ivm).isCovered() &&
+       !flags(ivm2).isCovered() &&
+       !flags(ivp).isCovered() &&
+       !flags(ivp2).isCovered())
+  {
+    return plm_slope(i, j, k, n, dir, q, flat, order);
+  }
+  // We have enough cells in the x-direction to do 2nd order slopes
+  // centered on (i,j,k)
+  else if ( (order > 1) &&
+            !flags(iv).isCovered() &&
+            !flags(ivm).isCovered() &&
+            !flags(ivp).isCovered())
+  {
+    return plm_slope(i, j, k, n, dir, q, flat, 2);
+  } else {
+    return 0.0;
+  }
+}
+
+AMREX_GPU_DEVICE
+AMREX_FORCE_INLINE
+amrex::Real
+plm_pslope_eb (
+  const int i,
+  const int j,
+  const int k,
+  const int n,
+  const int dir,
+  amrex::Array4<amrex::EBCellFlag const> const& flags,
+  amrex::Array4<const amrex::Real> const& q,
+  const amrex::Real dx,
+  amrex::Array4<const amrex::Real> const& srcQ,
+  amrex::Real flat,
+  const int order)
+{
+  const amrex::IntVect iv{AMREX_D_DECL(i, j, k)};
+  const amrex::IntVect dvec(amrex::IntVect::TheDimensionVector(dir));
+  const amrex::IntVect ivm2(iv - 2 * dvec);
+  const amrex::IntVect ivm(iv - dvec);
+  const amrex::IntVect ivp(iv + dvec);
+  const amrex::IntVect ivp2(iv + 2 * dvec);
+
+  // We have enough cells in the x-direction to do 4th order slopes
+  // centered on (i,j,k) with all values at cell centers
+  if ((order == 4) &&
+      !flags(iv).isCovered() &&
+      !flags(ivm).isCovered() &&
+      !flags(ivm2).isCovered() &&
+      !flags(ivp).isCovered() &&
+      !flags(ivp2).isCovered())
+  {
+    return plm_pslope(i, j, k, n, dir, q, dx, srcQ, flat, order);
+  }
+  // We have enough cells in the x-direction to do 2nd order slopes
+  // centered on (i,j,k) with all values at cell centers
+  else if ((order > 1) &&
+           !flags(iv).isCovered() &&
+           !flags(ivm).isCovered() &&
+           !flags(ivp).isCovered())
+  {
+    return plm_pslope(i, j, k, n, dir, q, dx, srcQ, flat, 2);
+  } else {
+    return 0.0;
+  }
+}
+#endif
+
+AMREX_GPU_DEVICE
+AMREX_FORCE_INLINE
+void
+hydro_plm_d_eb (
+  const int i,
+  const int j,
+  const int k,
+  const int dir,
+  amrex::Array4<amrex::Real> const& qm,
+  amrex::Array4<amrex::Real> const& qp,
+  amrex::Real const slope[QVAR],
+  amrex::Array4<const amrex::Real> const& q,
+  amrex::Real const cc,
+  amrex::Real const dx,
+  amrex::Real const dt,
+  amrex::Real const small_dens,
+  amrex::Real const small_pres,
+  PassMap const& pmap
+  , amrex::Array4<const amrex::Real> const& area = {}
+)
+{
+  amrex::ignore_unused(k);
+  const amrex::IntVect iv{AMREX_D_DECL(i, j, k)};
+  const amrex::IntVect ivp(iv + amrex::IntVect::TheDimensionVector(dir));
+  const amrex::GpuArray<const int, 3> bdim{{dir == 0, dir == 1, dir == 2}};
+  const amrex::GpuArray<const int, 3> l_idx{
+    {bdim[0] * 0 + bdim[1] * 1 + bdim[2] * 2,
+     bdim[0] * 1 + bdim[1] * 0 + bdim[2] * 0,
+     bdim[0] * 2 + bdim[1] * 2 + bdim[2] * 1}};
+
+  const amrex::Real dtdx = dt / dx;
+  const amrex::Real cs = cc * cc;
+  const amrex::Real rho = q(iv, QRHO);
+  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> vel{
+    AMREX_D_DECL(q(iv, QU), q(iv, QV), q(iv, QW))};
+  const amrex::Real p = q(iv, QPRES);
+  const amrex::Real rhoe = q(iv, QREINT);
+  const amrex::Real enth = ((rhoe + p) / rho) / cs;
+  const amrex::Real drho = slope[QRHO];
+  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>
+     dvel{AMREX_D_DECL(slope[QU], slope[QV], slope[QW])};
+  const amrex::Real dp = slope[QPRES];
+  const amrex::Real drhoe = slope[QREINT];
+  const amrex::Real alpham = 0.5 * (dp / (rho * cc) - dvel[dir]) * rho / cc;
+  const amrex::Real alphap = 0.5 * (dp / (rho * cc) + dvel[dir]) * rho / cc;
+  const amrex::Real alpha0r = drho - dp / cs;
+  const amrex::Real alpha0e = drhoe - dp * enth;
+  AMREX_D_TERM(, const amrex::Real alpha0v = dvel[l_idx[1]];
+               , const amrex::Real alpha0w = dvel[l_idx[2]];)
+  const amrex::GpuArray<amrex::Real, 3> wv = {
+    vel[dir] - cc, vel[dir], vel[dir] + cc};
+
+  //  Construct the right state on the i-1/2 interface
+  if (area(iv) > 0.0)
+  {
+  const amrex::Real rho_ref =
+    rho - 0.5 * (1.0 + dtdx * std::min(wv[0], 0.0)) * drho;
+  amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> vel_ref{
+    AMREX_D_DECL(0.0, 0.0, 0.0)};
+  for (int dim = 0; dim < AMREX_SPACEDIM; dim++) {
+    vel_ref[dim] =
+      vel[dim] -
+      0.5 * (1.0 + dtdx * std::min(wv[0], 0.0)) * dvel[dim];
+  }
+  const amrex::Real p_ref =
+    p - 0.5 * (1.0 + dtdx * std::min(wv[0], 0.0)) * dp;
+  const amrex::Real rhoe_ref =
+    rhoe - 0.5 * (1.0 + dtdx * std::min(wv[0], 0.0)) * drhoe;
+
+  const amrex::Real apright = 0.25 * dtdx * (wv[0] - wv[2]) *
+                              (1.0 - amrex::Math::copysign(1.0, wv[2])) *
+                              alphap;
+  const amrex::Real amright = 0.0;
+
+  const amrex::Real azrright = 0.25 * dtdx * (wv[0] - wv[1]) *
+                               (1.0 - amrex::Math::copysign(1.0, wv[1])) *
+                               alpha0r;
+  const amrex::Real azeright = 0.25 * dtdx * (wv[0] - wv[1]) *
+                               (1.0 - amrex::Math::copysign(1.0, wv[1])) *
+                               alpha0e;
+  AMREX_D_TERM(, const amrex::Real azv1rght =
+                   0.25 * dtdx * (wv[0] - wv[1]) *
+                   (1.0 - amrex::Math::copysign(1.0, wv[1])) * alpha0v;
+               , const amrex::Real azw1rght =
+                   0.25 * dtdx * (wv[0] - wv[1]) *
+                   (1.0 - amrex::Math::copysign(1.0, wv[1])) * alpha0w;)
+
+  qp(iv, QRHO) = rho_ref + apright + amright + azrright;
+  qp(iv, QRHO) = std::max( qp(iv, QRHO), small_dens);
+  AMREX_D_TERM(qp(iv, QU + l_idx[0]) =
+                 vel_ref[l_idx[0]] + (apright - amright) * cc / rho;
+               qp(iv, QU + l_idx[1]) = 0.; qp(iv, QU + l_idx[2]) = 0.;
+               , qp(iv, QU + l_idx[1]) = vel_ref[l_idx[1]] + azv1rght;
+               , qp(iv, QU + l_idx[2]) = vel_ref[l_idx[2]] + azw1rght;);
+  qp(iv, QPRES) = p_ref + (apright + amright) * cs;
+  qp(iv, QPRES) = std::max( qp(iv, QPRES), small_pres);
+  qp(iv, QREINT) = rhoe_ref + (apright + amright) * enth * cs + azeright;
+  }
+  else {
+    for(int n = 0; n < QVAR; n++){
+      qp(iv, n) = 0.0;
+    }
+  }
+
+  //  Construct the left state on the i+1/2 interface
+  if(area(ivp) > 0.0)
+  {
+  const amrex::Real rho_ref =
+    rho + 0.5 * (1.0 - dtdx * std::max(wv[2], 0.0)) * drho;
+  amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> vel_ref{
+    AMREX_D_DECL(0.0, 0.0, 0.0)};
+  for (int dim = 0; dim < AMREX_SPACEDIM; dim++) {
+    vel_ref[dim] =
+      vel[dim] +
+      0.5 * (1.0 - dtdx * std::max(wv[2], 0.0)) * dvel[dim];
+  }
+  const amrex::Real p_ref = p + 0.5 * (1.0 - dtdx * std::max(wv[2], 0.0)) * dp;
+
+  const amrex::Real rhoe_ref =
+    rhoe + 0.5 * (1.0 - dtdx * std::max(wv[2], 0.0)) * drhoe;
+
+  const amrex::Real apleft = 0.0;
+  const amrex::Real amleft = 0.25 * dtdx * (wv[2] - wv[0]) *
+                             (1.0 + amrex::Math::copysign(1.0, wv[0])) * alpham;
+
+  const amrex::Real azrleft = 0.25 * dtdx * (wv[2] - wv[1]) *
+                              (1.0 + amrex::Math::copysign(1.0, wv[1])) *
+                              alpha0r;
+  const amrex::Real azeleft = 0.25 * dtdx * (wv[2] - wv[1]) *
+                              (1.0 + amrex::Math::copysign(1.0, wv[1])) *
+                              alpha0e;
+  AMREX_D_TERM(, const amrex::Real azv1left =
+                   0.25 * dtdx * (wv[2] - wv[1]) *
+                   (1.0 + amrex::Math::copysign(1.0, wv[1])) * alpha0v;
+               , const amrex::Real azw1left =
+                   0.25 * dtdx * (wv[2] - wv[1]) *
+                   (1.0 + amrex::Math::copysign(1.0, wv[1])) * alpha0w;)
+  qm(ivp, QRHO) = rho_ref + apleft + amleft + azrleft;
+  qm(ivp, QRHO) = std::max( qm(ivp, QRHO), small_dens);
+  AMREX_D_TERM(qm(ivp, QU + l_idx[0]) =
+                 vel_ref[l_idx[0]] + (apleft - amleft) * cc / rho;
+               qm(ivp, QU + l_idx[1]) = 0.; qm(ivp, QU + l_idx[2]) = 0.;
+               , qm(ivp, QU + l_idx[1]) = vel_ref[l_idx[1]] + azv1left;
+               , qm(ivp, QU + l_idx[2]) = vel_ref[l_idx[2]] + azw1left;);
+  qm(ivp, QPRES) = p_ref + (apleft + amleft) * cs;
+  qm(ivp, QPRES) = std::max( qm(ivp, QPRES), small_pres);
+  qm(ivp, QREINT) = rhoe_ref + (apleft + amleft) * enth * cs + azeleft;
+  }
+  else {
+    for(int n = 0; n < QVAR; n++){
+      qm(ivp, n) = 0.0;
+    }
+  }
+
+  // Upwind the passive variables
+  for (int ipassive = 0; ipassive < NPASSIVE; ++ipassive) {
+    const int n = pmap.qpassMap[ipassive];
+    const amrex::Real vel_ad = q(iv, QU + l_idx[0]);
+    // Right state
+    if (area(iv) > 0.0)
+    {
+      const amrex::Real spzero = (vel_ad > 0) ? -1.0 : vel_ad * dtdx;
+      const amrex::Real acmprght = 0.5 * (-1.0 - spzero) * slope[n];
+      qp(iv, n) = q(iv, n) + acmprght;
+    }
+    else {
+      qp(iv, n) = 0.0;
+    }
+
+    if (area(iv) > 0.0)
+    {
+      const amrex::Real spzero = vel_ad >= 0 ? vel_ad * dtdx : 1.0;
+      const amrex::Real acmpleft = 0.5 * (1.0 - spzero) * slope[n];
+      qm(ivp, n) = q(iv, n) + acmpleft;
+    }
+    else {
+      qm(ivp, n) = 0.0;
+    }
+  }
+}
+#endif