m_global_parameters.fpp.f90

Go to the documentation of this file.

# 1 "/home/runner/work/MFC/MFC/src/post_process/m_global_parameters.fpp"
!>
!! @file
!! @brief Contains module m_global_parameters
 
# 1 "/home/runner/work/MFC/MFC/src/common/include/case.fpp" 1
! This file exists so that Fypp can be run without generating case.fpp files for
! each target. This is useful when generating documentation, for example. This
! should also let MFC be built with CMake directly, without invoking mfc.sh.
 
! For pre-process.
# 9 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
 
! For moving immersed boundaries in simulation
# 14 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
# 6 "/home/runner/work/MFC/MFC/src/post_process/m_global_parameters.fpp" 2
 
!> @brief Global parameters for the post-process: domain geometry, equation of state, and output database settings
module m_global_parameters
 
#ifdef MFC_MPI
    use mpi  !< message passing interface (mpi) module
#endif
 
    use m_derived_types
    use m_helper_basic
    use m_thermochem, only: num_species, species_names
 
    implicit none
 
    !> @name Logistics
    !> @{
    integer                 :: num_procs  !< Number of processors
    character(LEN=path_len) :: case_dir   !< Case folder location
    !> @}
 
    ! Computational Domain Parameters
 
    integer :: proc_rank  !< Rank of the local processor
    !> @name Number of cells in the x-, y- and z-coordinate directions
    !> @{
    integer :: m, m_root
    integer :: n
    integer :: p
    !> @}
 
    !> @name Max and min number of cells in a direction of each combination of x-,y-, and z-
    type(cell_num_bounds) :: cells_bounds
    integer(kind=8)       :: nglobal  !< Total number of cells in global domain
 
    !> @name Cylindrical coordinates (either axisymmetric or full 3D)
    !> @{
    logical :: cyl_coord
    integer :: grid_geometry
    !> @}
 
    !> @name Global number of cells in each direction
    !> @{
    integer :: m_glb, n_glb, p_glb
    !> @}
 
    integer :: num_dims  !< Number of spatial dimensions
    integer :: num_vels  !< Number of velocity components (different from num_dims for mhd)
    !> @name Cell-boundary locations in the x-, y- and z-coordinate directions
    !> @{
    real(wp), allocatable, dimension(:) :: x_cb, x_root_cb, y_cb, z_cb
    !> @}
 
    !> @name Cell-center locations in the x-, y- and z-coordinate directions
    !> @{
    real(wp), allocatable, dimension(:) :: x_cc, x_root_cc, y_cc, z_cc
    real(sp), allocatable, dimension(:) :: x_root_cc_s, x_cc_s
    !> @}
 
    !> Cell-width distributions in the x-, y- and z-coordinate directions
    !> @{
    real(wp), allocatable, dimension(:) :: dx, dy, dz
    !> @}
 
    integer :: buff_size     !< Number of ghost cells for boundary condition storage
    integer :: t_step_start  !< First time-step directory
    integer :: t_step_stop   !< Last time-step directory
    integer :: t_step_save   !< Interval between consecutive time-step directory
    !> @name IO options for adaptive time-stepping
    !> @{
    logical  :: cfl_adap_dt, cfl_const_dt, cfl_dt
    real(wp) :: t_save
    real(wp) :: t_stop
    real(wp) :: cfl_target
    integer  :: n_save
    integer  :: n_start
    !> @}
 
    ! NOTE: m_root, x_root_cb, x_root_cc = defragmented grid (1D only; equals m, x_cb, x_cc in serial)
 
    !> @name Simulation Algorithm Parameters
    !> @{
    integer            :: model_eqns                   !< Multicomponent flow model
    integer            :: num_fluids                   !< Number of different fluids present in the flow
    logical            :: relax                        !< phase change
    integer            :: relax_model                  !< Phase change relaxation model
    logical            :: mpp_lim                      !< Maximum volume fraction limiter
    integer            :: sys_size                     !< Number of unknowns in the system of equations
    integer            :: recon_type                   !< Which type of reconstruction to use
    integer            :: weno_order                   !< Order of accuracy for the WENO reconstruction
    integer            :: muscl_order                  !< Order of accuracy for the MUSCL reconstruction
    logical            :: mixture_err                  !< Mixture error limiter
    logical            :: alt_soundspeed               !< Alternate sound speed
    logical            :: mhd                          !< Magnetohydrodynamics
    logical            :: relativity                   !< Relativity for RMHD
    logical            :: hypoelasticity               !< Turn hypoelasticity on
    logical            :: hyperelasticity              !< Turn hyperelasticity on
    logical            :: elasticity                   !< elasticity modeling, true for hyper or hypo
    integer            :: b_size                       !< Number of components in the b tensor
    integer            :: tensor_size                  !< Number of components in the nonsymmetric tensor
    logical            :: cont_damage                  !< Continuum damage modeling
    logical            :: hyper_cleaning               !< Hyperbolic cleaning for MHD
    logical            :: igr                          !< enable IGR
    integer            :: igr_order                    !< IGR reconstruction order
    logical, parameter :: chemistry = .false.  !< Chemistry modeling
    !> @}
 
    integer :: avg_state  !< Average state evaluation method
    !> @name Annotations of the structure, i.e. the organization, of the state vectors
    !> @{
    type(eqn_idx_info)  :: eqn_idx   !< All conserved-variable equation index ranges and scalars.
    type(qbmm_idx_info) :: qbmm_idx  !< QBMM moment index mappings.
    integer             :: beta_idx  !< Index of lagrange bubbles beta
    !> @}
 
    ! Cell Indices for the (local) interior points (O-m, O-n, 0-p). Stands for "InDices With BUFFer".
    type(int_bounds_info) :: idwint(1:3)
 
    ! Cell indices (InDices With BUFFer): includes buffer in simulation only
    type(int_bounds_info) :: idwbuff(1:3)
    integer               :: num_bc_patches
    logical               :: bc_io
    !> @name Boundary conditions in the x-, y- and z-coordinate directions
    !> @{
    type(int_bounds_info) :: bc_x, bc_y, bc_z
    !> @}
 
    integer                            :: shear_num  !< Number of shear stress components
    integer, dimension(3)              :: shear_indices  !< Indices of the stress components that represent shear stress
    integer                            :: shear_bc_flip_num  !< Number of shear stress components to reflect for boundary conditions
    integer, dimension(3, 2)           :: shear_bc_flip_indices  !< Shear stress BC reflection indices (1:3, 1:shear_BC_flip_num)
    logical                            :: parallel_io  !< Format of the data files
    logical                            :: sim_data
    logical                            :: file_per_process  !< output format
    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
    integer                            :: num_ibs  !< Number of immersed boundaries
