m_global_parameters.f90 Source File

!>
!! @file m_global_parameters.f90
!! @brief Contains module m_global_parameters
!! @author S. Bryngelson, K. Schimdmayer, V. Coralic, J. Meng, K. Maeda, T. Colonius
!! @version 1.0
!! @date JUNE 06 2019

! This file was generated by MFC. It is only used if the --hard-code option is
! passed to ./mfc.sh run, enabling a GPU-oriented optimization that hard-codes
! some case parameters in order to improve runtime performance.


!> @brief The module contains all of the parameters describing the program
!!              logistics, the computational domain and the simulation algorithm.
!!              Additionally, for the volume fraction model, physical parameters
!!              of each of the fluids present in the flow are located here. They
!!              include stiffened gas equation of state parameters, the Reynolds
!!              numbers and the Weber numbers.
module m_global_parameters

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

    use m_derived_types        !< Definitions of the derived types

#ifdef _OPENACC
    use openacc
#endif

    ! ==========================================================================

    implicit none

    ! Logistics ================================================================
    integer                    :: num_procs             !< Number of processors
    integer, parameter :: num_stcls_min = 5     !< Mininum # of stencils
    integer, parameter :: path_len = 400   !< Maximum path length
    integer, parameter :: name_len = 50    !< Maximum name length
    character, parameter :: dflt_char = ' '   !< Default string value
    real(kind(0d0)), parameter :: dflt_real = -1d6  !< Default real value
    integer, parameter :: dflt_int = -100  !< Default integer value
    real(kind(0d0)), parameter :: sgm_eps = 1d-16 !< Segmentation tolerance
    integer, parameter :: fourier_rings = 5     !< Fourier filter ring limit
    character(LEN=path_len)  :: case_dir              !< Case folder location
    logical                    :: run_time_info         !< Run-time output flag
    logical                    :: debug                 !< Debug mode print statements
    integer                    :: t_step_old            !< Existing IC/grid folder
    real(kind(0d0)), PARAMETER :: small_alf = 1d-7 !< Small alf tolerance
    ! ==========================================================================
    ! Computational Domain Parameters ==========================================
    integer :: proc_rank !< Rank of the local processor

    !> @name Number of cells in the x-, y- and z-directions, respectively
    !> @{
    integer :: m, n, p
    !> @}

    !> @name Global number of cells in each direction
    !> @{
    integer :: m_glb, n_glb, p_glb
    !> @}

    !> @name Cylindrical coordinates (either axisymmetric or full 3D)
    !> @{
    logical :: cyl_coord
    integer :: grid_geometry
    !> @}
!$acc declare create(cyl_coord, grid_geometry)

    !> @name Cell-boundary (CB) locations in the x-, y- and z-directions, respectively
    !> @{
    real(kind(0d0)), target, allocatable, dimension(:) :: x_cb, y_cb, z_cb
    !> @}

    !> @name Cell-center (CC) locations in the x-, y- and z-directions, respectively
    !> @{
    real(kind(0d0)), target, allocatable, dimension(:) :: x_cc, y_cc, z_cc
    !> @}

    !> @name Cell-width distributions in the x-, y- and z-directions, respectively
    !> @{
    real(kind(0d0)), target, allocatable, dimension(:) :: dx, dy, dz
    !> @}

    real(kind(0d0)) :: dt !< Size of the time-step

!$acc declare create(x_cb, y_cb, z_cb, x_cc, y_cc, z_cc, dx, dy, dz, dt, m, n, p)

    !> @name Starting time-step iteration, stopping time-step iteration and the number
    !! of time-step iterations between successive solution backups, respectively
    !> @{
    integer :: t_step_start, t_step_stop, t_step_save
    !> @}

    ! ==========================================================================

    ! Simulation Algorithm Parameters ==========================================
    integer         :: model_eqns     !< Multicomponent flow model
    integer         :: num_dims       !< Number of spatial dimensions
    integer         :: num_fluids     !< Number of fluids in the flow
    logical         :: adv_alphan     !< Advection of the last volume fraction
    logical         :: mpp_lim        !< Mixture physical parameters (MPP) limits
    integer         :: time_stepper   !< Time-stepper algorithm
    integer         :: weno_vars      !< WENO-reconstructed state variables type
    integer         :: weno_order     !< Order of the WENO reconstruction

    integer         :: weno_polyn     !< Degree of the WENO polynomials (polyn)

    real(kind(0d0)) :: weno_eps       !< Binding for the WENO nonlinear weights
    logical         :: mapped_weno    !< WENO with mapping of nonlinear weights
    logical         :: mp_weno        !< Monotonicity preserving (MP) WENO
    logical         :: weno_Re_flux   !< WENO reconstruct velocity gradients for viscous stress tensor
    integer         :: riemann_solver !< Riemann solver algorithm
    integer         :: wave_speeds    !< Wave speeds estimation method
    integer         :: avg_state      !< Average state evaluation method
    logical         :: alt_soundspeed !< Alternate mixture sound speed
    logical         :: null_weights   !< Null undesired WENO weights
    logical         :: mixture_err    !< Mixture properties correction
    logical         :: cu_tensor   

    integer         :: cpu_start, cpu_end, cpu_rate

