m_phase_change.fpp.f90

Go to the documentation of this file.

# 1 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!>
!! @file
!! @brief Contains module m_phase_change
 
# 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"
# 6 "/home/runner/work/MFC/MFC/src/common/m_phase_change.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"
 
# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
! New line at end of file is required for FYPP
# 2 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
# 1 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 1
# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
! New line at end of file is required for FYPP
# 2 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 2
 
# 4 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 5 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 6 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 7 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 8 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 20 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 43 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 48 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 53 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 58 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 63 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 68 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 76 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 81 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 86 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 91 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 96 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 101 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 106 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 111 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 116 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 121 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 151 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 192 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 206 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 231 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 242 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 244 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 255 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 284 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 294 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 304 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 313 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 330 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 340 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 347 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 353 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 359 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 365 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 371 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 377 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
! New line at end of file is required for FYPP
# 3 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
# 1 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 1
# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
! New line at end of file is required for FYPP
# 2 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 2
 
# 7 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 17 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 22 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 27 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 32 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 37 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 42 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 47 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 52 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 57 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 62 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 73 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 78 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 83 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 88 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 103 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 131 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 160 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 175 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 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. NVIDIA unified memory page placement hint
# 57 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
! Allocate and create GPU device memory
# 77 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
! Free GPU device memory and deallocate
# 85 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
! Cray-specific GPU pointer setup for vector fields
# 109 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
! Cray-specific GPU pointer setup for scalar fields
# 125 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
! Cray-specific GPU pointer setup for acoustic source spatials
# 150 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 156 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 163 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
! New line at end of file is required for FYPP
# 7 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp" 2
 
!> @brief Phase transition relaxation solvers for liquid-vapor flows with cavitation and boiling
module m_phase_change
 
#ifndef MFC_POST_PROCESS
    use m_derived_types
    use m_global_parameters
    use m_mpi_proxy
    use m_variables_conversion
    use ieee_arithmetic
    use m_helper_basic
    use m_constants, only: model_eqns_6eq
 
    implicit none
 
    private
    public :: s_initialize_phasechange_module, s_relaxation_solver, s_infinite_relaxation_k, s_finalize_relaxation_solver_module
 
    !> @name Parameters for the first order transition phase change
    !> @{
    integer, parameter  :: max_iter = 1e8_wp            !< max # of iterations
    real(wp), parameter :: pcr = 4.94e7_wp              !< Critical pressure of water [Pa]
    real(wp), parameter :: tcr = 385.05_wp + 273.15_wp  !< Critical temperature of water [K]
    real(wp), parameter :: mixm = 1.0e-8_wp             !< Mixture mass fraction threshold for triggering phase change
    integer, parameter  :: lp = 1                       !< index for the liquid phase of the reacting fluid
    integer, parameter  :: vp = 2                       !< index for the vapor phase of the reacting fluid
    !> @}
 
    !> @name Gibbs free energy phase change parameters
    !> @{
    real(wp) :: a, b, c, d
    !> @}
 
 
# 40 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if defined(MFC_OpenACC)
# 40 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc declare create(A, B, C, D) 
# 40 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif defined(MFC_OpenMP)
# 40 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$omp declare target (A, B, C, D) 
# 40 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
 
contains
 
    !> Dispatch to the correct relaxation solver. Replaces the procedure pointer, which CCE is breaking on.
    impure subroutine s_relaxation_solver(q_cons_vf)
 
        type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
        ! This is empty because in current master the procedure pointer was never assigned
 
    if (.not. (.false.)) then
