m_mpi_common.fpp

!>
!! @file
!! @brief Contains module m_mpi_common

#:include 'case.fpp'
#:include 'macros.fpp'

!> @brief MPI communication layer: domain decomposition, halo exchange, reductions, and parallel I/O setup
module m_mpi_common

#ifdef MFC_MPI
    use mpi  !< Message passing interface (MPI) module
#endif

    use m_derived_types
    use m_global_parameters
    use m_helper
    use ieee_arithmetic
    use m_nvtx
    use m_constants, only: recon_type_weno

    implicit none

    private :: s_apply_decomposition_policies

    integer, private :: v_size
    $:GPU_DECLARE(create='[v_size]')

    real(wp), private, allocatable, dimension(:) :: buff_send  !< Primitive variable send buffer for halo exchange
    !> Primitive variable receive buffer for halo exchange Variables for EL bubbles communication
    real(wp), private, allocatable, dimension(:) :: buff_recv
    type(int_bounds_info)                        :: comm_coords(3)
    integer                                      :: comm_size(3)
    !> q_beta indices to communicate: 1=void fraction, 2=d(beta)/dt, 5=energy source
    integer :: beta_vars(1:3) = [1, 2, 5]
    $:GPU_DECLARE(create='[comm_coords, comm_size, beta_vars]')

#ifndef __NVCOMPILER_GPU_UNIFIED_MEM
    $:GPU_DECLARE(create='[buff_send, buff_recv]')
#endif

    integer(kind=8) :: halo_size
    $:GPU_DECLARE(create='[halo_size]')

    logical, private :: exchange_all_chemistry_temperatures = .false.
    logical, private :: use_rdma_transport = .false.

contains

    !> Initialize the module.
    impure subroutine s_initialize_mpi_common_module(exchange_all_chemistry_temperatures_in, use_rdma_transport_in)

        logical, intent(in) :: exchange_all_chemistry_temperatures_in
        logical, intent(in) :: use_rdma_transport_in

        exchange_all_chemistry_temperatures = exchange_all_chemistry_temperatures_in
        use_rdma_transport = use_rdma_transport_in

#ifdef MFC_MPI
        ! Allocating buff_send/recv and. Please note that for the sake of simplicity, both variables are provided sufficient storage
        ! to hold the largest buffer in the computational domain.

        if (qbmm .and. .not. polytropic) then
            v_size = sys_size + 2*nb*nnode
        else if (chemistry .and. (chem_params%diffusion .or. exchange_all_chemistry_temperatures)) then
            v_size = sys_size + 1
        else
            v_size = sys_size
        end if

        if (n > 0) then
            if (p > 0) then
                halo_size = nint(-1._wp + 1._wp*buff_size*(v_size)*(m + 2*buff_size + 1)*(n + 2*buff_size + 1)*(p + 2*buff_size &
                                 & + 1)/(cells_bounds%mnp_min + 2*buff_size + 1))
            else
                halo_size = -1 + buff_size*(v_size)*(cells_bounds%mn_max + 2*buff_size + 1)
            end if
        else
            halo_size = -1 + buff_size*(v_size)
        end if

        $:GPU_UPDATE(device='[halo_size, v_size]')

#ifndef __NVCOMPILER_GPU_UNIFIED_MEM
        @:ALLOCATE(buff_send(0:halo_size), buff_recv(0:halo_size))
#else
        allocate (buff_send(0:halo_size), buff_recv(0:halo_size))
        $:GPU_ENTER_DATA(create='[capture:buff_send]')
        $:GPU_ENTER_DATA(create='[capture:buff_recv]')
#endif
#endif

        $:GPU_UPDATE(device='[beta_vars]')

    end subroutine s_initialize_mpi_common_module

    !> Initialize the MPI execution environment and query the number of processors and local rank.
    impure subroutine s_mpi_initialize

#ifdef MFC_MPI
        integer :: ierr  !< Generic flag used to identify and report MPI errors

        call MPI_INIT(ierr)

        if (ierr /= MPI_SUCCESS) then
            print '(A)', 'Unable to initialize MPI environment. Exiting.'
            call MPI_ABORT(MPI_COMM_WORLD, 1, ierr)
        end if

        call MPI_COMM_SIZE(MPI_COMM_WORLD, num_procs, ierr)

        call MPI_COMM_RANK(MPI_COMM_WORLD, proc_rank, ierr)
#else
        num_procs = 1
        proc_rank = 0
#endif

        $:GPU_UPDATE(device='[num_procs, proc_rank]')

    end subroutine s_mpi_initialize

    !> Set up MPI I/O data views and variable pointers for parallel file output.
    impure subroutine s_initialize_mpi_data(q_cons_vf, ib_markers, ib_mpi_data, beta, qbmm_pb, qbmm_mv)

        type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf
        type(integer_field), optional, intent(in)           :: ib_markers
        type(mpi_io_ib_var), optional, intent(inout)        :: ib_mpi_data
        type(scalar_field), intent(in), optional            :: beta
        type(pres_field), intent(in), optional              :: qbmm_pb, qbmm_mv
        integer, dimension(num_dims)                        :: sizes_glb, sizes_loc

#ifdef MFC_MPI
        integer :: i, j
        integer :: ierr  !< Generic flag used to identify and report MPI errors
        integer :: alt_sys
        logical :: bind_qbmm_fields

        if (present(qbmm_pb) .neqv. present(qbmm_mv)) then
            call s_mpi_abort('QBMM MPI I/O requires both pressure and moment fields.')
        end if
        bind_qbmm_fields = qbmm .and. .not. polytropic .and. present(qbmm_pb) .and. present(qbmm_mv)

        if (present(beta)) then
            alt_sys = sys_size + 1
        else
            alt_sys = sys_size
        end if

        do i = 1, sys_size
            MPI_IO_DATA%var(i)%sf => q_cons_vf(i)%sf(0:m,0:n,0:p)
        end do

        if (present(beta)) then
            MPI_IO_DATA%var(alt_sys)%sf => beta%sf(0:m,0:n,0:p)
        end if

        ! Additional variables pb and mv for non-polytropic qbmm
        if (bind_qbmm_fields) then
            do i = 1, nb
                do j = 1, nnode
                    MPI_IO_DATA%var(sys_size + (i - 1)*nnode + j)%sf => qbmm_pb%sf(0:m,0:n,0:p,j, i)
                    MPI_IO_DATA%var(sys_size + (i - 1)*nnode + j + nb*nnode)%sf => qbmm_mv%sf(0:m,0:n,0:p,j, i)
                end do
            end do
        end if

        ! Define global(g) and local(l) sizes for flow variables
        sizes_glb(1) = m_glb + 1; sizes_loc(1) = m + 1
        if (n > 0) then
            sizes_glb(2) = n_glb + 1; sizes_loc(2) = n + 1
            if (p > 0) then
                sizes_glb(num_dims) = p_glb + 1; sizes_loc(num_dims) = p + 1
            end if
        end if

        ! Define the view for each variable
        do i = 1, alt_sys
            call MPI_TYPE_CREATE_SUBARRAY(num_dims, sizes_glb, sizes_loc, start_idx, MPI_ORDER_FORTRAN, mpi_p, &
                                          & MPI_IO_DATA%view(i), ierr)
            call MPI_TYPE_COMMIT(MPI_IO_DATA%view(i), ierr)
        end do

        if (bind_qbmm_fields) then
            do i = sys_size + 1, sys_size + 2*nb*nnode
                call MPI_TYPE_CREATE_SUBARRAY(num_dims, sizes_glb, sizes_loc, start_idx, MPI_ORDER_FORTRAN, mpi_p, &
                                              & MPI_IO_DATA%view(i), ierr)
                call MPI_TYPE_COMMIT(MPI_IO_DATA%view(i), ierr)
            end do
        end if

        if (present(ib_markers) .neqv. present(ib_mpi_data)) then
            call s_mpi_abort('Immersed-boundary MPI I/O requires both marker and descriptor fields.')
        end if

        if (present(ib_markers)) then
            ib_mpi_data%var%sf => ib_markers%sf(0:m,0:n,0:p)
            call MPI_TYPE_CREATE_SUBARRAY(num_dims, sizes_glb, sizes_loc, start_idx, MPI_ORDER_FORTRAN, MPI_INTEGER, &
                                          & ib_mpi_data%view, ierr)
            call MPI_TYPE_COMMIT(ib_mpi_data%view, ierr)
        end if
#endif

    end subroutine s_initialize_mpi_data

    !> Set up MPI I/O data views for downsampled (coarsened) parallel file output.
    subroutine s_initialize_mpi_data_ds(m_ds, n_ds, p_ds, q_cons_vf)

        integer, intent(in)                                           :: m_ds, n_ds, p_ds
        type(scalar_field), dimension(sys_size), intent(in), optional :: q_cons_vf
        integer, dimension(num_dims)                                  :: sizes_loc
        integer, dimension(3)                                         :: sf_start_idx

#ifdef MFC_MPI
        integer :: i, ierr

        sf_start_idx = (/0, 0, 0/)

        if (present(q_cons_vf)) then
            do i = 1, sys_size
                MPI_IO_DATA%var(i)%sf => q_cons_vf(i)%sf(-1:m_ds + 1,-1:n_ds + 1,-1:p_ds + 1)
            end do
        end if
        ! Define global(g) and local(l) sizes for flow variables
        sizes_loc(1) = m_ds + 3
        if (n > 0) then
            sizes_loc(2) = n_ds + 3
            if (p > 0) then
                sizes_loc(num_dims) = p_ds + 3
            end if
        end if

        ! Define the view for each variable
        do i = 1, sys_size
            call MPI_TYPE_CREATE_SUBARRAY(num_dims, sizes_loc, sizes_loc, sf_start_idx, MPI_ORDER_FORTRAN, mpi_p, &
                                          & MPI_IO_DATA%view(i), ierr)
            call MPI_TYPE_COMMIT(MPI_IO_DATA%view(i), ierr)
        end do
#endif

    end subroutine s_initialize_mpi_data_ds

    !> Gather variable-length real vectors from all MPI ranks onto the root process.
    impure subroutine s_mpi_gather_data(my_vector, counts, gathered_vector, root)

        integer, intent(in)                     :: counts              !< Array of vector lengths for each process
        real(wp), intent(in), dimension(counts) :: my_vector           !< Input vector on each process
        integer, intent(in)                     :: root                !< Rank of the root process
        real(wp), allocatable, intent(out)      :: gathered_vector(:)  !< Gathered vector on the root process
        integer                                 :: i
        integer                                 :: ierr                !< Generic flag used to identify and report MPI errors
        integer, allocatable                    :: recounts(:), displs(:)

