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"
 
# 167 "/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"
 
# 167 "/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"
 
# 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"
 
# 126 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 156 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 197 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 211 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 236 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 247 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 249 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 260 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 310 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 320 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 330 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 339 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 356 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 366 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 373 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 379 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 385 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 391 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 397 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 403 "/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"
 
# 167 "/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
# 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
 
    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 = 100000            !< max Newton iterations before accepting the last iterate
    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]
    integer, parameter  :: ptg_ls_max = 30              !< max backtracking-line-search halvings in the pTg solver
    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
    !> @}
 
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
# 43 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
        call s_mpi_abort("m_phase_change.fpp:43: " // "Assertion failed: .false.. " &
# 43 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
                         & // "s_relaxation_solver called but it currently does nothing")
# 43 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
    end if
 
    end subroutine s_relaxation_solver
 
    !> Initialize the phase change module (no module-level state to set up; the pT/pTg relaxation solvers are self-contained)
    impure subroutine s_initialize_phasechange_module
 
    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                    !< equilibrium pressure
        real(wp) :: ts                    !< equilibrium temperature
        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
 
# 65 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
            real(wp), dimension(num_fluids) :: p_infpt, sk, hk, gk, ek, rhok
# 67 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
        !> 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
#endif
#endif
 
        ! starting equilibrium solver
 
 
# 80 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 80 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if defined(MFC_OpenACC)
# 80 "/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_infpT, sk, hk, gk, ek, rhok, pS, TS, rhoe, dynE, rhos, rho, rM, m1, m2, MCT, TvF)
# 80 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif defined(MFC_OpenMP)
# 80 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 80 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 80 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 80 "/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) &
# 80 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$omp& private(i, j, k, l, p_infpT, sk, hk, gk, ek, rhok, pS, TS, rhoe, dynE, rhos, rho, rM, m1, m2, MCT, TvF)
# 80 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
# 82 "/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
 
# 86 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if defined(MFC_OpenACC)
# 86 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc loop seq 
# 86 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif defined(MFC_OpenMP)
# 86 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 86 "/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
 
# 110 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if defined(MFC_OpenACC)
# 110 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc loop seq 
# 110 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif defined(MFC_OpenMP)
# 110 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 110 "/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
                        ! Solve pTg-equilibrium directly on the actual reacting masses. The Newton solver projects
                        ! the liquid mass onto [0, mT], so it recovers the single-phase limits itself (ml -> 0 for
                        ! all-vapor, ml -> mT for all-liquid). The former overheated-vapor / subcooled-liquid pT
                        ! shortcuts were removed: their pT states differ O(1) from the pTg equilibrium, so the
                        ! sub-ULP shortcut/pTg branch decision flipped across backends (CPU vs GPU) near a phase
                        ! boundary and destroyed cross-backend reproducibility.
                        q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) = m1
                        q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) = m2
 
                        call s_infinite_ptg_relaxation_k(j, k, l, ps, rhoe, q_cons_vf, ts)
                    end if
 
                    ! Calculations AFTER equilibrium
 
 
# 142 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if defined(MFC_OpenACC)
# 142 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc loop seq 
# 142 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif defined(MFC_OpenMP)
# 142 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 142 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
                    do i = 1, num_fluids
                        ! entropy
                        sk(i) = cvs(i)*log((ts**isentrope_n(i))/((ps + isentrope_b(i))**(isentrope_n(i) - 1.0_wp))) + qvps(i)
 
                        ! enthalpy
                        hk(i) = isentrope_n(i)*cvs(i)*ts + qvs(i)
 
                        ! Gibbs-free energy
                        gk(i) = hk(i) - ts*sk(i)
 
                        ! densities
                        rhok(i) = f_sg_thermal(ps, ts, isentrope_n(i), isentrope_b(i), cvs(i))
 
                        ! internal energy
                        ek(i) = (ps + isentrope_n(i)*isentrope_b(i))/(ps + isentrope_b(i))*cvs(i)*ts + qvs(i)
                    end do
 
                    ! calculating volume fractions, internal energies, and total entropy
                    rhos = 0.0_wp
 
# 162 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if defined(MFC_OpenACC)
# 162 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc loop seq 
# 162 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif defined(MFC_OpenMP)
# 162 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 162 "/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
 
# 179 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if defined(MFC_OpenACC)
# 179 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc end parallel loop
# 179 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif defined(MFC_OpenMP)
# 179 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 179 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$omp end target teams loop
# 179 "/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)
 
 
# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#ifdef _CRAYFTN
# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if MFC_OpenACC
# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc routine seq 
# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif MFC_OpenMP
# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$omp declare target device_type(any) 
# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#else
# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!DIR$ NOINLINE s_infinite_pt_relaxation_k
# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif MFC_OpenACC
# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc routine seq 
# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif MFC_OpenMP
# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$omp declare target device_type(any) 
# 187 "/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
 
