diff --git a/Exec/Production/CounterFlow/PeleLMeX_PatchFlowVariables.cpp b/Exec/Production/CounterFlow/PeleLMeX_PatchFlowVariables.cpp
index 374c7ea3..cebcc6db 100644
--- a/Exec/Production/CounterFlow/PeleLMeX_PatchFlowVariables.cpp
+++ b/Exec/Production/CounterFlow/PeleLMeX_PatchFlowVariables.cpp
@@ -8,27 +8,26 @@ patchFlowVariables(
 {
 
   amrex::Print() << "\nPatching flow variables..";
-  
+
   const amrex::Real* prob_lo = geom.ProbLo();
   const amrex::Real* prob_hi = geom.ProbHi();
-  const amrex::Real* dx      = geom.CellSize();
+  const amrex::Real* dx = geom.CellSize();
 
   for (amrex::MFIter mfi(a_mf, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
     const amrex::Box& bx = mfi.tilebox();
-    
-	auto const& rho_arr = a_mf.array(mfi, DENSITY);
+
+    auto const& rho_arr = a_mf.array(mfi, DENSITY);
     auto const& rhoY_arr = a_mf.array(mfi, FIRSTSPEC);
     auto const& rhoH_arr = a_mf.array(mfi, RHOH);
-	auto const& temp_arr = a_mf.array(mfi, TEMP);
-	Real massfrac[NUM_SPECIES] = {0.0};
-    
-	amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
-      
-	  auto eos = pele::physics::PhysicsType::eos();
+    auto const& temp_arr = a_mf.array(mfi, TEMP);
+    Real massfrac[NUM_SPECIES] = {0.0};
+
+    amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+      auto eos = pele::physics::PhysicsType::eos();
       amrex::Real massfrac[NUM_SPECIES] = {0.0};
       amrex::Real sumYs = 0.0;
-	  
-	  for (int n = 0; n < NUM_SPECIES; n++) {
+
+      for (int n = 0; n < NUM_SPECIES; n++) {
         massfrac[n] = rhoY_arr(i, j, k, n);
 #ifdef N2_ID
         if (n != N2_ID) {
@@ -39,28 +38,29 @@ patchFlowVariables(
 #ifdef N2_ID
       massfrac[N2_ID] = 1.0 - sumYs;
 #endif
-       
-	  AMREX_D_TERM(const amrex::Real x = prob_lo[0] + (i+0.5)*dx[0];,
-               const amrex::Real y = prob_lo[1] + (j+0.5)*dx[1];,
-               const amrex::Real z = prob_lo[2] + (k+0.5)*dx[2];);
 
-      AMREX_D_TERM(const amrex::Real Lx = prob_hi[0] - prob_lo[0];,
-               const amrex::Real Ly = prob_hi[1] - prob_lo[1];,
-               const amrex::Real Lz = prob_hi[2] - prob_lo[2]);
+      AMREX_D_TERM(const amrex::Real x = prob_lo[0] + (i + 0.5) * dx[0];
+                   , const amrex::Real y = prob_lo[1] + (j + 0.5) * dx[1];
+                   , const amrex::Real z = prob_lo[2] + (k + 0.5) * dx[2];);
+
+      AMREX_D_TERM(const amrex::Real Lx = prob_hi[0] - prob_lo[0];
+                   , const amrex::Real Ly = prob_hi[1] - prob_lo[1];
+                   , const amrex::Real Lz = prob_hi[2] - prob_lo[2]);
+
+      AMREX_D_TERM(const amrex::Real xc = prob_lo[0] + Lx / 2.0;
+                   , const amrex::Real yc = prob_lo[1] + Ly / 2.0;
+                   , const amrex::Real zc = prob_lo[2] + Lz / 2.0);
+
+      amrex::Real radiusSq = AMREX_D_TERM(
+        (x - xc) * (x - xc), +(y - yc) * (y - yc), +(z - zc) * (z - zc));
 
-      AMREX_D_TERM(const amrex::Real xc = prob_lo[0] + Lx/2.0;,
-               const amrex::Real yc = prob_lo[1] + Ly/2.0;,
-               const amrex::Real zc = prob_lo[2] + Lz/2.0);
-	  
-	  amrex::Real radiusSq = AMREX_D_TERM(  (x-xc) * (x-xc),
-                                           + (y-yc) * (y-yc),
-                                           + (z-zc) * (z-zc));
-									   
       amrex::Real radius = std::sqrt(radiusSq);
-	  amrex::Real mixingWidth = 0.1*lprobparm.ignitSphereRad;
-	  amrex::Real mixingFunction = 0.5 * ( 1.0 + std::tanh((lprobparm.ignitSphereRad - radius)/mixingWidth));
-	  temp_arr(i, j, k) = mixingFunction * lprobparm.ignitT + (1.0 - mixingFunction) * lprobparm.T_inert;
-	  
+      amrex::Real mixingWidth = 0.1 * lprobparm.ignitSphereRad;
+      amrex::Real mixingFunction =
+        0.5 *
+        (1.0 + std::tanh((lprobparm.ignitSphereRad - radius) / mixingWidth));
+      temp_arr(i, j, k) = mixingFunction * lprobparm.ignitT +
+                          (1.0 - mixingFunction) * lprobparm.T_inert;
     });
   }
 
diff --git a/Exec/Production/CounterFlow/pelelmex_prob.H b/Exec/Production/CounterFlow/pelelmex_prob.H
index ba84b8f2..a83d8929 100644
--- a/Exec/Production/CounterFlow/pelelmex_prob.H
+++ b/Exec/Production/CounterFlow/pelelmex_prob.H
@@ -12,92 +12,105 @@
 
 AMREX_GPU_DEVICE
 AMREX_FORCE_INLINE
-void pelelmex_initdata(int i, int j, int k,
-                     int is_incompressible,
-                     amrex::Array4<amrex::Real> const& state,
-                     amrex::Array4<amrex::Real> const& aux,
-                     amrex::GeometryData const& geomdata,
-                     ProbParm const& prob_parm,
-                     pele::physics::PMF::PmfData::DataContainer const * pmf_data)
+void
+pelelmex_initdata(
+  int i,
+  int j,
+  int k,
+  int is_incompressible,
+  amrex::Array4<amrex::Real> const& state,
+  amrex::Array4<amrex::Real> const& aux,
+  amrex::GeometryData const& geomdata,
+  ProbParm const& prob_parm,
+  pele::physics::PMF::PmfData::DataContainer const* pmf_data)
 {
-    const amrex::Real* prob_lo = geomdata.ProbLo();
-    const amrex::Real* prob_hi = geomdata.ProbHi();
-    const amrex::Real* dx      = geomdata.CellSize();
-
-    AMREX_D_TERM(const amrex::Real x = prob_lo[0] + (i+0.5)*dx[0];,
-                 const amrex::Real y = prob_lo[1] + (j+0.5)*dx[1];,
-                 const amrex::Real z = prob_lo[2] + (k+0.5)*dx[2];);
-
-    AMREX_D_TERM(const amrex::Real Lx = prob_hi[0] - prob_lo[0];,
-                 const amrex::Real Ly = prob_hi[1] - prob_lo[1];,
-                 const amrex::Real Lz = prob_hi[2] - prob_lo[2]);
-
-    AMREX_D_TERM(const amrex::Real xc = prob_lo[0] + Lx/2.0;,
-                 const amrex::Real yc = prob_lo[1] + Ly/2.0;,
-                 const amrex::Real zc = prob_lo[2] + Lz/2.0);
-
-    constexpr amrex::Real Pi = 3.14159265358979323846264338327950288;
-
-    auto eos = pele::physics::PhysicsType::eos();
-    amrex::Real massfrac[NUM_SPECIES] = {0.0};
-    amrex::Real massfrac_cold[NUM_SPECIES] = {0.0};
-    amrex::Real massfrac_hot[NUM_SPECIES] = {0.0};
-
-    massfrac_cold[NC12H26_ID] = 0.05;
-	massfrac_cold[H2O_ID] = 0.0;
-	massfrac_cold[CO2_ID] = 0.0;
-	massfrac_cold[O2_ID] = 0.233 * (1.0 - massfrac_cold[NC12H26_ID]);;
-    massfrac_cold[N2_ID] = 1.0 - massfrac_cold[NC12H26_ID] - massfrac_cold[O2_ID];
-
-	state(i,j,k,TEMP) = prob_parm.T_inert;
-
-    if (prob_parm.do_ignit) {
-       amrex::Real radiusSq = AMREX_D_TERM(  (x-xc) * (x-xc),
-                                           + (y-yc) * (y-yc),
-                                           + (z-zc) * (z-zc));
-       amrex::Real radius = std::sqrt(radiusSq);
-       amrex::Real mixingWidth = 0.1*prob_parm.ignitSphereRad;
-       amrex::Real mixingFunction = 0.5 * ( 1.0 + std::tanh((prob_parm.ignitSphereRad - radius)/mixingWidth));
-	   state(i,j,k,TEMP) = mixingFunction * prob_parm.ignitT + (1.0 - mixingFunction) * prob_parm.T_inert;
-       massfrac_hot[NC12H26_ID] = 0.00;
-	   massfrac_hot[H2O_ID] = prob_parm.Y_H2O_ign;
-	   massfrac_hot[CO2_ID] = prob_parm.Y_CO2_ign;
-	   massfrac_hot[O2_ID] = 0.233 * (1.0 - massfrac_hot[H2O_ID] - massfrac_hot[CO2_ID]);
-	   massfrac_hot[N2_ID] = 1.0 - massfrac_hot[H2O_ID] - massfrac_hot[CO2_ID] - massfrac_hot[O2_ID];
-
-       massfrac[NC12H26_ID] = mixingFunction * massfrac_hot[NC12H26_ID] + (1.0-mixingFunction) * massfrac_cold[NC12H26_ID];
-       massfrac[H2O_ID] = mixingFunction * massfrac_hot[H2O_ID] + (1.0-mixingFunction) * massfrac_cold[H2O_ID];
-	   massfrac[CO2_ID] = mixingFunction * massfrac_hot[CO2_ID] + (1.0-mixingFunction) * massfrac_cold[CO2_ID];
-	   massfrac[O2_ID] = mixingFunction * massfrac_hot[O2_ID] + (1.0-mixingFunction) * massfrac_cold[O2_ID];
-       massfrac[N2_ID] = mixingFunction * massfrac_hot[N2_ID] + (1.0-mixingFunction) * massfrac_cold[N2_ID];
-    } else {
-       massfrac[NC12H26_ID] = massfrac_cold[NC12H26_ID];
-	   massfrac[H2O_ID] = massfrac_cold[H2O_ID];
-	   massfrac[CO2_ID] = massfrac_cold[CO2_ID];
-       massfrac[O2_ID] = massfrac_cold[O2_ID];
-       massfrac[N2_ID] = massfrac_cold[N2_ID];
-    }
-
-    AMREX_D_TERM(state(i,j,k,VELX) = 0.0;,
-                 state(i,j,k,VELY) = 0.0;,
-                 state(i,j,k,VELZ) = 0.0);
-
-    amrex::Real P_cgs = prob_parm.P_mean * 10.0;
-
-    // Density
-    amrex::Real rho_cgs = 0.0;
-    eos.PYT2R(P_cgs, massfrac, state(i,j,k,TEMP), rho_cgs);
-    state(i,j,k,DENSITY) = rho_cgs * 1.0e3;
-
-    // Enthalpy
-    amrex::Real h_cgs = 0.0;
-    eos.TY2H(state(i,j,k,TEMP), massfrac, h_cgs);
-    state(i,j,k,RHOH) = h_cgs * 1.0e-4 * state(i,j,k,DENSITY);
-
-    // Species mass
-    for (int n = 0; n < NUM_SPECIES; n++) {
-       state(i,j,k,FIRSTSPEC+n) = massfrac[n] * state(i,j,k,DENSITY);
-    }
+  const amrex::Real* prob_lo = geomdata.ProbLo();
+  const amrex::Real* prob_hi = geomdata.ProbHi();
+  const amrex::Real* dx = geomdata.CellSize();
+
+  AMREX_D_TERM(const amrex::Real x = prob_lo[0] + (i + 0.5) * dx[0];
+               , const amrex::Real y = prob_lo[1] + (j + 0.5) * dx[1];
+               , const amrex::Real z = prob_lo[2] + (k + 0.5) * dx[2];);
+
+  AMREX_D_TERM(const amrex::Real Lx = prob_hi[0] - prob_lo[0];
+               , const amrex::Real Ly = prob_hi[1] - prob_lo[1];
+               , const amrex::Real Lz = prob_hi[2] - prob_lo[2]);
+
+  AMREX_D_TERM(const amrex::Real xc = prob_lo[0] + Lx / 2.0;
+               , const amrex::Real yc = prob_lo[1] + Ly / 2.0;
+               , const amrex::Real zc = prob_lo[2] + Lz / 2.0);
+
+  constexpr amrex::Real Pi = 3.14159265358979323846264338327950288;
+
+  auto eos = pele::physics::PhysicsType::eos();
+  amrex::Real massfrac[NUM_SPECIES] = {0.0};
+  amrex::Real massfrac_cold[NUM_SPECIES] = {0.0};
+  amrex::Real massfrac_hot[NUM_SPECIES] = {0.0};
+
+  massfrac_cold[NC12H26_ID] = 0.05;
+  massfrac_cold[H2O_ID] = 0.0;
+  massfrac_cold[CO2_ID] = 0.0;
+  massfrac_cold[O2_ID] = 0.233 * (1.0 - massfrac_cold[NC12H26_ID]);
+  ;
+  massfrac_cold[N2_ID] = 1.0 - massfrac_cold[NC12H26_ID] - massfrac_cold[O2_ID];
+
+  state(i, j, k, TEMP) = prob_parm.T_inert;
+
+  if (prob_parm.do_ignit) {
+    amrex::Real radiusSq = AMREX_D_TERM(
+      (x - xc) * (x - xc), +(y - yc) * (y - yc), +(z - zc) * (z - zc));
+    amrex::Real radius = std::sqrt(radiusSq);
+    amrex::Real mixingWidth = 0.1 * prob_parm.ignitSphereRad;
+    amrex::Real mixingFunction =
+      0.5 *
+      (1.0 + std::tanh((prob_parm.ignitSphereRad - radius) / mixingWidth));
+    state(i, j, k, TEMP) = mixingFunction * prob_parm.ignitT +
+                           (1.0 - mixingFunction) * prob_parm.T_inert;
+    massfrac_hot[NC12H26_ID] = 0.00;
+    massfrac_hot[H2O_ID] = prob_parm.Y_H2O_ign;
+    massfrac_hot[CO2_ID] = prob_parm.Y_CO2_ign;
+    massfrac_hot[O2_ID] =
+      0.233 * (1.0 - massfrac_hot[H2O_ID] - massfrac_hot[CO2_ID]);
+    massfrac_hot[N2_ID] =
+      1.0 - massfrac_hot[H2O_ID] - massfrac_hot[CO2_ID] - massfrac_hot[O2_ID];
+
+    massfrac[NC12H26_ID] = mixingFunction * massfrac_hot[NC12H26_ID] +
+                           (1.0 - mixingFunction) * massfrac_cold[NC12H26_ID];
+    massfrac[H2O_ID] = mixingFunction * massfrac_hot[H2O_ID] +
+                       (1.0 - mixingFunction) * massfrac_cold[H2O_ID];
+    massfrac[CO2_ID] = mixingFunction * massfrac_hot[CO2_ID] +
+                       (1.0 - mixingFunction) * massfrac_cold[CO2_ID];
+    massfrac[O2_ID] = mixingFunction * massfrac_hot[O2_ID] +
+                      (1.0 - mixingFunction) * massfrac_cold[O2_ID];
+    massfrac[N2_ID] = mixingFunction * massfrac_hot[N2_ID] +
+                      (1.0 - mixingFunction) * massfrac_cold[N2_ID];
+  } else {
+    massfrac[NC12H26_ID] = massfrac_cold[NC12H26_ID];
+    massfrac[H2O_ID] = massfrac_cold[H2O_ID];
+    massfrac[CO2_ID] = massfrac_cold[CO2_ID];
+    massfrac[O2_ID] = massfrac_cold[O2_ID];
+    massfrac[N2_ID] = massfrac_cold[N2_ID];
+  }
+
+  AMREX_D_TERM(state(i, j, k, VELX) = 0.0;, state(i, j, k, VELY) = 0.0;
+               , state(i, j, k, VELZ) = 0.0);
+
+  amrex::Real P_cgs = prob_parm.P_mean * 10.0;
+
+  // Density
+  amrex::Real rho_cgs = 0.0;
+  eos.PYT2R(P_cgs, massfrac, state(i, j, k, TEMP), rho_cgs);
+  state(i, j, k, DENSITY) = rho_cgs * 1.0e3;
+
+  // Enthalpy
+  amrex::Real h_cgs = 0.0;
+  eos.TY2H(state(i, j, k, TEMP), massfrac, h_cgs);
+  state(i, j, k, RHOH) = h_cgs * 1.0e-4 * state(i, j, k, DENSITY);
+
+  // Species mass
+  for (int n = 0; n < NUM_SPECIES; n++) {
+    state(i, j, k, FIRSTSPEC + n) = massfrac[n] * state(i, j, k, DENSITY);
+  }
 }
 
 AMREX_GPU_DEVICE
@@ -112,117 +125,125 @@ bcnormal(
   const amrex::Real time,
   amrex::GeometryData const& geomdata,
   ProbParm const& prob_parm,
-  pele::physics::PMF::PmfData::DataContainer const * /*pmf_data*/)
+  pele::physics::PMF::PmfData::DataContainer const* /*pmf_data*/)
 {
-   const amrex::Real* prob_lo = geomdata.ProbLo();
-   const amrex::Real* prob_hi = geomdata.ProbHi();
-
-   AMREX_D_TERM(const amrex::Real Lx = prob_hi[0] - prob_lo[0];,
-                const amrex::Real Ly = prob_hi[1] - prob_lo[1];,
-                const amrex::Real Lz = prob_hi[2] - prob_lo[2]);
-
-   amrex::Real splitx = prob_lo[0] + 0.5 * Lx;
-   amrex::Real splity = prob_lo[1] + 0.5 * Ly;
-
-   amrex::Real massfrac_ox[NUM_SPECIES] = {0.0};
-   amrex::Real massfrac_fuel[NUM_SPECIES] = {0.0};
-   amrex::Real massfrac_inert[NUM_SPECIES] = {0.0};
-   amrex::Real massfrac_mix[NUM_SPECIES] = {0.0};
-
-   massfrac_ox[N2_ID] = prob_parm.Y_N2_ox;
-   massfrac_ox[O2_ID] = prob_parm.Y_O2_ox;
-   massfrac_fuel[NC12H26_ID] = prob_parm.Y_fuel; //Replace ID with gas fuel
-   massfrac_fuel[N2_ID] = prob_parm.Y_N2_fuel;
-   massfrac_inert[N2_ID] = 1.0;
-
-   auto eos = pele::physics::PhysicsType::eos();
-
-   for (int n = 0; n < NVAR; n++){
-      s_ext[n] = 0.0;
-   }
-
-// Get zone
-   int zone = -1;
-   if ( x[0] < splitx ) {
-      zone = 1; // Oxidizer side jet
-   } else {
-      zone = 2; // Fuel side jet
-   }
-   amrex::Real radius = std::abs(x[1]-splity);
-   amrex::Real mixingWidth = 0.1*prob_parm.jetRadius;
-   amrex::Real mixingFunction = 0.5 * ( 1.0 + std::tanh((prob_parm.jetRadius - radius)/mixingWidth));
-
-   amrex::Real P_cgs = prob_parm.P_mean * 10.0;
-   amrex::Real rho_cgs, RhoH_temp;
-   amrex::Real rhoOx = -1.0;
-   {
-      eos.PYT2R(P_cgs, massfrac_ox, prob_parm.T_ox, rhoOx);
-      rhoOx *= 1.0e3;
-   }
-   amrex::Real rhoFuel = -1.0;
-   {
-      eos.PYT2R(P_cgs, massfrac_fuel, prob_parm.T_fuel, rhoFuel);
-      rhoFuel *= 1.0e3;
-   }
-   if ( zone == 1 ) {
-	  s_ext[TEMP] = prob_parm.T_ox * mixingFunction +
-                    prob_parm.T_inert * (1.0 - mixingFunction);
-      for (int n = 0; n < NUM_SPECIES; n++) {
-         massfrac_mix[n] = massfrac_ox[n] * mixingFunction +
-                           massfrac_inert[n] * (1.0 - mixingFunction);
-      }
-      eos.PYT2R(P_cgs, massfrac_mix, s_ext[TEMP], rho_cgs);
-      s_ext[DENSITY] = rho_cgs * 1.0e3;
-
-      eos.TY2H(s_ext[TEMP], massfrac_mix, RhoH_temp);
-      s_ext[RHOH] = RhoH_temp * 1.0e-4 * s_ext[DENSITY];
-
-      for (int n = 0; n < NUM_SPECIES; n++) {
-         s_ext[FIRSTSPEC+n] = massfrac_mix[n] * s_ext[DENSITY];
-      }
-   } else if ( zone == 2 ) {
-	  s_ext[TEMP] = prob_parm.T_fuel * mixingFunction +
-                    prob_parm.T_inert * (1.0 - mixingFunction);
-      for (int n = 0; n < NUM_SPECIES; n++) {
-         massfrac_mix[n] = massfrac_fuel[n] * mixingFunction +
-                           massfrac_inert[n] * (1.0 - mixingFunction);
-      }
-      eos.PYT2R(P_cgs, massfrac_mix, s_ext[TEMP], rho_cgs);
-      s_ext[DENSITY] = rho_cgs * 1.0e3;
-
-      eos.TY2H(s_ext[TEMP], massfrac_mix, RhoH_temp);
-      s_ext[RHOH] = RhoH_temp * 1.0e-4 * s_ext[DENSITY];
-
-      for (int n = 0; n < NUM_SPECIES; n++) {
-         s_ext[FIRSTSPEC+n] = massfrac_mix[n] * s_ext[DENSITY];
-      }
-   }
-
-   s_ext[VELY] = 0.0;
-   if ( zone == 1 ) {
-	  amrex::Real vel = (prob_parm.mdot_ox / rhoOx) * mixingFunction +
-                        prob_parm.inertVel_ox /* (1.0 + (radius/0.025))*/ * (1.0 - mixingFunction);
-      s_ext[VELX] = vel;
-   } else if ( zone == 2 ) {
-	  amrex::Real vel = (prob_parm.mdot_fuel / rhoFuel) * mixingFunction +
-                        prob_parm.inertVel_fuel /* (1.0 + (radius/0.025))*/ * (1.0 - mixingFunction);
-      s_ext[VELX] = -vel;
-   }
+  const amrex::Real* prob_lo = geomdata.ProbLo();
+  const amrex::Real* prob_hi = geomdata.ProbHi();
+
+  AMREX_D_TERM(const amrex::Real Lx = prob_hi[0] - prob_lo[0];
+               , const amrex::Real Ly = prob_hi[1] - prob_lo[1];
+               , const amrex::Real Lz = prob_hi[2] - prob_lo[2]);
+
+  amrex::Real splitx = prob_lo[0] + 0.5 * Lx;
+  amrex::Real splity = prob_lo[1] + 0.5 * Ly;
+
+  amrex::Real massfrac_ox[NUM_SPECIES] = {0.0};
+  amrex::Real massfrac_fuel[NUM_SPECIES] = {0.0};
+  amrex::Real massfrac_inert[NUM_SPECIES] = {0.0};
+  amrex::Real massfrac_mix[NUM_SPECIES] = {0.0};
+
+  massfrac_ox[N2_ID] = prob_parm.Y_N2_ox;
+  massfrac_ox[O2_ID] = prob_parm.Y_O2_ox;
+  massfrac_fuel[NC12H26_ID] = prob_parm.Y_fuel; // Replace ID with gas fuel
+  massfrac_fuel[N2_ID] = prob_parm.Y_N2_fuel;
+  massfrac_inert[N2_ID] = 1.0;
+
+  auto eos = pele::physics::PhysicsType::eos();
+
+  for (int n = 0; n < NVAR; n++) {
+    s_ext[n] = 0.0;
+  }
+
+  // Get zone
+  int zone = -1;
+  if (x[0] < splitx) {
+    zone = 1; // Oxidizer side jet
+  } else {
+    zone = 2; // Fuel side jet
+  }
+  amrex::Real radius = std::abs(x[1] - splity);
+  amrex::Real mixingWidth = 0.1 * prob_parm.jetRadius;
+  amrex::Real mixingFunction =
+    0.5 * (1.0 + std::tanh((prob_parm.jetRadius - radius) / mixingWidth));
+
+  amrex::Real P_cgs = prob_parm.P_mean * 10.0;
+  amrex::Real rho_cgs, RhoH_temp;
+  amrex::Real rhoOx = -1.0;
+  {
+    eos.PYT2R(P_cgs, massfrac_ox, prob_parm.T_ox, rhoOx);
+    rhoOx *= 1.0e3;
+  }
+  amrex::Real rhoFuel = -1.0;
+  {
+    eos.PYT2R(P_cgs, massfrac_fuel, prob_parm.T_fuel, rhoFuel);
+    rhoFuel *= 1.0e3;
+  }
+  if (zone == 1) {
+    s_ext[TEMP] = prob_parm.T_ox * mixingFunction +
+                  prob_parm.T_inert * (1.0 - mixingFunction);
+    for (int n = 0; n < NUM_SPECIES; n++) {
+      massfrac_mix[n] = massfrac_ox[n] * mixingFunction +
+                        massfrac_inert[n] * (1.0 - mixingFunction);
+    }
+    eos.PYT2R(P_cgs, massfrac_mix, s_ext[TEMP], rho_cgs);
+    s_ext[DENSITY] = rho_cgs * 1.0e3;
+
+    eos.TY2H(s_ext[TEMP], massfrac_mix, RhoH_temp);
+    s_ext[RHOH] = RhoH_temp * 1.0e-4 * s_ext[DENSITY];
+
+    for (int n = 0; n < NUM_SPECIES; n++) {
+      s_ext[FIRSTSPEC + n] = massfrac_mix[n] * s_ext[DENSITY];
+    }
+  } else if (zone == 2) {
+    s_ext[TEMP] = prob_parm.T_fuel * mixingFunction +
+                  prob_parm.T_inert * (1.0 - mixingFunction);
+    for (int n = 0; n < NUM_SPECIES; n++) {
+      massfrac_mix[n] = massfrac_fuel[n] * mixingFunction +
+                        massfrac_inert[n] * (1.0 - mixingFunction);
+    }
+    eos.PYT2R(P_cgs, massfrac_mix, s_ext[TEMP], rho_cgs);
+    s_ext[DENSITY] = rho_cgs * 1.0e3;
+
+    eos.TY2H(s_ext[TEMP], massfrac_mix, RhoH_temp);
+    s_ext[RHOH] = RhoH_temp * 1.0e-4 * s_ext[DENSITY];
+
+    for (int n = 0; n < NUM_SPECIES; n++) {
+      s_ext[FIRSTSPEC + n] = massfrac_mix[n] * s_ext[DENSITY];
+    }
+  }
+
+  s_ext[VELY] = 0.0;
+  if (zone == 1) {
+    amrex::Real vel = (prob_parm.mdot_ox / rhoOx) * mixingFunction +
+                      prob_parm.inertVel_ox /* (1.0 + (radius/0.025))*/ *
+                        (1.0 - mixingFunction);
+    s_ext[VELX] = vel;
+  } else if (zone == 2) {
+    amrex::Real vel = (prob_parm.mdot_fuel / rhoFuel) * mixingFunction +
+                      prob_parm.inertVel_fuel /* (1.0 + (radius/0.025))*/ *
+                        (1.0 - mixingFunction);
+    s_ext[VELX] = -vel;
+  }
 }
 
 AMREX_GPU_DEVICE
 AMREX_FORCE_INLINE
 void
-zero_visc (int i, int j, int k,
-           amrex::Array4<amrex::Real> const& beta,
-           amrex::GeometryData const& geomdata,
-           amrex::Box const& domainBox,
-           const int  dir,
-           const int  beta_comp,
-           const int  nComp)
+zero_visc(
+  int i,
+  int j,
+  int k,
+  amrex::Array4<amrex::Real> const& beta,
+  amrex::GeometryData const& geomdata,
+  amrex::Box const& domainBox,
+  const int dir,
+  const int beta_comp,
+  const int nComp)
 {
-     amrex::ignore_unused(i,j,k,beta,geomdata,domainBox,dir,beta_comp,nComp);
-    // We treat species when beta_comp == 0 and nComp == NUM_SPECIES
-    // otherwise this routine could be called for other face diffusivity (Temp, velocity, ...)
+  amrex::ignore_unused(
+    i, j, k, beta, geomdata, domainBox, dir, beta_comp, nComp);
+  // We treat species when beta_comp == 0 and nComp == NUM_SPECIES
+  // otherwise this routine could be called for other face diffusivity (Temp,
+  // velocity, ...)
 }
 #endif
diff --git a/Exec/Production/CounterFlow/pelelmex_prob.cpp b/Exec/Production/CounterFlow/pelelmex_prob.cpp
index 94a3c6bf..538a5128 100644
--- a/Exec/Production/CounterFlow/pelelmex_prob.cpp
+++ b/Exec/Production/CounterFlow/pelelmex_prob.cpp
@@ -1,31 +1,32 @@
 #include <PeleLMeX.H>
 #include <AMReX_ParmParse.H>
 
-void PeleLM::readProbParm()
+void
+PeleLM::readProbParm()
 {
-   amrex::ParmParse pp("prob");
+  amrex::ParmParse pp("prob");
 
-   // CF params
-   pp.query("P_mean", PeleLM::prob_parm->P_mean);
-   pp.query("T_oxidizer", PeleLM::prob_parm->T_ox);
-   pp.query("T_fuel", PeleLM::prob_parm->T_fuel);
-   pp.query("T_inert", PeleLM::prob_parm->T_inert);
-   pp.query("Y_O2_ox", PeleLM::prob_parm->Y_O2_ox);
-   pp.query("Y_N2_ox", PeleLM::prob_parm->Y_N2_ox);
-   pp.query("Y_N2_fuel", PeleLM::prob_parm->Y_N2_fuel);
-   pp.query("Y_fuel", PeleLM::prob_parm->Y_fuel);
-   pp.query("Y_H2O_ign", PeleLM::prob_parm->Y_H2O_ign);
-   pp.query("Y_CO2_ign", PeleLM::prob_parm->Y_CO2_ign);
-   pp.query("massflow_ox", PeleLM::prob_parm->mdot_ox);
-   pp.query("massflow_fuel", PeleLM::prob_parm->mdot_fuel);
-   pp.query("pertmag",  PeleLM::prob_parm->pertmag);
-   pp.query("jet_radius",  PeleLM::prob_parm->jetRadius);
-   pp.query("inert_radius",  PeleLM::prob_parm->inertRadius);
-   pp.query("inert_velocity_ox",  PeleLM::prob_parm->inertVel_ox);
-   pp.query("inert_velocity_fuel",  PeleLM::prob_parm->inertVel_fuel);
+  // CF params
+  pp.query("P_mean", PeleLM::prob_parm->P_mean);
+  pp.query("T_oxidizer", PeleLM::prob_parm->T_ox);
+  pp.query("T_fuel", PeleLM::prob_parm->T_fuel);
+  pp.query("T_inert", PeleLM::prob_parm->T_inert);
+  pp.query("Y_O2_ox", PeleLM::prob_parm->Y_O2_ox);
+  pp.query("Y_N2_ox", PeleLM::prob_parm->Y_N2_ox);
+  pp.query("Y_N2_fuel", PeleLM::prob_parm->Y_N2_fuel);
+  pp.query("Y_fuel", PeleLM::prob_parm->Y_fuel);
+  pp.query("Y_H2O_ign", PeleLM::prob_parm->Y_H2O_ign);
+  pp.query("Y_CO2_ign", PeleLM::prob_parm->Y_CO2_ign);
+  pp.query("massflow_ox", PeleLM::prob_parm->mdot_ox);
+  pp.query("massflow_fuel", PeleLM::prob_parm->mdot_fuel);
+  pp.query("pertmag", PeleLM::prob_parm->pertmag);
+  pp.query("jet_radius", PeleLM::prob_parm->jetRadius);
+  pp.query("inert_radius", PeleLM::prob_parm->inertRadius);
+  pp.query("inert_velocity_ox", PeleLM::prob_parm->inertVel_ox);
+  pp.query("inert_velocity_fuel", PeleLM::prob_parm->inertVel_fuel);
 
-   // ignition params
-   pp.query("do_ignition",  PeleLM::prob_parm->do_ignit);
-   pp.query("ignition_SphRad",  PeleLM::prob_parm->ignitSphereRad);
-   pp.query("ignition_SphT",  PeleLM::prob_parm->ignitT);
+  // ignition params
+  pp.query("do_ignition", PeleLM::prob_parm->do_ignit);
+  pp.query("ignition_SphRad", PeleLM::prob_parm->ignitSphereRad);
+  pp.query("ignition_SphT", PeleLM::prob_parm->ignitT);
 }
diff --git a/Exec/Production/CounterFlow/pelelmex_prob_parm.H b/Exec/Production/CounterFlow/pelelmex_prob_parm.H
index efde759a..feb22e03 100644
--- a/Exec/Production/CounterFlow/pelelmex_prob_parm.H
+++ b/Exec/Production/CounterFlow/pelelmex_prob_parm.H
@@ -7,31 +7,31 @@ using namespace amrex::literals;
 
 struct ProbParm
 {
-   // Gas params
-   amrex::Real P_mean   = 101325.0_rt;
-   //amrex::Real T_ox     = 298.0_rt;
-   //amrex::Real T_fuel   = 298.0_rt;
-   //amrex::Real T_inert   = 298.0_rt;
-   amrex::Real T_ox     = 450.0_rt;
-   amrex::Real T_fuel   = 450.0_rt;
-   amrex::Real T_inert   = 450.0_rt;
-   amrex::Real Y_O2_ox   = 0.233_rt;
-   amrex::Real Y_N2_ox   = 0.767_rt;
-   amrex::Real Y_N2_fuel   = 0.0_rt;
-   amrex::Real Y_fuel   = 0.0_rt;
-   amrex::Real Y_H2O_ign   = 0.0_rt;
-   amrex::Real Y_CO2_ign   = 0.0_rt;
-   amrex::Real mdot_ox     = 0.01_rt;
-   amrex::Real mdot_fuel     = 0.01_rt;
-   amrex::Real pertmag  = 0.0004_rt;
-   amrex::Real inertVel_ox = 0.01_rt;
-   amrex::Real inertVel_fuel = 0.01_rt;
-   amrex::Real jetRadius = 0.004_rt;
-   amrex::Real inertRadius = 0.005_rt;
+  // Gas params
+  amrex::Real P_mean = 101325.0_rt;
+  // amrex::Real T_ox     = 298.0_rt;
+  // amrex::Real T_fuel   = 298.0_rt;
+  // amrex::Real T_inert   = 298.0_rt;
+  amrex::Real T_ox = 450.0_rt;
+  amrex::Real T_fuel = 450.0_rt;
+  amrex::Real T_inert = 450.0_rt;
+  amrex::Real Y_O2_ox = 0.233_rt;
+  amrex::Real Y_N2_ox = 0.767_rt;
+  amrex::Real Y_N2_fuel = 0.0_rt;
+  amrex::Real Y_fuel = 0.0_rt;
+  amrex::Real Y_H2O_ign = 0.0_rt;
+  amrex::Real Y_CO2_ign = 0.0_rt;
+  amrex::Real mdot_ox = 0.01_rt;
+  amrex::Real mdot_fuel = 0.01_rt;
+  amrex::Real pertmag = 0.0004_rt;
+  amrex::Real inertVel_ox = 0.01_rt;
+  amrex::Real inertVel_fuel = 0.01_rt;
+  amrex::Real jetRadius = 0.004_rt;
+  amrex::Real inertRadius = 0.005_rt;
 
-   // Misc
-   int do_ignit = 0;
-   amrex::Real ignitSphereRad = 0.003_rt;
-   amrex::Real ignitT  = 1800.0_rt;
+  // Misc
+  int do_ignit = 0;
+  amrex::Real ignitSphereRad = 0.003_rt;
+  amrex::Real ignitT = 1800.0_rt;
 };
 #endif