radiation_compute_dom Subroutine

   subroutine radiation_compute_dom(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 :: s_x(:), s_y(:), s_z(:), w(:)
      real(rk), allocatable :: diag(:)
      real(rk), allocatable :: coeff(:,:)
      integer, allocatable :: nb(:,:)
      real(rk), allocatable :: rhs(:)
      real(rk), allocatable :: I_m(:)
      real(rk), allocatable :: G_local(:)

      integer :: n_dirs, m, c, lf, f, neighbor, max_faces, max_iter
      real(rk) :: absorption, emissivity_input, I_b_val
      real(rk) :: u, area, T_w, rhs_bc
      real(rk) :: n_out(3)
      logical :: is_dir
      logical, save :: dom_bc_warned = .false.
      logical, save :: dom_neg_intensity_warned = .false.
      logical, save :: dom_emissivity_warned = .false.

      source%qrad = zero

      if (context%rad_size <= 0) then
         call fatal_error('dom', 'Invalid radiation angular communicator size')
      end if

      absorption = max(params%radiation_absorption_coeff, 1.0e-6_rk)
      emissivity_input = max(min(params%radiation_emissivity, 1.0_rk), 1.0e-6_rk)
      if (emissivity_input < 0.999999_rk .and. .not. dom_emissivity_warned) then
         dom_emissivity_warned = .true.
         write(*,'(A)') ' [WARNING] mod_radiation_dom: DOM currently assumes black diffuse boundaries; ' // &
                           'radiation_emissivity is ignored.'
      end if

      call get_quadrature(params%radiation_dom_quadrature, n_dirs, s_x, s_y, s_z, w)

      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(I_m(mesh%ncells))
      allocate(G_local(mesh%ncells))

      G_local = zero

      ! Solve each locally assigned angular direction.
      do m = 1, n_dirs
         if (mod(m - 1, context%rad_size) /= context%rad_rank) cycle

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

         ! Assemble direction-m upwind finite-volume transport system.
         do c = 1, mesh%ncells
            diag(c) = diag(c) + absorption * mesh%cells(c)%volume
            I_b_val = (sigma / pi) * (state%temperature(c)**4)
            rhs(c) = rhs(c) + absorption * I_b_val * mesh%cells(c)%volume

            do lf = 1, mesh%ncell_faces(c)
               f = mesh%cell_faces(lf, c)
               area = mesh%faces(f)%area

               if (mesh%faces(f)%owner == c) then
                  n_out = mesh%faces(f)%normal
               else
                  n_out = -mesh%faces(f)%normal
               end if

               u = s_x(m)*n_out(1) + s_y(m)*n_out(2) + s_z(m)*n_out(3)
               neighbor = face_effective_neighbor(mesh, bc, f, c)

               if (neighbor > 0) then
                  if (u >= zero) then
                     ! Outflow face: local cell contributes to diagonal.
                     diag(c) = diag(c) + u * area
                  else
                     ! Inflow face: upwind value is neighbor intensity.
                     nb(lf, c) = neighbor
                     coeff(lf, c) = -u * area
                  end if
               else
                  if (patch_type_for_face(mesh, bc, f) == bc_symmetry) cycle

                  if (u >= zero) then
                     diag(c) = diag(c) + u * area
                  else
                     ! Black diffuse wall/open-boundary inflow.
                     call boundary_temperature(mesh, bc, f, state%temperature(c), T_w, is_dir)
                     if (.not. is_dir .and. .not. dom_bc_warned) then
                        dom_bc_warned = .true.
                        write(*,'(A)') ' [WARNING] mod_radiation_dom: Boundary temperature missing; ' // &
                                          'falling back to cell temperature.'
                     end if

                     rhs_bc = -u * area * (sigma / pi) * (T_w**4)
                     rhs(c) = rhs(c) + rhs_bc
                  end if
               end if
            end do
         end do

         do c = 1, mesh%ncells
            I_m(c) = (sigma / pi) * (state%temperature(c)**4)
         end do

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

         if (minval(I_m) < -1.0e-10_rk .and. .not. dom_neg_intensity_warned) then
            dom_neg_intensity_warned = .true.
            write(*,'(A,ES12.5)') ' [WARNING] mod_radiation_dom: Negative DOM ordinate intensity detected: ', minval(I_m)
         end if

         G_local = G_local + w(m) * I_m
      end do

      ! Return only local kappa*G contribution. Driver reduces over angular ranks and
      ! subtracts 4*kappa*sigma*T^4.
      do c = 1, mesh%ncells
         source%qrad(c) = absorption * G_local(c)
      end do

      deallocate(diag, coeff, nb, rhs, I_m, G_local)
      deallocate(s_x, s_y, s_z, w)
   end subroutine radiation_compute_dom