# 50 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
        call s_mpi_abort("m_phase_change.fpp:50: " // "Assertion failed: .false.. " &
# 50 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
                         & // "s_relaxation_solver called but it currently does nothing")
# 50 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
    end if
 
    end subroutine s_relaxation_solver
 
    !> Initialize the phase change module by setting saturation curve coefficients for pT- or pTg-equilibrium
    impure subroutine s_initialize_phasechange_module
 
        ! Saturation curve coefficients via stiffened gas EOS. Saurel et al. JCP (2008), Le Metayer et al. JFE (2004)
        a = (gs_min(lp)*cvs(lp) - gs_min(vp)*cvs(vp) + qvps(vp) - qvps(lp))/((gs_min(vp) - 1.0_wp)*cvs(vp))
 
        b = (qvs(lp) - qvs(vp))/((gs_min(vp) - 1.0_wp)*cvs(vp))
 
        c = (gs_min(vp)*cvs(vp) - gs_min(lp)*cvs(lp))/((gs_min(vp) - 1.0_wp)*cvs(vp))
 
        d = ((gs_min(lp) - 1.0_wp)*cvs(lp))/((gs_min(vp) - 1.0_wp)*cvs(vp))
 
    end subroutine s_initialize_phasechange_module
 
    !> Apply pT- or pTg-equilibrium relaxation with mass depletion based on the incoming state conditions.
    subroutine s_infinite_relaxation_k(q_cons_vf)
 
        type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
        real(wp) :: ps, psov, pssl                  !< equilibrium pressure for mixture, overheated vapor, and subcooled liquid
        real(wp) :: ts, tsov, tssl, tsatov, tsatsl  !< Equilibrium and saturation temperatures
        real(wp) :: rhoe, dyne, rhos                !< total internal energy, kinetic energy, and total entropy
        real(wp) :: rho, rm, m1, m2, mct            !< total density, total reacting mass, individual reacting masses
        real(wp) :: tvf                             !< total volume fraction
        ! $:GPU_DECLARE(create='[pS,pSOV,pSSL,TS,TSOV,TSSL,TSatOV,TSatSL]')
        ! $:GPU_DECLARE(create='[rhoe,dynE,rhos,rho,rM,m1,m2,MCT,TvF]')
 
# 83 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
            real(wp), dimension(num_fluids) :: p_infov, p_infpt, p_infsl, sk, hk, gk, ek, rhok
# 85 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
        ! $:GPU_DECLARE(create='[p_infOV,p_infpT,p_infSL,sk,hk,gk,ek,rhok]')
 
        !> Generic loop iterators
        integer :: i, j, k, l
 
#ifdef _CRAYFTN
#ifdef MFC_OpenACC
        ! CCE 19 IPA workaround: prevent bring_routine_resident SIGSEGV DIR$ NOINLINE s_infinite_pt_relaxation_k DIR$ NOINLINE
        ! s_infinite_ptg_relaxation_k DIR$ NOINLINE s_correct_partial_densities DIR$ NOINLINE s_TSat
#endif
#endif
 
        ! starting equilibrium solver
 
 
# 99 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 99 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if defined(MFC_OpenACC)
# 99 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, l, p_infOV, p_infpT, p_infSL, sk, hk, gk, ek, rhok, pS, pSOV, pSSL, TS, TSOV, TSatOV, TSatSL, TSSL, rhoe, dynE, rhos, rho, rM, m1, m2, MCT, TvF) 
# 99 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif defined(MFC_OpenMP)
# 99 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 99 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 99 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 99 "/home/runner/work/MFC/MFC/src/common/m_phase_change.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, p_infOV, p_infpT, p_infSL, sk, hk, gk, ek, rhok, pS, pSOV, pSSL, TS, TSOV, TSatOV, TSatSL, TSSL, rhoe, dynE, rhos, rho, rM, m1, m2, MCT, TvF) 
# 99 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
# 101 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
        do j = 0, m
            do k = 0, n
                do l = 0, p
                    rho = 0.0_wp; tvf = 0.0_wp
 
# 105 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if defined(MFC_OpenACC)
# 105 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc loop seq 
# 105 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif defined(MFC_OpenMP)
# 105 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 105 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
                    do i = 1, num_fluids
                        ! Mixture density
                        rho = rho + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)
 
                        ! Total Volume Fraction
                        tvf = tvf + q_cons_vf(i + eqn_idx%adv%beg - 1)%sf(j, k, l)
                    end do
 
                    ! calculating the total reacting mass for the phase change process. By hypothesis, this should not change
                    ! throughout the phase-change process.
                    rm = q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) + q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l)
 
                    ! correcting negative (reacting) mass fraction values in case they happen
                    call s_correct_partial_densities(mct, q_cons_vf, rm, j, k, l)
 
                    ! fixing m1 and m2 AFTER correcting the partial densities. Note that these values must be stored for the phase
                    ! change process that will happen a posteriori
                    m1 = q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l)
 
                    m2 = q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l)
 
                    ! kinetic energy as an auxiliary variable to the calculation of the total internal energy
                    dyne = 0.0_wp
 
