m_bubbles.fpp.f90

Go to the documentation of this file.

# 1 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!>
!! @file
!! @brief Contains module m_bubbles
 
# 1 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 1
# 1 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 1
# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
! New line at end of file is required for FYPP
# 2 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
# 1 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 1
# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
! New line at end of file is required for FYPP
# 2 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 2
 
# 4 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 5 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 6 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 7 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 8 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 20 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 43 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 48 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 53 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 58 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 63 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 68 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 76 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 81 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 86 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 91 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 96 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 101 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 106 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 111 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 116 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 121 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 151 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 192 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 207 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 232 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 243 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 245 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 255 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 283 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 293 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 303 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 312 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 329 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 339 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 346 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 352 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 358 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 364 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 370 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 376 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
! New line at end of file is required for FYPP
# 3 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
# 1 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 1
# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
! New line at end of file is required for FYPP
# 2 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 2
 
# 7 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 17 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 22 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 27 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 32 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 37 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 42 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 47 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 52 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 57 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 62 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 73 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 78 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 83 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 88 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 103 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 131 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 160 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 175 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 192 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 213 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 241 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 256 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 266 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 275 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 291 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 301 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 308 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
! New line at end of file is required for FYPP
# 4 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
 
# 21 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 37 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 50 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 76 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 91 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 102 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 115 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 143 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 154 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 165 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 176 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 187 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 198 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 208 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 214 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 220 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 226 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 232 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 234 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
# 235 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
! New line at end of file is required for FYPP
# 2 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 2
 
# 14 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
! Caution:
! This macro requires the use of a binding script to set CUDA_VISIBLE_DEVICES, such that we have one GPU device per MPI rank.
! That's because for both cudaMemAdvise (preferred location) and cudaMemPrefetchAsync we use location = device_id = 0.
! For an example see misc/nvidia_uvm/bind.sh.
# 63 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 81 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 88 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 111 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 127 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 153 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 159 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 167 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
! New line at end of file is required for FYPP
# 6 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp" 2
 
!> @brief Shared bubble-dynamics procedures (radial acceleration, wall pressure, sound speed) for ensemble- and volume-averaged models
module m_bubbles
 
    use m_derived_types        !< definitions of the derived types
 
    use m_global_parameters    !< definitions of the global parameters
 
    use m_mpi_proxy            !< message passing interface (mpi) module proxy
 
    use m_variables_conversion !< state variables type conversion procedures
 
    use m_helper_basic         !< functions to compare floating point numbers
 
    implicit none
 
    real(wp) :: chi_vw  !< Bubble wall properties (Ando 2010)
    real(wp) :: k_mw    !< Bubble wall properties (Ando 2010)
    real(wp) :: rho_mw  !< Bubble wall properties (Ando 2010)
 
# 25 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#if defined(MFC_OpenACC)
# 25 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$acc declare create(chi_vw, k_mw, rho_mw) 
# 25 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#elif defined(MFC_OpenMP)
# 25 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$omp declare target (chi_vw, k_mw, rho_mw) 
# 25 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#endif
 
contains
 
    !> Function that computes the bubble radial acceleration based on bubble models
        !!  @param fRho Current density
        !!  @param fP Current driving pressure
        !!  @param fR Current bubble radius
        !!  @param fV Current bubble velocity
        !!  @param fR0 Equilibrium bubble radius
        !!  @param fpb Internal bubble pressure
        !!  @param fpbdot Time-derivative of internal bubble pressure
        !!  @param alf bubble volume fraction
        !!  @param fntait Tait EOS parameter
        !!  @param fBtait Tait EOS parameter
        !!  @param f_bub_adv_src Source for bubble volume fraction
        !!  @param f_divu Divergence of velocity
        !!  @param fCson Speed of sound from fP (EL)
    elemental function f_rddot(fRho, fP, fR, fV, fR0, fpb, fpbdot, alf, fntait, fBtait, f_bub_adv_src, f_divu, fCson)
 
# 44 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#if MFC_OpenACC
# 44 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$acc routine seq 
# 44 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#elif MFC_OpenMP
# 44 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
    
# 44 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
 
# 44 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$omp declare target device_type(any) 
# 44 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#endif
        real(wp), intent(in) :: frho, fp, fr, fv, fr0, fpb, fpbdot, alf
        real(wp), intent(in) :: fntait, fbtait, f_bub_adv_src, f_divu
        real(wp), intent(in) :: fcson
 
        real(wp) :: fcpbw, fcpinf, fcpinf_dot, fh, fhdot, c_gas, c_liquid
        real(wp) :: f_rddot
 
        if (bubble_model == 1) then
            ! Gilmore bubbles
            fcpinf = fp - eu
            fcpbw = f_cpbw(fr0, fr, fv, fpb)
            fh = f_h(fcpbw, fcpinf, fntait, fbtait)
            c_gas = f_cgas(fcpinf, fntait, fbtait, fh)
            fcpinf_dot = f_cpinfdot(frho, fp, alf, fntait, fbtait, f_bub_adv_src, f_divu)
            fhdot = f_hdot(fcpbw, fcpinf, fcpinf_dot, fntait, fbtait, fr, fv, fr0, fpbdot)
            f_rddot = f_rddot_g(fcpbw, fr, fv, fh, fhdot, c_gas, fntait, fbtait)
        else if (bubble_model == 2) then
            ! Keller-Miksis bubbles
            fcpinf = fp
            fcpbw = f_cpbw_km(fr0, fr, fv, fpb)
            if (bubbles_euler) then
                c_liquid = sqrt(fntait*(fp + fbtait)/(frho*(1._wp - alf)))
            else
                c_liquid = fcson
            end if
            f_rddot = f_rddot_km(fpbdot, fcpinf, fcpbw, frho, fr, fv, fr0, c_liquid)
        else if (bubble_model == 3) then
            ! Rayleigh-Plesset bubbles
            fcpbw = f_cpbw_km(fr0, fr, fv, fpb)
            f_rddot = f_rddot_rp(fp, frho, fr, fv, fcpbw)
        else
            ! Default: No bubble dynamics
            f_rddot = 0._wp
        end if
 
    end function f_rddot
 
    !>  Function that computes that bubble wall pressure for Gilmore bubbles
        !!  @param fR0 Equilibrium bubble radius
        !!  @param fR Current bubble radius
        !!  @param fV Current bubble velocity
        !!  @param fpb Internal bubble pressure
    elemental function f_cpbw(fR0, fR, fV, fpb)
 
