m_riemann_solver_lf.fpp.f90

Go to the documentation of this file.

# 1 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!>
!! @file
!! @brief Contains module m_riemann_solver_lf
 
!> @brief Lax-Friedrichs (Rusanov) approximate Riemann solver
# 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.
# 8 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
 
! For moving immersed boundaries in simulation
# 12 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
# 7 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp" 2
# 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"
 
# 12 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 1 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp" 1
! Default for Fypp runs outside the toolchain (e.g. documentation). Searched last, so the
! build/include/<target>/thermochem.fpp the toolchain writes from m_thermochem takes precedence.
# 4 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
# 5 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
# 13 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 2
 
# 16 "/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"
 
# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 174 "/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"
 
# 12 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 1 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp" 1
! Default for Fypp runs outside the toolchain (e.g. documentation). Searched last, so the
! build/include/<target>/thermochem.fpp the toolchain writes from m_thermochem takes precedence.
# 4 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
# 5 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
# 13 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 2
 
# 16 "/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"
 
# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 174 "/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"
 
# 81 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 95 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 100 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 105 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 110 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 115 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 120 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 125 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 130 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 135 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 165 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 219 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 233 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 258 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 269 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 271 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 282 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 332 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 342 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 352 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 361 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 378 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 388 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 395 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 401 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 407 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 413 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 419 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 425 "/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"
 
# 12 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 1 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp" 1
! Default for Fypp runs outside the toolchain (e.g. documentation). Searched last, so the
! build/include/<target>/thermochem.fpp the toolchain writes from m_thermochem takes precedence.
# 4 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
# 5 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
# 13 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 2
 
# 16 "/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"
 
# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 174 "/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"
 
# 193 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 215 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 244 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 259 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 269 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 278 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 294 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 304 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 311 "/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
 
! GPU parallel region (scalar reductions, maxval/minval)
# 23 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! GPU parallel loop over threads (most common GPU macro)
# 43 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Required closing for GPU_PARALLEL_LOOP
# 55 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Mark routine for device compilation
# 112 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Declare device-resident data
# 130 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Inner loop within a GPU parallel region
# 145 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Scoped GPU data region
# 164 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Host code with device pointers (for MPI with GPU buffers)
# 193 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Allocate device memory (unscoped)
# 207 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Free device memory
# 219 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Atomic operation on device
# 231 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! End atomic capture block
# 242 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Copy data between host and device
# 254 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Synchronization barrier
# 266 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Import GPU library module (openacc or omp_lib)
# 275 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Emit code only for AMD compiler
# 282 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Emit code for non-Cray compilers
# 289 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Emit code only for Cray compiler
# 296 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Emit code for non-NVIDIA compilers
# 303 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 305 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
# 306 "/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.
# 52 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
! Allocate and create GPU device memory
# 72 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
! Free GPU device memory and deallocate
# 80 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
! Cray-specific GPU pointer setup for vector fields
# 104 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
! Cray-specific GPU pointer setup for scalar fields
# 120 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
! Cray-specific GPU pointer setup for acoustic source spatials
# 145 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 151 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 158 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
! New line at end of file is required for FYPP
# 8 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp" 2
 
module m_riemann_solver_lf
 
    use m_derived_types
    use m_global_parameters
    use m_variables_conversion
    use m_constants, only: riemann_solver_hll, riemann_solver_hllc, riemann_solver_lax_friedrichs
    use m_thermochem, only: gas_constant, get_mixture_molecular_weight, molecular_weights
    use m_thermochem_state, only: get_mixture_caloric_state
    use m_riemann_state
 
    implicit none
 
contains
 
    !> Lax-Friedrichs (Rusanov) approximate Riemann solver
    subroutine s_lf_riemann_solver(qL_prim_rsx_vf, dqL_prim_dx_vf, dqL_prim_dy_vf, dqL_prim_dz_vf, qL_prim_vf, qR_prim_rsx_vf, &
                                   & dqR_prim_dx_vf, dqR_prim_dy_vf, dqR_prim_dz_vf, qR_prim_vf, q_prim_vf, flux_vf, flux_src_vf, &
                                   & flux_gsrc_vf, norm_dir, ix, iy, iz)
 
        real(wp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:), intent(inout) :: qL_prim_rsx_vf, qR_prim_rsx_vf
        type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
        type(scalar_field), allocatable, dimension(:), intent(inout) :: qL_prim_vf, qR_prim_vf
        type(scalar_field), allocatable, dimension(:), intent(inout) :: dqL_prim_dx_vf, dqR_prim_dx_vf, dqL_prim_dy_vf, &
             & dqR_prim_dy_vf, dqL_prim_dz_vf, dqR_prim_dz_vf
 
        ! Intercell fluxes
        type(scalar_field), dimension(sys_size), intent(inout) :: flux_vf, flux_src_vf, flux_gsrc_vf
        integer, intent(in)                                    :: norm_dir
        type(int_bounds_info), intent(in)                      :: ix, iy, iz
 
# 45 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
            real(wp), dimension(num_fluids) :: alpha_rho_L, alpha_rho_R
            real(wp), dimension(num_vels)   :: vel_L, vel_R
            real(wp), dimension(num_fluids) :: alpha_L, alpha_R
            !> Averaged velocity gradient tensor `d(vel_i)/d(coord_j)`.
            real(wp), dimension(num_dims, num_dims) :: vel_grad_L, vel_grad_R
# 51 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
        real(wp), dimension(1) :: Ys_L, Ys_R
# 52 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
        real(wp), dimension(1) :: Cp_iL, Cp_iR, Xs_L, Xs_R, Gamma_iL, Gamma_iR
# 53 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
        real(wp) :: rho_L, rho_R
        real(wp) :: pres_L, pres_R
        real(wp) :: E_L, E_R
        real(wp) :: T_L, T_R
        real(wp) :: Y_L, Y_R
        real(wp) :: MW_L, MW_R
        real(wp) :: R_gas_L, R_gas_R
        real(wp) :: Cp_L, Cp_R
        real(wp) :: Cv_L, Cv_R
        real(wp) :: Gamm_L, Gamm_R
        real(wp) :: gamma_L, gamma_R
        real(wp) :: pi_inf_L, pi_inf_R
        real(wp) :: qv_L, qv_R
        real(wp) :: c_L, c_R
        real(wp), dimension(2) :: Re_L, Re_R
        real(wp) :: s_L, s_R, s_M, s_P
        real(wp) :: ptilde_L, ptilde_R
        real(wp) :: vel_L_rms, vel_R_rms
        real(wp) :: alpha_L_sum, alpha_R_sum
        real(wp) :: pcorr  !< low Mach number correction
        integer :: i, j, k, l  !< Generic loop iterators
        integer :: Re_size_loc1, Re_size_loc2  !< host copies of Re_size; amdflang reads the declare-target original stale cross-TU
        integer, dimension(3) :: idx_right_phys  !< Physical (j,k,l) indices for right state.
        ! Populating the buffers of the left and right Riemann problem states variables, based on the choice of boundary conditions
 
        call s_populate_riemann_states_variables_buffers(ql_prim_rsx_vf, dql_prim_dx_vf, dql_prim_dy_vf, dql_prim_dz_vf, &
            & qr_prim_rsx_vf, dqr_prim_dx_vf, dqr_prim_dy_vf, dqr_prim_dz_vf, norm_dir, ix, iy, iz)
 
        ! Reshaping inputted data based on dimensional splitting direction
        call s_initialize_riemann_solver(flux_src_vf, norm_dir)
        re_size_loc1 = re_size(1); re_size_loc2 = re_size(2)
# 88 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
# 89 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
# 90 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
            if (norm_dir == 1) then
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, l, alpha_rho_L, alpha_rho_R, vel_L, vel_R, alpha_L, alpha_R, Re_L, Re_R, s_L, s_R, Ys_L, Ys_R, Cp_iL, Cp_iR, Xs_L, Xs_R, &
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc& Gamma_iL, Gamma_iR, pcorr, vel_grad_L, vel_grad_R, idx_right_phys, vel_L_rms, vel_R_rms, alpha_L_sum, alpha_R_sum, pres_L, pres_R, rho_L, rho_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, qv_L, &
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc& qv_R, c_L, c_R, Gamm_L, Gamm_R, E_L, E_R, ptilde_L, ptilde_R, s_M, s_P, Cp_L, Cp_R, Cv_L, Cv_R, R_gas_L, R_gas_R, MW_L, MW_R, T_L, T_R, Y_L, Y_R) firstprivate(Re_size_loc1, Re_size_loc2)
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, j, k, l, &
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$omp& alpha_rho_L, alpha_rho_R, vel_L, vel_R, alpha_L, alpha_R, Re_L, Re_R, s_L, s_R, Ys_L, Ys_R, Cp_iL, Cp_iR, Xs_L, Xs_R, Gamma_iL, Gamma_iR, pcorr, vel_grad_L, vel_grad_R, idx_right_phys, &
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$omp& vel_L_rms, vel_R_rms, alpha_L_sum, alpha_R_sum, pres_L, pres_R, rho_L, rho_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, qv_L, qv_R, c_L, c_R, Gamm_L, Gamm_R, E_L, E_R, ptilde_L, ptilde_R, s_M, s_P, &
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$omp& Cp_L, Cp_R, Cv_L, Cv_R, R_gas_L, R_gas_R, MW_L, MW_R, T_L, T_R, Y_L, Y_R) firstprivate(Re_size_loc1, Re_size_loc2)
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 97 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                do l = is3%beg, is3%end
# 98 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                    do k = is2%beg, is2%end
# 99 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                        do j = is1%beg, is1%end
# 100 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 100 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 100 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 100 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 100 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 100 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 101 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            do i = 1, eqn_idx%cont%end
                                alpha_rho_l(i) = ql_prim_rsx_vf(j, k, l, i)