#ifdef MFC_MPI
    type(mpi_io_var), public                      :: mpi_io_data
    type(mpi_io_ib_var), public                   :: mpi_io_ib_data
    type(mpi_io_levelset_var), public             :: mpi_io_levelset_data
    type(mpi_io_levelset_norm_var), public        :: mpi_io_levelsetnorm_data
    real(wp), allocatable, dimension(:,:), public :: mpi_io_data_lg_bubbles
#endif
 
    !> @name MPI info for parallel IO with Lustre file systems
    !> @{
    character(LEN=name_len) :: mpiiofs
    integer                 :: mpi_info_int
    !> @}
 
    type(physical_parameters), dimension(num_fluids_max) :: fluid_pp  !< Stiffened gas EOS parameters and Reynolds numbers per fluid
    ! Subgrid Bubble Parameters
    type(subgrid_bubble_physical_parameters) :: bub_pp
    real(wp), allocatable, dimension(:)      :: adv  !< Advection variables
    ! Formatted Database File(s) Structure Parameters
 
    integer               :: format                                    !< Format of the database file(s)
    integer               :: precision                                 !< Floating point precision of the database file(s)
    logical               :: down_sample                               !< down sampling of the database file(s)
    logical               :: output_partial_domain                     !< Specify portion of domain to output for post-processing
    type(bounds_info)     :: x_output, y_output, z_output              !< Portion of domain to output for post-processing
    type(int_bounds_info) :: x_output_idx, y_output_idx, z_output_idx  !< Indices of domain to output for post-processing
    !> @name Size of the ghost zone layer in the x-, y- and z-coordinate directions. The definition of the ghost zone layers is only
    !! necessary when using the Silo database file format in multidimensions. These zones provide VisIt with the subdomain
    !! connectivity information that it requires in order to produce smooth plots.
    !> @{
    type(int_bounds_info) :: offset_x, offset_y, offset_z
    !> @}
 
    !> @name The list of all possible flow variables that may be written to a database file. It includes partial densities, density,
    !! momentum, velocity, energy, pressure, volume fraction(s), specific heat ratio function, specific heat ratio, liquid stiffness
    !! function, liquid stiffness, primitive variables, conservative variables, speed of sound, the vorticity, and the numerical
    !! Schlieren function.
    !> @{
    logical, dimension(num_fluids_max) :: alpha_rho_wrt
    logical                            :: rho_wrt
    logical, dimension(3)              :: mom_wrt
    logical, dimension(3)              :: vel_wrt
    integer                            :: flux_lim
    logical, dimension(3)              :: flux_wrt
    logical                            :: e_wrt
    logical, dimension(num_fluids_max) :: alpha_rho_e_wrt
    logical                            :: fft_wrt
    logical                            :: dummy  !< AMDFlang workaround for case-optimization + GPU-kernel bug
    logical                            :: pres_wrt
    logical, dimension(num_fluids_max) :: alpha_wrt
    logical                            :: gamma_wrt
    logical                            :: heat_ratio_wrt
    logical                            :: pi_inf_wrt
    logical                            :: pres_inf_wrt
    logical                            :: prim_vars_wrt
    logical                            :: cons_vars_wrt
    logical                            :: c_wrt
    logical, dimension(3)              :: omega_wrt
    logical                            :: qm_wrt
    logical                            :: liutex_wrt
    logical                            :: schlieren_wrt
    logical                            :: cf_wrt
    logical                            :: ib
    logical                            :: ib_state_wrt
    logical                            :: chem_wrt_y(1:num_species)
    logical                            :: chem_wrt_t
    logical                            :: lag_header
    logical                            :: lag_txt_wrt
    logical                            :: lag_db_wrt
    logical                            :: lag_id_wrt
    logical                            :: lag_pos_wrt
    logical                            :: lag_pos_prev_wrt
    logical                            :: lag_vel_wrt
    logical                            :: lag_rad_wrt
    logical                            :: lag_rvel_wrt
    logical                            :: lag_r0_wrt
    logical                            :: lag_rmax_wrt
    logical                            :: lag_rmin_wrt
    logical                            :: lag_dphidt_wrt
    logical                            :: lag_pres_wrt
    logical                            :: lag_mv_wrt
    logical                            :: lag_mg_wrt
    logical                            :: lag_betat_wrt
    logical                            :: lag_betac_wrt
    !> @}
 
    real(wp), dimension(num_fluids_max) :: schlieren_alpha  !< Per-fluid Schlieren intensity amplitude coefficients
    integer                             :: fd_order         !< Finite-difference order for vorticity and Schlieren derivatives
    integer                             :: fd_number        !< Finite-difference half-stencil size: MAX(1, fd_order/2)
    !> @name Reference parameters for Tait EOS
    !> @{
    real(wp) :: rhoref, pref
    !> @}
 
    type(chemistry_parameters) :: chem_params
    !> @name Bubble modeling variables and parameters
    !> @{
    integer :: nb
    real(wp) :: eu, ca, web, re_inv
    real(wp), dimension(:), allocatable :: weight, r0
    logical :: bubbles_euler
    logical :: qbmm
    logical :: polytropic
    logical :: polydisperse
    logical :: adv_n
    integer :: thermal  !< 1 = adiabatic, 2 = isotherm, 3 = transfer
    real(wp) :: phi_vg, phi_gv, pe_c, tw, k_vl, k_gl
    real(wp) :: gam_m
    real(wp), dimension(:), allocatable :: pb0, mass_g0, mass_v0, pe_t, k_v, k_g
    real(wp), dimension(:), allocatable :: re_trans_t, re_trans_c, im_trans_t, im_trans_c, omegan
    real(wp) :: r0ref, p0ref, rho0ref, t0ref, ss, pv, vd, mu_l, mu_v, mu_g, gam_v, gam_g, m_v, m_g, cp_v, cp_g, r_v, r_g
    real(wp) :: g
    real(wp) :: poly_sigma
    real(wp) :: sigr
    integer :: nmom
    !> @}
 
    !> @name surface tension coefficient
    !> @{
    real(wp) :: sigma
    logical  :: surface_tension
    !> @}
 
    !> @name Lagrangian bubbles
    !> @{
    logical :: bubbles_lagrange
    !> @}
 
    real(wp) :: bx0                       !< Constant magnetic field in the x-direction (1D)
    real(wp) :: wall_time, wall_time_avg  !< Wall time measurements
 