# 88 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#if MFC_OpenACC
# 88 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$acc routine seq 
# 88 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#elif MFC_OpenMP
# 88 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
    
# 88 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
 
# 88 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$omp declare target device_type(any) 
# 88 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#endif
        real(wp), intent(in) :: fr0, fr, fv, fpb
 
        real(wp) :: f_cpbw
 
        if (polytropic) then
            f_cpbw = (ca + 2._wp/web/fr0)*((fr0/fr)**(3._wp*gam)) - ca - 4._wp*re_inv*fv/fr - 2._wp/(fr*web)
        else
            f_cpbw = fpb - 1._wp - 4._wp*re_inv*fv/fr - 2._wp/(fr*web)
        end if
 
    end function f_cpbw
 
    !>  Function that computes the bubble enthalpy
        !!  @param fCpbw Bubble wall pressure
        !!  @param fCpinf Driving bubble pressure
        !!  @param fntait Tait EOS parameter
        !!  @param fBtait Tait EOS parameter
    elemental function f_h(fCpbw, fCpinf, fntait, fBtait)
 
# 107 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#if MFC_OpenACC
# 107 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$acc routine seq 
# 107 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#elif MFC_OpenMP
# 107 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
    
# 107 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
 
# 107 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$omp declare target device_type(any) 
# 107 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#endif
        real(wp), intent(in) :: fcpbw, fcpinf, fntait, fbtait
 
        real(wp) :: tmp1, tmp2, tmp3
        real(wp) :: f_h
 
        tmp1 = (fntait - 1._wp)/fntait
        tmp2 = (fcpbw/(1._wp + fbtait) + 1._wp)**tmp1
        tmp3 = (fcpinf/(1._wp + fbtait) + 1._wp)**tmp1
 
        f_h = (tmp2 - tmp3)*fntait*(1._wp + fbtait)/(fntait - 1._wp)
 
    end function f_h
 
    !> Function that computes the sound speed for the bubble
        !! @param fCpinf Driving bubble pressure
        !! @param fntait Tait EOS parameter
        !! @param fBtait Tait EOS parameter
        !! @param fH Bubble enthalpy
    elemental function f_cgas(fCpinf, fntait, fBtait, fH)
 
# 127 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#if MFC_OpenACC
# 127 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$acc routine seq 
# 127 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#elif MFC_OpenMP
# 127 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
    
# 127 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
 
# 127 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$omp declare target device_type(any) 
# 127 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#endif
        real(wp), intent(in) :: fcpinf, fntait, fbtait, fh
 
        real(wp) :: tmp
        real(wp) :: f_cgas
 
        ! get sound speed for Gilmore equations "C" -> c_gas
        tmp = (fcpinf/(1._wp + fbtait) + 1._wp)**((fntait - 1._wp)/fntait)
        tmp = fntait*(1._wp + fbtait)*tmp
 
        f_cgas = sqrt(tmp + (fntait - 1._wp)*fh)
 
    end function f_cgas
 
    !>  Function that computes the time derivative of the driving pressure
        !!  @param fRho Local liquid density
        !!  @param fP Local pressure
        !!  @param falf Local void fraction
        !!  @param fntait Tait EOS parameter
        !!  @param fBtait Tait EOS parameter
        !!  @param advsrc Advection equation source term
        !!  @param divu Divergence of velocity
    elemental function f_cpinfdot(fRho, fP, falf, fntait, fBtait, advsrc, divu)
 
# 150 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#if MFC_OpenACC
# 150 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$acc routine seq 
# 150 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#elif MFC_OpenMP
# 150 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
    
# 150 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
 
# 150 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$omp declare target device_type(any) 
# 150 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#endif
        real(wp), intent(in) :: frho, fp, falf, fntait, fbtait, advsrc, divu
 
        real(wp) :: c2_liquid
        real(wp) :: f_cpinfdot
 
        ! get sound speed squared for liquid (only needed for pbdot)
        ! c_l^2 = gam (p+B) / (rho*(1-alf))
        if (mpp_lim) then
            c2_liquid = fntait*(fp + fbtait)/frho
        else
            c2_liquid = fntait*(fp + fbtait)/(frho*(1._wp - falf))
        end if
 
        ! \dot{Cp_inf} = rho sound^2 (alf_src - divu)
        f_cpinfdot = frho*c2_liquid*(advsrc - divu)
 
    end function f_cpinfdot
 
    !>  Function that computes the time derivative of the enthalpy
        !!  @param fCpbw Bubble wall pressure
        !!  @param fCpinf Driving bubble pressure
        !!  @param fCpinf_dot Time derivative of the driving pressure
        !!  @param fntait Tait EOS parameter
        !!  @param fBtait Tait EOS parameter
        !!  @param fR Current bubble radius
        !!  @param fV Current bubble velocity
        !!  @param fR0 Equilibrium bubble radius
        !!  @param fpbdot Time derivative of the internal bubble pressure
    elemental function f_hdot(fCpbw, fCpinf, fCpinf_dot, fntait, fBtait, fR, fV, fR0, fpbdot)
 
# 180 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#if MFC_OpenACC
# 180 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$acc routine seq 
# 180 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#elif MFC_OpenMP
# 180 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
    
# 180 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
 
