initialize_species Subroutine

   subroutine initialize_species(mesh, params, species)
      type(mesh_t), intent(in) :: mesh
      type(case_params_t), intent(in) :: params
      type(species_fields_t), intent(inout) :: species

      integer :: c, k, j
      real(rk) :: sum_Y
      real(rk) :: init_mixture(params%nspecies)
      character(len=name_len) :: target_name

      call finalize_species(species)

      species%nspecies = params%nspecies
      if (species%nspecies <= 0) return

      allocate(species%Y(species%nspecies, mesh%ncells))
      allocate(species%Y_old(species%nspecies, mesh%ncells))
      allocate(species%rhs_old(species%nspecies, mesh%ncells))
      allocate(species%names(species%nspecies))

      species%names = params%species_name(1:species%nspecies)

      ! Prefer named initial composition strings because Cantera may replace
      ! the input species ordering with the full mechanism species ordering.
      init_mixture = zero
      if (composition_string_is_set(params%initial_composition)) then
         call parse_named_composition(params%initial_composition, species%names, species%nspecies, &
                                      init_mixture, 'species initial_composition', &
                                      normalize=.true., require_sum_positive=.true.)
      else
         ! Legacy path: match namelist initial_Y entries by the species_name list
         ! supplied in &species_input.  This remains safer than assuming numeric
         ! indices are unchanged after Cantera species discovery.
         do j = 1, params%namelist_nspecies
            target_name = trim(lowercase(params%namelist_species_name(j)))
            if (len_trim(target_name) == 0) cycle
            do k = 1, species%nspecies
               if (trim(lowercase(species%names(k))) == target_name) then
                  init_mixture(k) = params%initial_Y(j)
                  exit
               end if
            end do
         end do
      end if

      ! Normalize the initial mixture vector
      sum_Y = sum(init_mixture)
      if (sum_Y > zero) then
         init_mixture = init_mixture / sum_Y
      else
         ! Fallback: If no IC specified, set the mixture to 100% of the first species.
         if (species%nspecies > 0) init_mixture(1) = one
      end if

      do c = 1, mesh%ncells
         species%Y(:, c) = init_mixture
      end do

      species%Y_old = species%Y
      species%rhs_old = zero
      species%rhs_old_valid = .false.
   end subroutine initialize_species