contains
 
    !> Assigns default values to user inputs prior to reading them in. This allows for an easier consistency check of these
    !! parameters once they are read from the input file.
    impure subroutine s_assign_default_values_to_user_inputs
 
        integer :: i  !< Generic loop iterator
        ! Logistics
 
        case_dir = '.'
 
        ! Computational domain parameters
        m = dflt_int; n = 0; p = 0
        call s_update_cell_bounds(cells_bounds, m, n, p)
 
        m_root = dflt_int
        cyl_coord = .false.
 
        t_step_start = dflt_int
        t_step_stop = dflt_int
        t_step_save = dflt_int
 
        cfl_adap_dt = .false.
        cfl_const_dt = .false.
        cfl_dt = .false.
        cfl_target = dflt_real
        t_save = dflt_real
        n_start = dflt_int
        t_stop = dflt_real
 
        ! Simulation algorithm parameters
        model_eqns = dflt_int
        num_fluids = dflt_int
        recon_type = weno_type
        weno_order = dflt_int
        muscl_order = dflt_int
        mixture_err = .false.
        alt_soundspeed = .false.
        relax = .false.
        relax_model = dflt_int
 
        mhd = .false.
        relativity = .false.
 
        hypoelasticity = .false.
        hyperelasticity = .false.
        elasticity = .false.
        b_size = dflt_int
        tensor_size = dflt_int
        cont_damage = .false.
        hyper_cleaning = .false.
        igr = .false.
 
        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
        bc_io = .false.
        num_bc_patches = dflt_int
 
# 333 "/home/runner/work/MFC/MFC/src/post_process/m_global_parameters.fpp"
# 334 "/home/runner/work/MFC/MFC/src/post_process/m_global_parameters.fpp"
                bc_x%vb1 = 0._wp
                bc_x%ve1 = 0._wp