# 129 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if defined(MFC_OpenACC)
# 129 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc loop seq 
# 129 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif defined(MFC_OpenMP)
# 129 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 129 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
                    do i = eqn_idx%mom%beg, eqn_idx%mom%end
                        dyne = dyne + 5.0e-1_wp*q_cons_vf(i)%sf(j, k, l)**2/rho
                    end do
 
                    ! calculating the total energy that MUST be preserved throughout the pT- and pTg-relaxation procedures at each
                    ! of the cells. The internal energy is calculated as the total energy minus the kinetic energy to preserved its
                    ! value at sharp interfaces
                    rhoe = q_cons_vf(eqn_idx%E)%sf(j, k, l) - dyne
 
                    ! Calling pT-equilibrium for either finishing phase-change module, or as an IC for the pTg-equilibrium for this
                    ! case, MFL cannot be either 0 or 1, so I chose it to be 2
                    call s_infinite_pt_relaxation_k(j, k, l, 2, ps, p_infpt, q_cons_vf, rhoe, ts)
 
                    ! Check if pTg-equilibrium needed; only partial densities require updating
                    if ((relax_model == 6) .and. ((q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, &
                        & l) > mixm*rm) .and. (q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, &
                        & l) > mixm*rm)) .and. (ps < pcr) .and. (ts < tcr)) then
                        ! Checking if phase change is needed, by checking whether the final solution is either subcoooled liquid or
                        ! overheated vapor.
 
                        ! overheated vapor case depleting the mass of liquid
                        q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) = mixm*rm
 
                        ! transferring the total mass to vapor
                        q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) = (1.0_wp - mixm)*rm
 
                        ! calling pT-equilibrium for overheated vapor, which is MFL = 0
                        call s_infinite_pt_relaxation_k(j, k, l, 0, psov, p_infov, q_cons_vf, rhoe, tsov)
 
                        ! calculating Saturation temperature
                        call s_tsat(psov, tsatov, tsov)
 
                        ! subcooled liquid case transferring the total mass to liquid
                        q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) = (1.0_wp - mixm)*rm
 
                        ! depleting the mass of vapor
                        q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) = mixm*rm
 
                        ! calling pT-equilibrium for subcooled liquid, which is MFL = 1
                        call s_infinite_pt_relaxation_k(j, k, l, 1, pssl, p_infsl, q_cons_vf, rhoe, tssl)
 
                        ! calculating Saturation temperature
                        call s_tsat(pssl, tsatsl, tssl)
 
                        ! checking the conditions for overheated vapor and subcooled liquide
                        if (tsov > tsatov) then
                            ! Assigning pressure
                            ps = psov
 
                            ! Assigning Temperature
                            ts = tsov
 
                            ! correcting the liquid partial density
                            q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) = mixm*rm
 
                            ! correcting the vapor partial density
                            q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) = (1.0_wp - mixm)*rm
                        else if (tssl < tsatsl) then
                            ! Assigning pressure
                            ps = pssl
 
                            ! Assigning Temperature
                            ts = tssl
 
                            ! correcting the liquid partial density
                            q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) = (1.0_wp - mixm)*rm
 
                            ! correcting the vapor partial density
                            q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) = mixm*rm
                        else
                            ! returning partial pressures to what they were from the homogeneous solver liquid
                            q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) = m1
 
                            ! vapor
                            q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) = m2
 
                            ! calling the pTg-equilibrium solver
                            call s_infinite_ptg_relaxation_k(j, k, l, ps, p_infpt, rhoe, q_cons_vf, ts)
                        end if
                    end if
 
                    ! Calculations AFTER equilibrium
 
 
# 213 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if defined(MFC_OpenACC)
# 213 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc loop seq 
# 213 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif defined(MFC_OpenMP)
# 213 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 213 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
                    do i = 1, num_fluids
                        ! entropy
                        sk(i) = cvs(i)*log((ts**gs_min(i))/((ps + ps_inf(i))**(gs_min(i) - 1.0_wp))) + qvps(i)
 
                        ! enthalpy
                        hk(i) = gs_min(i)*cvs(i)*ts + qvs(i)
 
                        ! Gibbs-free energy
                        gk(i) = hk(i) - ts*sk(i)
 
                        ! densities
                        rhok(i) = (ps + ps_inf(i))/((gs_min(i) - 1)*cvs(i)*ts)
 
                        ! internal energy
                        ek(i) = (ps + gs_min(i)*ps_inf(i))/(ps + ps_inf(i))*cvs(i)*ts + qvs(i)
                    end do
 
                    ! calculating volume fractions, internal energies, and total entropy
                    rhos = 0.0_wp
 