# 180 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$omp declare target device_type(any) 
# 180 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#endif
        real(wp), intent(in) :: fcpbw, fcpinf, fcpinf_dot, fntait, fbtait
        real(wp), intent(in) :: fr, fv, fr0, fpbdot
 
        real(wp) :: tmp1, tmp2
        real(wp) :: f_hdot
 
        if (polytropic) then
            tmp1 = (fr0/fr)**(3._wp*gam)
            tmp1 = -3._wp*gam*(ca + 2._wp/web/fr0)*tmp1*fv/fr
        else
            tmp1 = fpbdot
        end if
        tmp2 = (2._wp/web + 4._wp*re_inv*fv)*fv/(fr**2._wp)
 
        f_hdot = &
            (fcpbw/(1._wp + fbtait) + 1._wp)**(-1._wp/fntait)*(tmp1 + tmp2) &
            - (fcpinf/(1._wp + fbtait) + 1._wp)**(-1._wp/fntait)*fcpinf_dot
 
        ! Hdot = (Cpbw/(1+B) + 1)^(-1/n_tait)*(-3 gam)*(R0/R)^(3gam) V/R
        !f_Hdot = ((fCpbw/(1._wp+fBtait)+1._wp)**(-1._wp/fntait))*(-3._wp)*gam * &
        !            ( (fR0/fR)**(3._wp*gam ))*(fV/fR)
 
        ! Hdot = Hdot - (Cpinf/(1+B) + 1)^(-1/n_tait) Cpinfdot
        !f_Hdot = f_Hdot - ((fCpinf/(1._wp+fBtait)+1._wp)**(-1._wp/fntait))*fCpinf_dot
 
    end function f_hdot
 
    !>  Function that computes the bubble radial acceleration for Rayleigh-Plesset bubbles
        !!  @param fCp Driving pressure
        !!  @param fRho Current density
        !!  @param fR Current bubble radius
        !!  @param fV Current bubble velocity
        !!  @param fCpbw Boundary wall pressure
    elemental function f_rddot_rp(fCp, fRho, fR, fV, fCpbw)
 
# 215 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#if MFC_OpenACC
# 215 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$acc routine seq 
# 215 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#elif MFC_OpenMP
# 215 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
    
# 215 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
 
# 215 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$omp declare target device_type(any) 
# 215 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#endif
        real(wp), intent(in) :: fcp, frho, fr, fv, fcpbw
 
        real(wp) :: f_rddot_rp
 
            !! rddot = (1/r) (  -3/2 rdot^2 + ((r0/r)^3\gamma - Cp)/rho )
            !! rddot = (1/r) (  -3/2 rdot^2 + (tmp1 - Cp)/rho )
            !! rddot = (1/r) (  tmp2 )
 
        f_rddot_rp = (-1.5_wp*(fv**2._wp) + (fcpbw - fcp)/frho)/fr
 
    end function f_rddot_rp
 
    !>  Function that computes the bubble radial acceleration
        !!  @param fCpbw Bubble wall pressure
        !!  @param fR Current bubble radius
        !!  @param fV Current bubble velocity
        !!  @param fH Current enthalpy
        !!  @param fHdot Current time derivative of the enthalpy
        !!  @param fcgas Current gas sound speed
        !!  @param fntait Tait EOS parameter
        !!  @param fBtait Tait EOS parameter
    elemental function f_rddot_g(fCpbw, fR, fV, fH, fHdot, fcgas, fntait, fBtait)
 
# 238 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#if MFC_OpenACC
# 238 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$acc routine seq 
# 238 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#elif MFC_OpenMP
# 238 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
    
# 238 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
 
# 238 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$omp declare target device_type(any) 
# 238 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#endif
        real(wp), intent(in) :: fcpbw, fr, fv, fh, fhdot
        real(wp), intent(in) :: fcgas, fntait, fbtait
 
        real(wp) :: tmp1, tmp2, tmp3
        real(wp) :: f_rddot_g
 
        tmp1 = fv/fcgas
        tmp2 = 1._wp + 4._wp*re_inv/fcgas/fr*(fcpbw/(1._wp + fbtait) + 1._wp) &
               **(-1._wp/fntait)
        tmp3 = 1.5_wp*fv**2._wp*(tmp1/3._wp - 1._wp) + fh*(1._wp + tmp1) &
               + fr*fhdot*(1._wp - tmp1)/fcgas
 
        f_rddot_g = tmp3/(fr*(1._wp - tmp1)*tmp2)
 
    end function f_rddot_g
 
    !>  Function that computes the bubble wall pressure for Keller--Miksis bubbles
        !!  @param fR0 Equilibrium bubble radius
        !!  @param fR Current bubble radius
        !!  @param fV Current bubble velocity
        !!  @param fpb Internal bubble pressure
    elemental function f_cpbw_km(fR0, fR, fV, fpb)
 
# 261 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#if MFC_OpenACC
# 261 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$acc routine seq 
# 261 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#elif MFC_OpenMP
# 261 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
    
# 261 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
 
# 261 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$omp declare target device_type(any) 
# 261 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#endif
        real(wp), intent(in) :: fr0, fr, fv, fpb
        real(wp) :: f_cpbw_km
 
        if (polytropic) then
            f_cpbw_km = ca*((fr0/fr)**(3._wp*gam)) - ca + eu
            if (.not. f_is_default(web)) f_cpbw_km = f_cpbw_km + &
                                                     (2._wp/(web*fr0))*((fr0/fr)**(3._wp*gam))
        else
            f_cpbw_km = fpb
        end if
 
        if (.not. f_is_default(web)) f_cpbw_km = f_cpbw_km - 2._wp/(fr*web)
        if (.not. f_is_default(re_inv)) f_cpbw_km = f_cpbw_km - 4._wp*re_inv*fv/fr
 
    end function f_cpbw_km
 
    !>  Function that computes the bubble radial acceleration for Keller--Miksis bubbles
        !!  @param fpbdot Time-derivative of internal bubble pressure
        !!  @param fCp Driving pressure
        !!  @param fCpbw Bubble wall pressure
        !!  @param fRho Current density
        !!  @param fR Current bubble radius
        !!  @param fV Current bubble velocity
        !!  @param fR0 Equilibrium bubble radius
        !!  @param fC Current sound speed
    elemental function f_rddot_km(fpbdot, fCp, fCpbw, fRho, fR, fV, fR0, fC)
 
# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#if MFC_OpenACC
# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$acc routine seq 
# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#elif MFC_OpenMP
# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
    
# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
 
# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$omp declare target device_type(any) 
# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#endif
        real(wp), intent(in) :: fpbdot, fcp, fcpbw
        real(wp), intent(in) :: frho, fr, fv, fr0, fc
 
        real(wp) :: tmp1, tmp2, cdot_star
        real(wp) :: f_rddot_km
        if (polytropic) then
            cdot_star = -3._wp*gam*ca*((fr0/fr)**(3._wp*gam))*fv/fr
            if (.not. f_is_default(web)) cdot_star = cdot_star - &
                                                     3._wp*gam*(2._wp/(web*fr0))*((fr0/fr)**(3._wp*gam))*fv/fr
        else
            cdot_star = fpbdot
        end if
 
        if (.not. f_is_default(web)) cdot_star = cdot_star + (2._wp/web)*fv/(fr**2._wp)
        if (.not. f_is_default(re_inv)) cdot_star = cdot_star + 4._wp*re_inv*((fv/fr)**2._wp)
 
        tmp1 = fv/fc
        tmp2 = 1.5_wp*(fv**2._wp)*(tmp1/3._wp - 1._wp) + &
               (1._wp + tmp1)*(fcpbw - fcp)/frho + &
               cdot_star*fr/(frho*fc)
 
        if (f_is_default(re_inv)) then
            f_rddot_km = tmp2/(fr*(1._wp - tmp1))
        else
            f_rddot_km = tmp2/(fr*(1._wp - tmp1) + 4._wp*re_inv/(frho*fc))
        end if
 
    end function f_rddot_km
 
    !>  Subroutine that computes bubble wall properties for vapor bubbles
        !!  @param pb_in Internal bubble pressure
        !!  @param iR0 Current bubble size index
    elemental subroutine s_bwproperty(pb_in, iR0, chi_vw_out, k_mw_out, rho_mw_out)
 
# 322 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#if MFC_OpenACC
# 322 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$acc routine seq 
# 322 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#elif MFC_OpenMP
# 322 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
    
# 322 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
 
# 322 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$omp declare target device_type(any) 
# 322 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#endif
        real(wp), intent(in) :: pb_in
        integer, intent(in) :: ir0
        real(wp), intent(out) :: chi_vw_out  !< Bubble wall properties (Ando 2010)
        real(wp), intent(out) :: k_mw_out    !< Bubble wall properties (Ando 2010)
        real(wp), intent(out) :: rho_mw_out  !< Bubble wall properties (Ando 2010)
        real(wp) :: x_vw
 
        ! mass fraction of vapor
        chi_vw_out = 1._wp/(1._wp + r_v/r_g*(pb_in/pv - 1._wp))
        ! mole fraction of vapor & thermal conductivity of gas mixture
        x_vw = m_g*chi_vw_out/(m_v + (m_g - m_v)*chi_vw_out)
        k_mw_out = x_vw*k_v(ir0)/(x_vw + (1._wp - x_vw)*phi_vg) &
                   + (1._wp - x_vw)*k_g(ir0)/(x_vw*phi_gv + 1._wp - x_vw)
        ! gas mixture density
        rho_mw_out = pv/(chi_vw_out*r_v*tw)
 
    end subroutine s_bwproperty
 
    !>  Function that computes the vapour flux
        !!  @param fR Current bubble radius
        !!  @param fV Current bubble velocity
        !!  @param fpb
        !!  @param fmass_v Current mass of vapour
        !!  @param iR0 Bubble size index (EE) or bubble identifier (EL)
        !!  @param vflux Computed vapour flux
        !!  @param fmass_g Current gas mass (EL)
        !!  @param fbeta_c Mass transfer coefficient (EL)
        !!  @param fR_m Mixture gas constant (EL)
        !!  @param fgamma_m Mixture gamma (EL)
    elemental subroutine s_vflux(fR, fV, fpb, fmass_v, iR0, vflux, fmass_g, fbeta_c, fR_m, fgamma_m)
 
# 353 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#if MFC_OpenACC
# 353 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$acc routine seq 
# 353 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#elif MFC_OpenMP
# 353 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
    
# 353 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
 
# 353 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$omp declare target device_type(any) 
# 353 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#endif
        real(wp), intent(in) :: fr
        real(wp), intent(in) :: fv
        real(wp), intent(in) :: fpb
        real(wp), intent(in) :: fmass_v
        integer, intent(in) :: ir0
        real(wp), intent(out) :: vflux
        real(wp), intent(in), optional :: fmass_g, fbeta_c
        real(wp), intent(out), optional :: fr_m, fgamma_m
 
        real(wp) :: chi_bar
        real(wp) :: rho_mw_lag
        real(wp) :: grad_chi
        real(wp) :: conc_v
 
        if (thermal == 3) then !transfer
            ! constant transfer model
            if (bubbles_lagrange) then
                ! Mixture properties (gas+vapor) in the bubble
                conc_v = fmass_v/(fmass_v + fmass_g)
                if (lag_params%massTransfer_model) then
                    conc_v = 1._wp/(1._wp + (r_v/r_g)*(fpb/pv - 1._wp))
                end if
                fr_m = (fmass_g*r_g + fmass_v*r_v)
                fgamma_m = conc_v*gam_v + (1._wp - conc_v)*gam_g
 
                ! Vapor flux
                chi_bar = fmass_v/(fmass_v + fmass_g)
                grad_chi = (chi_bar - conc_v)
                rho_mw_lag = (fmass_g + fmass_v)/(4._wp/3._wp*pi*fr**3._wp)
                vflux = 0._wp
                if (lag_params%massTransfer_model) then
                    vflux = -fbeta_c*rho_mw_lag*grad_chi/(1._wp - conc_v)/fr
                end if
            else
                chi_bar = fmass_v/(fmass_v + mass_g0(ir0))
                grad_chi = -re_trans_c(ir0)*(chi_bar - chi_vw)
                vflux = rho_mw*grad_chi/pe_c/(1._wp - chi_vw)/fr
            end if
        else
            ! polytropic
            vflux = pv*fv/(r_v*tw)
        end if
 
    end subroutine s_vflux
 
    !>  Function that computes the time derivative of
        !!  the internal bubble pressure
        !!  @param fvflux Vapour flux
        !!  @param fR Current bubble radius
        !!  @param fV Current bubble velocity
        !!  @param fpb Current internal bubble pressure
        !!  @param fmass_v Current mass of vapour
        !!  @param iR0 Bubble size index (EE) or bubble identifier (EL)
        !!  @param fbeta_t Mass transfer coefficient (EL)
        !!  @param fR_m Mixture gas constant (EL)
        !!  @param fgamma_m Mixture gamma (EL)
    elemental function f_bpres_dot(fvflux, fR, fV, fpb, fmass_v, iR0, fbeta_t, fR_m, fgamma_m)
 