# 334 "/home/runner/work/MFC/MFC/src/post_process/m_global_parameters.fpp"
                bc_x%vb2 = 0._wp
                bc_x%ve2 = 0._wp
# 334 "/home/runner/work/MFC/MFC/src/post_process/m_global_parameters.fpp"
                bc_x%vb3 = 0._wp
                bc_x%ve3 = 0._wp
# 337 "/home/runner/work/MFC/MFC/src/post_process/m_global_parameters.fpp"
# 333 "/home/runner/work/MFC/MFC/src/post_process/m_global_parameters.fpp"
# 334 "/home/runner/work/MFC/MFC/src/post_process/m_global_parameters.fpp"
                bc_y%vb1 = 0._wp
                bc_y%ve1 = 0._wp
# 334 "/home/runner/work/MFC/MFC/src/post_process/m_global_parameters.fpp"
                bc_y%vb2 = 0._wp
                bc_y%ve2 = 0._wp
# 334 "/home/runner/work/MFC/MFC/src/post_process/m_global_parameters.fpp"
                bc_y%vb3 = 0._wp
                bc_y%ve3 = 0._wp
# 337 "/home/runner/work/MFC/MFC/src/post_process/m_global_parameters.fpp"
# 333 "/home/runner/work/MFC/MFC/src/post_process/m_global_parameters.fpp"
# 334 "/home/runner/work/MFC/MFC/src/post_process/m_global_parameters.fpp"
                bc_z%vb1 = 0._wp
                bc_z%ve1 = 0._wp
# 334 "/home/runner/work/MFC/MFC/src/post_process/m_global_parameters.fpp"
                bc_z%vb2 = 0._wp
                bc_z%ve2 = 0._wp
# 334 "/home/runner/work/MFC/MFC/src/post_process/m_global_parameters.fpp"
                bc_z%vb3 = 0._wp
                bc_z%ve3 = 0._wp
# 337 "/home/runner/work/MFC/MFC/src/post_process/m_global_parameters.fpp"
# 338 "/home/runner/work/MFC/MFC/src/post_process/m_global_parameters.fpp"
 
# 340 "/home/runner/work/MFC/MFC/src/post_process/m_global_parameters.fpp"
            bc_y%isothermal_in = .false.
            bc_y%isothermal_out = .false.
            bc_y%Twall_in = dflt_real
            bc_y%Twall_out = dflt_real
# 340 "/home/runner/work/MFC/MFC/src/post_process/m_global_parameters.fpp"
            bc_x%isothermal_in = .false.
            bc_x%isothermal_out = .false.
            bc_x%Twall_in = dflt_real
            bc_x%Twall_out = dflt_real
# 340 "/home/runner/work/MFC/MFC/src/post_process/m_global_parameters.fpp"
            bc_z%isothermal_in = .false.
            bc_z%isothermal_out = .false.
            bc_z%Twall_in = dflt_real
            bc_z%Twall_out = dflt_real