!$acc declare create(weno_polyn, mpp_lim, num_fluids, model_eqns, num_dims, mixture_err, alt_soundspeed, avg_state, mapped_weno, mp_weno, weno_eps)

    !> @name Boundary conditions (BC) in the x-, y- and z-directions, respectively
    !> @{
    type(int_bounds_info) :: bc_x, bc_y, bc_z
    type(int_bounds_info) :: bc_x_glb, bc_y_glb, bc_z_glb
    !> @}

    logical :: parallel_io !< Format of the data files
    integer :: precision !< Precision of output files

    integer, allocatable, dimension(:) :: proc_coords !<
    !! Processor coordinates in MPI_CART_COMM

    integer, allocatable, dimension(:) :: start_idx !<
    !! Starting cell-center index of local processor in global grid

    type(mpi_io_var), public :: MPI_IO_DATA

    !> @name MPI info for parallel IO with Lustre file systems
    !> @{
    character(LEN=name_len) :: mpiiofs
    integer :: mpi_info_int
    !> @}

    integer, private :: ierr

    !> @name Annotations of the structure of the state and flux vectors in terms of the
    !! size and the configuration of the system of equations to which they belong
    !> @{
    integer           :: sys_size                  !< Number of unknowns in system of eqns.
    type(int_bounds_info) :: cont_idx                  !< Indexes of first & last continuity eqns.
    type(int_bounds_info) :: mom_idx                   !< Indexes of first & last momentum eqns.
    integer           :: E_idx                     !< Index of energy equation
    type(int_bounds_info) :: adv_idx                   !< Indexes of first & last advection eqns.
    type(int_bounds_info) :: internalEnergies_idx      !< Indexes of first & last internal energy eqns.
    type(bub_bounds_info) :: bub_idx               !< Indexes of first & last bubble variable eqns.
    integer               :: alf_idx               !< Index of void fraction
    integer           :: gamma_idx                 !< Index of specific heat ratio func. eqn.
    integer           :: pi_inf_idx                !< Index of liquid stiffness func. eqn.
    !> @}

!$acc declare create(bub_idx)

    !> @name The number of fluids, along with their identifying indexes, respectively,
    !! for which viscous effects, e.g. the shear and/or the volume Reynolds (Re)
    !! numbers, will be non-negligible.
    !> @{
    integer, dimension(2)   :: Re_size
    integer, allocatable, dimension(:, :) :: Re_idx
    !> @}
!$acc declare create(Re_size, Re_idx)

    !> @name The coordinate direction indexes and flags (flg), respectively, for which
    !! the configurations will be determined with respect to a working direction
    !! and that will be used to isolate the contributions, in that direction, in
    !! the dimensionally split system of equations.
    !> @{
    integer, dimension(3) :: dir_idx
    real(kind(0d0)), dimension(3) :: dir_flg
    !> @}
!$acc declare create(dir_idx, dir_flg)

    integer :: buff_size !<
    !! The number of cells that are necessary to be able to store enough boundary
    !! conditions data to march the solution in the physical computational domain
    !! to the next time-step.

    integer :: startx, starty, startz

!$acc declare create(sys_size, buff_size, startx, starty, startz, E_idx, gamma_idx, pi_inf_idx, alf_idx)

    ! END: Simulation Algorithm Parameters =====================================

    ! Fluids Physical Parameters ===============================================

    type(physical_parameters), dimension(num_fluids_max) :: fluid_pp !<
    !! Database of the physical parameters of each of the fluids that is present
    !! in the flow. These include the stiffened gas equation of state parameters,
    !! the Reynolds numbers and the Weber numbers.


    ! ==========================================================================

    integer :: fd_order !<
    !! The order of the finite-difference (fd) approximations of the first-order
    !! derivatives that need to be evaluated when the CoM or flow probe data
    !! files are to be written at each time step

    integer :: fd_number !<
    !! The finite-difference number is given by MAX(1, fd_order/2). Essentially,
    !! it is a measure of the half-size of the finite-difference stencil for the
    !! selected order of accuracy.

    logical, dimension(num_fluids_max) :: com_wrt, cb_wrt
    logical :: probe_wrt
    logical :: integral_wrt
    integer :: num_probes
    integer :: num_integrals
    type(probe_parameters), dimension(num_probes_max) :: probe
    type(integral_parameters), dimension(num_probes_max) :: integral
    real(kind(0d0)), dimension(5) :: threshold_mf
    integer, dimension(5) :: moment_order

    !> @name Reference density and pressure for Tait EOS
    !> @{
    real(kind(0d0)) :: rhoref, pref
    !> @}
!$acc declare create(rhoref, pref)

    !> @name Bubble modeling
    !> @{
    integer            :: nb       !< Number of eq. bubble sizes

    real(kind(0d0)) :: R0ref    !< Reference bubble size
    real(kind(0d0)) :: Ca       !< Cavitation number
    real(kind(0d0)) :: Web      !< Weber number
    real(kind(0d0)) :: Re_inv   !< Inverse Reynolds number
    real(kind(0d0)), dimension(:), allocatable :: weight !< Simpson quadrature weights
    real(kind(0d0)), dimension(:), allocatable :: R0     !< Bubble sizes
    real(kind(0d0)), dimension(:), allocatable :: V0     !< Bubble velocities
    logical         :: bubbles      !< Bubbles on/off
    logical         :: polytropic   !< Polytropic  switch
    logical         :: polydisperse !< Polydisperse bubbles

    integer         :: bubble_model !< Gilmore or Keller--Miksis bubble model
    integer         :: thermal      !< Thermal behavior. 1 = adiabatic, 2 = isotherm, 3 = transfer
    real(kind(0d0)), allocatable, dimension(:, :, :) :: ptil  !< Pressure modification
    real(kind(0d0)) :: poly_sigma  !< log normal sigma for polydisperse PDF

    logical         :: qbmm !< Quadrature moment method
    integer         :: nmom !< Number of carried moments per R0 location
    integer         :: nnode !< Number of QBMM nodes
    integer         :: nmomsp !< Number of moments required by ensemble-averaging
    integer         :: nmomtot !< Total number of carried moments moments/transport equations
    integer         :: R0_type

!$acc declare create(nb, R0ref, Ca, Web, Re_inv, weight, R0, V0, bubbles, polytropic, polydisperse, qbmm, nmom, nnode, nmomsp, nmomtot, R0_type, ptil, bubble_model, thermal, poly_sigma)

    type(scalar_field), allocatable, dimension(:) :: mom_sp
    type(scalar_field), allocatable, dimension(:, :, :) :: mom_3d
    !> @}
