diff --git a/artisoptions_christinenonthermal.h b/artisoptions_christinenonthermal.h
index 4d5b917af..d918598d6 100644
--- a/artisoptions_christinenonthermal.h
+++ b/artisoptions_christinenonthermal.h
@@ -26,6 +26,8 @@ constexpr bool LEVEL_IS_NLTE(int element_z, int ionstage, int level) {
 
 constexpr bool LTEPOP_EXCITATION_USE_TJ = false;
 
+constexpr bool USE_TE_SOLVER = true; /// If not solving for Te then Te=TJ
+
 constexpr bool FORCE_SAHA_ION_BALANCE(int element_z) { return false; }
 
 constexpr bool single_level_top_ion = false;
diff --git a/artisoptions_classic.h b/artisoptions_classic.h
index 1900ec104..b5abf7e77 100644
--- a/artisoptions_classic.h
+++ b/artisoptions_classic.h
@@ -21,6 +21,8 @@ constexpr bool LEVEL_IS_NLTE(int element_z, int ionstage, int level) { return fa
 
 constexpr bool LTEPOP_EXCITATION_USE_TJ = true;
 
+constexpr bool USE_TE_SOLVER = true; /// If not solving for Te then Te=TJ
+
 constexpr bool FORCE_SAHA_ION_BALANCE(int element_z) { return false; }
 
 constexpr bool single_level_top_ion = true;
diff --git a/artisoptions_doc.md b/artisoptions_doc.md
index 3f02058bf..867fbdddc 100644
--- a/artisoptions_doc.md
+++ b/artisoptions_doc.md
@@ -24,6 +24,10 @@ constexpr bool LEVEL_IS_NLTE(int element_z, int ionstage, int level) { return fa
 // This is default on for classic, and off for nebularnlte, where it affects the super-level
 constexpr bool LTEPOP_EXCITATION_USE_TJ = false;
 
+//Switch to use Te solver to calculate Te from heating and cooling rates. Use Te=TJ instead when heating/cooling rates 
+// are too noisy to get a good Te solution (e.g. when only a few species are in NLTE)
+constexpr bool USE_TE_SOLVER = true; /// If not solving for Te then Te=TJ
+
 // Only include a single level for the highest ion stage
 constexpr bool single_level_top_ion;
 
diff --git a/artisoptions_kilonova_lte.h b/artisoptions_kilonova_lte.h
index 3d38dda29..66ac7aadf 100644
--- a/artisoptions_kilonova_lte.h
+++ b/artisoptions_kilonova_lte.h
@@ -21,6 +21,8 @@ constexpr bool LEVEL_IS_NLTE(int element_z, int ionstage, int level) { return fa
 
 constexpr bool LTEPOP_EXCITATION_USE_TJ = true;
 
+constexpr bool USE_TE_SOLVER = true; /// If not solving for Te then Te=TJ
+
 constexpr bool FORCE_SAHA_ION_BALANCE(int element_z) { return true; }
 
 constexpr bool single_level_top_ion = false;
diff --git a/artisoptions_nltenebular.h b/artisoptions_nltenebular.h
index 8194ad7db..223917bc4 100644
--- a/artisoptions_nltenebular.h
+++ b/artisoptions_nltenebular.h
@@ -26,6 +26,8 @@ constexpr bool LEVEL_IS_NLTE(int element_z, int ionstage, int level) {
 
 constexpr bool LTEPOP_EXCITATION_USE_TJ = false;
 
+constexpr bool USE_TE_SOLVER = true; /// If not solving for Te then Te=TJ
+
 constexpr bool FORCE_SAHA_ION_BALANCE(int element_z) { return false; }
 
 constexpr bool single_level_top_ion = false;
diff --git a/artisoptions_nltewithoutnonthermal.h b/artisoptions_nltewithoutnonthermal.h
index e144d015c..1c5b37ee6 100644
--- a/artisoptions_nltewithoutnonthermal.h
+++ b/artisoptions_nltewithoutnonthermal.h
@@ -26,6 +26,8 @@ constexpr bool LEVEL_IS_NLTE(int element_z, int ionstage, int level) {
 
 constexpr bool LTEPOP_EXCITATION_USE_TJ = false;
 
+constexpr bool USE_TE_SOLVER = true; /// If not solving for Te then Te=TJ
+
 constexpr bool FORCE_SAHA_ION_BALANCE(int element_z) { return false; }
 
 constexpr bool single_level_top_ion = true;
diff --git a/update_grid.cc b/update_grid.cc
index dd95df759..ec075caf3 100644
--- a/update_grid.cc
+++ b/update_grid.cc
@@ -714,11 +714,14 @@ void solve_Te_nltepops(const int mgi, const int nonemptymgi, const int nts, cons
 
     const double prev_T_e = grid::get_Te(mgi);
     const auto sys_time_start_Te = std::time(nullptr);
-    const int nts_for_te = (titer == 0) ? nts - 1 : nts;
 
-    /// Find T_e as solution for thermal balance
-    call_T_e_finder(mgi, nts, globals::timesteps[nts_for_te].mid, MINTEMP, MAXTEMP, heatingcoolingrates,
-                    bfheatingcoeffs);
+    if (USE_TE_SOLVER) {
+      const int nts_for_te = (titer == 0) ? nts - 1 : nts;
+
+      /// Find T_e as solution for thermal balance
+      call_T_e_finder(mgi, nts, globals::timesteps[nts_for_te].mid, MINTEMP, MAXTEMP, heatingcoolingrates,
+                      bfheatingcoeffs);
+    }
 
     const int duration_solve_T_e = std::time(nullptr) - sys_time_start_Te;
 
@@ -1024,6 +1027,11 @@ void update_grid_cell(const int mgi, const int nts, const int nts_prev, const in
       // full-spectrum and binned J and nuJ estimators
       radfield::fit_parameters(mgi, nts);
 
+      if (!USE_TE_SOLVER){
+        const double T_J = grid::get_TJ(mgi);
+        grid::set_Te(mgi, T_J); // instead of Te from heating/cooling rates
+      }
+
       solve_Te_nltepops(mgi, nonemptymgi, nts, titer, heatingcoolingrates);
     }
     printout("Temperature/NLTE solution for cell %d timestep %d took %ld seconds\n", mgi, nts,