# 411 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#if MFC_OpenACC
# 411 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$acc routine seq 
# 411 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#elif MFC_OpenMP
# 411 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
    
# 411 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
 
# 411 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$omp declare target device_type(any) 
# 411 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#endif
        real(wp), intent(in) :: fvflux
        real(wp), intent(in) :: fr
        real(wp), intent(in) :: fv
        real(wp), intent(in) :: fpb
        real(wp), intent(in) :: fmass_v
        integer, intent(in) :: ir0
        real(wp), intent(in), optional :: fbeta_t, fr_m, fgamma_m
 
        real(wp) :: t_bar
        real(wp) :: grad_t
        real(wp) :: f_bpres_dot
        real(wp) :: heatflux
 
        if (thermal == 3) then
            if (bubbles_lagrange) then
                t_bar = fpb*(4._wp/3._wp*pi*fr**3._wp)/fr_m
                grad_t = -fbeta_t*(t_bar - tw)
                heatflux = (fgamma_m - 1._wp)/fgamma_m*grad_t/fr
                f_bpres_dot = 3._wp*fgamma_m*(-fv*fpb + fvflux*r_v*tw &
                                              + heatflux)/fr
                return
            end if
            grad_t = -re_trans_t(ir0)*((fpb/pb0(ir0))*(fr/r0(ir0))**3 &
                                       *(mass_g0(ir0) + mass_v0(ir0))/(mass_g0(ir0) + fmass_v) - 1._wp)
            f_bpres_dot = 3._wp*gam_m*(-fv*fpb + fvflux*r_v*tw &
                                       + pb0(ir0)*k_mw*grad_t/pe_t(ir0)/fr)/fr
        else
            f_bpres_dot = -3._wp*gam_m*fv/fr*(fpb - pv)
        end if
 
    end function f_bpres_dot
 
    !> Adaptive time stepping routine for subgrid bubbles
        !!  (See Heirer, E. Hairer S.P.Nørsett G. Wanner, Solving Ordinary
        !!  Differential Equations I, Chapter II.4)
        !!  @param fRho Current density
        !!  @param fP Current driving pressure
        !!  @param fR Current bubble radius
        !!  @param fV Current bubble radial velocity
        !!  @param fR0 Equilibrium bubble radius
        !!  @param fpb Internal bubble pressure
        !!  @param fpbdot Time-derivative of internal bubble pressure
        !!  @param alf bubble volume fraction
        !!  @param fntait Tait EOS parameter
        !!  @param fBtait Tait EOS parameter
        !!  @param f_bub_adv_src Source for bubble volume fraction
        !!  @param f_divu Divergence of velocity
        !!  @param bub_id Bubble identifier (EL)
        !!  @param fmass_v Current mass of vapour (EL)
        !!  @param fmass_g Current mass of gas (EL)
        !!  @param fbeta_c Mass transfer coefficient (EL)
        !!  @param fbeta_t Heat transfer coefficient (EL)
        !!  @param fCson Speed of sound (EL)
        !!  @param adap_dt_stop Fail-safe exit if max iteration count reached
    subroutine s_advance_step(fRho, fP, fR, fV, fR0, fpb, fpbdot, alf, &
                              fntait, fBtait, f_bub_adv_src, f_divu, &
                              bub_id, fmass_v, fmass_g, fbeta_c, &
                              fbeta_t, fCson, adap_dt_stop)
 
# 470 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#ifdef _CRAYFTN
# 470 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!DIR$ INLINEALWAYS s_advance_step
# 470 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#elif MFC_OpenACC
# 470 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$acc routine seq 
# 470 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#elif MFC_OpenMP
# 470 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
    
# 470 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
 