#ifdef MFC_MPI
        allocate (recounts(num_procs))

        call MPI_GATHER(counts, 1, MPI_INTEGER, recounts, 1, MPI_INTEGER, root, MPI_COMM_WORLD, ierr)

        allocate (displs(size(recounts)))

        displs(1) = 0

        do i = 2, size(recounts)
            displs(i) = displs(i - 1) + recounts(i - 1)
        end do

        allocate (gathered_vector(sum(recounts)))
        call MPI_GATHERV(my_vector, counts, mpi_p, gathered_vector, recounts, displs, mpi_p, root, MPI_COMM_WORLD, ierr)
#endif

    end subroutine s_mpi_gather_data

    !> Gather per-rank time step wall-clock times onto rank 0 for performance reporting.
    impure subroutine mpi_bcast_time_step_values(proc_time, time_avg)

        real(wp), dimension(0:num_procs - 1), intent(inout) :: proc_time
        real(wp), intent(inout)                             :: time_avg

#ifdef MFC_MPI
        integer :: ierr  !< Generic flag used to identify and report MPI errors

        call MPI_GATHER(time_avg, 1, mpi_p, proc_time(0), 1, mpi_p, 0, MPI_COMM_WORLD, ierr)
#endif

    end subroutine mpi_bcast_time_step_values

    !> Print a case file error with the prohibited condition and message, then abort execution.
    impure subroutine s_prohibit_abort(condition, message)

        character(len=*), intent(in) :: condition, message

        print *, ""
        print *, "CASE FILE ERROR"
        print *, "  - Prohibited condition: ", trim(condition)
        if (len_trim(message) > 0) then
            print *, "  - Note: ", trim(message)
        end if
        print *, ""
        call s_mpi_abort(code=CASE_FILE_ERROR_CODE)

    end subroutine s_prohibit_abort

    !> The goal of this subroutine is to determine the global extrema of the stability criteria in the computational domain. This is
    !! performed by sifting through the local extrema of each stability criterion. Note that each of the local extrema is from a
    !! single process, within its assigned section of the computational domain. Finally, note that the global extrema values are
    !! only bookkeept on the rank 0 processor.
    impure subroutine s_mpi_reduce_stability_criteria_extrema(icfl_max_loc, vcfl_max_loc, Rc_min_loc, bubs_loc, icfl_max_glb, &
        & vcfl_max_glb, Rc_min_glb, bubs_glb, ccfl_max_loc, ccfl_max_glb)

        real(wp), intent(in)  :: icfl_max_loc
        real(wp), intent(in)  :: vcfl_max_loc
        real(wp), intent(in)  :: Rc_min_loc
        integer, intent(in)   :: bubs_loc
        real(wp), intent(out) :: icfl_max_glb
        real(wp), intent(out) :: vcfl_max_glb
        real(wp), intent(out) :: Rc_min_glb
        integer, intent(out)  :: bubs_glb
        real(wp), intent(in)  :: ccfl_max_loc
        real(wp), intent(out) :: ccfl_max_glb

        icfl_max_glb = icfl_max_loc
        vcfl_max_glb = vcfl_max_loc
        Rc_min_glb = Rc_min_loc
        ccfl_max_glb = ccfl_max_loc

#ifdef MFC_MPI
        block
            integer :: ierr

            bubs_glb = 0
            call MPI_REDUCE(icfl_max_loc, icfl_max_glb, 1, mpi_p, MPI_MAX, 0, MPI_COMM_WORLD, ierr)

            if (viscous) then
                call MPI_REDUCE(vcfl_max_loc, vcfl_max_glb, 1, mpi_p, MPI_MAX, 0, MPI_COMM_WORLD, ierr)
                call MPI_REDUCE(Rc_min_loc, Rc_min_glb, 1, mpi_p, MPI_MIN, 0, MPI_COMM_WORLD, ierr)
            end if

            if (surface_tension) then
                call MPI_REDUCE(ccfl_max_loc, ccfl_max_glb, 1, mpi_p, MPI_MAX, 0, MPI_COMM_WORLD, ierr)
            end if

            if (bubbles_lagrange) then
                call MPI_REDUCE(bubs_loc, bubs_glb, 1, MPI_INTEGER, MPI_SUM, 0, MPI_COMM_WORLD, ierr)
            end if
        end block
#else
        icfl_max_glb = icfl_max_loc
        bubs_glb = 0

        if (viscous) then
            vcfl_max_glb = vcfl_max_loc
            Rc_min_glb = Rc_min_loc
        end if

        if (surface_tension) then
            ccfl_max_glb = ccfl_max_loc
        end if

        if (bubbles_lagrange) bubs_glb = bubs_loc
#endif

    end subroutine s_mpi_reduce_stability_criteria_extrema

    !> Reduce a local integer value to its global sum across all MPI ranks.
    subroutine s_mpi_reduce_int_sum(var_loc, sum)

        integer, intent(in)  :: var_loc
        integer, intent(out) :: sum

#ifdef MFC_MPI
        integer :: ierr  !< Generic flag used to identify and report MPI errors

        call MPI_REDUCE(var_loc, sum, 1, MPI_INTEGER, MPI_SUM, 0, MPI_COMM_WORLD, ierr)
#else
        sum = var_loc
#endif

    end subroutine s_mpi_reduce_int_sum

    !> Reduce a local real value to its global sum across all MPI ranks.
    impure subroutine s_mpi_allreduce_sum(var_loc, var_glb)

        real(wp), intent(in)  :: var_loc
        real(wp), intent(out) :: var_glb

#ifdef MFC_MPI
        integer :: ierr  !< Generic flag used to identify and report MPI errors

        call MPI_ALLREDUCE(var_loc, var_glb, 1, mpi_p, MPI_SUM, MPI_COMM_WORLD, ierr)
#endif

    end subroutine s_mpi_allreduce_sum

    !> Reduce an array of vectors to their global sums across all MPI ranks.
    impure subroutine s_mpi_allreduce_vectors_sum(var_loc, var_glb, num_vectors, vector_length)

        integer, intent(in)                     :: num_vectors, vector_length
        real(wp), dimension(:,:), intent(in)    :: var_loc
        real(wp), dimension(:,:), intent(inout) :: var_glb

#ifdef MFC_MPI
        integer :: ierr  !< Generic flag used to identify and report MPI errors

        if (loc(var_loc) == loc(var_glb)) then
            call MPI_Allreduce(MPI_IN_PLACE, var_glb, num_vectors*vector_length, mpi_p, MPI_SUM, MPI_COMM_WORLD, ierr)
        else
            call MPI_Allreduce(var_loc, var_glb, num_vectors*vector_length, mpi_p, MPI_SUM, MPI_COMM_WORLD, ierr)
        end if
#else
        var_glb(1:num_vectors,1:vector_length) = var_loc(1:num_vectors,1:vector_length)
#endif

    end subroutine s_mpi_allreduce_vectors_sum

    !> Reduce a local integer value to its global sum across all MPI ranks.
    impure subroutine s_mpi_allreduce_integer_sum(var_loc, var_glb)

        integer(kind=8), intent(in)  :: var_loc
        integer(kind=8), intent(out) :: var_glb

#ifdef MFC_MPI
        integer :: ierr  !< Generic flag used to identify and report MPI errors

        call MPI_ALLREDUCE(var_loc, var_glb, 1, MPI_INTEGER8, MPI_SUM, MPI_COMM_WORLD, ierr)
#else
        var_glb = var_loc
#endif

    end subroutine s_mpi_allreduce_integer_sum

    !> Reduce a local real value to its global minimum across all MPI ranks.
    impure subroutine s_mpi_allreduce_min(var_loc, var_glb)

        real(wp), intent(in)  :: var_loc
        real(wp), intent(out) :: var_glb

#ifdef MFC_MPI
        integer :: ierr  !< Generic flag used to identify and report MPI errors

        call MPI_ALLREDUCE(var_loc, var_glb, 1, mpi_p, MPI_MIN, MPI_COMM_WORLD, ierr)
#endif

    end subroutine s_mpi_allreduce_min

    !> Reduce a local real vector to its elementwise global minimum across all MPI ranks.
    impure subroutine s_mpi_allreduce_min_vec(var_loc, var_glb)

        real(wp), dimension(:), intent(in)  :: var_loc
        real(wp), dimension(:), intent(out) :: var_glb

#ifdef MFC_MPI
        integer :: ierr  !< Generic flag used to identify and report MPI errors

        call MPI_ALLREDUCE(var_loc, var_glb, size(var_loc), mpi_p, MPI_MIN, MPI_COMM_WORLD, ierr)
#else
        var_glb = var_loc
#endif

    end subroutine s_mpi_allreduce_min_vec

    !> Reduce a local real value to its global maximum across all MPI ranks.
    impure subroutine s_mpi_allreduce_max(var_loc, var_glb)

        real(wp), intent(in)  :: var_loc
        real(wp), intent(out) :: var_glb

#ifdef MFC_MPI
        integer :: ierr  !< Generic flag used to identify and report MPI errors

        call MPI_ALLREDUCE(var_loc, var_glb, 1, mpi_p, MPI_MAX, MPI_COMM_WORLD, ierr)
#endif

    end subroutine s_mpi_allreduce_max

    !> Reduce a local real value to its global minimum across all ranks
    impure subroutine s_mpi_reduce_min(var_loc)

        real(wp), intent(inout) :: var_loc

#ifdef MFC_MPI
        integer  :: ierr  !< Generic flag used to identify and report MPI errors
        real(wp) :: var_glb

        call MPI_REDUCE(var_loc, var_glb, 1, mpi_p, MPI_MIN, 0, MPI_COMM_WORLD, ierr)

        call MPI_BCAST(var_glb, 1, mpi_p, 0, MPI_COMM_WORLD, ierr)

        var_loc = var_glb
#endif

    end subroutine s_mpi_reduce_min

    !> Reduce a 2-element variable to its global maximum value with the owning processor rank (MPI_MAXLOC).
    !> Reduce a local value to its global maximum with location (rank) across all ranks
    impure subroutine s_mpi_reduce_maxloc(var_loc)

        real(wp), dimension(2), intent(inout) :: var_loc

#ifdef MFC_MPI
        integer                :: ierr     !< Generic flag used to identify and report MPI errors
        real(wp), dimension(2) :: var_glb  !< Reduced (max value, rank) pair
        call MPI_REDUCE(var_loc, var_glb, 1, mpi_2p, MPI_MAXLOC, 0, MPI_COMM_WORLD, ierr)

        call MPI_BCAST(var_glb, 1, mpi_2p, 0, MPI_COMM_WORLD, ierr)

        var_loc = var_glb
#endif

    end subroutine s_mpi_reduce_maxloc

    !> The subroutine terminates the MPI execution environment.
    impure subroutine s_mpi_abort(prnt, code)

        character(len=*), intent(in), optional :: prnt
        integer, intent(in), optional          :: code

