m_bubbles_EL.fpp

!>
!! @file
!! @brief Contains module @ref m_bubbles_el "m_bubbles_EL"

#:include 'macros.fpp'

!> @brief Tracks Lagrangian bubbles and couples their dynamics to the Eulerian flow via volume averaging
module m_bubbles_EL

    use m_global_parameters
    use m_mpi_proxy
    use m_bubbles_EL_kernels
    use m_bubbles
    use m_variables_conversion
    use m_compile_specific
    use m_boundary_common
    use m_helper_basic
    use m_sim_helpers
    use m_helper
    use m_mpi_common
    use m_ibm
    use m_finite_differences
    use m_constants, only: time_stepper_rk1, time_stepper_rk2, time_stepper_rk3, precision_single

    implicit none

    ! (nBub)
    integer, allocatable, dimension(:,:) :: lag_id      !< Global and local IDs
    real(wp), allocatable, dimension(:)  :: bub_R0      !< Initial bubble radius
    real(wp), allocatable, dimension(:)  :: Rmax_stats  !< Maximum radius
    real(wp), allocatable, dimension(:)  :: Rmin_stats  !< Minimum radius
    $:GPU_DECLARE(create='[lag_id, bub_R0, Rmax_stats, Rmin_stats]')

    real(wp), allocatable, dimension(:) :: gas_mg      !< Bubble's gas mass
    real(wp), allocatable, dimension(:) :: gas_betaT   !< heatflux model (Preston et al., 2007)
    real(wp), allocatable, dimension(:) :: gas_betaC   !< massflux model (Preston et al., 2007)
    real(wp), allocatable, dimension(:) :: bub_dphidt  !< subgrid velocity potential (Maeda & Colonius, 2018)
    $:GPU_DECLARE(create='[gas_mg, gas_betaT, gas_betaC, bub_dphidt]')

    ! (nBub, 1 -> actual val or 2 -> temp val)
    real(wp), allocatable, dimension(:,:) :: gas_p      !< Pressure in the bubble
    real(wp), allocatable, dimension(:,:) :: gas_mv     !< Vapor mass in the bubble
    real(wp), allocatable, dimension(:,:) :: intfc_rad  !< Bubble radius
    real(wp), allocatable, dimension(:,:) :: intfc_vel  !< Velocity of the bubble interface
    $:GPU_DECLARE(create='[gas_p, gas_mv, intfc_rad, intfc_vel]')
    ! (nBub, 1-> x or 2->y or 3 ->z, 1 -> actual or 2 -> temporal val)
    real(wp), allocatable, dimension(:,:,:) :: mtn_pos      !< Bubble's position
    real(wp), allocatable, dimension(:,:,:) :: mtn_posPrev  !< Bubble's previous position
    real(wp), allocatable, dimension(:,:,:) :: mtn_vel      !< Bubble's velocity
    real(wp), allocatable, dimension(:,:,:) :: mtn_s        !< Bubble's computational cell position in real format
    $:GPU_DECLARE(create='[mtn_pos, mtn_posPrev, mtn_vel, mtn_s]')
    ! (nBub, 1-> x or 2->y or 3 ->z, time-stage)
    real(wp), allocatable, dimension(:,:)   :: intfc_draddt  !< Time derivative of bubble's radius
    real(wp), allocatable, dimension(:,:)   :: intfc_dveldt  !< Time derivative of bubble's interface velocity
    real(wp), allocatable, dimension(:,:)   :: gas_dpdt      !< Time derivative of gas pressure
    real(wp), allocatable, dimension(:,:)   :: gas_dmvdt     !< Time derivative of the vapor mass in the bubble
    real(wp), allocatable, dimension(:,:,:) :: mtn_dposdt    !< Time derivative of the bubble's position
    real(wp), allocatable, dimension(:,:,:) :: mtn_dveldt    !< Time derivative of the bubble's velocity
    $:GPU_DECLARE(create='[intfc_draddt, intfc_dveldt, gas_dpdt, gas_dmvdt, mtn_dposdt, mtn_dveldt]')

    integer, private :: lag_num_ts  !< Number of time stages in the time-stepping scheme
    $:GPU_DECLARE(create='[lag_num_ts]')

    real(wp) :: Rmax_glb, Rmin_glb  !< Maximum and minimum bubbe size in the local domain
    !> Projection of the lagrangian particles in the Eulerian framework
    type(scalar_field), dimension(:), allocatable :: q_beta
    type(scalar_field), dimension(:), allocatable :: kahan_comp  !< Kahan compensation for q_beta accumulation
    integer                                       :: q_beta_idx  !< Size of the q_beta vector field

    $:GPU_DECLARE(create='[Rmax_glb, Rmin_glb, q_beta, kahan_comp, q_beta_idx]')

    integer, parameter                   :: LAG_EVOL_ID = 11  ! File id for lag_bubbles_evol_*.dat
    integer, parameter                   :: LAG_STATS_ID = 12  ! File id for stats_lag_bubbles_*.dat
    integer, parameter                   :: LAG_VOID_ID = 13  ! File id for voidfraction.dat
    integer, allocatable, dimension(:)   :: keep_bubble
    integer, allocatable, dimension(:,:) :: wrap_bubble_loc, wrap_bubble_dir
    $:GPU_DECLARE(create='[keep_bubble]')
    $:GPU_DECLARE(create='[wrap_bubble_loc, wrap_bubble_dir]')