!$acc declare create(mom_sp, mom_3d)

    !> @name Physical bubble parameters (see Ando 2010, Preston 2007)
    !> @{
    real(kind(0d0)) :: R_n, R_v, phi_vn, phi_nv, Pe_c, Tw, pv, M_n, M_v
    real(kind(0d0)), dimension(:), allocatable :: k_n, k_v, pb0, mass_n0, mass_v0, Pe_T
    real(kind(0d0)), dimension(:), allocatable :: Re_trans_T, Re_trans_c, Im_trans_T, Im_trans_c, omegaN
    real(kind(0d0)) :: mul0, ss, gamma_v, mu_v
    real(kind(0d0)) :: gamma_m, gamma_n, mu_n
    real(kind(0d0)) :: gam
    !> @}
!$acc declare create(R_n, R_v, phi_vn, phi_nv, Pe_c, Tw, pv, M_n, M_v, k_n, k_v, pb0, mass_n0, mass_v0, Pe_T, Re_trans_T, Re_trans_c, Im_trans_T, Im_trans_c, omegaN , mul0, ss, gamma_v, mu_v, gamma_m, gamma_n, mu_n, gam)
    !> @name Acoustic monopole parameters
    !> @{
    logical         :: monopole !< Monopole switch
    type(mono_parameters), dimension(num_probes_max) :: mono !< Monopole parameters
    integer         :: num_mono !< Number of monopoles
    !> @}
!$acc declare create(monopole, mono, num_mono)


    real(kind(0d0)) :: mytime       !< Current simulation time
    real(kind(0d0)) :: finaltime    !< Final simulation time


    logical :: weno_flat, riemann_flat, cu_mpi

    ! ======================================================================

    ! Mathematical and Physical Constants ======================================
    real(kind(0d0)), PARAMETER :: pi = 3.141592653589793d0 !< Pi

    ! ==========================================================================

contains

    !> Assigns default values to the user inputs before reading
        !!  them in. This enables for an easier consistency check of
        !!  these parameters once they are read from the input file.
    subroutine s_assign_default_values_to_user_inputs() ! ------------------

        integer :: i, j !< Generic loop iterator

        ! Logistics
        case_dir = dflt_char
        run_time_info = .false.
        t_step_old = dflt_int

        debug = .false.

        ! Computational domain parameters
        m = dflt_int; n = dflt_int; p = dflt_int

        cyl_coord = .false.

        dt = dflt_real

        t_step_start = dflt_int
        t_step_stop = dflt_int
        t_step_save = dflt_int

        ! Simulation algorithm parameters
        model_eqns = dflt_int
        num_fluids = dflt_int
        adv_alphan = .false.
        mpp_lim = .false.
        time_stepper = dflt_int
        weno_vars = dflt_int
        weno_order = dflt_int
        weno_eps = dflt_real
        mapped_weno = .false.
        mp_weno = .false.
        weno_Re_flux = .false.
        riemann_solver = dflt_int
        wave_speeds = dflt_int
        avg_state = dflt_int
        alt_soundspeed = .false.
        null_weights = .false.
        mixture_err = .false.
        parallel_io = .false.
        precision = 2
        weno_flat = .true.
        riemann_flat = .true.
        cu_mpi = .true.

        bc_x%beg = dflt_int; bc_x%end = dflt_int
        bc_y%beg = dflt_int; bc_y%end = dflt_int
        bc_z%beg = dflt_int; bc_z%end = dflt_int

        ! Fluids physical parameters
        do i = 1, num_fluids_max
            fluid_pp(i)%gamma = dflt_real
            fluid_pp(i)%pi_inf = dflt_real
            fluid_pp(i)%Re(:) = dflt_real
            fluid_pp(i)%mul0 = dflt_real
            fluid_pp(i)%ss = dflt_real
            fluid_pp(i)%pv = dflt_real
            fluid_pp(i)%gamma_v = dflt_real
            fluid_pp(i)%M_v = dflt_real
            fluid_pp(i)%mu_v = dflt_real
            fluid_pp(i)%k_v = dflt_real
        end do

        ! Tait EOS
        rhoref = dflt_real
        pref = dflt_real

        ! Bubble modeling
        bubbles = .false.
        bubble_model = 1
        polytropic = .true.
        polydisperse = .false.
        thermal = dflt_int
        R0ref = dflt_real

        nb = dflt_int
        R0_type = dflt_int

        ! User inputs for qbmm for simulation code
        qbmm = .false.
        nnode = 1

        Ca = dflt_real
        Re_inv = dflt_real
        Web = dflt_real
        poly_sigma = dflt_real

        ! Monopole source
        monopole = .false.
        num_mono = 1

        cu_tensor = .false.

        do j = 1, num_probes_max
            do i = 1, 3
                mono(j)%loc(i) = dflt_real
            end do
            mono(j)%mag = dflt_real
            mono(j)%length = dflt_real
            mono(j)%delay = dflt_real
            mono(j)%dir = 1.d0
            mono(j)%npulse = 1.d0
            mono(j)%pulse = 1
            mono(j)%support = 1
        end do

        fd_order = dflt_int
        com_wrt = .false.
        cb_wrt = .false.
        probe_wrt = .false.
        integral_wrt = .false.
        num_probes = dflt_int
        num_integrals = dflt_int

        do i = 1, num_probes_max
            probe(i)%x = dflt_real
            probe(i)%y = dflt_real
            probe(i)%z = dflt_real
        end do

        do i = 1, num_probes_max
            integral(i)%xmin = dflt_real
            integral(i)%xmax = dflt_real
            integral(i)%ymin = dflt_real
            integral(i)%ymax = dflt_real
            integral(i)%ymin = dflt_real
            integral(i)%ymax = dflt_real
        end do

        do i = 1, 5
            threshold_mf(i) = dflt_real
            moment_order(i) = dflt_int
        end do

    end subroutine s_assign_default_values_to_user_inputs ! ----------------

    !>  The computation of parameters, the allocation of memory,
        !!      the association of pointers and/or the execution of any
        !!      other procedures that are necessary to setup the module.
    subroutine s_initialize_global_parameters_module() ! -------------------

        integer :: tmp_idx !< Temporary indexes storage
        integer :: i, j !< Generic loop iterators
        integer :: k !< Generic counter
        integer :: fac
        integer :: i1, i2, i3

        type(int_bounds_info) :: ix, iy, iz

        ! Determining the degree of the WENO polynomials
        weno_polyn = (weno_order - 1)/2
