diff --git a/Source/driver/Castro.H b/Source/driver/Castro.H
index 5d17812359..e9ddba5ceb 100644
--- a/Source/driver/Castro.H
+++ b/Source/driver/Castro.H
@@ -770,11 +770,11 @@ public:
 ///
     static void reset_internal_energy (const amrex::Box& bx,
 #ifdef MHD
-                                amrex::Array4<amrex::Real> const Bx,
-                                amrex::Array4<amrex::Real> const By,
-                                amrex::Array4<amrex::Real> const Bz,
+                                amrex::Array4<amrex::Real> const& Bx,
+                                amrex::Array4<amrex::Real> const& By,
+                                amrex::Array4<amrex::Real> const& Bz,
 #endif
-                                amrex::Array4<amrex::Real> const u);
+                                amrex::Array4<amrex::Real> const& u);
 
 
 #ifdef MHD
diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp
index f28d228b25..e524923df8 100644
--- a/Source/driver/Castro.cpp
+++ b/Source/driver/Castro.cpp
@@ -3739,9 +3739,9 @@ Castro::extern_init ()
 void
 Castro::reset_internal_energy(const Box& bx,
 #ifdef MHD
-                              Array4<Real> const Bx, Array4<Real> const By, Array4<Real> const Bz,
+                              Array4<Real> const& Bx, Array4<Real> const& By, Array4<Real> const& Bz,
 #endif
-                              Array4<Real> const u)
+                              Array4<Real> const& u)
 {
     BL_PROFILE("Castro::reset_internal_energy(Fab)");
 
@@ -3751,53 +3751,56 @@ Castro::reset_internal_energy(const Box& bx,
     amrex::ParallelFor(bx,
     [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
     {
-        Real rhoInv = 1.0_rt / u(i,j,k,URHO);
-        Real Up = u(i,j,k,UMX) * rhoInv;
-        Real Vp = u(i,j,k,UMY) * rhoInv;
-        Real Wp = u(i,j,k,UMZ) * rhoInv;
-        Real ke = 0.5_rt * (Up * Up + Vp * Vp + Wp * Wp);
+
+        auto U_cell = u.cellData(i, j, k);
+
+        const Real rhoInv = 1.0_rt / U_cell[URHO];
+        const Real Up = U_cell[UMX] * rhoInv;
+        const Real Vp = U_cell[UMY] * rhoInv;
+        const Real Wp = U_cell[UMZ] * rhoInv;
+        const Real ke = 0.5_rt * (Up * Up + Vp * Vp + Wp * Wp);
 
         eos_re_t eos_state;
 
-        eos_state.rho = u(i,j,k,URHO);
+        eos_state.rho = U_cell[URHO];
         eos_state.T   = lsmall_temp;
         for (int n = 0; n < NumSpec; ++n) {
-            eos_state.xn[n] = u(i,j,k,UFS+n) * rhoInv;
+            eos_state.xn[n] = U_cell[UFS+n] * rhoInv;
         }
 #if NAUX_NET > 0
         for (int n = 0; n < NumAux; ++n) {
-            eos_state.aux[n] = u(i,j,k,UFX+n) * rhoInv;
+            eos_state.aux[n] = U_cell[UFX+n] * rhoInv;
         }
 #endif
 
         eos(eos_input_rt, eos_state);
 
-        Real small_e = eos_state.e;
+        const Real small_e = eos_state.e;
 
 #ifdef MHD
-        Real bx_cell_c = 0.5_rt * (Bx(i,j,k) + Bx(i+1,j,k));
-        Real by_cell_c = 0.5_rt * (By(i,j,k) + By(i,j+1,k));
-        Real bz_cell_c = 0.5_rt * (Bz(i,j,k) + Bz(i,j,k+1));
+        const Real bx_cell_c = 0.5_rt * (Bx(i,j,k) + Bx(i+1,j,k));
+        const Real by_cell_c = 0.5_rt * (By(i,j,k) + By(i,j+1,k));
+        const Real bz_cell_c = 0.5_rt * (Bz(i,j,k) + Bz(i,j,k+1));
 
-        Real B_ener = 0.5_rt * (bx_cell_c*bx_cell_c +
-                                by_cell_c*by_cell_c +
-                                bz_cell_c*bz_cell_c);
+        const Real B_ener = 0.5_rt * (bx_cell_c*bx_cell_c +
+                                      by_cell_c*by_cell_c +
+                                      bz_cell_c*bz_cell_c);
 #else
-        Real B_ener = 0.0_rt;
+        const Real B_ener = 0.0_rt;
 #endif
 
         // Ensure the internal energy is at least as large as this minimum
         // from the EOS; the same holds true for the total energy.
 
-        u(i,j,k,UEINT) = amrex::max(u(i,j,k,UEINT), u(i,j,k,URHO) * small_e);
-        u(i,j,k,UEDEN) = amrex::max(u(i,j,k,UEDEN), u(i,j,k,URHO) * (small_e + ke) + B_ener);
+        U_cell[UEINT] = amrex::max(U_cell[UEINT], U_cell[URHO] * small_e);
+        U_cell[UEDEN] = amrex::max(U_cell[UEDEN], U_cell[URHO] * (small_e + ke) + B_ener);
 
         // Apply the dual energy criterion: get e from E if (E - K) > eta * E.
 
-        Real rho_eint = u(i,j,k,UEDEN) - u(i,j,k,URHO) * ke - B_ener;
+        const Real rho_eint = U_cell[UEDEN] - U_cell[URHO] * ke - B_ener;
 
-        if (rho_eint > ldual_energy_eta2 * u(i,j,k,UEDEN)) {
-            u(i,j,k,UEINT) = rho_eint;
+        if (rho_eint > ldual_energy_eta2 * U_cell[UEDEN]) {
+            U_cell[UEINT] = rho_eint;
         }
     });
 }
diff --git a/Source/reactions/sdc_cons_to_burn.H b/Source/reactions/sdc_cons_to_burn.H
index d4ef43fd8b..0ab2f404a1 100644
--- a/Source/reactions/sdc_cons_to_burn.H
+++ b/Source/reactions/sdc_cons_to_burn.H
@@ -9,23 +9,24 @@ copy_cons_to_burn_type(const int i, const int j, const int k,
                        Array4<const Real> const& state,
                        burn_t& burn_state) {
 
-    burn_state.y[SRHO] = state(i,j,k,URHO);
-    burn_state.y[SMX] = state(i,j,k,UMX);
-    burn_state.y[SMY] = state(i,j,k,UMY);
-    burn_state.y[SMZ] = state(i,j,k,UMZ);
-    burn_state.y[SEDEN] = state(i,j,k,UEDEN);
-    burn_state.y[SEINT] = state(i,j,k,UEINT);
-    for (int n = 0; n < NumSpec; n++) {
-        burn_state.y[SFS+n] = state(i,j,k,UFS+n);
-    }
+    const auto U_cell = state.cellData(i,j,k);
+    burn_state.y[SRHO] = U_cell[URHO];
+    burn_state.y[SMX] = U_cell[UMX];
+    burn_state.y[SMY] = U_cell[UMY];
+    burn_state.y[SMZ] = U_cell[UMZ];
+    burn_state.y[SEDEN] = U_cell[UEDEN];
+    burn_state.y[SEINT] = U_cell[UEINT];
+    amrex::constexpr_for<0, NumSpec>([&] (auto n) {
+        burn_state.y[SFS+n] = U_cell[UFS+n];
+    });
 #if NAUX_NET > 0
     for (int n = 0; n < NumAux; n++) {
         burn_state.y[SFX+n] = state(i,j,k,UFX+n);
     }
 #endif
 
-    burn_state.T = state(i,j,k,UTEMP);
-    burn_state.rho = state(i,j,k,URHO);
+    burn_state.T = U_cell[UTEMP];
+    burn_state.rho = burn_state.y[SRHO];
 
 }