contains

    !> Initializes the lagrangian subgrid bubble solver
    impure subroutine s_initialize_bubbles_EL_module(q_cons_vf, bc_type)

        type(scalar_field), dimension(sys_size), intent(inout)     :: q_cons_vf
        type(integer_field), dimension(1:num_dims,1:2), intent(in) :: bc_type
        integer                                                    :: nBubs_glb, i

        ! Setting number of time-stages for selected time-stepping scheme

        lag_num_ts = time_stepper

        ! Allocate space for the Eulerian fields needed to map the effect of the bubbles
        if (lag_params%solver_approach == 1) then
            ! One-way coupling
            q_beta_idx = 3
        else if (lag_params%solver_approach == 2) then
            ! Two-way coupling
            q_beta_idx = 4
            if (p == 0) then
                ! Subgrid noise model for 2D approximation
                q_beta_idx = 6
            end if
        else
            call s_mpi_abort('Please check the lag_params%solver_approach input')
        end if

        pcomm_coords(1)%beg = x_cb(-1)
        pcomm_coords(1)%end = x_cb(m)
        $:GPU_UPDATE(device='[pcomm_coords(1)]')
        if (n > 0) then
            pcomm_coords(2)%beg = y_cb(-1)
            pcomm_coords(2)%end = y_cb(n)
            $:GPU_UPDATE(device='[pcomm_coords(2)]')
            if (p > 0) then
                pcomm_coords(3)%beg = z_cb(-1)
                pcomm_coords(3)%end = z_cb(p)
                $:GPU_UPDATE(device='[pcomm_coords(3)]')
            end if
        end if

        $:GPU_UPDATE(device='[lag_num_ts, q_beta_idx]')

        @:ALLOCATE(q_beta(1:q_beta_idx))
        if (lag_params%kahan_summation) then
            @:ALLOCATE(kahan_comp(1:q_beta_idx))
        end if

        do i = 1, q_beta_idx
            @:ALLOCATE(q_beta(i)%sf(idwbuff(1)%beg:idwbuff(1)%end, idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end))
            @:ACC_SETUP_SFs(q_beta(i))
            if (lag_params%kahan_summation) then
                @:ALLOCATE(kahan_comp(i)%sf(idwbuff(1)%beg:idwbuff(1)%end, idwbuff(2)%beg:idwbuff(2)%end, &
                           & idwbuff(3)%beg:idwbuff(3)%end))
                @:ACC_SETUP_SFs(kahan_comp(i))
            end if
        end do

        ! Allocating space for lagrangian variables
        nBubs_glb = lag_params%nBubs_glb

        @:ALLOCATE(bub_R0(1:nBubs_glb))
        @:ALLOCATE(Rmax_stats(1:nBubs_glb))
        @:ALLOCATE(Rmin_stats(1:nBubs_glb))
        @:ALLOCATE(gas_mg(1:nBubs_glb))
        @:ALLOCATE(gas_betaT(1:nBubs_glb))
        @:ALLOCATE(gas_betaC(1:nBubs_glb))
        @:ALLOCATE(bub_dphidt(1:nBubs_glb))
        @:ALLOCATE(lag_id(1:nBubs_glb, 1:2))
        @:ALLOCATE(gas_p(1:nBubs_glb, 1:2))
        @:ALLOCATE(gas_mv(1:nBubs_glb, 1:2))
        @:ALLOCATE(intfc_rad(1:nBubs_glb, 1:2))
        @:ALLOCATE(intfc_vel(1:nBubs_glb, 1:2))
        @:ALLOCATE(mtn_pos(1:nBubs_glb, 1:3, 1:2))
        @:ALLOCATE(mtn_posPrev(1:nBubs_glb, 1:3, 1:2))
        @:ALLOCATE(mtn_vel(1:nBubs_glb, 1:3, 1:2))
        @:ALLOCATE(mtn_s(1:nBubs_glb, 1:3, 1:2))
        @:ALLOCATE(intfc_draddt(1:nBubs_glb, 1:lag_num_ts))
        @:ALLOCATE(intfc_dveldt(1:nBubs_glb, 1:lag_num_ts))
        @:ALLOCATE(gas_dpdt(1:nBubs_glb, 1:lag_num_ts))
        @:ALLOCATE(gas_dmvdt(1:nBubs_glb, 1:lag_num_ts))
        @:ALLOCATE(mtn_dposdt(1:nBubs_glb, 1:3, 1:lag_num_ts))
        @:ALLOCATE(mtn_dveldt(1:nBubs_glb, 1:3, 1:lag_num_ts))

        @:ALLOCATE(keep_bubble(1:nBubs_glb))
        @:ALLOCATE(wrap_bubble_loc(1:nBubs_glb, 1:num_dims), wrap_bubble_dir(1:nBubs_glb, 1:num_dims))

        if (adap_dt .and. f_is_default(adap_dt_tol)) adap_dt_tol = dflt_adap_dt_tol

        if (num_procs > 1) call s_initialize_particles_mpi(lag_num_ts)

        ! Starting bubbles
        if (lag_params%write_void_evol) call s_open_void_evol
        if (lag_params%write_bubbles) call s_open_lag_bubble_evol()
        if (lag_params%write_bubbles_stats) call s_open_lag_bubble_stats()

        if (lag_params%vel_model > 0) then
            moving_lag_bubbles = .true.
            lag_pressure_force = lag_params%pressure_force
            lag_gravity_force = lag_params%gravity_force
            lag_vel_model = lag_params%vel_model
            lag_drag_model = lag_params%drag_model
        end if
        $:GPU_UPDATE(device='[moving_lag_bubbles, lag_pressure_force, lag_gravity_force, lag_vel_model, lag_drag_model]')

        if (lag_params%vel_model > 0 .and. lag_params%pressure_force) then
            @:ALLOCATE(grad_p_x(0:m, 0:n, 0:p))
            ! s_compute_finite_difference_coefficients always extends fd_number beyond the interior on each side
            @:ALLOCATE(fd_coeff_x_pgrad(-fd_number:fd_number,-fd_number:m + fd_number))
            call s_compute_finite_difference_coefficients(m, x_cc, fd_coeff_x_pgrad, buff_size, fd_number, fd_order)
            $:GPU_UPDATE(device='[fd_coeff_x_pgrad]')
            if (n > 0) then
                @:ALLOCATE(grad_p_y(0:m, 0:n, 0:p))
                @:ALLOCATE(fd_coeff_y_pgrad(-fd_number:fd_number,-fd_number:n + fd_number))
                call s_compute_finite_difference_coefficients(n, y_cc, fd_coeff_y_pgrad, buff_size, fd_number, fd_order)
                $:GPU_UPDATE(device='[fd_coeff_y_pgrad]')
            end if
            if (p > 0) then
                @:ALLOCATE(grad_p_z(0:m, 0:n, 0:p))
                @:ALLOCATE(fd_coeff_z_pgrad(-fd_number:fd_number,-fd_number:p + fd_number))
                call s_compute_finite_difference_coefficients(p, z_cc, fd_coeff_z_pgrad, buff_size, fd_number, fd_order)
                $:GPU_UPDATE(device='[fd_coeff_z_pgrad]')
            end if
        end if

        @:ALLOCATE(cell_list_start(0:m, 0:n, 0:p))
        @:ALLOCATE(cell_list_count(0:m, 0:n, 0:p))
        @:ALLOCATE(cell_list_idx(1:lag_params%nBubs_glb))

        call s_read_input_bubbles(q_cons_vf, bc_type)

    end subroutine s_initialize_bubbles_EL_module

    !> Read initial bubble data from input files
    impure subroutine s_read_input_bubbles(q_cons_vf, bc_type)

        type(scalar_field), dimension(sys_size), intent(inout)     :: q_cons_vf
        type(integer_field), dimension(1:num_dims,1:2), intent(in) :: bc_type
        real(wp), dimension(8)                                     :: inputBubble
        real(wp)                                                   :: qtime
        integer                                                    :: id, bub_id, save_count
        integer                                                    :: i, ios
        logical                                                    :: file_exist, indomain
        integer, dimension(3)                                      :: cell
        character(LEN=path_len + 2*name_len)                       :: path_D_dir

        ! Initialize number of particles

        bub_id = 0
        id = 0

        ! Read the input lag_bubble file or restart point
        if (cfl_dt) then
            save_count = n_start
            qtime = n_start*t_save
        else
            save_count = t_step_start
            qtime = t_step_start*dt
        end if

        if (save_count == 0) then
            if (proc_rank == 0) print *, 'Reading lagrange bubbles input file.'
            call my_inquire(trim(lag_params%input_path), file_exist)
            if (file_exist) then
                open (94, file=trim(lag_params%input_path), form='formatted', iostat=ios)
                do while (ios == 0)
                    read (94, *, iostat=ios) (inputBubble(i), i=1, 8)
                    if (ios /= 0) cycle
                    indomain = particle_in_domain_physical(inputBubble(1:3))
                    id = id + 1
                    if (id > lag_params%nBubs_glb .and. proc_rank == 0) then
                        call s_mpi_abort("Current number of bubbles is larger than nBubs_glb")
                    end if
                    if (indomain) then
                        bub_id = bub_id + 1
                        call s_add_bubbles(inputBubble, q_cons_vf, bub_id)
                        lag_id(bub_id, 1) = id  ! global ID
                        lag_id(bub_id, 2) = bub_id  ! local ID
                        n_el_bubs_loc = bub_id  ! local number of bubbles
                    end if
                end do
                close (94)
            else
                call s_mpi_abort("Initialize the lagrange bubbles in " // trim(lag_params%input_path))
            end if
        else
            if (proc_rank == 0) print *, 'Restarting lagrange bubbles at save_count: ', save_count
            call s_restart_bubbles(bub_id, save_count)
        end if

        print *, " Lagrange bubbles running, in proc", proc_rank, "number:", bub_id, "/", id

        if (num_procs > 1) then
            call s_mpi_reduce_int_sum(n_el_bubs_loc, n_el_bubs_glb)
        else
            n_el_bubs_glb = n_el_bubs_loc
        end if

        if (proc_rank == 0) then
            if (n_el_bubs_glb == 0) call s_mpi_abort('No bubbles in the domain. Check ' // trim(lag_params%input_path))
        end if

        $:GPU_UPDATE(device='[bubbles_lagrange, lag_params]')

        $:GPU_UPDATE(device='[lag_id, bub_R0, Rmax_stats, Rmin_stats, gas_mg, gas_betaT, gas_betaC, bub_dphidt, gas_p, gas_mv, &
                     & intfc_rad, intfc_vel, mtn_pos, mtn_posPrev, mtn_vel, mtn_s, intfc_draddt, intfc_dveldt, gas_dpdt, &
                     & gas_dmvdt, mtn_dposdt, mtn_dveldt, n_el_bubs_loc]')

        Rmax_glb = min(dflt_real, -dflt_real)
        Rmin_glb = max(dflt_real, -dflt_real)
        $:GPU_UPDATE(device='[Rmax_glb, Rmin_glb]')

        $:GPU_UPDATE(device='[dx, dy, dz, x_cb, x_cc, y_cb, y_cc, z_cb, z_cc]')

        ! Populate temporal variables
        call s_transfer_data_to_tmp()
        call s_smear_voidfraction(bc_type)

        if (save_count == 0) then
            ! Create ./D directory
            if (proc_rank == 0) then
                write (path_D_dir, '(A,I0,A,I0)') trim(case_dir) // '/D'
                call my_inquire(trim(path_D_dir), file_exist)
                if (.not. file_exist) call s_create_directory(trim(path_D_dir))
            end if
            call s_mpi_barrier()
            call s_write_restart_lag_bubbles(save_count)  ! Needed for post_processing
            if (lag_params%write_void_evol) call s_write_void_evol(qtime)
        end if

        if (lag_params%write_bubbles) call s_write_lag_bubble_evol(qtime)

    end subroutine s_read_input_bubbles

    !> Add a new bubble from input data for a fresh start
    impure subroutine s_add_bubbles(inputBubble, q_cons_vf, bub_id)

        type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf
        real(wp), dimension(8), intent(in)                  :: inputBubble
        integer, intent(in)                                 :: bub_id
        integer                                             :: i
        real(wp)                                            :: pliq, volparticle, concvap, totalmass, kparticle, cpparticle
        real(wp)                                            :: omegaN_local, PeG, PeT, rhol, qv, gamma, pi_inf, dynP
        integer, dimension(3)                               :: cell
        real(wp), dimension(2)                              :: Re
        real(wp)                                            :: massflag, heatflag, Re_trans, Im_trans

        massflag = 0._wp
        heatflag = 0._wp
        if (lag_params%massTransfer_model) massflag = 1._wp
        if (lag_params%heatTransfer_model) heatflag = 1._wp

        bub_R0(bub_id) = inputBubble(7)
        Rmax_stats(bub_id) = min(dflt_real, -dflt_real)
        Rmin_stats(bub_id) = max(dflt_real, -dflt_real)
        bub_dphidt(bub_id) = 0._wp
        intfc_rad(bub_id, 1) = inputBubble(7)
        intfc_vel(bub_id, 1) = inputBubble(8)
        mtn_pos(bub_id,1:3,1) = inputBubble(1:3)
        mtn_posPrev(bub_id,1:3,1) = mtn_pos(bub_id,1:3,1)
        mtn_vel(bub_id,1:3,1) = inputBubble(4:6)

        if (cyl_coord .and. p == 0) then
            mtn_pos(bub_id, 2, 1) = sqrt(mtn_pos(bub_id, 2, 1)**2._wp + mtn_pos(bub_id, 3, 1)**2._wp)
            ! Storing azimuthal angle (-Pi to Pi)) into the third coordinate variable
            mtn_pos(bub_id, 3, 1) = atan2(inputBubble(3), inputBubble(2))
            mtn_posPrev(bub_id,1:3,1) = mtn_pos(bub_id,1:3,1)
        end if

        cell = fd_number - buff_size
        call s_locate_cell(mtn_pos(bub_id,1:3,1), cell, mtn_s(bub_id,1:3,1))

        ! Check if the bubble is located in the ghost cell of a symmetric, or wall boundary
        if ((any(bc_x%beg == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, &
            & BC_NO_SLIP_WALL/)) .and. cell(1) < 0) .or. (any(bc_x%end == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, &
            & BC_NO_SLIP_WALL/)) .and. cell(1) > m) .or. (any(bc_y%beg == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, &
            & BC_NO_SLIP_WALL/)) .and. cell(2) < 0) .or. (any(bc_y%end == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, &
            & BC_NO_SLIP_WALL/)) .and. cell(2) > n)) then
            call s_mpi_abort("Lagrange bubble is in the ghost cells of a symmetric or wall boundary.")
        end if

        if (p > 0) then
            if ((any(bc_z%beg == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, &
                & BC_NO_SLIP_WALL/)) .and. cell(3) < 0) .or. (any(bc_z%end == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, &
                & BC_NO_SLIP_WALL/)) .and. cell(3) > p)) then
                call s_mpi_abort("Lagrange bubble is in the ghost cells of a symmetric or wall boundary.")
            end if
        end if

        call s_convert_to_mixture_variables(q_cons_vf, cell(1), cell(2), cell(3), rhol, gamma, pi_inf, qv, Re)
        dynP = 0._wp
        do i = 1, num_dims
            dynP = dynP + 0.5_wp*q_cons_vf(eqn_idx%cont%end + i)%sf(cell(1), cell(2), cell(3))**2/rhol
        end do
        pliq = f_pressure(q_cons_vf(eqn_idx%E)%sf(cell(1), cell(2), cell(3)) - dynP, gamma, pi_inf, qv)
        if (pliq < 0) print *, "Negative pressure", proc_rank, q_cons_vf(eqn_idx%E)%sf(cell(1), cell(2), cell(3)), pi_inf, gamma, &
            & pliq, cell, dynP

        ! Initial particle pressure
        gas_p(bub_id, 1) = pliq + 2._wp*(1._wp/Web)/bub_R0(bub_id)

        ! Initial particle mass
        volparticle = 4._wp/3._wp*pi*bub_R0(bub_id)**3._wp  ! volume
        gas_mv(bub_id, 1) = pv*volparticle*(1._wp/(R_v*Tw))*(massflag)  ! vapermass
        gas_mg(bub_id) = (gas_p(bub_id, 1) - pv*(massflag))*volparticle*(1._wp/(R_g*Tw))  ! gasmass
        if (gas_mg(bub_id) <= 0._wp) then
            call s_mpi_abort("The initial mass of gas inside the bubble is negative. Check the initial conditions.")
        end if
        totalmass = gas_mg(bub_id) + gas_mv(bub_id, 1)  ! totalmass

        ! Bubble natural frequency
        concvap = gas_mv(bub_id, 1)/(gas_mv(bub_id, 1) + gas_mg(bub_id))
        omegaN_local = (3._wp*(gas_p(bub_id, 1) - pv*(massflag)) + 4._wp*(1._wp/Web)/bub_R0(bub_id))/rhol
        if (pv*(massflag) > gas_p(bub_id, 1)) then
            call s_mpi_abort("Lagrange bubble initially located in a region with pressure below the vapor pressure.")
        end if
        omegaN_local = sqrt(omegaN_local/bub_R0(bub_id)**2._wp)

        cpparticle = concvap*cp_v + (1._wp - concvap)*cp_g
        kparticle = concvap*k_vl + (1._wp - concvap)*k_gl

        ! Mass and heat transfer coefficients (based on Preston 2007)
        PeT = totalmass/volparticle*cpparticle*bub_R0(bub_id)**2._wp*omegaN_local/kparticle
        call s_transcoeff(1._wp, PeT, Re_trans, Im_trans)
        gas_betaT(bub_id) = Re_trans*(heatflag)*kparticle

        PeG = bub_R0(bub_id)**2._wp*omegaN_local/vd
        call s_transcoeff(1._wp, PeG, Re_trans, Im_trans)
        gas_betaC(bub_id) = Re_trans*(massflag)*vd

        if (gas_mg(bub_id) <= 0._wp) then
            call s_mpi_abort("Negative gas mass in the bubble, check if the bubble is in the domain.")
        end if

    end subroutine s_add_bubbles

    !> Restore bubble data from a restart file
    impure subroutine s_restart_bubbles(bub_id, save_count)

        integer, intent(inout)               :: bub_id, save_count
        character(LEN=path_len + 2*name_len) :: file_loc
        real(wp)                             :: file_time, file_dt
        integer                              :: file_num_procs, file_tot_part, tot_part

#ifdef MFC_MPI
        real(wp), dimension(20)                :: inputvals
        integer, dimension(MPI_STATUS_SIZE)    :: status
        integer(kind=MPI_OFFSET_KIND)          :: disp
        integer                                :: view
        integer, dimension(3)                  :: cell
        logical                                :: indomain, particle_file, file_exist
        integer, dimension(2)                  :: gsizes, lsizes, start_idx_part
        integer                                :: ifile, ierr, tot_data, id
        integer                                :: i
        integer, dimension(:), allocatable     :: proc_bubble_counts
        real(wp), dimension(1:1,1:lag_io_vars) :: dummy

        dummy = 0._wp

        ! Construct file path
        write (file_loc, '(A,I0,A)') 'lag_bubbles_', save_count, '.dat'
        file_loc = trim(case_dir) // '/restart_data' // trim(mpiiofs) // trim(file_loc)

        ! Check if file exists
        inquire (FILE=trim(file_loc), EXIST=file_exist)
        if (.not. file_exist) then
            call s_mpi_abort('Restart file ' // trim(file_loc) // ' does not exist!')
        end if

        if (.not. parallel_io) return

        if (proc_rank == 0) then
            call MPI_FILE_OPEN(MPI_COMM_SELF, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr)

            call MPI_FILE_READ(ifile, file_tot_part, 1, MPI_INTEGER, status, ierr)
            call MPI_FILE_READ(ifile, file_time, 1, mpi_p, status, ierr)
            call MPI_FILE_READ(ifile, file_dt, 1, mpi_p, status, ierr)
            call MPI_FILE_READ(ifile, file_num_procs, 1, MPI_INTEGER, status, ierr)

            call MPI_FILE_CLOSE(ifile, ierr)
        end if

        call MPI_BCAST(file_tot_part, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
        call MPI_BCAST(file_time, 1, mpi_p, 0, MPI_COMM_WORLD, ierr)
        call MPI_BCAST(file_dt, 1, mpi_p, 0, MPI_COMM_WORLD, ierr)
        call MPI_BCAST(file_num_procs, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)

        allocate (proc_bubble_counts(file_num_procs))

        if (proc_rank == 0) then
            call MPI_FILE_OPEN(MPI_COMM_SELF, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr)

            ! Skip to processor counts position
            disp = int(sizeof(file_tot_part) + 2*sizeof(file_time) + sizeof(file_num_procs), MPI_OFFSET_KIND)
            call MPI_FILE_SEEK(ifile, disp, MPI_SEEK_SET, ierr)
            call MPI_FILE_READ(ifile, proc_bubble_counts, file_num_procs, MPI_INTEGER, status, ierr)

            call MPI_FILE_CLOSE(ifile, ierr)
        end if

        call MPI_BCAST(proc_bubble_counts, file_num_procs, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)

        ! Set time variables from file
        mytime = file_time
        dt = file_dt

        bub_id = proc_bubble_counts(proc_rank + 1)

        start_idx_part(1) = 0
        do i = 1, proc_rank
            start_idx_part(1) = start_idx_part(1) + proc_bubble_counts(i)
        end do

        start_idx_part(2) = 0
        lsizes(1) = bub_id
        lsizes(2) = lag_io_vars

        gsizes(1) = file_tot_part
        gsizes(2) = lag_io_vars

        if (bub_id > 0) then
            allocate (MPI_IO_DATA_lag_bubbles(bub_id,1:lag_io_vars))

            call MPI_TYPE_CREATE_SUBARRAY(2, gsizes, lsizes, start_idx_part, MPI_ORDER_FORTRAN, mpi_p, view, ierr)
            call MPI_TYPE_COMMIT(view, ierr)

            call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr)

            ! Skip extended header
            disp = int(sizeof(file_tot_part) + 2*sizeof(file_time) + sizeof(file_num_procs) &
                       & + file_num_procs*sizeof(proc_bubble_counts(1)), MPI_OFFSET_KIND)
            call MPI_FILE_SET_VIEW(ifile, disp, mpi_p, view, 'native', mpi_info_int, ierr)

            call MPI_FILE_READ_ALL(ifile, MPI_IO_DATA_lag_bubbles, lag_io_vars*bub_id, mpi_p, status, ierr)

            call MPI_FILE_CLOSE(ifile, ierr)
            call MPI_TYPE_FREE(view, ierr)

            n_el_bubs_loc = bub_id

            do i = 1, bub_id
                lag_id(i, 1) = int(MPI_IO_DATA_lag_bubbles(i, 1))
                mtn_pos(i,1:3,1) = MPI_IO_DATA_lag_bubbles(i,2:4)
                mtn_posPrev(i,1:3,1) = MPI_IO_DATA_lag_bubbles(i,5:7)
                mtn_vel(i,1:3,1) = MPI_IO_DATA_lag_bubbles(i,8:10)
                intfc_rad(i, 1) = MPI_IO_DATA_lag_bubbles(i, 11)
                intfc_vel(i, 1) = MPI_IO_DATA_lag_bubbles(i, 12)
                bub_R0(i) = MPI_IO_DATA_lag_bubbles(i, 13)
                Rmax_stats(i) = MPI_IO_DATA_lag_bubbles(i, 14)
                Rmin_stats(i) = MPI_IO_DATA_lag_bubbles(i, 15)
                bub_dphidt(i) = MPI_IO_DATA_lag_bubbles(i, 16)
                gas_p(i, 1) = MPI_IO_DATA_lag_bubbles(i, 17)
                gas_mv(i, 1) = MPI_IO_DATA_lag_bubbles(i, 18)
                gas_mg(i) = MPI_IO_DATA_lag_bubbles(i, 19)
                gas_betaT(i) = MPI_IO_DATA_lag_bubbles(i, 20)
                gas_betaC(i) = MPI_IO_DATA_lag_bubbles(i, 21)
                cell = -buff_size
                call s_locate_cell(mtn_pos(i,1:3,1), cell, mtn_s(i,1:3,1))
            end do

            deallocate (MPI_IO_DATA_lag_bubbles)
        else
            n_el_bubs_loc = 0

            call MPI_TYPE_CONTIGUOUS(0, mpi_p, view, ierr)
            call MPI_TYPE_COMMIT(view, ierr)

            call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr)

            ! Skip extended header
            disp = int(sizeof(file_tot_part) + 2*sizeof(file_time) + sizeof(file_num_procs) &
                       & + file_num_procs*sizeof(proc_bubble_counts(1)), MPI_OFFSET_KIND)
            call MPI_FILE_SET_VIEW(ifile, disp, mpi_p, view, 'native', mpi_info_int, ierr)

            call MPI_FILE_READ_ALL(ifile, dummy, 0, mpi_p, status, ierr)

            call MPI_FILE_CLOSE(ifile, ierr)
            call MPI_TYPE_FREE(view, ierr)
        end if

        if (proc_rank == 0) then
            write (*, '(A,I0,A,I0)') 'Read ', file_tot_part, ' particles from restart file at t_step = ', save_count
            write (*, '(A,E15.7,A,E15.7)') 'Restart time = ', mytime, ', dt = ', dt
        end if

        deallocate (proc_bubble_counts)