# 233 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if defined(MFC_OpenACC)
# 233 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc loop seq 
# 233 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif defined(MFC_OpenMP)
# 233 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 233 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
                    do i = 1, num_fluids
                        ! volume fractions
                        q_cons_vf(i + eqn_idx%adv%beg - 1)%sf(j, k, l) = q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)/rhok(i)
 
                        ! alpha*rho*e
                        if (model_eqns == model_eqns_6eq) then
                            q_cons_vf(i + eqn_idx%int_en%beg - 1)%sf(j, k, l) = q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, &
                                      & l)*ek(i)
                        end if
 
                        ! Total entropy
                        rhos = rhos + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*sk(i)
                    end do
                end do
            end do
        end do
 
# 250 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if defined(MFC_OpenACC)
# 250 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc end parallel loop
# 250 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif defined(MFC_OpenMP)
# 250 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 250 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$omp end target teams loop
# 250 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
 
    end subroutine s_infinite_relaxation_k
 
    !> Apply pT-equilibrium relaxation for N fluids
    !! @param MFL flag: 0=gas, 1=liquid, 2=mixture
    subroutine s_infinite_pt_relaxation_k(j, k, l, MFL, pS, p_infpT, q_cons_vf, rhoe, TS)
 
 
# 258 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#ifdef _CRAYFTN
# 258 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if MFC_OpenACC
# 258 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc routine seq 
# 258 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif MFC_OpenMP
# 258 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 258 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 258 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$omp declare target device_type(any) 
# 258 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#else
# 258 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!DIR$ NOINLINE s_infinite_pt_relaxation_k
# 258 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
# 258 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif MFC_OpenACC
# 258 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc routine seq 
# 258 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif MFC_OpenMP
# 258 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 258 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 258 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$omp declare target device_type(any) 
# 258 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
 
        ! initializing variables
        integer, intent(in)                                 :: j, k, l, MFL
        real(wp), intent(out)                               :: pS
        real(wp), dimension(1:), intent(out)                :: p_infpT
        type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf
        real(wp), intent(in)                                :: rhoe
        real(wp), intent(out)                               :: TS
        real(wp)                                            :: gp, gpp, hp, pO, mCP, mQ  !< variables for the Newton Solver
        real(wp)                                            :: p_infpT_sum
        integer                                             :: i, ns                     !< generic loop iterators
        ! auxiliary variables for the pT-equilibrium solver
        mcp = 0.0_wp; mq = 0.0_wp; p_infpt_sum = 0._wp
 
# 272 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if defined(MFC_OpenACC)
# 272 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc loop seq 
# 272 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif defined(MFC_OpenMP)
# 272 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 272 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
        do i = 1, num_fluids
            p_infpt(i) = ps_inf(i)
            p_infpt_sum = p_infpt_sum + abs(p_infpt(i))
        end do
        ! Performing tests before initializing the pT-equilibrium
 
# 278 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if defined(MFC_OpenACC)
# 278 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc loop seq 
# 278 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif defined(MFC_OpenMP)
# 278 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 278 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
        do i = 1, num_fluids
            ! sum of the total alpha*rho*cp of the system
            mcp = mcp + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*cvs(i)*gs_min(i)
 
            ! sum of the total alpha*rho*q of the system
            mq = mq + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*qvs(i)
        end do
 
# 295 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
        ! Checking energy constraint
        if ((rhoe - mq - minval(p_infpt)) < 0.0_wp) then
            if ((mfl == 0) .or. (mfl == 1)) then
                ! Assigning zero values for mass depletion cases pressure
                ps = 0.0_wp
 
                ! temperature
                ts = 0.0_wp
 
                return
            end if
        end if
 
        ! calculating initial estimate for pressure in the pT-relaxation procedure. I will also use this variable to iterate over
        ! the Newton's solver
        po = 0.0_wp
 
        ! Maybe improve this condition afterwards. As long as the initial guess is in between -min(ps_inf) and infinity, a solution
        ! should be able to be found.
        ps = 1.0e4_wp
 
        ! Newton Solver for the pT-equilibrium
        ns = 0
        ! change this relative error metric. 1.e4_wp is just arbitrary
        do while ((abs(ps - po) > palpha_eps) .and. (abs((ps - po)/po) > palpha_eps/1.e4_wp) .or. (ns == 0))
            ! increasing counter
            ns = ns + 1
 
            ! updating old pressure
            po = ps
 
            ! updating functions used in the Newton's solver
            gpp = 0.0_wp; gp = 0.0_wp; hp = 0.0_wp
 
