Update nltepop.cc

nltepop.cc

@@ -600,9 +600,10 @@ void nltepop_matrix_add_nt_ionisation(const int nonemptymgi, const int element,
       for (int level = 0; level < nlevels; level++) {
         const int lower_index = get_nlte_vector_index(element, ion, level);
 
-        *gsl_matrix_ptr(rate_matrix_ntcoll_bf, lower_index, lower_index) -= Y_nt_thisupperion * s_renorm[level];
-        *gsl_matrix_ptr(rate_matrix_ntcoll_bf, upper_groundstate_index, lower_index) +=
-            Y_nt_thisupperion * s_renorm[level];
+        atomicadd(*gsl_matrix_ptr(rate_matrix_ntcoll_bf, lower_index, lower_index),
+                  -Y_nt_thisupperion * s_renorm[level]);
+        atomicadd(*gsl_matrix_ptr(rate_matrix_ntcoll_bf, upper_groundstate_index, lower_index),
+                  Y_nt_thisupperion * s_renorm[level]);
       }
     }
   }
@@ -899,7 +900,8 @@ void solve_nlte_pops_element(const int element, const int nonemptymgi, const int
   }
 
   // printout("  Adding rates for ion stages:");
-  for (int ion = 0; ion < nions; ion++) {
+  const auto ions = std::ranges::iota_view{0, nions};
+  std::for_each(EXEC_PAR ions.begin(), ions.end(), [&](const auto ion) {
     // const int ionstage = get_ionstage(element, ion);
     // printout(" %d", ionstage);
 
@@ -916,14 +918,14 @@ void solve_nlte_pops_element(const int element, const int nonemptymgi, const int
     nltepop_matrix_add_boundbound(nonemptymgi, element, ion, t_mid, s_renorm, &rate_matrix_rad_bb, &rate_matrix_coll_bb,
                                   &rate_matrix_ntcoll_bb);
 
-    if (ion < nions - 1) {
+    if (ion < (nions - 1)) {
       // this is the slowest component
       nltepop_matrix_add_ionisation(modelgridindex, element, ion, s_renorm, &rate_matrix_rad_bf, &rate_matrix_coll_bf);
       if (NT_ON) {
         nltepop_matrix_add_nt_ionisation(nonemptymgi, element, ion, s_renorm, &rate_matrix_ntcoll_bf);
       }
     }
-  }
+  });
   // printout("\n");
 
   if (individual_process_matrices) {