# 201 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if defined(MFC_OpenACC)
# 201 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc loop seq 
# 201 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif defined(MFC_OpenMP)
# 201 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 201 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
        do i = 1, num_fluids
            p_infpt(i) = isentrope_b(i)
            p_infpt_sum = p_infpt_sum + abs(p_infpt(i))
        end do
        ! Performing tests before initializing the pT-equilibrium
 
# 207 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if defined(MFC_OpenACC)
# 207 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc loop seq 
# 207 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif defined(MFC_OpenMP)
# 207 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 207 "/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)*isentrope_n(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
 
# 224 "/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(isentrope_B) 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
        ! Relative criterion written in multiply form to avoid dividing by pO (pO = 0 on the first pass).
        do while ((abs(ps - po) > palpha_eps) .and. (abs(ps - po) > (palpha_eps/1.e4_wp)*abs(po)) .or. (ns == 0))
            ! increasing counter
            ns = ns + 1
            ! guard against non-convergence: accept the last iterate rather than looping forever
            if (ns >= max_iter) exit
 
            ! 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
 
# 262 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if defined(MFC_OpenACC)
# 262 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc loop seq 
# 262 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif defined(MFC_OpenMP)
# 262 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 262 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
            do i = 1, num_fluids
                gp = gp + (isentrope_n(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 + (isentrope_n(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
 
    !> Evaluate the pTg-equilibrium residual R2D and temperature TS at a trial state (ml, pS) WITHOUT mutating q_cons_vf, so the
    !! Newton driver can line-search. The total reacting mass mT is conserved, so the reacting masses are (ml, mT - ml) and only the
    !! inert fluids are read from q_cons_vf. Also returns the mixture sums the Jacobian and the final temperature need.
    subroutine s_compute_ptg_residual(ml, mT, pS, j, k, l, q_cons_vf, rhoe, R2D, TS, mCP, mQ, mCVGP, mCVGP2, mCPD)
 
 
# 287 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if MFC_OpenACC
# 287 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc routine seq 
# 287 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif MFC_OpenMP
# 287 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 287 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 287 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$omp declare target device_type(any) 
# 287 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
 
        real(wp), intent(in)                                :: ml, mT, pS, rhoe
        integer, intent(in)                                 :: j, k, l
        type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf
        real(wp), dimension(2), intent(out)                 :: R2D
        real(wp), intent(out)                               :: TS, mCP, mQ, mCVGP, mCVGP2, mCPD
        real(wp)                                            :: mQD
        integer                                             :: i
 
        ! reacting fluids contribute via (ml, mT - ml); inert fluids are summed from q_cons_vf
        mcp = ml*cvs(lp)*isentrope_n(lp) + (mt - ml)*cvs(vp)*isentrope_n(vp)
        mq = ml*qvs(lp) + (mt - ml)*qvs(vp)
        mcvgp = 0.0_wp; mcvgp2 = 0.0_wp; mcpd = 0.0_wp; mqd = 0.0_wp
 
# 301 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if defined(MFC_OpenACC)
# 301 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc loop seq 
# 301 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif defined(MFC_OpenMP)
# 301 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 301 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
        do i = 1, num_fluids
            if ((i /= lp) .and. (i /= vp)) then
                mcp = mcp + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*cvs(i)*isentrope_n(i)
                mq = mq + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*qvs(i)
                mcvgp = mcvgp + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*cvs(i)*(isentrope_n(i) - 1)/(ps + isentrope_b(i))
                mcvgp2 = mcvgp2 + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, &
                                            & l)*cvs(i)*(isentrope_n(i) - 1)/((ps + isentrope_b(i))**2)
                mqd = mqd + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*qvs(i)
                mcpd = mcpd + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*cvs(i)*isentrope_n(i)
            end if
        end do
 
        ts = 1.0_wp/(mt*cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp)) + ml*(cvs(lp)*(isentrope_n(lp) - 1)/(ps &
                     & + isentrope_b(lp)) - cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp))) + mcvgp)
 
        ! (i) Gibbs free-energy equality
        r2d(1) = ts*((cvs(lp)*isentrope_n(lp) - cvs(vp)*isentrope_n(vp))*(1 - log(ts)) - (qvps(lp) - qvps(vp)) + cvs(lp) &
            & *(isentrope_n(lp) - 1)*log(ps + isentrope_b(lp)) - cvs(vp)*(isentrope_n(vp) - 1)*log(ps + isentrope_b(vp))) &
            & + qvs(lp) - qvs(vp)
 
        ! (ii) constant-energy condition
        r2d(2) = rhoe + ps + ml*(qvs(vp) - qvs(lp)) - mt*qvs(vp) - mqd + (ml*(isentrope_n(vp)*cvs(vp) - isentrope_n(lp)*cvs(lp)) &
            & - mt*isentrope_n(vp)*cvs(vp) - mcpd)/(ml*(cvs(lp)*(isentrope_n(lp) - 1)/(ps + isentrope_b(lp)) - cvs(vp) &
            & *(isentrope_n(vp) - 1)/(ps + isentrope_b(vp))) + mt*cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp)) + mcvgp)
 
    end subroutine s_compute_ptg_residual
 
    !> Apply pTg-equilibrium relaxation: a damped (backtracking line search) Newton solve for the reacting liquid mass ml and
    !! pressure pS enforcing Gibbs equality and energy conservation, converging on the residual norm (absolute ptgalpha_eps, or the
    !! rhoe-relative branch). Every step is projected onto the physical bounds 0 <= ml <= mT, pS > pmin. This converges in a handful
    !! of iterations with a bounded, uniform count (no GPU warp divergence), unlike the former fixed 1e-3 underrelaxation that
    !! stalled far from the root.
    subroutine s_infinite_ptg_relaxation_k(j, k, l, pS, rhoe, q_cons_vf, TS)
 
 
# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#ifdef _CRAYFTN
# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if MFC_OpenACC
# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc routine seq 
# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif MFC_OpenMP
# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$omp declare target device_type(any) 
# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#else
# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!DIR$ NOINLINE s_infinite_ptg_relaxation_k
# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif MFC_OpenACC
# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc routine seq 
# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif MFC_OpenMP
# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$omp declare target device_type(any) 
# 336 "/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), intent(in)                                   :: rhoe
        type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
        real(wp), intent(inout)                                :: TS
        real(wp), dimension(2, 2)                              :: Jac, InvJac
        real(wp), dimension(2)                                 :: R2D, R2D_try, DeltamP
        real(wp)                                               :: mCP, mCPD, mCVGP, mCVGP2, mQ
        real(wp)                                               :: ml, ml_try, mT, pS_try, pmin, lambda, resnorm, resnorm_try
        real(wp)                                               :: dFdT, dTdm, dTdp, detJ
        integer                                                :: ns, ls
 
        ! total reacting mass is conserved; the liquid mass ml is the primary unknown, vapor mass = mT - ml
        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)
        ml = q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l)
 
        ! recover a physical pressure guess when the incoming pS is non-physical
        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 - isentrope_n(lp)*isentrope_b(lp)/(isentrope_n(lp) - 1))/qvs(lp)))) .or. ((ps >= 0.0_wp) &
            & .and. (ps < 1.0e-1_wp))) then
            ps = 1.0e4_wp
        end if
 
        ! pressure floor (stiffened gas requires pS + isentrope_B > 0 for both phases)
        pmin = -min(isentrope_b(lp), isentrope_b(vp)) + 1.0_wp
 
        call s_compute_ptg_residual(ml, mt, ps, j, k, l, q_cons_vf, rhoe, r2d, ts, mcp, mq, mcvgp, mcvgp2, mcpd)
        resnorm = sqrt(r2d(1)**2 + r2d(2)**2)
 
        do ns = 1, max_iter
            ! converged on the absolute residual, or on the rhoe-relative residual (multiply form, rhoe > 0)
            if ((resnorm <= ptgalpha_eps) .or. (resnorm <= (ptgalpha_eps/1.e6_wp)*rhoe)) exit
 
            ! 2x2 Jacobian of (Gibbs equality, energy) with respect to (ml, pS) at the current state
            dfdt = -(cvs(lp)*isentrope_n(lp) - cvs(vp)*isentrope_n(vp))*log(ts) - (qvps(lp) - qvps(vp)) + cvs(lp)*(isentrope_n(lp) &
                     & - 1)*log(ps + isentrope_b(lp)) - cvs(vp)*(isentrope_n(vp) - 1)*log(ps + isentrope_b(vp))
            dtdm = -(cvs(lp)*(isentrope_n(lp) - 1)/(ps + isentrope_b(lp)) - cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp))) &
                     & *ts**2
            dtdp = (mt*cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp))**2 + ml*(cvs(lp)*(isentrope_n(lp) - 1)/(ps &
                    & + isentrope_b(lp))**2 - cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp))**2) + mcvgp2)*ts**2
 
            jac(1, 1) = dfdt*dtdm
            jac(1, &
                & 2) = dfdt*dtdp + ts*(cvs(lp)*(isentrope_n(lp) - 1)/(ps + isentrope_b(lp)) - cvs(vp)*(isentrope_n(vp) - 1)/(ps &
                & + isentrope_b(vp)))
            jac(2, &
                & 1) = qvs(vp) - qvs(lp) + (cvs(vp)*isentrope_n(vp) - cvs(lp)*isentrope_n(lp))/(ml*(cvs(lp)*(isentrope_n(lp) - 1) &
                & /(ps + isentrope_b(lp)) - cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp))) + mt*cvs(vp)*(isentrope_n(vp) &
                & - 1)/(ps + isentrope_b(vp)) + mcvgp) - (ml*(cvs(vp)*isentrope_n(vp) - cvs(lp)*isentrope_n(lp)) - mt*cvs(vp) &
                & *isentrope_n(vp) - mcpd)*(cvs(lp)*(isentrope_n(lp) - 1)/(ps + isentrope_b(lp)) - cvs(vp)*(isentrope_n(vp) - 1) &
                & /(ps + isentrope_b(vp)))/((ml*(cvs(lp)*(isentrope_n(lp) - 1)/(ps + isentrope_b(lp)) - cvs(vp)*(isentrope_n(vp) &
                & - 1)/(ps + isentrope_b(vp))) + mt*cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp)) + mcvgp)**2)
            jac(2, &
                & 2) = 1 + (ml*(cvs(vp)*isentrope_n(vp) - cvs(lp)*isentrope_n(lp)) - mt*cvs(vp)*isentrope_n(vp) - mcpd) &
                & *(ml*(cvs(lp)*(isentrope_n(lp) - 1)/(ps + isentrope_b(lp))**2 - cvs(vp)*(isentrope_n(vp) - 1)/(ps &
                & + isentrope_b(vp))**2) + mt*cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp))**2 + mcvgp2)/(ml*(cvs(lp) &
                & *(isentrope_n(lp) - 1)/(ps + isentrope_b(lp)) - cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp))) &
                & + mt*cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp)) + mcvgp)**2
 
            detj = jac(1, 1)*jac(2, 2) - jac(1, 2)*jac(2, 1)
            ! singular Jacobian: no usable Newton direction, accept the current (best) state
            if (detj == 0.0_wp) exit
 
            invjac(1, 1) = jac(2, 2)/detj
            invjac(1, 2) = -jac(1, 2)/detj
            invjac(2, 1) = -jac(2, 1)/detj
            invjac(2, 2) = jac(1, 1)/detj
 
            deltamp(1) = -(invjac(1, 1)*r2d(1) + invjac(1, 2)*r2d(2))
            deltamp(2) = -(invjac(2, 1)*r2d(1) + invjac(2, 2)*r2d(2))
 
            ! backtracking line search: halve the step until the residual decreases, keeping the state
            ! physical (0 <= ml <= mT, pS above the stiffened-gas floor)
            lambda = 1.0_wp
            do ls = 1, ptg_ls_max
                ml_try = min(max(ml + lambda*deltamp(1), 0.0_wp), mt)
                ps_try = max(ps + lambda*deltamp(2), pmin)
                call s_compute_ptg_residual(ml_try, mt, ps_try, j, k, l, q_cons_vf, rhoe, r2d_try, ts, mcp, mq, mcvgp, mcvgp2, mcpd)
                resnorm_try = sqrt(r2d_try(1)**2 + r2d_try(2)**2)
                if ((resnorm_try < resnorm) .or. (ls == ptg_ls_max)) exit
                lambda = 0.5_wp*lambda
            end do
 
            ! accept the trial state (TS, mCP, mQ, mCVGP, mCVGP2, mCPD already set to it by the last call)
            ml = ml_try; ps = ps_try; r2d = r2d_try; resnorm = resnorm_try
        end do
 
        ! commit the reacting masses (mT conserved) and set the common temperature
        q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) = ml
        q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) = mt - ml
 
        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)
 
 
# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#ifdef _CRAYFTN
# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if MFC_OpenACC
# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc routine seq 
# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif MFC_OpenMP
# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$omp declare target device_type(any) 
# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#else
# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!DIR$ NOINLINE s_correct_partial_densities
# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif MFC_OpenACC
# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc routine seq 
# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif MFC_OpenMP
# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$omp declare target device_type(any) 
# 438 "/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
 
    !> Finalize the phase change module
    impure subroutine s_finalize_relaxation_solver_module
 
    end subroutine s_finalize_relaxation_solver_module
 
end module m_phase_change