# 329 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if defined(MFC_OpenACC)
# 329 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc loop seq 
# 329 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif defined(MFC_OpenMP)
# 329 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 329 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
            do i = 1, num_fluids
                gp = gp + (gs_min(i) - 1.0_wp)*q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, &
                           & l)*cvs(i)*(rhoe + ps - mq)/(mcp*(ps + p_infpt(i)))
 
                gpp = gpp + (gs_min(i) - 1.0_wp)*q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, &
                             & l)*cvs(i)*(p_infpt(i) - rhoe + mq)/(mcp*(ps + p_infpt(i))**2)
            end do
 
            hp = 1.0_wp/(rhoe + ps - mq) + 1.0_wp/(ps + minval(p_infpt))
 
            ! updating common pressure for the newton solver
            ps = po + ((1.0_wp - gp)/gpp)/(1.0_wp - (1.0_wp - gp + abs(1.0_wp - gp))/(2.0_wp*gpp)*hp)
        end do
 
        ! common temperature
        ts = (rhoe + ps - mq)/mcp
 
    end subroutine s_infinite_pt_relaxation_k
 
    !> Apply pTg-equilibrium relaxation for N fluids under pT and 2 fluids under pTg-equilibrium. There is a final common p and T
    !! during relaxation
    subroutine s_infinite_ptg_relaxation_k(j, k, l, pS, p_infpT, rhoe, q_cons_vf, TS)
 
 
# 353 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#ifdef _CRAYFTN
# 353 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if MFC_OpenACC
# 353 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc routine seq 
# 353 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif MFC_OpenMP
# 353 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 353 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 353 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$omp declare target device_type(any) 
# 353 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#else
# 353 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!DIR$ NOINLINE s_infinite_ptg_relaxation_k
# 353 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
# 353 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif MFC_OpenACC
# 353 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc routine seq 
# 353 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif MFC_OpenMP
# 353 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 353 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 353 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$omp declare target device_type(any) 
# 353 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
 
        integer, intent(in)                                    :: j, k, l
        real(wp), intent(inout)                                :: pS
        real(wp), dimension(1:), intent(in)                    :: p_infpT
        real(wp), intent(in)                                   :: rhoe
        type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
        real(wp), intent(inout)                                :: TS
# 364 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
            real(wp), dimension(num_fluids) :: p_infpTg  !< stiffness for the participating fluids for pTg-equilibrium
# 366 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
        real(wp), dimension(2, 2) :: Jac, InvJac, TJac                  !< matrices for the Newton Solver
        real(wp), dimension(2)    :: R2D, DeltamP                       !< residual and correction array
        real(wp)                  :: Om                                 !< underrelaxation factor
        real(wp)                  :: mCP, mCPD, mCVGP, mCVGP2, mQ, mQD  !< auxiliary variables for the pTg-solver
        real(wp)                  :: ml, mT, dFdT, dTdm, dTdp
 
        !> Generic loop iterators
        integer :: i, ns
        ! pTg-equilibrium solution procedure Newton Solver parameters counter
        ns = 0
 
        ! Relaxation factor
        om = 1.0e-3_wp
 
        p_infptg = p_infpt
 
        if (((ps < 0.0_wp) .and. ((q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) + q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, &
            & k, &
            & l)) > ((rhoe - gs_min(lp)*ps_inf(lp)/(gs_min(lp) - 1))/qvs(lp)))) .or. ((ps >= 0.0_wp) .and. (ps < 1.0e-1_wp))) then
            ! improve this initial condition
            ps = 1.0e4_wp
        end if
 
        ! Loop until the solution for F(X) is satisfied Check whether I need to use both absolute and relative values for the
        ! residual, and how to do it adequately. Dummy guess to start the pTg-equilibrium problem. improve this initial condition
        r2d(1) = 0.0_wp; r2d(2) = 0.0_wp
        deltamp(1) = 0.0_wp; deltamp(2) = 0.0_wp
        do while (((sqrt(r2d(1)**2 + r2d(2)**2) > ptgalpha_eps) .and. ((sqrt(r2d(1)**2 + r2d(2)**2)/rhoe) > (ptgalpha_eps/1.e6_wp) &
                  & )) .or. (ns == 0))
 
            ! Updating counter for the iterative procedure
            ns = ns + 1
 
            ! Auxiliary variables to help in the calculation of the residue
            mcp = 0.0_wp; mcpd = 0.0_wp; mcvgp = 0.0_wp; mcvgp2 = 0.0_wp; mq = 0.0_wp; mqd = 0.0_wp
            ! Those must be updated through the iterations, as they either depend on the partial masses for all fluids, or on the
            ! equilibrium pressure
 