#ifdef MFC_MPI
        integer :: ierr  !< Generic flag used to identify and report MPI errors
#endif

        if (present(prnt)) then
            print *, prnt
            call flush (6)
        end if

#ifndef MFC_MPI
        if (present(code)) then
            stop code
        else
            stop 1
        end if
#else
        if (present(code)) then
            call MPI_ABORT(MPI_COMM_WORLD, code, ierr)
        else
            call MPI_ABORT(MPI_COMM_WORLD, 1, ierr)
        end if
#endif

    end subroutine s_mpi_abort

    !> Halts all processes until all have reached barrier.
    impure subroutine s_mpi_barrier

#ifdef MFC_MPI
        integer :: ierr  !< Generic flag used to identify and report MPI errors

        call MPI_BARRIER(MPI_COMM_WORLD, ierr)
#endif

    end subroutine s_mpi_barrier

    !> The subroutine finalizes the MPI execution environment.
    impure subroutine s_mpi_finalize

#ifdef MFC_MPI
        integer :: ierr  !< Generic flag used to identify and report MPI errors

        call MPI_FINALIZE(ierr)
#endif

    end subroutine s_mpi_finalize

    !> The goal of this procedure is to populate the buffers of the cell-average conservative variables by communicating with the
    !! neighboring processors.
    subroutine s_mpi_sendrecv_variables_buffers(q_comm, mpi_dir, pbc_loc, nVar, pb_in, mv_in, q_T_sf)

        type(scalar_field), dimension(1:), intent(inout) :: q_comm
        real(stp), optional, dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), intent(inout) :: pb_in, mv_in
        integer, intent(in) :: mpi_dir, pbc_loc, nVar
        integer :: i, j, k, l, r, q  !< Generic loop iterators
        integer :: buffer_counts(1:3), buffer_count
        type(int_bounds_info) :: boundary_conditions(1:3)
        integer :: beg_end(1:2), grid_dims(1:3)
        integer :: dst_proc, src_proc, recv_tag, send_tag
        logical :: beg_end_geq_0, qbmm_comm, chem_diff_comm
        integer :: pack_offset, unpack_offset
        type(scalar_field), optional, intent(inout) :: q_T_sf

#ifdef MFC_MPI
        integer :: ierr  !< Generic flag used to identify and report MPI errors

        call nvtxStartRange("RHS-COMM-PACKBUF")

        qbmm_comm = .false.
        chem_diff_comm = .false.

        if (present(pb_in) .and. present(mv_in) .and. qbmm .and. .not. polytropic) then
            qbmm_comm = .true.
            v_size = nVar + 2*nb*nnode
            buffer_counts = (/buff_size*v_size*(n + 1)*(p + 1), buff_size*v_size*(m + 2*buff_size + 1)*(p + 1), &
                             & buff_size*v_size*(m + 2*buff_size + 1)*(n + 2*buff_size + 1)/)
        else if (present(q_T_sf) .and. chemistry .and. (chem_params%diffusion .or. exchange_all_chemistry_temperatures)) then
            ! Consumers that convert over ghost-inclusive bounds request temperature exchange for every chemistry run.
            ! The temperature Newton guess must be valid at rank seams even when diffusion is disabled:
            ! an unexchanged seam ghost is an uninitialized guess -> NaN T/pres/c in the output
            chem_diff_comm = .true.
            v_size = nVar + 1
            buffer_counts = (/buff_size*v_size*(n + 1)*(p + 1), buff_size*v_size*(m + 2*buff_size + 1)*(p + 1), &
                             & buff_size*v_size*(m + 2*buff_size + 1)*(n + 2*buff_size + 1)/)
        else
            v_size = nVar
            buffer_counts = (/buff_size*v_size*(n + 1)*(p + 1), buff_size*v_size*(m + 2*buff_size + 1)*(p + 1), &
                             & buff_size*v_size*(m + 2*buff_size + 1)*(n + 2*buff_size + 1)/)
        end if

        $:GPU_UPDATE(device='[v_size]')

        buffer_count = buffer_counts(mpi_dir)
        boundary_conditions = (/bc_x, bc_y, bc_z/)
        beg_end = (/boundary_conditions(mpi_dir)%beg, boundary_conditions(mpi_dir)%end/)
        beg_end_geq_0 = beg_end(max(pbc_loc, 0) - pbc_loc + 1) >= 0

        ! Implements: pbc_loc bc_x >= 0 -> [send/recv]_tag [dst/src]_proc -1 (=0) 0 -> [1,0] [0,0] | 0 0 [1,0] [beg,beg] -1 (=0) 1
        ! -> [0,0] [1,0] | 0 1 [0,0] [end,beg] +1 (=1) 0 -> [0,1] [1,1] | 1 0 [0,1] [end,end] +1 (=1) 1 -> [1,1] [0,1] | 1 1 [1,1]
        ! [beg,end]

        send_tag = f_logical_to_int(.not. f_xor(beg_end_geq_0, pbc_loc == 1))
        recv_tag = f_logical_to_int(pbc_loc == 1)

        dst_proc = beg_end(1 + f_logical_to_int(f_xor(pbc_loc == 1, beg_end_geq_0)))
        src_proc = beg_end(1 + f_logical_to_int(pbc_loc == 1))

        grid_dims = (/m, n, p/)

        pack_offset = 0
        if (f_xor(pbc_loc == 1, beg_end_geq_0)) then
            pack_offset = grid_dims(mpi_dir) - buff_size + 1
        end if

        unpack_offset = 0
        if (pbc_loc == 1) then
            unpack_offset = grid_dims(mpi_dir) + buff_size + 1
        end if

        ! Pack Buffer to Send
        #:for mpi_dir in [1, 2, 3]
            if (mpi_dir == ${mpi_dir}$) then
                #:if mpi_dir == 1
                    $:GPU_PARALLEL_LOOP(collapse=4,private='[r]')
                    do l = 0, p
                        do k = 0, n
                            do j = 0, buff_size - 1
                                do i = 1, nVar
                                    r = (i - 1) + v_size*(j + buff_size*(k + (n + 1)*l))
                                    buff_send(r) = real(q_comm(i)%sf(j + pack_offset, k, l), kind=wp)
                                end do
                            end do
                        end do
                    end do
                    $:END_GPU_PARALLEL_LOOP()

                    if (chem_diff_comm) then
                        $:GPU_PARALLEL_LOOP(collapse=3,private='[r]')
                        do l = 0, p
                            do k = 0, n
                                do j = 0, buff_size - 1
                                    r = nVar + v_size*(j + buff_size*(k + (n + 1)*l))
                                    buff_send(r) = real(q_T_sf%sf(j + pack_offset, k, l), kind=wp)
                                end do
                            end do
                        end do
                        $:END_GPU_PARALLEL_LOOP()
                    end if

                    if (qbmm_comm) then
                        $:GPU_PARALLEL_LOOP(collapse=4,private='[r]')
                        do l = 0, p
                            do k = 0, n
                                do j = 0, buff_size - 1
                                    do i = nVar + 1, nVar + nnode
                                        do q = 1, nb
                                            r = (i - 1) + (q - 1)*nnode + v_size*(j + buff_size*(k + (n + 1)*l))
                                            buff_send(r) = real(pb_in(j + pack_offset, k, l, i - nVar, q), kind=wp)
                                        end do
                                    end do
                                end do
                            end do
                        end do
                        $:END_GPU_PARALLEL_LOOP()

                        $:GPU_PARALLEL_LOOP(collapse=5,private='[r]')
                        do l = 0, p
                            do k = 0, n
                                do j = 0, buff_size - 1
                                    do i = nVar + 1, nVar + nnode
                                        do q = 1, nb
                                            r = (i - 1) + (q - 1)*nnode + nb*nnode + v_size*(j + buff_size*(k + (n + 1)*l))
                                            buff_send(r) = real(mv_in(j + pack_offset, k, l, i - nVar, q), kind=wp)
                                        end do
                                    end do
                                end do
                            end do
                        end do
                        $:END_GPU_PARALLEL_LOOP()
                    end if
                #:elif mpi_dir == 2
                    $:GPU_PARALLEL_LOOP(collapse=4,private='[r]')
                    do i = 1, nVar
                        do l = 0, p
                            do k = 0, buff_size - 1
                                do j = -buff_size, m + buff_size
                                    r = (i - 1) + v_size*((j + buff_size) + (m + 2*buff_size + 1)*(k + buff_size*l))
                                    buff_send(r) = real(q_comm(i)%sf(j, k + pack_offset, l), kind=wp)
                                end do
                            end do
                        end do
                    end do
                    $:END_GPU_PARALLEL_LOOP()

                    if (chem_diff_comm) then
                        $:GPU_PARALLEL_LOOP(collapse=3,private='[r]')
                        do l = 0, p
                            do k = 0, buff_size - 1
                                do j = -buff_size, m + buff_size
                                    r = nVar + v_size*((j + buff_size) + (m + 2*buff_size + 1)*(k + buff_size*l))
                                    buff_send(r) = real(q_T_sf%sf(j, k + pack_offset, l), kind=wp)
                                end do
                            end do
                        end do
                        $:END_GPU_PARALLEL_LOOP()
                    end if

                    if (qbmm_comm) then
                        $:GPU_PARALLEL_LOOP(collapse=5,private='[r]')
                        do i = nVar + 1, nVar + nnode
                            do l = 0, p
                                do k = 0, buff_size - 1
                                    do j = -buff_size, m + buff_size
                                        do q = 1, nb
                                            r = (i - 1) + (q - 1)*nnode + v_size*((j + buff_size) + (m + 2*buff_size + 1)*(k &
                                                 & + buff_size*l))
                                            buff_send(r) = real(pb_in(j, k + pack_offset, l, i - nVar, q), kind=wp)
                                        end do
                                    end do
                                end do
                            end do
                        end do
                        $:END_GPU_PARALLEL_LOOP()

                        $:GPU_PARALLEL_LOOP(collapse=5,private='[r]')
                        do i = nVar + 1, nVar + nnode
                            do l = 0, p
                                do k = 0, buff_size - 1
                                    do j = -buff_size, m + buff_size
                                        do q = 1, nb
                                            r = (i - 1) + (q - 1)*nnode + nb*nnode + v_size*((j + buff_size) + (m + 2*buff_size &
                                                 & + 1)*(k + buff_size*l))
                                            buff_send(r) = real(mv_in(j, k + pack_offset, l, i - nVar, q), kind=wp)
                                        end do
                                    end do
                                end do
                            end do
                        end do
                        $:END_GPU_PARALLEL_LOOP()
                    end if
                #:else
                    $:GPU_PARALLEL_LOOP(collapse=4,private='[r]')
                    do i = 1, nVar
                        do l = 0, buff_size - 1
                            do k = -buff_size, n + buff_size
                                do j = -buff_size, m + buff_size
                                    r = (i - 1) + v_size*((j + buff_size) + (m + 2*buff_size + 1)*((k + buff_size) + (n &
                                         & + 2*buff_size + 1)*l))
                                    buff_send(r) = real(q_comm(i)%sf(j, k, l + pack_offset), kind=wp)
                                end do
                            end do
                        end do
                    end do
                    $:END_GPU_PARALLEL_LOOP()

                    if (chem_diff_comm) then
                        $:GPU_PARALLEL_LOOP(collapse=3,private='[r]')
                        do l = 0, buff_size - 1
                            do k = -buff_size, n + buff_size
                                do j = -buff_size, m + buff_size
                                    r = nVar + v_size*((j + buff_size) + (m + 2*buff_size + 1)*((k + buff_size) + (n &
                                                       & + 2*buff_size + 1)*l))
                                    buff_send(r) = real(q_T_sf%sf(j, k, l + pack_offset), kind=wp)
                                end do
                            end do
                        end do
                        $:END_GPU_PARALLEL_LOOP()
                    end if

                    if (qbmm_comm) then
                        $:GPU_PARALLEL_LOOP(collapse=5,private='[r]')
                        do i = nVar + 1, nVar + nnode
                            do l = 0, buff_size - 1
                                do k = -buff_size, n + buff_size
                                    do j = -buff_size, m + buff_size
                                        do q = 1, nb
                                            r = (i - 1) + (q - 1)*nnode + v_size*((j + buff_size) + (m + 2*buff_size + 1)*((k &
                                                 & + buff_size) + (n + 2*buff_size + 1)*l))
                                            buff_send(r) = real(pb_in(j, k, l + pack_offset, i - nVar, q), kind=wp)
                                        end do
                                    end do
                                end do
                            end do
                        end do
                        $:END_GPU_PARALLEL_LOOP()

                        $:GPU_PARALLEL_LOOP(collapse=5,private='[r]')
                        do i = nVar + 1, nVar + nnode
                            do l = 0, buff_size - 1
                                do k = -buff_size, n + buff_size
                                    do j = -buff_size, m + buff_size
                                        do q = 1, nb
                                            r = (i - 1) + (q - 1)*nnode + nb*nnode + v_size*((j + buff_size) + (m + 2*buff_size &
                                                 & + 1)*((k + buff_size) + (n + 2*buff_size + 1)*l))
                                            buff_send(r) = real(mv_in(j, k, l + pack_offset, i - nVar, q), kind=wp)
                                        end do
                                    end do
                                end do
                            end do
                        end do
                        $:END_GPU_PARALLEL_LOOP()
                    end if
                #:endif
            end if
        #:endfor
        call nvtxEndRange  ! Packbuf

        ! Send/Recv
        #:for rdma_mpi in [False, True]
            if (use_rdma_transport .eqv. ${'.true.' if rdma_mpi else '.false.'}$) then
                #:if rdma_mpi
                    #:call GPU_HOST_DATA(use_device_addr='[buff_send, buff_recv]')
                        call nvtxStartRange("RHS-COMM-SENDRECV-RDMA")

                        call MPI_SENDRECV(buff_send, buffer_count, mpi_p, dst_proc, send_tag, buff_recv, buffer_count, mpi_p, &
                                          & src_proc, recv_tag, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr)

                        call nvtxEndRange  ! RHS-MPI-SENDRECV-(NO)-RDMA
                    #:endcall GPU_HOST_DATA
                    $:GPU_WAIT()
                #:else
                    call nvtxStartRange("RHS-COMM-DEV2HOST")
                    $:GPU_UPDATE(host='[buff_send]')
                    call nvtxEndRange
                    call nvtxStartRange("RHS-COMM-SENDRECV-NO-RMDA")

                    call MPI_SENDRECV(buff_send, buffer_count, mpi_p, dst_proc, send_tag, buff_recv, buffer_count, mpi_p, &
                                      & src_proc, recv_tag, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr)

                    call nvtxEndRange  ! RHS-MPI-SENDRECV-(NO)-RDMA

                    call nvtxStartRange("RHS-COMM-HOST2DEV")
                    $:GPU_UPDATE(device='[buff_recv]')
                    call nvtxEndRange
                #:endif
            end if
        #:endfor

        ! Unpack Received Buffer
        call nvtxStartRange("RHS-COMM-UNPACKBUF")
        #:for mpi_dir in [1, 2, 3]
            if (mpi_dir == ${mpi_dir}$) then
                #:if mpi_dir == 1
                    $:GPU_PARALLEL_LOOP(collapse=4,private='[r]')
                    do l = 0, p
                        do k = 0, n
                            do j = -buff_size, -1
                                do i = 1, nVar
                                    r = (i - 1) + v_size*(j + buff_size*((k + 1) + (n + 1)*l))
                                    q_comm(i)%sf(j + unpack_offset, k, l) = real(buff_recv(r), kind=stp)
