repair_cell_mass_fractions Subroutine

   subroutine repair_cell_mass_fractions(Y_cell, nspecies)
      real(rk), intent(inout) :: Y_cell(:) !< Current mass fractions \(Y_k\)
      integer, intent(in) :: nspecies !< Number of species \(N_s\)

      integer :: k !< Species index loop counter
      integer :: k_best !< Index of the species with the highest mass fraction
      integer :: iter !< Loop iteration counter for convergence safety
      real(rk) :: sum_Y !< Sum of mass fractions
      real(rk) :: defect !< The discrepancy from unity constraint: 1 - sum(Y)
      real(rk) :: amount !< The incremental correction applied to the dominant species
      real(rk) :: best_value !< The maximum species mass fraction found
      real(rk), parameter :: simplex_tol = 1.0e-12_rk !< Numerical tolerance for constraint convergence


      if (nspecies <= 0) return

      do k = 1, nspecies
         Y_cell(k) = max(zero, min(one, Y_cell(k)))
      end do

      sum_Y = sum(Y_cell(1:nspecies))
      if (sum_Y <= tiny_safe) then
         Y_cell(1:nspecies) = zero
         Y_cell(1) = one
         return
      end if

      defect = one - sum_Y
      iter = 0
      do while (abs(defect) > simplex_tol .and. iter < 2*nspecies)
         iter = iter + 1
         k_best = 0
         best_value = -one
         
         ! Find the dominant species in this cell (highest mass fraction)
         ! to absorb both surplus and deficit, preserving trace species.
         do k = 1, nspecies
            if (Y_cell(k) > best_value) then
               best_value = Y_cell(k)
               k_best = k
            end if
         end do
         
         if (k_best <= 0) exit
         
         if (defect > zero) then
            amount = min(defect, one - Y_cell(k_best))
            Y_cell(k_best) = Y_cell(k_best) + amount
         else
            amount = min(-defect, Y_cell(k_best))
            Y_cell(k_best) = Y_cell(k_best) - amount
         end if
         
         defect = one - sum(Y_cell(1:nspecies))
      end do

      if (abs(defect) > 10.0_rk*simplex_tol) then
         sum_Y = sum(Y_cell(1:nspecies))
         if (sum_Y > tiny_safe) then
            Y_cell(1:nspecies) = Y_cell(1:nspecies) / sum_Y
         else
            Y_cell(1:nspecies) = zero
            Y_cell(1) = one
         end if
      end if
   end subroutine repair_cell_mass_fractions