# 403 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if defined(MFC_OpenACC)
# 403 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc loop seq 
# 403 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif defined(MFC_OpenMP)
# 403 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 403 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
            do i = 1, num_fluids
                ! sum of the total alpha*rho*cp of the system
                mcp = mcp + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*cvs(i)*gs_min(i)
 
                ! sum of the total alpha*rho*q of the system
                mq = mq + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*qvs(i)
 
                ! These auxiliary variables now need to be updated, as the partial densities now vary at every iteration
                if ((i /= lp) .and. (i /= vp)) then
                    mcvgp = mcvgp + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*cvs(i)*(gs_min(i) - 1)/(ps + ps_inf(i))
 
                    mcvgp2 = mcvgp2 + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*cvs(i)*(gs_min(i) - 1)/((ps + ps_inf(i))**2)
 
                    mqd = mqd + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*qvs(i)
 
                    ! sum of the total alpha*rho*cp of the system
                    mcpd = mcpd + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*cvs(i)*gs_min(i)
                end if
            end do
 
            ! calculating the (2D) Jacobian Matrix used in the solution of the pTg-quilibrium model
 
            ! mass of the reacting liquid
            ml = q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l)
 
            ! mass of the two participating fluids
            mt = q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) + q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l)
 
            ts = 1/(mt*cvs(vp)*(gs_min(vp) - 1)/(ps + ps_inf(vp)) + ml*(cvs(lp)*(gs_min(lp) - 1)/(ps + ps_inf(lp)) - cvs(vp) &
                    & *(gs_min(vp) - 1)/(ps + ps_inf(vp))) + mcvgp)
 
            dfdt = -(cvs(lp)*gs_min(lp) - cvs(vp)*gs_min(vp))*log(ts) - (qvps(lp) - qvps(vp)) + cvs(lp)*(gs_min(lp) - 1)*log(ps &
                     & + ps_inf(lp)) - cvs(vp)*(gs_min(vp) - 1)*log(ps + ps_inf(vp))
 
            dtdm = -(cvs(lp)*(gs_min(lp) - 1)/(ps + ps_inf(lp)) - cvs(vp)*(gs_min(vp) - 1)/(ps + ps_inf(vp)))*ts**2
 
            dtdp = (mt*cvs(vp)*(gs_min(vp) - 1)/(ps + ps_inf(vp))**2 + ml*(cvs(lp)*(gs_min(lp) - 1)/(ps + ps_inf(lp))**2 - cvs(vp) &
                    & *(gs_min(vp) - 1)/(ps + ps_inf(vp))**2) + mcvgp2)*ts**2
 
            ! F = (F1,F2) is the function whose roots we are looking for x = (m1, p) are the independent variables. m1 = mass of the
            ! first participant fluid, p = pressure F1 = 0 is the Gibbs free energy quality F2 = 0 is the enforcement of the
            ! thermodynamic (total - kinectic) energy dF1dm
            jac(1, 1) = dfdt*dtdm
 
            ! dF1dp
            jac(1, 2) = dfdt*dtdp + ts*(cvs(lp)*(gs_min(lp) - 1)/(ps + ps_inf(lp)) - cvs(vp)*(gs_min(vp) - 1)/(ps + ps_inf(vp)))
 
            ! dF2dm
            jac(2, &
                & 1) = (qvs(vp) - qvs(lp) + (cvs(vp)*gs_min(vp) - cvs(lp)*gs_min(lp))/(ml*(cvs(lp)*(gs_min(lp) - 1)/(ps &
                & + ps_inf(lp)) - cvs(vp)*(gs_min(vp) - 1)/(ps + ps_inf(vp))) + mt*cvs(vp)*(gs_min(vp) - 1)/(ps + ps_inf(vp)) &
                & + mcvgp) - (ml*(cvs(vp)*gs_min(vp) - cvs(lp)*gs_min(lp)) - mt*cvs(vp)*gs_min(vp) - mcpd)*(cvs(lp)*(gs_min(lp) &
                & - 1)/(ps + ps_inf(lp)) - cvs(vp)*(gs_min(vp) - 1)/(ps + ps_inf(vp)))/((ml*(cvs(lp)*(gs_min(lp) - 1)/(ps &
                & + ps_inf(lp)) - cvs(vp)*(gs_min(vp) - 1)/(ps + ps_inf(vp))) + mt*cvs(vp)*(gs_min(vp) - 1)/(ps + ps_inf(vp)) &
                & + mcvgp)**2))/1
            ! dF2dp
            jac(2, &
                & 2) = (1 + (ml*(cvs(vp)*gs_min(vp) - cvs(lp)*gs_min(lp)) - mt*cvs(vp)*gs_min(vp) - mcpd)*(ml*(cvs(lp)*(gs_min(lp) &
                & - 1)/(ps + ps_inf(lp))**2 - cvs(vp)*(gs_min(vp) - 1)/(ps + ps_inf(vp))**2) + mt*cvs(vp)*(gs_min(vp) - 1)/(ps &
                & + ps_inf(vp))**2 + mcvgp2)/(ml*(cvs(lp)*(gs_min(lp) - 1)/(ps + ps_inf(lp)) - cvs(vp)*(gs_min(vp) - 1)/(ps &
                & + ps_inf(vp))) + mt*cvs(vp)*(gs_min(vp) - 1)/(ps + ps_inf(vp)) + mcvgp)**2)/1
 
            ! intermediate elements of J^{-1}
            invjac(1, 1) = jac(2, 2)
            invjac(1, 2) = -1.0_wp*jac(1, 2)
            invjac(2, 1) = -1.0_wp*jac(2, 1)
            invjac(2, 2) = jac(1, 1)
 
            ! elements of J^{T}
            tjac(1, 1) = jac(1, 1)
            tjac(1, 2) = jac(2, 1)
            tjac(2, 1) = jac(1, 2)
            tjac(2, 2) = jac(2, 2)
 
            ! dividing by det(J)
            invjac = invjac/(jac(1, 1)*jac(2, 2) - jac(1, 2)*jac(2, 1))
 
            ! calculating correction array for Newton's method
            deltamp(1) = -1.0_wp*(invjac(1, 1)*r2d(1) + invjac(1, 2)*r2d(2))
            deltamp(2) = -1.0_wp*(invjac(2, 1)*r2d(1) + invjac(2, 2)*r2d(2))
 
            ! updating two reacting 'masses'. Recall that inert 'masses' do not change during the phase change liquid
            q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) = q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) + om*deltamp(1)
 
            ! gas
            q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) = q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) - om*deltamp(1)
 
            ! updating pressure
            ps = ps + om*deltamp(2)
 
            ! calculating residuals, which are (i) the difference between the Gibbs Free energy of the gas and the liquid and (ii)
            ! the energy before and after the phase-change process.
 
            ! mass of the reacting liquid
            ml = q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l)
 
            ! mass of the two participating fluids
            mt = q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) + q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l)
 
            ts = 1/(mt*cvs(vp)*(gs_min(vp) - 1)/(ps + ps_inf(vp)) + ml*(cvs(lp)*(gs_min(lp) - 1)/(ps + ps_inf(lp)) - cvs(vp) &
                    & *(gs_min(vp) - 1)/(ps + ps_inf(vp))) + mcvgp)
 
            ! Gibbs Free Energy Equality condition (DG)
            r2d(1) = ts*((cvs(lp)*gs_min(lp) - cvs(vp)*gs_min(vp))*(1 - log(ts)) - (qvps(lp) - qvps(vp)) + cvs(lp)*(gs_min(lp) &
                & - 1)*log(ps + ps_inf(lp)) - cvs(vp)*(gs_min(vp) - 1)*log(ps + ps_inf(vp))) + qvs(lp) - qvs(vp)
 
            ! Constant Energy Process condition (DE)
            r2d(2) = (rhoe + ps + ml*(qvs(vp) - qvs(lp)) - mt*qvs(vp) - mqd + (ml*(gs_min(vp)*cvs(vp) - gs_min(lp)*cvs(lp)) &
                & - mt*gs_min(vp)*cvs(vp) - mcpd)/(ml*(cvs(lp)*(gs_min(lp) - 1)/(ps + ps_inf(lp)) - cvs(vp)*(gs_min(vp) - 1)/(ps &
                & + ps_inf(vp))) + mt*cvs(vp)*(gs_min(vp) - 1)/(ps + ps_inf(vp)) + mcvgp))/1
        end do
 
        ! common temperature
        ts = (rhoe + ps - mq)/mcp
 
    end subroutine s_infinite_ptg_relaxation_k
 
    !> Correct the partial densities of the reacting fluids in case one of them is negative but their sum is positive. Inert phases
    !! are not corrected at this moment
    subroutine s_correct_partial_densities(MCT, q_cons_vf, rM, j, k, l)
 
 