#endif

    end subroutine s_restart_bubbles

    !> Contains the bubble dynamics subroutines.
    subroutine s_compute_bubble_EL_dynamics(q_prim_vf, bc_type, stage)

        type(scalar_field), dimension(sys_size), intent(inout)     :: q_prim_vf
        type(integer_field), dimension(1:num_dims,1:2), intent(in) :: bc_type
        integer, intent(in)                                        :: stage
        real(wp)                                                   :: myVapFlux
        real(wp)                                                   :: preterm1, term2, paux, pint, Romega, term1_fac
        real(wp)                                                   :: myR_m, mygamma_m, myPb, myMass_n, myMass_v
        real(wp)                                                   :: myR, myV, myBeta_c, myBeta_t, myR0, myPbdot, myMvdot
        real(wp)                                                   :: myPinf, aux1, aux2, myCson, myRho
        real(wp), dimension(3)                                     :: myPos, myVel
        real(wp)                                                   :: gamma, pi_inf, qv, f_b, myRe

        #:if not MFC_CASE_OPTIMIZATION and USING_AMD
            real(wp), dimension(3) :: myalpha_rho, myalpha
        #:else
            real(wp), dimension(num_fluids) :: myalpha_rho, myalpha
        #:endif
        real(wp), dimension(2) :: Re
        integer, dimension(3) :: cell
        integer :: adap_dt_stop_sum, adap_dt_stop                     !< Fail-safe exit if max iteration count reached
        real(wp) :: dmalf, dmntait, dmBtait, dm_bub_adv_src, dm_divu  !< Dummy variables for unified subgrid bubble subroutines
        integer :: k, l

        ! Subgrid p_inf model based on Maeda and Colonius (2018).
        if (lag_params%pressure_corrector) then
            call nvtxStartRange("LAGRANGE-BUBBLE-PINF-CORRECTION")
            ! Calculate velocity potentials (valid for one bubble per cell)
            $:GPU_PARALLEL_LOOP(private='[k, cell, paux, preterm1, term2, Romega, myR0, myR, myV, myPb, pint, term1_fac]')
            do k = 1, n_el_bubs_loc
                call s_get_pinf(k, q_prim_vf, 2, paux, cell, preterm1, term2, Romega)
                myR0 = bub_R0(k)
                myR = intfc_rad(k, 2)
                myV = intfc_vel(k, 2)
                myPb = gas_p(k, 2)
                pint = f_cpbw_KM(myR0, myR, myV, myPb)
                pint = pint + 0.5_wp*myV**2._wp
                if (lag_params%cluster_type == 2) then
                    bub_dphidt(k) = (paux - pint) + term2
                    ! Accounting for the potential induced by the bubble averaged over the control volume Note that this is based on
                    ! the incompressible flow assumption near the bubble.
                    term1_fac = 3._wp/2._wp*(myR*(Romega**2._wp - myR**2._wp))/(Romega**3._wp - myR**3._wp)
                    bub_dphidt(k) = bub_dphidt(k)/(1._wp - term1_fac)
                end if
            end do
            $:END_GPU_PARALLEL_LOOP()
            call nvtxEndRange()
        end if

        ! Precompute cell-centered pressure gradients for translational motion
        if (moving_lag_bubbles .and. lag_pressure_force) then
            call nvtxStartRange("LAGRANGE-BUBBLE-PRESSURE-GRADIENT")
            call s_compute_pressure_gradients(q_prim_vf)
            call nvtxEndRange()
        end if

        call nvtxStartRange("LAGRANGE-BUBBLE-DYNAMICS")
        ! Radial motion model
        adap_dt_stop_sum = 0
        $:GPU_PARALLEL_LOOP(private='[k, myalpha_rho, myalpha, Re, cell, myVapFlux, preterm1, term2, paux, pint, Romega, &
                            & term1_fac, myR_m, mygamma_m, myPb, myMass_n, myMass_v, myR, myV, myBeta_c, myBeta_t, myR0, myPbdot, &
                            & myMvdot, myPinf, aux1, aux2, myCson, myRho, gamma, pi_inf, qv, dmalf, dmntait, dmBtait, &
                            & dm_bub_adv_src, dm_divu, adap_dt_stop, myPos, myVel]', copy='[adap_dt_stop_sum]',copyin='[stage]')
        do k = 1, n_el_bubs_loc
            ! Keller-Miksis model

            ! Current bubble state
            myPb = gas_p(k, 2)
            myMass_n = gas_mg(k)
            myMass_v = gas_mv(k, 2)
            myR = intfc_rad(k, 2)
            myV = intfc_vel(k, 2)
            myBeta_c = gas_betaC(k)
            myBeta_t = gas_betaT(k)
            myR0 = bub_R0(k)
            myPos = mtn_pos(k,:,2)
            myVel = mtn_vel(k,:,2)

            ! Vapor and heat fluxes
            call s_vflux(myR, myV, myPb, myMass_v, k, myVapFlux, myMass_n, myBeta_c, myR_m, mygamma_m)
            myPbdot = f_bpres_dot(myVapFlux, myR, myV, myPb, myMass_v, k, myBeta_t, myR_m, mygamma_m)
            myMvdot = 4._wp*pi*myR**2._wp*myVapFlux

            ! Obtaining driving pressure
            call s_get_pinf(k, q_prim_vf, 1, myPinf, cell, aux1, aux2)

            ! Obtain liquid density and computing speed of sound from pinf
            call s_compute_species_fraction(q_prim_vf, cell(1), cell(2), cell(3), myalpha_rho, myalpha)
            call s_convert_species_to_mixture_variables_kernel(myRho, gamma, pi_inf, qv, myalpha, myalpha_rho, Re)
            myCson = sqrt(f_bulk_modulus(myPinf, gamma, pi_inf)/myRho)

            ! Adaptive time stepping
            adap_dt_stop = 0

            if (adap_dt) then
                mtn_posPrev(k,:,1) = myPos

                myRe = Re(1)
                adap_dt_stop = f_advance_step(myRho, myPinf, myR, myV, myR0, myPb, myPbdot, dmalf, dmntait, dmBtait, &
                                              & dm_bub_adv_src, dm_divu, k, myMass_v, myMass_n, myBeta_c, myBeta_t, myCson, myRe, &
                                              & myPos, myVel, cell, q_prim_vf)

                ! Update bubble state
                intfc_rad(k, 1) = myR
                intfc_vel(k, 1) = myV
                gas_p(k, 1) = myPb
                gas_mv(k, 1) = myMass_v
                mtn_pos(k,:,1) = myPos
                mtn_vel(k,:,1) = myVel
            else
                ! Radial acceleration from bubble models
                intfc_dveldt(k, stage) = f_rddot(myRho, myPinf, myR, myV, myR0, myPb, myPbdot, dmalf, dmntait, dmBtait, &
                             & dm_bub_adv_src, dm_divu, myCson)

                intfc_draddt(k, stage) = myV
                gas_dmvdt(k, stage) = myMvdot
                gas_dpdt(k, stage) = myPbdot

                if (moving_lag_bubbles) then
                    do l = 1, num_dims
                        mtn_dposdt(k, l, stage) = 0._wp
                        mtn_dveldt(k, l, stage) = 0._wp
                        select case (lag_vel_model)
                        case (1)
                            mtn_dposdt(k, l, stage) = f_interpolate_velocity(myPos(l), cell, l, q_prim_vf)
                            mtn_dveldt(k, l, stage) = 0._wp
                            mtn_vel(k, l, 1) = mtn_dposdt(k, l, stage)
                        case (2)
                            mtn_dposdt(k, l, stage) = myVel(l)
                            f_b = f_get_bubble_force(myPos(l), myR, myV, myVel(l), myMass_n, myMass_v, Re(1), myRho, cell, l, &
                                                     & q_prim_vf)
                            mtn_dveldt(k, l, stage) = f_b/(myMass_n + myMass_v)
                        case default
                            mtn_dposdt(k, l, stage) = 0._wp
                            mtn_dveldt(k, l, stage) = 0._wp
                        end select
                    end do
                end if
            end if

            $:GPU_ATOMIC(atomic='update')
            adap_dt_stop_sum = adap_dt_stop_sum + adap_dt_stop
        end do
        $:END_GPU_PARALLEL_LOOP()
        call nvtxEndRange

        if (adap_dt .and. adap_dt_stop_sum > 0) call s_mpi_abort("Adaptive time stepping failed to converge.")

        if (adap_dt) then
            call s_transfer_data_to_tmp()
            if (moving_lag_bubbles) call s_enforce_EL_bubbles_boundary_conditions(q_prim_vf)
            call s_smear_voidfraction(bc_type)
        end if

    end subroutine s_compute_bubble_EL_dynamics

    !> Compute the Lagrangian bubble source terms and add them to the RHS
    subroutine s_compute_bubbles_EL_source(q_cons_vf, q_prim_vf, rhs_vf)

        type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
        type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
        type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf
        integer                                                :: i, j, k, l

        call nvtxStartRange("LAGRANGE-BUBBLE-EL-SOURCE")
        ! (q / (1 - beta)) * d(beta)/dt source
        if (lag_params%cluster_type >= 4) then
            $:GPU_PARALLEL_LOOP(private='[i, j, k, l]', collapse=4)
            do k = idwint(3)%beg, idwint(3)%end
                do j = idwint(2)%beg, idwint(2)%end
                    do i = idwint(1)%beg, idwint(1)%end
                        do l = 1, eqn_idx%E
                            if (q_beta(1)%sf(i, j, k) > (1._wp - lag_params%valmaxvoid)) then
                                rhs_vf(l)%sf(i, j, k) = rhs_vf(l)%sf(i, j, k) + q_cons_vf(l)%sf(i, j, k)*(q_beta(2)%sf(i, j, &
                                       & k) + q_beta(5)%sf(i, j, k))
                            end if
                        end do
                    end do
                end do
            end do
            $:END_GPU_PARALLEL_LOOP()
        else
            $:GPU_PARALLEL_LOOP(private='[i, j, k, l]', collapse=4)
            do k = idwint(3)%beg, idwint(3)%end
                do j = idwint(2)%beg, idwint(2)%end
                    do i = idwint(1)%beg, idwint(1)%end
                        do l = 1, eqn_idx%E
                            if (q_beta(1)%sf(i, j, k) > (1._wp - lag_params%valmaxvoid)) then
                                rhs_vf(l)%sf(i, j, k) = rhs_vf(l)%sf(i, j, k) + (q_cons_vf(l)%sf(i, j, k)/q_beta(1)%sf(i, j, &
                                       & k))*q_beta(2)%sf(i, j, k)
                            end if
                        end do
                    end do
                end do
            end do
            $:END_GPU_PARALLEL_LOOP()
        end if

        do l = 1, num_dims
            call s_gradient_dir(q_prim_vf(eqn_idx%E)%sf, q_beta(3)%sf, l)

            ! (q / (1 - beta)) * d(beta)/dt source
            $:GPU_PARALLEL_LOOP(private='[i, j, k]', collapse=3)
            do k = idwint(3)%beg, idwint(3)%end
                do j = idwint(2)%beg, idwint(2)%end
                    do i = idwint(1)%beg, idwint(1)%end
                        if (q_beta(1)%sf(i, j, k) > (1._wp - lag_params%valmaxvoid)) then
                            rhs_vf(eqn_idx%cont%end + l)%sf(i, j, k) = rhs_vf(eqn_idx%cont%end + l)%sf(i, j, &
                                   & k) - (1._wp - q_beta(1)%sf(i, j, k))/q_beta(1)%sf(i, j, k)*q_beta(3)%sf(i, j, k)
                        end if
                    end do
                end do
            end do
            $:END_GPU_PARALLEL_LOOP()

            ! source in energy
            $:GPU_PARALLEL_LOOP(private='[i, j, k]', collapse=3)
            do k = idwbuff(3)%beg, idwbuff(3)%end
                do j = idwbuff(2)%beg, idwbuff(2)%end
                    do i = idwbuff(1)%beg, idwbuff(1)%end
                        q_beta(3)%sf(i, j, k) = q_prim_vf(eqn_idx%E)%sf(i, j, k)*q_prim_vf(eqn_idx%cont%end + l)%sf(i, j, k)
                    end do
                end do
            end do
            $:END_GPU_PARALLEL_LOOP()

            call s_gradient_dir(q_beta(3)%sf, q_beta(4)%sf, l)

            ! (beta / (1 - beta)) * d(Pu)/dl source
            $:GPU_PARALLEL_LOOP(private='[i, j, k]', collapse=3)
            do k = idwint(3)%beg, idwint(3)%end
                do j = idwint(2)%beg, idwint(2)%end
                    do i = idwint(1)%beg, idwint(1)%end
                        if (q_beta(1)%sf(i, j, k) > (1._wp - lag_params%valmaxvoid)) then
                            rhs_vf(eqn_idx%E)%sf(i, j, k) = rhs_vf(eqn_idx%E)%sf(i, j, k) - q_beta(4)%sf(i, j, &
                                   & k)*(1._wp - q_beta(1)%sf(i, j, k))/q_beta(1)%sf(i, j, k)
                        end if
                    end do
                end do
            end do
            $:END_GPU_PARALLEL_LOOP()
        end do
        call nvtxEndRange

    end subroutine s_compute_bubbles_EL_source

    !> Smear the bubble effects onto the Eulerian grid
    subroutine s_smear_voidfraction(bc_type)

        type(integer_field), dimension(1:num_dims,1:2), intent(in) :: bc_type
        integer                                                    :: i, j, k, l

        call nvtxStartRange("BUBBLES-LAGRANGE-SMEARING")
        $:GPU_PARALLEL_LOOP(private='[i, j, k, l]', collapse=4)
        do i = 1, q_beta_idx
            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_beta(i)%sf(j, k, l) = 0._wp
                        if (lag_params%kahan_summation) kahan_comp(i)%sf(j, k, l) = 0._wp
                    end do
                end do
            end do
        end do
        $:END_GPU_PARALLEL_LOOP()

        ! Build cell list before smearing (CPU-side counting sort)
        call s_build_cell_list(n_el_bubs_loc, mtn_s)

        call s_smoothfunction(n_el_bubs_loc, intfc_rad, intfc_vel, mtn_s, mtn_pos, q_beta, kahan_comp)

        call nvtxStartRange("BUBBLES-LAGRANGE-BETA-COMM")
        if (lag_params%cluster_type >= 4) then
            call s_populate_beta_buffers(q_beta, kahan_comp, bc_type, 3)
        else
            call s_populate_beta_buffers(q_beta, kahan_comp, bc_type, 2)
        end if
        call nvtxEndRange

        ! Store 1-beta
        $:GPU_PARALLEL_LOOP(private='[j, k, l]', collapse=3)
        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_beta(1)%sf(j, k, l) = 1._wp - q_beta(1)%sf(j, k, l)
                    ! Limiting void fraction given max value
                    q_beta(1)%sf(j, k, l) = max(q_beta(1)%sf(j, k, l), 1._wp - lag_params%valmaxvoid)
                end do
            end do
        end do
        $:END_GPU_PARALLEL_LOOP()
        call nvtxEndRange  ! BUBBLES-LAGRANGE-SMEARING

    end subroutine s_smear_voidfraction

    !> Compute the bubble driving pressure p_inf
    subroutine s_get_pinf(bub_id, q_prim_vf, ptype, f_pinfl, cell, preterm1, term2, Romega)

        $:GPU_ROUTINE(function_name='s_get_pinf',parallelism='[seq]', cray_inline=True)

        integer, intent(in)                                 :: bub_id, ptype
        type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
        real(wp), intent(out)                               :: f_pinfl
        integer, dimension(3), intent(out)                  :: cell
        real(wp), intent(out), optional                     :: preterm1, term2, Romega
        real(wp), dimension(3)                              :: scoord, psi_pos, psi_x, psi_y, psi_z
        real(wp)                                            :: xi, eta, zeta
        real(wp)                                            :: dc, vol, aux
        real(wp)                                            :: volgas, term1, Rbeq, denom
        real(wp)                                            :: charvol, charpres, charvol2, charpres2
        integer, dimension(3)                               :: cellaux
        integer                                             :: i, j, k
        integer                                             :: smearGrid, smearGridz

        f_pinfl = 0._wp

        if (moving_lag_bubbles) then
            cell = fd_number - buff_size
            call s_locate_cell(mtn_pos(bub_id,1:3,2), cell, mtn_s(bub_id,1:3,2))
            scoord = mtn_s(bub_id,1:3,2)
        else
            scoord = mtn_s(bub_id,1:3,2)
            cell(:) = int(scoord(:))
            $:GPU_LOOP(parallelism='[seq]')
            do i = 1, num_dims
                if (scoord(i) < 0._wp) cell(i) = cell(i) - 1
            end do
        end if

        if ((lag_params%cluster_type == 1)) then
            !> Getting p_cell in terms of only the current cell by interpolation

            if (fd_order == 2) then  ! Bilinear interpolation
                if (p > 0) then
                    vol = dx(cell(1))*dy(cell(2))*dz(cell(3))
                else
                    if (cyl_coord) then
                        vol = dx(cell(1))*dy(cell(2))*y_cc(cell(2))*2._wp*pi
                    else
                        vol = dx(cell(1))*dy(cell(2))*lag_params%charwidth
                    end if
                end if

                !> Obtain bilinear interpolation coefficients, based on the current location of the bubble.
                psi_pos(1) = (scoord(1) - real(cell(1)))*dx(cell(1)) + x_cb(cell(1) - 1)
                psi_pos(1) = abs((psi_pos(1) - x_cc(cell(1)))/(x_cc(cell(1) + 1) - x_cc(cell(1))))

                psi_pos(2) = (scoord(2) - real(cell(2)))*dy(cell(2)) + y_cb(cell(2) - 1)
                psi_pos(2) = abs((psi_pos(2) - y_cc(cell(2)))/(y_cc(cell(2) + 1) - y_cc(cell(2))))

                if (p > 0) then
                    psi_pos(3) = (scoord(3) - real(cell(3)))*dz(cell(3)) + z_cb(cell(3) - 1)
                    psi_pos(3) = abs((psi_pos(3) - z_cc(cell(3)))/(z_cc(cell(3) + 1) - z_cc(cell(3))))
                else
                    psi_pos(3) = 0._wp
                end if

                ! Calculate bilinear basis functions for each direction For normalized coordinate xi in [0, 1], the two basis
                ! functions are: phi_0(xi) = 1 - xi, phi_1(xi) = xi

                ! X-direction basis functions
                psi_x(1) = 1._wp - psi_pos(1)  ! Left basis function
                psi_x(2) = psi_pos(1)  ! Right basis function

                ! Y-direction basis functions
                psi_y(1) = 1._wp - psi_pos(2)  ! Left basis function
                psi_y(2) = psi_pos(2)  ! Right basis function

                if (p > 0) then
                    ! Z-direction basis functions
                    psi_z(1) = 1._wp - psi_pos(3)  ! Left basis function
                    psi_z(2) = psi_pos(3)  ! Right basis function
                else  ! 3D
                    psi_z(1) = 1._wp
                    psi_z(2) = 0._wp
                end if

                !> Perform bilinear interpolation
                f_pinfl = 0._wp

                if (p == 0) then
                    do j = 1, 2
                        do i = 1, 2
                            f_pinfl = f_pinfl + q_prim_vf(eqn_idx%E)%sf(cell(1) + i - 1, cell(2) + j - 1, cell(3))*psi_x(i)*psi_y(j)
                        end do
                    end do
                else
                    do k = 1, 2
                        do j = 1, 2
                            do i = 1, 2
                                f_pinfl = f_pinfl + q_prim_vf(eqn_idx%E)%sf(cell(1) + i - 1, cell(2) + j - 1, &
                                                              & cell(3) + k - 1)*psi_x(i)*psi_y(j)*psi_z(k)
                            end do
                        end do
                    end do
                end if
            else if (fd_order == 4) then  ! Biquadratic interpolation
                if (p > 0) then
                    vol = dx(cell(1))*dy(cell(2))*dz(cell(3))
                else
                    if (cyl_coord) then
                        vol = dx(cell(1))*dy(cell(2))*y_cc(cell(2))*2._wp*pi
                    else
                        vol = dx(cell(1))*dy(cell(2))*lag_params%charwidth
                    end if
                end if

                !> Obtain biquadratic interpolation coefficients, based on the current location of the bubble.
                ! For biquadratic interpolation, we need coefficients for 3 points in each direction
                psi_pos(1) = (scoord(1) - real(cell(1)))*dx(cell(1)) + x_cb(cell(1) - 1)
                psi_pos(1) = (psi_pos(1) - x_cc(cell(1)))/(x_cc(cell(1) + 1) - x_cc(cell(1)))

                psi_pos(2) = (scoord(2) - real(cell(2)))*dy(cell(2)) + y_cb(cell(2) - 1)
                psi_pos(2) = (psi_pos(2) - y_cc(cell(2)))/(y_cc(cell(2) + 1) - y_cc(cell(2)))

                if (p > 0) then
                    psi_pos(3) = (scoord(3) - real(cell(3)))*dz(cell(3)) + z_cb(cell(3) - 1)
                    psi_pos(3) = (psi_pos(3) - z_cc(cell(3)))/(z_cc(cell(3) + 1) - z_cc(cell(3)))
                else
                    psi_pos(3) = 0._wp
                end if

                ! Calculate biquadratic basis functions for each direction For normalized coordinate xi in [-1, 1], the three basis
                ! functions are: phi_0(xi) = xi*(xi-1)/2, phi_1(xi) = (1-xi)*(1+xi), phi_2(xi) = xi*(xi+1)/2

                ! X-direction basis functions
                xi = 2._wp*psi_pos(1) - 1._wp  ! Convert to [-1, 1] range
                psi_x(1) = xi*(xi - 1._wp)/2._wp  ! Left basis function
                psi_x(2) = (1._wp - xi)*(1._wp + xi)  ! Center basis function
                psi_x(3) = xi*(xi + 1._wp)/2._wp  ! Right basis function

                ! Y-direction basis functions
                eta = 2._wp*psi_pos(2) - 1._wp  ! Convert to [-1, 1] range
                psi_y(1) = eta*(eta - 1._wp)/2._wp  ! Left basis function
                psi_y(2) = (1._wp - eta)*(1._wp + eta)  ! Center basis function
                psi_y(3) = eta*(eta + 1._wp)/2._wp  ! Right basis function

                if (p > 0) then
                    ! Z-direction basis functions
                    zeta = 2._wp*psi_pos(3) - 1._wp  ! Convert to [-1, 1] range
                    psi_z(1) = zeta*(zeta - 1._wp)/2._wp  ! Left basis function
                    psi_z(2) = (1._wp - zeta)*(1._wp + zeta)  ! Center basis function
                    psi_z(3) = zeta*(zeta + 1._wp)/2._wp  ! Right basis function
                else
                    psi_z(1) = 0._wp
                    psi_z(2) = 1._wp
                    psi_z(3) = 0._wp
                end if

                !> Perform biquadratic interpolation
                f_pinfl = 0._wp

                if (p == 0) then  ! 2D
                    do j = 1, 3
                        do i = 1, 3
                            f_pinfl = f_pinfl + q_prim_vf(eqn_idx%E)%sf(cell(1) + i - 2, cell(2) + j - 2, cell(3))*psi_x(i)*psi_y(j)
                        end do
                    end do
                else
                    do k = 1, 3
                        do j = 1, 3
                            do i = 1, 3
                                f_pinfl = f_pinfl + q_prim_vf(eqn_idx%E)%sf(cell(1) + i - 2, cell(2) + j - 2, &
                                                              & cell(3) + k - 2)*psi_x(i)*psi_y(j)*psi_z(k)
                            end do
                        end do
                    end do
                end if
            end if

            ! R_Omega
            dc = (3._wp*vol/(4._wp*pi))**(1._wp/3._wp)
        else if (lag_params%cluster_type >= 2) then
            ! Bubble dynamic closure from Maeda and Colonius (2018)

            ! Include the cell that contains the bubble (mapCells+1+mapCells)
            smearGrid = mapCells - (-mapCells) + 1
            smearGridz = smearGrid
            if (p == 0) smearGridz = 1

            charvol = 0._wp
            charpres = 0._wp
            charvol2 = 0._wp
            charpres2 = 0._wp
            vol = 0._wp

            $:GPU_LOOP(parallelism='[seq]')
            do i = 1, smearGrid
                $:GPU_LOOP(parallelism='[seq]')
                do j = 1, smearGrid
                    $:GPU_LOOP(parallelism='[seq]')
                    do k = 1, smearGridz
                        cellaux(1) = cell(1) + i - (mapCells + 1)
                        cellaux(2) = cell(2) + j - (mapCells + 1)
                        cellaux(3) = cell(3) + k - (mapCells + 1)
                        if (p == 0) cellaux(3) = 0

                        !> Obtaining the cell volume
                        if (p > 0) then
                            vol = dx(cellaux(1))*dy(cellaux(2))*dz(cellaux(3))
                        else
                            if (cyl_coord) then
                                vol = dx(cellaux(1))*dy(cellaux(2))*y_cc(cellaux(2))*2._wp*pi
                            else
                                vol = dx(cellaux(1))*dy(cellaux(2))*lag_params%charwidth
                            end if
                        end if

                        !> Update values
                        charvol = charvol + vol
                        charpres = charpres + q_prim_vf(eqn_idx%E)%sf(cellaux(1), cellaux(2), cellaux(3))*vol
                        charvol2 = charvol2 + vol*q_beta(1)%sf(cellaux(1), cellaux(2), cellaux(3))
                        charpres2 = charpres2 + q_prim_vf(eqn_idx%E)%sf(cellaux(1), cellaux(2), &
                                                          & cellaux(3))*vol*q_beta(1)%sf(cellaux(1), cellaux(2), cellaux(3))
                    end do
                end do
            end do
            f_pinfl = charpres2/charvol2
            vol = charvol
            dc = (3._wp*abs(vol)/(4._wp*pi))**(1._wp/3._wp)
        end if

        if (lag_params%pressure_corrector) then
            ! Valid if only one bubble exists per cell
            volgas = intfc_rad(bub_id, 2)**3._wp
            denom = intfc_rad(bub_id, 2)**2._wp
            term1 = bub_dphidt(bub_id)*intfc_rad(bub_id, 2)**2._wp
            term2 = intfc_vel(bub_id, 2)*intfc_rad(bub_id, 2)**2._wp

            Rbeq = volgas**(1._wp/3._wp)  ! surrogate bubble radius
            aux = dc**3._wp - Rbeq**3._wp
            term2 = term2/denom
            term2 = 3._wp/2._wp*term2**2._wp*Rbeq**3._wp*(1._wp - Rbeq/dc)/aux
            preterm1 = 3._wp/2._wp*Rbeq*(dc**2._wp - Rbeq**2._wp)/(aux*denom)

            ! Control volume radius
            if (ptype == 2) Romega = dc

            ! Getting p_inf
            if (ptype == 1) then
                f_pinfl = f_pinfl + preterm1*term1 + term2
            end if
        end if

    end subroutine s_get_pinf

    !> Update Lagrangian bubble variables using TVD Runge-Kutta time stepping
    impure subroutine s_update_lagrange_tdv_rk(q_prim_vf, bc_type, stage)

        type(scalar_field), dimension(sys_size), intent(in)        :: q_prim_vf
        type(integer_field), dimension(1:num_dims,1:2), intent(in) :: bc_type
        integer, intent(in)                                        :: stage
        integer                                                    :: k, q

        if (time_stepper == time_stepper_rk1) then  ! 1st order TVD RK
            $:GPU_PARALLEL_LOOP(private='[k]')
            do k = 1, n_el_bubs_loc
                ! u{1} = u{n} +  dt * RHS{n}
                intfc_rad(k, 1) = intfc_rad(k, 1) + dt*intfc_draddt(k, 1)
                intfc_vel(k, 1) = intfc_vel(k, 1) + dt*intfc_dveldt(k, 1)
                gas_p(k, 1) = gas_p(k, 1) + dt*gas_dpdt(k, 1)
                gas_mv(k, 1) = gas_mv(k, 1) + dt*gas_dmvdt(k, 1)
                if (moving_lag_bubbles) then
                    do q = 1, num_dims
                        mtn_posPrev(k, q, 1) = mtn_pos(k, q, 1)
                        mtn_pos(k, q, 1) = mtn_pos(k, q, 1) + dt*mtn_dposdt(k, q, 1)
                        mtn_vel(k, q, 1) = mtn_vel(k, q, 1) + dt*mtn_dveldt(k, q, 1)
                    end do
                end if
            end do
            $:END_GPU_PARALLEL_LOOP()

            call s_transfer_data_to_tmp()
            if (moving_lag_bubbles) call s_enforce_EL_bubbles_boundary_conditions(q_prim_vf)
            call s_smear_voidfraction(bc_type)
            if (lag_params%write_void_evol) call s_write_void_evol(mytime)
            if (lag_params%write_bubbles_stats) call s_calculate_lag_bubble_stats()
            if (lag_params%write_bubbles) then
                $:GPU_UPDATE(host='[gas_p, gas_mv, intfc_rad, intfc_vel]')
                call s_write_lag_bubble_evol(mytime)
            end if
        else if (time_stepper == time_stepper_rk2) then  ! 2nd order TVD RK
            if (stage == 1) then
                $:GPU_PARALLEL_LOOP(private='[k]')
                do k = 1, n_el_bubs_loc
                    ! u{1} = u{n} +  dt * RHS{n}
                    intfc_rad(k, 2) = intfc_rad(k, 1) + dt*intfc_draddt(k, 1)
                    intfc_vel(k, 2) = intfc_vel(k, 1) + dt*intfc_dveldt(k, 1)
                    gas_p(k, 2) = gas_p(k, 1) + dt*gas_dpdt(k, 1)
                    gas_mv(k, 2) = gas_mv(k, 1) + dt*gas_dmvdt(k, 1)
                    if (moving_lag_bubbles) then
                        do q = 1, num_dims
                            mtn_posPrev(k, q, 2) = mtn_pos(k, q, 1)
                            mtn_pos(k, q, 2) = mtn_pos(k, q, 1) + dt*mtn_dposdt(k, q, 1)
                            mtn_vel(k, q, 2) = mtn_vel(k, q, 1) + dt*mtn_dveldt(k, q, 1)
                        end do
                    end if
                end do
                $:END_GPU_PARALLEL_LOOP()

                if (moving_lag_bubbles) call s_enforce_EL_bubbles_boundary_conditions(q_prim_vf)
                call s_smear_voidfraction(bc_type)
            else if (stage == 2) then
                $:GPU_PARALLEL_LOOP(private='[k]')
                do k = 1, n_el_bubs_loc
                    ! u{1} = u{n} + (1/2) * dt * (RHS{n} + RHS{1})
                    intfc_rad(k, 1) = intfc_rad(k, 1) + dt*(intfc_draddt(k, 1) + intfc_draddt(k, 2))/2._wp
                    intfc_vel(k, 1) = intfc_vel(k, 1) + dt*(intfc_dveldt(k, 1) + intfc_dveldt(k, 2))/2._wp
                    gas_p(k, 1) = gas_p(k, 1) + dt*(gas_dpdt(k, 1) + gas_dpdt(k, 2))/2._wp
                    gas_mv(k, 1) = gas_mv(k, 1) + dt*(gas_dmvdt(k, 1) + gas_dmvdt(k, 2))/2._wp
                    if (moving_lag_bubbles) then
                        do q = 1, num_dims
                            mtn_posPrev(k, q, 1) = mtn_pos(k, q, 2)
                            mtn_pos(k, q, 1) = mtn_pos(k, q, 1) + dt*(mtn_dposdt(k, q, 1) + mtn_dposdt(k, q, 2))/2._wp
                            mtn_vel(k, q, 1) = mtn_vel(k, q, 1) + dt*(mtn_dveldt(k, q, 1) + mtn_dveldt(k, q, 2))/2._wp
                        end do
                    end if
                end do
                $:END_GPU_PARALLEL_LOOP()

                call s_transfer_data_to_tmp()
                if (moving_lag_bubbles) call s_enforce_EL_bubbles_boundary_conditions(q_prim_vf)
                call s_smear_voidfraction(bc_type)
                if (lag_params%write_void_evol) call s_write_void_evol(mytime)
                if (lag_params%write_bubbles_stats) call s_calculate_lag_bubble_stats()
                if (lag_params%write_bubbles) then
                    $:GPU_UPDATE(host='[gas_p, gas_mv, intfc_rad, intfc_vel]')
                    call s_write_lag_bubble_evol(mytime)
                end if
            end if
        else if (time_stepper == time_stepper_rk3) then  ! 3rd order TVD RK
            if (stage == 1) then
                $:GPU_PARALLEL_LOOP(private='[k]')
                do k = 1, n_el_bubs_loc
                    ! u{1} = u{n} +  dt * RHS{n}
                    intfc_rad(k, 2) = intfc_rad(k, 1) + dt*intfc_draddt(k, 1)
                    intfc_vel(k, 2) = intfc_vel(k, 1) + dt*intfc_dveldt(k, 1)
                    gas_p(k, 2) = gas_p(k, 1) + dt*gas_dpdt(k, 1)
                    gas_mv(k, 2) = gas_mv(k, 1) + dt*gas_dmvdt(k, 1)
                    if (moving_lag_bubbles) then
                        do q = 1, num_dims
                            mtn_posPrev(k, q, 2) = mtn_pos(k, q, 1)
                            mtn_pos(k, q, 2) = mtn_pos(k, q, 1) + dt*mtn_dposdt(k, q, 1)
                            mtn_vel(k, q, 2) = mtn_vel(k, q, 1) + dt*mtn_dveldt(k, q, 1)
                        end do
                    end if
                end do
                $:END_GPU_PARALLEL_LOOP()

                if (moving_lag_bubbles) call s_enforce_EL_bubbles_boundary_conditions(q_prim_vf)
                call s_smear_voidfraction(bc_type)
            else if (stage == 2) then
                $:GPU_PARALLEL_LOOP(private='[k]')
                do k = 1, n_el_bubs_loc
                    ! u{2} = u{n} + (1/4) * dt * [RHS{n} + RHS{1}]
                    intfc_rad(k, 2) = intfc_rad(k, 1) + dt*(intfc_draddt(k, 1) + intfc_draddt(k, 2))/4._wp
                    intfc_vel(k, 2) = intfc_vel(k, 1) + dt*(intfc_dveldt(k, 1) + intfc_dveldt(k, 2))/4._wp
                    gas_p(k, 2) = gas_p(k, 1) + dt*(gas_dpdt(k, 1) + gas_dpdt(k, 2))/4._wp
                    gas_mv(k, 2) = gas_mv(k, 1) + dt*(gas_dmvdt(k, 1) + gas_dmvdt(k, 2))/4._wp
                    if (moving_lag_bubbles) then
                        do q = 1, num_dims
                            mtn_posPrev(k, q, 2) = mtn_pos(k, q, 2)
                            mtn_pos(k, q, 2) = mtn_pos(k, q, 1) + dt*(mtn_dposdt(k, q, 1) + mtn_dposdt(k, q, 2))/4._wp
                            mtn_vel(k, q, 2) = mtn_vel(k, q, 1) + dt*(mtn_dveldt(k, q, 1) + mtn_dveldt(k, q, 2))/4._wp
                        end do
                    end if
                end do
                $:END_GPU_PARALLEL_LOOP()

                if (moving_lag_bubbles) call s_enforce_EL_bubbles_boundary_conditions(q_prim_vf)
                call s_smear_voidfraction(bc_type)
            else if (stage == 3) then
                $:GPU_PARALLEL_LOOP(private='[k]')
                do k = 1, n_el_bubs_loc
                    ! u{n+1} = u{n} + (2/3) * dt * [(1/4)* RHS{n} + (1/4)* RHS{1} + RHS{2}]
                    intfc_rad(k, 1) = intfc_rad(k, 1) + (2._wp/3._wp)*dt*(intfc_draddt(k, 1)/4._wp + intfc_draddt(k, &
                              & 2)/4._wp + intfc_draddt(k, 3))
                    intfc_vel(k, 1) = intfc_vel(k, 1) + (2._wp/3._wp)*dt*(intfc_dveldt(k, 1)/4._wp + intfc_dveldt(k, &
                              & 2)/4._wp + intfc_dveldt(k, 3))
                    gas_p(k, 1) = gas_p(k, 1) + (2._wp/3._wp)*dt*(gas_dpdt(k, 1)/4._wp + gas_dpdt(k, 2)/4._wp + gas_dpdt(k, 3))
                    gas_mv(k, 1) = gas_mv(k, 1) + (2._wp/3._wp)*dt*(gas_dmvdt(k, 1)/4._wp + gas_dmvdt(k, 2)/4._wp + gas_dmvdt(k, 3))
                    if (moving_lag_bubbles) then
                        do q = 1, num_dims
                            mtn_posPrev(k, q, 1) = mtn_pos(k, q, 2)
                            mtn_pos(k, q, 1) = mtn_pos(k, q, 1) + (2._wp/3._wp)*dt*(mtn_dposdt(k, q, 1)/4._wp + mtn_dposdt(k, q, &
                                    & 2)/4._wp + mtn_dposdt(k, q, 3))
                            mtn_vel(k, q, 1) = mtn_vel(k, q, 1) + (2._wp/3._wp)*dt*(mtn_dveldt(k, q, 1)/4._wp + mtn_dveldt(k, q, &
                                    & 2)/4._wp + mtn_dveldt(k, q, 3))
                        end do
                    end if
                end do
                $:END_GPU_PARALLEL_LOOP()

                call s_transfer_data_to_tmp()
                if (moving_lag_bubbles) call s_enforce_EL_bubbles_boundary_conditions(q_prim_vf)
                call s_smear_voidfraction(bc_type)
                if (lag_params%write_void_evol) call s_write_void_evol(mytime)
                if (lag_params%write_bubbles_stats) call s_calculate_lag_bubble_stats()
                if (lag_params%write_bubbles) then
                    $:GPU_UPDATE(host='[gas_p, gas_mv, gas_mg, intfc_rad, intfc_vel]')
                    call s_write_lag_bubble_evol(mytime)
                end if
            end if
        end if

    end subroutine s_update_lagrange_tdv_rk

    impure subroutine s_enforce_EL_bubbles_boundary_conditions(q_prim_vf)

        type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
        integer                                             :: k, i, q
        integer                                             :: patch_id, newBubs, new_idx
        real(wp)                                            :: offset
        integer, dimension(3)                               :: cell

        call nvtxStartRange("LAG-BC")
        call nvtxStartRange("LAG-BC-DEV2HOST")
        $:GPU_UPDATE(host='[bub_R0, Rmax_stats, Rmin_stats, gas_mg, gas_betaT, gas_betaC, bub_dphidt, lag_id, gas_p, gas_mv, &
                     & intfc_rad, intfc_vel, mtn_pos, mtn_posPrev, mtn_vel, mtn_s, intfc_draddt, intfc_dveldt, gas_dpdt, &
                     & gas_dmvdt, mtn_dposdt, mtn_dveldt, keep_bubble, n_el_bubs_loc, wrap_bubble_dir, wrap_bubble_loc]')
        call nvtxEndRange

        ! Handle MPI transfer of bubbles going to another processor's local domain
        if (num_procs > 1) then
            call nvtxStartRange("LAG-BC-TRANSFER-LIST")
            call s_add_particles_to_transfer_list(n_el_bubs_loc, mtn_pos(:,:,2), mtn_posPrev(:,:,2))
            call nvtxEndRange

            call nvtxStartRange("LAG-BC-SENDRECV")
            call s_mpi_sendrecv_particles(bub_R0, Rmax_stats, Rmin_stats, gas_mg, gas_betaT, gas_betaC, bub_dphidt, lag_id, &
                                          & gas_p, gas_mv, intfc_rad, intfc_vel, mtn_pos, mtn_posPrev, mtn_vel, mtn_s, &
                                          & intfc_draddt, intfc_dveldt, gas_dpdt, gas_dmvdt, mtn_dposdt, mtn_dveldt, lag_num_ts, &
                                          & n_el_bubs_loc, 2)
            call nvtxEndRange
        end if

        call nvtxStartRange("LAG-BC-HOST2DEV")
        $:GPU_UPDATE(device='[bub_R0, Rmax_stats, Rmin_stats, gas_mg, gas_betaT, gas_betaC, bub_dphidt, lag_id, gas_p, gas_mv, &
                     & intfc_rad, intfc_vel, mtn_pos, mtn_posPrev, mtn_vel, mtn_s, intfc_draddt, intfc_dveldt, gas_dpdt, &
                     & gas_dmvdt, mtn_dposdt, mtn_dveldt, n_el_bubs_loc]')
        call nvtxEndRange

        $:GPU_PARALLEL_LOOP(private='[k, cell]')
        do k = 1, n_el_bubs_loc
            keep_bubble(k) = 1
            wrap_bubble_loc(k,:) = 0
            wrap_bubble_dir(k,:) = 0

            ! Relocate bubbles at solid boundaries and delete bubbles that leave buffer regions
            if (any(bc_x%beg == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/)) .and. mtn_pos(k, 1, &
                & 2) < x_cb(-1) + intfc_rad(k, 2)) then
                mtn_pos(k, 1, 2) = x_cb(-1) + intfc_rad(k, 2)
            else if (any(bc_x%end == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/)) .and. mtn_pos(k, 1, &
                     & 2) > x_cb(m) - intfc_rad(k, 2)) then
                mtn_pos(k, 1, 2) = x_cb(m) - intfc_rad(k, 2)
            else if (bc_x%beg == BC_PERIODIC .and. mtn_pos(k, 1, 2) < pcomm_coords(1)%beg .and. mtn_posPrev(k, 1, &
                     & 2) >= pcomm_coords(1)%beg) then
                wrap_bubble_dir(k, 1) = 1
                wrap_bubble_loc(k, 1) = -1
            else if (bc_x%end == BC_PERIODIC .and. mtn_pos(k, 1, 2) > pcomm_coords(1)%end .and. mtn_posPrev(k, 1, &
                     & 2) <= pcomm_coords(1)%end) then
                wrap_bubble_dir(k, 1) = 1
                wrap_bubble_loc(k, 1) = 1
            else if (.not. (mtn_pos(k, 1, 2) >= x_cb(-1) .and. mtn_pos(k, 1, 2) < x_cb(m))) then
                keep_bubble(k) = 0
            end if

            if (any(bc_y%beg == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/)) .and. mtn_pos(k, 2, &
                & 2) < y_cb(-1) + intfc_rad(k, 2)) then
                mtn_pos(k, 2, 2) = y_cb(-1) + intfc_rad(k, 2)
            else if (any(bc_y%end == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/)) .and. mtn_pos(k, 2, &
                     & 2) > y_cb(n) - intfc_rad(k, 2)) then
                mtn_pos(k, 2, 2) = y_cb(n) - intfc_rad(k, 2)
            else if (bc_y%beg == BC_PERIODIC .and. mtn_pos(k, 2, 2) < pcomm_coords(2)%beg .and. mtn_posPrev(k, 2, &
                     & 2) >= pcomm_coords(2)%beg) then
                wrap_bubble_dir(k, 2) = 1
                wrap_bubble_loc(k, 2) = -1
            else if (bc_y%end == BC_PERIODIC .and. mtn_pos(k, 2, 2) > pcomm_coords(2)%end .and. mtn_posPrev(k, 2, &
                     & 2) <= pcomm_coords(2)%end) then
                wrap_bubble_dir(k, 2) = 1
                wrap_bubble_loc(k, 2) = 1
            else if (.not. (mtn_pos(k, 2, 2) >= y_cb(-1) .and. mtn_pos(k, 2, 2) < y_cb(n))) then
                keep_bubble(k) = 0
            end if

            if (p > 0) then
                if (any(bc_z%beg == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/)) .and. mtn_pos(k, 3, &
                    & 2) < z_cb(-1) + intfc_rad(k, 2)) then
                    mtn_pos(k, 3, 2) = z_cb(-1) + intfc_rad(k, 2)
                else if (any(bc_z%end == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/)) .and. mtn_pos(k, 3, &
                         & 2) > z_cb(p) - intfc_rad(k, 2)) then
                    mtn_pos(k, 3, 2) = z_cb(p) - intfc_rad(k, 2)
                else if (bc_z%beg == BC_PERIODIC .and. mtn_pos(k, 3, 2) < pcomm_coords(3)%beg .and. mtn_posPrev(k, 3, &
                         & 2) >= pcomm_coords(3)%beg) then
                    wrap_bubble_dir(k, 3) = 1
                    wrap_bubble_loc(k, 3) = -1
                else if (bc_z%end == BC_PERIODIC .and. mtn_pos(k, 3, 2) > pcomm_coords(3)%end .and. mtn_posPrev(k, 3, &
                         & 2) <= pcomm_coords(3)%end) then
                    wrap_bubble_dir(k, 3) = 1
                    wrap_bubble_loc(k, 3) = 1
                else if (.not. (mtn_pos(k, 3, 2) >= z_cb(-1) .and. mtn_pos(k, 3, 2) < z_cb(p))) then
                    keep_bubble(k) = 0
                end if
            end if

            if (keep_bubble(k) == 1) then
                ! Remove bubbles that are no longer in a liquid
                cell = fd_number - buff_size
                call s_locate_cell(mtn_pos(k,1:3,2), cell, mtn_s(k,1:3,2))

                if (q_prim_vf(eqn_idx%adv%beg)%sf(cell(1), cell(2), cell(3)) < (1._wp - lag_params%valmaxvoid)) then
                    keep_bubble(k) = 0
                end if
            end if
        end do
        $:END_GPU_PARALLEL_LOOP()

        if (n_el_bubs_loc > 0) then
            call nvtxStartRange("LAG-BC-DEV2HOST")
            $:GPU_UPDATE(host='[bub_R0, Rmax_stats, Rmin_stats, gas_mg, gas_betaT, gas_betaC, bub_dphidt, lag_id, gas_p, gas_mv, &
                         & intfc_rad, intfc_vel, mtn_pos, mtn_posPrev, mtn_vel, mtn_s, intfc_draddt, intfc_dveldt, gas_dpdt, &
                         & gas_dmvdt, mtn_dposdt, mtn_dveldt, keep_bubble, n_el_bubs_loc, wrap_bubble_dir, wrap_bubble_loc]')
            call nvtxEndRange

            newBubs = 0
            do k = 1, n_el_bubs_loc
                if (keep_bubble(k) == 1) then
                    newBubs = newBubs + 1
                    if (newBubs /= k) then
                        call s_copy_lag_bubble(newBubs, k)
                        wrap_bubble_dir(newBubs,:) = wrap_bubble_dir(k,:)
                        wrap_bubble_loc(newBubs,:) = wrap_bubble_loc(k,:)
                    end if
                end if
            end do
            n_el_bubs_loc = newBubs

            ! Handle periodic wrapping of bubbles on same processor
            do k = 1, n_el_bubs_loc
                if (any(wrap_bubble_dir(k,:) == 1)) then
                    do i = 1, num_dims
                        if (wrap_bubble_dir(k, i) == 1) then
                            offset = glb_bounds(i)%end - glb_bounds(i)%beg
                            if (wrap_bubble_loc(k, i) == 1) then
                                do q = 1, 2
                                    mtn_pos(k, i, q) = mtn_pos(k, i, q) - offset
                                    mtn_posPrev(k, i, q) = mtn_posPrev(k, i, q) - offset
                                end do
                            else if (wrap_bubble_loc(k, i) == -1) then
                                do q = 1, 2
                                    mtn_pos(k, i, q) = mtn_pos(k, i, q) + offset
                                    mtn_posPrev(k, i, q) = mtn_posPrev(k, i, q) + offset
                                end do
                            end if
                        end if
                    end do
                end if
            end do
            call nvtxStartRange("LAG-BC-HOST2DEV")
            $:GPU_UPDATE(device='[bub_R0, Rmax_stats, Rmin_stats, gas_mg, gas_betaT, gas_betaC, bub_dphidt, lag_id, gas_p, &
                         & gas_mv, intfc_rad, intfc_vel, mtn_pos, mtn_posPrev, mtn_vel, mtn_s, intfc_draddt, intfc_dveldt, &
                         & gas_dpdt, gas_dmvdt, mtn_dposdt, mtn_dveldt, n_el_bubs_loc]')
            call nvtxEndRange
        end if

        $:GPU_PARALLEL_LOOP(private='[cell]')
        do k = 1, n_el_bubs_loc
            cell = fd_number - buff_size
            call s_locate_cell(mtn_pos(k,1:3,2), cell, mtn_s(k,1:3,2))
        end do
        $:END_GPU_PARALLEL_LOOP()

        call nvtxEndRange  ! LAG-BC

    end subroutine s_enforce_EL_bubbles_boundary_conditions

    !> Locate the cell index for a given physical position
    subroutine s_locate_cell(pos, cell, scoord)

        $:GPU_ROUTINE(function_name='s_locate_cell',parallelism='[seq]', cray_inline=True)

        real(wp), dimension(3), intent(in)   :: pos
        real(wp), dimension(3), intent(out)  :: scoord
        integer, dimension(3), intent(inout) :: cell
        integer                              :: i

        do while (pos(1) < x_cb(cell(1) - 1))
            cell(1) = cell(1) - 1
        end do

        do while (pos(1) >= x_cb(cell(1)))
            cell(1) = cell(1) + 1
        end do

        do while (pos(2) < y_cb(cell(2) - 1))
            cell(2) = cell(2) - 1
        end do

        do while (pos(2) >= y_cb(cell(2)))
            cell(2) = cell(2) + 1
        end do

        if (p > 0) then
            do while (pos(3) < z_cb(cell(3) - 1))
                cell(3) = cell(3) - 1
            end do
            do while (pos(3) >= z_cb(cell(3)))
                cell(3) = cell(3) + 1
            end do
        else
            cell(3) = 0
        end if

        ! The numbering of the cell of which left boundary is the domain boundary is 0. if comp.coord of the pos is s, the real
        ! coordinate of s is (the coordinate of the left boundary of the Floor(s)-th cell) + (s-(int(s))*(cell-width). In other
        ! words, the coordinate of the center of the cell is x_cc(cell).

        ! coordinates in computational space
        scoord(1) = cell(1) + (pos(1) - x_cb(cell(1) - 1))/dx(cell(1))
        scoord(2) = cell(2) + (pos(2) - y_cb(cell(2) - 1))/dy(cell(2))
        scoord(3) = 0._wp
        if (p > 0) scoord(3) = cell(3) + (pos(3) - z_cb(cell(3) - 1))/dz(cell(3))

    end subroutine s_locate_cell

    !> Transfer data into the temporal variables
    impure subroutine s_transfer_data_to_tmp()

        integer :: k

        $:GPU_PARALLEL_LOOP(private='[k]')
        do k = 1, n_el_bubs_loc
            gas_p(k, 2) = gas_p(k, 1)
            gas_mv(k, 2) = gas_mv(k, 1)
            intfc_rad(k, 2) = intfc_rad(k, 1)
            intfc_vel(k, 2) = intfc_vel(k, 1)
            mtn_pos(k,1:3,2) = mtn_pos(k,1:3,1)
            mtn_posPrev(k,1:3,2) = mtn_posPrev(k,1:3,1)
            mtn_vel(k,1:3,2) = mtn_vel(k,1:3,1)
            mtn_s(k,1:3,2) = mtn_s(k,1:3,1)
        end do
        $:END_GPU_PARALLEL_LOOP()

    end subroutine s_transfer_data_to_tmp

    !> Determine if a bubble position lies within the current MPI subdomain including ghost cells
    function particle_in_domain(pos_part)

        logical                            :: particle_in_domain
        real(wp), dimension(3), intent(in) :: pos_part

        ! 2D

        if (p == 0 .and. cyl_coord .neqv. .true.) then
            ! Defining a virtual z-axis that has the same dimensions as y-axis defined in the input file
            particle_in_domain = ((pos_part(1) < x_cb(m + buff_size - fd_number)) .and. (pos_part(1) >= x_cb(fd_number &
                                  & - buff_size - 1)) .and. (pos_part(2) < y_cb(n + buff_size - fd_number)) .and. (pos_part(2) &
                                  & >= y_cb(fd_number - buff_size - 1)) .and. (pos_part(3) < lag_params%charwidth/2._wp) &
                                  & .and. (pos_part(3) > -lag_params%charwidth/2._wp))
        else
            ! cyl_coord
            particle_in_domain = ((pos_part(1) < x_cb(m + buff_size - fd_number)) .and. (pos_part(1) >= x_cb(fd_number &
                                  & - buff_size - 1)) .and. (abs(pos_part(2)) < y_cb(n + buff_size - fd_number)) &
                                  & .and. (abs(pos_part(2)) >= max(y_cb(fd_number - buff_size - 1), 0._wp)))
        end if

        ! 3D
        if (p > 0) then
            particle_in_domain = ((pos_part(1) < x_cb(m + buff_size - fd_number)) .and. (pos_part(1) >= x_cb(fd_number &
                                  & - buff_size - 1)) .and. (pos_part(2) < y_cb(n + buff_size - fd_number)) .and. (pos_part(2) &
                                  & >= y_cb(fd_number - buff_size - 1)) .and. (pos_part(3) < z_cb(p + buff_size - fd_number)) &
                                  & .and. (pos_part(3) >= z_cb(fd_number - buff_size - 1)))
        end if

        ! For symmetric and wall boundary condition
        if (any(bc_x%beg == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/))) then
            particle_in_domain = (particle_in_domain .and. (pos_part(1) >= x_cb(-1)))
        end if
        if (any(bc_x%end == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/))) then
            particle_in_domain = (particle_in_domain .and. (pos_part(1) < x_cb(m)))
        end if
        if (any(bc_y%beg == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/)) .and. (.not. cyl_coord)) then
            particle_in_domain = (particle_in_domain .and. (pos_part(2) >= y_cb(-1)))
        end if
        if (any(bc_y%end == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/)) .and. (.not. cyl_coord)) then
            particle_in_domain = (particle_in_domain .and. (pos_part(2) < y_cb(n)))
        end if
        if (p > 0) then
            if (any(bc_z%beg == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/))) then
                particle_in_domain = (particle_in_domain .and. (pos_part(3) >= z_cb(-1)))
            end if
            if (any(bc_z%end == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/))) then
                particle_in_domain = (particle_in_domain .and. (pos_part(3) < z_cb(p)))
            end if
        end if

    end function particle_in_domain

    !> Determine if a Lagrangian bubble is within the physical domain excluding ghost cells
    function particle_in_domain_physical(pos_part)

        logical                            :: particle_in_domain_physical
        real(wp), dimension(3), intent(in) :: pos_part

        particle_in_domain_physical = ((pos_part(1) < x_cb(m)) .and. (pos_part(1) >= x_cb(-1)) .and. (pos_part(2) < y_cb(n)) &
                                       & .and. (pos_part(2) >= y_cb(-1)))

        if (p > 0) then
            particle_in_domain_physical = (particle_in_domain_physical .and. (pos_part(3) < z_cb(p)) .and. (pos_part(3) &
                                           & >= z_cb(-1)))
        end if

    end function particle_in_domain_physical

    !> Compute the gradient of a scalar field using second-order central differences on a non-uniform grid
    subroutine s_gradient_dir(q, dq, dir)

        real(stp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:), intent(inout) :: q, dq
        integer, intent(in)                                                                  :: dir
        integer                                                                              :: i, j, k

        if (dir == 1) then
            ! Gradient in x dir.
            $:GPU_PARALLEL_LOOP(private='[i, j, k]', collapse=3)
            do k = 0, p
                do j = 0, n
                    do i = 0, m
                        dq(i, j, k) = q(i, j, k)*(dx(i + 1) - dx(i - 1)) + q(i + 1, j, k)*(dx(i) + dx(i - 1)) - q(i - 1, j, &
                           & k)*(dx(i) + dx(i + 1))
                        dq(i, j, k) = dq(i, j, k)/((dx(i) + dx(i - 1))*(dx(i) + dx(i + 1)))
                    end do
                end do
            end do
            $:END_GPU_PARALLEL_LOOP()
        else if (dir == 2) then
            ! Gradient in y dir.
            $:GPU_PARALLEL_LOOP(private='[i, j, k]', collapse=3)
            do k = 0, p
                do j = 0, n
                    do i = 0, m
                        dq(i, j, k) = q(i, j, k)*(dy(j + 1) - dy(j - 1)) + q(i, j + 1, k)*(dy(j) + dy(j - 1)) - q(i, j - 1, &
                           & k)*(dy(j) + dy(j + 1))
                        dq(i, j, k) = dq(i, j, k)/((dy(j) + dy(j - 1))*(dy(j) + dy(j + 1)))
                    end do
                end do
            end do
            $:END_GPU_PARALLEL_LOOP()
        else if (dir == 3) then
            ! Gradient in z dir.
            $:GPU_PARALLEL_LOOP(private='[i, j, k]', collapse=3)
            do k = 0, p
                do j = 0, n
                    do i = 0, m
                        dq(i, j, k) = q(i, j, k)*(dz(k + 1) - dz(k - 1)) + q(i, j, k + 1)*(dz(k) + dz(k - 1)) - q(i, j, &
                           & k - 1)*(dz(k) + dz(k + 1))
                        dq(i, j, k) = dq(i, j, k)/((dz(k) + dz(k - 1))*(dz(k) + dz(k + 1)))
                    end do
                end do
            end do
            $:END_GPU_PARALLEL_LOOP()
        end if

    end subroutine s_gradient_dir

    !> Write Lagrangian bubble state data at each time step
    impure subroutine s_write_lag_particles(qtime)

        real(wp), intent(in)                 :: qtime
        integer                              :: k
        character(LEN=path_len + 2*name_len) :: file_loc
        logical                              :: file_exist
        character(LEN=25)                    :: FMT

        write (file_loc, '(A,I0,A)') 'lag_bubble_evol_', proc_rank, '.dat'
        file_loc = trim(case_dir) // '/D/' // trim(file_loc)
        call my_inquire(trim(file_loc), file_exist)

        if (precision == precision_single) then
            FMT = "(A16,A14,8A16)"
        else
            FMT = "(A24,A14,8A24)"
        end if

        if (.not. file_exist) then
            open (LAG_EVOL_ID, FILE=trim(file_loc), form='formatted', position='rewind')
            write (LAG_EVOL_ID, FMT) 'currentTime', 'particleID', 'x', 'y', 'z', 'coreVaporMass', 'coreVaporConcentration', &
                   & 'radius', 'interfaceVelocity', 'corePressure'
        else
            open (LAG_EVOL_ID, FILE=trim(file_loc), form='formatted', position='append')
        end if

    end subroutine s_write_lag_particles

    !> Open the file to write the evolution of the lagrangian bubbles on each time step.
    impure subroutine s_open_lag_bubble_evol()

        character(LEN=path_len + 2*name_len) :: file_loc
        logical                              :: file_exist
        character(LEN=25)                    :: FMT

        write (file_loc, '(A,I0,A)') 'lag_bubble_evol_', proc_rank, '.dat'
        file_loc = trim(case_dir) // '/D/' // trim(file_loc)
        call my_inquire(trim(file_loc), file_exist)

        if (precision == precision_single) then
            FMT = "(A16,A14,8A16)"
        else
            FMT = "(A24,A14,8A24)"
        end if

        if (.not. file_exist) then
            open (LAG_EVOL_ID, FILE=trim(file_loc), form='formatted', position='rewind')
            write (LAG_EVOL_ID, FMT) 'currentTime', 'particleID', 'x', 'y', 'z', 'coreVaporMass', 'coreVaporConcentration', &
                   & 'radius', 'interfaceVelocity', 'corePressure'
        else
            open (LAG_EVOL_ID, FILE=trim(file_loc), form='formatted', position='append')
        end if

    end subroutine s_open_lag_bubble_evol

    !> Write on each time step the changes of the lagrangian bubbles.
    impure subroutine s_write_lag_bubble_evol(qtime)

        real(wp), intent(in)                 :: qtime
        integer                              :: k, ios
        character(LEN=25)                    :: FMT
        character(LEN=path_len + 2*name_len) :: file_loc, path
        logical                              :: file_exist

        if (precision == precision_single) then
            FMT = "(F16.8,I14,8F16.8)"
        else
            FMT = "(F24.16,I14,8F24.16)"
        end if

        ! Cycle through list
        do k = 1, n_el_bubs_loc
            write (LAG_EVOL_ID, FMT) qtime, lag_id(k, 1), mtn_pos(k, 1, 1), mtn_pos(k, 2, 1), mtn_pos(k, 3, 1), gas_mv(k, 1), &
                   & gas_mv(k, 1)/(gas_mv(k, 1) + gas_mg(k)), intfc_rad(k, 1), intfc_vel(k, 1), gas_p(k, 1)
        end do

    end subroutine s_write_lag_bubble_evol

    impure subroutine s_close_lag_bubble_evol

        close (LAG_EVOL_ID)

    end subroutine s_close_lag_bubble_evol

    subroutine s_open_void_evol

        character(LEN=path_len + 2*name_len) :: file_loc
        logical                              :: file_exist

        if (proc_rank == 0) then
            write (file_loc, '(A)') 'voidfraction.dat'
            file_loc = trim(case_dir) // '/D/' // trim(file_loc)
            call my_inquire(trim(file_loc), file_exist)
            if (.not. file_exist) then
                open (LAG_VOID_ID, FILE=trim(file_loc), form='formatted', position='rewind')
            else
                open (LAG_VOID_ID, FILE=trim(file_loc), form='formatted', position='append')
            end if
        end if

    end subroutine s_open_void_evol

    !> Write some useful statistics related to the volume fraction of the particles (void fraction)
    impure subroutine s_write_void_evol(qtime)

        real(wp), intent(in)                 :: qtime
        real(wp)                             :: volcell, voltot
        real(wp)                             :: lag_void_max, lag_void_avg, lag_vol
        real(wp)                             :: void_max_glb, void_avg_glb, vol_glb
        integer                              :: i, j, k
        character(LEN=path_len + 2*name_len) :: file_loc
        logical                              :: file_exist

        lag_void_max = 0._wp
        lag_void_avg = 0._wp
        lag_vol = 0._wp
        $:GPU_PARALLEL_LOOP(private='[volcell]', collapse=3, reduction='[[lag_vol, lag_void_avg], [lag_void_max]]', &
                            & reductionOp='[+, MAX]', copy='[lag_vol, lag_void_avg, lag_void_max]')
        do k = 0, p
            do j = 0, n
                do i = 0, m
                    lag_void_max = max(lag_void_max, 1._wp - q_beta(1)%sf(i, j, k))
                    call s_get_char_vol(i, j, k, volcell)
                    if ((1._wp - q_beta(1)%sf(i, j, k)) > 5.0e-11_wp) then
                        lag_void_avg = lag_void_avg + (1._wp - q_beta(1)%sf(i, j, k))*volcell
                        lag_vol = lag_vol + volcell
                    end if
                end do
            end do
        end do
        $:END_GPU_PARALLEL_LOOP()