#if defined(__INTEL_COMPILER)
                                    if (ieee_is_nan(q_comm(i)%sf(j + unpack_offset, k, l))) then
                                        print *, "Error", j, k, l, i
                                        call s_mpi_abort("NaN(s) in recv")
                                    end if
#endif
                                end do
                            end do
                        end do
                    end do
                    $:END_GPU_PARALLEL_LOOP()

                    if (chem_diff_comm) then
                        $:GPU_PARALLEL_LOOP(collapse=3,private='[r]')
                        do l = 0, p
                            do k = 0, n
                                do j = -buff_size, -1
                                    r = nVar + v_size*(j + buff_size*((k + 1) + (n + 1)*l))
                                    q_T_sf%sf(j + unpack_offset, k, l) = real(buff_recv(r), kind=stp)
#if defined(__INTEL_COMPILER)
                                    if (ieee_is_nan(q_T_sf%sf(j + unpack_offset, k, l))) then
                                        print *, "Error", j, k, l
                                        call s_mpi_abort("NaN(s) in recv")
                                    end if
#endif
                                end do
                            end do
                        end do
                        $:END_GPU_PARALLEL_LOOP()
                    end if

                    if (qbmm_comm) then
                        $:GPU_PARALLEL_LOOP(collapse=5,private='[r]')
                        do l = 0, p
                            do k = 0, n
                                do j = -buff_size, -1
                                    do i = nVar + 1, nVar + nnode
                                        do q = 1, nb
                                            r = (i - 1) + (q - 1)*nnode + v_size*(j + buff_size*((k + 1) + (n + 1)*l))
                                            pb_in(j + unpack_offset, k, l, i - nVar, q) = real(buff_recv(r), kind=stp)
                                        end do
                                    end do
                                end do
                            end do
                        end do
                        $:END_GPU_PARALLEL_LOOP()

                        $:GPU_PARALLEL_LOOP(collapse=5,private='[r]')
                        do l = 0, p
                            do k = 0, n
                                do j = -buff_size, -1
                                    do i = nVar + 1, nVar + nnode
                                        do q = 1, nb
                                            r = (i - 1) + (q - 1)*nnode + nb*nnode + v_size*(j + buff_size*((k + 1) + (n + 1)*l))
                                            mv_in(j + unpack_offset, k, l, i - nVar, q) = real(buff_recv(r), kind=stp)
                                        end do
                                    end do
                                end do
                            end do
                        end do
                        $:END_GPU_PARALLEL_LOOP()
                    end if
                #:elif mpi_dir == 2
                    $:GPU_PARALLEL_LOOP(collapse=4,private='[r]')
                    do i = 1, nVar
                        do l = 0, p
                            do k = -buff_size, -1
                                do j = -buff_size, m + buff_size
                                    r = (i - 1) + v_size*((j + buff_size) + (m + 2*buff_size + 1)*((k + buff_size) + buff_size*l))
                                    q_comm(i)%sf(j, k + unpack_offset, l) = real(buff_recv(r), kind=stp)
#if defined(__INTEL_COMPILER)
                                    if (ieee_is_nan(q_comm(i)%sf(j, k + unpack_offset, l))) then
                                        print *, "Error", j, k, l, i
                                        call s_mpi_abort("NaN(s) in recv")
                                    end if
#endif
                                end do
                            end do
                        end do
                    end do
                    $:END_GPU_PARALLEL_LOOP()

                    if (chem_diff_comm) then
                        $:GPU_PARALLEL_LOOP(collapse=3,private='[r]')
                        do l = 0, p
                            do k = -buff_size, -1
                                do j = -buff_size, m + buff_size
                                    r = nVar + v_size*((j + buff_size) + (m + 2*buff_size + 1)*((k + buff_size) + buff_size*l))
                                    q_T_sf%sf(j, k + unpack_offset, l) = real(buff_recv(r), kind=stp)
#if defined(__INTEL_COMPILER)
                                    if (ieee_is_nan(q_T_sf%sf(j, k + unpack_offset, l))) then
                                        print *, "Error", j, k, l
                                        call s_mpi_abort("NaN(s) in recv")
                                    end if
#endif
                                end do
                            end do
                        end do
                        $:END_GPU_PARALLEL_LOOP()
                    end if

                    if (qbmm_comm) then
                        $:GPU_PARALLEL_LOOP(collapse=5,private='[r]')
                        do i = nVar + 1, nVar + nnode
                            do l = 0, p
                                do k = -buff_size, -1
                                    do j = -buff_size, m + buff_size
                                        do q = 1, nb
                                            r = (i - 1) + (q - 1)*nnode + v_size*((j + buff_size) + (m + 2*buff_size + 1)*((k &
                                                 & + buff_size) + buff_size*l))
                                            pb_in(j, k + unpack_offset, l, i - nVar, q) = real(buff_recv(r), kind=stp)
                                        end do
                                    end do
                                end do
                            end do
                        end do
                        $:END_GPU_PARALLEL_LOOP()

                        $:GPU_PARALLEL_LOOP(collapse=5,private='[r]')
                        do i = nVar + 1, nVar + nnode
                            do l = 0, p
                                do k = -buff_size, -1
                                    do j = -buff_size, m + buff_size
                                        do q = 1, nb
                                            r = (i - 1) + (q - 1)*nnode + nb*nnode + v_size*((j + buff_size) + (m + 2*buff_size &
                                                 & + 1)*((k + buff_size) + buff_size*l))
                                            mv_in(j, k + unpack_offset, l, i - nVar, q) = real(buff_recv(r), kind=stp)
                                        end do
                                    end do
                                end do
                            end do
                        end do
                        $:END_GPU_PARALLEL_LOOP()
                    end if
                #:else
                    $:GPU_PARALLEL_LOOP(collapse=4,private='[r]')
                    do i = 1, nVar
                        do l = -buff_size, -1
                            do k = -buff_size, n + buff_size
                                do j = -buff_size, m + buff_size
                                    r = (i - 1) + v_size*((j + buff_size) + (m + 2*buff_size + 1)*((k + buff_size) + (n &
                                         & + 2*buff_size + 1)*(l + buff_size)))
                                    q_comm(i)%sf(j, k, l + unpack_offset) = real(buff_recv(r), kind=stp)