# 103 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                alpha_rho_r(i) = qr_prim_rsx_vf(j + 1, k, l, i)
# 104 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            end do
 
                            vel_l_rms = 0._wp; vel_r_rms = 0._wp
 
 
# 108 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 108 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 108 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 108 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 108 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 109 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            do i = 1, num_vels
                                vel_l(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%cont%end + i)
# 111 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                vel_r(i) = qr_prim_rsx_vf(j + 1, k, l, eqn_idx%cont%end + i)
# 112 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                vel_l_rms = vel_l_rms + vel_l(i)**2._wp
                                vel_r_rms = vel_r_rms + vel_r(i)**2._wp
                            end do
 
 
# 116 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 116 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 116 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 116 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 116 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 117 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            do i = 1, num_fluids
                                alpha_l(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%E + i)
# 119 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                alpha_r(i) = qr_prim_rsx_vf(j + 1, k, l, eqn_idx%E + i)
# 120 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            end do
 
                            pres_l = ql_prim_rsx_vf(j, k, l, eqn_idx%E)
# 123 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            pres_r = qr_prim_rsx_vf(j + 1, k, l, eqn_idx%E)
# 124 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
                            rho_l = 0._wp
                            gamma_l = 0._wp
                            pi_inf_l = 0._wp
                            qv_l = 0._wp
 
                            rho_r = 0._wp
                            gamma_r = 0._wp
                            pi_inf_r = 0._wp
                            qv_r = 0._wp
 
                            alpha_l_sum = 0._wp
                            alpha_r_sum = 0._wp
 
                            if (mpp_lim) then
 
# 139 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 139 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 139 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 139 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 139 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 140 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                do i = 1, num_fluids
                                    alpha_rho_l(i) = max(0._wp, alpha_rho_l(i))
                                    alpha_l(i) = min(max(0._wp, alpha_l(i)), 1._wp)
                                    alpha_l_sum = alpha_l_sum + alpha_l(i)
                                    alpha_rho_r(i) = max(0._wp, alpha_rho_r(i))
                                    alpha_r(i) = min(max(0._wp, alpha_r(i)), 1._wp)
                                    alpha_r_sum = alpha_r_sum + alpha_r(i)
                                end do
 
                                alpha_l = alpha_l/max(alpha_l_sum, sgm_eps)
                                alpha_r = alpha_r/max(alpha_r_sum, sgm_eps)
                            end if
 
                            call s_compute_mixture_coefficients(alpha_rho_l, alpha_l, rho_l, gamma_l, pi_inf_l, qv_l)
                            call s_compute_mixture_coefficients(alpha_rho_r, alpha_r, rho_r, gamma_r, pi_inf_r, qv_r)
 
                            if (viscous) then
                                call s_compute_interface_reynolds(alpha_l, re_l, re_size_loc1, re_size_loc2)
                                call s_compute_interface_reynolds(alpha_r, re_r, re_size_loc1, re_size_loc2)
                            end if
 
                            if (chemistry) then
 
# 162 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 162 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 162 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 162 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 162 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 163 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                do i = eqn_idx%species%beg, eqn_idx%species%end
                                    ys_l(i - eqn_idx%species%beg + 1) = ql_prim_rsx_vf(j, k, l, i)
# 165 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                    ys_r(i - eqn_idx%species%beg + 1) = qr_prim_rsx_vf(j + 1, k, l, i)
# 166 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                end do
 
                                call get_mixture_molecular_weight(ys_l, mw_l)
                                call get_mixture_molecular_weight(ys_r, mw_r)
 
                                xs_l(1:num_species) = ys_l(1:num_species)*mw_l/molecular_weights(:)
                                xs_r(1:num_species) = ys_r(1:num_species)*mw_r/molecular_weights(:)
 
                                r_gas_l = gas_constant/mw_l
                                r_gas_r = gas_constant/mw_r
                                t_l = pres_l/rho_l/r_gas_l
                                t_r = pres_r/rho_r/r_gas_r
 
                                call get_mixture_caloric_state(t_l, ys_l, cp_il, cp_l, cv_l, e_l)
                                call get_mixture_caloric_state(t_r, ys_r, cp_ir, cp_r, cv_r, e_r)
 
                                if (chem_params%gamma_method == 1) then
                                    ! gamma_method = 1: Ref. Section 2.3.1 Formulation of doi:10.7907/ZKW8-ES97.
                                    gamma_il(1:num_species) = cp_il(1:num_species)/(cp_il(1:num_species) - 1.0_wp)
                                    gamma_ir(1:num_species) = cp_ir(1:num_species)/(cp_ir(1:num_species) - 1.0_wp)
 
                                    gamma_l = sum(xs_l(1:num_species)/(gamma_il(1:num_species) - 1.0_wp))
                                    gamma_r = sum(xs_r(1:num_species)/(gamma_ir(1:num_species) - 1.0_wp))
                                else if (chem_params%gamma_method == 2) then
                                    ! gamma_method = 2: c_p / c_v where c_p, c_v are specific heats.
                                    gamm_l = cp_l/cv_l
                                    gamma_l = 1.0_wp/(gamm_l - 1.0_wp)
                                    gamm_r = cp_r/cv_r
                                    gamma_r = 1.0_wp/(gamm_r - 1.0_wp)
                                end if
 
                                e_l = rho_l*e_l + 5.e-1*rho_l*vel_l_rms
                                e_r = rho_r*e_r + 5.e-1*rho_r*vel_r_rms
                            else
                                call s_compute_energy(pres_l, alpha_rho_l, alpha_l, vel_l_rms, e_l)
                                call s_compute_energy(pres_r, alpha_rho_r, alpha_r, vel_r_rms, e_r)
                            end if
 
                            call s_compute_speed_of_sound(pres_l, rho_l, gamma_l, pi_inf_l, alpha_l, c_l, alpha_rho_l)
 
                            call s_compute_speed_of_sound(pres_r, rho_r, gamma_r, pi_inf_r, alpha_r, c_r, alpha_rho_r)
 
                            s_l = 0._wp; s_r = 0._wp
 
 
# 210 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 210 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 210 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 210 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 210 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 211 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            do i = 1, num_dims
                                s_l = s_l + vel_l(i)**2._wp
                                s_r = s_r + vel_r(i)**2._wp
                            end do
 
                            s_l = sqrt(s_l)
                            s_r = sqrt(s_r)
 
                            s_p = max(s_l, s_r) + max(c_l, c_r)
                            s_m = -s_p
 
                            s_l = s_m
                            s_r = s_p
 
                            ! Low Mach correction
                            pcorr = f_low_mach_pcorr_hll(vel_l_rms, vel_r_rms, c_l, c_r, rho_l, rho_r, s_m, s_p)
 
                            ! Mass
 
# 229 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 229 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 229 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 229 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 229 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 230 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            do i = 1, eqn_idx%cont%end
                                flux_rsx_vf(j, k, l, &
# 232 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                            & i) = (s_m*alpha_rho_r(i)*vel_r(norm_dir) - s_p*alpha_rho_l(i)*vel_l(norm_dir) &
                                            & + s_m*s_p*(alpha_rho_l(i) - alpha_rho_r(i)))/(s_m - s_p)
                            end do
 
                            ! Momentum
                            if (bubbles_euler) then
 
# 238 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 238 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 238 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 238 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 238 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 239 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                do i = 1, num_vels
                                    flux_rsx_vf(j, k, l, &
# 241 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                                & eqn_idx%cont%end + dir_idx(i)) = (s_m*(rho_r*vel_r(dir_idx(1))*vel_r(dir_idx(i)) &
                                                & + dir_flg(dir_idx(i))*(pres_r - ptilde_r)) - s_p*(rho_l*vel_l(dir_idx(1)) &
                                                & *vel_l(dir_idx(i)) + dir_flg(dir_idx(i))*(pres_l - ptilde_l)) &
                                                & + s_m*s_p*(rho_l*vel_l(dir_idx(i)) - rho_r*vel_r(dir_idx(i))))/(s_m - s_p) &
                                                & + (s_m/s_l)*(s_p/s_r)*pcorr*(vel_r(dir_idx(i)) - vel_l(dir_idx(i)))
                                end do
                            else
 