#ifdef MFC_MPI
        if (num_procs > 1) then
            call s_mpi_allreduce_max(lag_void_max, void_max_glb)
            lag_void_max = void_max_glb
            call s_mpi_allreduce_sum(lag_vol, vol_glb)
            lag_vol = vol_glb
            call s_mpi_allreduce_sum(lag_void_avg, void_avg_glb)
            lag_void_avg = void_avg_glb
        end if
#endif
        voltot = lag_void_avg
        ! This voidavg value does not reflect the real void fraction in the cloud since the cell which does not have bubbles are not
        ! accounted
        if (lag_vol > 0._wp) lag_void_avg = lag_void_avg/lag_vol

        if (proc_rank == 0) then
            write (LAG_VOID_ID, '(6X,4e24.8)') qtime, lag_void_avg, lag_void_max, voltot
        end if

    end subroutine s_write_void_evol

    subroutine s_close_void_evol

        if (proc_rank == 0) close (LAG_VOID_ID)

    end subroutine s_close_void_evol

    !> Write restart files for the Lagrangian bubble solver
    impure subroutine s_write_restart_lag_bubbles(t_step)

        ! Generic string used to store the address of a particular file
        integer, intent(in)                  :: t_step
        character(LEN=path_len + 2*name_len) :: file_loc
        logical                              :: file_exist
        integer                              :: bub_id, tot_part
        integer                              :: i, k