!$acc update device(weno_polyn)

!$acc update device(nb)


        ! Initializing the number of fluids for which viscous effects will
        ! be non-negligible, the number of distinctive material interfaces
        ! for which surface tension will be important and also, the number
        ! of fluids for which the physical and geometric curvatures of the
        ! interfaces will be computed
        Re_size = 0

        ! Gamma/Pi_inf Model ===============================================
        if (model_eqns == 1) then

            ! Setting number of fluids
            num_fluids = 1

            ! Annotating structure of the state and flux vectors belonging
            ! to the system of equations defined by the selected number of
            ! spatial dimensions and the gamma/pi_inf model
            cont_idx%beg = 1
            cont_idx%end = cont_idx%beg
            mom_idx%beg = cont_idx%end + 1
            mom_idx%end = cont_idx%end + num_dims
            E_idx = mom_idx%end + 1
            adv_idx%beg = E_idx + 1
            adv_idx%end = adv_idx%beg + 1
            gamma_idx = adv_idx%beg
            pi_inf_idx = adv_idx%end
            sys_size = adv_idx%end

            ! ==================================================================

            ! Volume Fraction Model ============================================
        else

            ! Annotating structure of the state and flux vectors belonging
            ! to the system of equations defined by the selected number of
            ! spatial dimensions and the volume fraction model
            if (model_eqns == 2) then
                cont_idx%beg = 1
                cont_idx%end = num_fluids
                mom_idx%beg = cont_idx%end + 1
                mom_idx%end = cont_idx%end + num_dims
                E_idx = mom_idx%end + 1
                adv_idx%beg = E_idx + 1
                adv_idx%end = E_idx + num_fluids


                sys_size = adv_idx%end

                if (bubbles) then
                    alf_idx = adv_idx%end
                else
                    alf_idx = 0
                end if

                if (bubbles) then
                    bub_idx%beg = sys_size + 1
                    if (qbmm) then
                        nmomsp = 4 !number of special moments
                        if (nnode == 4) then
                            nmom = 6
                            nmomtot = nmom*nb
                        end if
                        bub_idx%end = adv_idx%end + nb*nmom
                    else
                        if (.not. polytropic) then
                            bub_idx%end = sys_size + 4*nb
                        else
                            bub_idx%end = sys_size + 2*nb
                        end if
                    end if
                    sys_size = bub_idx%end
                    ! print*, 'alf idx', alf_idx
                    ! print*, 'bub -idx beg end', bub_idx%beg, bub_idx%end

                    allocate (weight(nb), R0(nb), V0(nb))
                    allocate (bub_idx%rs(nb), bub_idx%vs(nb))
                    allocate (bub_idx%ps(nb), bub_idx%ms(nb))

                    if (num_fluids == 1) then
                        gam = 1.d0/fluid_pp(num_fluids + 1)%gamma + 1.d0
                    else
                        gam = 1.d0/fluid_pp(num_fluids)%gamma + 1.d0
                    end if

                    if (qbmm) then
                        allocate (bub_idx%moms(nb, nmom))
                        do i = 1, nb
                            do j = 1, nmom
                                bub_idx%moms(i, j) = bub_idx%beg + (j - 1) + (i - 1)*nmom
                            end do
                            bub_idx%rs(i) = bub_idx%moms(i, 2)
                            bub_idx%vs(i) = bub_idx%moms(i, 3)
                        end do
                    else
                        do i = 1, nb
                            if (.not. polytropic) then
                                fac = 4
                            else
                                fac = 2
                            end if

                            bub_idx%rs(i) = bub_idx%beg + (i - 1)*fac
                            bub_idx%vs(i) = bub_idx%rs(i) + 1

                            if (.not. polytropic) then
                                bub_idx%ps(i) = bub_idx%vs(i) + 1
                                bub_idx%ms(i) = bub_idx%ps(i) + 1
                            end if
                        end do
                    end if


                    if (nb == 1) then
                        weight(:) = 1d0
                        R0(:) = 1d0
                        V0(:) = 1d0
                    else if (nb > 1) then
                        if (R0_type == 1) then
                            call s_simpson
                        else 
                            print*, 'Invalid R0 type - abort'
                            stop
                        end if
                        V0(:) = 1d0
                    else
                        stop 'Invalid value of nb'
                    end if

                    print *, 'R0 weights: ', weight(:)
                    print *, 'R0 abscissas: ', R0(:)

                    if (.not. polytropic) then
                        call s_initialize_nonpoly
                    else
                        rhoref = 1.d0
                        pref = 1.d0
                    end if
                end if


            else if (model_eqns == 3) then
                cont_idx%beg = 1
                cont_idx%end = num_fluids
                mom_idx%beg = cont_idx%end + 1
                mom_idx%end = cont_idx%end + num_dims
                E_idx = mom_idx%end + 1
                adv_idx%beg = E_idx + 1
                adv_idx%end = E_idx + num_fluids
                internalEnergies_idx%beg = adv_idx%end + 1
                internalEnergies_idx%end = adv_idx%end + num_fluids
                sys_size = internalEnergies_idx%end
            else if (model_eqns == 4) then
                cont_idx%beg = 1 ! one continuity equation
                cont_idx%end = 1 !num_fluids
                mom_idx%beg = cont_idx%end + 1 ! one momentum equation in each direction
                mom_idx%end = cont_idx%end + num_dims
                E_idx = mom_idx%end + 1 ! one energy equation
                adv_idx%beg = E_idx + 1
                adv_idx%end = adv_idx%beg !one volume advection equation
                alf_idx = adv_idx%end
                sys_size = alf_idx !adv_idx%end

                if (bubbles) then
                    bub_idx%beg = sys_size + 1
                    bub_idx%end = sys_size + 2*nb
                    if (.not. polytropic) then
                        bub_idx%end = sys_size + 4*nb
                    end if
                    sys_size = bub_idx%end

                    allocate (bub_idx%rs(nb), bub_idx%vs(nb))
                    allocate (bub_idx%ps(nb), bub_idx%ms(nb))
                    allocate (weight(nb), R0(nb), V0(nb))

                    do i = 1, nb
                        if (polytropic) then
                            fac = 2
                        else
                            fac = 4
                        end if

                        bub_idx%rs(i) = bub_idx%beg + (i - 1)*fac
                        bub_idx%vs(i) = bub_idx%rs(i) + 1

                        if (.not. polytropic) then
                            bub_idx%ps(i) = bub_idx%vs(i) + 1
                            bub_idx%ms(i) = bub_idx%ps(i) + 1
                        end if
                    end do
                    if (nb == 1) then
                        weight(:) = 1d0
                        R0(:) = 1d0
                        V0(:) = 0d0
                    else if (nb > 1) then
                        if (R0_type == 1) then
                            call s_simpson
                        else 
                            print*, 'Invalid R0 type - abort'
                            stop
                        end if
                        V0(:) = 1d0
                    else
                        stop 'Invalid value of nb'
                    end if

                    if (.not. polytropic) then
                        call s_initialize_nonpoly
                    else
                        rhoref = 1.d0
                        pref = 1.d0
                    end if
                end if
            end if

            ! Determining the number of fluids for which the shear and the
            ! volume Reynolds numbers, e.g. viscous effects, are important
            do i = 1, num_fluids
                if (fluid_pp(i)%Re(1) > 0) Re_size(1) = Re_size(1) + 1
                if (fluid_pp(i)%Re(2) > 0) Re_size(2) = Re_size(2) + 1
            end do

            ! Bookkeeping the indexes of any viscous fluids and any pairs of
            ! fluids whose interface will support effects of surface tension
            if (any(Re_size > 0)) then

                allocate (Re_idx(1:2, 1:maxval(Re_size)))

                k = 0
                do i = 1, num_fluids
                    if (fluid_pp(i)%Re(1) > 0) then
                        k = k + 1; Re_idx(1, k) = i
                    end if
                end do

                k = 0
                do i = 1, num_fluids
                    if (fluid_pp(i)%Re(2) > 0) then
                        k = k + 1; Re_idx(2, k) = i
                    end if
                end do

            end if

        end if
        ! END: Volume Fraction Model =======================================

        allocate (MPI_IO_DATA%view(1:sys_size))
        allocate (MPI_IO_DATA%var(1:sys_size))

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