# 248 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 248 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 248 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 248 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 248 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 249 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                do i = 1, num_vels
                                    flux_rsx_vf(j, k, l, &
# 251 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                                & eqn_idx%cont%end + dir_idx(i)) = (s_m*(rho_r*vel_r(dir_idx(1))*vel_r(dir_idx(i)) &
                                                & + dir_flg(dir_idx(i))*pres_r) - s_p*(rho_l*vel_l(dir_idx(1))*vel_l(dir_idx(i)) &
                                                & + dir_flg(dir_idx(i))*pres_l) + s_m*s_p*(rho_l*vel_l(dir_idx(i)) &
                                                & - rho_r*vel_r(dir_idx(i))))/(s_m - s_p) + (s_m/s_l)*(s_p/s_r) &
                                                & *pcorr*(vel_r(dir_idx(i)) - vel_l(dir_idx(i)))
                                end do
                            end if
 
                            ! Energy
                            if (bubbles_euler) then
                                flux_rsx_vf(j, k, l, &
# 262 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                            & eqn_idx%E) = (s_m*vel_r(dir_idx(1))*(e_r + pres_r - ptilde_r) - s_p*vel_l(dir_idx(1) &
                                            & )*(e_l + pres_l - ptilde_l) + s_m*s_p*(e_l - e_r))/(s_m - s_p) + (s_m/s_l)*(s_p/s_r) &
                                            & *pcorr*(vel_r_rms - vel_l_rms)/2._wp
                            else
                                flux_rsx_vf(j, k, l, &
# 267 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                            & eqn_idx%E) = (s_m*vel_r(dir_idx(1))*(e_r + pres_r) - s_p*vel_l(dir_idx(1))*(e_l &
                                            & + pres_l) + s_m*s_p*(e_l - e_r))/(s_m - s_p) + (s_m/s_l)*(s_p/s_r)*pcorr*(vel_r_rms &
                                            & - vel_l_rms)/2._wp
                            end if
 
                            ! Advection flux and source: interface velocity for volume fraction transport
 
# 273 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 273 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 273 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 273 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 273 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 274 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            do i = eqn_idx%adv%beg, eqn_idx%adv%end
                                flux_rsx_vf(j, k, l, i) = (ql_prim_rsx_vf(j, k, l, i) - qr_prim_rsx_vf(j + 1, k, l, &
# 276 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                            & i))*s_m*s_p/(s_m - s_p)
                                flux_src_rsx_vf(j, k, l, i) = (s_m*qr_prim_rsx_vf(j + 1, k, l, &
# 278 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                                & i) - s_p*ql_prim_rsx_vf(j, k, l, i))/(s_m - s_p)
# 279 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            end do
 
                            if (bubbles_euler) then
                                ! From HLLC: Kills mass transport @ bubble gas density
                                if (num_fluids > 1) then
                                    flux_rsx_vf(j, k, l, eqn_idx%cont%end) = 0._wp
# 285 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                end if
                            end if
 
                            if (chemistry) then
 
# 289 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 289 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 289 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 289 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 289 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 290 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                do i = eqn_idx%species%beg, eqn_idx%species%end
                                    y_l = ql_prim_rsx_vf(j, k, l, i)
# 292 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                    y_r = qr_prim_rsx_vf(j + 1, k, l, i)
# 293 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
                                    flux_rsx_vf(j, k, l, &
# 295 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                                & i) = (s_m*y_r*rho_r*vel_r(dir_idx(1)) - s_p*y_l*rho_l*vel_l(dir_idx(1)) &
                                                & + s_m*s_p*(y_l*rho_l - y_r*rho_r))/(s_m - s_p)
                                    flux_src_rsx_vf(j, k, l, i) = 0._wp
# 298 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                end do
                            end if
 
# 318 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                        end do
                    end do
                end do
 
# 321 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 321 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc end parallel loop
# 321 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 321 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 321 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$omp end target teams loop
# 321 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 322 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
            end if
# 88 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
# 89 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
# 90 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
            if (norm_dir == 2) then
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, l, alpha_rho_L, alpha_rho_R, vel_L, vel_R, alpha_L, alpha_R, Re_L, Re_R, s_L, s_R, Ys_L, Ys_R, Cp_iL, Cp_iR, Xs_L, Xs_R, &
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc& Gamma_iL, Gamma_iR, pcorr, vel_grad_L, vel_grad_R, idx_right_phys, vel_L_rms, vel_R_rms, alpha_L_sum, alpha_R_sum, pres_L, pres_R, rho_L, rho_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, qv_L, &
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc& qv_R, c_L, c_R, Gamm_L, Gamm_R, E_L, E_R, ptilde_L, ptilde_R, s_M, s_P, Cp_L, Cp_R, Cv_L, Cv_R, R_gas_L, R_gas_R, MW_L, MW_R, T_L, T_R, Y_L, Y_R) firstprivate(Re_size_loc1, Re_size_loc2)
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, j, k, l, &
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$omp& alpha_rho_L, alpha_rho_R, vel_L, vel_R, alpha_L, alpha_R, Re_L, Re_R, s_L, s_R, Ys_L, Ys_R, Cp_iL, Cp_iR, Xs_L, Xs_R, Gamma_iL, Gamma_iR, pcorr, vel_grad_L, vel_grad_R, idx_right_phys, &
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$omp& vel_L_rms, vel_R_rms, alpha_L_sum, alpha_R_sum, pres_L, pres_R, rho_L, rho_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, qv_L, qv_R, c_L, c_R, Gamm_L, Gamm_R, E_L, E_R, ptilde_L, ptilde_R, s_M, s_P, &
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$omp& Cp_L, Cp_R, Cv_L, Cv_R, R_gas_L, R_gas_R, MW_L, MW_R, T_L, T_R, Y_L, Y_R) firstprivate(Re_size_loc1, Re_size_loc2)
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 97 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                do l = is3%beg, is3%end
# 98 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                    do k = is1%beg, is1%end
# 99 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                        do j = is2%beg, is2%end
# 100 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 100 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 100 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 100 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 100 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 100 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 101 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            do i = 1, eqn_idx%cont%end
                                alpha_rho_l(i) = ql_prim_rsx_vf(j, k, l, i)
# 103 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                alpha_rho_r(i) = qr_prim_rsx_vf(j, k + 1, l, i)
# 104 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            end do
 
                            vel_l_rms = 0._wp; vel_r_rms = 0._wp
 
 
# 108 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 108 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 108 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 108 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 108 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 109 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            do i = 1, num_vels
                                vel_l(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%cont%end + i)
# 111 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                vel_r(i) = qr_prim_rsx_vf(j, k + 1, l, eqn_idx%cont%end + i)
# 112 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                vel_l_rms = vel_l_rms + vel_l(i)**2._wp
                                vel_r_rms = vel_r_rms + vel_r(i)**2._wp
                            end do
 
 
# 116 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 116 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 116 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 116 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 116 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 117 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            do i = 1, num_fluids
                                alpha_l(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%E + i)
# 119 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                alpha_r(i) = qr_prim_rsx_vf(j, k + 1, l, eqn_idx%E + i)
# 120 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            end do
 
                            pres_l = ql_prim_rsx_vf(j, k, l, eqn_idx%E)
# 123 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            pres_r = qr_prim_rsx_vf(j, k + 1, l, eqn_idx%E)
# 124 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
                            rho_l = 0._wp
                            gamma_l = 0._wp
                            pi_inf_l = 0._wp
                            qv_l = 0._wp
 
                            rho_r = 0._wp
                            gamma_r = 0._wp
                            pi_inf_r = 0._wp
                            qv_r = 0._wp
 
                            alpha_l_sum = 0._wp
                            alpha_r_sum = 0._wp
 
                            if (mpp_lim) then
 
# 139 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 139 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 139 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 139 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 139 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 140 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                do i = 1, num_fluids
                                    alpha_rho_l(i) = max(0._wp, alpha_rho_l(i))
                                    alpha_l(i) = min(max(0._wp, alpha_l(i)), 1._wp)
                                    alpha_l_sum = alpha_l_sum + alpha_l(i)
                                    alpha_rho_r(i) = max(0._wp, alpha_rho_r(i))
                                    alpha_r(i) = min(max(0._wp, alpha_r(i)), 1._wp)
                                    alpha_r_sum = alpha_r_sum + alpha_r(i)
                                end do
 
                                alpha_l = alpha_l/max(alpha_l_sum, sgm_eps)
                                alpha_r = alpha_r/max(alpha_r_sum, sgm_eps)
                            end if
 
                            call s_compute_mixture_coefficients(alpha_rho_l, alpha_l, rho_l, gamma_l, pi_inf_l, qv_l)
                            call s_compute_mixture_coefficients(alpha_rho_r, alpha_r, rho_r, gamma_r, pi_inf_r, qv_r)
 
                            if (viscous) then
                                call s_compute_interface_reynolds(alpha_l, re_l, re_size_loc1, re_size_loc2)
                                call s_compute_interface_reynolds(alpha_r, re_r, re_size_loc1, re_size_loc2)
                            end if
 
                            if (chemistry) then
 
# 162 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 162 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 162 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 162 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 162 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 163 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                do i = eqn_idx%species%beg, eqn_idx%species%end
                                    ys_l(i - eqn_idx%species%beg + 1) = ql_prim_rsx_vf(j, k, l, i)
# 165 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                    ys_r(i - eqn_idx%species%beg + 1) = qr_prim_rsx_vf(j, k + 1, l, i)
# 166 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                end do
 
                                call get_mixture_molecular_weight(ys_l, mw_l)
                                call get_mixture_molecular_weight(ys_r, mw_r)
 
                                xs_l(1:num_species) = ys_l(1:num_species)*mw_l/molecular_weights(:)
                                xs_r(1:num_species) = ys_r(1:num_species)*mw_r/molecular_weights(:)
 
                                r_gas_l = gas_constant/mw_l
                                r_gas_r = gas_constant/mw_r
                                t_l = pres_l/rho_l/r_gas_l
                                t_r = pres_r/rho_r/r_gas_r
 
                                call get_mixture_caloric_state(t_l, ys_l, cp_il, cp_l, cv_l, e_l)
                                call get_mixture_caloric_state(t_r, ys_r, cp_ir, cp_r, cv_r, e_r)
 
                                if (chem_params%gamma_method == 1) then
                                    ! gamma_method = 1: Ref. Section 2.3.1 Formulation of doi:10.7907/ZKW8-ES97.
                                    gamma_il(1:num_species) = cp_il(1:num_species)/(cp_il(1:num_species) - 1.0_wp)
                                    gamma_ir(1:num_species) = cp_ir(1:num_species)/(cp_ir(1:num_species) - 1.0_wp)
 
                                    gamma_l = sum(xs_l(1:num_species)/(gamma_il(1:num_species) - 1.0_wp))
                                    gamma_r = sum(xs_r(1:num_species)/(gamma_ir(1:num_species) - 1.0_wp))
                                else if (chem_params%gamma_method == 2) then
                                    ! gamma_method = 2: c_p / c_v where c_p, c_v are specific heats.
                                    gamm_l = cp_l/cv_l
                                    gamma_l = 1.0_wp/(gamm_l - 1.0_wp)
                                    gamm_r = cp_r/cv_r
                                    gamma_r = 1.0_wp/(gamm_r - 1.0_wp)
                                end if
 
                                e_l = rho_l*e_l + 5.e-1*rho_l*vel_l_rms
                                e_r = rho_r*e_r + 5.e-1*rho_r*vel_r_rms
                            else
                                call s_compute_energy(pres_l, alpha_rho_l, alpha_l, vel_l_rms, e_l)
                                call s_compute_energy(pres_r, alpha_rho_r, alpha_r, vel_r_rms, e_r)
                            end if
 
                            call s_compute_speed_of_sound(pres_l, rho_l, gamma_l, pi_inf_l, alpha_l, c_l, alpha_rho_l)
 
                            call s_compute_speed_of_sound(pres_r, rho_r, gamma_r, pi_inf_r, alpha_r, c_r, alpha_rho_r)
 
                            s_l = 0._wp; s_r = 0._wp
 
 
# 210 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 210 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 210 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 210 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 210 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 211 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            do i = 1, num_dims
                                s_l = s_l + vel_l(i)**2._wp
                                s_r = s_r + vel_r(i)**2._wp
                            end do
 
                            s_l = sqrt(s_l)
                            s_r = sqrt(s_r)
 
                            s_p = max(s_l, s_r) + max(c_l, c_r)
                            s_m = -s_p
 
                            s_l = s_m
                            s_r = s_p
 
                            ! Low Mach correction
                            pcorr = f_low_mach_pcorr_hll(vel_l_rms, vel_r_rms, c_l, c_r, rho_l, rho_r, s_m, s_p)
 
                            ! Mass
 
# 229 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 229 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 229 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 229 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 229 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 230 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            do i = 1, eqn_idx%cont%end
                                flux_rsx_vf(j, k, l, &
# 232 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                            & i) = (s_m*alpha_rho_r(i)*vel_r(norm_dir) - s_p*alpha_rho_l(i)*vel_l(norm_dir) &
                                            & + s_m*s_p*(alpha_rho_l(i) - alpha_rho_r(i)))/(s_m - s_p)
                            end do
 
                            ! Momentum
                            if (bubbles_euler) then
 
# 238 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 238 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 238 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 238 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 238 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 239 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                do i = 1, num_vels
                                    flux_rsx_vf(j, k, l, &
# 241 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                                & eqn_idx%cont%end + dir_idx(i)) = (s_m*(rho_r*vel_r(dir_idx(1))*vel_r(dir_idx(i)) &
                                                & + dir_flg(dir_idx(i))*(pres_r - ptilde_r)) - s_p*(rho_l*vel_l(dir_idx(1)) &
                                                & *vel_l(dir_idx(i)) + dir_flg(dir_idx(i))*(pres_l - ptilde_l)) &
                                                & + s_m*s_p*(rho_l*vel_l(dir_idx(i)) - rho_r*vel_r(dir_idx(i))))/(s_m - s_p) &
                                                & + (s_m/s_l)*(s_p/s_r)*pcorr*(vel_r(dir_idx(i)) - vel_l(dir_idx(i)))
                                end do
                            else
 
# 248 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 248 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 248 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 248 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 248 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 249 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                do i = 1, num_vels
                                    flux_rsx_vf(j, k, l, &
# 251 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                                & eqn_idx%cont%end + dir_idx(i)) = (s_m*(rho_r*vel_r(dir_idx(1))*vel_r(dir_idx(i)) &
                                                & + dir_flg(dir_idx(i))*pres_r) - s_p*(rho_l*vel_l(dir_idx(1))*vel_l(dir_idx(i)) &
                                                & + dir_flg(dir_idx(i))*pres_l) + s_m*s_p*(rho_l*vel_l(dir_idx(i)) &
                                                & - rho_r*vel_r(dir_idx(i))))/(s_m - s_p) + (s_m/s_l)*(s_p/s_r) &
                                                & *pcorr*(vel_r(dir_idx(i)) - vel_l(dir_idx(i)))
                                end do
                            end if
 
                            ! Energy
                            if (bubbles_euler) then
                                flux_rsx_vf(j, k, l, &
# 262 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                            & eqn_idx%E) = (s_m*vel_r(dir_idx(1))*(e_r + pres_r - ptilde_r) - s_p*vel_l(dir_idx(1) &
                                            & )*(e_l + pres_l - ptilde_l) + s_m*s_p*(e_l - e_r))/(s_m - s_p) + (s_m/s_l)*(s_p/s_r) &
                                            & *pcorr*(vel_r_rms - vel_l_rms)/2._wp
                            else
                                flux_rsx_vf(j, k, l, &
# 267 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                            & eqn_idx%E) = (s_m*vel_r(dir_idx(1))*(e_r + pres_r) - s_p*vel_l(dir_idx(1))*(e_l &
                                            & + pres_l) + s_m*s_p*(e_l - e_r))/(s_m - s_p) + (s_m/s_l)*(s_p/s_r)*pcorr*(vel_r_rms &
                                            & - vel_l_rms)/2._wp
                            end if
 
                            ! Advection flux and source: interface velocity for volume fraction transport
 