#ifdef MFC_MPI
        ! For Parallel I/O
        integer                                :: ifile, ierr
        integer, dimension(MPI_STATUS_SIZE)    :: status
        integer(KIND=MPI_OFFSET_KIND)          :: disp
        integer                                :: view
        integer, dimension(2)                  :: gsizes, lsizes, start_idx_part
        integer, allocatable                   :: proc_bubble_counts(:)
        real(wp), dimension(1:1,1:lag_io_vars) :: dummy
        dummy = 0._wp

        bub_id = 0
        if (n_el_bubs_loc /= 0) then
            do k = 1, n_el_bubs_loc
                if (particle_in_domain_physical(mtn_pos(k,1:3,1))) then
                    bub_id = bub_id + 1
                end if
            end do
        end if

        if (.not. parallel_io) return

        allocate (proc_bubble_counts(num_procs))

        lsizes(1) = bub_id
        lsizes(2) = lag_io_vars

        ! Total number of particles
        call MPI_ALLREDUCE(bub_id, tot_part, 1, MPI_integer, MPI_SUM, MPI_COMM_WORLD, ierr)

        call MPI_ALLGATHER(bub_id, 1, MPI_INTEGER, proc_bubble_counts, 1, MPI_INTEGER, MPI_COMM_WORLD, ierr)

        ! Calculate starting index for this processor's particles
        call MPI_EXSCAN(lsizes(1), start_idx_part(1), 1, MPI_INTEGER, MPI_SUM, MPI_COMM_WORLD, ierr)
        if (proc_rank == 0) start_idx_part(1) = 0
        start_idx_part(2) = 0

        gsizes(1) = tot_part
        gsizes(2) = lag_io_vars

        write (file_loc, '(A,I0,A)') 'lag_bubbles_', t_step, '.dat'
        file_loc = trim(case_dir) // '/restart_data' // trim(mpiiofs) // trim(file_loc)

        ! Clean up existing file
        if (proc_rank == 0) then
            inquire (FILE=trim(file_loc), EXIST=file_exist)
            if (file_exist) then
                call MPI_FILE_DELETE(file_loc, mpi_info_int, ierr)
            end if
        end if

        call MPI_BARRIER(MPI_COMM_WORLD, ierr)

        if (proc_rank == 0) then
            call MPI_FILE_OPEN(MPI_COMM_SELF, file_loc, ior(MPI_MODE_WRONLY, MPI_MODE_CREATE), mpi_info_int, ifile, ierr)

            ! Write header using MPI I/O for consistency
            call MPI_FILE_WRITE(ifile, tot_part, 1, MPI_INTEGER, status, ierr)
            call MPI_FILE_WRITE(ifile, mytime, 1, mpi_p, status, ierr)
            call MPI_FILE_WRITE(ifile, dt, 1, mpi_p, status, ierr)
            call MPI_FILE_WRITE(ifile, num_procs, 1, MPI_INTEGER, status, ierr)
            call MPI_FILE_WRITE(ifile, proc_bubble_counts, num_procs, MPI_INTEGER, status, ierr)

            call MPI_FILE_CLOSE(ifile, ierr)
        end if

        call MPI_BARRIER(MPI_COMM_WORLD, ierr)

        if (bub_id > 0) then
            allocate (MPI_IO_DATA_lag_bubbles(max(1, bub_id),1:lag_io_vars))

            i = 0
            do k = 1, n_el_bubs_loc
                if (.not. particle_in_domain_physical(mtn_pos(k,1:3,1))) cycle
                i = i + 1
                MPI_IO_DATA_lag_bubbles(i, 1) = real(lag_id(k, 1))
                MPI_IO_DATA_lag_bubbles(i,2:4) = mtn_pos(k,1:3,1)
                MPI_IO_DATA_lag_bubbles(i,5:7) = mtn_posPrev(k,1:3,1)
                MPI_IO_DATA_lag_bubbles(i,8:10) = mtn_vel(k,1:3,1)
                MPI_IO_DATA_lag_bubbles(i, 11) = intfc_rad(k, 1)
                MPI_IO_DATA_lag_bubbles(i, 12) = intfc_vel(k, 1)
                MPI_IO_DATA_lag_bubbles(i, 13) = bub_R0(k)
                MPI_IO_DATA_lag_bubbles(i, 14) = Rmax_stats(k)
                MPI_IO_DATA_lag_bubbles(i, 15) = Rmin_stats(k)
                MPI_IO_DATA_lag_bubbles(i, 16) = bub_dphidt(k)
                MPI_IO_DATA_lag_bubbles(i, 17) = gas_p(k, 1)
                MPI_IO_DATA_lag_bubbles(i, 18) = gas_mv(k, 1)
                MPI_IO_DATA_lag_bubbles(i, 19) = gas_mg(k)
                MPI_IO_DATA_lag_bubbles(i, 20) = gas_betaT(k)
                MPI_IO_DATA_lag_bubbles(i, 21) = gas_betaC(k)
            end do

            call MPI_TYPE_CREATE_SUBARRAY(2, gsizes, lsizes, start_idx_part, MPI_ORDER_FORTRAN, mpi_p, view, ierr)
            call MPI_TYPE_COMMIT(view, ierr)

            call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, ior(MPI_MODE_WRONLY, MPI_MODE_CREATE), mpi_info_int, ifile, ierr)

            ! Skip header (written by rank 0)
            disp = int(sizeof(tot_part) + 2*sizeof(mytime) + sizeof(num_procs) + num_procs*sizeof(proc_bubble_counts(1)), &
                       & MPI_OFFSET_KIND)
            call MPI_FILE_SET_VIEW(ifile, disp, mpi_p, view, 'native', mpi_info_int, ierr)

            call MPI_FILE_WRITE_ALL(ifile, MPI_IO_DATA_lag_bubbles, lag_io_vars*bub_id, mpi_p, status, ierr)

            call MPI_FILE_CLOSE(ifile, ierr)

            deallocate (MPI_IO_DATA_lag_bubbles)
        else
            call MPI_TYPE_CONTIGUOUS(0, mpi_p, view, ierr)
            call MPI_TYPE_COMMIT(view, ierr)

            call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, ior(MPI_MODE_WRONLY, MPI_MODE_CREATE), mpi_info_int, ifile, ierr)

            ! Skip header (written by rank 0)
            disp = int(sizeof(tot_part) + 2*sizeof(mytime) + sizeof(num_procs) + num_procs*sizeof(proc_bubble_counts(1)), &
                       & MPI_OFFSET_KIND)
            call MPI_FILE_SET_VIEW(ifile, disp, mpi_p, view, 'native', mpi_info_int, ierr)

            call MPI_FILE_WRITE_ALL(ifile, dummy, 0, mpi_p, status, ierr)

            call MPI_FILE_CLOSE(ifile, ierr)
        end if

        deallocate (proc_bubble_counts)