!$acc update device(Re_size)
        ! Determining the number of cells that are needed in order to store
        ! sufficient boundary conditions data as to iterate the solution in
        ! the physical computational domain from one time-step iteration to
        ! the next one
        if (any(Re_size > 0)) then
            buff_size = 2*weno_polyn + 2
        else
            buff_size = weno_polyn + 2
        end if





        ! Configuring Coordinate Direction Indexes =========================
        if (bubbles) then
            ix%beg = -buff_size; iy%beg = 0; iz%beg = 0
            if (n > 0) iy%beg = -buff_size; if (p > 0) iz%beg = -buff_size
            ix%end = m - ix%beg; iy%end = n - iy%beg; iz%end = p - iz%beg
            allocate (ptil(ix%beg:ix%end, &
                           iy%beg:iy%end, &
                           iz%beg:iz%end))
        end if

        if (probe_wrt) then
            fd_number = max(1, fd_order/2)
            buff_size = buff_size + fd_number
        end if


        startx = -buff_size
        starty = 0
        startz = 0
        if(n > 0) then
            starty = -buff_size
        end if
        if(p > 0) then
            startz = -buff_size
        end if

!$acc update device(startx, starty, startz)

        if (cyl_coord .neqv. .true.) then ! Cartesian grid
            grid_geometry = 1
        elseif (cyl_coord .and. p == 0) then ! Axisymmetric cylindrical grid
            grid_geometry = 2
        else ! Fully 3D cylindrical grid
            grid_geometry = 3
        end if

        ! Allocating grid variables for the x-, y- and z-directions
        allocate (x_cb(-1 - buff_size:m + buff_size))
        allocate (x_cc(-buff_size:m + buff_size))
        allocate (dx(-buff_size:m + buff_size))

        if (n == 0) return; allocate (y_cb(-1 - buff_size:n + buff_size))
        allocate (y_cc(-buff_size:n + buff_size))
        allocate (dy(-buff_size:n + buff_size))

        if (p == 0) return; allocate (z_cb(-1 - buff_size:p + buff_size))
        allocate (z_cc(-buff_size:p + buff_size))
        allocate (dz(-buff_size:p + buff_size))

    end subroutine s_initialize_global_parameters_module ! -----------------

    !> Initializes non-polydisperse bubble modeling
    subroutine s_initialize_nonpoly
        integer :: ir
        real(kind(0.d0)) :: rhol0
        real(kind(0.d0)) :: pl0
        real(kind(0.d0)) :: uu
        real(kind(0.d0)) :: D_m
        real(kind(0.d0)) :: temp
        real(kind(0.d0)) :: omega_ref
        real(kind(0.d0)), dimension(Nb) :: chi_vw0
        real(kind(0.d0)), dimension(Nb) :: cp_m0
        real(kind(0.d0)), dimension(Nb) :: k_m0
        real(kind(0.d0)), dimension(Nb) :: rho_m0
        real(kind(0.d0)), dimension(Nb) :: x_vw

        real(kind(0.d0)), parameter :: k_poly = 1.d0 !<
            !! polytropic index used to compute isothermal natural frequency

        real(kind(0.d0)), parameter :: Ru = 8314.d0 !<
            !! universal gas constant

        rhol0 = rhoref
        pl0 = pref

        allocate (pb0(nb), mass_n0(nb), mass_v0(nb), Pe_T(nb))
        allocate (k_n(nb), k_v(nb), omegaN(nb))
        allocate (Re_trans_T(nb), Re_trans_c(nb), Im_trans_T(nb), Im_trans_c(nb))

        pb0(:) = dflt_real
        mass_n0(:) = dflt_real
        mass_v0(:) = dflt_real
        Pe_T(:) = dflt_real
        omegaN(:) = dflt_real

        mul0 = fluid_pp(1)%mul0
        ss = fluid_pp(1)%ss
        pv = fluid_pp(1)%pv
        gamma_v = fluid_pp(1)%gamma_v
        M_v = fluid_pp(1)%M_v
        mu_v = fluid_pp(1)%mu_v
        k_v(:) = fluid_pp(1)%k_v

        gamma_n = fluid_pp(2)%gamma_v
        M_n = fluid_pp(2)%M_v
        mu_n = fluid_pp(2)%mu_v
        k_n(:) = fluid_pp(2)%k_v

        gamma_m = gamma_n
        if (thermal == 2) gamma_m = 1.d0

        temp = 293.15d0
        D_m = 0.242d-4
        uu = DSQRT(pl0/rhol0)

        omega_ref = 3.d0*k_poly*Ca + 2.d0*(3.d0*k_poly - 1.d0)/Web

            !!! thermal properties !!!
        ! gas constants
        R_n = Ru/M_n
        R_v = Ru/M_v
        ! phi_vn & phi_nv (phi_nn = phi_vv = 1)
        phi_vn = (1.d0 + DSQRT(mu_v/mu_n)*(M_n/M_v)**(0.25d0))**2 &
                 /(DSQRT(8.d0)*DSQRT(1.d0 + M_v/M_n))
        phi_nv = (1.d0 + DSQRT(mu_n/mu_v)*(M_v/M_n)**(0.25d0))**2 &
                 /(DSQRT(8.d0)*DSQRT(1.d0 + M_n/M_v))
        ! internal bubble pressure
        pb0 = pl0 + 2.d0*ss/(R0ref*R0)

        ! mass fraction of vapor
        chi_vw0 = 1.d0/(1.d0 + R_v/R_n*(pb0/pv - 1.d0))
        ! specific heat for gas/vapor mixture
        cp_m0 = chi_vw0*R_v*gamma_v/(gamma_v - 1.d0) &
                + (1.d0 - chi_vw0)*R_n*gamma_n/(gamma_n - 1.d0)
        ! mole fraction of vapor
        x_vw = M_n*chi_vw0/(M_v + (M_n - M_v)*chi_vw0)
        ! thermal conductivity for gas/vapor mixture
        k_m0 = x_vw*k_v/(x_vw + (1.d0 - x_vw)*phi_vn) &
               + (1.d0 - x_vw)*k_n/(x_vw*phi_nv + 1.d0 - x_vw)
        ! mixture density
        rho_m0 = pv/(chi_vw0*R_v*temp)

        ! mass of gas/vapor computed using dimensional quantities
        mass_n0 = 4.d0*(pb0 - pv)*pi/(3.d0*R_n*temp*rhol0)*R0**3
        mass_v0 = 4.d0*pv*pi/(3.d0*R_v*temp*rhol0)*R0**3
        ! Peclet numbers
        Pe_T = rho_m0*cp_m0*uu*R0ref/k_m0
        Pe_c = uu*R0ref/D_m

        ! nondimensional properties
        R_n = rhol0*R_n*temp/pl0
        R_v = rhol0*R_v*temp/pl0
        k_n = k_n/k_m0
        k_v = k_v/k_m0
        pb0 = pb0/pl0
        pv = pv/pl0

        ! bubble wall temperature, normalized by T0, in the liquid
        ! keeps a constant (cold liquid assumption)
        Tw = 1.d0
        ! natural frequencies
        omegaN = DSQRT(3.d0*k_poly*Ca + 2.d0*(3.d0*k_poly - 1.d0)/(Web*R0))/R0

        pl0 = 1.d0
        do ir = 1, Nb
            call s_transcoeff(omegaN(ir)*R0(ir), Pe_T(ir)*R0(ir), &
                              Re_trans_T(ir), Im_trans_T(ir))
            call s_transcoeff(omegaN(ir)*R0(ir), Pe_c*R0(ir), &
                              Re_trans_c(ir), Im_trans_c(ir))
        end do
        Im_trans_T = 0d0
        Im_trans_c = 0d0

        rhoref = 1.d0
        pref = 1.d0
    end subroutine s_initialize_nonpoly

    !>  Computes transfer coefficient for non-polydisperse bubble modeling (Preston 2007)
        !!  @param omega Frequency
        !!  @param peclet Peclet number
        !!  @param Re_trans Real part of transfer coefficient
        !!  @param Im_trans Imaginary part of transfer coefficient
    subroutine s_transcoeff(omega, peclet, Re_trans, Im_trans)

        real(kind(0.d0)), intent(IN) :: omega
        real(kind(0.d0)), intent(IN) :: peclet
        real(kind(0.d0)), intent(OUT) :: Re_trans
        real(kind(0.d0)), intent(OUT) :: Im_trans
        complex :: trans, c1, c2, c3
        complex :: imag = (0., 1.)
        real(kind(0.d0)) :: f_transcoeff

        c1 = imag*omega*peclet
        c2 = CSQRT(c1)
        c3 = (CEXP(c2) - CEXP(-c2))/(CEXP(c2) + CEXP(-c2)) ! TANH(c2)
        trans = ((c2/c3 - 1.d0)**(-1) - 3.d0/c1)**(-1) ! transfer function

        Re_trans = dble(trans)
        Im_trans = aimag(trans)

    end subroutine s_transcoeff

    !> Initializes parallel infrastructure
    subroutine s_initialize_parallel_io() ! --------------------------------

        num_dims = 1 + min(1, n) + min(1, p)

        allocate (proc_coords(1:num_dims))

        if (parallel_io .neqv. .true.) return

