prepare_radiation_state Subroutine

   subroutine prepare_radiation_state(mesh, flow, context, params, energy, species_Y, step, time, dt, state, &
                                      debug_T_diff, debug_Y_diff)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(inout) :: flow
      type(radiation_context_t), intent(in) :: context
      type(case_params_t), intent(in) :: params
      type(energy_fields_t), intent(in) :: energy
      real(rk), intent(in), optional :: species_Y(:,:)
      integer, intent(in) :: step
      real(rk), intent(in) :: time, dt
      type(radiation_state_t), intent(inout) :: state
      real(rk), intent(out) :: debug_T_diff, debug_Y_diff

      real(rk), allocatable :: selected_local(:,:)
      integer :: s, c

      state%step = step
      state%time = time
      state%dt = dt
      if (params%enable_chemistry_load_balancing) then
         state%temperature = energy%T
      else
         call flow_allgather_owned_scalar(flow, energy%T, state%temperature)
      end if

      ! Populate radiation thermodynamic pressure from the selected source.
      ! Note: in this low-Mach solver fields%p is a projection correction
      ! potential, NOT an absolute thermodynamic pressure.  Both 'background'
      ! and 'system' therefore resolve to params%background_press (the uniform
      ! thermodynamic p0 passed to Cantera and the energy module).
      select case (trim(context%pressure_source))
      case ('background', 'system')
         state%pressure = params%background_press
      case default
         call fatal_error('radiation', 'unsupported radiation_pressure_source: ' // &
                          trim(context%pressure_source))
      end select

      if (context%n_species > 0) then
         if (.not. present(species_Y)) call fatal_error('radiation', 'radiation selected species but species_Y was not passed')
         if (params%enable_chemistry_load_balancing) then
            do s = 1, context%n_species
               do c = 1, mesh%ncells
                  state%Y(s, c) = species_Y(context%species_index(s), c)
               end do
            end do
         else
            allocate(selected_local(context%n_species, mesh%ncells))
            selected_local = zero
            do s = 1, context%n_species
               do c = flow%first_cell, flow%last_cell
                  selected_local(s, c) = species_Y(context%species_index(s), c)
               end do
            end do
            call flow_allgather_owned_matrix(flow, selected_local, state%Y)
            deallocate(selected_local)
         end if
      end if

      ! Generic scalar gathering is intentionally reserved for the next scalar
      ! registry patch.  The state allocation keeps the API stable.
      if (allocated(state%scalars)) state%scalars = zero

      call compute_gather_debug(flow, context, state, debug_T_diff, debug_Y_diff)
   end subroutine prepare_radiation_state