write_radiation_setup Subroutine

   subroutine write_radiation_setup(mesh, rad, params, context)
      type(mesh_t), intent(in) :: mesh
      type(radiation_mpi_t), intent(in) :: rad
      type(case_params_t), intent(in) :: params
      type(radiation_context_t), intent(inout) :: context
      integer :: unit_id, r, first, last, nlocal, s

      if (.not. context%enabled) return
      if (.not. context%write_diagnostics) return
      if (context%rad_rank /= 0) return

      open(newunit=unit_id, file=trim(context%setup_file), status='replace', action='write')
      write(unit_id,'(a)') 'Radiation setup'
      write(unit_id,'(a)') '---------------'
      write(unit_id,'(a,l1)') 'enabled: ', context%enabled
      write(unit_id,'(a,a)') 'model: ', trim(context%model)
      write(unit_id,'(a,i0)') 'radiation communicator size: ', rad%nprocs
      write(unit_id,'(a,i0)') 'n_wavenumbers: ', context%n_wavenumbers
      write(unit_id,'(a,i0)') 'ncells: ', mesh%ncells
      write(unit_id,'(a,i0)') 'nfaces: ', mesh%nfaces
      write(unit_id,'(a,a)') 'pressure_source: ', trim(context%pressure_source)
      write(unit_id,'(a,l1)') 'mesh_cached: ', context%mesh_cached
      write(unit_id,'(a,l1)') 'debug_checksums: ', context%debug
      write(unit_id,'(a)') ''
      write(unit_id,'(a)') 'Wavenumber decomposition:'
      do r = 0, max(0, rad%nprocs-1)
         call local_wavenumber_bounds(context%n_wavenumbers, rad%nprocs, r, first, last, nlocal)
         write(unit_id,'(a,i0,a,i0,a,i0,a,i0)') '  rank ', r, ': ', first, '-', last, ' count=', nlocal
      end do
      write(unit_id,'(a)') ''
      write(unit_id,'(a)') 'Selected radiation species:'
      if (context%n_species == 0) then
         write(unit_id,'(a)') '  <none>'
      else
         do s = 1, context%n_species
            write(unit_id,'(a,a,a,i0)') '  ', trim(context%species_name(s)), ' -> solver species index ', context%species_index(s)
         end do
      end if
      write(unit_id,'(a)') ''
      write(unit_id,'(a)') 'Selected radiation scalars:'
      if (context%n_scalars == 0) then
         write(unit_id,'(a)') '  <none currently gathered; scalar registry pending>'
      else
         do s = 1, context%n_scalars
            write(unit_id,'(a,a)') '  ', trim(context%scalar_name(s))
         end do
      end if
      write(unit_id,'(a)') ''
      write(unit_id,'(a)') 'Developer contract:'
      write(unit_id,'(a)') '  Implement physics kernels under src/radiation/models/mod_radiation_<modelname>.f90.'
      write(unit_id,'(a)') '  Models receive full mesh/state and assigned wavenumber range; do not call MPI in model files.'
      write(unit_id,'(a,es16.8)') 'background_press: ', params%background_press
      close(unit_id)
      context%setup_written = .true.
   end subroutine write_radiation_setup