#endif

    end subroutine s_write_restart_lag_bubbles

    !> Compute the maximum and minimum radius of each bubble
    subroutine s_calculate_lag_bubble_stats()

        integer :: k

        $:GPU_PARALLEL_LOOP(private='[k]', reduction='[[Rmax_glb], [Rmin_glb]]', reductionOp='[MAX, MIN]', &
                            & copy='[Rmax_glb, Rmin_glb]')
        do k = 1, n_el_bubs_loc
            Rmax_glb = max(Rmax_glb, intfc_rad(k, 1)/bub_R0(k))
            Rmin_glb = min(Rmin_glb, intfc_rad(k, 1)/bub_R0(k))
            Rmax_stats(k) = max(Rmax_stats(k), intfc_rad(k, 1)/bub_R0(k))
            Rmin_stats(k) = min(Rmin_stats(k), intfc_rad(k, 1)/bub_R0(k))
        end do
        $:END_GPU_PARALLEL_LOOP()

    end subroutine s_calculate_lag_bubble_stats

    impure subroutine s_open_lag_bubble_stats()

        character(LEN=path_len + 2*name_len) :: file_loc
        character(LEN=20)                    :: FMT
        logical                              :: file_exist

        write (file_loc, '(A,I0,A)') 'stats_lag_bubbles_', proc_rank, '.dat'
        file_loc = trim(case_dir) // '/D/' // trim(file_loc)
        call my_inquire(trim(file_loc), file_exist)

        if (precision == precision_single) then
            FMT = "(A10,A14,5A16)"
        else
            FMT = "(A10,A14,5A24)"
        end if

        if (.not. file_exist) then
            open (LAG_STATS_ID, FILE=trim(file_loc), form='formatted', position='rewind')
            write (LAG_STATS_ID, *) 'proc_rank, particleID, x, y, z, Rmax_glb, Rmin_glb'
        else
            open (LAG_STATS_ID, FILE=trim(file_loc), form='formatted', position='append')
        end if

    end subroutine s_open_lag_bubble_stats

    !> Write the maximum and minimum radius of each bubble.
    impure subroutine s_write_lag_bubble_stats()

        integer                              :: k
        character(LEN=path_len + 2*name_len) :: file_loc
        character(LEN=20)                    :: FMT

        $:GPU_UPDATE(host='[Rmax_glb, Rmin_glb]')

        if (precision == precision_single) then
            FMT = "(I10,I14,5F16.8)"
        else
            FMT = "(I10,I14,5F24.16)"
        end if

        do k = 1, n_el_bubs_loc
            write (LAG_STATS_ID, FMT) proc_rank, lag_id(k, 1), mtn_pos(k, 1, 1), mtn_pos(k, 2, 1), mtn_pos(k, 3, 1), &
                   & Rmax_stats(k), Rmin_stats(k)
        end do

    end subroutine s_write_lag_bubble_stats

    subroutine s_close_lag_bubble_stats

        close (LAG_STATS_ID)

    end subroutine s_close_lag_bubble_stats

    !> Copy a lag bubble from a src index to a destination index
    impure subroutine s_copy_lag_bubble(dest, src)

        integer, intent(in) :: src, dest

        bub_R0(dest) = bub_R0(src)
        Rmax_stats(dest) = Rmax_stats(src)
        Rmin_stats(dest) = Rmin_stats(src)
        gas_mg(dest) = gas_mg(src)
        gas_betaT(dest) = gas_betaT(src)
        gas_betaC(dest) = gas_betaC(src)
        bub_dphidt(dest) = bub_dphidt(src)
        lag_id(dest, 1) = lag_id(src, 1)
        gas_p(dest,1:2) = gas_p(src,1:2)
        gas_mv(dest,1:2) = gas_mv(src,1:2)
        intfc_rad(dest,1:2) = intfc_rad(src,1:2)
        intfc_vel(dest,1:2) = intfc_vel(src,1:2)
        mtn_vel(dest,1:3,1:2) = mtn_vel(src,1:3,1:2)
        mtn_s(dest,1:3,1:2) = mtn_s(src,1:3,1:2)
        mtn_pos(dest,1:3,1:2) = mtn_pos(src,1:3,1:2)
        mtn_posPrev(dest,1:3,1:2) = mtn_posPrev(src,1:3,1:2)
        intfc_draddt(dest,1:lag_num_ts) = intfc_draddt(src,1:lag_num_ts)
        intfc_dveldt(dest,1:lag_num_ts) = intfc_dveldt(src,1:lag_num_ts)
        gas_dpdt(dest,1:lag_num_ts) = gas_dpdt(src,1:lag_num_ts)
        gas_dmvdt(dest,1:lag_num_ts) = gas_dmvdt(src,1:lag_num_ts)
        mtn_dposdt(dest,1:3,1:lag_num_ts) = mtn_dposdt(src,1:3,1:lag_num_ts)
        mtn_dveldt(dest,1:3,1:lag_num_ts) = mtn_dveldt(src,1:3,1:lag_num_ts)

    end subroutine s_copy_lag_bubble

    impure subroutine s_finalize_lagrangian_solver()

        integer :: i

        if (lag_params%write_void_evol) call s_close_void_evol
        if (lag_params%write_bubbles) call s_close_lag_bubble_evol()
        if (lag_params%write_bubbles_stats) call s_close_lag_bubble_stats()

        do i = 1, q_beta_idx
            @:DEALLOCATE(q_beta(i)%sf)
            if (lag_params%kahan_summation) then
                @:DEALLOCATE(kahan_comp(i)%sf)
            end if
        end do
        @:DEALLOCATE(q_beta)
        if (lag_params%kahan_summation) then
            @:DEALLOCATE(kahan_comp)
        end if

        ! Deallocating space
        @:DEALLOCATE(lag_id)
        @:DEALLOCATE(bub_R0)
        @:DEALLOCATE(Rmax_stats)
        @:DEALLOCATE(Rmin_stats)
        @:DEALLOCATE(gas_mg)
        @:DEALLOCATE(gas_betaT)
        @:DEALLOCATE(gas_betaC)
        @:DEALLOCATE(bub_dphidt)
        @:DEALLOCATE(gas_p)
        @:DEALLOCATE(gas_mv)
        @:DEALLOCATE(intfc_rad)
        @:DEALLOCATE(intfc_vel)
        @:DEALLOCATE(mtn_pos)
        @:DEALLOCATE(mtn_posPrev)
        @:DEALLOCATE(mtn_vel)
        @:DEALLOCATE(mtn_s)
        @:DEALLOCATE(intfc_draddt)
        @:DEALLOCATE(intfc_dveldt)
        @:DEALLOCATE(gas_dpdt)
        @:DEALLOCATE(gas_dmvdt)
        @:DEALLOCATE(mtn_dposdt)
        @:DEALLOCATE(mtn_dveldt)

        @:DEALLOCATE(keep_bubble)
        @:DEALLOCATE(wrap_bubble_loc, wrap_bubble_dir)

        if (lag_params%vel_model > 0 .and. lag_params%pressure_force) then
            @:DEALLOCATE(grad_p_x)
            @:DEALLOCATE(fd_coeff_x_pgrad)
            if (n > 0) then
                @:DEALLOCATE(grad_p_y)
                @:DEALLOCATE(fd_coeff_y_pgrad)
                if (p > 0) then
                    @:DEALLOCATE(grad_p_z)
                    @:DEALLOCATE(fd_coeff_z_pgrad)
                end if
            end if
        end if

        @:DEALLOCATE(cell_list_start)
        @:DEALLOCATE(cell_list_count)
        @:DEALLOCATE(cell_list_idx)

    end subroutine s_finalize_lagrangian_solver

end module m_bubbles_EL