# 525 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#ifdef _CRAYFTN
# 525 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if MFC_OpenACC
# 525 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc routine seq 
# 525 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif MFC_OpenMP
# 525 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 525 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 525 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$omp declare target device_type(any) 
# 525 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#else
# 525 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!DIR$ NOINLINE s_correct_partial_densities
# 525 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
# 525 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif MFC_OpenACC
# 525 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc routine seq 
# 525 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif MFC_OpenMP
# 525 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 525 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 525 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$omp declare target device_type(any) 
# 525 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
 
        !> @name variables for the correction of the reacting partial densities
        !> @{
        real(wp), intent(out)                                  :: mct
        type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
        real(wp), intent(inout)                                :: rm
        integer, intent(in)                                    :: j, k, l
        !> @}
        if (rm < 0.0_wp) then
            if ((q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, &
                & l) >= -1.0_wp*mixm) .and. (q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) >= -1.0_wp*mixm)) then
                q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) = 0.0_wp
 
                q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) = 0.0_wp
 
                rm = q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) + q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l)
            end if
        end if
 
        ! TODO: Consider partitioning partial densities instead of absolute-value correction
        mct = 2*mixm
 
        ! correcting the partial densities of the reacting fluids. What to do for the nonreacting ones?
        if (q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) < 0.0_wp) then
            q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) = mct*rm
 
            q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) = (1.0_wp - mct)*rm
        else if (q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) < 0.0_wp) then
            q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) = (1.0_wp - mct)*rm
 
            q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) = mct*rm
        end if
 
    end subroutine s_correct_partial_densities
 
    !> Find the saturation temperature for a given saturation pressure using a Newton solver
    elemental subroutine s_tsat(pSat, TSat, TSIn)
 
 