# 470 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$omp declare target device_type(any) 
# 470 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#endif
# 472 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
 
        real(wp), intent(inout) :: fR, fV, fpb, fmass_v
        real(wp), intent(in) :: fRho, fP, fR0, fpbdot, alf
        real(wp), intent(in) :: fntait, fBtait, f_bub_adv_src, f_divu
        integer, intent(in) :: bub_id
        real(wp), intent(in) :: fmass_g, fbeta_c, fbeta_t, fCson
        integer, intent(inout) :: adap_dt_stop
 
        real(wp), dimension(5) :: err !< Error estimates for adaptive time stepping
        real(wp) :: t_new !< Updated time step size
        real(wp) :: h0, h !< Time step size
        real(wp), dimension(4) :: myR_tmp1, myV_tmp1, myR_tmp2, myV_tmp2 !< Bubble radius, radial velocity, and radial acceleration for the inner loop
        real(wp), dimension(4) :: myPb_tmp1, myMv_tmp1, myPb_tmp2, myMv_tmp2 !< Gas pressure and vapor mass for the inner loop (EL)
        real(wp) :: fR2, fV2, fpb2, fmass_v2
        integer :: iter_count
 
        call s_initial_substep_h(frho, fp, fr, fv, fr0, fpb, fpbdot, alf, &
                                 fntait, fbtait, f_bub_adv_src, f_divu, fcson, h0)
        h = h0
        ! Advancing one step
        t_new = 0._wp
        iter_count = 0
        adap_dt_stop = 0
 
        do
            if (t_new + h > 0.5_wp*dt) then
                h = 0.5_wp*dt - t_new
            end if
 
            ! Advancing one sub-step
            do while (iter_count < adap_dt_max_iters)
 
                iter_count = iter_count + 1
 
                ! Advance one sub-step
                call s_advance_substep(err(1), &
                                       frho, fp, fr, fv, fr0, fpb, fpbdot, alf, &
                                       fntait, fbtait, f_bub_adv_src, f_divu, &
                                       bub_id, fmass_v, fmass_g, fbeta_c, &
                                       fbeta_t, fcson, h, &
                                       myr_tmp1, myv_tmp1, mypb_tmp1, mymv_tmp1)
                if (err(1) > adap_dt_tol) then
                    h = 0.25_wp*h
                    cycle
                end if
 
                ! Advance one sub-step by advancing two half steps
                call s_advance_substep(err(2), &
                                       frho, fp, fr, fv, fr0, fpb, fpbdot, alf, &
                                       fntait, fbtait, f_bub_adv_src, f_divu, &
                                       bub_id, fmass_v, fmass_g, fbeta_c, &
                                       fbeta_t, fcson, 0.5_wp*h, &
                                       myr_tmp2, myv_tmp2, mypb_tmp2, mymv_tmp2)
                if (err(2) > adap_dt_tol) then
                    h = 0.25_wp*h
                    cycle
                end if
 
                fr2 = myr_tmp2(4); fv2 = myv_tmp2(4)
                fpb2 = mypb_tmp2(4); fmass_v2 = mymv_tmp2(4)
 
                call s_advance_substep(err(3), &
                                       frho, fp, fr2, fv2, fr0, fpb2, fpbdot, alf, &
                                       fntait, fbtait, f_bub_adv_src, f_divu, &
                                       bub_id, fmass_v2, fmass_g, fbeta_c, &
                                       fbeta_t, fcson, 0.5_wp*h, &
                                       myr_tmp2, myv_tmp2, mypb_tmp2, mymv_tmp2)
                if (err(3) > adap_dt_tol) then
                    h = 0.5_wp*h
                    cycle
                end if
 
                err(4) = abs((myr_tmp1(4) - myr_tmp2(4))/myr_tmp1(4))
                err(5) = abs((myv_tmp1(4) - myv_tmp2(4))/myv_tmp1(4))
                if (abs(myv_tmp1(4)) < verysmall) err(5) = 0._wp
 
                ! Determine acceptance/rejection and update step size
                !   Rule 1: err1, err2, err3 < tol
                !   Rule 2: myR_tmp1(4) > 0._wp
                !   Rule 3: abs((myR_tmp1(4) - myR_tmp2(4))/fR) < tol
                !   Rule 4: abs((myV_tmp1(4) - myV_tmp2(4))/fV) < tol
                if ((err(1) <= adap_dt_tol) .and. (err(2) <= adap_dt_tol) .and. &
                    (err(3) <= adap_dt_tol) .and. (err(4) <= adap_dt_tol) .and. &
                    (err(5) <= adap_dt_tol) .and. myr_tmp1(4) > 0._wp) then
 
                    ! Accepted. Finalize the sub-step
                    t_new = t_new + h
 
                    ! Update R and V
                    fr = myr_tmp1(4)
                    fv = myv_tmp1(4)
 
                    if (bubbles_lagrange) then
                        ! Update pb and mass_v
                        fpb = mypb_tmp1(4)
                        fmass_v = mymv_tmp1(4)
                    end if
 
                    ! Update step size for the next sub-step
                    h = h*min(2._wp, max(0.5_wp, (adap_dt_tol/err(1))**(1._wp/3._wp)))
 
                    exit
                else
                    ! Rejected. Update step size for the next try on sub-step
                    if (err(2) <= adap_dt_tol) then
                        h = 0.5_wp*h
                    else
                        h = 0.25_wp*h
                    end if
                end if
            end do
 
            ! Exit the loop if the final time reached dt
            if (f_approx_equal(t_new, 0.5_wp*dt) .or. iter_count >= adap_dt_max_iters) exit
 
        end do
 
        if (iter_count >= adap_dt_max_iters) adap_dt_stop = 1
 
    end subroutine s_advance_step
 
    !> Choose the initial time step size for the adaptive time stepping routine
        !!  (See Heirer, E. Hairer S.P.Nørsett G. Wanner, Solving Ordinary
        !!  Differential Equations I, Chapter II.4)
        !!  @param fRho Current density
        !!  @param fP Current driving pressure
        !!  @param fR Current bubble radius
        !!  @param fV Current bubble velocity
        !!  @param fR0 Equilibrium bubble radius
        !!  @param fpb Internal bubble pressure
        !!  @param fpbdot Time-derivative of internal bubble pressure
        !!  @param alf bubble volume fraction
        !!  @param fntait Tait EOS parameter
        !!  @param fBtait Tait EOS parameter
        !!  @param f_bub_adv_src Source for bubble volume fraction
        !!  @param f_divu Divergence of velocity
        !!  @param fCson Speed of sound (EL)
        !!  @param h Time step size
    subroutine s_initial_substep_h(fRho, fP, fR, fV, fR0, fpb, fpbdot, alf, &
                                   fntait, fBtait, f_bub_adv_src, f_divu, &
                                   fCson, h)
 
# 613 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#ifdef _CRAYFTN
# 613 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!DIR$ INLINEALWAYS s_initial_substep_h
# 613 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#elif MFC_OpenACC
# 613 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$acc routine seq 
# 613 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#elif MFC_OpenMP
# 613 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
    
# 613 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
 
