diff --git a/ltepop.cc b/ltepop.cc
index 4ab837c77..e81ae5347 100644
--- a/ltepop.cc
+++ b/ltepop.cc
@@ -451,7 +451,7 @@ static auto find_uppermost_ion(const int modelgridindex, const int element, cons
   int uppermost_ion = 0;
 
   uppermost_ion = nions - 1;
-  if (!force_lte) {
+  if (!use_lte) {
     for (int ion = 0; ion < nions - 1; ion++) {
       if (iongamma_is_zero(nonemptymgi, element, ion) &&
           (!NT_ON || ((globals::dep_estimator_gamma[nonemptymgi] == 0.) &&
@@ -504,8 +504,10 @@ void set_groundlevelpops(const int modelgridindex, const int element, const floa
   /// calculate number density of the current element (abundances are given by mass)
   const double nnelement = grid::get_elem_numberdens(modelgridindex, element);
 
+  const bool use_phi_lte = force_lte || FORCE_SAHA_ION_BALANCE(get_atomicnumber(element));
+
   const auto ionfractions =
-      (nnelement > 0) ? calculate_ionfractions(element, modelgridindex, nne, force_lte) : std::vector<double>();
+      (nnelement > 0) ? calculate_ionfractions(element, modelgridindex, nne, use_phi_lte) : std::vector<double>();
 
   const int uppermost_ion = static_cast<int>(ionfractions.size() - 1);
 
diff --git a/nltepop.cc b/nltepop.cc
index cbfdc8429..ad5971262 100644
--- a/nltepop.cc
+++ b/nltepop.cc
@@ -1108,6 +1108,9 @@ void nltepop_write_to_file(const int modelgridindex, const int timestep) {
   //         timestep, n, grid::get_TR(n), grid::get_Te(n), grid::get_W(n), grid::get_TJ(n), grid::get_nne(n));
 
   for (int element = 0; element < get_nelements(); element++) {
+    if (!elem_has_nlte_levels(element)) {
+      continue;
+    }
     const int nions = get_nions(element);
     const int atomic_number = get_atomicnumber(element);
 
diff --git a/update_grid.cc b/update_grid.cc
index 7aa7add90..63c581cf3 100644
--- a/update_grid.cc
+++ b/update_grid.cc
@@ -733,7 +733,7 @@ void solve_Te_nltepops(const int mgi, const int nonemptymgi, const int nts, cons
       // population change))
       double fracdiff_nne = 0.;
       for (int element = 0; element < get_nelements(); element++) {
-        if (get_nions(element) > 0) {
+        if (get_nions(element) > 0 && elem_has_nlte_levels(element)) {
           solve_nlte_pops_element(element, mgi, nts, nlte_iter);
           calculate_cellpartfuncts(mgi, element);
         }