# 564 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#ifdef _CRAYFTN
# 564 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if MFC_OpenACC
# 564 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc routine seq 
# 564 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif MFC_OpenMP
# 564 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 564 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 564 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$omp declare target device_type(any) 
# 564 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#else
# 564 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!DIR$ NOINLINE s_TSat
# 564 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
# 564 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif MFC_OpenACC
# 564 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc routine seq 
# 564 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif MFC_OpenMP
# 564 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 564 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 564 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$omp declare target device_type(any) 
# 564 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
 
        real(wp), intent(in)  :: psat
        real(wp), intent(out) :: tsat
        real(wp), intent(in)  :: tsin
        real(wp)              :: dfdt, ft, om  !< auxiliary variables
        ! Generic loop iterators
        integer :: ns
 
        if ((f_approx_equal(psat, 0.0_wp)) .and. (f_approx_equal(tsin, 0.0_wp))) then
            ! assigning Saturation temperature
            tsat = 0.0_wp
        else
            ! calculating initial estimate for temperature in the TSat procedure. I will also use this variable to iterate over the
            ! Newton's solver
            tsat = tsin
 
            ! iteration counter
            ns = 0
 
            ! underrelaxation factor
            om = 1.0e-3_wp
 
            ! FT must be initialized before the do while condition is evaluated. Fortran .or. is not short-circuit: abs(FT) is
            ! always evaluated even when ns == 0, so FT must have a defined value here.
            ft = huge(1.0_wp)
 
            do while ((abs(ft) > ptgalpha_eps) .or. (ns == 0))
                ! increasing counter
                ns = ns + 1
 
                ! calculating residual
                ft = tsat*((cvs(lp)*gs_min(lp) - cvs(vp)*gs_min(vp))*(1 - log(tsat)) - (qvps(lp) - qvps(vp)) + cvs(lp)*(gs_min(lp) &
                           & - 1)*log(psat + ps_inf(lp)) - cvs(vp)*(gs_min(vp) - 1)*log(psat + ps_inf(vp))) + qvs(lp) - qvs(vp)
 
                ! calculating the jacobian
                dfdt = -(cvs(lp)*gs_min(lp) - cvs(vp)*gs_min(vp))*log(tsat) - (qvps(lp) - qvps(vp)) + cvs(lp)*(gs_min(lp) - 1) &
                         & *log(psat + ps_inf(lp)) - cvs(vp)*(gs_min(vp) - 1)*log(psat + ps_inf(vp))
 
                ! updating saturation temperature
                tsat = tsat - om*ft/dfdt
 
                if (abs(ft) <= ptgalpha_eps) exit
            end do
        end if
 
    end subroutine s_tsat
 
    !> Finalize the phase change module
    impure subroutine s_finalize_relaxation_solver_module
 
    end subroutine s_finalize_relaxation_solver_module
#endif
end module m_phase_change