diff --git a/external/parthenon b/external/parthenon
index eb16b502..01cd49ec 160000
--- a/external/parthenon
+++ b/external/parthenon
@@ -1 +1 @@
-Subproject commit eb16b502dc03d5288ef7ea06fb3399b480e17374
+Subproject commit 01cd49ec65506421a1a0b86dd9ae122de68a663d
diff --git a/src/fixup/fixup_radc2p.cpp b/src/fixup/fixup_radc2p.cpp
index 63d6c1da..3cfa6962 100644
--- a/src/fixup/fixup_radc2p.cpp
+++ b/src/fixup/fixup_radc2p.cpp
@@ -41,14 +41,11 @@ using robust::ratio;
 
 namespace fixup {
 
-template <typename T, class CLOSURE>
+template <typename T, class CLOSURE, radiation::RadEnergyMoment EMOMENT>
 TaskStatus RadConservedToPrimitiveFixupImpl(T *rc) {
   namespace p = fluid_prim;
   namespace c = fluid_cons;
   namespace impl = internal_variables;
-  namespace ir = radmoment_internal;
-  namespace pr = radmoment_prim;
-  namespace cr = radmoment_cons;
 
   auto *pmb = rc->GetParentPointer().get();
   IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior);
@@ -62,9 +59,10 @@ TaskStatus RadConservedToPrimitiveFixupImpl(T *rc) {
   bool enable_c2p_fixup = fix_pkg->Param<bool>("enable_c2p_fixup");
   bool update_rad = rad_pkg->Param<bool>("active");
   if (!enable_c2p_fixup || !update_rad) return TaskStatus::complete;
-
-  const std::vector<std::string> vars({p::velocity, p::ye, c::ye, pr::J, pr::H, cr::E,
-                                       cr::F, ir::tilPi, ir::c2pfail, impl::fail});
+  
+  auto var = radiation::RadiationVariableNames::GetByMomentType(EMOMENT);
+  const std::vector<std::string> vars({p::velocity, p::ye, c::ye, var.J, var.H, var.E,
+                                       var.F, var.tilPi, var.c2pfail, impl::fail});
 
   PackIndexMap imap;
   auto v = rc->PackVariables(vars, imap);
@@ -72,13 +70,13 @@ TaskStatus RadConservedToPrimitiveFixupImpl(T *rc) {
   auto idx_pvel = imap.GetFlatIdx(p::velocity);
   int pye = imap[p::ye].second; // negative if not present
   int cye = imap[c::ye].second;
-  auto idx_J = imap.GetFlatIdx(pr::J, false);
-  auto idx_H = imap.GetFlatIdx(pr::H, false);
-  auto idx_E = imap.GetFlatIdx(cr::E, false);
-  auto idx_F = imap.GetFlatIdx(cr::F, false);
+  auto idx_J = imap.GetFlatIdx(var.J, false);
+  auto idx_H = imap.GetFlatIdx(var.H, false);
+  auto idx_E = imap.GetFlatIdx(var.E, false);
+  auto idx_F = imap.GetFlatIdx(var.F, false);
   int ifluidfail = imap[impl::fail].first;
-  int iradfail = imap[ir::c2pfail].first;
-  auto iTilPi = imap.GetFlatIdx(ir::tilPi, false);
+  int iradfail = imap[var.c2pfail].first;
+  auto iTilPi = imap.GetFlatIdx(var.tilPi, false);
 
   bool report_c2p_fails = fix_pkg->Param<bool>("report_c2p_fails");
   if (report_c2p_fails) {
@@ -227,22 +225,28 @@ TaskStatus RadConservedToPrimitiveFixup(T *rc) {
   StateDescriptor *rad_pkg = pm->packages.Get("radiation").get();
   StateDescriptor *fix_pkg = pm->packages.Get("fixup").get();
   const bool enable_rad_floors = fix_pkg->Param<bool>("enable_rad_floors");
+  auto do_num = rad_pkg->Param<bool>("do_number_evolution");
   std::string method;
   if (enable_rad_floors) {
     method = rad_pkg->Param<std::string>("method");
   } else {
     return TaskStatus::complete;
   }
-
+  using namespace radiation;
   // TODO(BRR) share these settings somewhere else. Set at configure time?
-  using settings =
+  using settings_N =
+      ClosureSettings<ClosureEquation::number_conserve, ClosureVerbosity::quiet>;
+  using settings_E =
       ClosureSettings<ClosureEquation::energy_conserve, ClosureVerbosity::quiet>;
   if (method == "moment_m1") {
-    return RadConservedToPrimitiveFixupImpl<T, radiation::ClosureM1<settings>>(rc);
+    if (do_num) RadConservedToPrimitiveFixupImpl<T, ClosureM1<settings_N>, RadEnergyMoment::Number>(rc);
+    return RadConservedToPrimitiveFixupImpl<T, ClosureM1<settings_E>, RadEnergyMoment::Energy>(rc);
   } else if (method == "moment_eddington") {
-    return RadConservedToPrimitiveFixupImpl<T, radiation::ClosureEdd<settings>>(rc);
+    if (do_num) RadConservedToPrimitiveFixupImpl<T, ClosureEdd<settings_N>, RadEnergyMoment::Number>(rc);
+    return RadConservedToPrimitiveFixupImpl<T, ClosureEdd<settings_E>, RadEnergyMoment::Energy>(rc);
   } else if (method == "mocmc") {
-    return RadConservedToPrimitiveFixupImpl<T, radiation::ClosureMOCMC<settings>>(rc);
+    if (do_num) RadConservedToPrimitiveFixupImpl<T, ClosureMOCMC<settings_N>, RadEnergyMoment::Number>(rc);
+    return RadConservedToPrimitiveFixupImpl<T, ClosureMOCMC<settings_E>, RadEnergyMoment::Energy>(rc);
   }
   return TaskStatus::fail;
 }
diff --git a/src/phoebus_driver.cpp b/src/phoebus_driver.cpp
index eaff0cbc..eb5dfa0a 100644
--- a/src/phoebus_driver.cpp
+++ b/src/phoebus_driver.cpp
@@ -54,7 +54,7 @@ namespace phoebus {
 PhoebusDriver::PhoebusDriver(ParameterInput *pin, ApplicationInput *app_in, Mesh *pm,
                              const bool is_restart)
     : EvolutionDriver(pin, app_in, pm),
-      integrator(std::make_unique<StagedIntegrator>(pin)), is_restart_(is_restart) {
+      integrator(std::make_unique<parthenon::LowStorageIntegrator>(pin)), is_restart_(is_restart) {
 
   // fail if these are not specified in the input file
   pin->CheckRequired("parthenon/mesh", "ix1_bc");
@@ -192,7 +192,12 @@ TaskCollection PhoebusDriver::RungeKuttaStage(const int stage) {
   bool rad_mocmc_active = false;
   if (rad_active) {
     rad_mocmc_active = (rad->Param<std::string>("method") == "mocmc");
+  } 
+  bool rad_moment_number_evolve = false;
+  if (rad_moments_active) {
+    rad_moment_number_evolve = rad->Param<bool>("do_number_evolution");
   }
+
   const auto monopole_enabled = monopole->Param<bool>("enable_monopole_gr");
   // Force static here means monopole only called at initialization.
   // and source terms are disabled
@@ -312,13 +317,22 @@ TaskCollection PhoebusDriver::RungeKuttaStage(const int stage) {
 
     if (rad_moments_active) {
       using MDT = std::remove_pointer<decltype(sc0.get())>::type;
-      auto moment_recon =
-          tl.AddTask(none, radiation::ReconstructEdgeStates<MDT>, sc0.get());
-      // TODO(BRR) Remove from task list if not using due to MOCMC
       auto get_opacities =
-          tl.AddTask(moment_recon, radiation::MomentCalculateOpacities<MDT>, sc0.get());
-      auto moment_flux =
-          tl.AddTask(get_opacities, radiation::CalculateFluxes<MDT>, sc0.get());
+          tl.AddTask(none, radiation::MomentCalculateOpacities<MDT>, sc0.get());
+      
+      auto moment_recon_E =
+          tl.AddTask(get_opacities, radiation::ReconstructEdgeStates<MDT>, sc0.get());
+      // TODO(BRR) Remove from task list if not using due to MOCMC
+      auto moment_flux_E =
+          tl.AddTask(moment_recon_E, radiation::CalculateFluxes<MDT>, sc0.get());
+      auto moment_recon = moment_recon_E;
+      auto moment_flux = moment_flux_E; 
+      if (rad_moment_number_evolve) { 
+        moment_recon =
+            tl.AddTask(moment_flux_E, radiation::ReconstructEdgeStates<MDT, radiation::RadEnergyMoment::Number>, sc0.get());
+        moment_flux =
+            tl.AddTask(moment_recon, radiation::CalculateFluxes<MDT, radiation::RadEnergyMoment::Number>, sc0.get());
+      }
       auto moment_geom_src = tl.AddTask(none, radiation::CalculateGeometricSource<MDT>,
                                         sc0.get(), gsrc.get());
       sndrcv_flux_depend = sndrcv_flux_depend | moment_flux;
diff --git a/src/phoebus_driver.hpp b/src/phoebus_driver.hpp
index e38105bf..da0d20e9 100644
--- a/src/phoebus_driver.hpp
+++ b/src/phoebus_driver.hpp
@@ -38,7 +38,7 @@ class PhoebusDriver : public EvolutionDriver {
   TaskListStatus Step();
 
  private:
-  std::unique_ptr<StagedIntegrator> integrator;
+  std::unique_ptr<parthenon::LowStorageIntegrator> integrator;
   const bool is_restart_;
   Real dt_init, dt_init_fact;
 
diff --git a/src/phoebus_utils/variables.hpp b/src/phoebus_utils/variables.hpp
index 5b5f1ef7..8d32a39b 100644
--- a/src/phoebus_utils/variables.hpp
+++ b/src/phoebus_utils/variables.hpp
@@ -36,11 +36,19 @@ constexpr char ye[] = "c.ye";
 namespace radmoment_prim {
 constexpr char J[] = "r.p.J";
 constexpr char H[] = "r.p.H";
+  namespace num {
+  constexpr char J[] = "r.p.n.J";
+  constexpr char H[] = "r.p.n.H";
+  } // namespace num
 } // namespace radmoment_prim
 
 namespace radmoment_cons {
 constexpr char E[] = "r.c.E";
 constexpr char F[] = "r.c.F";
+  namespace num {
+  constexpr char E[] = "r.c.n.E";
+  constexpr char F[] = "r.c.n.F";
+  } // namespace radmoment_cons
 } // namespace radmoment_cons
 
 namespace radmoment_internal {
@@ -58,6 +66,16 @@ constexpr char tilPi[] = "r.i.tilPi";
 constexpr char kappaH_mean[] = "r.i.kappaH_mean";
 constexpr char c2pfail[] = "r.i.c2p_fail";
 constexpr char srcfail[] = "r.i.src_fail";
+  namespace num { 
+  constexpr char xi[] = "r.i.n.xi";
+  constexpr char phi[] = "r.i.n.phi";
+  constexpr char dJ[] = "r.i.n.dJ";
+  constexpr char kappaJ[] = "r.i.n.kappaJ";
+  constexpr char kappaH[] = "r.i.n.kappaH";
+  constexpr char kappaH_mean[] = "r.i.n.kappaH_mean";
+  constexpr char JBB[] = "r.i.n.JBB";
+  constexpr char tilPi[] = "r.i.n.tilPi";
+  }
 } // namespace radmoment_internal
 
 namespace mocmc_internal {
diff --git a/src/radiation/closure.hpp b/src/radiation/closure.hpp
index aee022c5..fa52e29d 100644
--- a/src/radiation/closure.hpp
+++ b/src/radiation/closure.hpp
@@ -48,8 +48,8 @@ enum class ClosureVerbosity { quiet, v1, v2 };
 template <ClosureEquation EQ = ClosureEquation::energy_conserve,
           ClosureVerbosity VB = ClosureVerbosity::quiet>
 struct ClosureSettings {
-  static const ClosureEquation eqn_type = EQ;
-  static const ClosureVerbosity verbosity = VB;
+  static constexpr ClosureEquation eqn_type = EQ;
+  static constexpr ClosureVerbosity verbosity = VB;
 };
 
 enum class ClosureCon2PrimStrategy {
@@ -66,7 +66,7 @@ template <class SET = ClosureSettings<>>
 class ClosureEdd {
  public:
   using LocalGeometryType = LocalThreeGeometry;
-
+  static constexpr ClosureEquation eqn_type = SET::eqn_type;
   //-------------------------------------------------------------------------------------
   /// Constructor just calculates the inverse 3-metric, covariant three-velocity, and the
   /// Lorentz factor for the given background state.
diff --git a/src/radiation/moments.cpp b/src/radiation/moments.cpp
index 1a0859db..775e96fa 100644
--- a/src/radiation/moments.cpp
+++ b/src/radiation/moments.cpp
@@ -57,12 +57,8 @@ class ReconstructionIndexer {
   const int block_;
 };
 
-template <class T, class CLOSURE, bool STORE_GUESS>
+template <class T, class CLOSURE, bool STORE_GUESS, RadEnergyMoment EMOMENT>
 TaskStatus MomentCon2PrimImpl(T *rc) {
-  namespace cr = radmoment_cons;
-  namespace pr = radmoment_prim;
-  namespace ir = radmoment_internal;
-
   auto *pm = rc->GetParentPointer().get();
   StateDescriptor *rad = pm->packages.Get("radiation").get();
 
@@ -70,25 +66,27 @@ TaskStatus MomentCon2PrimImpl(T *rc) {
   IndexRange jb = pm->cellbounds.GetBoundsJ(IndexDomain::entire);
   IndexRange kb = pm->cellbounds.GetBoundsK(IndexDomain::entire);
 
+  auto var = RadiationVariableNames::GetByMomentType(EMOMENT);  
+  
   std::vector<std::string> variables{
-      cr::E,  cr::F,   pr::J,       pr::H,    fluid_prim::velocity,
-      ir::xi, ir::phi, ir::c2pfail, ir::tilPi};
+      var.E, var.F, var.J, var.H, fluid_prim::velocity,
+      var.xi, var.phi, var.c2pfail, var.tilPi};
   PackIndexMap imap;
   auto v = rc->PackVariables(variables, imap);
 
-  auto cE = imap.GetFlatIdx(cr::E);
-  auto pJ = imap.GetFlatIdx(pr::J);
-  auto cF = imap.GetFlatIdx(cr::F);
-  auto pH = imap.GetFlatIdx(pr::H);
+  auto cE = imap.GetFlatIdx(var.E);
+  auto pJ = imap.GetFlatIdx(var.J);
+  auto cF = imap.GetFlatIdx(var.F);
+  auto pH = imap.GetFlatIdx(var.H);
   auto pv = imap.GetFlatIdx(fluid_prim::velocity);
-  auto iTilPi = imap.GetFlatIdx(ir::tilPi, false);
+  auto iTilPi = imap.GetFlatIdx(var.tilPi, false);
   auto specB = cE.GetBounds(1);
   auto dirB = pH.GetBounds(2);
 
-  auto iXi = imap.GetFlatIdx(ir::xi);
-  auto iPhi = imap.GetFlatIdx(ir::phi);
+  auto iXi = imap.GetFlatIdx(var.xi);
+  auto iPhi = imap.GetFlatIdx(var.phi);
 
-  auto ifail = imap[ir::c2pfail].first;
+  auto ifail = imap[var.c2pfail].first;
 
   auto geom = Geometry::GetCoordinateSystem(rc);
   const Real pi = acos(-1);
@@ -124,7 +122,6 @@ TaskStatus MomentCon2PrimImpl(T *rc) {
         Vec covF = {{v(b, cF(ispec, 0), k, j, i) * isdetgam,
                      v(b, cF(ispec, 1), k, j, i) * isdetgam,
                      v(b, cF(ispec, 2), k, j, i) * isdetgam}};
-
         if (iTilPi.IsValid()) {
           SPACELOOP2(ii, jj) { conTilPi(ii, jj) = v(b, iTilPi(ispec, ii, jj), k, j, i); }
         } else {
@@ -171,15 +168,21 @@ TaskStatus MomentCon2Prim(T *rc) {
   auto *pm = rc->GetParentPointer().get();
   StateDescriptor *rad = pm->packages.Get("radiation").get();
   auto method = rad->Param<std::string>("method");
+  auto do_num = rad->Param<bool>("do_number_evolution");
 
-  using settings =
+  using settings_E =
       ClosureSettings<ClosureEquation::energy_conserve, ClosureVerbosity::quiet>;
+  using settings_N =
+      ClosureSettings<ClosureEquation::number_conserve, ClosureVerbosity::quiet>;
   if (method == "moment_m1") {
-    return MomentCon2PrimImpl<T, ClosureM1<settings>, true>(rc);
+    if (do_num) MomentCon2PrimImpl<T, ClosureM1<settings_N>, true, RadEnergyMoment::Number>(rc);
+    return MomentCon2PrimImpl<T, ClosureM1<settings_E>, true, RadEnergyMoment::Energy>(rc);
   } else if (method == "moment_eddington") {
-    return MomentCon2PrimImpl<T, ClosureEdd<settings>, false>(rc);
+    if (do_num) MomentCon2PrimImpl<T, ClosureEdd<settings_N>, false, RadEnergyMoment::Number>(rc);
+    return MomentCon2PrimImpl<T, ClosureEdd<settings_E>, false, RadEnergyMoment::Energy>(rc);
   } else if (method == "mocmc") {
-    return MomentCon2PrimImpl<T, ClosureMOCMC<settings>, false>(rc);
+    if (do_num) MomentCon2PrimImpl<T, ClosureMOCMC<settings_N>, false, RadEnergyMoment::Number>(rc);
+    return MomentCon2PrimImpl<T, ClosureMOCMC<settings_E>, false, RadEnergyMoment::Energy>(rc);
   } else {
     PARTHENON_FAIL("Radiation method unknown");
   }
@@ -188,31 +191,28 @@ TaskStatus MomentCon2Prim(T *rc) {
 // template TaskStatus MomentCon2Prim<MeshData<Real>>(MeshData<Real> *);
 template TaskStatus MomentCon2Prim<MeshBlockData<Real>>(MeshBlockData<Real> *);
 
-template <class T, class CLOSURE>
+template <class T, class CLOSURE, RadEnergyMoment EMOMENT>
 TaskStatus MomentPrim2ConImpl(T *rc, IndexDomain domain) {
-  namespace cr = radmoment_cons;
-  namespace pr = radmoment_prim;
-  namespace ir = radmoment_internal;
-
   auto *pm = rc->GetParentPointer().get();
 
   IndexRange ib = pm->cellbounds.GetBoundsI(domain);
   IndexRange jb = pm->cellbounds.GetBoundsJ(domain);
   IndexRange kb = pm->cellbounds.GetBoundsK(domain);
-
-  std::vector<std::string> variables{cr::E, cr::F, pr::J, pr::H, fluid_prim::velocity};
+  
+  auto var = RadiationVariableNames::GetByMomentType(EMOMENT);  
+  std::vector<std::string> variables{var.E, var.F, var.J, var.H, fluid_prim::velocity};
   if (programming::is_specialization_of<CLOSURE, ClosureMOCMC>::value) {
-    variables.push_back(ir::tilPi);
+    variables.push_back(var.tilPi);
   }
   PackIndexMap imap;
   auto v = rc->PackVariables(variables, imap);
 
-  auto cE = imap.GetFlatIdx(cr::E);
-  auto pJ = imap.GetFlatIdx(pr::J);
-  auto cF = imap.GetFlatIdx(cr::F);
-  auto pH = imap.GetFlatIdx(pr::H);
+  auto cE = imap.GetFlatIdx(var.E);
+  auto pJ = imap.GetFlatIdx(var.J);
+  auto cF = imap.GetFlatIdx(var.F);
+  auto pH = imap.GetFlatIdx(var.H);
   auto pv = imap.GetFlatIdx(fluid_prim::velocity);
-  auto iTilPi = imap.GetFlatIdx(ir::tilPi, false);
+  auto iTilPi = imap.GetFlatIdx(var.tilPi, false);
 
   auto specB = cE.GetBounds(1);
   auto dirB = pH.GetBounds(2);
@@ -256,7 +256,6 @@ TaskStatus MomentPrim2ConImpl(T *rc, IndexDomain domain) {
         PARTHENON_DEBUG_REQUIRE(!std::isnan(J), "NAN J in rad P2C!");
 
         c.Prim2Con(J, covH, conTilPi, &E, &covF);
-
         v(b, cE(ispec), k, j, i) = sdetgam * E;
         for (int idir = dirB.s; idir <= dirB.e; ++idir) {
           v(b, cF(ispec, idir), k, j, i) = sdetgam * covF(idir);
@@ -271,14 +270,21 @@ TaskStatus MomentPrim2Con(T *rc, IndexDomain domain) {
   auto *pm = rc->GetParentPointer().get();
   StateDescriptor *rad = pm->packages.Get("radiation").get();
   auto method = rad->Param<std::string>("method");
-  using settings =
+  auto do_num = rad->Param<bool>("do_number_evolution");
+
+  using settings_N =
+      ClosureSettings<ClosureEquation::number_conserve, ClosureVerbosity::quiet>;
+  using settings_E =
       ClosureSettings<ClosureEquation::energy_conserve, ClosureVerbosity::quiet>;
   if (method == "moment_m1") {
-    return MomentPrim2ConImpl<T, ClosureM1<settings>>(rc, domain);
+    if (do_num) MomentPrim2ConImpl<T, ClosureM1<settings_N>, RadEnergyMoment::Number>(rc, domain);
+    return MomentPrim2ConImpl<T, ClosureM1<settings_E>, RadEnergyMoment::Energy>(rc, domain);
   } else if (method == "moment_eddington") {
-    return MomentPrim2ConImpl<T, ClosureEdd<settings>>(rc, domain);
+    if (do_num) MomentPrim2ConImpl<T, ClosureEdd<settings_N>, RadEnergyMoment::Number>(rc, domain);
+    return MomentPrim2ConImpl<T, ClosureEdd<settings_E>, RadEnergyMoment::Energy>(rc, domain);
   } else if (method == "mocmc") {
-    return MomentPrim2ConImpl<T, ClosureMOCMC<settings>>(rc, domain);
+    if (do_num) MomentPrim2ConImpl<T, ClosureMOCMC<settings_N>, RadEnergyMoment::Number>(rc, domain);
+    return MomentPrim2ConImpl<T, ClosureMOCMC<settings_E>, RadEnergyMoment::Energy>(rc, domain);
   } else {
     PARTHENON_FAIL("Radiation method unknown!");
   }
@@ -288,7 +294,7 @@ TaskStatus MomentPrim2Con(T *rc, IndexDomain domain) {
 template TaskStatus MomentPrim2Con<MeshBlockData<Real>>(MeshBlockData<Real> *,
                                                         IndexDomain);
 
-template <class T>
+template <class T, RadEnergyMoment EMOMENT>
 TaskStatus ReconstructEdgeStates(T *rc) {
   using namespace PhoebusReconstruction;
 
@@ -304,27 +310,24 @@ TaskStatus ReconstructEdgeStates(T *rc) {
   IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior);
   const int dk = (pmb->pmy_mesh->ndim > 2 ? 1 : 0);
 
-  namespace cr = radmoment_cons;
-  namespace pr = radmoment_prim;
-  namespace ir = radmoment_internal;
-
+  auto var = RadiationVariableNames::GetByMomentType(EMOMENT);  
   PackIndexMap imap_ql, imap_qr, imap;
   VariablePack<Real> ql_base =
-      rc->PackVariables(std::vector<std::string>{ir::ql}, imap_ql);
+      rc->PackVariables(std::vector<std::string>{var.ql}, imap_ql);
   VariablePack<Real> qr_base =
-      rc->PackVariables(std::vector<std::string>{ir::qr}, imap_qr);
-  std::vector<std::string> variables = {pr::J, pr::H, ir::tilPi};
-  variables.push_back(ir::dJ);
+      rc->PackVariables(std::vector<std::string>{var.qr}, imap_qr);
+  
+  std::vector<std::string> variables = {var.J, var.H, var.tilPi, var.dJ};
   VariablePack<Real> v = rc->PackVariables(variables, imap);
-  auto idx_J = imap.GetFlatIdx(pr::J);
-  auto idx_dJ = imap.GetFlatIdx(ir::dJ);
-  auto iTilPi = imap.GetFlatIdx(ir::tilPi, false);
+  auto idx_J = imap.GetFlatIdx(var.J);
+  auto idx_dJ = imap.GetFlatIdx(var.dJ);
+  auto iTilPi = imap.GetFlatIdx(var.tilPi, false);
 
-  ParArrayND<Real> ql_v = rc->Get(ir::ql_v).data;
-  ParArrayND<Real> qr_v = rc->Get(ir::qr_v).data;
+  ParArrayND<Real> ql_v = rc->Get(var.ql_v).data;
+  ParArrayND<Real> qr_v = rc->Get(var.qr_v).data;
   VariablePack<Real> v_vel =
       rc->PackVariables(std::vector<std::string>{fluid_prim::velocity});
-  auto qIdx = imap_ql.GetFlatIdx(ir::ql);
+  auto qIdx = imap_ql.GetFlatIdx(var.ql);
 
   const int nspec = qIdx.DimSize(1);
   int nrecon = 4 * nspec;
@@ -332,7 +335,7 @@ TaskStatus ReconstructEdgeStates(T *rc) {
     nrecon = (4 + 9) * nspec; // TODO(BRR) 6 instead of 9 for conTilPi by symmetry
   }
 
-  const int offset = imap_ql[ir::ql].first;
+  const int offset = imap_ql[var.ql].first;
 
   const int nblock = ql_base.GetDim(5);
   const int ndim = pmb->pmy_mesh->ndim;
@@ -485,11 +488,12 @@ TaskStatus ReconstructEdgeStates(T *rc) {
 
   return TaskStatus::complete;
 }
-template TaskStatus ReconstructEdgeStates<MeshBlockData<Real>>(MeshBlockData<Real> *);
+template TaskStatus ReconstructEdgeStates<MeshBlockData<Real>, RadEnergyMoment::Number>(MeshBlockData<Real> *);
+template TaskStatus ReconstructEdgeStates<MeshBlockData<Real>, RadEnergyMoment::Energy>(MeshBlockData<Real> *);
 
 // This really only works for MeshBlockData right now since fluxes don't have a block
 // index
-template <class T, class CLOSURE>
+template <class T, class CLOSURE, RadEnergyMoment EMOMENT>
 TaskStatus CalculateFluxesImpl(T *rc) {
   auto *pmb = rc->GetParentPointer().get();
   StateDescriptor *rad_pkg = pmb->packages.Get("radiation").get();
@@ -504,25 +508,23 @@ TaskStatus CalculateFluxesImpl(T *rc) {
   IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior);
   const int dk = (pmb->pmy_mesh->ndim > 2 ? 1 : 0);
 
-  namespace cr = radmoment_cons;
-  namespace pr = radmoment_prim;
-  namespace ir = radmoment_internal;
+  auto var = RadiationVariableNames::GetByMomentType(EMOMENT); 
 
   PackIndexMap imap_ql, imap_qr, imap;
-  std::vector<std::string> vars{ir::ql, ir::qr, ir::ql_v, ir::qr_v, ir::dJ, ir::kappaH};
-  std::vector<std::string> flxs{cr::E, cr::F};
+  std::vector<std::string> vars{var.ql, var.qr, var.ql_v, var.qr_v, var.dJ, var.kappaH};
+  std::vector<std::string> flxs{var.E, var.F};
 
   auto v = rc->PackVariablesAndFluxes(vars, flxs, imap);
 
-  auto idx_qlv = imap.GetFlatIdx(ir::ql_v);
-  auto idx_qrv = imap.GetFlatIdx(ir::qr_v);
-  auto idx_ql = imap.GetFlatIdx(ir::ql);
-  auto idx_qr = imap.GetFlatIdx(ir::qr);
-  auto idx_dJ = imap.GetFlatIdx(ir::dJ);
-  auto idx_kappaH = imap.GetFlatIdx(ir::kappaH);
+  auto idx_qlv = imap.GetFlatIdx(var.ql_v);
+  auto idx_qrv = imap.GetFlatIdx(var.qr_v);
+  auto idx_ql = imap.GetFlatIdx(var.ql);
+  auto idx_qr = imap.GetFlatIdx(var.qr);
+  auto idx_dJ = imap.GetFlatIdx(var.dJ);
+  auto idx_kappaH = imap.GetFlatIdx(var.kappaH);
 
-  auto idx_Ff = imap.GetFlatIdx(cr::F);
-  auto idx_Ef = imap.GetFlatIdx(cr::E);
+  auto idx_Ff = imap.GetFlatIdx(var.F);
+  auto idx_Ef = imap.GetFlatIdx(var.E);
 
   auto num_species = rad_pkg->Param<int>("num_species");
 
@@ -746,25 +748,27 @@ TaskStatus CalculateFluxesImpl(T *rc) {
   return TaskStatus::complete;
 }
 
-template <class T>
+template <class T, RadEnergyMoment EMOMENT>
 TaskStatus CalculateFluxes(T *rc) {
   auto *pm = rc->GetParentPointer().get();
   StateDescriptor *rad = pm->packages.Get("radiation").get();
   auto method = rad->Param<std::string>("method");
-  using settings =
-      ClosureSettings<ClosureEquation::energy_conserve, ClosureVerbosity::quiet>;
+  using settings = typename std::conditional<EMOMENT == RadEnergyMoment::Number, 
+                                    ClosureSettings<ClosureEquation::number_conserve, ClosureVerbosity::quiet>,
+                                    ClosureSettings<ClosureEquation::energy_conserve, ClosureVerbosity::quiet>>::type;
   if (method == "moment_m1") {
-    return CalculateFluxesImpl<T, ClosureM1<settings>>(rc);
+    return CalculateFluxesImpl<T, ClosureM1<settings>, EMOMENT>(rc);
   } else if (method == "moment_eddington") {
-    return CalculateFluxesImpl<T, ClosureEdd<settings>>(rc);
+    return CalculateFluxesImpl<T, ClosureEdd<settings>, EMOMENT>(rc);
   } else if (method == "mocmc") {
-    return CalculateFluxesImpl<T, ClosureMOCMC<settings>>(rc);
+    return CalculateFluxesImpl<T, ClosureMOCMC<settings>, EMOMENT>(rc);
   } else {
     PARTHENON_FAIL("Radiation method unknown!");
   }
   return TaskStatus::fail;
 }
-template TaskStatus CalculateFluxes<MeshBlockData<Real>>(MeshBlockData<Real> *);
+template TaskStatus CalculateFluxes<MeshBlockData<Real>, RadEnergyMoment::Number>(MeshBlockData<Real> *);
+template TaskStatus CalculateFluxes<MeshBlockData<Real>, RadEnergyMoment::Energy>(MeshBlockData<Real> *);
 
 template <class T, class CLOSURE>
 TaskStatus CalculateGeometricSourceImpl(T *rc, T *rc_src) {
@@ -872,17 +876,20 @@ TaskStatus CalculateGeometricSourceImpl(T *rc, T *rc_src) {
             }
           }
           SPACELOOP2(ii, jj) { con_T[ii + 1][jj + 1] += conP(ii, jj); }
-
-          Real TGam = 0.0;
-          SPACETIMELOOP2(m, n) {
-            Real gam0 = 0.;
-            SPACETIMELOOP(r) { gam0 += con_g[0][r] * Gamma[r][m][n]; }
-            TGam += con_T[m][n] * gam0;
+          
+          if constexpr(CLOSURE::eqn_type == ClosureEquation::energy_conserve) {
+            Real TGam = 0.0;
+            SPACETIMELOOP2(m, n) {
+              Real gam0 = 0.;
+              SPACETIMELOOP(r) { gam0 += con_g[0][r] * Gamma[r][m][n]; }
+              TGam += con_T[m][n] * gam0;
+            }
+            Real Ta = 0.0;
+            SPACETIMELOOP(m) { Ta += con_T[m][0] * dlnalp[m]; }
+            v_src(iblock, idx_E_src(ispec), k, j, i) = sdetgam * alp * alp * (Ta - TGam);
+          } else { 
+            v_src(iblock, idx_E_src(ispec), k, j, i) = 0.0;
           }
-          Real Ta = 0.0;
-          SPACETIMELOOP(m) { Ta += con_T[m][0] * dlnalp[m]; }
-          v_src(iblock, idx_E_src(ispec), k, j, i) = sdetgam * alp * alp * (Ta - TGam);
-
           SPACELOOP(l) {
             Real src_mom = 0.0;
             SPACETIMELOOP2(m, n) {
diff --git a/src/radiation/radiation.cpp b/src/radiation/radiation.cpp
index 0ad1327c..89becfc8 100644
--- a/src/radiation/radiation.cpp
+++ b/src/radiation/radiation.cpp
@@ -424,13 +424,27 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
 
     physics->AddField(c::E, mspecies_scalar_cons);
     physics->AddField(c::F, mspecies_three_vector_cons);
-
+    
     physics->AddField(p::J, mspecies_scalar);
     physics->AddField(p::H, mspecies_three_vector);
 
     // Fields for saving guesses for NR iteration in the radiation Con2Prim type solve
     physics->AddField(i::xi, mspecies_scalar);
     physics->AddField(i::phi, mspecies_scalar);
+    
+    const bool do_num_evol =
+        pin->GetOrAddBoolean("radiation", "do_number_evolution", true);
+    params.Add("do_number_evolution", do_num_evol);
+    if (do_num_evol) { 
+      physics->AddField(c::num::E, mspecies_scalar_cons);
+      physics->AddField(c::num::F, mspecies_three_vector_cons);
+    
+      physics->AddField(p::num::J, mspecies_scalar);
+      physics->AddField(p::num::H, mspecies_three_vector);
+
+      physics->AddField(i::num::xi, mspecies_scalar);
+      physics->AddField(i::num::phi, mspecies_scalar);
+    }
 
     // Fields for cell edge reconstruction
     /// TODO: (LFR) The amount of storage can likely be reduced, but maybe at the expense
@@ -452,6 +466,7 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
     Metadata mdJ = Metadata({Metadata::Cell, Metadata::Derived, Metadata::OneCopy},
                             std::vector<int>{num_species, ndim, ndim});
     physics->AddField(i::dJ, mdJ);
+    if (do_num_evol) physics->AddField(i::num::dJ, mdJ); 
 
     // Add variables for source functions
     Metadata mSourceVar = Metadata({Metadata::Cell, Metadata::Derived, Metadata::OneCopy},
@@ -459,7 +474,11 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
     physics->AddField(i::kappaJ, mSourceVar);
     physics->AddField(i::kappaH, mSourceVar);
     physics->AddField(i::JBB, mSourceVar);
-
+    if (do_num_evol) {
+      physics->AddField(i::num::kappaJ, mSourceVar);
+      physics->AddField(i::num::kappaH, mSourceVar);
+      physics->AddField(i::num::JBB, mSourceVar);
+    }
     // this fail flag should really be an enum or something
     // but parthenon doesn't yet support that kind of thing
     Metadata m_scalar = Metadata({Metadata::Cell, Metadata::OneCopy, Metadata::Derived,
@@ -474,6 +493,7 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
           std::vector<int>{num_species, 3, 3});
 
       physics->AddField(i::tilPi, mspecies_three_tensor);
+      if (do_num_evol) physics->AddField(i::num::tilPi, mspecies_three_tensor);
       physics->AddField(mocmc_internal::dnsamp, mscalar);
     }
 
diff --git a/src/radiation/radiation.hpp b/src/radiation/radiation.hpp
index 6b02f758..247de959 100644
--- a/src/radiation/radiation.hpp
+++ b/src/radiation/radiation.hpp
@@ -123,16 +123,18 @@ TaskStatus InitializeCommunicationMesh(const std::string swarmName,
 Real EstimateTimestepBlock(MeshBlockData<Real> *rc);
 
 // Moment tasks
+enum class RadEnergyMoment {Number, Energy};
+
 template <class T>
 TaskStatus MomentCon2Prim(T *rc);
 
 template <class T>
 TaskStatus MomentPrim2Con(T *rc, IndexDomain domain = IndexDomain::entire);
 
-template <class T>
+template <class T, RadEnergyMoment EMOMENT = RadEnergyMoment::Energy>
 TaskStatus ReconstructEdgeStates(T *rc);
 
-template <class T>
+template <class T, RadEnergyMoment EMOMENT = RadEnergyMoment::Energy>
 TaskStatus CalculateFluxes(T *rc);
 
 template <class T>
@@ -165,6 +167,52 @@ TaskStatus MOCMCFluidSource(T *rc, const Real dt, const bool update_fluid);
 
 TaskStatus MOCMCUpdateParticleCount(Mesh *pmesh, std::vector<Real> *resolution);
 
+struct RadiationVariableNames {
+  static RadiationVariableNames GetByMomentType(RadEnergyMoment type) {
+    namespace cr = radmoment_cons;
+    namespace pr = radmoment_prim;
+    namespace ir = radmoment_internal; 
+    if (type == RadEnergyMoment::Number) { 
+      return RadiationVariableNames{cr::num::E, cr::num::F, 
+                                    pr::num::J, pr::num::H, 
+                                    ir::num::xi, ir::num::phi,
+                                    ir::ql, ir::qr,
+                                    ir::ql_v, ir::qr_v, 
+                                    ir::num::dJ, ir::num::kappaJ, ir::num::kappaH, 
+                                    ir::num::JBB, ir::num::tilPi, ir::num::kappaH_mean,
+                                    ir::c2pfail, ir::srcfail};
+    } else /*if (type == RadEnergyMoment::Energy)*/ { 
+      return RadiationVariableNames{cr::E, cr::F, 
+                                    pr::J, pr::H, 
+                                    ir::xi, ir::phi, 
+                                    ir::ql, ir::qr,
+                                    ir::ql_v, ir::qr_v, 
+                                    ir::dJ, ir::kappaJ, ir::kappaH, 
+                                    ir::JBB, ir::tilPi, ir::kappaH_mean,
+                                    ir::c2pfail, ir::srcfail};
+    }
+  }
+  const std::string& E; 
+  const std::string& F; 
+  const std::string& J;
+  const std::string& H; 
+  
+  const std::string& xi; 
+  const std::string& phi; 
+  const std::string& ql; 
+  const std::string& qr; 
+  const std::string& ql_v; 
+  const std::string& qr_v; 
+  const std::string& dJ; 
+  const std::string& kappaJ; 
+  const std::string& kappaH; 
+  const std::string& JBB; 
+  const std::string& tilPi;
+  const std::string& kappaH_mean; 
+  const std::string& c2pfail; 
+  const std::string& srcfail; 
+};
+
 } // namespace radiation
 
 #endif // RADIATION_HPP_