m_start_up.fpp

!>
!! @file
!! @brief Contains module m_start_up

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

!> @brief Reads input files, loads initial conditions and grid data, and orchestrates solver initialization and finalization
module m_start_up

    use m_derived_types
    use m_global_parameters
    use m_mpi_proxy
    use m_mpi_common
    use m_variables_conversion
    use m_weno
    use m_muscl
    use m_thinc
    use m_riemann_solvers
    use m_cbc
    use m_boundary_common
    use m_boundary_io
    use m_acoustic_src
    use m_rhs
    use m_pressure_relaxation, only: s_report_pressure_relaxation
    use m_chemistry
    use m_data_output
    use m_time_steppers
    use m_qbmm
    use m_derived_variables
    use m_hypoelastic
    use m_phase_change
    use m_viscous
    use m_bubbles_EE
    use m_bubbles_EL
    use ieee_arithmetic
    use m_helper_basic
    use m_helper

    $:USE_GPU_MODULE()

    use m_nvtx
    use m_ibm
    use m_ib_patches
    use m_model
    use m_particle_cloud
    use m_collisions
    use m_compile_specific
    use m_checker_common
    use m_checker
    use m_surface_tension
    use m_body_forces
    use m_sim_helpers
    use m_igr
    use m_constants, only: model_eqns_6eq, time_stepper_rk1, time_stepper_rk2, time_stepper_rk3, recon_type_weno, recon_type_muscl

    implicit none

    private; public :: s_read_input_file, s_check_input_file, s_read_data_files, s_read_serial_data_files, &
        & s_read_parallel_data_files, s_initialize_internal_energy_equations, s_initialize_modules, s_initialize_gpu_vars, &
        & s_initialize_mpi_domain, s_finalize_modules, s_perform_time_step, s_save_data, s_save_performance_metrics

    type(scalar_field), allocatable, dimension(:) :: q_cons_temp
    real(wp)                                      :: dt_init

