!> Fractional-step projection solver for constant- and variable-density modes.
!!
!! This module implements the core hydrodynamic solver for the 
!! LowMachReact-Hex framework. Constant-density mode uses `params%rho` and
!! enforces \( \nabla \cdot u = 0 \). Guarded variable-density low-Mach mode
!! uses active density `transport%rho` and enforces
!! \( \nabla \cdot u = S \), where \(S\) is stored in
!! `fields%projection_divergence_source` for the projection time level.
!!
!! 1. **Predictor Step**: An intermediate velocity \(\mathbf{u}^*\) is
!!    calculated by advancing the momentum equation explicitly, excluding 
!!     the new pressure gradient.
!!    $$\frac{\mathbf{u}^* - \mathbf{u}^n}{\Delta t} = -(\mathbf{u} \cdot \nabla) \mathbf{u} + \nu \nabla^2 \mathbf{u} - \frac{1}{\rho} \nabla p^n + \mathbf{f}$$
!!
!! 2. **Poisson Solve**: A pressure correction potential \(\phi = p^{n+1} - p^n\)
!!    is found by solving the Poisson equation derived from the active
!!    continuity constraint.
!!    $$\nabla^2 \phi = \frac{\rho}{\Delta t} \nabla \cdot \mathbf{u}^*$$
!!
!! 3. **Corrector Step**: The final velocity \(\mathbf{u}^{n+1}\) and
!!    pressure \(p^{n+1}\) are updated using the potential gradient.
!!    $$\mathbf{u}^{n+1} = \mathbf{u}^* - \frac{\Delta t}{\rho} \nabla \phi$$
!!    $$p^{n+1} = p^n + \phi$$
!!
!! The linear system for the Poisson equation is solved using a 
!! **Preconditioned Conjugate Gradient (PCG)** method with a 
!! diagonal (Jacobi) preconditioner.
module mod_flow_projection
   use mpi_f08
   use mod_precision, only : rk, zero, one, half, tiny_safe, fatal_error
   use mod_input, only : case_params_t, max_species
   use mod_mesh, only : mesh_t
   use mod_mpi_flow, only : flow_mpi_t, flow_exchange_cell_scalar, &
                            flow_exchange_cell_matrix, flow_exchange_face_scalar, &
                            flow_global_dot_owned, flow_global_two_dots_owned, &
                            flow_global_max_owned
   use mod_bc, only : bc_set_t, bc_periodic, bc_neumann, bc_dirichlet, &
                      patch_type_for_face, boundary_velocity, &
                      face_effective_neighbor, boundary_pressure_type, &
                      boundary_pressure, boundary_species, boundary_temperature
   use mod_profiling, only : profiler_start, profiler_stop
   use mod_flow_fields, only : flow_fields_t
   use mod_cantera, only : transport_properties_t, evaluate_cantera_boundary_state
   use mod_linear_solver, only : solver_stats_t, &
                                 projection_work, &
                                 solve_pressure_correction, &
                                 subtract_owned_mean_scalar, &
                                 finalize_flow_projection_workspace, &
                                 ensure_projection_workspace, &
                                 ensure_pressure_operator_cache, &
                                 face_normal_distance, &
                                 outward_normal, &
                                 pressure_cache_has_dirichlet_pressure, &
                                 pressure_cache_has_neumann_outlet, &
                                 pressure_face_inverse_density

   implicit none

   private

   public :: advance_projection_step, compute_flow_diagnostics
   public :: compute_and_update_cfl, update_lowmach_divergence_source_from_density
   public :: correct_projection_to_lowmach_source
   public :: finalize_flow_projection_workspace, face_normal_distance, solver_stats_t