# 613 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$omp declare target device_type(any) 
# 613 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#endif
# 615 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
 
        real(wp), intent(IN) :: fRho, fP, fR, fV, fR0, fpb, fpbdot, alf
        real(wp), intent(IN) :: fntait, fBtait, f_bub_adv_src, f_divu
        real(wp), intent(IN) :: fCson
        real(wp), intent(OUT) :: h
 
        real(wp), dimension(2) :: h_size !< Time step size (h0, h1)
        real(wp), dimension(3) :: d_norms !< norms (d_0, d_1, d_2)
        real(wp), dimension(2) :: myR_tmp, myV_tmp, myA_tmp !< Bubble radius, radial velocity, and radial acceleration
 
        ! Determine the starting time step
        ! Evaluate f(x0,y0)
        myr_tmp(1) = fr
        myv_tmp(1) = fv
        mya_tmp(1) = f_rddot(frho, fp, myr_tmp(1), myv_tmp(1), fr0, &
                             fpb, fpbdot, alf, fntait, fbtait, &
                             f_bub_adv_src, f_divu, &
                             fcson)
 
        ! Compute d_0 = ||y0|| and d_1 = ||f(x0,y0)||
        d_norms(1) = sqrt((myr_tmp(1)**2._wp + myv_tmp(1)**2._wp)/2._wp)
        d_norms(2) = sqrt((myv_tmp(1)**2._wp + mya_tmp(1)**2._wp)/2._wp)
        if (d_norms(1) < threshold_first_guess .or. d_norms(2) < threshold_first_guess) then
            h_size(1) = small_guess
        else
            h_size(1) = scale_guess*(d_norms(1)/d_norms(2))
        end if
 
        ! Evaluate f(x0+h0,y0+h0*f(x0,y0))
        myr_tmp(2) = myr_tmp(1) + h_size(1)*myv_tmp(1)
        myv_tmp(2) = myv_tmp(1) + h_size(1)*mya_tmp(1)
        mya_tmp(2) = f_rddot(frho, fp, myr_tmp(2), myv_tmp(2), fr0, &
                             fpb, fpbdot, alf, fntait, fbtait, &
                             f_bub_adv_src, f_divu, &
                             fcson)
 
        ! Compute d_2 = ||f(x0+h0,y0+h0*f(x0,y0))-f(x0,y0)||/h0
        d_norms(3) = sqrt(((myv_tmp(2) - myv_tmp(1))**2._wp + (mya_tmp(2) - mya_tmp(1))**2._wp)/2._wp)/h_size(1)
 
        ! Set h1 = (0.01/max(d_1,d_2))^{1/(p+1)}
        !      if max(d_1,d_2) < 1.e-15_wp, h_size(2) = max(1.e-6_wp, h0*1.e-3_wp)
        if (max(d_norms(2), d_norms(3)) < threshold_second_guess) then
            h_size(2) = max(small_guess, h_size(1)*scale_first_guess)
        else
            h_size(2) = (scale_guess/max(d_norms(2), d_norms(3)))**(1._wp/3._wp)
        end if
 
        h = min(h_size(1)/scale_guess, h_size(2))
 
    end subroutine s_initial_substep_h
 
    !>  Integrate bubble variables over the given time step size, h, using a
        !!      third-order accurate embedded Runge–Kutta scheme.
        !!  @param err Estimated error
        !!  @param fRho Current density
        !!  @param fP Current driving pressure
        !!  @param fR Current bubble radius
        !!  @param fV Current bubble velocity
        !!  @param fR0 Equilibrium bubble radius
        !!  @param fpb Internal bubble pressure
        !!  @param fpbdot Time-derivative of internal bubble pressure
        !!  @param alf bubble volume fraction
        !!  @param fntait Tait EOS parameter
        !!  @param fBtait Tait EOS parameter
        !!  @param f_bub_adv_src Source for bubble volume fraction
        !!  @param f_divu Divergence of velocity
        !!  @param bub_id Bubble identifier (EL)
        !!  @param fmass_v Current mass of vapour (EL)
        !!  @param fmass_g Current mass of gas (EL)
        !!  @param fbeta_c Mass transfer coefficient (EL)
        !!  @param fbeta_t Heat transfer coefficient (EL)
        !!  @param fCson Speed of sound (EL)
        !!  @param h Time step size
        !!  @param myR_tmp Bubble radius at each stage
        !!  @param myV_tmp Bubble radial velocity at each stage
        !!  @param myPb_tmp Internal bubble pressure at each stage (EL)
        !!  @param myMv_tmp Mass of vapor in the bubble at each stage (EL)
    subroutine s_advance_substep(err, fRho, fP, fR, fV, fR0, fpb, fpbdot, alf, &
                                 fntait, fBtait, f_bub_adv_src, f_divu, &
                                 bub_id, fmass_v, fmass_g, fbeta_c, &
                                 fbeta_t, fCson, h, &
                                 myR_tmp, myV_tmp, myPb_tmp, myMv_tmp)
 
# 697 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#ifdef _CRAYFTN
# 697 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!DIR$ INLINEALWAYS s_advance_substep
# 697 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#elif MFC_OpenACC
# 697 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$acc routine seq 
# 697 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#elif MFC_OpenMP
# 697 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
    
# 697 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
 