#if defined(__INTEL_COMPILER)
                                    if (ieee_is_nan(q_comm(i)%sf(j, k, l + unpack_offset))) then
                                        print *, "Error", j, k, l, i
                                        call s_mpi_abort("NaN(s) in recv")
                                    end if
#endif
                                end do
                            end do
                        end do
                    end do
                    $:END_GPU_PARALLEL_LOOP()

                    if (chem_diff_comm) then
                        $:GPU_PARALLEL_LOOP(collapse=3,private='[r]')
                        do l = -buff_size, -1
                            do k = -buff_size, n + buff_size
                                do j = -buff_size, m + buff_size
                                    r = nVar + v_size*((j + buff_size) + (m + 2*buff_size + 1)*((k + buff_size) + (n &
                                                       & + 2*buff_size + 1)*(l + buff_size)))
                                    q_T_sf%sf(j, k, l + unpack_offset) = real(buff_recv(r), kind=stp)
#if defined(__INTEL_COMPILER)
                                    if (ieee_is_nan(q_T_sf%sf(j, k, l + unpack_offset))) then
                                        print *, "Error", j, k, l
                                        call s_mpi_abort("NaN(s) in recv")
                                    end if
#endif
                                end do
                            end do
                        end do
                        $:END_GPU_PARALLEL_LOOP()
                    end if

                    if (qbmm_comm) then
                        $:GPU_PARALLEL_LOOP(collapse=5,private='[r]')
                        do i = nVar + 1, nVar + nnode
                            do l = -buff_size, -1
                                do k = -buff_size, n + buff_size
                                    do j = -buff_size, m + buff_size
                                        do q = 1, nb
                                            r = (i - 1) + (q - 1)*nnode + v_size*((j + buff_size) + (m + 2*buff_size + 1)*((k &
                                                 & + buff_size) + (n + 2*buff_size + 1)*(l + buff_size)))
                                            pb_in(j, k, l + unpack_offset, i - nVar, q) = real(buff_recv(r), kind=stp)
                                        end do
                                    end do
                                end do
                            end do
                        end do
                        $:END_GPU_PARALLEL_LOOP()

                        $:GPU_PARALLEL_LOOP(collapse=5,private='[r]')
                        do i = nVar + 1, nVar + nnode
                            do l = -buff_size, -1
                                do k = -buff_size, n + buff_size
                                    do j = -buff_size, m + buff_size
                                        do q = 1, nb
                                            r = (i - 1) + (q - 1)*nnode + nb*nnode + v_size*((j + buff_size) + (m + 2*buff_size &
                                                 & + 1)*((k + buff_size) + (n + 2*buff_size + 1)*(l + buff_size)))
                                            mv_in(j, k, l + unpack_offset, i - nVar, q) = real(buff_recv(r), kind=stp)
                                        end do
                                    end do
                                end do
                            end do
                        end do
                        $:END_GPU_PARALLEL_LOOP()
                    end if
                #:endif
            end if
        #:endfor
        call nvtxEndRange
#endif

    end subroutine s_mpi_sendrecv_variables_buffers

    !> The goal of this procedure is to populate the buffers of the cell-average conservative variables by communicating with the
    !! neighboring processors.
    !! @param q_cons_vf Cell-average conservative variables
    !! @param mpi_dir MPI communication coordinate direction
    !! @param pbc_loc Processor boundary condition (PBC) location
    subroutine s_mpi_reduce_beta_variables_buffers(q_comm, kahan_comp, mpi_dir, pbc_loc, nVar)

        type(scalar_field), dimension(1:), intent(inout) :: q_comm
        type(scalar_field), dimension(1:), intent(inout) :: kahan_comp
        integer, intent(in)                              :: mpi_dir, pbc_loc, nVar
        integer                                          :: i, j, k, l, r, q  !< Generic loop iterators
        integer                                          :: lb_size
        integer                                          :: buffer_counts(1:3), buffer_count
        type(int_bounds_info)                            :: boundary_conditions(1:3)
        integer                                          :: beg_end(1:2), grid_dims(1:3)
        integer                                          :: dst_proc, src_proc, recv_tag, send_tag
        logical                                          :: replace_buff
        integer                                          :: pack_offset, unpack_offset
        real(wp)                                         :: y_kahan, t_kahan