# 345 "/home/runner/work/MFC/MFC/src/post_process/m_global_parameters.fpp"
 
        chem_params%gamma_method = 1
        chem_params%transport_model = 1
 
        ! 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)%cv = 0._wp
            fluid_pp(i)%qv = 0._wp
            fluid_pp(i)%qvp = 0._wp
            fluid_pp(i)%G = dflt_real
        end do
 
        ! Subgrid bubble parameters
        bub_pp%R0ref = dflt_real; r0ref = dflt_real
        bub_pp%p0ref = dflt_real; p0ref = dflt_real
        bub_pp%rho0ref = dflt_real; rho0ref = dflt_real
        bub_pp%T0ref = dflt_real; t0ref = dflt_real
        bub_pp%ss = dflt_real; ss = dflt_real
        bub_pp%pv = dflt_real; pv = dflt_real
        bub_pp%vd = dflt_real; vd = dflt_real
        bub_pp%mu_l = dflt_real; mu_l = dflt_real
        bub_pp%mu_v = dflt_real; mu_v = dflt_real
        bub_pp%mu_g = dflt_real; mu_g = dflt_real
        bub_pp%gam_v = dflt_real; gam_v = dflt_real
        bub_pp%gam_g = dflt_real; gam_g = dflt_real
        bub_pp%M_v = dflt_real; m_v = dflt_real
        bub_pp%M_g = dflt_real; m_g = dflt_real
        bub_pp%k_v = dflt_real
        bub_pp%k_g = dflt_real
        bub_pp%cp_v = dflt_real; cp_v = dflt_real
        bub_pp%cp_g = dflt_real; cp_g = dflt_real
        bub_pp%R_v = dflt_real; r_v = dflt_real
        bub_pp%R_g = dflt_real; r_g = dflt_real
 
        ! Formatted database file(s) structure parameters
        format = dflt_int
 
        precision = dflt_int
        down_sample = .false.
 
        alpha_rho_wrt = .false.
        alpha_rho_e_wrt = .false.
        rho_wrt = .false.
        mom_wrt = .false.
        vel_wrt = .false.
        chem_wrt_y = .false.
        chem_wrt_t = .false.
        flux_lim = dflt_int
        flux_wrt = .false.
        parallel_io = .false.
        file_per_process = .false.
        e_wrt = .false.
        fft_wrt = .false.
        dummy = .false.
        pres_wrt = .false.
        alpha_wrt = .false.
        gamma_wrt = .false.
        heat_ratio_wrt = .false.
        pi_inf_wrt = .false.
        pres_inf_wrt = .false.
        prim_vars_wrt = .false.
        cons_vars_wrt = .false.
        c_wrt = .false.
        omega_wrt = .false.
        qm_wrt = .false.
        liutex_wrt = .false.
        schlieren_wrt = .false.
        sim_data = .false.
        cf_wrt = .false.
        ib = .false.
        ib_state_wrt = .false.
        lag_txt_wrt = .false.
        lag_header = .true.
        lag_db_wrt = .false.
        lag_id_wrt = .true.
        lag_pos_wrt = .true.
        lag_pos_prev_wrt = .false.
        lag_vel_wrt = .true.
        lag_rad_wrt = .true.
        lag_rvel_wrt = .false.
        lag_r0_wrt = .false.
        lag_rmax_wrt = .false.
        lag_rmin_wrt = .false.
        lag_dphidt_wrt = .false.
        lag_pres_wrt = .false.
        lag_mv_wrt = .false.
        lag_mg_wrt = .false.
        lag_betat_wrt = .false.
        lag_betac_wrt = .false.
 
        schlieren_alpha = dflt_real
 
        fd_order = dflt_int
        avg_state = dflt_int
 
        ! Tait EOS
        rhoref = dflt_real
        pref = dflt_real
 
        ! Bubble modeling
        bubbles_euler = .false.
        qbmm = .false.
        r0ref = dflt_real
        nb = dflt_int
        polydisperse = .false.
        poly_sigma = dflt_real
        sigr = dflt_real
        sigma = dflt_real
        surface_tension = .false.
        adv_n = .false.
 
        ! Lagrangian bubbles modeling
        bubbles_lagrange = .false.
 
        ! IBM
        num_ibs = dflt_int
 
        ! Output partial domain
        output_partial_domain = .false.
        x_output%beg = dflt_real
        x_output%end = dflt_real
        y_output%beg = dflt_real
        y_output%end = dflt_real
        z_output%beg = dflt_real
        z_output%end = dflt_real
 
        ! MHD
        bx0 = dflt_real
 
    end subroutine s_assign_default_values_to_user_inputs
 
    !> Computation of parameters, allocation procedures, and/or any other tasks needed to properly setup the module
    impure subroutine s_initialize_global_parameters_module
 
        integer :: i, j, fac
 
        ! Setting m_root equal to m in the case of a 1D serial simulation
 
        if (n == 0) m_root = m_glb
 
        ! 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
            eqn_idx%cont%beg = 1
            eqn_idx%cont%end = eqn_idx%cont%beg
            eqn_idx%mom%beg = eqn_idx%cont%end + 1
            eqn_idx%mom%end = eqn_idx%cont%end + num_vels
            eqn_idx%E = eqn_idx%mom%end + 1
            eqn_idx%adv%beg = eqn_idx%E + 1
            eqn_idx%adv%end = eqn_idx%adv%beg + 1
            eqn_idx%gamma = eqn_idx%adv%beg
            eqn_idx%pi_inf = eqn_idx%adv%end
            sys_size = eqn_idx%adv%end
 
            ! Volume Fraction Model (5-equation model)
        else if (model_eqns == 2) then
            ! 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
            eqn_idx%cont%beg = 1
            eqn_idx%cont%end = num_fluids
            eqn_idx%mom%beg = eqn_idx%cont%end + 1
            eqn_idx%mom%end = eqn_idx%cont%end + num_vels
            eqn_idx%E = eqn_idx%mom%end + 1
 
            if (igr) then
                ! Volume fractions are stored in the indices immediately following the energy equation. IGR tracks a total of (N-1)
                ! volume fractions for N fluids, hence the "-1" in eqn_idx%adv%end. If num_fluids = 1 then eqn_idx%adv%end <
                ! eqn_idx%adv%beg, which skips all loops over the volume fractions since there is no volume fraction to track
                eqn_idx%adv%beg = eqn_idx%E + 1  ! Alpha for fluid 1
                eqn_idx%adv%end = eqn_idx%E + num_fluids - 1
            else
                ! Volume fractions are stored in the indices immediately following the energy equation. WENO/MUSCL + Riemann tracks
                ! a total of (N) volume fractions for N fluids, hence the lack of "-1" in eqn_idx%adv%end
                eqn_idx%adv%beg = eqn_idx%E + 1
                eqn_idx%adv%end = eqn_idx%E + num_fluids
            end if
 
            sys_size = eqn_idx%adv%end
 
            if (bubbles_euler) then
                eqn_idx%alf = eqn_idx%adv%end
            else
                eqn_idx%alf = 1
            end if
 
            if (qbmm) then
                nmom = 6
            end if
 
            if (bubbles_euler) then
                eqn_idx%bub%beg = sys_size + 1
                if (qbmm) then
                    eqn_idx%bub%end = eqn_idx%adv%end + nb*nmom
                else
                    if (.not. polytropic) then
                        eqn_idx%bub%end = sys_size + 4*nb
                    else
                        eqn_idx%bub%end = sys_size + 2*nb
                    end if
                end if
                sys_size = eqn_idx%bub%end
 
                if (adv_n) then
                    eqn_idx%n = eqn_idx%bub%end + 1
                    sys_size = eqn_idx%n
                end if
 
                allocate (qbmm_idx%rs(nb), qbmm_idx%vs(nb))
                allocate (qbmm_idx%ps(nb), qbmm_idx%ms(nb))
 
                if (qbmm) then
                    allocate (qbmm_idx%moms(nb, nmom))
                    do i = 1, nb
                        do j = 1, nmom
                            qbmm_idx%moms(i, j) = eqn_idx%bub%beg + (j - 1) + (i - 1)*nmom
                        end do
                        qbmm_idx%rs(i) = qbmm_idx%moms(i, 2)
                        qbmm_idx%vs(i) = qbmm_idx%moms(i, 3)
                    end do
                else
                    do i = 1, nb
                        if (polytropic .neqv. .true.) then
                            fac = 4
                        else
                            fac = 2
                        end if
 
                        qbmm_idx%rs(i) = eqn_idx%bub%beg + (i - 1)*fac
                        qbmm_idx%vs(i) = qbmm_idx%rs(i) + 1
 
                        if (polytropic .neqv. .true.) then
                            qbmm_idx%ps(i) = qbmm_idx%vs(i) + 1
                            qbmm_idx%ms(i) = qbmm_idx%ps(i) + 1
                        end if
                    end do
                end if
            end if
 
            if (bubbles_lagrange) then
                beta_idx = sys_size + 1
                sys_size = beta_idx
            end if
 
            if (mhd) then
                eqn_idx%B%beg = sys_size + 1
                if (n == 0) then
                    eqn_idx%B%end = sys_size + 2  ! 1D: By, Bz
                else
                    eqn_idx%B%end = sys_size + 3  ! 2D/3D: Bx, By, Bz
                end if
                sys_size = eqn_idx%B%end
            end if
 
            ! Volume Fraction Model (6-equation model)
        else if (model_eqns == 3) then
            ! 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
            eqn_idx%cont%beg = 1
            eqn_idx%cont%end = num_fluids
            eqn_idx%mom%beg = eqn_idx%cont%end + 1
            eqn_idx%mom%end = eqn_idx%cont%end + num_vels
            eqn_idx%E = eqn_idx%mom%end + 1
            eqn_idx%adv%beg = eqn_idx%E + 1
            eqn_idx%adv%end = eqn_idx%E + num_fluids
            eqn_idx%int_en%beg = eqn_idx%adv%end + 1
            eqn_idx%int_en%end = eqn_idx%adv%end + num_fluids
            sys_size = eqn_idx%int_en%end
            eqn_idx%alf = 1  ! dummy, cannot actually have a void fraction
        else if (model_eqns == 4) then
            eqn_idx%cont%beg = 1  ! one continuity equation
            eqn_idx%cont%end = 1  ! num_fluids
            eqn_idx%mom%beg = eqn_idx%cont%end + 1  ! one momentum equation in each
            eqn_idx%mom%end = eqn_idx%cont%end + num_vels
            eqn_idx%E = eqn_idx%mom%end + 1  ! one energy equation
            eqn_idx%adv%beg = eqn_idx%E + 1
            eqn_idx%adv%end = eqn_idx%adv%beg  ! one volume advection equation
            eqn_idx%alf = eqn_idx%adv%end
            sys_size = eqn_idx%alf  ! eqn_idx%adv%end
 
            if (bubbles_euler) then
                eqn_idx%bub%beg = sys_size + 1
                eqn_idx%bub%end = sys_size + 2*nb
                if (polytropic .neqv. .true.) then
                    eqn_idx%bub%end = sys_size + 4*nb
                end if
                sys_size = eqn_idx%bub%end
 
                allocate (qbmm_idx%rs(nb), qbmm_idx%vs(nb))
                allocate (qbmm_idx%ps(nb), qbmm_idx%ms(nb))
                allocate (weight(nb), r0(nb))
 
                do i = 1, nb
                    if (polytropic .neqv. .true.) then
                        fac = 4
                    else
                        fac = 2
                    end if
 
                    qbmm_idx%rs(i) = eqn_idx%bub%beg + (i - 1)*fac
                    qbmm_idx%vs(i) = qbmm_idx%rs(i) + 1
 
                    if (polytropic .neqv. .true.) then
                        qbmm_idx%ps(i) = qbmm_idx%vs(i) + 1
                        qbmm_idx%ms(i) = qbmm_idx%ps(i) + 1
                    end if
                end do
 
                if (nb == 1) then
                    weight(:) = 1._wp
                    r0(:) = 1._wp
                else if (nb < 1) then
                    stop 'Invalid value of nb'
                end if
 
                if (polytropic) then
                    rhoref = 1._wp
                    pref = 1._wp
                end if
            end if
        end if
 
        if (model_eqns == 2 .or. model_eqns == 3) then
            if (hypoelasticity .or. hyperelasticity) then
                elasticity = .true.
                eqn_idx%stress%beg = sys_size + 1
                eqn_idx%stress%end = sys_size + (num_dims*(num_dims + 1))/2
                if (cyl_coord) eqn_idx%stress%end = eqn_idx%stress%end + 1
                ! number of stresses is 1 in 1D, 3 in 2D, 4 in 2D-Axisym, 6 in 3D
                sys_size = eqn_idx%stress%end
 
                ! shear stress index is 2 for 2D and 2,4,5 for 3D
                if (num_dims == 1) then
                    shear_num = 0
                else if (num_dims == 2) then
                    shear_num = 1
                    shear_indices(1) = eqn_idx%stress%beg - 1 + 2
                    shear_bc_flip_num = 1
                    shear_bc_flip_indices(1:2,1) = shear_indices(1)
                    ! Both x-dir and y-dir: flip tau_xy only
                else if (num_dims == 3) then
                    shear_num = 3
                    shear_indices(1:3) = eqn_idx%stress%beg - 1 + (/2, 4, 5/)
                    shear_bc_flip_num = 2
                    shear_bc_flip_indices(1,1:2) = shear_indices((/1, 2/))
                    shear_bc_flip_indices(2,1:2) = shear_indices((/1, 3/))
                    shear_bc_flip_indices(3,1:2) = shear_indices((/2, 3/))
                    ! x-dir: flip tau_xy and tau_xz y-dir: flip tau_xy and tau_yz z-dir: flip tau_xz and tau_yz
                end if
            end if
 
            if (hyperelasticity) then
                eqn_idx%xi%beg = sys_size + 1
                eqn_idx%xi%end = sys_size + num_dims
                ! adding three more equations for the \xi field and the elastic energy
                sys_size = eqn_idx%xi%end + 1
                ! number of entries in the symmetric btensor plus the jacobian
                b_size = (num_dims*(num_dims + 1))/2 + 1
                tensor_size = num_dims**2 + 1
            end if
 
            if (surface_tension) then
                eqn_idx%c = sys_size + 1
                sys_size = eqn_idx%c
            end if
 
            if (cont_damage) then
                eqn_idx%damage = sys_size + 1
                sys_size = eqn_idx%damage
            else
                eqn_idx%damage = dflt_int
            end if
 
            if (hyper_cleaning) then
                eqn_idx%psi = sys_size + 1
                sys_size = eqn_idx%psi
            else
                eqn_idx%psi = dflt_int
            end if
        end if
 
        if (chemistry) then
            eqn_idx%species%beg = sys_size + 1
            eqn_idx%species%end = sys_size + num_species
            sys_size = eqn_idx%species%end
        else
            eqn_idx%species%beg = 1
            eqn_idx%species%end = 1
        end if
 
        if (output_partial_domain) then
            x_output_idx%beg = 0
            x_output_idx%end = 0
            y_output_idx%beg = 0
            y_output_idx%end = 0
            z_output_idx%beg = 0
            z_output_idx%end = 0
        end if
 