contains

   !> Orchestrates one full fractional-step iteration.
   !!
   !! This routine implements the high-level logic of the projection method. 
   !! It coordinates the explicit momentum update, the elliptic pressure 
   !! solve, and the final divergence-free correction.
   !!
   !! @param mesh The computational mesh.
   !! @param flow MPI decomposition context.
   !! @param bc Boundary condition set.
   !! @param params Case parameters (dt, rho, etc).
   !! @param transport Physical property fields.
   !! @param fields Flow field containers to be updated.
   !! @param stats Diagnostic stats to be populated.
   !> Advances the low-Mach pressure projection system by one fractional step.
   !!
   !! Solves the fractional-step projection flow equations:
   !! 1. Predictor: Evaluates the momentum advection and diffusion terms using
   !!    explicit second-order Adams-Bashforth time integration to compute
   !!    the intermediate velocity field \( u^* \).
   !! 2. Volumetric Flux Predictor: Interpolates \( u^* \) to cell faces to obtain
   !!    face-centered predicted volumetric fluxes.
   !! 3. Neumann Boundary Balancing: Adjusts Neumann outlet fluxes globally to
   !!    enforce strict numerical conservation constraints and prevent Poisson drift.
   !! 4. Pressure Poisson Solve: Resolves the pressure potential \( \phi \) that
   !!    enforces \( \nabla \cdot u = S_u \) where \( S_u \) is the low-Mach divergence source.
   !! 5. Corrector Step: Updates cell velocities, face fluxes, and thermodynamic pressure.
   !!
   !! @param mesh Computational mesh context.
   !! @param flow MPI parallelization data structure.
   !! @param bc Boundary condition patches.
   !! @param params Configuration input parameters.
   !! @param transport Physical and thermodynamic properties (e.g. density).
   !! @param fields Flow field arrays (u, u_star, p, phi, face_flux, mass_flux) updated in-place.
   !! @param stats Solver diagnostics statistics container.
   subroutine advance_projection_step(mesh, flow, bc, params, transport, fields, stats)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(inout) :: flow
      type(bc_set_t), intent(in) :: bc
      type(case_params_t), intent(in) :: params
      type(transport_properties_t), intent(inout) :: transport
      type(flow_fields_t), intent(inout) :: fields
      type(solver_stats_t), intent(inout) :: stats

      integer :: c                                     !< Cell index counter.
      integer, save :: projection_step_counter = 0     !< Persistent timestep counter to schedule diagnostics.
      logical :: do_projection_diagnostics             !< True if diagnostics/VTK output is scheduled this step.
      real(rk) :: source_value                         !< Resolved divergence source for the current cell.

      projection_step_counter = projection_step_counter + 1
      do_projection_diagnostics = (mod(projection_step_counter, params%output_interval) == 0) .or. &
                                  (projection_step_counter == params%nsteps)

      call ensure_projection_workspace(mesh)
      call ensure_pressure_operator_cache(mesh, flow, bc)

      associate( &
         local_vec => projection_work%local_vec, &
         local_scalar => projection_work%local_scalar, &
         rhs_poisson => projection_work%rhs_poisson, &
         local_face_flux => projection_work%local_face_flux, &
         local_vec_star => projection_work%local_vec_star, &
         predicted_face_flux => projection_work%predicted_face_flux)

         ! 1. Predictor: Compute Explicit momentum RHS and Advance intermediate velocity u*
      ! Constant-density path has S = 0. Variable-density mode preserves the
      ! source prepared after the previous thermo-density sync.
      if (.not. params%enable_variable_density) then
         if (allocated(fields%divergence_source)) fields%divergence_source = zero
      end if

      ! Preserve the source time level actually consumed by this projection.
      ! Later energy/thermo updates may advance fields%divergence_source before
      ! diagnostics are written, so projection residuals compare against
      ! fields%projection_divergence_source.
      if (allocated(fields%projection_divergence_source)) then
         fields%projection_divergence_source = zero
         if (allocated(fields%divergence_source)) then
            fields%projection_divergence_source = fields%divergence_source
         end if
      end if

         call profiler_start('Projection_Momentum_RHS')
         call compute_momentum_rhs(mesh, flow, bc, params, transport, fields%u, fields%p, local_vec)
         call profiler_stop('Projection_Momentum_RHS')
         
         ! 2. Advance intermediate velocity u* locally
         ! This uses the PREVIOUS fields%rhs_old and the CURRENT local_vec
         call profiler_start('Projection_AB2')
         call advance_ab2(mesh, flow, params, fields, local_vec, local_vec_star)
         call profiler_stop('Projection_AB2')
         
         ! Store current RHS for next step (AB2)
         ! rhs_old is effectively partitioned; each rank only needs its owned cell values.
         fields%rhs_old = local_vec
         fields%rhs_old_valid = .true.

         fields%u_star = local_vec_star

         call profiler_start('Projection_Predict_Flux')

         call profiler_start('PredictFlux_Exchange_UStar')
         call flow_exchange_cell_matrix(flow, fields%u_star)
         call profiler_stop('PredictFlux_Exchange_UStar')

         ! 3. Interpolate predicted cell velocity to intermediate face fluxes
         call profiler_start('PredictFlux_Compute')
         call compute_predicted_face_flux(mesh, flow, bc, params, transport, fields%u_star, predicted_face_flux)
         call profiler_stop('PredictFlux_Compute')

         ! 3b. Balance mass at open boundaries to prevent drift in closed-loop systems
         call profiler_start('PredictFlux_Balance')
         call balance_neumann_outlet_flux(mesh, flow, bc, params, fields, predicted_face_flux)
         call profiler_stop('PredictFlux_Balance')

         call profiler_start('PredictFlux_Exchange_Face')
         call flow_exchange_face_scalar(flow, predicted_face_flux)
         call profiler_stop('PredictFlux_Exchange_Face')

         call profiler_stop('Projection_Predict_Flux')

         ! 4. Construct the Poisson RHS for div(u)=S.  Constant-density mode
         ! scales by rho/dt; variable-density mode solves the variable-
         ! coefficient projection in flux form with the source volume integral.
         call profiler_start('Projection_Poisson_RHS')
         call compute_flux_divergence(mesh, flow, predicted_face_flux, local_scalar)
         fields%div = local_scalar

         rhs_poisson = zero
         do c = flow%first_cell, flow%last_cell
            source_value = zero
            if (allocated(fields%divergence_source)) then
               source_value = fields%divergence_source(c)
            end if

            if (params%enable_variable_density) then
               rhs_poisson(c) = -(fields%div(c) - source_value) * &
                                mesh%cells(c)%volume / params%dt
            else
               rhs_poisson(c) = -params%rho / params%dt * &
                                (fields%div(c) - source_value) * mesh%cells(c)%volume
            end if
         end do

         ! Floating pressure handle.
         ! Constant-density mode keeps the historical cell-1 pin. Experimental
         ! variable-density mode keeps every continuity equation active and
         ! removes the pure-Neumann nullspace with a zero-mean RHS gauge.
         if (.not. pressure_cache_has_dirichlet_pressure()) then
            if (params%enable_variable_density) then
               call subtract_owned_mean_scalar(flow, rhs_poisson)
            else
               rhs_poisson(1) = zero
            end if
         end if
         fields%phi = zero
         call profiler_stop('Projection_Poisson_RHS')

         ! Solve Laplacian system for potential phi
         call profiler_start('Projection_PCG')
         call solve_pressure_correction(mesh, flow, bc, params, transport, rhs_poisson, fields%phi, stats)
         call profiler_stop('Projection_PCG')

         ! 5. Update Pressure: p(n+1) = p(n) + phi
         call profiler_start('Projection_Pressure_Update')
         do c = flow%first_cell, flow%last_cell
            fields%p(c) = fields%p(c) + fields%phi(c)
         end do
         if (.not. pressure_cache_has_dirichlet_pressure()) then
            if (params%enable_variable_density) then
               call subtract_owned_mean_scalar(flow, fields%p)
            else
               fields%p(1) = zero
            end if
         end if
         call flow_exchange_cell_scalar(flow, fields%p)
         call profiler_stop('Projection_Pressure_Update')

         ! 6. Correct face fluxes and cell-centered velocity locally
         call profiler_start('Projection_Correction')
         call correct_face_flux(mesh, flow, bc, params, transport, predicted_face_flux, fields%phi, local_face_flux)
         fields%face_flux = local_face_flux
         call flow_exchange_face_scalar(flow, fields%face_flux)

         ! Keep density-weighted face flux available for scalar/enthalpy
         ! transport. In constant-density mode this is exactly
         ! params%rho * face_flux; in variable-density mode rho_f is
         ! interpolated from active density.
         call compute_face_mass_flux(mesh, flow, bc, transport, fields%face_flux, fields%mass_flux)
         call flow_exchange_face_scalar(flow, fields%mass_flux)

         call correct_cell_velocity(mesh, flow, bc, params, transport, fields%u_star, fields%phi, local_vec)

         ! 7. Keep velocity ghosts current for the next predictor/species step.
         fields%u = local_vec
         call flow_exchange_cell_matrix(flow, fields%u)

         ! Divergence is only needed for diagnostics/output, so avoid computing
         ! and exchanging it on non-output steps.
         if (do_projection_diagnostics) then
            call compute_flux_divergence(mesh, flow, fields%face_flux, local_scalar)
            fields%div = local_scalar
            call flow_exchange_cell_scalar(flow, fields%div)
         end if
         call profiler_stop('Projection_Correction')

         if (do_projection_diagnostics) then
            call profiler_start('Projection_Diagnostics')
            call compute_flow_diagnostics(mesh, flow, bc, params, fields, stats)
            call profiler_stop('Projection_Diagnostics')
         end if

      end associate
   end subroutine advance_projection_step


   !> Advances velocity using the 2nd-order Adams-Bashforth scheme.
   !!
   !! Falls back to 1st-order Forward Euler on the very first timestep 
   !! when previous RHS data is unavailable.
   subroutine advance_ab2(mesh, flow, params, fields, rhs, local_ustar)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(in) :: flow
      type(case_params_t), intent(in) :: params
      type(flow_fields_t), intent(in) :: fields
      real(rk), intent(in) :: rhs(:,:)
      real(rk), intent(out) :: local_ustar(:,:)

      integer :: c
      real(rk) :: c_curr, c_old

      local_ustar = zero

      if (fields%rhs_old_valid .and. params%dt_old > tiny_safe) then
         c_curr = 1.0_rk + 0.5_rk * (params%dt / params%dt_old)
         c_old = -0.5_rk * (params%dt / params%dt_old)
         do c = flow%first_cell, flow%last_cell
            local_ustar(:, c) = fields%u(:, c) + params%dt * &
               (c_curr * rhs(:, c) + c_old * fields%rhs_old(:, c))
         end do
      else
         do c = flow%first_cell, flow%last_cell
            local_ustar(:, c) = fields%u(:, c) + params%dt * rhs(:, c)
         end do
      end if

      associate(dummy => mesh%ncells)
      end associate
   end subroutine advance_ab2


   !> Evaluates the advective, diffusive, and pressure terms of the momentum equation.
   !!
   !! Supports both **Upwind** (stable) and **Central** (high-accuracy) 
   !! convection schemes. Advection is computed using a flux-form 
   !! discretization.
   subroutine compute_momentum_rhs(mesh, flow, bc, params, transport, u, p, local_rhs)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(in) :: flow
      type(bc_set_t), intent(in) :: bc
      type(case_params_t), intent(in) :: params
      type(transport_properties_t), intent(in) :: transport
      real(rk), intent(in) :: u(:,:)
      real(rk), intent(in) :: p(:)
      real(rk), intent(out) :: local_rhs(:,:)

      integer :: c, lf, f, nb
      real(rk) :: nvec(3), nhat(3)
      real(rk) :: uf(3), ub(3), advected(3)
      real(rk) :: un_area
      real(rk) :: conv(3), diff(3), gradp(3)
      real(rk) :: dist, diff_face, rho_cell
      logical :: use_central

      use_central = .false.

      select case (trim(params%momentum_convection_scheme))
      case ('central', 'central_difference', 'central-difference')
         use_central = .true.
      case ('upwind', 'first_order_upwind', 'first-order-upwind')
         use_central = .false.
      case default
         call fatal_error('flow', 'unknown momentum_convection_scheme '//trim(params%momentum_convection_scheme))
      end select

      local_rhs = zero

      do c = flow%first_cell, flow%last_cell
         conv = zero
         diff = zero

         ! Compute local pressure gradient at cell center
         call pressure_gradient_cell(mesh, bc, p, c, gradp)

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

            nvec = outward_normal(mesh, f, c)
            nhat = nvec

            nb = face_effective_neighbor(mesh, bc, f, c)

            if (nb > 0) then
               uf = face_linear_vector(mesh, bc, f, c, nb, u(:, c), u(:, nb))
               ub = u(:, nb)
               dist = face_normal_distance(mesh, bc, f, c, nb)
            else
               call boundary_velocity(mesh, bc, f, u(:, c), ub, &
                                      boundary_density=boundary_density_for_cell(params, transport, c))
               uf = ub
               dist = face_normal_distance(mesh, bc, f, c, 0)
            end if

            un_area = dot_product(uf, nhat) * mesh%faces(f)%area

            ! Scheme selection for advection
            if (use_central) then
               advected = uf
            else
               if (un_area >= zero) then
                  advected = u(:, c)
               else
                  advected = ub
               end if
            end if

            conv = conv + un_area * advected

            ! Viscous diffusion Term
            if (nb > 0) then
               diff_face = half * (transport%nu(c) + transport%nu(nb))
            else
               diff_face = transport%nu(c)
            end if

            diff = diff + diff_face * mesh%faces(f)%area * (ub - u(:, c)) / dist
         end do

         ! Assemble total RHS. Constant-density mode uses params%rho; guarded
         ! variable-density mode uses the active cell density.
         rho_cell = params%rho
         if (params%enable_variable_density) rho_cell = max(transport%rho(c), tiny_safe)

         local_rhs(:, c) = -conv / mesh%cells(c)%volume + &
                            diff / mesh%cells(c)%volume + &
                            params%body_force - gradp / rho_cell
      end do
   end subroutine compute_momentum_rhs


   real(rk) function boundary_density_for_cell(params, transport, cell_id) result(rho_b)
      type(case_params_t), intent(in) :: params
      type(transport_properties_t), intent(in) :: transport
      integer, intent(in) :: cell_id

      rho_b = params%rho
      if (params%enable_variable_density) then
         if (allocated(transport%rho)) then
            if (cell_id >= 1 .and. cell_id <= size(transport%rho)) then
               rho_b = max(transport%rho(cell_id), tiny_safe)
            end if
         end if
      end if
   end function boundary_density_for_cell


   !> Calculates the pressure gradient at a cell center using Gauss's Theorem.
   !!
   !! $$\nabla p = \frac{1}{V_c} \int_S p \mathbf{n} dS$$
   subroutine pressure_gradient_cell(mesh, bc, p, cell_id, gradp)
      type(mesh_t), intent(in) :: mesh
      type(bc_set_t), intent(in) :: bc
      real(rk), intent(in) :: p(:)
      integer, intent(in) :: cell_id
      real(rk), intent(out) :: gradp(3)

      integer :: lf, f, nb, patch_id
      real(rk) :: nvec(3)
      real(rk) :: pf, dist
      logical :: is_dirichlet

      gradp = zero

      do lf = 1, mesh%ncell_faces(cell_id)
         f = mesh%cell_faces(lf, cell_id)
         nvec = outward_normal(mesh, f, cell_id)

         nb = face_effective_neighbor(mesh, bc, f, cell_id)

         if (nb > 0) then
            pf = face_linear_scalar(mesh, bc, f, cell_id, nb, p(cell_id), p(nb))
         else
            if (boundary_pressure_type(mesh, bc, f) == bc_dirichlet) then
               call boundary_pressure(mesh, bc, f, p(cell_id), pf, is_dirichlet)
            else
               patch_id = mesh%faces(f)%patch
               if (patch_id > 0) then
                  dist = face_normal_distance(mesh, bc, f, cell_id, 0)
                  pf = p(cell_id) + bc%patches(patch_id)%dpdn * dist
               else
                  pf = p(cell_id)
               end if
            end if
         end if

         gradp = gradp + pf * nvec * mesh%faces(f)%area
      end do

      gradp = gradp / mesh%cells(cell_id)%volume
   end subroutine pressure_gradient_cell


   !> Linearly interpolates cell-centered intermediate velocity to mesh faces.
   subroutine compute_predicted_face_flux(mesh, flow, bc, params, transport, ustar, face_flux)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(in) :: flow
      type(bc_set_t), intent(in) :: bc
      type(case_params_t), intent(in) :: params
      type(transport_properties_t), intent(in) :: transport
      real(rk), intent(in) :: ustar(:,:)
      real(rk), intent(out) :: face_flux(:)

      integer :: i, f, owner, nb
      real(rk) :: uf(3), ub(3)

      face_flux = zero

      do i = 1, size(flow%owned_faces)
         f = flow%owned_faces(i)
         owner = mesh%faces(f)%owner

         nb = face_effective_neighbor(mesh, bc, f, owner)

         if (nb > 0) then
            uf = face_linear_vector(mesh, bc, f, owner, nb, ustar(:, owner), ustar(:, nb))
         else
            call boundary_velocity(mesh, bc, f, ustar(:, owner), ub, &
                                   boundary_density=boundary_density_for_cell(params, transport, owner))
            uf = ub
         end if

         face_flux(f) = dot_product(uf, mesh%faces(f)%normal) * mesh%faces(f)%area
      end do
   end subroutine compute_predicted_face_flux


   !> Corrects face fluxes using the pressure potential gradient.
   subroutine correct_face_flux(mesh, flow, bc, params, transport, predicted_flux, phi, face_flux)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(in) :: flow
      type(bc_set_t), intent(in) :: bc
      type(case_params_t), intent(in) :: params
      type(transport_properties_t), intent(in) :: transport
      real(rk), intent(in) :: predicted_flux(:)
      real(rk), intent(in) :: phi(:)
      real(rk), intent(out) :: face_flux(:)

      integer :: i, f, owner, nb
      real(rk) :: dist, beta_face

      face_flux = zero

      do i = 1, size(flow%owned_faces)
         f = flow%owned_faces(i)
         owner = mesh%faces(f)%owner

         nb = face_effective_neighbor(mesh, bc, f, owner)

         face_flux(f) = predicted_flux(f)

         if (nb > 0) then
            dist = face_normal_distance(mesh, bc, f, owner, nb)

            if (params%enable_variable_density) then
               beta_face = pressure_face_inverse_density(transport, owner, nb)
               face_flux(f) = predicted_flux(f) - params%dt * beta_face * &
                               mesh%faces(f)%area * (phi(nb) - phi(owner)) / dist
            else
               face_flux(f) = predicted_flux(f) - params%dt / params%rho * &
                               mesh%faces(f)%area * (phi(nb) - phi(owner)) / dist
            end if
         else
            if (boundary_pressure_type(mesh, bc, f) == bc_dirichlet) then
               dist = face_normal_distance(mesh, bc, f, owner, 0)
               if (params%enable_variable_density) then
                  beta_face = pressure_face_inverse_density(transport, owner, 0)
                  face_flux(f) = predicted_flux(f) - params%dt * beta_face * &
                                  mesh%faces(f)%area * (zero - phi(owner)) / dist
               else
                  face_flux(f) = predicted_flux(f) - params%dt / params%rho * &
                                  mesh%faces(f)%area * (zero - phi(owner)) / dist
               end if
            end if
         end if
      end do
   end subroutine correct_face_flux


   !> Computes density-weighted face mass flux from volumetric flux.
   !!
   !! This keeps a globally valid face-centered mass flux for scalar and
   !! enthalpy transport.  `face_flux` is volumetric flux in \( \mathrm{m^3/s} \);
   !! `mass_flux` is \( \rho_f \) times that flux in \( \mathrm{kg/s} \).  Both
   !! use owner-to-neighbor orientation.
   !!
   !! In constant-density mode \(\rho_f = \rho_{\text{params}}\). In
   !! variable-density mode \(\rho_f\) is the arithmetic face density from
   !! `transport%rho` for interior faces and the owner density on physical
   !! boundaries.
   subroutine compute_face_mass_flux(mesh, flow, bc, transport, face_flux, mass_flux)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(in) :: flow
      type(bc_set_t), intent(in) :: bc
      type(transport_properties_t), intent(in) :: transport
      real(rk), intent(in) :: face_flux(:)
      real(rk), intent(out) :: mass_flux(:)

      integer :: i, f, owner, nb
      real(rk) :: rho_face

      if (.not. allocated(transport%rho)) then
         call fatal_error('flow', 'transport rho must be allocated before mass-flux update')
      end if

      mass_flux = zero

      do i = 1, size(flow%owned_faces)
         f = flow%owned_faces(i)
         owner = mesh%faces(f)%owner
         nb = face_effective_neighbor(mesh, bc, f, owner)

         if (nb > 0) then
            rho_face = 0.5_rk * (transport%rho(owner) + transport%rho(nb))
         else
            rho_face = transport%rho(owner)
         end if

         mass_flux(f) = rho_face * face_flux(f)
      end do
   end subroutine compute_face_mass_flux


   !> Computes the discrete divergence of face fluxes for each cell.
   subroutine compute_flux_divergence(mesh, flow, face_flux, local_div)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(in) :: flow
      real(rk), intent(in) :: face_flux(:)
      real(rk), intent(out) :: local_div(:)

      integer :: c, lf, f
      real(rk) :: sgn

      local_div = zero

      do c = flow%first_cell, flow%last_cell
         do lf = 1, mesh%ncell_faces(c)
            f = mesh%cell_faces(lf, c)

            if (mesh%faces(f)%owner == c) then
               sgn = one
            else
               sgn = -one
            end if

            local_div(c) = local_div(c) + sgn * face_flux(f) / mesh%cells(c)%volume
         end do
      end do
   end subroutine compute_flux_divergence


   !> Updates cell-centered velocity using the pressure potential gradient.
   subroutine correct_cell_velocity(mesh, flow, bc, params, transport, ustar, phi, local_u)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(in) :: flow
      type(bc_set_t), intent(in) :: bc
      type(case_params_t), intent(in) :: params
      type(transport_properties_t), intent(in) :: transport
      real(rk), intent(in) :: ustar(:,:)
      real(rk), intent(in) :: phi(:)
      real(rk), intent(out) :: local_u(:,:)

      integer :: c, lf, f, nb
      real(rk) :: nvec(3), grad(3), phi_face, rho_cell

      local_u = zero

      do c = flow%first_cell, flow%last_cell
         grad = zero

         do lf = 1, mesh%ncell_faces(c)
            f = mesh%cell_faces(lf, c)
            nvec = outward_normal(mesh, f, c)

            nb = face_effective_neighbor(mesh, bc, f, c)

            if (nb > 0) then
               phi_face = face_linear_scalar(mesh, bc, f, c, nb, phi(c), phi(nb))
            else
               if (boundary_pressure_type(mesh, bc, f) == bc_dirichlet) then
                  phi_face = zero
               else
                  phi_face = phi(c)
               end if
            end if

            grad = grad + phi_face * nvec * mesh%faces(f)%area
         end do

         grad = grad / mesh%cells(c)%volume

         rho_cell = params%rho
         if (params%enable_variable_density) rho_cell = max(transport%rho(c), tiny_safe)

         local_u(:, c) = ustar(:, c) - params%dt / rho_cell * grad
      end do
   end subroutine correct_cell_velocity


   !> Aggregates global raw-divergence, kinetic-energy, and boundary-flux data.
   !!
   !! `stats%max_divergence` and `stats%rms_divergence` are raw
   !! \( \nabla \cdot u \) diagnostics.  They are primary projection metrics for
   !! constant-density runs, but variable-density validation should use
   !! `div(u)-S_projection` diagnostics because the target is \( \nabla \cdot u=S \).
   subroutine compute_flow_diagnostics(mesh, flow, bc, params, fields, stats)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(in) :: flow
      type(bc_set_t), intent(in) :: bc
      type(case_params_t), intent(in) :: params
      type(flow_fields_t), intent(in) :: fields
      type(solver_stats_t), intent(inout) :: stats

      real(rk) :: local_ke, global_ke
      real(rk) :: local_rms, global_rms
      real(rk) :: local_flux, global_flux
      real(rk) :: local_sum_vec(3), global_sum_vec(3)
      integer :: c, ierr

      stats%max_divergence = flow_global_max_owned(flow, fields%div)

      local_rms = zero
      local_ke = zero

      do c = flow%first_cell, flow%last_cell
         local_rms = local_rms + fields%div(c) * fields%div(c)

         local_ke = local_ke + half * params%rho * mesh%cells(c)%volume * &
                    dot_product(fields%u(:, c), fields%u(:, c))
      end do

      call compute_boundary_flux(mesh, flow, bc, fields%face_flux, local_flux)

      local_sum_vec = [local_rms, local_ke, local_flux]
      call MPI_Allreduce(local_sum_vec, global_sum_vec, 3, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      call check_mpi(ierr, 'diagnostic flow sums')

      global_rms = global_sum_vec(1)
      global_ke = global_sum_vec(2)
      global_flux = global_sum_vec(3)

      stats%rms_divergence = sqrt(global_rms / max(real(mesh%ncells, rk), one))
      stats%kinetic_energy = global_ke
      stats%net_boundary_flux = global_flux
   end subroutine compute_flow_diagnostics


   !> Adjusts flux at Neumann outlets to ensure strict global mass balance.
   !!
   !! This is critical for solvers without a pressure-pinned cell, 
   !! as numerical drift in the Poisson RHS can lead to a singular system.
   subroutine balance_neumann_outlet_flux(mesh, flow, bc, params, fields, face_flux)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(in) :: flow
      type(bc_set_t), intent(in) :: bc
      type(case_params_t), intent(in) :: params
      type(flow_fields_t), intent(in) :: fields
      real(rk), intent(inout) :: face_flux(:)

      integer :: c, f, owner, btype, ierr
      real(rk) :: local_net_flux
      real(rk) :: global_net_flux
      real(rk) :: local_target_flux
      real(rk) :: global_target_flux
      real(rk) :: local_outlet_area
      real(rk) :: global_outlet_area
      real(rk) :: correction_per_area
      real(rk) :: local_sum_vec(3), global_sum_vec(3)

      if (pressure_cache_has_dirichlet_pressure()) return

      ! If the case has no Neumann outlet faces, there is nothing to balance.
      ! This avoids two unnecessary MPI_Allreduce calls per timestep for closed
      ! domains such as lid-driven cavity.
      if (.not. pressure_cache_has_neumann_outlet()) return

      local_net_flux = zero
      local_target_flux = zero
      local_outlet_area = zero

      ! Low-Mach compatibility target.  For incompressible/constant-density
      ! flow the target net volume flux is zero.  For variable-density
      ! low-Mach flow, the boundary volume flux must equal integral(S)dV.
      if (params%enable_variable_density .and. allocated(fields%divergence_source)) then
         do c = flow%first_cell, flow%last_cell
            local_target_flux = local_target_flux + fields%divergence_source(c) * mesh%cells(c)%volume
         end do
      end if

      do f = 1, mesh%nfaces
         if (mesh%faces(f)%neighbor /= 0) cycle

         btype = patch_type_for_face(mesh, bc, f)

         if (btype == bc_periodic) cycle

         owner = mesh%faces(f)%owner
         if (.not. flow%owned(owner)) cycle

         local_net_flux = local_net_flux + face_flux(f)

         if (btype == bc_neumann) then
            local_outlet_area = local_outlet_area + mesh%faces(f)%area
         end if
      end do

      local_sum_vec = [local_net_flux, local_outlet_area, local_target_flux]
      call profiler_start('MPI_Communication')
      call MPI_Allreduce(local_sum_vec, global_sum_vec, 3, &
                         MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      call check_mpi(ierr, 'outlet balance sums')
      call profiler_stop('MPI_Communication')

      global_net_flux = global_sum_vec(1)
      global_outlet_area = global_sum_vec(2)
      global_target_flux = global_sum_vec(3)

      if (global_outlet_area <= tiny_safe) then
         if (abs(global_net_flux - global_target_flux) > tiny_safe) then
            call fatal_error('flow', 'boundary flux incompatible with low-Mach target and no neumann outlet faces found')
         end if
         return
      end if

      correction_per_area = (global_target_flux - global_net_flux) / global_outlet_area

      do f = 1, mesh%nfaces
         if (mesh%faces(f)%neighbor /= 0) cycle

         btype = patch_type_for_face(mesh, bc, f)
         if (btype /= bc_neumann) cycle
         owner = mesh%faces(f)%owner
         if (.not. flow%owned(owner)) cycle

         face_flux(f) = face_flux(f) + correction_per_area * mesh%faces(f)%area
      end do
   end subroutine balance_neumann_outlet_flux


   !> Helper to integrated boundary flux on MPI-owned partitions.
   subroutine compute_boundary_flux(mesh, flow, bc, face_flux, local_flux)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(in) :: flow
      type(bc_set_t), intent(in) :: bc
      real(rk), intent(in) :: face_flux(:)
      real(rk), intent(out) :: local_flux

      integer :: f, owner, btype

      local_flux = zero

      do f = 1, mesh%nfaces
         if (mesh%faces(f)%neighbor /= 0) cycle

         btype = patch_type_for_face(mesh, bc, f)

         if (btype == bc_periodic) cycle

         owner = mesh%faces(f)%owner

         if (.not. flow%owned(owner)) cycle

         local_flux = local_flux + face_flux(f)
      end do
   end subroutine compute_boundary_flux


   !> Allocates and resets temporary solver vectors.



   !> Calculates the normal distance between cell centers or cell-to-face.
   !!
   !! Handles periodic boundary logic by accounting for the face pair offset.
   !!
   !! @param mesh The mesh.
   !! @param bc Boundary conditions.
   !! @param face_id Face index.
   !! @param cell_id Source cell index.
   !! @param nb Neighbor cell index (or 0 for boundaries).



   !> Computes the linear interpolation weight for a neighbor cell.
   function face_neighbor_weight(mesh, bc, face_id, cell_id, nb) result(w_nb)
      type(mesh_t), intent(in) :: mesh
      type(bc_set_t), intent(in) :: bc
      integer, intent(in) :: face_id
      integer, intent(in) :: cell_id
      integer, intent(in) :: nb
      real(rk) :: w_nb

      integer :: pair_face
      integer :: btype
      real(rk) :: nvec(3)
      real(rk) :: d_owner
      real(rk) :: d_nb
      real(rk) :: d_total

      if (nb <= 0) then
         w_nb = zero
         return
      end if

      nvec = outward_normal(mesh, face_id, cell_id)

      d_owner = abs(dot_product(mesh%faces(face_id)%center - &
                                mesh%cells(cell_id)%center, nvec))

      if (mesh%faces(face_id)%neighbor == 0) then
         btype = patch_type_for_face(mesh, bc, face_id)

         if (btype == bc_periodic) then
            pair_face = mesh%faces(face_id)%periodic_face

            if (pair_face <= 0) then
               call fatal_error('flow', 'periodic face has no paired face')
            end if

            d_nb = abs(dot_product(mesh%cells(nb)%center - &
                                   mesh%faces(pair_face)%center, nvec))
         else
            d_nb = abs(dot_product(mesh%cells(nb)%center - &
                                   mesh%faces(face_id)%center, nvec))
         end if
      else
         d_nb = abs(dot_product(mesh%cells(nb)%center - &
                                mesh%faces(face_id)%center, nvec))
      end if

      d_total = max(d_owner + d_nb, tiny_safe)

      w_nb = d_owner / d_total
   end function face_neighbor_weight


   !> Linearly interpolates a scalar field to a face.
   function face_linear_scalar(mesh, bc, face_id, cell_id, nb, owner_value, neighbor_value) result(face_value)
      type(mesh_t), intent(in) :: mesh
      type(bc_set_t), intent(in) :: bc
      integer, intent(in) :: face_id
      integer, intent(in) :: cell_id
      integer, intent(in) :: nb
      real(rk), intent(in) :: owner_value, neighbor_value
      real(rk) :: face_value

      real(rk) :: w_nb

      w_nb = face_neighbor_weight(mesh, bc, face_id, cell_id, nb)

      face_value = (one - w_nb) * owner_value + w_nb * neighbor_value
   end function face_linear_scalar


   !> Linearly interpolates a vector field to a face.
   function face_linear_vector(mesh, bc, face_id, cell_id, nb, owner_value, neighbor_value) result(face_value)
      type(mesh_t), intent(in) :: mesh
      type(bc_set_t), intent(in) :: bc
      integer, intent(in) :: face_id
      integer, intent(in) :: cell_id
      integer, intent(in) :: nb
      real(rk), intent(in) :: owner_value(3), neighbor_value(3)
      real(rk) :: face_value(3)

      real(rk) :: w_nb

      w_nb = face_neighbor_weight(mesh, bc, face_id, cell_id, nb)

      face_value = (one - w_nb) * owner_value + w_nb * neighbor_value
   end function face_linear_vector





   !> Internal helper for MPI error checking.
   subroutine check_mpi(ierr, where)
      integer, intent(in) :: ierr
      character(len=*), intent(in) :: where

      if (ierr /= MPI_SUCCESS) call fatal_error('flow', trim(where)//' MPI failure')
   end subroutine check_mpi



   !> Update the guarded variable-density low-Mach divergence source.
   !!
   !! Constant-density mode returns without changing the zero source. In
   !! variable-density mode the next projection targets \( \nabla \cdot u=S \)
   !! with the conservative source
   !! $$
   !! S =
   !! \frac{\rho^{old} - \rho}{\rho \Delta t}
   !! -
   !! \frac{u \cdot \nabla \rho}{\rho}.
   !! $$
   !! The advective term is evaluated from corrected volumetric face fluxes as
   !! a cell-local outward sum,
   !! \( u\cdot\nabla\rho \approx V_c^{-1}\sum_f F_{f,out}(\rho_f-\rho_c) \).
   !! Interior \(\rho_f\) values are centered; physical boundary faces use a
   !! zero-normal-density-gradient assumption because no explicit density
   !! boundary state exists yet. The routine updates `fields%divergence_source`
   !! for the next projection and advances `transport%rho_old` after doing so.
   !!
   !! @param mesh Computational mesh context.
   !! @param flow MPI parallelization communicator data.
   !! @param bc Boundary condition set.
   !! @param params Global simulation configuration parameters.
   !! @param transport Thermodynamic property array containing density.
   !! @param fields Flow fields where the divergence source will be populated.
   !! @param species_Y Optional species mass fraction array for boundary density resolution.
   !! @param T_state Optional temperature field for boundary density resolution.
   subroutine update_lowmach_divergence_source_from_density(mesh, flow, bc, params, transport, fields, species_Y, T_state)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(inout) :: flow
      type(bc_set_t), intent(in) :: bc
      type(case_params_t), intent(in) :: params
      type(transport_properties_t), intent(inout) :: transport
      type(flow_fields_t), intent(inout) :: fields
      real(rk), intent(in), optional :: species_Y(:,:)
      real(rk), intent(in), optional :: T_state(:)

      integer :: c                                     !< Cell index counter.
      integer :: lf                                    !< Local cell-to-face loop counter.
      integer :: fid                                   !< Global face index.
      integer :: nb                                    !< Neighbor cell index across the face.
      integer :: first_cell                            !< Starting cell index for the MPI partition.
      integer :: last_cell                             !< Ending cell index for the MPI partition.
      real(rk) :: vol                                  !< Cell volume [m^3].
      real(rk) :: rho_c                                !< Active cell density at current time level [kg/m^3].
      real(rk) :: rho_old_c                            !< Cell density at the previous time level [kg/m^3].
      real(rk) :: rho_nb                               !< Neighbor cell density [kg/m^3].
      real(rk) :: rho_f                                !< Face density resolved by upwind/interpolation [kg/m^3].
      real(rk) :: flux_out                             !< Outward volumetric face flux [m^3/s].
      real(rk) :: history_source                       !< Temporal history source contribution (drho/dt) / rho.
      real(rk) :: advective_density_term               !< Sum of outward convective density flux [kg/s].
      real(rk) :: advective_source                     !< Advective source contribution -(u . grad(rho)) / rho.
      real(rk) :: rho_floor                            !< Small positive density safety floor to prevent division-by-zero.
      logical :: have_face_flux                        !< Flag indicating whether the face volumetric flux field is allocated.

      ! Local boundary-evaluation variables
      real(rk) :: T_b                                  !< Interpolated boundary temperature [K].
      real(rk) :: Y_b(max_species)                     !< Interpolated boundary species mass fraction vector.
      real(rk) :: h_sens_b                             !< Boundary enthalpy placeholder.
      real(rk) :: rho_b                                !< Resolved thermodynamic density at boundary [kg/m^3].
      real(rk) :: cp_b                                 !< Boundary specific heat placeholder.
      real(rk) :: lambda_b                             !< Boundary conductivity placeholder.
      real(rk) :: mu_b                                 !< Boundary dynamic viscosity placeholder.
      logical :: is_dirichlet_T                        !< True if temperature boundary is Dirichlet.
      logical :: is_dirichlet_Y                        !< True if species boundaries are Dirichlet.
      integer :: k                                     !< Species index loop counter.

      if (.not. params%enable_variable_density) return
      if (.not. allocated(transport%rho)) return

      if (.not. allocated(fields%divergence_source)) then
         allocate(fields%divergence_source(mesh%ncells))
         fields%divergence_source = 0.0_rk
      else if (size(fields%divergence_source) /= mesh%ncells) then
         deallocate(fields%divergence_source)
         allocate(fields%divergence_source(mesh%ncells))
         fields%divergence_source = 0.0_rk
      end if

      if (.not. allocated(transport%rho_old)) then
         allocate(transport%rho_old(mesh%ncells))
         transport%rho_old = transport%rho
      else if (size(transport%rho_old) /= mesh%ncells) then
         deallocate(transport%rho_old)
         allocate(transport%rho_old(mesh%ncells))
         transport%rho_old = transport%rho
      end if

      if (size(transport%rho) < mesh%ncells) return

      rho_floor = tiny(1.0_rk)
      have_face_flux = allocated(fields%face_flux)
      if (have_face_flux) have_face_flux = size(fields%face_flux) >= mesh%nfaces

      first_cell = max(1, flow%first_cell)
      last_cell = min(mesh%ncells, flow%last_cell)
      if (last_cell < first_cell) return

      do c = first_cell, last_cell
         vol = mesh%cells(c)%volume
         rho_c = max(transport%rho(c), rho_floor)
         rho_old_c = transport%rho_old(c)

         if (params%dt > rho_floor) then
            history_source = (rho_old_c - rho_c) / (rho_c * params%dt)
         else
            history_source = 0.0_rk
         end if

         advective_density_term = 0.0_rk
         if (have_face_flux .and. vol > rho_floor) then
            do lf = 1, mesh%ncell_faces(c)
               fid = mesh%cell_faces(lf, c)

               if (mesh%faces(fid)%owner == c) then
                  flux_out = fields%face_flux(fid)
                  nb = mesh%faces(fid)%neighbor
                  if (nb <= 0) nb = mesh%faces(fid)%periodic_neighbor
               else
                  flux_out = -fields%face_flux(fid)
                  nb = mesh%faces(fid)%owner
               end if

               if (nb > 0 .and. nb <= mesh%ncells) then
                  rho_nb = max(transport%rho(nb), rho_floor)
                  if (params%upwind_dilatation_density) then
                     if (flux_out >= 0.0_rk) then
                        rho_f = rho_c
                     else
                        rho_f = rho_nb
                     end if
                  else
                     rho_f = 0.5_rk * (rho_c + rho_nb)
                  end if
               else
                  ! Physical boundary face. If variable density and Cantera are enabled,
                  ! and species/energy are present, evaluate the boundary density using Cantera.
                  rho_f = rho_c
                  if (params%enable_variable_density .and. present(species_Y) .and. present(T_state)) then
                     if (allocated(transport%rho) .and. size(T_state) >= c) then
                        ! 1. Retrieve the boundary temperature
                        call boundary_temperature(mesh, bc, fid, T_state(c), T_b, is_dirichlet_T)
                        
                        ! 2. Retrieve the boundary species mass fractions
                        Y_b = zero
                        do k = 1, params%nspecies
                           call boundary_species(mesh, bc, fid, k, species_Y(k, c), Y_b(k), is_dirichlet_Y)
                        end do
                        
                        ! 3. Evaluate the boundary state density using Cantera
                        call evaluate_cantera_boundary_state(params, T_b, Y_b(1:params%nspecies), &
                                                             h_sens_b, rho_b, cp_b, lambda_b, mu_b)
                        
                        if (params%upwind_dilatation_density) then
                           if (flux_out >= 0.0_rk) then
                              rho_f = rho_c
                           else
                              rho_f = rho_b
                           end if
                        else
                           rho_f = 0.5_rk * (rho_c + rho_b)
                        end if
                     end if
                  end if
               end if

               advective_density_term = advective_density_term + flux_out * (rho_f - rho_c)
            end do

            advective_density_term = advective_density_term / vol
         end if

         advective_source = -advective_density_term / rho_c
         fields%divergence_source(c) = history_source + advective_source

         ! Advance density history only after storing the source that will be
         ! used by the next projection.
         transport%rho_old(c) = transport%rho(c)
      end do

      call flow_exchange_cell_scalar(flow, fields%divergence_source)
      call flow_exchange_cell_scalar(flow, transport%rho_old)

   end subroutine update_lowmach_divergence_source_from_density



   !> Calculates domain-wide CFL and optionally scales the timestep size.
   !!
   !! For stability in reacting flows, the timestep growth is capped 
   !! to 2% per step when `use_dynamic_dt` is active.
   subroutine compute_and_update_cfl(mesh, flow, params, fields, stats)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(in) :: flow
      type(case_params_t), intent(inout) :: params
      type(flow_fields_t), intent(in) :: fields
      type(solver_stats_t), intent(inout) :: stats

      integer :: c, lf, f
      real(rk) :: local_cfl_rate, max_cfl_rate
      real(rk) :: cell_outward_flux
      integer :: ierr

      local_cfl_rate = zero

      do c = flow%first_cell, flow%last_cell
         cell_outward_flux = zero
         do lf = 1, mesh%ncell_faces(c)
            f = mesh%cell_faces(lf, c)
            cell_outward_flux = cell_outward_flux + abs(fields%face_flux(f))
         end do
         cell_outward_flux = half * cell_outward_flux / mesh%cells(c)%volume
         if (cell_outward_flux > local_cfl_rate) then
            local_cfl_rate = cell_outward_flux
         end if
      end do

      call MPI_Allreduce(local_cfl_rate, max_cfl_rate, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr)
      call check_mpi(ierr, 'cfl max rate')

      if (params%use_dynamic_dt) then
         if (max_cfl_rate > tiny_safe) then
            params%dt = min(params%max_cfl / max_cfl_rate, &
                            params%dt * params%dt_growth_limit, &
                            params%max_dt)
         else
            params%dt = min(params%dt * params%dt_growth_limit, params%max_dt)
         end if
         params%dt = max(params%min_dt, min(params%dt, params%max_dt))
      end if

      stats%cfl = max_cfl_rate * params%dt
   end subroutine compute_and_update_cfl


   !> Performs a post-chemistry corrective pressure projection step.
   !!
   !! Re-projects the velocity field onto the divergence manifold after stiff chemical
   !! integration has altered the local density, and thus the divergence source terms.
   !! This step resolves the new pressure correction potential \( \phi \) via Poisson
   !! solver and updates cell velocities, face volumetric fluxes, and mass fluxes to
   !! guarantee strict low-Mach mass conservation.
   !!
   !! @param mesh Computational mesh context.
   !! @param flow MPI parallelization data structure.
   !! @param bc Boundary condition patches.
   !! @param params Configuration input parameters.
   !! @param transport Physical and thermodynamic properties (e.g. density).
   !! @param fields Flow field arrays corrected in-place.
   !! @param stats Solver diagnostics statistics container.
   subroutine correct_projection_to_lowmach_source(mesh, flow, bc, params, transport, fields, stats)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(inout) :: flow
      type(bc_set_t), intent(in) :: bc
      type(case_params_t), intent(in) :: params
      type(transport_properties_t), intent(inout) :: transport
      type(flow_fields_t), intent(inout) :: fields
      type(solver_stats_t), intent(inout) :: stats

      integer :: c                                     !< Cell index loop counter.
      real(rk) :: source_value                         !< Local divergence source term.

      if (.not. params%enable_density_corrector_projection) return
      if (.not. params%enable_variable_density) return

      call ensure_projection_workspace(mesh)
      call ensure_pressure_operator_cache(mesh, flow, bc)

      associate(rhs_poisson => projection_work%rhs_poisson, &
                local_scalar => projection_work%local_scalar, &
                local_face_flux => projection_work%local_face_flux, &
                local_vec => projection_work%local_vec)

         ! Compute current face-flux divergence using the live face_flux.
         call compute_flux_divergence(mesh, flow, fields%face_flux, local_scalar)

         ! Build Poisson RHS for the residual div(u) - S.
         rhs_poisson = zero
         do c = flow%first_cell, flow%last_cell
            source_value = zero
            if (allocated(fields%divergence_source)) then
               source_value = fields%divergence_source(c)
            end if

            if (params%enable_variable_density) then
               rhs_poisson(c) = -(local_scalar(c) - source_value) * &
                                mesh%cells(c)%volume / params%dt
            else
               rhs_poisson(c) = -params%rho / params%dt * &
                                (local_scalar(c) - source_value) * mesh%cells(c)%volume
            end if
         end do

         ! Remove pure-Neumann nullspace (zero-mean gauge for VD mode).
         if (.not. pressure_cache_has_dirichlet_pressure()) then
            if (params%enable_variable_density) then
               call subtract_owned_mean_scalar(flow, rhs_poisson)
            else
               rhs_poisson(1) = zero
            end if
         end if

         fields%phi = zero
         call solve_pressure_correction(mesh, flow, bc, params, transport, rhs_poisson, fields%phi, stats)

         ! Pressure update.
         do c = flow%first_cell, flow%last_cell
            fields%p(c) = fields%p(c) + fields%phi(c)
         end do
         if (.not. pressure_cache_has_dirichlet_pressure()) then
            if (params%enable_variable_density) then
               call subtract_owned_mean_scalar(flow, fields%p)
            else
               fields%p(1) = zero
            end if
         end if
         call flow_exchange_cell_scalar(flow, fields%p)

         ! Correct face fluxes.
         call correct_face_flux(mesh, flow, bc, params, transport, fields%face_flux, fields%phi, local_face_flux)
         fields%face_flux = local_face_flux
         call flow_exchange_face_scalar(flow, fields%face_flux)

         call compute_face_mass_flux(mesh, flow, bc, transport, fields%face_flux, fields%mass_flux)
         call flow_exchange_face_scalar(flow, fields%mass_flux)

         ! Correct cell-centred velocity.
         call correct_cell_velocity(mesh, flow, bc, params, transport, fields%u, fields%phi, local_vec)
         fields%u = local_vec
         call flow_exchange_cell_matrix(flow, fields%u)

      end associate
   end subroutine correct_projection_to_lowmach_source


end module mod_flow_projection