#ifdef MFC_MPI
        integer :: ierr  !< Generic flag used to identify and report MPI errors

        call nvtxStartRange("BETA-COMM-PACKBUF")

        ! Set bounds for each dimension Always include the full buffer range for each existing dimension. The Gaussian smearing
        ! kernel writes to buffer cells even at physical boundaries, and these contributions must be communicated to neighbors in
        ! other directions via ADD operations.
        comm_coords(1)%beg = -mapcells - 1
        comm_coords(1)%end = m + mapcells + 1
        comm_coords(2)%beg = merge(-mapcells - 1, 0, n > 0)
        comm_coords(2)%end = merge(n + mapcells + 1, n, n > 0)
        comm_coords(3)%beg = merge(-mapcells - 1, 0, p > 0)
        comm_coords(3)%end = merge(p + mapcells + 1, p, p > 0)

        ! Compute sizes
        comm_size(1) = comm_coords(1)%end - comm_coords(1)%beg + 1
        comm_size(2) = comm_coords(2)%end - comm_coords(2)%beg + 1
        comm_size(3) = comm_coords(3)%end - comm_coords(3)%beg + 1

        ! Buffer counts using the conditional sizes
        v_size = nVar
        lb_size = 2*(mapcells + 1)  ! Size of the buffer region for beta variables (-mapcells - 1, mapcells)
        buffer_counts = (/lb_size*v_size*comm_size(2)*comm_size(3), lb_size*v_size*comm_size(1)*comm_size(3), &
                         & lb_size*v_size*comm_size(1)*comm_size(2)/)

        $:GPU_UPDATE(device='[v_size, comm_coords, comm_size]')

        buffer_count = buffer_counts(mpi_dir)
        boundary_conditions = (/bc_x, bc_y, bc_z/)
        beg_end = (/boundary_conditions(mpi_dir)%beg, boundary_conditions(mpi_dir)%end/)
        grid_dims = (/m, n, p/)

        if (pbc_loc == -1) then  ! PBC at the beginning
            ! Phase 1: Rightward accumulation Send END buffer to right neighbor, recv from left into BEG, ADD
            pack_offset = grid_dims(mpi_dir) + 1
            unpack_offset = 0
            dst_proc = merge(beg_end(2), MPI_PROC_NULL, beg_end(2) >= 0)
            src_proc = merge(beg_end(1), MPI_PROC_NULL, beg_end(1) >= 0)
            send_tag = 0
            recv_tag = 0
            replace_buff = .false.
        else
            ! Phase 2: Leftward distribution Send BEG buffer to left neighbor, recv from right into END, REPLACE
            pack_offset = 0
            unpack_offset = grid_dims(mpi_dir) + 1
            dst_proc = merge(beg_end(1), MPI_PROC_NULL, beg_end(1) >= 0)
            src_proc = merge(beg_end(2), MPI_PROC_NULL, beg_end(2) >= 0)
            send_tag = 1
            recv_tag = 1
            replace_buff = .true.
        end if

        ! Pack Buffer to Send
        #:for mpi_dir in [1, 2, 3]
            if (mpi_dir == ${mpi_dir}$) then
                #:if mpi_dir == 1
                    $:GPU_PARALLEL_LOOP(collapse=4,private='[r]')
                    do l = comm_coords(3)%beg, comm_coords(3)%end
                        do k = comm_coords(2)%beg, comm_coords(2)%end
                            do j = -mapcells - 1, mapcells
                                do i = 1, v_size
                                    r = (i - 1) + v_size*((j + mapcells + 1) + lb_size*((k - comm_coords(2)%beg) + comm_size(2) &
                                         & *(l - comm_coords(3)%beg)))
                                    buff_send(r) = real(q_comm(beta_vars(i))%sf(j + pack_offset, k, l), &
                                              & kind=wp) - real(kahan_comp(beta_vars(i))%sf(j + pack_offset, k, l), kind=wp)
                                end do
                            end do
                        end do
                    end do
                    $:END_GPU_PARALLEL_LOOP()
                #:elif mpi_dir == 2
                    $:GPU_PARALLEL_LOOP(collapse=4,private='[r]')
                    do i = 1, v_size
                        do l = comm_coords(3)%beg, comm_coords(3)%end
                            do k = -mapcells - 1, mapcells
                                do j = comm_coords(1)%beg, comm_coords(1)%end
                                    r = (i - 1) + v_size*((j - comm_coords(1)%beg) + comm_size(1)*((k + mapcells + 1) &
                                         & + lb_size*(l - comm_coords(3)%beg)))
                                    buff_send(r) = real(q_comm(beta_vars(i))%sf(j, k + pack_offset, l), &
                                              & kind=wp) - real(kahan_comp(beta_vars(i))%sf(j, k + pack_offset, l), kind=wp)
                                end do
                            end do
                        end do
                    end do
                    $:END_GPU_PARALLEL_LOOP()
                #:else
                    $:GPU_PARALLEL_LOOP(collapse=4,private='[r]')
                    do i = 1, v_size
                        do l = -mapcells - 1, mapcells
                            do k = comm_coords(2)%beg, comm_coords(2)%end
                                do j = comm_coords(1)%beg, comm_coords(1)%end
                                    r = (i - 1) + v_size*((j - comm_coords(1)%beg) + comm_size(1)*((k - comm_coords(2)%beg) &
                                         & + comm_size(2)*(l + mapcells + 1)))
                                    buff_send(r) = real(q_comm(beta_vars(i))%sf(j, k, l + pack_offset), &
                                              & kind=wp) - real(kahan_comp(beta_vars(i))%sf(j, k, l + pack_offset), kind=wp)
                                end do
                            end do
                        end do
                    end do
                    $:END_GPU_PARALLEL_LOOP()
                #:endif
            end if
        #:endfor
        call nvtxEndRange  ! Packbuf

        ! Send/Recv
        #:for rdma_mpi in [False, True]
            if (use_rdma_transport .eqv. ${'.true.' if rdma_mpi else '.false.'}$) then
                #:if rdma_mpi
                    #:call GPU_HOST_DATA(use_device_addr='[buff_send, buff_recv]')
                        call nvtxStartRange("BETA-COMM-SENDRECV-RDMA")

                        call MPI_SENDRECV(buff_send, buffer_count, mpi_p, dst_proc, send_tag, buff_recv, buffer_count, mpi_p, &
                                          & src_proc, recv_tag, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr)

                        call nvtxEndRange  ! BETA-MPI-SENDRECV-(NO)-RDMA
                    #:endcall GPU_HOST_DATA
                    $:GPU_WAIT()
                #:else
                    call nvtxStartRange("BETA-COMM-DEV2HOST")
                    $:GPU_UPDATE(host='[buff_send]')
                    call nvtxEndRange
                    call nvtxStartRange("BETA-COMM-SENDRECV-NO-RMDA")

                    call MPI_SENDRECV(buff_send, buffer_count, mpi_p, dst_proc, send_tag, buff_recv, buffer_count, mpi_p, &
                                      & src_proc, recv_tag, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr)

                    call nvtxEndRange  ! BETA-MPI-SENDRECV-(NO)-RDMA

                    call nvtxStartRange("BETA-COMM-HOST2DEV")
                    $:GPU_UPDATE(device='[buff_recv]')
                    call nvtxEndRange
                #:endif
            end if
        #:endfor

        ! Unpack Received Buffer (skip if no source rank)
        call nvtxStartRange("BETA-COMM-UNPACKBUF")
        if (src_proc /= MPI_PROC_NULL) then
            #:for mpi_dir in [1, 2, 3]
                if (mpi_dir == ${mpi_dir}$) then
                    #:if mpi_dir == 1
                        $:GPU_PARALLEL_LOOP(collapse=4,private='[r, y_kahan, t_kahan]',copyin='[replace_buff]')
                        do l = comm_coords(3)%beg, comm_coords(3)%end
                            do k = comm_coords(2)%beg, comm_coords(2)%end
                                do j = -mapcells - 1, mapcells
                                    do i = 1, v_size
                                        r = (i - 1) + v_size*((j + mapcells + 1) + lb_size*((k - comm_coords(2)%beg) &
                                             & + comm_size(2)*(l - comm_coords(3)%beg)))
                                        if (replace_buff) then
                                            q_comm(beta_vars(i))%sf(j + unpack_offset, k, l) = real(buff_recv(r), kind=stp)
                                            kahan_comp(beta_vars(i))%sf(j + unpack_offset, k, &
                                                       & l) = real(q_comm(beta_vars(i))%sf(j + unpack_offset, k, l), &
                                                       & kind=wp) - buff_recv(r)
                                        else
                                            y_kahan = buff_recv(r) - real(kahan_comp(beta_vars(i))%sf(j + unpack_offset, k, l), &
                                                                & kind=wp)
                                            t_kahan = real(q_comm(beta_vars(i))%sf(j + unpack_offset, k, l), kind=wp) + y_kahan
                                            kahan_comp(beta_vars(i))%sf(j + unpack_offset, k, &
                                                       & l) = (t_kahan - q_comm(beta_vars(i))%sf(j + unpack_offset, k, l)) - y_kahan
                                            q_comm(beta_vars(i))%sf(j + unpack_offset, k, l) = t_kahan
                                        end if
                                    end do
                                end do
                            end do
                        end do
                        $:END_GPU_PARALLEL_LOOP()
                    #:elif mpi_dir == 2
                        $:GPU_PARALLEL_LOOP(collapse=4,private='[r, y_kahan, t_kahan]',copyin='[replace_buff]')
                        do i = 1, v_size
                            do l = comm_coords(3)%beg, comm_coords(3)%end
                                do k = -mapcells - 1, mapcells
                                    do j = comm_coords(1)%beg, comm_coords(1)%end
                                        r = (i - 1) + v_size*((j - comm_coords(1)%beg) + comm_size(1)*((k + mapcells + 1) &
                                             & + lb_size*(l - comm_coords(3)%beg)))
                                        if (replace_buff) then
                                            q_comm(beta_vars(i))%sf(j, k + unpack_offset, l) = real(buff_recv(r), kind=stp)
                                            kahan_comp(beta_vars(i))%sf(j, k + unpack_offset, &
                                                       & l) = real(q_comm(beta_vars(i))%sf(j, k + unpack_offset, l), &
                                                       & kind=wp) - buff_recv(r)
                                        else
                                            y_kahan = buff_recv(r) - real(kahan_comp(beta_vars(i))%sf(j, k + unpack_offset, l), &
                                                                & kind=wp)
                                            t_kahan = real(q_comm(beta_vars(i))%sf(j, k + unpack_offset, l), kind=wp) + y_kahan
                                            kahan_comp(beta_vars(i))%sf(j, k + unpack_offset, &
                                                       & l) = (t_kahan - q_comm(beta_vars(i))%sf(j, k + unpack_offset, l)) - y_kahan
                                            q_comm(beta_vars(i))%sf(j, k + unpack_offset, l) = t_kahan
                                        end if
                                    end do
                                end do
                            end do
                        end do
                        $:END_GPU_PARALLEL_LOOP()
                    #:else
                        $:GPU_PARALLEL_LOOP(collapse=4,private='[r, y_kahan, t_kahan]',copyin='[replace_buff]')
                        do i = 1, v_size
                            do l = -mapcells - 1, mapcells
                                do k = comm_coords(2)%beg, comm_coords(2)%end
                                    do j = comm_coords(1)%beg, comm_coords(1)%end
                                        r = (i - 1) + v_size*((j - comm_coords(1)%beg) + comm_size(1)*((k - comm_coords(2)%beg) &
                                             & + comm_size(2)*(l + mapcells + 1)))
                                        if (replace_buff) then
                                            q_comm(beta_vars(i))%sf(j, k, l + unpack_offset) = real(buff_recv(r), kind=stp)
                                            kahan_comp(beta_vars(i))%sf(j, k, &
                                                       & l + unpack_offset) = real(q_comm(beta_vars(i))%sf(j, k, &
                                                       & l + unpack_offset), kind=wp) - buff_recv(r)
                                        else
                                            y_kahan = buff_recv(r) - real(kahan_comp(beta_vars(i))%sf(j, k, l + unpack_offset), &
                                                                & kind=wp)
                                            t_kahan = real(q_comm(beta_vars(i))%sf(j, k, l + unpack_offset), kind=wp) + y_kahan
                                            kahan_comp(beta_vars(i))%sf(j, k, &
                                                       & l + unpack_offset) = (t_kahan - q_comm(beta_vars(i))%sf(j, k, &
                                                       & l + unpack_offset)) - y_kahan
                                            q_comm(beta_vars(i))%sf(j, k, l + unpack_offset) = t_kahan
                                        end if
                                    end do
                                end do
                            end do
                        end do
                        $:END_GPU_PARALLEL_LOOP()
                    #:endif
                end if
            #:endfor
        end if
        call nvtxEndRange
#endif

    end subroutine s_mpi_reduce_beta_variables_buffers

    !> The purpose of this procedure is to optimally decompose the computational domain among the available processors. This is
    !! performed by attempting to award each processor, in each of the coordinate directions, approximately the same number of
    !! cells, and then recomputing the affected global parameters.
    subroutine s_mpi_decompose_computational_domain(write_silo_ghost_offsets, adjust_local_domains, output_offsets, local_domains)

        logical, intent(in)                                          :: write_silo_ghost_offsets
        logical, intent(in)                                          :: adjust_local_domains
        type(int_bounds_info), dimension(3), intent(inout), optional :: output_offsets
        type(bounds_info), dimension(3), intent(inout), optional     :: local_domains

