diff --git a/src/m_fluid.f90 b/src/m_fluid.f90
index 0eebf3f8..ff9a0eb3 100644
--- a/src/m_fluid.f90
+++ b/src/m_fluid.f90
@@ -445,16 +445,21 @@ subroutine add_source_terms(box, dt, n_vars, i_cc, s_deriv, s_out, n_dt, dt_lim,
           box%cc(IJK, i_electron_energy+s_out) = box%cc(IJK, i_electron_energy+s_out) &
                + dt * (gain - loss_rate * box%cc(IJK, i_electron+s_out))
        end if
+    end do; CLOSE_DO
 
-       do n = n_gas_species+1, n_species
-          iv = species_itree(n)
+    do n = n_gas_species+1, n_species
+       ix = 0
+       iv = species_itree(n)
+       do KJI_DO(1,nc)
+          ix = ix + 1
+          if (.not. mask(IJK)) cycle
           box%cc(IJK, iv+s_out) = box%cc(IJK, iv+s_out) + dt * derivs(ix, n)
-       end do
-    end do; CLOSE_DO
+       end do; CLOSE_DO
+    end do
 
-    tmp = maxval(mean_energies)
-    if (tmp > 0) then
+    if (model_has_energy_equation) then
        ! Set time step restriction for energy loss
+       tmp = maxval(mean_energies)
        dt_lim(2) = tmp/LT_get_col(td_ee_tbl, td_ee_loss, tmp)
     end if