balance_neumann_outlet_flux Subroutine

   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