#ifdef MFC_MPI
        !> Non-optimal number of processors in the x-, y- and z-directions
        real(wp) :: tmp_num_procs_x, tmp_num_procs_y, tmp_num_procs_z
        real(wp) :: fct_min        !< Processor factorization (fct) minimization parameter
        integer  :: MPI_COMM_CART  !< Cartesian processor topology communicator
        integer  :: rem_cells      !< Remaining cells after distribution among processors
        integer  :: rem_cells_by_dim(3)
        integer  :: recon_order    !< WENO or MUSCL reconstruction order
        integer  :: i, j, k        !< Generic loop iterators
        integer  :: ierr           !< Generic flag used to identify and report MPI errors

        ! temp array to store neighbor rank coordinates
        integer, dimension(1:num_dims) :: neighbor_coords

        rem_cells_by_dim = 0

        ! Zeroing out communication needs for moving EL bubbles/particles
        nidx(1)%beg = 0; nidx(1)%end = 0
        nidx(2)%beg = 0; nidx(2)%end = 0
        nidx(3)%beg = 0; nidx(3)%end = 0

        if (recon_type == recon_type_weno) then
            recon_order = weno_order
        else
            recon_order = muscl_order
        end if

        if (num_procs == 1 .and. parallel_io) then
            do i = 1, num_dims
                start_idx(i) = 0
            end do
            return
        end if

        if (igr) then
            recon_order = igr_order
        end if

        ! 3D Cartesian Processor Topology
        if (n > 0) then
            if (p > 0) then
                if (fft_wrt) then
                    ! Initial estimate of optimal processor topology
                    num_procs_x = 1
                    num_procs_y = 1
                    num_procs_z = num_procs
                    ierr = -1

                    ! Benchmarking the quality of this initial guess
                    tmp_num_procs_y = num_procs_y
                    tmp_num_procs_z = num_procs_z
                    fct_min = 10._wp*abs((n + 1)/tmp_num_procs_y - (p + 1)/tmp_num_procs_z)

                    ! Optimization of the initial processor topology
                    do i = 1, num_procs
                        if (mod(num_procs, i) == 0 .and. (n + 1)/i >= num_stcls_min*recon_order) then
                            tmp_num_procs_y = i
                            tmp_num_procs_z = num_procs/i

                            if (fct_min >= abs((n + 1)/tmp_num_procs_y - (p + 1)/tmp_num_procs_z) .and. (p + 1) &
                                & /tmp_num_procs_z >= num_stcls_min*recon_order) then
                                num_procs_y = i
                                num_procs_z = num_procs/i
                                fct_min = abs((n + 1)/tmp_num_procs_y - (p + 1)/tmp_num_procs_z)
                                ierr = 0
                            end if
                        end if
                    end do
                else
                    if (cyl_coord .and. p > 0) then
                        ! Pencil blocking for cylindrical coordinates (Fourier filter near axis)

                        ! Initial values of the processor factorization optimization
                        num_procs_x = 1
                        num_procs_y = num_procs
                        num_procs_z = 1
                        ierr = -1

                        ! Computing minimization variable for these initial values
                        tmp_num_procs_x = num_procs_x
                        tmp_num_procs_y = num_procs_y
                        tmp_num_procs_z = num_procs_z
                        fct_min = 10._wp*abs((m + 1)/tmp_num_procs_x - (n + 1)/tmp_num_procs_y)

                        ! Searching for optimal computational domain distribution
                        do i = 1, num_procs
                            if (mod(num_procs, i) == 0 .and. (m + 1)/i >= num_stcls_min*recon_order) then
                                tmp_num_procs_x = i
                                tmp_num_procs_y = num_procs/i

                                if (fct_min >= abs((m + 1)/tmp_num_procs_x - (n + 1)/tmp_num_procs_y) .and. (n + 1) &
                                    & /tmp_num_procs_y >= num_stcls_min*recon_order) then
                                    num_procs_x = i
                                    num_procs_y = num_procs/i
                                    fct_min = abs((m + 1)/tmp_num_procs_x - (n + 1)/tmp_num_procs_y)
                                    ierr = 0
                                end if
                            end if
                        end do
                    else
                        ! Initial estimate of optimal processor topology
                        num_procs_x = 1
                        num_procs_y = 1
                        num_procs_z = num_procs
                        ierr = -1

                        ! Benchmarking the quality of this initial guess
                        tmp_num_procs_x = num_procs_x
                        tmp_num_procs_y = num_procs_y
                        tmp_num_procs_z = num_procs_z
                        fct_min = 10._wp*abs((m + 1)/tmp_num_procs_x - (n + 1)/tmp_num_procs_y) + 10._wp*abs((n + 1) &
                                             & /tmp_num_procs_y - (p + 1)/tmp_num_procs_z)

                        ! Optimization of the initial processor topology
                        do i = 1, num_procs
                            if (mod(num_procs, i) == 0 .and. (m + 1)/i >= num_stcls_min*recon_order) then
                                do j = 1, num_procs/i
                                    if (mod(num_procs/i, j) == 0 .and. (n + 1)/j >= num_stcls_min*recon_order) then
                                        tmp_num_procs_x = i
                                        tmp_num_procs_y = j
                                        tmp_num_procs_z = num_procs/(i*j)

                                        if (fct_min >= abs((m + 1)/tmp_num_procs_x - (n + 1)/tmp_num_procs_y) + abs((n + 1) &
                                            & /tmp_num_procs_y - (p + 1)/tmp_num_procs_z) .and. (p + 1) &
                                            & /tmp_num_procs_z >= num_stcls_min*recon_order) then
                                            num_procs_x = i
                                            num_procs_y = j
                                            num_procs_z = num_procs/(i*j)
                                            fct_min = abs((m + 1)/tmp_num_procs_x - (n + 1)/tmp_num_procs_y) + abs((n + 1) &
                                                          & /tmp_num_procs_y - (p + 1)/tmp_num_procs_z)
                                            ierr = 0
                                        end if
                                    end if
                                end do
                            end if
                        end do
                    end if
                end if

                ! Verifying that a valid decomposition of the computational domain has been established. If not, the simulation
                ! exits.
                if (proc_rank == 0 .and. ierr == -1) then
                    call s_mpi_abort('Unsupported combination of values ' // 'of num_procs, m, n, p and ' &
                                     & // 'weno/muscl/igr_order. Exiting.')
                end if

                ! Creating new communicator using the Cartesian topology
                call MPI_CART_CREATE(MPI_COMM_WORLD, 3, (/num_procs_x, num_procs_y, num_procs_z/), (/.true., .true., .true./), &
                                     & .false., MPI_COMM_CART, ierr)

                ! Finding the Cartesian coordinates of the local process
                call MPI_CART_COORDS(MPI_COMM_CART, proc_rank, 3, proc_coords, ierr)

                ! Global Parameters for z-direction

                ! Number of remaining cells
                rem_cells = mod(p + 1, num_procs_z)
                rem_cells_by_dim(3) = rem_cells

                ! Optimal number of cells per processor
                p = (p + 1)/num_procs_z - 1

                ! Distributing the remaining cells
                do i = 1, rem_cells
                    if (proc_coords(3) == i - 1) then
                        p = p + 1; exit
                    end if
                end do

                ! Boundary condition at the beginning
                if (proc_coords(3) > 0 .or. (bc_z%beg == BC_PERIODIC .and. num_procs_z > 1)) then
                    proc_coords(3) = proc_coords(3) - 1
                    call MPI_CART_RANK(MPI_COMM_CART, proc_coords, bc_z%beg, ierr)
                    proc_coords(3) = proc_coords(3) + 1
                    nidx(3)%beg = -1
                end if

                ! Boundary condition at the end
                if (proc_coords(3) < num_procs_z - 1 .or. (bc_z%end == BC_PERIODIC .and. num_procs_z > 1)) then
                    proc_coords(3) = proc_coords(3) + 1
                    call MPI_CART_RANK(MPI_COMM_CART, proc_coords, bc_z%end, ierr)
                    proc_coords(3) = proc_coords(3) - 1
                    nidx(3)%end = 1
                end if

                ! Beginning and end sub-domain boundary locations
                if (parallel_io) then
                    if (proc_coords(3) < rem_cells) then
                        start_idx(3) = (p + 1)*proc_coords(3)
                    else
                        start_idx(3) = (p + 1)*proc_coords(3) + rem_cells
                    end if
                end if

                ! 2D Cartesian Processor Topology
            else
                ! Initial estimate of optimal processor topology
                num_procs_x = 1
                num_procs_y = num_procs
                ierr = -1

                ! Benchmarking the quality of this initial guess
                tmp_num_procs_x = num_procs_x
                tmp_num_procs_y = num_procs_y
                fct_min = 10._wp*abs((m + 1)/tmp_num_procs_x - (n + 1)/tmp_num_procs_y)

                ! Optimization of the initial processor topology
                do i = 1, num_procs
                    if (mod(num_procs, i) == 0 .and. (m + 1)/i >= num_stcls_min*recon_order) then
                        tmp_num_procs_x = i
                        tmp_num_procs_y = num_procs/i

                        if (fct_min >= abs((m + 1)/tmp_num_procs_x - (n + 1)/tmp_num_procs_y) .and. (n + 1) &
                            & /tmp_num_procs_y >= num_stcls_min*recon_order) then
                            num_procs_x = i
                            num_procs_y = num_procs/i
                            fct_min = abs((m + 1)/tmp_num_procs_x - (n + 1)/tmp_num_procs_y)
                            ierr = 0
                        end if
                    end if
                end do

                ! Verifying that a valid decomposition of the computational domain has been established. If not, the simulation
                ! exits.
                if (proc_rank == 0 .and. ierr == -1) then
                    call s_mpi_abort('Unsupported combination of values ' // 'of num_procs, m, n and ' &
                                     & // 'weno/muscl/igr_order. Exiting.')
                end if

                ! Creating new communicator using the Cartesian topology
                call MPI_CART_CREATE(MPI_COMM_WORLD, 2, (/num_procs_x, num_procs_y/), (/.true., .true./), .false., MPI_COMM_CART, &
                                     & ierr)

                ! Finding the Cartesian coordinates of the local process
                call MPI_CART_COORDS(MPI_COMM_CART, proc_rank, 2, proc_coords, ierr)
            end if

            ! Global Parameters for y-direction

            ! Number of remaining cells
            rem_cells = mod(n + 1, num_procs_y)
            rem_cells_by_dim(2) = rem_cells

            ! Optimal number of cells per processor
            n = (n + 1)/num_procs_y - 1

            ! Distributing the remaining cells
            do i = 1, rem_cells
                if (proc_coords(2) == i - 1) then
                    n = n + 1; exit
                end if
            end do

            ! Boundary condition at the beginning
            if (proc_coords(2) > 0 .or. (bc_y%beg == BC_PERIODIC .and. num_procs_y > 1)) then
                proc_coords(2) = proc_coords(2) - 1
                call MPI_CART_RANK(MPI_COMM_CART, proc_coords, bc_y%beg, ierr)
                proc_coords(2) = proc_coords(2) + 1
                nidx(2)%beg = -1
            end if

            ! Boundary condition at the end
            if (proc_coords(2) < num_procs_y - 1 .or. (bc_y%end == BC_PERIODIC .and. num_procs_y > 1)) then
                proc_coords(2) = proc_coords(2) + 1
                call MPI_CART_RANK(MPI_COMM_CART, proc_coords, bc_y%end, ierr)
                proc_coords(2) = proc_coords(2) - 1
                nidx(2)%end = 1
            end if

            ! Beginning and end sub-domain boundary locations
            if (parallel_io) then
                if (proc_coords(2) < rem_cells) then
                    start_idx(2) = (n + 1)*proc_coords(2)
                else
                    start_idx(2) = (n + 1)*proc_coords(2) + rem_cells
                end if
            end if

            ! 1D Cartesian Processor Topology
        else
            ! Optimal processor topology
            num_procs_x = num_procs

            ! Creating new communicator using the Cartesian topology
            call MPI_CART_CREATE(MPI_COMM_WORLD, 1, (/num_procs_x/), (/.true./), .false., MPI_COMM_CART, ierr)

            ! Finding the Cartesian coordinates of the local process
            call MPI_CART_COORDS(MPI_COMM_CART, proc_rank, 1, proc_coords, ierr)
        end if

        ! Global Parameters for x-direction

        ! Number of remaining cells
        rem_cells = mod(m + 1, num_procs_x)
        rem_cells_by_dim(1) = rem_cells

        ! Optimal number of cells per processor
        m = (m + 1)/num_procs_x - 1

        ! Distributing the remaining cells
        do i = 1, rem_cells
            if (proc_coords(1) == i - 1) then
                m = m + 1; exit
            end if
        end do

        call s_update_cell_bounds(cells_bounds, m, n, p)

        ! Boundary condition at the beginning
        if (proc_coords(1) > 0 .or. (bc_x%beg == BC_PERIODIC .and. num_procs_x > 1)) then
            proc_coords(1) = proc_coords(1) - 1
            call MPI_CART_RANK(MPI_COMM_CART, proc_coords, bc_x%beg, ierr)
            proc_coords(1) = proc_coords(1) + 1
            nidx(1)%beg = -1
        end if

        ! Boundary condition at the end
        if (proc_coords(1) < num_procs_x - 1 .or. (bc_x%end == BC_PERIODIC .and. num_procs_x > 1)) then
            proc_coords(1) = proc_coords(1) + 1
            call MPI_CART_RANK(MPI_COMM_CART, proc_coords, bc_x%end, ierr)
            proc_coords(1) = proc_coords(1) - 1
            nidx(1)%end = 1
        end if

        ! Beginning and end sub-domain boundary locations
        if (parallel_io) then
            if (proc_coords(1) < rem_cells) then
                start_idx(1) = (m + 1)*proc_coords(1)
            else
                start_idx(1) = (m + 1)*proc_coords(1) + rem_cells
            end if
        end if

        call s_apply_decomposition_policies((/num_procs_x, num_procs_y, num_procs_z/), rem_cells_by_dim, (/m, n, p/), (/m_glb, &
                                            & n_glb, p_glb/), write_silo_ghost_offsets, &
                                            & adjust_local_domains .and. (.not. parallel_io), output_offsets, local_domains)

        @:ALLOCATE(neighbor_ranks(nidx(1)%beg:nidx(1)%end, nidx(2)%beg:nidx(2)%end, nidx(3)%beg:nidx(3)%end))
        do k = nidx(3)%beg, nidx(3)%end
            do j = nidx(2)%beg, nidx(2)%end
                do i = nidx(1)%beg, nidx(1)%end
                    if (abs(i) + abs(j) + abs(k) > 0) then
                        neighbor_coords(1) = proc_coords(1) + i
                        if (num_dims > 1) neighbor_coords(2) = proc_coords(2) + j
                        if (num_dims > 2) neighbor_coords(3) = proc_coords(3) + k
                        call MPI_CART_RANK(MPI_COMM_CART, neighbor_coords, neighbor_ranks(i, j, k), ierr)
                    end if
                end do
            end do
        end do
#endif

    end subroutine s_mpi_decompose_computational_domain

    !> Apply executable-configured output and local-domain policies after the shared Cartesian decomposition.
    subroutine s_apply_decomposition_policies(proc_counts, remainders, local_cells, global_cells, write_silo_ghost_offsets, &
        & adjust_local_domains, output_offsets, local_domains)

        integer, dimension(3), intent(in)                            :: proc_counts, remainders, local_cells, global_cells
        logical, intent(in)                                          :: write_silo_ghost_offsets, adjust_local_domains
        type(int_bounds_info), dimension(3), intent(inout), optional :: output_offsets
        type(bounds_info), dimension(3), intent(inout), optional     :: local_domains
        integer                                                      :: dim
        real(wp)                                                     :: domain_beg, domain_end, spacing

        if (write_silo_ghost_offsets .and. .not. present(output_offsets)) then
            call s_mpi_abort('Silo ghost-offset policy requires output offset storage.')
        end if
        if (adjust_local_domains .and. .not. present(local_domains)) then
            call s_mpi_abort('Local-domain adjustment policy requires domain storage.')
        end if

        if (present(output_offsets)) then
            do dim = 1, 3
                output_offsets(dim)%beg = 0
                output_offsets(dim)%end = 0
            end do

            if (write_silo_ghost_offsets) then
                do dim = 1, num_dims
                    if (proc_coords(dim) > 0) output_offsets(dim)%beg = 2
                    if (proc_coords(dim) < proc_counts(dim) - 1) output_offsets(dim)%end = 2
                end do
            end if
        end if

        if (adjust_local_domains) then
            do dim = 1, num_dims
                domain_beg = local_domains(dim)%beg
                domain_end = local_domains(dim)%end
                spacing = (domain_end - domain_beg)/real(global_cells(dim) + 1, wp)

                select case (dim)
                case (1)
                    dx_min = spacing
                case (2)
                    dy_min = spacing
                case (3)
                    dz_min = spacing
                end select

                if (proc_coords(dim) < remainders(dim)) then
                    local_domains(dim)%beg = domain_beg + spacing*real((local_cells(dim) + 1)*proc_coords(dim), wp)
                    local_domains(dim)%end = domain_end - spacing*real((local_cells(dim) + 1)*(proc_counts(dim) - proc_coords(dim) &
                                  & - 1) - (proc_counts(dim) - remainders(dim)), wp)
                else
                    local_domains(dim)%beg = domain_beg + spacing*real((local_cells(dim) + 1)*proc_coords(dim) + remainders(dim), &
                                  & wp)
                    local_domains(dim)%end = domain_end - spacing*real((local_cells(dim) + 1)*(proc_counts(dim) - proc_coords(dim) &
                                  & - 1), wp)
                end if
            end do
        end if

    end subroutine s_apply_decomposition_policies

    !> The goal of this procedure is to populate the buffers of the grid variables by communicating with the neighboring processors.
    !! Note that only the buffers of the cell-width distributions are handled in such a way. This is because the buffers of
    !! cell-boundary locations may be calculated directly from those of the cell-width distributions.
    subroutine s_mpi_sendrecv_grid_variable_buffer(cell_boundaries, cell_centers, cell_widths, num_cells, bc_bounds, pbc_loc, &
        & offset)

        integer, intent(in)               :: num_cells, pbc_loc
        type(int_bounds_info), intent(in) :: bc_bounds, offset
        ! Contiguous so that passing an element to MPI is a plain address, with no descriptor
        ! or copy-in/copy-out. Every actual argument is a whole module array. The attribute must
        ! be on ALL THREE, and on every frame that forwards them: CCE 19 IPA drops the stores to
        ! the non-contiguous dummies when one call mixes contiguous and non-contiguous arrays.
        real(wp), contiguous, intent(inout) :: cell_boundaries(-1 - offset%beg:)
        real(wp), contiguous, intent(inout) :: cell_centers(-buff_size:)
        real(wp), contiguous, intent(inout) :: cell_widths(-buff_size:)

#ifdef MFC_MPI
        integer :: ierr
        integer :: i

        if (pbc_loc == -1) then
            if (bc_bounds%end >= 0) then
                call MPI_SENDRECV(cell_widths(num_cells - buff_size + 1), buff_size, mpi_p, bc_bounds%end, 0, &
                                  & cell_widths(-buff_size), buff_size, mpi_p, bc_bounds%beg, 0, MPI_COMM_WORLD, &
                                  & MPI_STATUS_IGNORE, ierr)
            else
                call MPI_SENDRECV(cell_widths(0), buff_size, mpi_p, bc_bounds%beg, 1, cell_widths(-buff_size), buff_size, mpi_p, &
                                  & bc_bounds%beg, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr)
            end if
            do i = 1, offset%beg
                cell_boundaries(-1 - i) = cell_boundaries(-i) - cell_widths(-i)
            end do
            do i = 1, buff_size
                cell_centers(-i) = cell_centers(1 - i) - (cell_widths(1 - i) + cell_widths(-i))/2._wp
            end do
        else
            if (bc_bounds%beg >= 0) then
                call MPI_SENDRECV(cell_widths(0), buff_size, mpi_p, bc_bounds%beg, 1, cell_widths(num_cells + 1), buff_size, &
                                  & mpi_p, bc_bounds%end, 1, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr)
            else
                call MPI_SENDRECV(cell_widths(num_cells - buff_size + 1), buff_size, mpi_p, bc_bounds%end, 0, &
                                  & cell_widths(num_cells + 1), buff_size, mpi_p, bc_bounds%end, 1, MPI_COMM_WORLD, &
                                  & MPI_STATUS_IGNORE, ierr)
            end if
            do i = 1, offset%end
                cell_boundaries(num_cells + i) = cell_boundaries(num_cells + i - 1) + cell_widths(num_cells + i)
            end do
            do i = 1, buff_size
                cell_centers(num_cells + i) = cell_centers(num_cells + i - 1) + (cell_widths(num_cells + i - 1) &
                             & + cell_widths(num_cells + i))/2._wp
            end do
        end if
#endif

    end subroutine s_mpi_sendrecv_grid_variable_buffer

    !> Populate the local cell-boundary, cell-center, and cell-width arrays in one direction directly from the global cell-boundary
    !! array. This guarantees that every rank sees bitwise-identical values at any shared physical cell or boundary
    subroutine s_apply_grid_from_global_dim(x_cb_glb, m_dim_glb, m_dim, sidx, bc_beg, bc_end, cb_lo, cb_hi, cw_lo, cw_hi, &
                                            & x_cb_loc, x_cc_loc, dx_loc)

        integer, intent(in)     :: m_dim_glb, m_dim, sidx, bc_beg, bc_end
        integer, intent(in)     :: cb_lo, cb_hi, cw_lo, cw_hi
        real(wp), intent(in)    :: x_cb_glb(-1:m_dim_glb)
        real(wp), intent(inout) :: x_cb_loc(-1 - cb_lo:m_dim + cb_hi)
        real(wp), intent(inout) :: x_cc_loc(-cw_lo:m_dim + cw_hi)
        real(wp), intent(inout) :: dx_loc(-cw_lo:m_dim + cw_hi)
        real(wp)                :: domain_len
        integer                 :: i, gidx, lo, hi

        domain_len = x_cb_glb(m_dim_glb) - x_cb_glb(-1)

        ! Interior cell boundaries sliced directly from the global list
        do i = -1, m_dim
            x_cb_loc(i) = x_cb_glb(sidx + i)
        end do

        ! Left ghost cell boundaries
        if (bc_beg >= 0) then
            if (sidx == 0) then
                ! Leftmost rank with a neighbor -> periodic+multirank, so wrap from the global right end
                do i = 1, cb_lo
                    x_cb_loc(-1 - i) = x_cb_glb(m_dim_glb - i) - domain_len
                end do
            else
                do i = 1, cb_lo
                    gidx = sidx - 1 - i
                    if (gidx >= -1) then
                        x_cb_loc(-1 - i) = x_cb_glb(gidx)
                    else
                        x_cb_loc(-1 - i) = x_cb_glb(m_dim_glb + 1 + gidx) - domain_len
                    end if
                end do
            end if
        end if

        ! Right ghost cell boundaries
        if (bc_end >= 0) then
            if (sidx + m_dim == m_dim_glb) then
                ! Rightmost rank with a neighbor -> periodic+multirank, wrap from the global left end
                do i = 1, cb_hi
                    x_cb_loc(m_dim + i) = x_cb_glb(i - 1) + domain_len
                end do
            else
                do i = 1, cb_hi
                    gidx = sidx + m_dim + i
                    if (gidx <= m_dim_glb) then
                        x_cb_loc(m_dim + i) = x_cb_glb(gidx)
                    else
                        x_cb_loc(m_dim + i) = x_cb_glb(gidx - m_dim_glb - 1) + domain_len
                    end if
                end do
            end if
        end if

        ! Recompute dx and x_cc over the range where x_cb is now valid using one formula so values are bitwise-identical
        if (bc_beg >= 0) then
            lo = -min(cw_lo, cb_lo)
        else
            lo = 0
        end if

        if (bc_end >= 0) then
            hi = m_dim + min(cw_hi, cb_hi)
        else
            hi = m_dim
        end if

        do i = lo, hi
            dx_loc(i) = x_cb_loc(i) - x_cb_loc(i - 1)
            x_cc_loc(i) = (x_cb_loc(i) + x_cb_loc(i - 1))/2._wp
        end do

    end subroutine s_apply_grid_from_global_dim

    !> Module deallocation and/or disassociation procedures
    impure subroutine s_finalize_mpi_common_module

#ifdef MFC_MPI
#ifndef __NVCOMPILER_GPU_UNIFIED_MEM
        @:DEALLOCATE(buff_send, buff_recv)
#else
        $:GPU_EXIT_DATA(delete='[buff_send, buff_recv]')
        deallocate (buff_send, buff_recv)
#endif
#endif

    end subroutine s_finalize_mpi_common_module

end module m_mpi_common