diff --git a/integration/nse_update_simplified_sdc.H b/integration/nse_update_simplified_sdc.H
index 19cdd82ef6..1293fbe73b 100644
--- a/integration/nse_update_simplified_sdc.H
+++ b/integration/nse_update_simplified_sdc.H
@@ -26,6 +26,106 @@
 using namespace amrex;
 
 #if defined(NSE_TABLE)
+
+///
+/// this acts as an explicit Euler step for the system (rho e, rho aux)
+/// on input, *_source are the reactive sources at time t0 and on output
+/// they are the sources at time t0+dt
+///
+template <bool skip_X_fill>
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+void evolve(const Real rho0, const Real T0, const Real rhoe0, const Real *rhoaux0,
+            const Real dt, const Real *ydot_a,
+            Real& rhoe_new, Real* rhoaux_new, Real* rhoX_new,
+            Real& rhoe_source, Real* rhoaux_source, const Real T_fixed) {
+
+    // integrate rho to the new time
+
+    Real rho_new = rho0 + dt * ydot_a[SRHO];
+
+    // integrate rhoe and rhoaux using the input sources
+
+    rhoe_new = rhoe0 + dt * ydot_a[SEINT] + dt * rhoe_source;
+    for (int n = 0; n < NumAux; n++) {
+        rhoaux_new[n] = rhoaux0[n] + dt * ydot_a[SFX+n] + dt * rhoaux_source[n];
+    }
+
+    // find the temperature at t0 + dt
+
+    eos_t eos_state;
+    eos_state.rho = rho_new;
+    eos_state.e = rhoe_new / rho_new;
+    eos_state.T = T0;
+    for (int n = 0; n < NumAux; n++) {
+        eos_state.aux[n] = rhoaux_new[n] / rho_new;
+    }
+
+    Real T_new;
+    if (T_fixed > 0) {
+        T_new = T_fixed;
+    } else {
+        eos(eos_input_re, eos_state);
+        T_new = eos_state.T;
+    }
+
+    // compute the plasma neutrino losses at t0 + dt
+
+    Real snu{0.0};
+    Real dsnudt{0.0};
+    Real dsnudd{0.0};
+    Real dsnuda{0.0};
+    Real dsnudz{0.0};
+
+    Real abar = eos_state.aux[iabar];
+    Real zbar = abar * eos_state.aux[iye];
+
+    constexpr int do_derivatives = 0;
+    sneut5<do_derivatives>(T_new, rho_new, abar, zbar,
+                           snu, dsnudt, dsnudd, dsnuda, dsnudz);
+
+    // call the NSE table
+
+    Real abar_out;
+    Real bea_out;
+    Real dyedt;
+    Real dabardt;
+    Real dbeadt;
+    Real enu;
+
+    nse_interp(T_new, rho_new, eos_state.aux[iye],
+               abar_out, bea_out, dyedt, dabardt, dbeadt, enu, rhoX_new,
+               skip_X_fill);
+
+    if constexpr (!skip_X_fill) {
+        for (int n = 0; n < NumSpec; ++n) {
+            rhoX_new[n] *= rho_new;
+        }
+    }
+
+    // construct the sources at time t0+dt
+
+    // note that abar, (B/A) here have already seen advection
+    // implicitly
+
+    Real rho_bea_tilde = rho_new * bea_out - dt * ydot_a[SFX+ibea];
+    Real rho_dBEA = rho_bea_tilde - rhoaux0[ibea]; // this is MeV / nucleon * g / cm**3
+
+    rhoe_source = rho_dBEA * C::MeV2eV * C::ev2erg * C::n_A / dt;
+    rhoe_source -= rho_new * (enu + snu);
+
+    rhoaux_source[iabar] = 0.0;
+    rhoaux_source[iye] = rho_new * dyedt;
+    rhoaux_source[ibea] = rho_dBEA / dt;
+
+    // reset the new <B/A> -- we want to ensure that if we use the
+    // updated T(rho, e, Y_e, Abar) that we get exactly the same B/A,
+    // so there is no instantaneous change in the binding energy
+
+    rhoaux_new[ibea] = rho_new * bea_out;
+    rhoaux_new[iabar] = rho_new * abar_out;
+
+}
+
 ///
 /// update the state due to NSE changes for the simplified-SDC case
 /// this version works with the tabulated NSE and requires AUX_THERMO
