write_radiation_diagnostics Subroutine

   subroutine write_radiation_diagnostics(mesh, rad, context, state, source, step, time, dt, &
                                          gather_time, compute_time, reduce_time, debug_T_diff, debug_Y_diff)
      type(mesh_t), intent(in) :: mesh
      type(radiation_mpi_t), intent(in) :: rad
      type(radiation_context_t), intent(inout) :: context
      type(radiation_state_t), intent(in) :: state
      type(radiation_source_t), intent(in) :: source
      integer, intent(in) :: step
      real(rk), intent(in) :: time, dt, gather_time, compute_time, reduce_time, debug_T_diff, debug_Y_diff
      integer :: unit_id, s
      real(rk) :: vol, T_avg, P_avg, q_avg, q_sum
      real(rk) :: ymin, ymax, yavg

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

      if (.not. context%diagnostics_initialized) then
         open(newunit=unit_id, file=trim(context%diagnostics_file), status='replace', action='write')
         write(unit_id,'(a)', advance='no') 'step,time,dt,rad_size,n_wavenumbers,wn_first,wn_last,ncells,T_min,T_max,T_avg,P_min,P_max,P_avg,qrad_min,qrad_max,qrad_avg,qrad_sum,gather_time,compute_time,reduce_time,debug_T_sum_diff,debug_Y_sum_maxdiff'
         do s = 1, context%n_species
            write(unit_id,'(a,a,a)', advance='no') ',Y_', trim(context%species_name(s)), '_min'
            write(unit_id,'(a,a,a)', advance='no') ',Y_', trim(context%species_name(s)), '_max'
            write(unit_id,'(a,a,a)', advance='no') ',Y_', trim(context%species_name(s)), '_avg'
         end do
         write(unit_id,*)
         context%diagnostics_initialized = .true.
      else
         open(newunit=unit_id, file=trim(context%diagnostics_file), status='old', position='append', action='write')
      end if

      vol = real(max(1, mesh%ncells), rk)
      T_avg = sum(state%temperature) / vol
      P_avg = sum(state%pressure) / vol
      q_avg = sum(source%qrad) / vol
      q_sum = zero
      do s = 1, mesh%ncells
         q_sum = q_sum + source%qrad(s) * mesh%cells(s)%volume
      end do

      write(unit_id,'(i0,a,es16.8,a,es16.8,a,i0,a,i0,a,i0,a,i0,a,i0,a,es16.8,a,es16.8,a,es16.8,a,es16.8,a,es16.8,a,es16.8,a,es16.8,a,es16.8,a,es16.8,a,es16.8,a,es16.8,a,es16.8,a,es16.8,a,es16.8,a,es16.8)', advance='no') &
         step, ',', time, ',', dt, ',', rad%nprocs, ',', context%n_wavenumbers, ',', context%wn_first, ',', context%wn_last, ',', mesh%ncells, ',', &
         minval(state%temperature), ',', maxval(state%temperature), ',', T_avg, ',', &
         minval(state%pressure), ',', maxval(state%pressure), ',', P_avg, ',', &
         minval(source%qrad), ',', maxval(source%qrad), ',', q_avg, ',', q_sum, ',', &
         gather_time, ',', compute_time, ',', reduce_time, ',', debug_T_diff, ',', debug_Y_diff

      if (allocated(state%Y)) then
         do s = 1, context%n_species
            ymin = minval(state%Y(s, :))
            ymax = maxval(state%Y(s, :))
            yavg = sum(state%Y(s, :)) / vol
            write(unit_id,'(a,es16.8,a,es16.8,a,es16.8)', advance='no') ',', ymin, ',', ymax, ',', yavg
         end do
      end if
      write(unit_id,*)
      close(unit_id)
   end subroutine write_radiation_diagnostics