#ifndef MFC_MPI

        print '(A)', '[m_global_parameters] s_initialize_parallel_io not supported without MPI.'

#else

    ! Option for Lustre file system (Darter/Comet/Stampede)
    write (mpiiofs, '(A)') '/lustre_'
    mpiiofs = trim(mpiiofs)


    call MPI_INFO_CREATE(mpi_info_int, ierr)
    call MPI_INFO_SET(mpi_info_int, 'romio_ds_write', 'disable', ierr)

        ! Option for UNIX file system (Hooke/Thomson)
        ! WRITE(mpiiofs, '(A)') '/ufs_'
        ! mpiiofs = TRIM(mpiiofs)
        ! mpi_info_int = MPI_INFO_NULL

        allocate (start_idx(1:num_dims))

#endif

    end subroutine s_initialize_parallel_io ! ------------------------------

    !> Module deallocation and/or disassociation procedures
    subroutine s_finalize_global_parameters_module() ! ---------------------

        integer :: i

        ! Deallocating the variables bookkeeping the indexes of any viscous
        ! fluids and any pairs of fluids whose interfaces supported effects
        ! of surface tension
        if (any(Re_size > 0)) deallocate (Re_idx)

        ! Deallocating grid variables for the x-, y- and z-directions
        deallocate (x_cb, x_cc, dx)

        if (n == 0) return; deallocate (y_cb, y_cc, dy)

        if (p == 0) return; deallocate (z_cb, z_cc, dz)

        deallocate (proc_coords)
        if (parallel_io) then
            deallocate (start_idx)
            do i = 1, sys_size
                MPI_IO_DATA%var(i)%sf => null()
            end do

            deallocate (MPI_IO_DATA%var)
            deallocate (MPI_IO_DATA%view)
        end if

    end subroutine s_finalize_global_parameters_module ! -------------------

    !>  Computes the bubble number density n from the conservative variables
        !!  \f$ n = \sqrt{ \frac{4 \pi}{3} } \frac{ nR^3}{\alpha} \f$
        !! @param vftmp is the void fraction
        !! @param nRtmp is the bubble number  density times the bubble radii
        !! @param ntmp is the output number bubble density
    subroutine s_comp_n_from_cons(vftmp, nRtmp, ntmp)