@@ -40,141 +140,228 @@ void sdc_nse_burn(BurnT& state, const Real dt) {
     state.n_rhs = 0;
     state.n_jac = 0;
 
-    // call the NSE table to get (dYe/dt)^n
-    Real abar_out;
-    Real dq_out;
-    Real dyedt;
-    Real dabardt;
-    Real dbeadt;
-    Real enu;
-    Real X[NumSpec];
+    // we need the initial Ye to get the NSE state.  We also
+    // want the initial rho <B/A> since we may not be entering
+    // this routine already in NSE.  This means that there will
+    // be an energy release from us instantaneously adjusting
+    // into NSE (the first call to nse_interp) + an energy
+    // release from the evolution of NSE over the timestep
+
+    // we will write the energy update as:
+    //
+    // e^new = e^old + dbea_instantantous + dbea_over_dt
+    //
+    // where dbea_over_dt will be an integral over dt using
+    // a simple RK2 method
 
     Real ye_in = state.y[SFX+iye] / state.rho;
+    Real rho_bea_orig = state.y[SFX+ibea];
+    Real rhoe_orig = state.y[SEINT];
+
+    // density and momentum have no reactive sources
+    Real rho_old = state.y[SRHO];
+
+    state.y[SRHO] += dt * state.ydot_a[SRHO];
+    state.y[SMX] += dt * state.ydot_a[SMX];
+    state.y[SMY] += dt * state.ydot_a[SMY];
+    state.y[SMZ] += dt * state.ydot_a[SMZ];
 
     // if we are doing drive_initial_convection, we want to use
     // the temperature that comes in through T_fixed
 
     Real T_in = state.T_fixed > 0.0_rt ? state.T_fixed : state.T;
 
-    // get the current NSE state from the table
+    // get the instantaneous change in B/A, as well as the source
+    // terms for the new evolution.
 
-    nse_interp(T_in, state.rho, ye_in,
-               abar_out, dq_out, dyedt, dabardt, dbeadt, enu, X);
+    Real abar_out{};
+    Real bea_out{};
+    Real dyedt{};
+    Real dabardt{};
+    Real dbeadt{};
+    Real enu{};
+    Real X[NumSpec] = {};
 
-    Real dyedt_old = dyedt;
+    {
+        constexpr bool skip_X_fill{true};
 
-    // density and momentum have no reactive sources
-    Real rho_old = state.y[SRHO];
-    Real rho_bea_old = state.y[SFX+ibea];
+        nse_interp(T_in, state.rho, ye_in,
+                   abar_out, bea_out, dyedt, dabardt, dbeadt, enu, X, skip_X_fill);
+    }
 
-    state.y[SRHO] += dt * state.ydot_a[SRHO];
-    state.y[SMX] += dt * state.ydot_a[SMX];
-    state.y[SMY] += dt * state.ydot_a[SMY];
-    state.y[SMZ] += dt * state.ydot_a[SMZ];
+    // the starting state for integration uses the state that we just
+    // adjusted instantaneously into NSE
 
-    // predict the U^{n+1,*} state with only estimates of the aux
-    // reaction terms and advection to dt
+    Real rhoaux_old[NumAux];
+    rhoaux_old[iye] = rho_old * ye_in;
+    rhoaux_old[iabar] = rho_old * abar_out;
+    rhoaux_old[ibea] = rho_old * bea_out;
 
-    eos_re_t eos_state;
-    eos_state.T = T_in; // initial guess
+    Real rhoenuc = (rho_old * bea_out - rho_bea_orig) * C::MeV2eV * C::ev2erg * C::n_A;
+    Real rhoe_old = state.y[SEINT] + rhoenuc;
 
-    // initial aux_sources
-    Real aux_source[NumAux] = {0.0_rt};
-    aux_source[iye] = rho_old * dyedt;
 
-    Real rho_aux_new[NumAux];
+    // get the starting temperature
+
+    eos_t eos_state;
+    eos_state.rho = rho_old;
+    eos_state.e = rhoe_old / rho_old;
+    for (int n = 0; n < NumAux; ++n) {
+        eos_state.aux[n] = rhoaux_old[n] / rho_old;
+    }
+    eos_state.T = T_in;
+
+    Real T_old;
+    if (state.T_fixed > 0) {
+        T_old = state.T_fixed;
+    } else {
+        eos(eos_input_re, eos_state);
+        T_old = eos_state.T;
+    }
+
+    // start of integration
+
+    // we want to integrate from t^n to t^{n+1/2}
+    // start by computing the sources
+
+    // the aux sources are just due to weak rates
+    Real rhoaux_source[NumAux];
+    rhoaux_source[iye] = rho_old * dyedt;
+    rhoaux_source[ibea] = rho_old * dbeadt;
+    rhoaux_source[iabar] = 0.0;
+
+    // the energy source is due to plasma neutrinos + weak rate neutrinos
+
+    Real snu{};
+    Real dsnudt{};
+    Real dsnudd{};
+    Real dsnuda{};
+    Real dsnudz{};
+    constexpr int do_derivatives = 0;
+
+    Real zbar = abar_out * ye_in;
+    sneut5<do_derivatives>(T_old, rho_old, abar_out, zbar,
+                           snu, dsnudt, dsnudd, dsnuda, dsnudz);
+
+    Real rhoe_source = -rho_old * (snu + enu);
 
     Real rhoe_new;
-    Real rho_enucdot = 0.0_rt;
+    Real rhoaux_new[NumAux];
+    Real rhoX_new[NumSpec];
 
-    Real rho_half = 0.5_rt * (rho_old + state.y[SRHO]);
+    // evolve for eps dt -- the sources will be overwritten with the
+    // sources at the end point of this integration
 
-    // predict the new aux for the first iteration -- this is really
-    // just including the advection bits
+    {
+        constexpr bool skip_X_fill{true};
 
-    for (int n = 0; n < NumAux; n++) {
-        rho_aux_new[n] = state.y[SFX+n] + dt * state.ydot_a[SFX+n] + dt * aux_source[n];
+        evolve<skip_X_fill>(rho_old, T_old, rhoe_old, rhoaux_old,
+                            0.05_rt*dt, state.ydot_a,
+                            rhoe_new, rhoaux_new, rhoX_new,
+                            rhoe_source, rhoaux_source, state.T_fixed);
     }
 
-    for (int iter = 0; iter < integrator_rp::nse_iters; iter++) {
+    //std::cout << "starting integration, dt = " << dt << std::endl;
+    //std::cout << " rhoe: " << rhoe_old << " " << rhoe_source << std::endl;
+    //std::cout << " rhoYe: " << rhoaux_old[iye] << " " << rhoaux_source[iye] << std::endl;
+    //std::cout << " rhoAbar: " << rhoaux_old[iabar] << " " << rhoaux_source[iabar] << std::endl;
+    //std::cout << " rhoBA: " << rhoaux_old[ibea] << " " << rhoaux_source[ibea] << std::endl;
 
-        // update (rho e)^{n+1} based on the new energy generation rate
-        rhoe_new = state.y[SEINT] + dt * state.ydot_a[SEINT] + dt * rho_enucdot;
+    // evolve for dt/2 -- the sources will be overwritten with the
+    // sources at the end point of this integration
 
-        // call the EOS to get the updated T*
+    {
+        constexpr bool skip_X_fill{true};
 
-        Real T_new;
-        eos_state.rho = state.y[SRHO];
-        eos_state.e = rhoe_new / state.y[SRHO];
-        for (int n = 0; n < NumAux; n++) {
-            eos_state.aux[n] = rho_aux_new[n] / state.y[SRHO];
-        }
+        evolve<skip_X_fill>(rho_old, T_old, rhoe_old, rhoaux_old,
+                            0.5_rt*dt, state.ydot_a,
+                            rhoe_new, rhoaux_new, rhoX_new,
+                            rhoe_source, rhoaux_source, state.T_fixed);
+    }
 
-        if (state.T_fixed > 0) {
-            T_new = state.T_fixed;
-        } else {
-            eos(eos_input_re, eos_state);
-            T_new = eos_state.T;
-        }
+    // evolve from dt/2 to dt/2 + eps -- the sources will be overwritten with the
+    // sources at the end point of this integration
 
-        // call the NSE table using the * state to get the t^{n+1}
-        // source estimates
+    {
+        constexpr bool skip_X_fill{true};
 
-        nse_interp(T_new, eos_state.rho, eos_state.aux[iye],
-                   abar_out, dq_out, dyedt, dabardt, dbeadt, enu, X);
+        Real rho_tmp = rho_old + 0.5*dt * state.ydot_a[SRHO];
+        Real rhoe_tmp{rhoe_new};
+        Real rhoaux_tmp[NumAux] = {rhoaux_new[0], rhoaux_new[1], rhoaux_new[2]};
 
-        // note that the abar, dq (B/A), and X here have all already
-        // seen advection implicitly
+        evolve<skip_X_fill>(rho_tmp, T_old, rhoe_tmp, rhoaux_tmp,
+                            0.05_rt*dt, state.ydot_a,
+                            rhoe_new, rhoaux_new, rhoX_new,
+                            rhoe_source, rhoaux_source, state.T_fixed);
+    }
 
+    // this gives us the sources at n+1/2.  Now we do the full dt step,
+    // starting again from the beginning, but using the time-centered
+    // sources
 
-        // compute the energy release -- we need to remove the advective part
 
-        Real rho_bea_tilde = state.y[SRHO] * dq_out - dt * state.ydot_a[SFX+ibea];
-        Real rho_dBEA = rho_bea_tilde - rho_bea_old;  // this is MeV / nucleon * g / cm**3
+    //std::cout << "second half of integration, dt = " << dt << std::endl;
+    //std::cout << " rhoe: " << rhoe_old << " " << rhoe_source << std::endl;
+    //std::cout << " rhoYe: " << rhoaux_old[iye] << " " << rhoaux_source[iye] << std::endl;
+    //std::cout << " rhoAbar: " << rhoaux_old[iabar] << " " << rhoaux_source[iabar] << std::endl;
+    //std::cout << " rhoBA: " << rhoaux_old[ibea] << " " << rhoaux_source[ibea] << std::endl;
 
-        // convert the energy to erg / cm**3
-        rho_enucdot  = rho_dBEA * C::MeV2eV * C::ev2erg * C::n_A / dt;
 
-        // update the new state for the next pass
+    {
+        constexpr bool skip_X_fill{false};
 
-        aux_source[iye] = 0.5_rt * rho_half * (dyedt_old + dyedt);
-        rho_aux_new[iye] = state.y[SFX+iye] + dt * state.ydot_a[SFX+iye] + dt * aux_source[iye];
+        evolve<skip_X_fill>(rho_old, T_old, rhoe_old, rhoaux_old,
+                            dt, state.ydot_a,
+                            rhoe_new, rhoaux_new, rhoX_new,
+                            rhoe_source, rhoaux_source, state.T_fixed);
+    }
 
-        rho_aux_new[iabar] = state.y[SRHO] * abar_out;
-        rho_aux_new[ibea] = state.y[SRHO] * dq_out;
+    // done with the integration.  This gave us a second-order
+    // accurate update to the rhoe and Ye, and the aux data and mass
+    // fractions are consistent with the NSE table at the end of the
+    // integration
 
-    }
+    //std::cout << "done" << std::endl;
+    //std::cout << " rhoe: " << rhoe_new << std::endl;
+    //std::cout << " rhoYe: " << rhoaux_new[iye] << std::endl;
+    //std::cout << " rhoAbar: " << rhoaux_new[iabar] << std::endl;
+    //std::cout << " rhoBA: " << rhoaux_new[ibea] << std::endl;
+
+    //amrex::Error("stop");
 
     // now update the aux quantities
 
     // the new mass fractions are just those that come from the table
     // make sure they are normalized
-    Real sum_X{0.0_rt};
-    for (auto & xn : X) {
-        xn = amrex::max(small_x, amrex::min(1.0_rt, xn));
-        sum_X += xn;
+    Real sum_rhoX{0.0_rt};
+    for (auto & rx : rhoX_new) {
+        rx = amrex::max(state.y[SRHO] * small_x, amrex::min(state.y[SRHO], rx));
+        sum_rhoX += rx;
     }
 
-    for (auto & xn : X) {
-        xn /= sum_X;
+    for (auto & rx : rhoX_new) {
+        rx /= sum_rhoX;
     }
 
     for (int n = 0; n < NumSpec; ++n) {
-        state.y[SFS+n] = state.y[SRHO] * X[n];
+        state.y[SFS+n] = rhoX_new[n];
     }
 
-    // aux data comes from the table (except Ye, which we predicted)
+    // aux data comes from the integration
 
-    state.y[SFX+iye] = eos_state.rho * eos_state.aux[iye];
-    state.y[SFX+iabar] = eos_state.rho * abar_out;
-    state.y[SFX+ibea] = eos_state.rho * dq_out;
+    for (int n = 0; n < NumAux; ++n) {
+        state.y[SFX+n] = rhoaux_new[n];
+    }
 
 
     // density and momenta have already been updated
 
-    // update the total and internal energy now
+    // get the energy release from the change in rhoe over the timestep
+    // excluding any advection, and use that to update the total energy
 
-    state.y[SEINT] += dt * state.ydot_a[SEINT] + dt * rho_enucdot;
+    Real rho_enucdot = (rhoe_new - rhoe_orig) / dt - state.ydot_a[SEINT];
+
+    state.y[SEINT] = rhoe_new;
     state.y[SEDEN] += dt * state.ydot_a[SEDEN] + dt * rho_enucdot;
 
 }
@@ -192,53 +379,62 @@ void sdc_nse_burn(BurnT& state, const Real dt) {
     state.n_rhs = 0;
     state.n_jac = 0;
 
-    // call the NSE table to get (dYe/dt)^n
-    Real abar_out;
-    Real dq_out;
-    Real dyedt;
-    Real X[NumSpec];
+    // store the initial mass fractions -- we will need these
+    // to compute the energy release.
+
+    Real X_old[NumSpec];
+
+    for (int n = 0; n < NumSpec; ++n) {
+        X_old[n] = state.y[SFS+n] / state.y[SRHO];
+    }
+
+    // density and momentum have no reactive sources
+    Real rho_old = state.y[SRHO];
 
-    // TODO: are state.y[SFS:] and state.xn[:] synced?
+    state.y[SRHO] += dt * state.ydot_a[SRHO];
+    state.y[SMX] += dt * state.ydot_a[SMX];
+    state.y[SMY] += dt * state.ydot_a[SMY];
+    state.y[SMZ] += dt * state.ydot_a[SMZ];
 
     // if we are doing drive_initial_convection, we want to use
     // the temperature that comes in through T_fixed
 
     Real T_in = state.T_fixed > 0.0_rt ? state.T_fixed : state.T;
 
-    // We will get initial NSE prediction using the input X's
-    BurnT nse_state_in{state};
-    nse_state_in.T = T_in;
+    // get the neutrino loss term -- we want to use the state that we
+    // came in here with, so the original Abar and Zbar
+    Real snu{0.0};
+    Real dsnudt{0.0};
+    Real dsnudd{0.0};
+    Real dsnuda{0.0};
+    Real dsnudz{0.0};
 
-    // solve for the NSE state directly -- we'll have it compute
-    // ye from the input Xs
-    auto nse_state = get_actual_nse_state(nse_state_in);
-
-    // compute the new binding energy (B/A) and dyedt
-    dq_out = 0.0;
+    Real abar{0.0};
+    Real zbar{0.0};
     for (int n = 0; n < NumSpec; ++n) {
-        dq_out += nse_state.xn[n] * network::bion(n+1) * aion_inv[n];
+        abar += X_old[n] * aion_inv[n];
+        zbar += X_old[n] * zion[n] * aion_inv[n];
     }
+    abar = 1.0 / abar;
+    zbar *= abar;
 
-    dyedt = 0.0;   // we can update this in the future by calling actual_rhs()
+    constexpr int do_derivatives = 0;
+    sneut5<do_derivatives>(T_in, rho_old, abar, zbar,
+                           snu, dsnudt, dsnudd, dsnuda, dsnudz);
 
-    Real dyedt_old = dyedt;
+    Real snu_old = snu;
 
-    // density and momentum have no reactive sources
-    Real rho_old = state.y[SRHO];
 
-    // compute the original rho (B/A) of the composition
-    Real rho_bea_old = 0.0;
-    for (int n = 0; n < NumSpec; ++n) {
-        rho_bea_old += state.y[SFS+n] * network::bion(n+1) * aion_inv[n];
-    }
+    // if our network could return the evolution of Ye due to the
+    // weak interactions, we would evaluate the NSE state here and
+    // get dYe/dt.
+
+    Real dyedt_old = 0.0;
 
-    state.y[SRHO] += dt * state.ydot_a[SRHO];
-    state.y[SMX] += dt * state.ydot_a[SMX];
-    state.y[SMY] += dt * state.ydot_a[SMY];
-    state.y[SMZ] += dt * state.ydot_a[SMZ];
 
     // predict the U^{n+1,*} state with only estimates of the X
-    // updates with advection to dt
+    // updates with advection to dt and the neutrino loss term in
+    // energy
 
     BurnT burn_state;
     burn_state.T = T_in; // initial guess
@@ -250,7 +446,7 @@ void sdc_nse_burn(BurnT& state, const Real dt) {
     Real rhoX_new[NumSpec];
 
     Real rhoe_new;
-    Real rho_enucdot = 0.0_rt;
+    Real rho_enucdot = -rho_old * snu;
 
     Real rho_half = 0.5_rt * (rho_old + state.y[SRHO]);
 
@@ -261,6 +457,8 @@ void sdc_nse_burn(BurnT& state, const Real dt) {
         rhoX_new[n] = state.y[SFS+n] + dt * state.ydot_a[SFS+n] + dt * rhoX_source[n];
     }
 
+    burn_t nse_state;
+
     for (int iter = 0; iter < integrator_rp::nse_iters; iter++) {
 
         // update (rho e)^{n+1} based on the new energy generation rate
@@ -290,27 +488,39 @@ void sdc_nse_burn(BurnT& state, const Real dt) {
 
         nse_state = get_actual_nse_state(burn_state);
 
-        // compute (B/A) and dyedt
-        dq_out = 0.0;
+        // compute the energy release.  The mass fractions in nse_state.xn[]
+        // include the advective parts, so first we need to remove that.
+
+        Real rhoX_tilde[NumSpec];
         for (int n = 0; n < NumSpec; ++n) {
-            dq_out += nse_state.xn[n] * network::bion(n+1) * aion_inv[n];
+            rhoX_tilde[n] = state.y[SRHO] * nse_state.xn[n] - dt * state.ydot_a[SFS+n];
         }
 
-        dyedt = 0.0_rt;   // we can update this in the future by calling actual_rhs()
+        Real dyedt = 0.0_rt;   // we can update this in the future by calling actual_rhs()
+
+        // we want to compute (rho eps) = - N_A c^2 sum{m_i (rhoX_tilde - rhoX_old) / A_i}
+        rho_enucdot = 0.0;
+        for (int n = 0; n < NumSpec; ++n) {
+            rho_enucdot += (rhoX_tilde[n] - rho_old * X_old[n]) *
+                network::mion(n+1) * aion_inv[n];
+        }
+        rho_enucdot *= C::Legacy::enuc_conv2;
 
-        // compute the energy release -- we need to remove the
-        // advective part.  To do this, we must first compute the
-        // advective update for B/A from those of the mass fractions
-        Real ydot_a_BA = 0.0_rt;
+        // now get the updated neutrino term
+        abar = 0.0;
+        zbar = 0.0;
         for (int n = 0; n < NumSpec; ++n) {
-            ydot_a_BA += network::bion(n+1) * aion_inv[n] * state.ydot_a[SFS+n];
+            abar += nse_state.xn[n] * aion_inv[n];
+            zbar += nse_state.xn[n] * zion[n] * aion_inv[n];
         }
+        abar = 1.0 / abar;
+        zbar *= abar;
 
-        Real rho_bea_tilde = state.y[SRHO] * dq_out - dt * ydot_a_BA;
-        Real rho_dBEA = rho_bea_tilde - rho_bea_old;  // this is MeV / nucleon * g / cm**3
+        constexpr int do_derivatives = 0;
+        sneut5<do_derivatives>(T_new, state.y[SRHO], abar, zbar,
+                               snu, dsnudt, dsnudd, dsnuda, dsnudz);
 
-        // convert the energy to erg / cm**3
-        rho_enucdot  = rho_dBEA * C::MeV2eV * C::ev2erg * C::n_A / dt;
+        rho_enucdot -= 0.5_rt * rho_half * (snu_old + snu);
 
         // update the new state for the next pass -- this should
         // already implicitly have the advective portion included,
@@ -327,6 +537,7 @@ void sdc_nse_burn(BurnT& state, const Real dt) {
     // the new mass fractions are just those that come from the table
     // make sure they are normalized
     Real sum_X{0.0_rt};
+    Real X[NumSpec] = {0.0_rt};
     for (int n = 0; n < NumSpec; ++n) {
         X[n] = amrex::max(small_x, amrex::min(1.0_rt, nse_state.xn[n]));
         sum_X += X[n];