# 273 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 273 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 273 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 273 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 273 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 274 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            do i = eqn_idx%adv%beg, eqn_idx%adv%end
                                flux_rsx_vf(j, k, l, i) = (ql_prim_rsx_vf(j, k, l, i) - qr_prim_rsx_vf(j, k + 1, l, &
# 276 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                            & i))*s_m*s_p/(s_m - s_p)
                                flux_src_rsx_vf(j, k, l, i) = (s_m*qr_prim_rsx_vf(j, k + 1, l, &
# 278 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                                & i) - s_p*ql_prim_rsx_vf(j, k, l, i))/(s_m - s_p)
# 279 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            end do
 
                            if (bubbles_euler) then
                                ! From HLLC: Kills mass transport @ bubble gas density
                                if (num_fluids > 1) then
                                    flux_rsx_vf(j, k, l, eqn_idx%cont%end) = 0._wp
# 285 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                end if
                            end if
 
                            if (chemistry) then
 
# 289 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 289 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 289 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 289 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 289 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 290 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                do i = eqn_idx%species%beg, eqn_idx%species%end
                                    y_l = ql_prim_rsx_vf(j, k, l, i)
# 292 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                    y_r = qr_prim_rsx_vf(j, k + 1, l, i)
# 293 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
                                    flux_rsx_vf(j, k, l, &
# 295 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                                & i) = (s_m*y_r*rho_r*vel_r(dir_idx(1)) - s_p*y_l*rho_l*vel_l(dir_idx(1)) &
                                                & + s_m*s_p*(y_l*rho_l - y_r*rho_r))/(s_m - s_p)
                                    flux_src_rsx_vf(j, k, l, i) = 0._wp
# 298 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                end do
                            end if
 
# 302 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                if (cyl_coord) then
                                    ! Substituting the advective flux into the inviscid geometrical source flux
 
# 304 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 304 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 304 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 304 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 304 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 305 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                    do i = 1, eqn_idx%E
                                        flux_gsrc_rsx_vf(j, k, l, i) = flux_rsx_vf(j, k, l, i)
# 307 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                    end do
                                    ! Recalculating the radial momentum geometric source flux
                                    flux_gsrc_rsx_vf(j, k, l, eqn_idx%cont%end + 2) = flux_rsx_vf(j, k, l, &
# 310 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                                     & eqn_idx%cont%end + 2) - (s_m*pres_r - s_p*pres_l)/(s_m - s_p)
                                    ! Geometrical source of the void fraction(s) is zero
 
# 312 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 312 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 312 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 312 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 312 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 313 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                    do i = eqn_idx%adv%beg, eqn_idx%adv%end
                                        flux_gsrc_rsx_vf(j, k, l, i) = flux_rsx_vf(j, k, l, i)
# 315 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                    end do
                                end if
# 318 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                        end do
                    end do
                end do
 
# 321 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 321 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc end parallel loop
# 321 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 321 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 321 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$omp end target teams loop
# 321 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 322 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
            end if
# 88 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
# 89 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
# 90 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
            if (norm_dir == 3) then
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, l, alpha_rho_L, alpha_rho_R, vel_L, vel_R, alpha_L, alpha_R, Re_L, Re_R, s_L, s_R, Ys_L, Ys_R, Cp_iL, Cp_iR, Xs_L, Xs_R, &
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc& Gamma_iL, Gamma_iR, pcorr, vel_grad_L, vel_grad_R, idx_right_phys, vel_L_rms, vel_R_rms, alpha_L_sum, alpha_R_sum, pres_L, pres_R, rho_L, rho_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, qv_L, &
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc& qv_R, c_L, c_R, Gamm_L, Gamm_R, E_L, E_R, ptilde_L, ptilde_R, s_M, s_P, Cp_L, Cp_R, Cv_L, Cv_R, R_gas_L, R_gas_R, MW_L, MW_R, T_L, T_R, Y_L, Y_R) firstprivate(Re_size_loc1, Re_size_loc2)
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, j, k, l, &
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$omp& alpha_rho_L, alpha_rho_R, vel_L, vel_R, alpha_L, alpha_R, Re_L, Re_R, s_L, s_R, Ys_L, Ys_R, Cp_iL, Cp_iR, Xs_L, Xs_R, Gamma_iL, Gamma_iR, pcorr, vel_grad_L, vel_grad_R, idx_right_phys, &
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$omp& vel_L_rms, vel_R_rms, alpha_L_sum, alpha_R_sum, pres_L, pres_R, rho_L, rho_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, qv_L, qv_R, c_L, c_R, Gamm_L, Gamm_R, E_L, E_R, ptilde_L, ptilde_R, s_M, s_P, &
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$omp& Cp_L, Cp_R, Cv_L, Cv_R, R_gas_L, R_gas_R, MW_L, MW_R, T_L, T_R, Y_L, Y_R) firstprivate(Re_size_loc1, Re_size_loc2)
# 91 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 97 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                do l = is1%beg, is1%end
# 98 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                    do k = is2%beg, is2%end
# 99 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                        do j = is3%beg, is3%end
# 100 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 100 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 100 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 100 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 100 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 100 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 101 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            do i = 1, eqn_idx%cont%end
                                alpha_rho_l(i) = ql_prim_rsx_vf(j, k, l, i)
# 103 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                alpha_rho_r(i) = qr_prim_rsx_vf(j, k, l + 1, i)
# 104 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            end do
 
                            vel_l_rms = 0._wp; vel_r_rms = 0._wp
 
 
# 108 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 108 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 108 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 108 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 108 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 109 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            do i = 1, num_vels
                                vel_l(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%cont%end + i)
# 111 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                vel_r(i) = qr_prim_rsx_vf(j, k, l + 1, eqn_idx%cont%end + i)
# 112 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                vel_l_rms = vel_l_rms + vel_l(i)**2._wp
                                vel_r_rms = vel_r_rms + vel_r(i)**2._wp
                            end do
 
 
# 116 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 116 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 116 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 116 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 116 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 117 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            do i = 1, num_fluids
                                alpha_l(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%E + i)
# 119 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                alpha_r(i) = qr_prim_rsx_vf(j, k, l + 1, eqn_idx%E + i)
# 120 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            end do
 
                            pres_l = ql_prim_rsx_vf(j, k, l, eqn_idx%E)
# 123 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            pres_r = qr_prim_rsx_vf(j, k, l + 1, eqn_idx%E)
# 124 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
                            rho_l = 0._wp
                            gamma_l = 0._wp
                            pi_inf_l = 0._wp
                            qv_l = 0._wp
 
                            rho_r = 0._wp
                            gamma_r = 0._wp
                            pi_inf_r = 0._wp
                            qv_r = 0._wp
 
                            alpha_l_sum = 0._wp
                            alpha_r_sum = 0._wp
 
                            if (mpp_lim) then
 
# 139 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 139 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 139 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 139 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 139 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 140 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                do i = 1, num_fluids
                                    alpha_rho_l(i) = max(0._wp, alpha_rho_l(i))
                                    alpha_l(i) = min(max(0._wp, alpha_l(i)), 1._wp)
                                    alpha_l_sum = alpha_l_sum + alpha_l(i)
                                    alpha_rho_r(i) = max(0._wp, alpha_rho_r(i))
                                    alpha_r(i) = min(max(0._wp, alpha_r(i)), 1._wp)
                                    alpha_r_sum = alpha_r_sum + alpha_r(i)
                                end do
 
                                alpha_l = alpha_l/max(alpha_l_sum, sgm_eps)
                                alpha_r = alpha_r/max(alpha_r_sum, sgm_eps)
                            end if
 
                            call s_compute_mixture_coefficients(alpha_rho_l, alpha_l, rho_l, gamma_l, pi_inf_l, qv_l)
                            call s_compute_mixture_coefficients(alpha_rho_r, alpha_r, rho_r, gamma_r, pi_inf_r, qv_r)
 
                            if (viscous) then
                                call s_compute_interface_reynolds(alpha_l, re_l, re_size_loc1, re_size_loc2)
                                call s_compute_interface_reynolds(alpha_r, re_r, re_size_loc1, re_size_loc2)
                            end if
 
                            if (chemistry) then
 
# 162 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 162 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 162 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 162 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 162 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 163 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                do i = eqn_idx%species%beg, eqn_idx%species%end
                                    ys_l(i - eqn_idx%species%beg + 1) = ql_prim_rsx_vf(j, k, l, i)
# 165 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                    ys_r(i - eqn_idx%species%beg + 1) = qr_prim_rsx_vf(j, k, l + 1, i)
# 166 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                end do
 
                                call get_mixture_molecular_weight(ys_l, mw_l)
                                call get_mixture_molecular_weight(ys_r, mw_r)
 
                                xs_l(1:num_species) = ys_l(1:num_species)*mw_l/molecular_weights(:)
                                xs_r(1:num_species) = ys_r(1:num_species)*mw_r/molecular_weights(:)
 
                                r_gas_l = gas_constant/mw_l
                                r_gas_r = gas_constant/mw_r
                                t_l = pres_l/rho_l/r_gas_l
                                t_r = pres_r/rho_r/r_gas_r
 
                                call get_mixture_caloric_state(t_l, ys_l, cp_il, cp_l, cv_l, e_l)
                                call get_mixture_caloric_state(t_r, ys_r, cp_ir, cp_r, cv_r, e_r)
 
                                if (chem_params%gamma_method == 1) then
                                    ! gamma_method = 1: Ref. Section 2.3.1 Formulation of doi:10.7907/ZKW8-ES97.
                                    gamma_il(1:num_species) = cp_il(1:num_species)/(cp_il(1:num_species) - 1.0_wp)
                                    gamma_ir(1:num_species) = cp_ir(1:num_species)/(cp_ir(1:num_species) - 1.0_wp)
 
                                    gamma_l = sum(xs_l(1:num_species)/(gamma_il(1:num_species) - 1.0_wp))
                                    gamma_r = sum(xs_r(1:num_species)/(gamma_ir(1:num_species) - 1.0_wp))
                                else if (chem_params%gamma_method == 2) then
                                    ! gamma_method = 2: c_p / c_v where c_p, c_v are specific heats.
                                    gamm_l = cp_l/cv_l
                                    gamma_l = 1.0_wp/(gamm_l - 1.0_wp)
                                    gamm_r = cp_r/cv_r
                                    gamma_r = 1.0_wp/(gamm_r - 1.0_wp)
                                end if
 
                                e_l = rho_l*e_l + 5.e-1*rho_l*vel_l_rms
                                e_r = rho_r*e_r + 5.e-1*rho_r*vel_r_rms
                            else
                                call s_compute_energy(pres_l, alpha_rho_l, alpha_l, vel_l_rms, e_l)
                                call s_compute_energy(pres_r, alpha_rho_r, alpha_r, vel_r_rms, e_r)
                            end if
 
                            call s_compute_speed_of_sound(pres_l, rho_l, gamma_l, pi_inf_l, alpha_l, c_l, alpha_rho_l)
 
                            call s_compute_speed_of_sound(pres_r, rho_r, gamma_r, pi_inf_r, alpha_r, c_r, alpha_rho_r)
 
                            s_l = 0._wp; s_r = 0._wp
 
 
# 210 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 210 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 210 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 210 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 210 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 211 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            do i = 1, num_dims
                                s_l = s_l + vel_l(i)**2._wp
                                s_r = s_r + vel_r(i)**2._wp
                            end do
 
                            s_l = sqrt(s_l)
                            s_r = sqrt(s_r)
 
                            s_p = max(s_l, s_r) + max(c_l, c_r)
                            s_m = -s_p
 
                            s_l = s_m
                            s_r = s_p
 
                            ! Low Mach correction
                            pcorr = f_low_mach_pcorr_hll(vel_l_rms, vel_r_rms, c_l, c_r, rho_l, rho_r, s_m, s_p)
 
                            ! Mass
 
# 229 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 229 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 229 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 229 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 229 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 230 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            do i = 1, eqn_idx%cont%end
                                flux_rsx_vf(j, k, l, &
# 232 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                            & i) = (s_m*alpha_rho_r(i)*vel_r(norm_dir) - s_p*alpha_rho_l(i)*vel_l(norm_dir) &
                                            & + s_m*s_p*(alpha_rho_l(i) - alpha_rho_r(i)))/(s_m - s_p)
                            end do
 
                            ! Momentum
                            if (bubbles_euler) then
 
# 238 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 238 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 238 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 238 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 238 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 239 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                do i = 1, num_vels
                                    flux_rsx_vf(j, k, l, &
# 241 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                                & eqn_idx%cont%end + dir_idx(i)) = (s_m*(rho_r*vel_r(dir_idx(1))*vel_r(dir_idx(i)) &
                                                & + dir_flg(dir_idx(i))*(pres_r - ptilde_r)) - s_p*(rho_l*vel_l(dir_idx(1)) &
                                                & *vel_l(dir_idx(i)) + dir_flg(dir_idx(i))*(pres_l - ptilde_l)) &
                                                & + s_m*s_p*(rho_l*vel_l(dir_idx(i)) - rho_r*vel_r(dir_idx(i))))/(s_m - s_p) &
                                                & + (s_m/s_l)*(s_p/s_r)*pcorr*(vel_r(dir_idx(i)) - vel_l(dir_idx(i)))
                                end do
                            else
 
# 248 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 248 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 248 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 248 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 248 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 249 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                do i = 1, num_vels
                                    flux_rsx_vf(j, k, l, &
# 251 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                                & eqn_idx%cont%end + dir_idx(i)) = (s_m*(rho_r*vel_r(dir_idx(1))*vel_r(dir_idx(i)) &
                                                & + dir_flg(dir_idx(i))*pres_r) - s_p*(rho_l*vel_l(dir_idx(1))*vel_l(dir_idx(i)) &
                                                & + dir_flg(dir_idx(i))*pres_l) + s_m*s_p*(rho_l*vel_l(dir_idx(i)) &
                                                & - rho_r*vel_r(dir_idx(i))))/(s_m - s_p) + (s_m/s_l)*(s_p/s_r) &
                                                & *pcorr*(vel_r(dir_idx(i)) - vel_l(dir_idx(i)))
                                end do
                            end if
 
                            ! Energy
                            if (bubbles_euler) then
                                flux_rsx_vf(j, k, l, &
# 262 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                            & eqn_idx%E) = (s_m*vel_r(dir_idx(1))*(e_r + pres_r - ptilde_r) - s_p*vel_l(dir_idx(1) &
                                            & )*(e_l + pres_l - ptilde_l) + s_m*s_p*(e_l - e_r))/(s_m - s_p) + (s_m/s_l)*(s_p/s_r) &
                                            & *pcorr*(vel_r_rms - vel_l_rms)/2._wp
                            else
                                flux_rsx_vf(j, k, l, &
# 267 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                            & eqn_idx%E) = (s_m*vel_r(dir_idx(1))*(e_r + pres_r) - s_p*vel_l(dir_idx(1))*(e_l &
                                            & + pres_l) + s_m*s_p*(e_l - e_r))/(s_m - s_p) + (s_m/s_l)*(s_p/s_r)*pcorr*(vel_r_rms &
                                            & - vel_l_rms)/2._wp
                            end if
 
                            ! Advection flux and source: interface velocity for volume fraction transport
 