!$acc routine seq

        real(kind(0d0)), intent(IN) :: vftmp
        real(kind(0d0)), dimension(:), intent(IN) :: nRtmp
        real(kind(0d0)), intent(OUT) :: ntmp
        real(kind(0d0)) :: nR3
        integer :: i

        nR3 = 0d0
        do i = 1, nb
            nR3 = nR3 + weight(i)*(nRtmp(i)**3d0)
        end do


        !if (nR3 < 0d0) then
            ! DO i = 1,nb
            ! IF (nRtmp(i) < small_alf) THEN
            ! nRtmp(i) = small_alf
            ! END IF
            ! END DO
            ! nR3 = 1.d-12
            !print *, vftmp, nR3, nRtmp(:)
         !   stop 'nR3 is negative'
        !end if
        !if (vftmp < 0d0) then
            ! vftmp = small_alf
            ! ntmp = DSQRT( (4.d0*pi/3.d0)*nR3/1.d-12 )
            !print *, vftmp, nR3, nRtmp(:)
         !   stop 'vf negative'
        !end if

        ntmp = DSQRT((4.d0*pi/3.d0)*nR3/vftmp)

    end subroutine s_comp_n_from_cons

    !> Computes the bubble number density n from the primitive variables
        !!  \f$ n = \sqrt{ \frac{3}{4 \pi} } \frac{ \alpha }{ R^3} \f$
        !! @param vftmp is the void fraction
        !! @param Rtmp is the  bubble radii
        !! @param ntmp is the output number bubble density
    subroutine s_comp_n_from_prim(vftmp, Rtmp, ntmp)