#ifdef MFC_MPI
        if (qbmm .and. .not. polytropic) then
            allocate (mpi_io_data%view(1:sys_size + 2*nb*nnode))
            allocate (mpi_io_data%var(1:sys_size + 2*nb*nnode))
        else
            allocate (mpi_io_data%view(1:sys_size))
            allocate (mpi_io_data%var(1:sys_size))
        end if
 
        do i = 1, sys_size
            if (down_sample) then
                allocate (mpi_io_data%var(i)%sf(-1:m + 1,-1:n + 1,-1:p + 1))
            else
                allocate (mpi_io_data%var(i)%sf(0:m,0:n,0:p))
            end if
            mpi_io_data%var(i)%sf => null()
        end do
        if (qbmm .and. .not. polytropic) then
            do i = sys_size + 1, sys_size + 2*nb*nnode
                allocate (mpi_io_data%var(i)%sf(0:m,0:n,0:p))
                mpi_io_data%var(i)%sf => null()
            end do
        end if
 
        if (ib) allocate (mpi_io_ib_data%var%sf(0:m,0:n,0:p))
#endif
 
        ! Size of the ghost zone layer is non-zero only when post-processing the raw simulation data of a parallel multidimensional
        ! computation in the Silo-HDF5 format. If this is the case, one must also verify whether the raw simulation data is 2D or
        ! 3D. In the 2D case, size of the z-coordinate direction ghost zone layer must be zeroed out.
        if (num_procs == 1 .or. format /= 1) then
            offset_x%beg = 0
            offset_x%end = 0
            offset_y%beg = 0
            offset_y%end = 0
            offset_z%beg = 0
            offset_z%end = 0
        else if (n == 0) then
            offset_y%beg = 0
            offset_y%end = 0
            offset_z%beg = 0
            offset_z%end = 0
        else if (p == 0) then
            offset_z%beg = 0
            offset_z%end = 0
        end if
 
        ! Determining the finite-difference number and the buffer size. Note that the size of the buffer is unrelated to the order
        ! of the WENO scheme. Rather, it is directly dependent on maximum size of ghost zone layers and possibly the order of the
        ! finite difference scheme used for the computation of vorticity and/or numerical Schlieren function.
        buff_size = max(offset_x%beg, offset_x%end, offset_y%beg, offset_y%end, offset_z%beg, offset_z%end)
 
        if (any(omega_wrt) .or. schlieren_wrt .or. qm_wrt .or. liutex_wrt) then
            fd_number = max(1, fd_order/2)
            buff_size = buff_size + fd_number
        end if
 
        ! Configuring Coordinate Direction Indexes
        idwint(1)%beg = 0; idwint(2)%beg = 0; idwint(3)%beg = 0
        idwint(1)%end = m; idwint(2)%end = n; idwint(3)%end = p
 
        idwbuff(1)%beg = -buff_size
        if (num_dims > 1) then; idwbuff(2)%beg = -buff_size; else; idwbuff(2)%beg = 0; end if
        if (num_dims > 2) then; idwbuff(3)%beg = -buff_size; else; idwbuff(3)%beg = 0; end if
 
        idwbuff(1)%end = idwint(1)%end - idwbuff(1)%beg
        idwbuff(2)%end = idwint(2)%end - idwbuff(2)%beg
        idwbuff(3)%end = idwint(3)%end - idwbuff(3)%beg
 
        ! Allocating single precision grid variables if needed
        allocate (x_cc_s(-buff_size:m + buff_size))
 
        ! Allocating the grid variables in the x-coordinate direction
        allocate (x_cb(-1 - offset_x%beg:m + offset_x%end))
        allocate (x_cc(-buff_size:m + buff_size))
        allocate (dx(-buff_size:m + buff_size))
 
        ! Allocating grid variables in the y- and z-coordinate directions
        if (n > 0) then
            allocate (y_cb(-1 - offset_y%beg:n + offset_y%end))
            allocate (y_cc(-buff_size:n + buff_size))
            allocate (dy(-buff_size:n + buff_size))
 
            if (p > 0) then
                allocate (z_cb(-1 - offset_z%beg:p + offset_z%end))
                allocate (z_cc(-buff_size:p + buff_size))
                allocate (dz(-buff_size:p + buff_size))
            end if
 
            ! Allocating the grid variables, only used for the 1D simulations, and containing the defragmented computational domain
            ! grid data
        else
            allocate (x_root_cb(-1:m_root))
            allocate (x_root_cc(0:m_root))
 
            if (precision == 1) then
                allocate (x_root_cc_s(0:m_root))
            end if
        end if
 
        allocate (adv(num_fluids))
 
        if (cyl_coord .neqv. .true.) then  ! Cartesian grid
            grid_geometry = 1
        else if (cyl_coord .and. p == 0) then  ! Axisymmetric cylindrical grid
            grid_geometry = 2
        else  ! Fully 3D cylindrical grid
            grid_geometry = 3
        end if
 
    end subroutine s_initialize_global_parameters_module
 
    !> Subroutine to initialize parallel infrastructure
    impure subroutine s_initialize_parallel_io
 