# 697 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$omp declare target device_type(any) 
# 697 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#endif
# 699 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
 
        real(wp), intent(OUT) :: err
        real(wp), intent(IN) :: fRho, fP, fR, fV, fR0, fpb, fpbdot, alf
        real(wp), intent(IN) :: fntait, fBtait, f_bub_adv_src, f_divu, h
        integer, intent(IN) :: bub_id
        real(wp), intent(IN) :: fmass_v, fmass_g, fbeta_c, fbeta_t, fCson
        real(wp), dimension(4), intent(OUT) :: myR_tmp, myV_tmp, myPb_tmp, myMv_tmp
 
        real(wp), dimension(4) :: myA_tmp, mydPbdt_tmp, mydMvdt_tmp
        real(wp) :: err_R, err_V
 
        mypb_tmp(1:4) = fpb
        mydpbdt_tmp(1:4) = fpbdot
 
        ! Stage 0
        myr_tmp(1) = fr
        myv_tmp(1) = fv
        if (bubbles_lagrange) then
            mypb_tmp(1) = fpb
            mymv_tmp(1) = fmass_v
            call s_advance_el(myr_tmp(1), myv_tmp(1), mypb_tmp(1), mymv_tmp(1), bub_id, &
                              fmass_g, fbeta_c, fbeta_t, mydpbdt_tmp(1), mydmvdt_tmp(1))
        end if
        mya_tmp(1) = f_rddot(frho, fp, myr_tmp(1), myv_tmp(1), fr0, &
                             mypb_tmp(1), mydpbdt_tmp(1), alf, fntait, fbtait, &
                             f_bub_adv_src, f_divu, &
                             fcson)
 
        ! Stage 1
        myr_tmp(2) = myr_tmp(1) + h*myv_tmp(1)
        if (myr_tmp(2) < 0._wp) then
            err = adap_dt_tol + 1._wp; return
        end if
        myv_tmp(2) = myv_tmp(1) + h*mya_tmp(1)
        if (bubbles_lagrange) then
            mypb_tmp(2) = mypb_tmp(1) + h*mydpbdt_tmp(1)
            mymv_tmp(2) = mymv_tmp(1) + h*mydmvdt_tmp(1)
            call s_advance_el(myr_tmp(2), myv_tmp(2), mypb_tmp(2), mymv_tmp(2), &
                              bub_id, fmass_g, fbeta_c, fbeta_t, mydpbdt_tmp(2), mydmvdt_tmp(2))
        end if
        mya_tmp(2) = f_rddot(frho, fp, myr_tmp(2), myv_tmp(2), fr0, &
                             mypb_tmp(2), mydpbdt_tmp(2), alf, fntait, fbtait, &
                             f_bub_adv_src, f_divu, &
                             fcson)
 
        ! Stage 2
        myr_tmp(3) = myr_tmp(1) + (h/4._wp)*(myv_tmp(1) + myv_tmp(2))
        if (myr_tmp(3) < 0._wp) then
            err = adap_dt_tol + 1._wp; return
        end if
        myv_tmp(3) = myv_tmp(1) + (h/4._wp)*(mya_tmp(1) + mya_tmp(2))
        if (bubbles_lagrange) then
            mypb_tmp(3) = mypb_tmp(1) + (h/4._wp)*(mydpbdt_tmp(1) + mydpbdt_tmp(2))
            mymv_tmp(3) = mymv_tmp(1) + (h/4._wp)*(mydmvdt_tmp(1) + mydmvdt_tmp(2))
            call s_advance_el(myr_tmp(3), myv_tmp(3), mypb_tmp(3), mymv_tmp(3), &
                              bub_id, fmass_g, fbeta_c, fbeta_t, mydpbdt_tmp(3), mydmvdt_tmp(3))
        end if
        mya_tmp(3) = f_rddot(frho, fp, myr_tmp(3), myv_tmp(3), fr0, &
                             mypb_tmp(3), mydpbdt_tmp(3), alf, fntait, fbtait, &
                             f_bub_adv_src, f_divu, &
                             fcson)
 
        ! Stage 3
        myr_tmp(4) = myr_tmp(1) + (h/6._wp)*(myv_tmp(1) + myv_tmp(2) + 4._wp*myv_tmp(3))
        if (myr_tmp(4) < 0._wp) then
            err = adap_dt_tol + 1._wp; return
        end if
        myv_tmp(4) = myv_tmp(1) + (h/6._wp)*(mya_tmp(1) + mya_tmp(2) + 4._wp*mya_tmp(3))
        if (bubbles_lagrange) then
            mypb_tmp(4) = mypb_tmp(1) + (h/6._wp)*(mydpbdt_tmp(1) + mydpbdt_tmp(2) + 4._wp*mydpbdt_tmp(3))
            mymv_tmp(4) = mymv_tmp(1) + (h/6._wp)*(mydmvdt_tmp(1) + mydmvdt_tmp(2) + 4._wp*mydmvdt_tmp(3))
            call s_advance_el(myr_tmp(4), myv_tmp(4), mypb_tmp(4), mymv_tmp(4), &
                              bub_id, fmass_g, fbeta_c, fbeta_t, mydpbdt_tmp(4), mydmvdt_tmp(4))
        end if
        mya_tmp(4) = f_rddot(frho, fp, myr_tmp(4), myv_tmp(4), fr0, &
                             mypb_tmp(4), mydpbdt_tmp(4), alf, fntait, fbtait, &
                             f_bub_adv_src, f_divu, &
                             fcson)
 
        ! Estimate error
        err_r = (-5._wp*h/24._wp)*(myv_tmp(2) + myv_tmp(3) - 2._wp*myv_tmp(4)) &
                /max(abs(myr_tmp(1)), abs(myr_tmp(4)))
        err_v = (-5._wp*h/24._wp)*(mya_tmp(2) + mya_tmp(3) - 2._wp*mya_tmp(4)) &
                /max(abs(myv_tmp(1)), abs(myv_tmp(4)))
        ! Error correction for non-oscillating bubbles
        if (max(abs(myv_tmp(1)), abs(myv_tmp(4))) < 1.e-12_wp) then
            err_v = 0._wp
        end if
        if (bubbles_lagrange .and. f_approx_equal(mya_tmp(1), 0._wp) .and. f_approx_equal(mya_tmp(2), 0._wp) .and. &
            f_approx_equal(mya_tmp(3), 0._wp) .and. f_approx_equal(mya_tmp(4), 0._wp)) then
            err_v = 0._wp
        end if
        err = sqrt((err_r**2._wp + err_v**2._wp)/2._wp)
 
    end subroutine s_advance_substep
 
    !>  Changes of pressure and vapor mass in the lagrange bubbles.
        !!  @param fR_tmp Bubble radius
        !!  @param fV_tmp Bubble radial velocity
        !!  @param fPb_tmp Internal bubble pressure
        !!  @param fMv_tmp Mass of vapor in the bubble
        !!  @param bub_id Bubble identifier
        !!  @param fmass_g Current mass of gas
        !!  @param fbeta_c Mass transfer coefficient
        !!  @param fbeta_t Heat transfer coefficient
        !!  @param fdPbdt_tmp Rate of change of the internal bubble pressure
        !!  @param advance_EL Rate of change of the mass of vapor in the bubble
    elemental subroutine s_advance_el(fR_tmp, fV_tmp, fPb_tmp, fMv_tmp, bub_id, &
                                      fmass_g, fbeta_c, fbeta_t, fdPbdt_tmp, advance_EL)
 
# 808 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#if MFC_OpenACC
# 808 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$acc routine seq 
# 808 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#elif MFC_OpenMP
# 808 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
    
# 808 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
 
# 808 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
!$omp declare target device_type(any) 
# 808 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles.fpp"
#endif
        real(wp), intent(IN) :: fr_tmp, fv_tmp, fpb_tmp, fmv_tmp
        real(wp), intent(IN) :: fmass_g, fbeta_c, fbeta_t
        integer, intent(IN) :: bub_id
        real(wp), intent(INOUT) :: fdpbdt_tmp
        real(wp), intent(out) :: advance_el
        real(wp) :: fvapflux, myr_m, mygamma_m
 
        call s_vflux(fr_tmp, fv_tmp, fpb_tmp, fmv_tmp, bub_id, fvapflux, fmass_g, fbeta_c, myr_m, mygamma_m)
        fdpbdt_tmp = f_bpres_dot(fvapflux, fr_tmp, fv_tmp, fpb_tmp, fmv_tmp, bub_id, fbeta_t, myr_m, mygamma_m)
        advance_el = 4._wp*pi*fr_tmp**2._wp*fvapflux
 
    end subroutine s_advance_el
 
end module m_bubbles