!$acc routine seq

        real(kind(0.d0)), intent(IN) :: vftmp
        real(kind(0.d0)), dimension(:), intent(IN) :: Rtmp
        real(kind(0.d0)), intent(OUT) :: ntmp
        real(kind(0.d0)) :: R3
        integer :: i

        R3 = 0d0
        do i = 1, nb
            R3 = R3 + weight(i)*(Rtmp(i)**3d0)
        end do

        IF ( R3 < 0d0 ) THEN
            !PRINT*, vftmp, R3, Rtmp(:)
            STOP 'R3 is negative'
        END IF
        IF (vftmp < 0d0) THEN
            !PRINT*, vftmp, R3, Rtmp(:)
            STOP 'vf negative'
        END IF

        ntmp = (3.d0/(4.d0*pi))*vftmp/R3

    end subroutine s_comp_n_from_prim

    !> Computes the quadrature for polydisperse bubble populations
        !! @param func is the bubble dynamic variables for each bin
        !! @param mom is the computed moment
    subroutine s_quad(func, mom)
!$acc routine seq

        real(kind(0.d0)), dimension(:), intent(IN) :: func
        real(kind(0.d0)), intent(OUT) :: mom
        integer :: i

        mom = 0d0
        do i = 1, nb
            mom = mom + weight(i)*func(i)
        end do


    end subroutine s_quad

    !> Computes the Simpson weights for quadrature
    subroutine s_simpson

        integer :: ir
        real(kind(0.d0)) :: R0mn
        real(kind(0.d0)) :: R0mx
        real(kind(0.d0)) :: dphi
        real(kind(0.d0)) :: tmp
        real(kind(0.d0)) :: sd
        real(kind(0.d0)), dimension(nb) :: phi

        ! nondiml. min. & max. initial radii for numerical quadrature
        !sd   = 0.05D0
        !R0mn = 0.75D0
        !R0mx = 1.3D0

        !sd   = 0.3D0
        !R0mn = 0.3D0
        !R0mx = 6.D0

        !sd   = 0.7D0
        !R0mn = 0.12D0
        !R0mx = 150.D0

        sd = poly_sigma
        R0mn = 0.8d0*DEXP(-2.8d0*sd)
        R0mx = 0.2d0*DEXP(9.5d0*sd) + 1.d0

        ! phi = ln( R0 ) & return R0
        do ir = 1, nb
            phi(ir) = DLOG(R0mn) &
                      + dble(ir - 1)*DLOG(R0mx/R0mn)/dble(nb - 1)
            R0(ir) = DEXP(phi(ir))
        end do
        dphi = phi(2) - phi(1)

        ! weights for quadrature using Simpson's rule
        do ir = 2, nb - 1
            ! Gaussian
            tmp = DEXP(-0.5d0*(phi(ir)/sd)**2)/DSQRT(2.d0*pi)/sd
            if (mod(ir, 2) == 0) then
                weight(ir) = tmp*4.d0*dphi/3.d0
            else
                weight(ir) = tmp*2.d0*dphi/3.d0
            end if
        end do

        tmp = DEXP(-0.5d0*(phi(1)/sd)**2)/DSQRT(2.d0*pi)/sd
        weight(1) = tmp*dphi/3.d0
        tmp = DEXP(-0.5d0*(phi(nb)/sd)**2)/DSQRT(2.d0*pi)/sd
        weight(nb) = tmp*dphi/3.d0
    end subroutine s_simpson

end module m_global_parameters