# 273 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 273 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 273 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 273 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 273 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 274 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            do i = eqn_idx%adv%beg, eqn_idx%adv%end
                                flux_rsx_vf(j, k, l, i) = (ql_prim_rsx_vf(j, k, l, i) - qr_prim_rsx_vf(j, k, l + 1, &
# 276 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                            & i))*s_m*s_p/(s_m - s_p)
                                flux_src_rsx_vf(j, k, l, i) = (s_m*qr_prim_rsx_vf(j, k, l + 1, &
# 278 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                                & i) - s_p*ql_prim_rsx_vf(j, k, l, i))/(s_m - s_p)
# 279 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            end do
 
                            if (bubbles_euler) then
                                ! From HLLC: Kills mass transport @ bubble gas density
                                if (num_fluids > 1) then
                                    flux_rsx_vf(j, k, l, eqn_idx%cont%end) = 0._wp
# 285 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                end if
                            end if
 
                            if (chemistry) then
 
# 289 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 289 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 289 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 289 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 289 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 290 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                do i = eqn_idx%species%beg, eqn_idx%species%end
                                    y_l = ql_prim_rsx_vf(j, k, l, i)
# 292 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                    y_r = qr_prim_rsx_vf(j, k, l + 1, i)
# 293 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
                                    flux_rsx_vf(j, k, l, &
# 295 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                                & i) = (s_m*y_r*rho_r*vel_r(dir_idx(1)) - s_p*y_l*rho_l*vel_l(dir_idx(1)) &
                                                & + s_m*s_p*(y_l*rho_l - y_r*rho_r))/(s_m - s_p)
                                    flux_src_rsx_vf(j, k, l, i) = 0._wp
# 298 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                end do
                            end if
 
# 318 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                        end do
                    end do
                end do
 
# 321 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 321 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc end parallel loop
# 321 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 321 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 321 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$omp end target teams loop
# 321 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 322 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
            end if
# 324 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
        if (viscous) then
 
# 326 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 326 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 326 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, l, idx_right_phys, vel_grad_L, vel_grad_R, alpha_L, alpha_R, vel_L, vel_R, Re_L, Re_R) &
# 326 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc& firstprivate(Re_size_loc1, Re_size_loc2) copyin(norm_dir)
# 326 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 326 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 326 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 326 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 326 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
# 326 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$omp& private(i, j, k, l, idx_right_phys, vel_grad_L, vel_grad_R, alpha_L, alpha_R, vel_L, vel_R, Re_L, Re_R) firstprivate(Re_size_loc1, Re_size_loc2) map(to:norm_dir)
# 326 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 328 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
            do l = isz%beg, isz%end
                do k = isy%beg, isy%end
                    do j = isx%beg, isx%end
                        idx_right_phys(1) = j
                        idx_right_phys(2) = k
                        idx_right_phys(3) = l
                        idx_right_phys(norm_dir) = idx_right_phys(norm_dir) + 1
 
                        if (norm_dir == 1) then
 
# 337 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 337 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 337 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 337 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 337 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 338 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            do i = 1, num_fluids
                                alpha_l(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%E + i)
                                alpha_r(i) = qr_prim_rsx_vf(j + 1, k, l, eqn_idx%E + i)
                            end do
 
 
# 343 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 343 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 343 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 343 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 343 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 344 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            do i = 1, num_dims
                                vel_l(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%mom%beg + i - 1)
                                vel_r(i) = qr_prim_rsx_vf(j + 1, k, l, eqn_idx%mom%beg + i - 1)
                            end do
                        else if (norm_dir == 2) then
 
# 349 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 349 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 349 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 349 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 349 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 350 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            do i = 1, num_fluids
                                alpha_l(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%E + i)
                                alpha_r(i) = qr_prim_rsx_vf(j, k + 1, l, eqn_idx%E + i)
                            end do
 
# 354 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 354 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 354 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 354 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 354 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 355 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            do i = 1, num_dims
                                vel_l(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%mom%beg + i - 1)
                                vel_r(i) = qr_prim_rsx_vf(j, k + 1, l, eqn_idx%mom%beg + i - 1)
                            end do
                        else
 
# 360 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 360 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 360 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 360 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 360 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 361 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            do i = 1, num_fluids
                                alpha_l(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%E + i)
                                alpha_r(i) = qr_prim_rsx_vf(j, k, l + 1, eqn_idx%E + i)
                            end do
 
 
# 366 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 366 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 366 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 366 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 366 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 367 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            do i = 1, num_dims
                                vel_l(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%mom%beg + i - 1)
                                vel_r(i) = qr_prim_rsx_vf(j, k, l + 1, eqn_idx%mom%beg + i - 1)
                            end do
                        end if
 
                        call s_compute_interface_reynolds(alpha_l, re_l, re_size_loc1, re_size_loc2)
                        call s_compute_interface_reynolds(alpha_r, re_r, re_size_loc1, re_size_loc2)
 
                        if (shear_stress) then
 
# 377 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 377 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 377 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 377 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 377 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 378 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            do i = 1, num_dims
                                vel_grad_l(i, 1) = (dql_prim_dx_vf(eqn_idx%mom%beg + i - 1)%sf(j, k, l)/re_l(1))
                                vel_grad_r(i, 1) = (dqr_prim_dx_vf(eqn_idx%mom%beg + i - 1)%sf(idx_right_phys(1), &
                                           & idx_right_phys(2), idx_right_phys(3))/re_r(1))
# 383 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                    if (num_dims > 1) then
                                        vel_grad_l(i, 2) = (dql_prim_dy_vf(eqn_idx%mom%beg + i - 1)%sf(j, k, l)/re_l(1))
                                        vel_grad_r(i, 2) = (dqr_prim_dy_vf(eqn_idx%mom%beg + i - 1)%sf(idx_right_phys(1), &
                                                   & idx_right_phys(2), idx_right_phys(3))/re_r(1))
                                    end if
# 389 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                        if (num_dims > 2) then
                                            vel_grad_l(i, 3) = (dql_prim_dz_vf(eqn_idx%mom%beg + i - 1)%sf(j, k, l)/re_l(1))
                                            vel_grad_r(i, 3) = (dqr_prim_dz_vf(eqn_idx%mom%beg + i - 1)%sf(idx_right_phys(1), &
                                                       & idx_right_phys(2), idx_right_phys(3))/re_r(1))
                                        end if
# 395 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
# 396 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            end do
 
                            if (norm_dir == 1) then
                                flux_src_vf(eqn_idx%mom%beg)%sf(j, k, l) = flux_src_vf(eqn_idx%mom%beg)%sf(j, k, &
                                            & l) - (4._wp/3._wp)*0.5_wp*(vel_grad_l(1, 1) + vel_grad_r(1, 1))
                                flux_src_vf(eqn_idx%E)%sf(j, k, l) = flux_src_vf(eqn_idx%E)%sf(j, k, &
                                            & l) - (4._wp/3._wp)*0.5_wp*(vel_grad_l(1, 1)*vel_l(1) + vel_grad_r(1, 1)*vel_r(1))
# 404 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                    if (num_dims > 1) then
                                        flux_src_vf(eqn_idx%mom%beg)%sf(j, k, l) = flux_src_vf(eqn_idx%mom%beg)%sf(j, k, &
                                                    & l) - (-2._wp/3._wp)*0.5_wp*(vel_grad_l(2, 2) + vel_grad_r(2, 2))
                                        flux_src_vf(eqn_idx%E)%sf(j, k, l) = flux_src_vf(eqn_idx%E)%sf(j, k, &
                                                    & l) - (-2._wp/3._wp)*0.5_wp*(vel_grad_l(2, 2)*vel_l(1) + vel_grad_r(2, &
                                                    & 2)*vel_r(1))
 
                                        flux_src_vf(eqn_idx%mom%beg + 1)%sf(j, k, l) = flux_src_vf(eqn_idx%mom%beg + 1)%sf(j, k, &
                                                    & l) - 0.5_wp*(vel_grad_l(1, 2) + vel_grad_r(1, 2)) - 0.5_wp*(vel_grad_l(2, &
                                                    & 1) + vel_grad_r(2, 1))
                                        flux_src_vf(eqn_idx%E)%sf(j, k, l) = flux_src_vf(eqn_idx%E)%sf(j, k, &
                                                    & l) - 0.5_wp*(vel_grad_l(1, 2)*vel_l(2) + vel_grad_r(1, &
                                                    & 2)*vel_r(2)) - 0.5_wp*(vel_grad_l(2, 1)*vel_l(2) + vel_grad_r(2, 1)*vel_r(2))
# 418 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                            if (num_dims > 2) then
                                                flux_src_vf(eqn_idx%mom%beg)%sf(j, k, l) = flux_src_vf(eqn_idx%mom%beg)%sf(j, k, &
                                                            & l) - (-2._wp/3._wp)*0.5_wp*(vel_grad_l(3, 3) + vel_grad_r(3, 3))
                                                flux_src_vf(eqn_idx%E)%sf(j, k, l) = flux_src_vf(eqn_idx%E)%sf(j, k, &
                                                            & l) - (-2._wp/3._wp)*0.5_wp*(vel_grad_l(3, &
                                                            & 3)*vel_l(1) + vel_grad_r(3, 3)*vel_r(1))
 
                                                flux_src_vf(eqn_idx%mom%beg + 2)%sf(j, k, &
                                                            & l) = flux_src_vf(eqn_idx%mom%beg + 2)%sf(j, k, &
                                                            & l) - 0.5_wp*(vel_grad_l(1, 3) + vel_grad_r(1, &
                                                            & 3)) - 0.5_wp*(vel_grad_l(3, 1) + vel_grad_r(3, 1))
                                                flux_src_vf(eqn_idx%E)%sf(j, k, l) = flux_src_vf(eqn_idx%E)%sf(j, k, &
                                                            & l) - 0.5_wp*(vel_grad_l(1, 3)*vel_l(3) + vel_grad_r(1, &
                                                            & 3)*vel_r(3)) - 0.5_wp*(vel_grad_l(3, 1)*vel_l(3) + vel_grad_r(3, &
                                                            & 1)*vel_r(3))
                                            end if
# 435 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                    end if
# 437 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            else if (norm_dir == 2) then
# 439 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                    flux_src_vf(eqn_idx%mom%beg + 1)%sf(j, k, l) = flux_src_vf(eqn_idx%mom%beg + 1)%sf(j, k, &
                                                & l) - (-2._wp/3._wp)*0.5_wp*(vel_grad_l(1, 1) + vel_grad_r(1, 1))
                                    flux_src_vf(eqn_idx%E)%sf(j, k, l) = flux_src_vf(eqn_idx%E)%sf(j, k, &
                                                & l) - (-2._wp/3._wp)*0.5_wp*(vel_grad_l(1, 1)*vel_l(2) + vel_grad_r(1, 1)*vel_r(2))
 
                                    flux_src_vf(eqn_idx%mom%beg + 1)%sf(j, k, l) = flux_src_vf(eqn_idx%mom%beg + 1)%sf(j, k, &
                                                & l) - (4._wp/3._wp)*0.5_wp*(vel_grad_l(2, 2) + vel_grad_r(2, 2))
                                    flux_src_vf(eqn_idx%E)%sf(j, k, l) = flux_src_vf(eqn_idx%E)%sf(j, k, &
                                                & l) - (4._wp/3._wp)*0.5_wp*(vel_grad_l(2, 2)*vel_l(2) + vel_grad_r(2, 2)*vel_r(2))
 
                                    flux_src_vf(eqn_idx%mom%beg)%sf(j, k, l) = flux_src_vf(eqn_idx%mom%beg)%sf(j, k, &
                                                & l) - 0.5_wp*(vel_grad_l(1, 2) + vel_grad_r(1, 2)) - 0.5_wp*(vel_grad_l(2, &
                                                & 1) + vel_grad_r(2, 1))
                                    flux_src_vf(eqn_idx%E)%sf(j, k, l) = flux_src_vf(eqn_idx%E)%sf(j, k, &
                                                & l) - 0.5_wp*(vel_grad_l(1, 2)*vel_l(1) + vel_grad_r(1, &
                                                & 2)*vel_r(1)) - 0.5_wp*(vel_grad_l(2, 1)*vel_l(1) + vel_grad_r(2, 1)*vel_r(1))
# 456 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                        if (num_dims > 2) then
                                            flux_src_vf(eqn_idx%mom%beg + 1)%sf(j, k, l) = flux_src_vf(eqn_idx%mom%beg + 1)%sf(j, &
                                                        & k, l) - (-2._wp/3._wp)*0.5_wp*(vel_grad_l(3, 3) + vel_grad_r(3, 3))
                                            flux_src_vf(eqn_idx%E)%sf(j, k, l) = flux_src_vf(eqn_idx%E)%sf(j, k, &
                                                        & l) - (-2._wp/3._wp)*0.5_wp*(vel_grad_l(3, 3)*vel_l(2) + vel_grad_r(3, &
                                                        & 3)*vel_r(2))
 
                                            flux_src_vf(eqn_idx%mom%beg + 2)%sf(j, k, l) = flux_src_vf(eqn_idx%mom%beg + 2)%sf(j, &
                                                        & k, l) - 0.5_wp*(vel_grad_l(2, 3) + vel_grad_r(2, &
                                                        & 3)) - 0.5_wp*(vel_grad_l(3, 2) + vel_grad_r(3, 2))
                                            flux_src_vf(eqn_idx%E)%sf(j, k, l) = flux_src_vf(eqn_idx%E)%sf(j, k, &
                                                        & l) - 0.5_wp*(vel_grad_l(2, 3)*vel_l(3) + vel_grad_r(2, &
                                                        & 3)*vel_r(3)) - 0.5_wp*(vel_grad_l(3, 2)*vel_l(3) + vel_grad_r(3, &
                                                        & 2)*vel_r(3))
                                        end if
# 472 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
# 473 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            else
# 475 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                    flux_src_vf(eqn_idx%mom%beg + 2)%sf(j, k, l) = flux_src_vf(eqn_idx%mom%beg + 2)%sf(j, k, &
                                                & l) - (-2._wp/3._wp)*0.5_wp*(vel_grad_l(1, 1) + vel_grad_r(1, 1))
                                    flux_src_vf(eqn_idx%E)%sf(j, k, l) = flux_src_vf(eqn_idx%E)%sf(j, k, &
                                                & l) - (-2._wp/3._wp)*0.5_wp*(vel_grad_l(1, 1)*vel_l(3) + vel_grad_r(1, 1)*vel_r(3))
 
                                    flux_src_vf(eqn_idx%mom%beg + 2)%sf(j, k, l) = flux_src_vf(eqn_idx%mom%beg + 2)%sf(j, k, &
                                                & l) - (-2._wp/3._wp)*0.5_wp*(vel_grad_l(2, 2) + vel_grad_r(2, 2))
                                    flux_src_vf(eqn_idx%E)%sf(j, k, l) = flux_src_vf(eqn_idx%E)%sf(j, k, &
                                                & l) - (-2._wp/3._wp)*0.5_wp*(vel_grad_l(2, 2)*vel_l(3) + vel_grad_r(2, 2)*vel_r(3))
 
                                    flux_src_vf(eqn_idx%mom%beg)%sf(j, k, l) = flux_src_vf(eqn_idx%mom%beg)%sf(j, k, &
                                                & l) - 0.5_wp*(vel_grad_l(1, 3) + vel_grad_r(1, 3)) - 0.5_wp*(vel_grad_l(3, &
                                                & 1) + vel_grad_r(3, 1))
                                    flux_src_vf(eqn_idx%E)%sf(j, k, l) = flux_src_vf(eqn_idx%E)%sf(j, k, &
                                                & l) - 0.5_wp*(vel_grad_l(1, 3)*vel_l(1) + vel_grad_r(1, &
                                                & 3)*vel_r(1)) - 0.5_wp*(vel_grad_l(3, 1)*vel_l(1) + vel_grad_r(3, 1)*vel_r(1))
 
                                    flux_src_vf(eqn_idx%mom%beg + 2)%sf(j, k, l) = flux_src_vf(eqn_idx%mom%beg + 2)%sf(j, k, &
                                                & l) - (4._wp/3._wp)*0.5_wp*(vel_grad_l(3, 3) + vel_grad_r(3, 3))
                                    flux_src_vf(eqn_idx%E)%sf(j, k, l) = flux_src_vf(eqn_idx%E)%sf(j, k, &
                                                & l) - (4._wp/3._wp)*0.5_wp*(vel_grad_l(3, 3)*vel_l(3) + vel_grad_r(3, 3)*vel_r(3))
 
                                    flux_src_vf(eqn_idx%mom%beg + 1)%sf(j, k, l) = flux_src_vf(eqn_idx%mom%beg + 1)%sf(j, k, &
                                                & l) - 0.5_wp*(vel_grad_l(2, 3) + vel_grad_r(2, 3)) - 0.5_wp*(vel_grad_l(3, &
                                                & 2) + vel_grad_r(3, 2))
                                    flux_src_vf(eqn_idx%E)%sf(j, k, l) = flux_src_vf(eqn_idx%E)%sf(j, k, &
                                                & l) - 0.5_wp*(vel_grad_l(2, 3)*vel_l(2) + vel_grad_r(2, &
                                                & 3)*vel_r(2)) - 0.5_wp*(vel_grad_l(3, 2)*vel_l(2) + vel_grad_r(3, 2)*vel_r(2))
# 504 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            end if
                        end if
 
                        if (bulk_stress) then
 
# 508 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 508 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc loop seq 
# 508 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 508 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 508 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 509 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            do i = 1, num_dims
                                vel_grad_l(i, 1) = (dql_prim_dx_vf(eqn_idx%mom%beg + i - 1)%sf(j, k, l)/re_l(2))
                                vel_grad_r(i, 1) = (dqr_prim_dx_vf(eqn_idx%mom%beg + i - 1)%sf(idx_right_phys(1), &
                                           & idx_right_phys(2), idx_right_phys(3))/re_r(2))
# 514 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                    if (num_dims > 1) then
                                        vel_grad_l(i, 2) = (dql_prim_dy_vf(eqn_idx%mom%beg + i - 1)%sf(j, k, l)/re_l(2))
                                        vel_grad_r(i, 2) = (dqr_prim_dy_vf(eqn_idx%mom%beg + i - 1)%sf(idx_right_phys(1), &
                                                   & idx_right_phys(2), idx_right_phys(3))/re_r(2))
                                    end if
# 520 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
# 521 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                    if (num_dims > 2) then
                                        vel_grad_l(i, 3) = (dql_prim_dz_vf(eqn_idx%mom%beg + i - 1)%sf(j, k, l)/re_l(2))
                                        vel_grad_r(i, 3) = (dqr_prim_dz_vf(eqn_idx%mom%beg + i - 1)%sf(idx_right_phys(1), &
                                                   & idx_right_phys(2), idx_right_phys(3))/re_r(2))
                                    end if
# 527 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            end do
 
                            if (norm_dir == 1) then
                                flux_src_vf(eqn_idx%mom%beg)%sf(j, k, l) = flux_src_vf(eqn_idx%mom%beg)%sf(j, k, &
                                            & l) - 0.5_wp*(vel_grad_l(1, 1) + vel_grad_r(1, 1))
                                flux_src_vf(eqn_idx%E)%sf(j, k, l) = flux_src_vf(eqn_idx%E)%sf(j, k, l) - 0.5_wp*(vel_grad_l(1, &
                                            & 1)*vel_l(1) + vel_grad_r(1, 1)*vel_r(1))
# 535 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                    if (num_dims > 1) then
                                        flux_src_vf(eqn_idx%mom%beg)%sf(j, k, l) = flux_src_vf(eqn_idx%mom%beg)%sf(j, k, &
                                                    & l) - 0.5_wp*(vel_grad_l(2, 2) + vel_grad_r(2, 2))
                                        flux_src_vf(eqn_idx%E)%sf(j, k, l) = flux_src_vf(eqn_idx%E)%sf(j, k, &
                                                    & l) - 0.5_wp*(vel_grad_l(2, 2)*vel_l(1) + vel_grad_r(2, 2)*vel_r(1))
 
# 542 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                            if (num_dims > 2) then
                                                flux_src_vf(eqn_idx%mom%beg)%sf(j, k, l) = flux_src_vf(eqn_idx%mom%beg)%sf(j, k, &
                                                            & l) - 0.5_wp*(vel_grad_l(3, 3) + vel_grad_r(3, 3))
                                                flux_src_vf(eqn_idx%E)%sf(j, k, l) = flux_src_vf(eqn_idx%E)%sf(j, k, &
                                                            & l) - 0.5_wp*(vel_grad_l(3, 3)*vel_l(1) + vel_grad_r(3, 3)*vel_r(1))
                                            end if
# 549 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                    end if
# 551 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            else if (norm_dir == 2) then
# 553 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                    flux_src_vf(eqn_idx%mom%beg + 1)%sf(j, k, l) = flux_src_vf(eqn_idx%mom%beg + 1)%sf(j, k, &
                                                & l) - 0.5_wp*(vel_grad_l(1, 1) + vel_grad_r(1, 1))
                                    flux_src_vf(eqn_idx%E)%sf(j, k, l) = flux_src_vf(eqn_idx%E)%sf(j, k, &
                                                & l) - 0.5_wp*(vel_grad_l(1, 1)*vel_l(2) + vel_grad_r(1, 1)*vel_r(2))
 
                                    flux_src_vf(eqn_idx%mom%beg + 1)%sf(j, k, l) = flux_src_vf(eqn_idx%mom%beg + 1)%sf(j, k, &
                                                & l) - 0.5_wp*(vel_grad_l(2, 2) + vel_grad_r(2, 2))
                                    flux_src_vf(eqn_idx%E)%sf(j, k, l) = flux_src_vf(eqn_idx%E)%sf(j, k, &
                                                & l) - 0.5_wp*(vel_grad_l(2, 2)*vel_l(2) + vel_grad_r(2, 2)*vel_r(2))
 
# 564 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                        if (num_dims > 2) then
                                            flux_src_vf(eqn_idx%mom%beg + 1)%sf(j, k, l) = flux_src_vf(eqn_idx%mom%beg + 1)%sf(j, &
                                                        & k, l) - 0.5_wp*(vel_grad_l(3, 3) + vel_grad_r(3, 3))
                                            flux_src_vf(eqn_idx%E)%sf(j, k, l) = flux_src_vf(eqn_idx%E)%sf(j, k, &
                                                        & l) - 0.5_wp*(vel_grad_l(3, 3)*vel_l(2) + vel_grad_r(3, 3)*vel_r(2))
                                        end if
# 571 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
# 572 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            else
# 574 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                                    flux_src_vf(eqn_idx%mom%beg + 2)%sf(j, k, l) = flux_src_vf(eqn_idx%mom%beg + 2)%sf(j, k, &
                                                & l) - 0.5_wp*(vel_grad_l(1, 1) + vel_grad_r(1, 1))
                                    flux_src_vf(eqn_idx%E)%sf(j, k, l) = flux_src_vf(eqn_idx%E)%sf(j, k, &
                                                & l) - 0.5_wp*(vel_grad_l(1, 1)*vel_l(3) + vel_grad_r(1, 1)*vel_r(3))
 
                                    flux_src_vf(eqn_idx%mom%beg + 2)%sf(j, k, l) = flux_src_vf(eqn_idx%mom%beg + 2)%sf(j, k, &
                                                & l) - 0.5_wp*(vel_grad_l(2, 2) + vel_grad_r(2, 2))
                                    flux_src_vf(eqn_idx%E)%sf(j, k, l) = flux_src_vf(eqn_idx%E)%sf(j, k, &
                                                & l) - 0.5_wp*(vel_grad_l(2, 2)*vel_l(3) + vel_grad_r(2, 2)*vel_r(3))
 
                                    flux_src_vf(eqn_idx%mom%beg + 2)%sf(j, k, l) = flux_src_vf(eqn_idx%mom%beg + 2)%sf(j, k, &
                                                & l) - 0.5_wp*(vel_grad_l(3, 3) + vel_grad_r(3, 3))
                                    flux_src_vf(eqn_idx%E)%sf(j, k, l) = flux_src_vf(eqn_idx%E)%sf(j, k, &
                                                & l) - 0.5_wp*(vel_grad_l(3, 3)*vel_l(3) + vel_grad_r(3, 3)*vel_r(3))
# 589 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
                            end if
                        end if
                    end do
                end do
            end do
 
# 594 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#if defined(MFC_OpenACC)
# 594 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$acc end parallel loop
# 594 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#elif defined(MFC_OpenMP)
# 594 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
 
# 594 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
!$omp end target teams loop
# 594 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
#endif
# 595 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_lf.fpp"
        end if
 
        call s_finalize_riemann_solver(flux_vf, flux_src_vf, flux_gsrc_vf, norm_dir)
 
    end subroutine s_lf_riemann_solver
 
end module m_riemann_solver_lf