contains

    !> Read data files. Dispatch subroutine that replaces procedure pointer.
    impure subroutine s_read_data_files(q_cons_vf)

        type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf

        if (.not. parallel_io) then
            call s_read_serial_data_files(q_cons_vf)
        else
            call s_read_parallel_data_files(q_cons_vf)
        end if

    end subroutine s_read_data_files

    !> Verify the input file exists and read it
    impure subroutine s_read_input_file

        character(LEN=name_len), parameter :: file_path = './simulation.inp'
        logical                            :: file_exist  !< Logical used to check the existence of the input file
        integer                            :: iostatus
        ! Integer to check iostat of file read

        character(len=1000) :: line

        #:include 'generated_namelist.fpp'

        inquire (FILE=trim(file_path), EXIST=file_exist)

        if (file_exist) then
            open (1, FILE=trim(file_path), form='formatted', ACTION='read', STATUS='old')
            read (1, NML=user_inputs, iostat=iostatus)

            if (iostatus /= 0) then
                backspace (1)
                read (1, fmt='(A)') line
                print *, 'Invalid line in namelist: ' // trim(line)
                call s_mpi_abort('Invalid line in simulation.inp. It is ' // 'likely due to a datatype mismatch. Exiting.')
            end if

            close (1)

            if ((bf_x) .or. (bf_y) .or. (bf_z) .or. (bf_spatial_support)) then
                bodyForces = .true.
            end if

            m_glb = m
            n_glb = n
            p_glb = p

            call s_update_cell_bounds(cells_bounds, m, n, p)

            if (cfl_adap_dt .or. cfl_const_dt) cfl_dt = .true.

            if (any((/bc_x%beg, bc_x%end, bc_y%beg, bc_y%end, bc_z%beg, bc_z%end/) == -17) .or. num_bc_patches > 0) then
                bc_io = .true.
            end if

            if (bc_x%beg == BC_PERIODIC .and. bc_x%end == BC_PERIODIC) periodic_bc(1) = .true.
            if (bc_y%beg == BC_PERIODIC .and. bc_y%end == BC_PERIODIC) periodic_bc(2) = .true.
            if (bc_z%beg == BC_PERIODIC .and. bc_z%end == BC_PERIODIC) periodic_bc(3) = .true.
        else
            call s_mpi_abort(trim(file_path) // ' is missing. Exiting.')
        end if

    end subroutine s_read_input_file

    !> Validate that all user-provided inputs form a consistent simulation configuration
    impure subroutine s_check_input_file

        character(LEN=path_len) :: file_path
        logical                 :: file_exist

        file_path = trim(case_dir) // '/.'

        call my_inquire(file_path, file_exist)

        if (file_exist .neqv. .true.) then
            call s_mpi_abort(trim(file_path) // ' is missing. Exiting.')
        end if

        call s_check_inputs_common(check_total_cells=.false., n_global=0_8)
        call s_check_inputs()

    end subroutine s_check_input_file

    !> Read serial initial condition and grid data files and compute cell-width distributions
    impure subroutine s_read_serial_data_files(q_cons_vf)

        type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
        character(LEN=path_len + 2*name_len) :: t_step_dir  !< Relative path to the starting time-step directory
        character(LEN=path_len + 3*name_len) :: file_path   !< Relative path to the grid and conservative variables data files
        logical :: file_exist                               !< Logical used to check the existence of the input file
        integer :: i, r

        if (cfl_dt) then
            write (t_step_dir, '(A,I0,A,I0)') trim(case_dir) // '/p_all/p', proc_rank, '/', n_start
        else
            write (t_step_dir, '(A,I0,A,I0)') trim(case_dir) // '/p_all/p', proc_rank, '/', t_step_start
        end if

        file_path = trim(t_step_dir) // '/.'
        call my_inquire(file_path, file_exist)

        if (file_exist .neqv. .true.) then
            call s_mpi_abort(trim(file_path) // ' is missing. Exiting.')
        end if

        if (bc_io) then
            call s_read_serial_boundary_condition_files(t_step_dir, bc_type)
        else
            call s_assign_default_bc_type(bc_type)
        end if

        file_path = trim(t_step_dir) // '/x_cb.dat'

        inquire (FILE=trim(file_path), EXIST=file_exist)

        if (file_exist) then
            open (2, FILE=trim(file_path), form='unformatted', ACTION='read', STATUS='old')
            read (2) x_cb(-1:m); close (2)
        else
            call s_mpi_abort(trim(file_path) // ' is missing. Exiting.')
        end if

        dx(0:m) = x_cb(0:m) - x_cb(-1:m - 1)
        x_cc(0:m) = x_cb(-1:m - 1) + dx(0:m)/2._wp

        if (n > 0) then
            file_path = trim(t_step_dir) // '/y_cb.dat'

            inquire (FILE=trim(file_path), EXIST=file_exist)

            if (file_exist) then
                open (2, FILE=trim(file_path), form='unformatted', ACTION='read', STATUS='old')
                read (2) y_cb(-1:n); close (2)
            else
                call s_mpi_abort(trim(file_path) // ' is missing. Exiting.')
            end if

            dy(0:n) = y_cb(0:n) - y_cb(-1:n - 1)
            y_cc(0:n) = y_cb(-1:n - 1) + dy(0:n)/2._wp
        end if

        if (p > 0) then
            file_path = trim(t_step_dir) // '/z_cb.dat'

            inquire (FILE=trim(file_path), EXIST=file_exist)

            if (file_exist) then
                open (2, FILE=trim(file_path), form='unformatted', ACTION='read', STATUS='old')
                read (2) z_cb(-1:p); close (2)
            else
                call s_mpi_abort(trim(file_path) // ' is missing. Exiting.')
            end if

            dz(0:p) = z_cb(0:p) - z_cb(-1:p - 1)
            z_cc(0:p) = z_cb(-1:p - 1) + dz(0:p)/2._wp
        end if

        do i = 1, sys_size
            write (file_path, '(A,I0,A)') trim(t_step_dir) // '/q_cons_vf', i, '.dat'
            inquire (FILE=trim(file_path), EXIST=file_exist)
            if (file_exist) then
                open (2, FILE=trim(file_path), form='unformatted', ACTION='read', STATUS='old')
                read (2) q_cons_vf(i)%sf(0:m,0:n,0:p); close (2)
            else
                call s_mpi_abort(trim(file_path) // ' is missing. Exiting.')
            end if
        end do

        if (bubbles_euler .or. hypoelasticity) then
            ! Read pb and mv for non-polytropic qbmm
            if (qbmm .and. .not. polytropic) then
                do i = 1, nb
                    do r = 1, nnode
                        write (file_path, '(A,I0,A)') trim(t_step_dir) // '/pb', sys_size + (i - 1)*nnode + r, '.dat'
                        inquire (FILE=trim(file_path), EXIST=file_exist)
                        if (file_exist) then
                            open (2, FILE=trim(file_path), form='unformatted', ACTION='read', STATUS='old')
                            read (2) pb_ts(1)%sf(0:m,0:n,0:p,r, i); close (2)
                        else
                            call s_mpi_abort(trim(file_path) // ' is missing. Exiting.')
                        end if
                    end do
                end do
                do i = 1, nb
                    do r = 1, nnode
                        write (file_path, '(A,I0,A)') trim(t_step_dir) // '/mv', sys_size + (i - 1)*nnode + r, '.dat'
                        inquire (FILE=trim(file_path), EXIST=file_exist)
                        if (file_exist) then
                            open (2, FILE=trim(file_path), form='unformatted', ACTION='read', STATUS='old')
                            read (2) mv_ts(1)%sf(0:m,0:n,0:p,r, i); close (2)
                        else
                            call s_mpi_abort(trim(file_path) // ' is missing. Exiting.')
                        end if
                    end do
                end do
            end if
        end if

    end subroutine s_read_serial_data_files

    !> Read parallel initial condition and grid data files via MPI I/O
    impure subroutine s_read_parallel_data_files(q_cons_vf)

        type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf

#ifdef MFC_MPI
        real(wp), allocatable, dimension(:)  :: x_cb_glb, y_cb_glb, z_cb_glb
        integer                              :: ifile, ierr, data_size
        integer, dimension(MPI_STATUS_SIZE)  :: status
        integer(KIND=MPI_OFFSET_KIND)        :: disp
        integer(KIND=MPI_OFFSET_KIND)        :: m_MOK, n_MOK, p_MOK
        integer(KIND=MPI_OFFSET_KIND)        :: WP_MOK, var_MOK
        integer(KIND=MPI_OFFSET_KIND)        :: MOK
        character(LEN=path_len + 2*name_len) :: file_loc
        logical                              :: file_exist
        character(len=10)                    :: t_step_start_string
        integer                              :: i, j

        ! Downsampled data variables
        integer :: m_ds, n_ds, p_ds
        integer :: m_glb_ds, n_glb_ds, p_glb_ds
        integer :: m_glb_read, n_glb_read, p_glb_read  !< data size of read

        allocate (x_cb_glb(-1:m_glb))
        allocate (y_cb_glb(-1:n_glb))
        allocate (z_cb_glb(-1:p_glb))

        file_loc = trim(case_dir) // '/restart_data' // trim(mpiiofs) // 'x_cb.dat'
        inquire (FILE=trim(file_loc), EXIST=file_exist)

        if (down_sample) then
            m_ds = int((m + 1)/3) - 1
            n_ds = int((n + 1)/3) - 1
            p_ds = int((p + 1)/3) - 1

            m_glb_ds = int((m_glb + 1)/3) - 1
            n_glb_ds = int((n_glb + 1)/3) - 1
            p_glb_ds = int((p_glb + 1)/3) - 1
        end if

        if (file_exist) then
            data_size = m_glb + 2
            call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr)
            call MPI_FILE_READ(ifile, x_cb_glb, data_size, mpi_p, status, ierr)
            call MPI_FILE_CLOSE(ifile, ierr)
        else
            call s_mpi_abort('File ' // trim(file_loc) // ' is missing. Exiting.')
        end if

        call s_apply_grid_from_global_dim(x_cb_glb, m_glb, m, start_idx(1), bc_x%beg, bc_x%end, buff_size, buff_size, buff_size, &
                                          & buff_size, x_cb, x_cc, dx)

        if (n > 0) then
            file_loc = trim(case_dir) // '/restart_data' // trim(mpiiofs) // 'y_cb.dat'
            inquire (FILE=trim(file_loc), EXIST=file_exist)

            if (file_exist) then
                data_size = n_glb + 2
                call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr)
                call MPI_FILE_READ(ifile, y_cb_glb, data_size, mpi_p, status, ierr)
                call MPI_FILE_CLOSE(ifile, ierr)
            else
                call s_mpi_abort('File ' // trim(file_loc) // ' is missing. Exiting.')
            end if

            call s_apply_grid_from_global_dim(y_cb_glb, n_glb, n, start_idx(2), bc_y%beg, bc_y%end, buff_size, buff_size, &
                                              & buff_size, buff_size, y_cb, y_cc, dy)

            if (p > 0) then
                file_loc = trim(case_dir) // '/restart_data' // trim(mpiiofs) // 'z_cb.dat'
                inquire (FILE=trim(file_loc), EXIST=file_exist)

                if (file_exist) then
                    data_size = p_glb + 2
                    call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr)
                    call MPI_FILE_READ(ifile, z_cb_glb, data_size, mpi_p, status, ierr)
                    call MPI_FILE_CLOSE(ifile, ierr)
                else
                    call s_mpi_abort('File ' // trim(file_loc) // 'is missing. Exiting.')
                end if

                call s_apply_grid_from_global_dim(z_cb_glb, p_glb, p, start_idx(3), bc_z%beg, bc_z%end, buff_size, buff_size, &
                                                  & buff_size, buff_size, z_cb, z_cc, dz)
            end if
        end if

        if (file_per_process) then
            if (cfl_dt) then
                call s_int_to_str(n_start, t_step_start_string)
                write (file_loc, '(I0,A1,I7.7,A)') n_start, '_', proc_rank, '.dat'
            else
                call s_int_to_str(t_step_start, t_step_start_string)
                write (file_loc, '(I0,A1,I7.7,A)') t_step_start, '_', proc_rank, '.dat'
            end if
            file_loc = trim(case_dir) // '/restart_data/lustre_' // trim(t_step_start_string) // trim(mpiiofs) // trim(file_loc)
            inquire (FILE=trim(file_loc), EXIST=file_exist)

            if (file_exist) then
                call MPI_FILE_OPEN(MPI_COMM_SELF, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr)

                if (down_sample) then
                    call s_initialize_mpi_data_ds(m_ds, n_ds, p_ds)
                else
                    if (ib) then
                        call s_initialize_mpi_data(q_cons_vf, ib_markers=ib_markers, ib_mpi_data=MPI_IO_IB_DATA, &
                                                   & qbmm_pb=pb_ts(1), qbmm_mv=mv_ts(1))
                    else
                        call s_initialize_mpi_data(q_cons_vf, qbmm_pb=pb_ts(1), qbmm_mv=mv_ts(1))
                    end if
                end if

                if (down_sample) then
                    data_size = (m_ds + 3)*(n_ds + 3)*(p_ds + 3)
                    m_glb_read = m_glb_ds + 1
                    n_glb_read = n_glb_ds + 1
                    p_glb_read = p_glb_ds + 1
                else
                    data_size = (m + 1)*(n + 1)*(p + 1)
                    m_glb_read = m_glb + 1
                    n_glb_read = n_glb + 1
                    p_glb_read = p_glb + 1
                end if

                m_MOK = int(m_glb_read + 1, MPI_OFFSET_KIND)
                n_MOK = int(m_glb_read + 1, MPI_OFFSET_KIND)
                p_MOK = int(m_glb_read + 1, MPI_OFFSET_KIND)
                WP_MOK = int(storage_size(0._stp)/8, MPI_OFFSET_KIND)
                MOK = int(1._wp, MPI_OFFSET_KIND)

                if (bubbles_euler .or. hypoelasticity) then
                    do i = 1, sys_size
                        var_MOK = int(i, MPI_OFFSET_KIND)

                        call MPI_FILE_READ(ifile, MPI_IO_DATA%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
                    end do
                    ! Read pb and mv for non-polytropic qbmm
                    if (qbmm .and. .not. polytropic) then
                        do i = sys_size + 1, sys_size + 2*nb*nnode
                            var_MOK = int(i, MPI_OFFSET_KIND)

                            call MPI_FILE_READ(ifile, MPI_IO_DATA%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
                        end do
                    end if
                else
                    if (down_sample) then
                        do i = 1, sys_size
                            var_MOK = int(i, MPI_OFFSET_KIND)

                            call MPI_FILE_READ(ifile, q_cons_temp(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
                        end do
                    else
                        do i = 1, sys_size
                            var_MOK = int(i, MPI_OFFSET_KIND)

                            call MPI_FILE_READ(ifile, MPI_IO_DATA%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
                        end do
                    end if
                end if

                call s_mpi_barrier()

                call MPI_FILE_CLOSE(ifile, ierr)
            else
                call s_mpi_abort('File ' // trim(file_loc) // ' is missing. Exiting.')
            end if
        else
            if (cfl_dt) then
                write (file_loc, '(I0,A)') n_start, '.dat'
            else
                write (file_loc, '(I0,A)') t_step_start, '.dat'
            end if
            file_loc = trim(case_dir) // '/restart_data' // trim(mpiiofs) // trim(file_loc)
            inquire (FILE=trim(file_loc), EXIST=file_exist)

            if (file_exist) then
                call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr)

                if (ib) then
                    call s_initialize_mpi_data(q_cons_vf, ib_markers=ib_markers, ib_mpi_data=MPI_IO_IB_DATA, qbmm_pb=pb_ts(1), &
                                               & qbmm_mv=mv_ts(1))
                else
                    call s_initialize_mpi_data(q_cons_vf, qbmm_pb=pb_ts(1), qbmm_mv=mv_ts(1))
                end if

                data_size = (m + 1)*(n + 1)*(p + 1)

                m_MOK = int(m_glb + 1, MPI_OFFSET_KIND)
                n_MOK = int(n_glb + 1, MPI_OFFSET_KIND)
                p_MOK = int(p_glb + 1, MPI_OFFSET_KIND)
                WP_MOK = int(storage_size(0._stp)/8, MPI_OFFSET_KIND)
                MOK = int(1._wp, MPI_OFFSET_KIND)

                if (bubbles_euler .or. hypoelasticity) then
                    do i = 1, sys_size
                        var_MOK = int(i, MPI_OFFSET_KIND)
                        disp = m_MOK*max(MOK, n_MOK)*max(MOK, p_MOK)*WP_MOK*(var_MOK - 1)

                        call MPI_FILE_SET_VIEW(ifile, disp, mpi_io_p, MPI_IO_DATA%view(i), 'native', mpi_info_int, ierr)
                        call MPI_FILE_READ(ifile, MPI_IO_DATA%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
                    end do
                    ! Read pb and mv for non-polytropic qbmm
                    if (qbmm .and. .not. polytropic) then
                        do i = sys_size + 1, sys_size + 2*nb*nnode
                            var_MOK = int(i, MPI_OFFSET_KIND)
                            disp = m_MOK*max(MOK, n_MOK)*max(MOK, p_MOK)*WP_MOK*(var_MOK - 1)

                            call MPI_FILE_SET_VIEW(ifile, disp, mpi_io_p, MPI_IO_DATA%view(i), 'native', mpi_info_int, ierr)
                            call MPI_FILE_READ(ifile, MPI_IO_DATA%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
                        end do
                    end if
                else
                    do i = 1, sys_size
                        var_MOK = int(i, MPI_OFFSET_KIND)

                        disp = m_MOK*max(MOK, n_MOK)*max(MOK, p_MOK)*WP_MOK*(var_MOK - 1)

                        call MPI_FILE_SET_VIEW(ifile, disp, mpi_io_p, MPI_IO_DATA%view(i), 'native', mpi_info_int, ierr)
                        call MPI_FILE_READ_ALL(ifile, MPI_IO_DATA%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
                    end do
                end if

                call s_mpi_barrier()

                call MPI_FILE_CLOSE(ifile, ierr)
            else
                call s_mpi_abort('File ' // trim(file_loc) // ' is missing. Exiting.')
            end if
        end if

        deallocate (x_cb_glb, y_cb_glb, z_cb_glb)

        if (bc_io) then
            call s_read_parallel_boundary_condition_files(bc_type)
        else
            call s_assign_default_bc_type(bc_type)
        end if
#endif

    end subroutine s_read_parallel_data_files

    !> Initialize internal-energy equations from phase mass, mixture momentum, and total energy
    subroutine s_initialize_internal_energy_equations(v_vf)

        type(scalar_field), dimension(sys_size), intent(inout) :: v_vf
        real(wp)                                               :: rho
        real(wp)                                               :: dyn_pres
        real(wp)                                               :: gamma
        real(wp)                                               :: pi_inf
        real(wp)                                               :: qv
        real(wp), dimension(2)                                 :: Re
        real(wp)                                               :: pres, T
        real(wp)                                               :: alpha_i, alpha_rho_i, e_i
        integer                                                :: i, j, k, l, c
        real(wp), dimension(num_species)                       :: rhoYks
        real(wp)                                               :: pres_mag

        pres_mag = 0._wp

        T = dflt_T_guess

        do j = 0, m
            do k = 0, n
                do l = 0, p
                    call s_convert_to_mixture_variables(v_vf, j, k, l, rho, gamma, pi_inf, qv, Re)

                    dyn_pres = 0._wp
                    do i = eqn_idx%mom%beg, eqn_idx%mom%end
                        dyn_pres = dyn_pres + 5.e-1_wp*v_vf(i)%sf(j, k, l)*v_vf(i)%sf(j, k, l)/max(rho, sgm_eps)
                    end do

                    if (chemistry) then
                        do c = 1, num_species
                            rhoYks(c) = v_vf(eqn_idx%species%beg + c - 1)%sf(j, k, l)
                        end do
                    end if

                    if (mhd) then
                        if (n == 0) then
                            pres_mag = 0.5_wp*(Bx0**2 + v_vf(eqn_idx%B%beg)%sf(j, k, l)**2 + v_vf(eqn_idx%B%beg + 1)%sf(j, k, l)**2)
                        else
                            pres_mag = 0.5_wp*(v_vf(eqn_idx%B%beg)%sf(j, k, l)**2 + v_vf(eqn_idx%B%beg + 1)%sf(j, k, &
                                               & l)**2 + v_vf(eqn_idx%B%beg + 2)%sf(j, k, l)**2)
                        end if
                    end if

                    call s_compute_pressure(v_vf(eqn_idx%E)%sf(j, k, l), 0._stp, dyn_pres, pi_inf, gamma, rho, qv, rhoYks, pres, &
                                            & T, pres_mag=pres_mag)

                    do i = 1, num_fluids
                        alpha_i = v_vf(i + eqn_idx%adv%beg - 1)%sf(j, k, l)
                        alpha_rho_i = v_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)
                        call s_phase_internal_energy(pres, alpha_i, alpha_rho_i, i, e_i)
                        v_vf(i + eqn_idx%int_en%beg - 1)%sf(j, k, l) = e_i
                    end do
                end do
            end do
        end do

    end subroutine s_initialize_internal_energy_equations

    !> Advance the simulation by one time step, handling CFL-based dt and time-stepper dispatch
    impure subroutine s_perform_time_step(t_step, time_avg)

        integer, intent(inout)  :: t_step
        real(wp), intent(inout) :: time_avg
        integer                 :: i, eta_hh, eta_mm, eta_ss
        real(wp)                :: eta_sec
        real(wp)                :: dt_floor
        character(len=8)        :: lim_str  !< Time-step limiter tag, e.g. ' (ICFL)'

        if (cfl_dt) then
            if (cfl_const_dt .and. t_step == 0) call s_compute_dt()

            if (cfl_adap_dt) call s_compute_dt()

            if (t_step == 0) dt_init = dt

            ! the collision restriction deliberately drops dt to collision_time/collision_temporal_resolution, so lower the
            ! runaway-dt abort threshold below that cap when it is enabled
            dt_floor = 1.e-3_wp*dt_init
            if (collision_model > 0 .and. collision_temporal_resolution > 0) then
                dt_floor = min(dt_floor, 1.e-3_wp*collision_time/real(collision_temporal_resolution, wp))
            end if

            if (dt < dt_floor .and. cfl_adap_dt .and. proc_rank == 0) then
                print *, "Delta t = ", dt
                call s_mpi_abort("Delta t has become too small")
            end if
        end if

        if (cfl_dt) then
            if ((mytime + dt) >= t_stop) then
                dt = t_stop - mytime
                $:GPU_UPDATE(device='[dt]')
            end if
        else
            if ((mytime + dt) >= finaltime) then
                dt = finaltime - mytime
                $:GPU_UPDATE(device='[dt]')
            end if
        end if

        if (cfl_dt) then
            if (proc_rank == 0 .and. mod(t_step - t_step_start, t_step_print) == 0) then
                eta_sec = wall_time_avg*(t_stop - mytime)/max(dt, tiny(dt))
                eta_hh = int(eta_sec)/3600
                eta_mm = mod(int(eta_sec), 3600)/60
                eta_ss = mod(int(eta_sec), 60)
                lim_str = ''
                if (cfl_adap_dt) lim_str = ' (' // dt_limiter // ')'
                print '(" [", I3, "%] t = ", ES11.4, " dt = ", ES11.4, A, " @ step ", I0, " t/step ", ES9.2, "s (avg ", ES9.2, "s) ETA ", I0, ":", I2.2, ":", I2.2)', &
                    & int(ceiling(100._wp*(mytime/t_stop))), mytime, dt, trim(lim_str), t_step, wall_time, wall_time_avg, eta_hh, &
                    & eta_mm, eta_ss
            end if
        else
            if (proc_rank == 0 .and. mod(t_step - t_step_start, t_step_print) == 0) then
                eta_sec = wall_time_avg*real(t_step_stop - t_step, wp)
                eta_hh = int(eta_sec)/3600
                eta_mm = mod(int(eta_sec), 3600)/60
                eta_ss = mod(int(eta_sec), 60)
                print '(" [", I3, "%] step ", I0, " of ", I0, " (t_step ", I0, ") t/step ", ES9.2, "s (avg ", ES9.2, "s) ETA ", I0, ":", I2.2, ":", I2.2)', &
                    & int(ceiling(100._wp*(real(t_step - t_step_start)/(t_step_stop - t_step_start + 1)))), &
                    & t_step - t_step_start + 1, t_step_stop - t_step_start + 1, t_step, wall_time, wall_time_avg, eta_hh, &
                    & eta_mm, eta_ss
            end if
        end if

        if (probe_wrt) then
            do i = 1, sys_size
                $:GPU_UPDATE(host='[q_cons_ts(1)%vf(i)%sf]')
            end do
            if (bubbles_euler) then
                $:GPU_UPDATE(host='[ptil]')
            end if
        end if

        ! Total-variation-diminishing (TVD) Runge-Kutta (RK) time-steppers
        if (any(time_stepper == (/time_stepper_rk1, time_stepper_rk2, time_stepper_rk3/))) then
            call s_tvd_rk(t_step, time_avg, time_stepper)
        end if

        ! Advance time after RK so source terms see current-step time
        mytime = mytime + dt

        if (relax) call s_infinite_relaxation_k(q_cons_ts(1)%vf)

        ! Time-stepping loop controls
        t_step = t_step + 1

    end subroutine s_perform_time_step

    !> Collect per-process wall-clock times and write aggregate performance metrics to file
    impure subroutine s_save_performance_metrics(time_avg, time_final, io_time_avg, io_time_final, proc_time, io_proc_time, &
        & file_exists)

        real(wp), intent(inout)               :: time_avg, time_final
        real(wp), intent(inout)               :: io_time_avg, io_time_final
        real(wp), dimension(:), intent(inout) :: proc_time
        real(wp), dimension(:), intent(inout) :: io_proc_time
        logical, intent(inout)                :: file_exists
        real(wp)                              :: grind_time

        call s_mpi_barrier()

        if (num_procs > 1) then
            call mpi_bcast_time_step_values(proc_time, time_avg)

            call mpi_bcast_time_step_values(io_proc_time, io_time_avg)
        end if

        if (proc_rank == 0) then
            time_final = 0._wp
            io_time_final = 0._wp
            if (num_procs == 1) then
                time_final = time_avg
                io_time_final = io_time_avg
            else
                time_final = maxval(proc_time)
                io_time_final = maxval(io_proc_time)
            end if

            grind_time = time_final*1.0e9_wp/(real(sys_size, wp)*real(maxval((/1, m_glb/)), wp)*real(maxval((/1, n_glb/)), &
                                              & wp)*real(maxval((/1, p_glb/)), wp))

            print *, "Performance:", grind_time, "ns/gp/eq/rhs"
            inquire (FILE='time_data.dat', EXIST=file_exists)
            if (file_exists) then
                open (1, file='time_data.dat', position='append', status='old')
            else
                open (1, file='time_data.dat', status='new')
                write (1, '(A10, A15, A15)') "Ranks", "s/step", "ns/gp/eq/rhs"
            end if

            write (1, '(I10, 2(F15.8))') num_procs, time_final, grind_time

            close (1)

            inquire (FILE='io_time_data.dat', EXIST=file_exists)
            if (file_exists) then
                open (1, file='io_time_data.dat', position='append', status='old')
            else
                open (1, file='io_time_data.dat', status='new')
                write (1, '(A10, A15)') "Ranks", "s/step"
            end if

            write (1, '(I10, F15.8)') num_procs, io_time_final
            close (1)
        end if

    end subroutine s_save_performance_metrics

    !> Save conservative variable data to disk at the current time step
    impure subroutine s_save_data(t_step, start, finish, io_time_avg, nt)

        integer, intent(inout)  :: t_step
        real(wp), intent(inout) :: start, finish, io_time_avg
        integer, intent(inout)  :: nt
        integer(kind=8)         :: i, j, k, l
        integer                 :: stor
        integer                 :: save_count

        if (down_sample) then
            call s_populate_variables_buffers(bc_type, q_cons_ts(1)%vf)
        end if

        stor = 1

        if (time_stepper /= time_stepper_rk1) then
            $:GPU_PARALLEL_LOOP(collapse=4, copyin='[idwbuff]')
            do i = 1, sys_size
                do l = idwbuff(3)%beg, idwbuff(3)%end
                    do k = idwbuff(2)%beg, idwbuff(2)%end
                        do j = idwbuff(1)%beg, idwbuff(1)%end
                            q_cons_ts(2)%vf(i)%sf(j, k, l) = q_cons_ts(1)%vf(i)%sf(j, k, l)
                        end do
                    end do
                end do
            end do
            $:END_GPU_PARALLEL_LOOP()
            stor = 2
        end if

        call cpu_time(start)
        call nvtxStartRange("SAVE-DATA")
        do i = 1, sys_size
#ifndef FRONTIER_UNIFIED
            $:GPU_UPDATE(host='[q_cons_ts(stor)%vf(i)%sf]')
#endif
            do l = 0, p
                do k = 0, n
                    do j = 0, m
                        if (ieee_is_nan(real(q_cons_ts(stor)%vf(i)%sf(j, k, l), kind=wp))) then
                            print *, "NaN(s) in timestep output.", j, k, l, i, proc_rank, t_step, m, n, p
                            call s_mpi_abort("NaN(s) in timestep output.")
                        end if
                    end do
                end do
            end do
        end do

        if (qbmm .and. .not. polytropic) then
            $:GPU_UPDATE(host='[pb_ts(1)%sf]')
            $:GPU_UPDATE(host='[mv_ts(1)%sf]')
        end if

        if (cfl_dt) then
            save_count = int(mytime/t_save)
        else
            save_count = t_step
        end if

        if (bubbles_lagrange) then
            $:GPU_UPDATE(host='[lag_id, mtn_pos, mtn_posPrev, mtn_vel, intfc_rad, intfc_vel, bub_R0, Rmax_stats, Rmin_stats, &
                         & bub_dphidt, gas_p, gas_mv, gas_mg, gas_betaT, gas_betaC]')
            do i = 1, n_el_bubs_loc
                if (ieee_is_nan(intfc_rad(i, 1)) .or. intfc_rad(i, 1) <= 0._wp) then
                    call s_mpi_abort("Bubble radius is negative or NaN, please reduce dt.")
                end if
            end do

            $:GPU_UPDATE(host='[q_beta(1)%sf]')
            call s_write_data_files(q_cons_ts(stor)%vf, q_T_sf, q_prim_vf, save_count, bc_type, q_beta(1))
            $:GPU_UPDATE(host='[Rmax_stats, Rmin_stats, gas_p, gas_mv, intfc_vel]')
            call s_write_restart_lag_bubbles(save_count)  ! parallel
            if (lag_params%write_bubbles_stats) call s_write_lag_bubble_stats()
        else
            call s_write_data_files(q_cons_ts(stor)%vf, q_T_sf, q_prim_vf, save_count, bc_type)
        end if

        ! Write IB kinematic state for restart
        if (ib) call s_write_ib_state_file(save_count)

        call nvtxEndRange
        call cpu_time(finish)
        if (cfl_dt) then
            nt = mytime/t_save
        else
            nt = int((t_step - t_step_start)/(t_step_save))
        end if

        if (nt == 1) then
            io_time_avg = abs(finish - start)
        else
            io_time_avg = (abs(finish - start) + io_time_avg*(nt - 1))/nt
        end if

    end subroutine s_save_data

    !> Initialize all simulation sub-modules in the required dependency order
    impure subroutine s_initialize_modules

        integer :: m_ds, n_ds, p_ds
        integer :: i

        call s_initialize_global_parameters_module()
        #:if USING_AMD
            #:for BC in [-5, -6, -7, -8, -9, -10, -11, -12, -13]
                @:PROHIBIT(any((/bc_x%beg, bc_x%end, bc_y%beg, bc_y%end, bc_z%beg, &
                           & bc_z%end/) == ${BC}$) .and. eqn_idx%adv%end > 20 .and. (.not. chemistry), &
                           & "CBC module with AMD compiler requires eqn_idx%adv%end <= 20 when case optimization is turned off")
                @:PROHIBIT(any((/bc_x%beg, bc_x%end, bc_y%beg, bc_y%end, bc_z%beg, &
                           & bc_z%end/) == ${BC}$) .and. sys_size > 20 .and. (chemistry), &
                           & "CBC module with AMD compiler and chemistry requires sys_size <= 20 when case optimization is turned off")
            #:endfor
        #:endif
        if (bubbles_euler .or. bubbles_lagrange) then
            call s_initialize_bubbles_model()
        end if
        call s_initialize_mpi_common_module(exchange_all_chemistry_temperatures_in=.false., use_rdma_transport_in=rdma_mpi)
        call s_initialize_mpi_proxy_module()
        call s_initialize_variables_conversion_module(enforce_density_floor=.true., preserve_qbmm_number=.true.)
        if (grid_geometry == 3) call s_initialize_fftw_module()

        if (bubbles_euler) call s_initialize_bubbles_EE_module()
        if (ib) then
            call s_initialize_ibm_module()
        end if
        if (qbmm) call s_initialize_qbmm_module()

        if (acoustic_source) then
            call s_initialize_acoustic_src()
        end if

        if (viscous .and. (.not. igr)) then
            call s_initialize_viscous_module()
        end if

        call s_initialize_rhs_module()

        if (surface_tension) call s_initialize_surface_tension_module()

        if (relax) call s_initialize_phasechange_module()

        call s_initialize_data_output_module()
        call s_initialize_derived_variables_module()
        call s_initialize_time_steppers_module()

        call s_initialize_boundary_common_module(use_dirichlet_buffers=.true.)

        if (down_sample) then
            m_ds = int((m + 1)/3) - 1
            n_ds = int((n + 1)/3) - 1
            p_ds = int((p + 1)/3) - 1

            allocate (q_cons_temp(1:sys_size))
            do i = 1, sys_size
                allocate (q_cons_temp(i)%sf(-1:m_ds + 1,-1:n_ds + 1,-1:p_ds + 1))
            end do
        end if

        if (down_sample) then
            call s_read_data_files(q_cons_temp)
            call s_upsample_data(q_cons_ts(1)%vf, q_cons_temp)
            do i = 1, sys_size
                $:GPU_UPDATE(device='[q_cons_ts(1)%vf(i)%sf]')
            end do
            do i = 1, sys_size
                deallocate (q_cons_temp(i)%sf)
            end do
            deallocate (q_cons_temp)
        else
            call s_read_data_files(q_cons_ts(1)%vf)
        end if

        block
            type(int_bounds_info), dimension(3) :: grid_offsets

            grid_offsets(:)%beg = buff_size
            grid_offsets(:)%end = buff_size
            if (n == 0) then
                call s_populate_grid_variables_buffers(x_cb, x_cc, dx, grid_offsets(1), grid_offsets(2), grid_offsets(3), &
                                                       & global_bounds=glb_bounds)
            else if (p == 0) then
                call s_populate_grid_variables_buffers(x_cb, x_cc, dx, grid_offsets(1), grid_offsets(2), grid_offsets(3), y_cb, &
                                                       & y_cc, dy, global_bounds=glb_bounds)
            else
                call s_populate_grid_variables_buffers(x_cb, x_cc, dx, grid_offsets(1), grid_offsets(2), grid_offsets(3), y_cb, &
                                                       & y_cc, dy, z_cb, z_cc, dz, glb_bounds)
            end if
        end block
        $:GPU_UPDATE(device='[glb_bounds]')
        dx_min = minval(dx)
        if (n > 0) dy_min = minval(dy)
        if (p > 0) dz_min = minval(dz)

        if (model_eqns == model_eqns_6eq) call s_initialize_internal_energy_equations(q_cons_ts(1)%vf)
        if (ib) then
            block
                type(ib_patch_parameters), allocatable :: particle_cloud_ibs(:)
                integer                                :: num_particle_cloud_ibs

                call s_instantiate_STL_models()
                call s_initialize_ib_airfoils()
                call s_get_neighbor_bounds()

                if (cfl_dt .and. n_start > 0) then
                    call s_read_ib_restart_data(n_start)
                    allocate (particle_cloud_ibs(0))
                    num_particle_cloud_ibs = 0
                else if (t_step_start > 0) then
                    call s_read_ib_restart_data(t_step_start)
                    allocate (particle_cloud_ibs(0))
                    num_particle_cloud_ibs = 0
                else
                    call s_generate_particle_clouds(particle_cloud_ibs, num_particle_cloud_ibs)
                end if
                call s_reduce_ib_patch_array(particle_cloud_ibs, num_particle_cloud_ibs)
                deallocate (particle_cloud_ibs)
            end block
            call s_ibm_setup()
            if (t_step_start == 0 .or. (cfl_dt .and. n_start == 0)) then
                call s_write_ib_data_file(0)
                call s_write_ib_state_file(0)
            end if
        end if
        if (bodyForces .or. synthetic_turbulence) call s_initialize_body_forces_module()
        if (acoustic_source) call s_precalculate_acoustic_spatial_sources()

        ! Initialize the Temperature cache.
        if (chemistry) call s_compute_q_T_sf(q_T_sf, q_cons_ts(1)%vf, idwint)

        ! Computation of parameters, allocation of memory, association of pointers, and/or execution of any other tasks that are
        ! needed to properly configure the modules. The preparations below DO DEPEND on the grid being complete.
        if (igr) then
            call s_initialize_igr_module()
        end if
        if (.not. igr) then
            if (recon_type == recon_type_weno) then
                call s_initialize_weno_module()
            else if (recon_type == recon_type_muscl) then
                call s_initialize_muscl_module()
            end if
            call s_initialize_cbc_module()
            call s_initialize_riemann_solvers_module()
        end if
        if (int_comp > 0) call s_initialize_thinc_module()
        call s_initialize_derived_variables()
        if (bubbles_lagrange) call s_initialize_bubbles_EL_module(q_cons_ts(1)%vf, bc_type)

        if (hypoelasticity) call s_initialize_hypoelastic_module()

    end subroutine s_initialize_modules

    !> Set up the MPI execution environment, bind GPUs, and decompose the computational domain
    impure subroutine s_initialize_mpi_domain

        integer :: ierr

#ifdef MFC_GPU
        real(wp) :: starttime, endtime
        integer  :: num_devices, local_size, num_nodes, ppn, my_device_num
        integer  :: dev, devNum, local_rank
#ifdef MFC_MPI
        integer :: local_comm
#endif
#if defined(MFC_OpenACC)
        integer(acc_device_kind) :: devtype
#endif
#endif

        call s_mpi_initialize()

#ifdef MFC_GPU
#ifndef MFC_MPI
        local_size = 1
        local_rank = 0
#else
        call MPI_Comm_split_type(MPI_COMM_WORLD, MPI_COMM_TYPE_SHARED, 0, MPI_INFO_NULL, local_comm, ierr)
        call MPI_Comm_size(local_comm, local_size, ierr)
        call MPI_Comm_rank(local_comm, local_rank, ierr)
#endif
#if defined(MFC_OpenACC)
        devtype = acc_get_device_type()
        devNum = acc_get_num_devices(devtype)
        dev = mod(local_rank, devNum)

        call acc_set_device_num(dev, devtype)
#elif defined(MFC_OpenMP)
        devNum = omp_get_num_devices()
        dev = mod(local_rank, devNum)
        call omp_set_default_device(dev)
#endif
#endif

        if (proc_rank == 0) then
            call s_assign_default_values_to_user_inputs()
            call s_read_input_file()
            call s_check_input_file()

            print '(" Simulating a ", A, " ", I0, "x", I0, "x", I0, " case on ", I0, " rank(s) ", A, ".")', &
            #:if not MFC_CASE_OPTIMIZATION
                "regular", &
            #:else
                "case-optimized", &
            #:endif
            m, n, p, num_procs, &
#if defined(MFC_OpenACC)
            "with OpenACC offloading"
#elif defined(MFC_OpenMP)
            "with OpenMP offloading"
#else
            "on CPUs"
#endif
        end if

        call s_mpi_bcast_user_inputs()

        ! Save original BCs before decomposition overwrites them with MPI neighbor ranks
        ib_bc_x = bc_x
        ib_bc_y = bc_y
        ib_bc_z = bc_z

        call s_initialize_parallel_io()

        call s_mpi_decompose_computational_domain(write_silo_ghost_offsets=.false., adjust_local_domains=.false.)

        bc = bc_xyz_info(bc_x, bc_y, bc_z)

    end subroutine s_initialize_mpi_domain

    !> Transfer initial conservative variable and model parameter data to the GPU device
    subroutine s_initialize_gpu_vars

        integer :: i

        if (.not. down_sample) then
            do i = 1, sys_size
                $:GPU_UPDATE(device='[q_cons_ts(1)%vf(i)%sf]')
            end do
        end if

        if (qbmm .and. .not. polytropic) then
            $:GPU_UPDATE(device='[pb_ts(1)%sf, mv_ts(1)%sf]')
        end if
        if (chemistry) then
            $:GPU_UPDATE(device='[q_T_sf%sf]')
        end if

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

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

        $:GPU_UPDATE(device='[R0ref, p0ref, rho0ref, ss, pv, vd, mu_l, mu_v, mu_g, gam_v, gam_g, M_v, M_g, R_v, R_g, Tw, cp_v, &
                     & cp_g, k_vl, k_gl, gam, gam_m, Eu, Ca, Web, Re_inv, Pe_c, phi_vg, phi_gv, omegaN, bubbles_euler, &
                     & polytropic, polydisperse, qbmm, ptil, bubble_model, thermal, poly_sigma, adv_n, adap_dt, adap_dt_tol, &
                     & adap_dt_max_iters, eqn_idx%n, pi_fac, low_Mach]')

        if (bubbles_euler) then
            $:GPU_UPDATE(device='[weight, R0]')
            if (.not. polytropic) then
                $:GPU_UPDATE(device='[pb0, Pe_T, k_g, k_v, mass_g0, mass_v0, Re_trans_T, Re_trans_c, Im_trans_T, Im_trans_c]')
            else if (qbmm) then
                $:GPU_UPDATE(device='[pb0]')
            end if
        end if

        $:GPU_UPDATE(device='[adv_n, adap_dt, adap_dt_tol, adap_dt_max_iters, pi_fac, low_Mach]')

        $:GPU_UPDATE(device='[acoustic_source, num_source]')
        $:GPU_UPDATE(device='[sigma, surface_tension]')

        $:GPU_UPDATE(device='[dx, dy, dz, x_cb, x_cc, y_cb, y_cc, z_cb, z_cc]')
        $:GPU_UPDATE(device='[bc_x%beg, bc_x%end, bc_y%beg, bc_y%end, bc_z%beg, bc_z%end]')
        $:GPU_UPDATE(device='[bc_x%vb1, bc_x%vb2, bc_x%vb3, bc_x%ve1, bc_x%ve2, bc_x%ve3]')
        $:GPU_UPDATE(device='[bc_y%vb1, bc_y%vb2, bc_y%vb3, bc_y%ve1, bc_y%ve2, bc_y%ve3]')
        $:GPU_UPDATE(device='[bc_z%vb1, bc_z%vb2, bc_z%vb3, bc_z%ve1, bc_z%ve2, bc_z%ve3]')

        $:GPU_UPDATE(device='[bc_x%grcbc_in, bc_x%grcbc_out, bc_x%grcbc_vel_out]')
        $:GPU_UPDATE(device='[bc_y%grcbc_in, bc_y%grcbc_out, bc_y%grcbc_vel_out]')
        $:GPU_UPDATE(device='[bc_z%grcbc_in, bc_z%grcbc_out, bc_z%grcbc_vel_out]')

        $:GPU_UPDATE(device='[bc_x%isothermal_in, bc_x%isothermal_out]')
        $:GPU_UPDATE(device='[bc_y%isothermal_in, bc_y%isothermal_out]')
        $:GPU_UPDATE(device='[bc_z%isothermal_in, bc_z%isothermal_out]')
        $:GPU_UPDATE(device='[bc_x%Twall_in, bc_x%Twall_out, bc_y%Twall_in, bc_y%Twall_out, bc_z%Twall_in, bc_z%Twall_out]')

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

        $:GPU_UPDATE(device='[relax, relax_model]')
        if (relax) then
            $:GPU_UPDATE(device='[palpha_eps, ptgalpha_eps]')
        end if

        if (ib) then
            $:GPU_UPDATE(device='[ib_markers%sf]')
        end if
        #:if not MFC_CASE_OPTIMIZATION
            $:GPU_UPDATE(device='[igr, nb, igr_order]')
        #:endif

    end subroutine s_initialize_gpu_vars

    !> Finalize and deallocate all simulation sub-modules in reverse initialization order
    impure subroutine s_finalize_modules

        if (model_eqns == model_eqns_6eq) call s_report_pressure_relaxation()

        call s_finalize_time_steppers_module()
        if (hypoelasticity) call s_finalize_hypoelastic_module()
        call s_finalize_derived_variables_module()
        call s_finalize_data_output_module()
        call s_finalize_rhs_module()
        if (igr) then
            call s_finalize_igr_module()
        else
            call s_finalize_cbc_module()
            call s_finalize_riemann_solvers_module()
            if (recon_type == recon_type_weno) then
                call s_finalize_weno_module()
            else if (recon_type == recon_type_muscl) then
                call s_finalize_muscl_module()
            end if
        end if
        if (int_comp > 0) call s_finalize_thinc_module()
        call s_finalize_variables_conversion_module()
        if (grid_geometry == 3) call s_finalize_fftw_module
        call s_finalize_mpi_common_module()
        call s_finalize_global_parameters_module()
        call s_finalize_boundary_common_module()
        if (relax) call s_finalize_relaxation_solver_module()
        if (bubbles_lagrange) call s_finalize_lagrangian_solver()
        if (viscous .and. (.not. igr)) then
            call s_finalize_viscous_module()
        end if
        call s_finalize_mpi_proxy_module()

        if (surface_tension) call s_finalize_surface_tension_module()
        if (bodyForces .or. synthetic_turbulence) call s_finalize_body_forces_module()
        if (ib) call s_finalize_ibm_module()

        call s_mpi_finalize()

    end subroutine s_finalize_modules

    !> @brief Reads IB kinematic state from restart_data/ib_state.dat on restart. Rank 0 reads the last num_ibs records and
    !! broadcasts to all ranks. Overwrites patch_ib vel, angular_vel, angles, and centroid.
    impure subroutine s_read_ib_restart_data(t_step)

        integer, intent(in)                  :: t_step
        character(len=path_len + 2*name_len) :: file_loc
        integer                              :: i, ios, file_unit, ierr
        integer                              :: r, nlocal, gbl_id
        integer, parameter                   :: NFIELDS_PER_IB = 20
        real(wp)                             :: ib_buf(NFIELDS_PER_IB)
        logical                              :: file_exist
        character(len=10)                    :: t_step_string

        if (file_per_process) then
            call s_int_to_str(t_step, t_step_string)

            do r = 0, num_procs - 1
                write (file_loc, '(A,I0,A,i7.7,A)') 'ib_state_', t_step, '_', r, '.dat'
                file_loc = trim(case_dir) // '/restart_data/lustre_' // trim(t_step_string) // '/' // trim(file_loc)

                inquire (FILE=trim(file_loc), EXIST=file_exist)
                if (.not. file_exist) call s_mpi_abort('Cannot open IB state file for restart: ' // trim(file_loc))

                open (newunit=file_unit, file=trim(file_loc), form='unformatted', access='stream', status='old', iostat=ios)
                if (ios /= 0) call s_mpi_abort('Error opening IB state restart file: ' // trim(file_loc))

                read (file_unit, iostat=ios) nlocal
                if (ios /= 0) call s_mpi_abort('Error reading IB state file header: ' // trim(file_loc))

                do i = 1, nlocal
                    read (file_unit, iostat=ios) gbl_id
                    if (ios /= 0) call s_mpi_abort('Error reading IB patch ID: ' // trim(file_loc))
                    read (file_unit, iostat=ios) ib_buf
                    if (ios /= 0) call s_mpi_abort('Error reading IB state data: ' // trim(file_loc))

                    patch_ib(gbl_id)%vel = ib_buf(8:10)
                    patch_ib(gbl_id)%angular_vel = ib_buf(11:13)
                    patch_ib(gbl_id)%angles = ib_buf(14:16)
                    patch_ib(gbl_id)%x_centroid = ib_buf(17)
                    patch_ib(gbl_id)%y_centroid = ib_buf(18)
                    patch_ib(gbl_id)%z_centroid = ib_buf(19)
                end do

                close (file_unit)
            end do
        else
            write (file_loc, '(A,I0,A)') '/restart_data/ib_state_', t_step, '.dat'
            file_loc = trim(case_dir) // trim(file_loc)

            if (proc_rank == 0) then
                inquire (FILE=trim(file_loc), EXIST=file_exist)
                if (.not. file_exist) then
                    call s_mpi_abort('Cannot open IB state file for restart: ' // trim(file_loc))
                end if

                open (newunit=file_unit, file=trim(file_loc), form='unformatted', access='stream', status='old', iostat=ios)
                if (ios /= 0) call s_mpi_abort('Error opening IB state restart file: ' // trim(file_loc))

                do i = 1, num_ibs
                    read (file_unit, iostat=ios) ib_buf
                    if (ios /= 0) call s_mpi_abort('Error reading IB state restart file')

                    patch_ib(i)%vel = ib_buf(8:10)
                    patch_ib(i)%angular_vel = ib_buf(11:13)
                    patch_ib(i)%angles = ib_buf(14:16)
                    patch_ib(i)%x_centroid = ib_buf(17)
                    patch_ib(i)%y_centroid = ib_buf(18)
                    patch_ib(i)%z_centroid = ib_buf(19)
                end do

                close (file_unit)
            end if

#ifdef MFC_MPI
            do i = 1, num_ibs
                call MPI_BCAST(patch_ib(i)%vel, 3, mpi_p, 0, MPI_COMM_WORLD, ierr)
                call MPI_BCAST(patch_ib(i)%angular_vel, 3, mpi_p, 0, MPI_COMM_WORLD, ierr)
                call MPI_BCAST(patch_ib(i)%angles, 3, mpi_p, 0, MPI_COMM_WORLD, ierr)
                call MPI_BCAST(patch_ib(i)%x_centroid, 1, mpi_p, 0, MPI_COMM_WORLD, ierr)
                call MPI_BCAST(patch_ib(i)%y_centroid, 1, mpi_p, 0, MPI_COMM_WORLD, ierr)
                call MPI_BCAST(patch_ib(i)%z_centroid, 1, mpi_p, 0, MPI_COMM_WORLD, ierr)
            end do
#endif
        end if

    end subroutine s_read_ib_restart_data

    !> @brief Merges patch_ib (namelist patches, fixed at num_ib_patches_max_namelist) with particle_cloud_ibs (already filtered by
    !! s_generate_particle_clouds to this rank's IB neighborhood, each entry already tagged with its final, absolute gbl_patch_id)
    !! and reduces to only the patches in or near the local computational domain. patch_ib is never reallocated; the local subset is
    !! written in-place from the front. particle_cloud_ibs is owned by the caller and freed there after this returns.
    !! num_particle_cloud_ibs is the number of entries s_generate_particle_clouds actually wrote into particle_cloud_ibs - it may be
    !! allocated to a larger worst-case capacity, so size() of it must never be used as the valid-entry count.
    subroutine s_reduce_ib_patch_array(particle_cloud_ibs, num_particle_cloud_ibs)

        type(ib_patch_parameters), intent(in), dimension(:) :: particle_cloud_ibs
        integer, intent(in)                                 :: num_particle_cloud_ibs
        real(wp), dimension(3)                              :: centroid
        integer                                             :: i
        integer                                             :: num_namelist_ibs, num_bed_ibs

        num_namelist_ibs = num_ibs
        num_bed_ibs = 0
        do i = 1, num_particle_clouds
            num_bed_ibs = num_bed_ibs + particle_cloud(i)%num_particles
        end do

        ! Check for moving IBs across both namelist and particle cloud patches.
        moving_immersed_boundary_flag = .false.
        do i = 1, num_namelist_ibs
            if (patch_ib(i)%moving_ibm /= 0) then
                moving_immersed_boundary_flag = .true.
                exit
            end if
        end do
        if (.not. moving_immersed_boundary_flag) then
            do i = 1, num_particle_clouds
                if (particle_cloud(i)%moving_ibm /= 0) then
                    moving_immersed_boundary_flag = .true.
                    exit
                end if
            end do
        end if

        call s_compute_ib_neighbor_ranks()

#ifdef MFC_MPI
        if (num_procs == 1) then
            ! single-rank: all patches are local; append particle bed entries directly into patch_ib.
            do i = 1, num_particle_cloud_ibs
                patch_ib(num_namelist_ibs + i) = particle_cloud_ibs(i)
            end do
            num_gbl_ibs = num_namelist_ibs + num_particle_cloud_ibs
            @:PROHIBIT(num_gbl_ibs > num_ib_patches_max_namelist, &
                       & "Total IB count exceeds patch_ib capacity. Increase num_ib_patches_max_namelist.")
            num_ibs = num_gbl_ibs
            num_local_ibs = num_gbl_ibs
            do i = 1, num_gbl_ibs
                local_ib_patch_ids(i) = i
            end do
        else
            ! multi-rank: compact namelist patches in-place (write_idx <= read_idx, no aliasing), then append local particle beds.
            num_ibs = 0
            num_local_ibs = 0
            num_gbl_ibs = num_namelist_ibs + num_bed_ibs
            do i = 1, num_namelist_ibs
                centroid = [patch_ib(i)%x_centroid, patch_ib(i)%y_centroid, 0._wp]
                if (num_dims == 3) centroid(3) = patch_ib(i)%z_centroid
                if (f_neighborhood_ranks_own_location(centroid)) then
                    num_ibs = num_ibs + 1
                    patch_ib(num_ibs) = patch_ib(i)
                    patch_ib(num_ibs)%gbl_patch_id = i
                    if (f_local_rank_owns_location(centroid)) then
                        num_local_ibs = num_local_ibs + 1
                        @:PROHIBIT(num_local_ibs > num_local_ibs_max, &
                                   & "Too many IBs on a single processor rank. Modify case file or increase limit of num_local_ibs_max to resolve.")
                        local_ib_patch_ids(num_local_ibs) = num_ibs
                    end if
                end if
            end do
            ! particle_cloud_ibs entries already passed the neighborhood check at generation time so no need to recheck it here.
            do i = 1, num_particle_cloud_ibs
                centroid = [particle_cloud_ibs(i)%x_centroid, particle_cloud_ibs(i)%y_centroid, 0._wp]
                if (num_dims == 3) centroid(3) = particle_cloud_ibs(i)%z_centroid
                num_ibs = num_ibs + 1
                @:PROHIBIT(num_ibs > num_ib_patches_max_namelist, &
                           & "Local IB count exceeds patch_ib capacity. Increase num_ib_patches_max_namelist.")
                patch_ib(num_ibs) = particle_cloud_ibs(i)
                if (f_local_rank_owns_location(centroid)) then
                    num_local_ibs = num_local_ibs + 1
                    @:PROHIBIT(num_local_ibs > num_local_ibs_max, &
                               & "Too many IBs on a single processor rank. Modify case file or increase limit of num_local_ibs_max to resolve.")
                    local_ib_patch_ids(num_local_ibs) = num_ibs
                end if
            end do
        end if
#else
        ! no-MPI: all patches are local; append particle bed entries directly into patch_ib.
        do i = 1, num_particle_cloud_ibs
            patch_ib(num_namelist_ibs + i) = particle_cloud_ibs(i)
        end do
        num_gbl_ibs = num_namelist_ibs + num_particle_cloud_ibs
        @:PROHIBIT(num_gbl_ibs > num_ib_patches_max_namelist, &
                   & "Total IB count exceeds patch_ib capacity. Increase num_ib_patches_max_namelist.")
        num_ibs = num_gbl_ibs
        num_local_ibs = num_gbl_ibs
        do i = 1, num_gbl_ibs
            local_ib_patch_ids(i) = i
        end do
#endif

        @:ALLOCATE(ib_gbl_idx_lookup(1:num_gbl_ibs))

    end subroutine s_reduce_ib_patch_array

    !> Build ib_neighbor_ranks(-1:1,-1:1,-1:1): MPI ranks of all neighbor domains. Uses two rounds of MPI_SENDRECV cascades - face
    !! neighbors are known from bc_*, edge neighbors are obtained in round 1, and (3D) corner neighbors in round 2.
    subroutine s_compute_ib_neighbor_ranks()

        integer                :: ax, k, nbr_idx, nreqs, sx, sy, sz, dx, dy, dz
        integer, allocatable   :: send_table(:,:,:), recv_tables(:,:,:,:)
        integer, dimension(52) :: requests

#ifdef MFC_MPI
        integer               :: ierr
        integer, dimension(4) :: buf4, rbuf4
        integer, dimension(2) :: buf2, rbuf2

        ax = ib_neighborhood_radius

        if (allocated(ib_neighbor_ranks)) deallocate (ib_neighbor_ranks)
        allocate (ib_neighbor_ranks(-ax:ax,-ax:ax,-ax:ax))
        ib_neighbor_ranks = MPI_PROC_NULL
        ib_neighbor_ranks(0, 0, 0) = proc_rank

        ! Fill radius-1 entries: face neighbors are known from domain decomposition
        ib_neighbor_ranks(-1, 0, 0) = bc_x%beg
        ib_neighbor_ranks(+1, 0, 0) = bc_x%end
        if (num_dims >= 2) then
            ib_neighbor_ranks(0, -1, 0) = bc_y%beg
            ib_neighbor_ranks(0, +1, 0) = bc_y%end
        end if
        if (num_dims == 3) then
            ib_neighbor_ranks(0, 0, -1) = bc_z%beg
            ib_neighbor_ranks(0, 0, +1) = bc_z%end
        end if

        if (num_dims >= 2) then
            ! Round 1a: exchange y/z face ranks with +/-x face neighbors -> xy and xz edge ranks
            buf4 = [bc_y%beg, bc_y%end, bc_z%beg, bc_z%end]

            ! Send to -x, receive from +x -> edges (+1,+/-1,0) and (+1,0,+/-1)
            call MPI_SENDRECV(buf4, 4, MPI_INTEGER, merge(bc_x%beg, MPI_PROC_NULL, bc_x%beg >= 0), 310, rbuf4, 4, MPI_INTEGER, &
                              & merge(bc_x%end, MPI_PROC_NULL, bc_x%end >= 0), 310, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr)
            if (bc_x%end >= 0) then
                ib_neighbor_ranks(+1, -1, 0) = rbuf4(1)
                ib_neighbor_ranks(+1, +1, 0) = rbuf4(2)
                ib_neighbor_ranks(+1, 0, -1) = rbuf4(3)
                ib_neighbor_ranks(+1, 0, +1) = rbuf4(4)
            end if

            call MPI_SENDRECV(buf4, 4, MPI_INTEGER, merge(bc_x%end, MPI_PROC_NULL, bc_x%end >= 0), 311, rbuf4, 4, MPI_INTEGER, &
                              & merge(bc_x%beg, MPI_PROC_NULL, bc_x%beg >= 0), 311, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr)
            if (bc_x%beg >= 0) then
                ib_neighbor_ranks(-1, -1, 0) = rbuf4(1)
                ib_neighbor_ranks(-1, +1, 0) = rbuf4(2)
                ib_neighbor_ranks(-1, 0, -1) = rbuf4(3)
                ib_neighbor_ranks(-1, 0, +1) = rbuf4(4)
            end if
        end if

        if (num_dims == 3) then
            ! Round 1b: exchange z face ranks with +/-y face neighbors -> yz edge ranks
            buf2 = [bc_z%beg, bc_z%end]

            call MPI_SENDRECV(buf2, 2, MPI_INTEGER, merge(bc_y%beg, MPI_PROC_NULL, bc_y%beg >= 0), 312, rbuf2, 2, MPI_INTEGER, &
                              & merge(bc_y%end, MPI_PROC_NULL, bc_y%end >= 0), 312, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr)
            if (bc_y%end >= 0) then
                ib_neighbor_ranks(0, +1, -1) = rbuf2(1)
                ib_neighbor_ranks(0, +1, +1) = rbuf2(2)
            end if

            call MPI_SENDRECV(buf2, 2, MPI_INTEGER, merge(bc_y%end, MPI_PROC_NULL, bc_y%end >= 0), 313, rbuf2, 2, MPI_INTEGER, &
                              & merge(bc_y%beg, MPI_PROC_NULL, bc_y%beg >= 0), 313, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr)
            if (bc_y%beg >= 0) then
                ib_neighbor_ranks(0, -1, -1) = rbuf2(1)
                ib_neighbor_ranks(0, -1, +1) = rbuf2(2)
            end if

            ! Round 2: exchange z face ranks with xy-diagonal edge neighbors -> corner ranks. Each of the 4 xy diagonals gives 2
            ! corners (the +/-z variants). Pattern: send buf2 to mirror diagonal, receive from this diagonal -> that edge's z face
            ! ranks.
            #:for DX, DY, MDX, MDY, TAG in [(1,1,-1,-1,320), (1,-1,-1,1,321), (-1,1,1,-1,322), (-1,-1,1,1,323)]
                call MPI_SENDRECV(buf2, 2, MPI_INTEGER, merge(ib_neighbor_ranks(${MDX}$, ${MDY}$, 0), MPI_PROC_NULL, &
                                  & ib_neighbor_ranks(${MDX}$, ${MDY}$, 0) >= 0), ${TAG}$, rbuf2, 2, MPI_INTEGER, &
                                  & merge(ib_neighbor_ranks(${DX}$, ${DY}$, 0), MPI_PROC_NULL, ib_neighbor_ranks(${DX}$, ${DY}$, &
                                  & 0) >= 0), ${TAG}$, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr)
                if (ib_neighbor_ranks(${DX}$, ${DY}$, 0) >= 0) then
                    ib_neighbor_ranks(${DX}$, ${DY}$, -1) = rbuf2(1)
                    ib_neighbor_ranks(${DX}$, ${DY}$, +1) = rbuf2(2)
                end if
            #:endfor
        end if

        ! For radius > 1: extend the table by iterative 26-neighbor full-table exchanges. In each round, every rank broadcasts its
        ! current table to all 26 immediate neighbors. Their entry at offset (dx,dy,dz) from them = our entry at
        ! (dx+sx,dy+sy,dz+sz). One extension round fills the entire next shell, so ax-1 rounds suffice.
        if (ax > 1) then
            allocate (send_table(-ax:ax,-ax:ax,-ax:ax))
            allocate (recv_tables(-ax:ax,-ax:ax,-ax:ax,1:26))

            do k = 2, ax
                send_table = ib_neighbor_ranks

                nreqs = 0
                nbr_idx = 0
                do sz = -1, 1
                    do sy = -1, 1
                        do sx = -1, 1
                            if (sx == 0 .and. sy == 0 .and. sz == 0) cycle
                            nbr_idx = nbr_idx + 1
                            if (ib_neighbor_ranks(sx, sy, sz) < 0) cycle
                            nreqs = nreqs + 1
                            call MPI_IRECV(recv_tables(:,:,:,nbr_idx), (2*ax + 1)**3, MPI_INTEGER, ib_neighbor_ranks(sx, sy, sz), &
                                           & 400, MPI_COMM_WORLD, requests(nreqs), ierr)
                        end do
                    end do
                end do

                do sz = -1, 1
                    do sy = -1, 1
                        do sx = -1, 1
                            if (sx == 0 .and. sy == 0 .and. sz == 0) cycle
                            if (ib_neighbor_ranks(sx, sy, sz) < 0) cycle
                            nreqs = nreqs + 1
                            call MPI_ISEND(send_table, (2*ax + 1)**3, MPI_INTEGER, ib_neighbor_ranks(sx, sy, sz), 400, &
                                           & MPI_COMM_WORLD, requests(nreqs), ierr)
                        end do
                    end do
                end do

                call MPI_WAITALL(nreqs, requests, MPI_STATUSES_IGNORE, ierr)

                nbr_idx = 0
                do sz = -1, 1
                    do sy = -1, 1
                        do sx = -1, 1
                            if (sx == 0 .and. sy == 0 .and. sz == 0) cycle
                            nbr_idx = nbr_idx + 1
                            if (ib_neighbor_ranks(sx, sy, sz) < 0) cycle
                            do dz = -ax, ax
                                do dy = -ax, ax
                                    do dx = -ax, ax
                                        if (recv_tables(dx, dy, dz, nbr_idx) == MPI_PROC_NULL) cycle
                                        if (dx + sx < -ax .or. dx + sx > ax) cycle
                                        if (dy + sy < -ax .or. dy + sy > ax) cycle
                                        if (dz + sz < -ax .or. dz + sz > ax) cycle
                                        if (ib_neighbor_ranks(dx + sx, dy + sy, dz + sz) /= MPI_PROC_NULL) cycle
                                        ib_neighbor_ranks(dx + sx, dy + sy, dz + sz) = recv_tables(dx, dy, dz, nbr_idx)
                                    end do
                                end do
                            end do
                        end do
                    end do
                end do
            end do

            deallocate (send_table, recv_tables)
        end if
#endif

    end subroutine s_compute_ib_neighbor_ranks

    subroutine s_get_neighbor_bounds()

        real(wp) :: beg_val, end_val, recv_val, bound, max_ib_bound, local_rank_width, min_rank_width
        integer  :: k, send_neighbor, recv_neighbor, ierr, temporary_radius

        ! Default: unbounded in all directions (covers single-rank and no-MPI cases)

        neighbor_domain_x%beg = -huge(0._wp)
        neighbor_domain_x%end = huge(0._wp)
        neighbor_domain_y%beg = -huge(0._wp)
        neighbor_domain_y%end = huge(0._wp)
        neighbor_domain_z%beg = -huge(0._wp)
        neighbor_domain_z%end = huge(0._wp)

#ifdef MFC_MPI
        ! perform setup if we are doing automatic radius checking
        if (ib_neighborhood_radius < 1) then
            ib_neighborhood_radius = 0  ! ensure we are starting with 0 neighborhood radius

            ! determine the maximum length of space that needs to be contained by the neighborhood
            max_ib_bound = -1._wp
            do k = 1, num_ibs
                call s_get_ib_bound(patch_ib(k), bound)
                max_ib_bound = max(max_ib_bound, bound)
            end do
            do k = 1, num_particle_clouds
                max_ib_bound = max(max_ib_bound, particle_cloud(k)%radius)
            end do

            ! determine the upper bound on the size
            local_rank_width = -1._wp
            #:for X, ID, DIM in [('x', 1, 'm'), ('y', 2, 'n'), ('z', 3, 'p')]
                if (num_dims >= ${ID}$) local_rank_width = max(local_rank_width, abs(${X}$_cb(${DIM}$) - ${X}$_cb(-1)))
            #:endfor
            call s_mpi_allreduce_min(local_rank_width, min_rank_width)

            ! approximate the size of the neighborhood with a local 1.1x fudge factor for safety, lower bound of 1
            ib_neighborhood_radius = max(1, ceiling(1.1_wp*max_ib_bound/(min_rank_width)))
            if (proc_rank == 0) print *, "Automatic choice of ib_neighborhood_radius selected: ", ib_neighborhood_radius
        end if

        ! For each direction, propagate the left/right boundary edges outward ib_neighborhood_radius hops. After k rounds: beg_val =
        ! left edge of the rank k hops to the left; end_val = right edge of the rank k hops to the right.
        #:for X, ID, TAG, DIM in [('x', 1, 100, 'm'), ('y', 2, 102, 'n'), ('z', 3, 104, 'p')]
            if (num_dims >= ${ID}$) then
                beg_val = ${X}$_cb(-1)
                end_val = ${X}$_cb(${DIM}$)
                do k = 1, ib_neighborhood_radius
                    send_neighbor = merge(bc_${X}$%end, MPI_PROC_NULL, bc_${X}$%end >= 0)
                    recv_neighbor = merge(bc_${X}$%beg, MPI_PROC_NULL, bc_${X}$%beg >= 0)
                    recv_val = -huge(0._wp)
                    call MPI_SENDRECV(beg_val, 1, mpi_p, send_neighbor, ${TAG}$, recv_val, 1, mpi_p, recv_neighbor, ${TAG}$, &
                                      & MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr)
                    beg_val = recv_val

                    send_neighbor = merge(bc_${X}$%beg, MPI_PROC_NULL, bc_${X}$%beg >= 0)
                    recv_neighbor = merge(bc_${X}$%end, MPI_PROC_NULL, bc_${X}$%end >= 0)
                    recv_val = huge(0._wp)
                    call MPI_SENDRECV(end_val, 1, mpi_p, send_neighbor, ${TAG + 1}$, recv_val, 1, mpi_p, recv_neighbor, &
                                      & ${TAG + 1}$, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr)
                    end_val = recv_val

                    ! protect from looping back around on yourself multiple times
                    if (f_approx_equal(beg_val, ${X}$_cb(${DIM}$)) .or. f_approx_equal(end_val, ${X}$_cb(-1))) then
                        beg_val = -huge(0._wp)
                        end_val = huge(0._wp)
                        exit
                    end if
                end do
                neighbor_domain_${X}$%beg = beg_val
                neighbor_domain_${X}$%end = end_val
            end if
        #:endfor
#endif

    end subroutine s_get_neighbor_bounds

end module m_start_up