diff --git a/src/analysis/history.hpp b/src/analysis/history.hpp
index 645dcda1..62c2a610 100644
--- a/src/analysis/history.hpp
+++ b/src/analysis/history.hpp
@@ -112,6 +112,7 @@ Real ReduceInGain(MeshData<Real> *md, const std::string &varname, int idx = 0) {
       std::is_same<Reducer_t, Kokkos::Sum<Real, HostExecSpace>>::value ||
       std::is_same<Reducer_t, Kokkos::Sum<Real, DevExecSpace>>::value ||
       std::is_same<Reducer_t, Kokkos::Sum<Real, Kokkos::DefaultExecutionSpace>>::value;
+  auto geom = Geometry::GetCoordinateSystem(md);
 
   parthenon::par_reduce(
       parthenon::LoopPatternMDRange(), "Phoebus History for " + varname, DevExecSpace(),
@@ -120,10 +121,12 @@ Real ReduceInGain(MeshData<Real> *md, const std::string &varname, int idx = 0) {
         // join is a Kokkos construct
         // that automatically does the
         // reduction operation locally
+        Real gdet = geom.DetGamma(CellLocation::Cent, 0, k, j, i);
         bool is_netheat = (pack(b, iheat, k, j, i) - pack(b, icool, k, j, i) > 0);
         const auto &coords = pack.GetCoords(b);
         const Real vol = volume_weighting ? coords.CellVolume(k, j, i) : 1.0;
-        reducer.join(lresult, pack(b, ivar.first + idx, k, j, i) * vol * is_netheat);
+        reducer.join(lresult,
+                     pack(b, ivar.first + idx, k, j, i) * gdet * vol * is_netheat);
       },
       reducer);
   return result;
diff --git a/src/fluid/fluid.cpp b/src/fluid/fluid.cpp
index 83462166..6fd709aa 100644
--- a/src/fluid/fluid.cpp
+++ b/src/fluid/fluid.cpp
@@ -311,6 +311,7 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   // Reductions
   // By default compute integrated value of scalar conserved vars
   auto HstSum = parthenon::UserHistoryOperation::sum;
+  auto HstMax = parthenon::UserHistoryOperation::max;
   using History::ReduceInGain;
   using History::ReduceOneVar;
   using parthenon::HistoryOutputVar;
@@ -322,7 +323,7 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   auto ReduceEn = [](MeshData<Real> *md) {
     return ReduceOneVar<Kokkos::Sum<Real>>(md, fluid_cons::energy::name(), 0);
   };
-  auto CentralDensitySN = [](MeshData<Real> *md) {
+  auto MaxDensitySN = [](MeshData<Real> *md) {
     History::ReduceCentralDensitySN(md);
     return ReduceOneVar<Kokkos::Max<Real>>(md, diag::central_density::name(), 0);
   };
@@ -330,9 +331,16 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
     return ReduceInGain<Kokkos::Sum<Real, HostExecSpace>>(md, fluid_prim::density::name(),
                                                           0);
   };
+  auto Qgain = [](MeshData<Real> *md) {
+    return ReduceInGain<Kokkos::Sum<Real, HostExecSpace>>(
+               md, internal_variables::GcovHeat::name(), 0) -
+           ReduceInGain<Kokkos::Sum<Real, HostExecSpace>>(
+               md, internal_variables::GcovCool::name(), 0);
+  };
 
-  hst_vars.emplace_back(HistoryOutputVar(HstSum, CentralDensitySN, "central density SN"));
+  hst_vars.emplace_back(HistoryOutputVar(HstMax, MaxDensitySN, "maximum density SN"));
   hst_vars.emplace_back(HistoryOutputVar(HstSum, Mgain, "Mgain"));
+  hst_vars.emplace_back(HistoryOutputVar(HstSum, Qgain, "total net heat"));
   hst_vars.emplace_back(HistoryOutputVar(HstSum, ReduceMass, "total baryon number"));
   hst_vars.emplace_back(HistoryOutputVar(HstSum, ReduceEn, "total conserved energy tau"));