#ifdef MFC_MPI
        integer :: ierr  !< Generic flag used to identify and report MPI errors
#endif
 
        num_dims = 1 + min(1, n) + min(1, p)
 
        if (mhd) then
            num_vels = 3
        else
            num_vels = num_dims
        end if
 
        allocate (proc_coords(1:num_dims))
 
        if (parallel_io .neqv. .true.) return
 
#ifdef MFC_MPI
        ! 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
 
    !> Deallocation procedures for the module
    impure subroutine s_finalize_global_parameters_module
 
        integer :: i
 
        if (bubbles_euler) then
            deallocate (qbmm_idx%rs, qbmm_idx%vs, qbmm_idx%ps, qbmm_idx%ms)
            if (qbmm) deallocate (qbmm_idx%moms)
        end if
 
        ! Deallocating the grid variables for the x-coordinate direction
 
        deallocate (x_cc, x_cb, dx)
 
        ! Deallocating grid variables for the y- and z-coordinate directions
        if (n > 0) then
            deallocate (y_cc, y_cb, dy)
            if (p > 0) then
                deallocate (z_cc, z_cb, dz)
            end if
        else
            ! Deallocating the grid variables, only used for the 1D simulations, and containing the defragmented computational
            ! domain grid data
            deallocate (x_root_cb, x_root_cc)
        end if
 
        deallocate (proc_coords)
 
        deallocate (adv)
 
#ifdef MFC_MPI
        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
 
        if (ib) mpi_io_ib_data%var%sf => null()
#endif
 
    end subroutine s_finalize_global_parameters_module
 
end module m_global_parameters