radiation_compute_p1 Subroutine

   subroutine radiation_compute_p1(mesh, bc, context, state, source, params)
      type(mesh_t), intent(in) :: mesh
      type(bc_set_t), intent(in) :: bc
      type(radiation_context_t), intent(in) :: context
      type(radiation_state_t), intent(in) :: state
      type(radiation_source_t), intent(inout) :: source
      type(case_params_t), intent(in) :: params

      real(rk), allocatable :: diag(:)
      real(rk), allocatable :: coeff(:,:)
      integer, allocatable :: nb(:,:)
      real(rk), allocatable :: rhs(:)
      real(rk), allocatable :: G(:)

      integer :: c, lf, f, neighbor, max_faces, max_iter
      real(rk) :: absorption, gamma, emissivity
      real(rk) :: dist, coeff_val, denom, T_w, rhs_bc
      logical :: is_dir
      logical, save :: p1_bc_warned = .false.
      logical, save :: p1_neg_g_warned = .false.

      ! Initialize the output partial source term.
      source%qrad = zero

      ! In the parallel spectral framework, only the rank assigned to wavenumber 1
      ! computes this single grey band.
      if (context%wn_first > 1 .or. context%wn_last < 1) return

      max_faces = size(mesh%cell_faces, 1)

      allocate(diag(mesh%ncells))
      allocate(coeff(max_faces, mesh%ncells))
      allocate(nb(max_faces, mesh%ncells))
      allocate(rhs(mesh%ncells))
      allocate(G(mesh%ncells))

      diag = zero
      coeff = zero
      nb = 0
      rhs = zero

      absorption = max(params%radiation_absorption_coeff, 1.0e-6_rk)
      emissivity = max(min(params%radiation_emissivity, 1.0_rk), 1.0e-6_rk)
      gamma = 1.0_rk / (3.0_rk * absorption)

      ! Assemble the P-1 finite-volume system.
      do c = 1, mesh%ncells
         diag(c) = diag(c) + absorption * mesh%cells(c)%volume
         rhs(c) = rhs(c) + 4.0_rk * absorption * sigma * &
                            (state%temperature(c)**4) * mesh%cells(c)%volume

         do lf = 1, mesh%ncell_faces(c)
            f = mesh%cell_faces(lf, c)
            neighbor = face_effective_neighbor(mesh, bc, f, c)

            if (neighbor > 0) then
               ! Internal or periodic face.
               dist = face_normal_distance(mesh, bc, f, c, neighbor)
               if (dist <= tiny_safe) then
                  call fatal_error('p1', 'Non-positive face-normal distance in internal P-1 face assembly')
               end if

               coeff_val = gamma * mesh%faces(f)%area / dist
               nb(lf, c) = neighbor
               coeff(lf, c) = coeff_val
               diag(c) = diag(c) + coeff_val
            else
               ! Boundary face.
               if (patch_type_for_face(mesh, bc, f) == bc_symmetry) cycle

               ! Non-periodic wall or open boundary: apply Marshak BC.
               ! The denominator represents boundary resistance:
               !   denom = Delta_x/Gamma + 2*(2 - epsilon)/epsilon
               dist = face_normal_distance(mesh, bc, f, c, 0)
               if (dist <= tiny_safe) then
                  call fatal_error('p1', 'Non-positive face-normal distance in boundary P-1 face assembly')
               end if

               denom = 2.0_rk * (2.0_rk - emissivity) / emissivity + dist / gamma
               if (denom <= tiny_safe) then
                  call fatal_error('p1', 'Non-positive Marshak boundary denominator')
               end if

               coeff_val = mesh%faces(f)%area / denom
               diag(c) = diag(c) + coeff_val

               call boundary_temperature(mesh, bc, f, state%temperature(c), T_w, is_dir)
               if (.not. is_dir .and. .not. p1_bc_warned) then
                  p1_bc_warned = .true.
                  write(*,'(A)') ' [WARNING] mod_radiation_p1: Radiative boundary faces have no prescribed temperature; ' // &
                                    'falling back to cell temperature.'
               end if

               rhs_bc = coeff_val * 4.0_rk * sigma * (T_w**4)
               rhs(c) = rhs(c) + rhs_bc
            end if
         end do
      end do

      ! Initial guess: local blackbody incident radiation.
      do c = 1, mesh%ncells
         G(c) = 4.0_rk * sigma * (state%temperature(c)**4)
      end do

      max_iter = 1000
      call solve_pcg(mesh, diag, coeff, nb, rhs, G, max_iter, 1.0e-8_rk)

      ! Incident radiation should be nonnegative. Warn once if the linear solve produces
      ! an unphysical value, but do not clip silently.
      if (minval(G) < -1.0e-10_rk .and. .not. p1_neg_g_warned) then
         p1_neg_g_warned = .true.
         write(*,'(A,ES12.5)') ' [WARNING] mod_radiation_p1: Negative incident radiation detected: ', minval(G)
      end if

      ! Return only the partial absorption-like contribution. The driver subtracts
      ! 4*kappa*sigma*T^4 after global synchronization.
      do c = 1, mesh%ncells
         source%qrad(c) = absorption * G(c)
      end do

      deallocate(diag, coeff, nb, rhs, G)
   end subroutine radiation_compute_p1