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.
# 9 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
 
! For moving immersed boundaries in simulation
# 14 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
# 6 "/home/runner/work/MFC/MFC/src/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"
 
# 207 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 232 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 243 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 245 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 255 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 283 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 293 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 303 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 312 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 329 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 339 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 346 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 352 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 358 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 364 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 370 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 376 "/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"
 
# 192 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 213 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 241 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 256 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 266 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 275 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 291 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 301 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 308 "/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
 
# 21 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 37 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 50 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 76 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 91 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 102 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 115 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 143 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 154 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 165 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 176 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 186 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 197 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 208 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 218 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 224 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 230 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 236 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 242 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 244 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
# 245 "/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.
# 63 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 81 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 88 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 111 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 127 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 153 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 159 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 167 "/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        !< definitions of the derived types
 
    use m_global_parameters    !< definitions of the global parameters
 
    use m_mpi_proxy            !< message passing interface (mpi) module proxy
 
    use m_variables_conversion !< state variables type conversion procedures
 
    use ieee_arithmetic
 
    use m_helper_basic         !< functions to compare floating point numbers
 
    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 water pressure
    real(wp), parameter :: tcr = 385.05_wp + 273.15_wp  !< Critical water temperature
    real(wp), parameter :: mixm = 1.0e-8_wp !< threshold for 'mixture cell'. If Y < mixM, phase change does not happen
    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
    !> @}
 
 
# 48 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if defined(MFC_OpenACC)
# 48 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc declare create(A, B, C, D) 
# 48 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif defined(MFC_OpenMP)
# 48 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$omp declare target (A, B, C, D) 
# 48 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
 
contains
 
    !> This subroutine should dispatch to the correct relaxation solver based
        !!      some parameter. It 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
# 59 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
        call s_mpi_abort("m_phase_change.fpp:59: "// &
# 59 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
                         "Assertion failed: .false.. " &
# 59 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
                         //"s_relaxation_solver called but it currently does nothing")
# 59 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
    end if
    end subroutine s_relaxation_solver
 
    !>  The purpose of this subroutine is to initialize the phase change module
        !!      by setting the parameters needed for phase change and
        !!      selecting the phase change module that will be used
        !!      (pT- or pTg-equilibrium)
    impure subroutine s_initialize_phasechange_module
        ! variables used in the calculation of the saturation curves for fluids 1 and 2
        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
 
    !>  This subroutine is created to activate either the pT- (N fluids) or the
        !!      pTg-equilibrium (2 fluids for g-equilibrium)
        !!      model, also considering mass depletion, depending on the incoming
        !!      state conditions.
        !!  @param q_cons_vf Cell-average conservative variables
    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 temperature for mixture, overheated vapor, and subcooled liquid. Saturation Temperatures at overheated vapor and subcooled liquid
        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]')
# 100 "/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
# 102 "/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
 
        ! starting equilibrium solver
 
# 108 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 108 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if defined(MFC_OpenACC)
# 108 "/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) 
# 108 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif defined(MFC_OpenMP)
# 108 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
    
# 108 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 108 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 108 "/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) 
# 108 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
# 108 "/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
 
# 114 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if defined(MFC_OpenACC)
# 114 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc loop seq 
# 114 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif defined(MFC_OpenMP)
# 114 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 114 "/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 + contxb - 1)%sf(j, k, l)
 
                        ! Total Volume Fraction
                        tvf = tvf + q_cons_vf(i + advxb - 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 + contxb - 1)%sf(j, k, l) + q_cons_vf(vp + contxb - 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 + contxb - 1)%sf(j, k, l)
 
                    m2 = q_cons_vf(vp + contxb - 1)%sf(j, k, l)
 
                    ! kinetic energy as an auxiliary variable to the calculation of the total internal energy
                    dyne = 0.0_wp
 
# 140 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if defined(MFC_OpenACC)
# 140 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc loop seq 
# 140 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif defined(MFC_OpenMP)
# 140 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 140 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
                    do i = momxb, momxe
 
                        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(e_idx)%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 is required
                    ! NOTE that NOTHING else needs to be updated OTHER than the individual partial densities
                    ! given the outputs from the pT- and pTg-equilibrium solvers are just p and one of the partial masses
                    ! (pTg- case)
                    if ((relax_model == 6) .and. ((q_cons_vf(lp + contxb - 1)%sf(j, k, l) > mixm*rm) &
                                                  .and. (q_cons_vf(vp + contxb - 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 + contxb - 1)%sf(j, k, l) = mixm*rm
 
                        ! transferring the total mass to vapor
                        q_cons_vf(vp + contxb - 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 + contxb - 1)%sf(j, k, l) = (1.0_wp - mixm)*rm
 
                        ! depleting the mass of vapor
                        q_cons_vf(vp + contxb - 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 + contxb - 1)%sf(j, k, l) = mixm*rm
 
                            ! correcting the vapor partial density
                            q_cons_vf(vp + contxb - 1)%sf(j, k, l) = (1.0_wp - mixm)*rm
 
                        elseif (tssl < tsatsl) then
 
                            ! Assigning pressure
                            ps = pssl
 
                            ! Assigning Temperature
                            ts = tssl
 
                            ! correcting the liquid partial density
                            q_cons_vf(lp + contxb - 1)%sf(j, k, l) = (1.0_wp - mixm)*rm
 
                            ! correcting the vapor partial density
                            q_cons_vf(vp + contxb - 1)%sf(j, k, l) = mixm*rm
 
                        else
 
                            ! returning partial pressures to what they were from the homogeneous solver
                            ! liquid
                            q_cons_vf(lp + contxb - 1)%sf(j, k, l) = m1
 
                            ! vapor
                            q_cons_vf(vp + contxb - 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
 
 
# 240 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if defined(MFC_OpenACC)
# 240 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc loop seq 
# 240 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif defined(MFC_OpenMP)
# 240 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 240 "/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
 
# 265 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if defined(MFC_OpenACC)
# 265 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc loop seq 
# 265 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif defined(MFC_OpenMP)
# 265 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 265 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
                    do i = 1, num_fluids
 
                        ! volume fractions
                        q_cons_vf(i + advxb - 1)%sf(j, k, l) = q_cons_vf(i + contxb - 1)%sf(j, k, l)/rhok(i)
 
                        ! alpha*rho*e
                        if (model_eqns == 3) then
                            q_cons_vf(i + intxb - 1)%sf(j, k, l) = q_cons_vf(i + contxb - 1)%sf(j, k, l)*ek(i)
                        end if
 
                        ! Total entropy
                        rhos = rhos + q_cons_vf(i + contxb - 1)%sf(j, k, l)*sk(i)
 
                    end do
                end do
            end do
        end do
 
# 283 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 283 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if defined(MFC_OpenACC)
# 283 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc end parallel loop
# 283 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif defined(MFC_OpenMP)
# 283 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
    
# 283 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 283 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$omp end target teams loop
# 283 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
# 283 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
 
    end subroutine s_infinite_relaxation_k
 
    !>  This auxiliary subroutine is created to activate the pT-equilibrium for N fluids
        !!  @param j generic loop iterator for x direction
        !!  @param k generic loop iterator for y direction
        !!  @param l generic loop iterator for z direction
        !!  @param MFL flag that tells whether the fluid is gas (0), liquid (1), or a mixture (2)
        !!  @param pS equilibrium pressure at the interface
        !!  @param p_infpT stiffness for the participating fluids under pT-equilibrium
        !!  @param q_cons_vf Cell-average conservative variables
        !!  @param rhoe mixture energy
        !!  @param TS equilibrium temperature at the interface
    subroutine s_infinite_pt_relaxation_k(j, k, l, MFL, pS, p_infpT, q_cons_vf, rhoe, TS)
 
# 298 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#ifdef _CRAYFTN
# 298 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!DIR$ INLINEALWAYS s_infinite_pt_relaxation_k
# 298 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif MFC_OpenACC
# 298 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc routine seq 
# 298 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif MFC_OpenMP
# 298 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
    
# 298 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 298 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$omp declare target device_type(any) 
# 298 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
# 300 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
        ! 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
 
# 315 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if defined(MFC_OpenACC)
# 315 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc loop seq 
# 315 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif defined(MFC_OpenMP)
# 315 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 315 "/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
 
# 321 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if defined(MFC_OpenACC)
# 321 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc loop seq 
# 321 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif defined(MFC_OpenMP)
# 321 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 321 "/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 + contxb - 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 + contxb - 1)%sf(j, k, l)*qvs(i)
 
        end do
 
# 340 "/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
 
# 379 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if defined(MFC_OpenACC)
# 379 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc loop seq 
# 379 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif defined(MFC_OpenMP)
# 379 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 379 "/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 + contxb - 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 + contxb - 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
 
    !>  This auxiliary subroutine is created to activate the pTg-equilibrium for N fluids under pT
        !!      and 2 fluids under pTg-equilibrium. There is a final common p and T during relaxation
        !!  @param j generic loop iterator for x direction
        !!  @param k generic loop iterator for y direction
        !!  @param l generic loop iterator for z direction
        !!  @param pS equilibrium pressure at the interface
        !!  @param p_infpT stiffness for the participating fluids under pT-equilibrium
        !!  @param rhoe mixture energy
        !!  @param q_cons_vf Cell-average conservative variables
        !!  @param TS equilibrium temperature at the interface
    subroutine s_infinite_ptg_relaxation_k(j, k, l, pS, p_infpT, rhoe, q_cons_vf, TS)
 
# 413 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#ifdef _CRAYFTN
# 413 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!DIR$ INLINEALWAYS s_infinite_ptg_relaxation_k
# 413 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif MFC_OpenACC
# 413 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc routine seq 
# 413 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif MFC_OpenMP
# 413 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
    
# 413 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 413 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$omp declare target device_type(any) 
# 413 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
# 415 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
        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
# 425 "/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
# 427 "/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 + contxb - 1)%sf(j, k, l) &
                                   + q_cons_vf(vp + contxb - 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
 
# 473 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#if defined(MFC_OpenACC)
# 473 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc loop seq 
# 473 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif defined(MFC_OpenMP)
# 473 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 473 "/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 + contxb - 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 + contxb - 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 + contxb - 1)%sf(j, k, l) &
                            *cvs(i)*(gs_min(i) - 1)/(ps + ps_inf(i))
 
                    mcvgp2 = mcvgp2 + q_cons_vf(i + contxb - 1)%sf(j, k, l) &
                             *cvs(i)*(gs_min(i) - 1)/((ps + ps_inf(i))**2)
 
                    mqd = mqd + q_cons_vf(i + contxb - 1)%sf(j, k, l)*qvs(i)
 
                    ! sum of the total alpha*rho*cp of the system
                    mcpd = mcpd + q_cons_vf(i + contxb - 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 + contxb - 1)%sf(j, k, l)
 
            ! mass of the two participating fluids
            mt = q_cons_vf(lp + contxb - 1)%sf(j, k, l) &
                 + q_cons_vf(vp + contxb - 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.0_wp*(matmul(invjac, r2d))
 
            ! updating two reacting 'masses'. Recall that inert 'masses' do not change during the phase change
            ! liquid
            q_cons_vf(lp + contxb - 1)%sf(j, k, l) = q_cons_vf(lp + contxb - 1)%sf(j, k, l) + om*deltamp(1)
 
            ! gas
            q_cons_vf(vp + contxb - 1)%sf(j, k, l) = q_cons_vf(vp + contxb - 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 + contxb - 1)%sf(j, k, l)
 
            ! mass of the two participating fluids
            mt = q_cons_vf(lp + contxb - 1)%sf(j, k, l) &
                 + q_cons_vf(vp + contxb - 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
 
    !>  This auxiliary subroutine corrects 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
        !!  @param MCT partial density correction parameter
        !!  @param q_cons_vf Cell-average conservative variables
        !!  @param rM sum of the reacting masses
        !!  @param j generic loop iterator for x direction
        !!  @param k generic loop iterator for y direction
        !!  @param l generic loop iterator for z direction
    subroutine s_correct_partial_densities(MCT, q_cons_vf, rM, j, k, l)
 
# 640 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#ifdef _CRAYFTN
# 640 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!DIR$ INLINEALWAYS s_correct_partial_densities
# 640 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif MFC_OpenACC
# 640 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc routine seq 
# 640 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif MFC_OpenMP
# 640 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
    
# 640 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 640 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$omp declare target device_type(any) 
# 640 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
# 642 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
        !> @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 + contxb - 1)%sf(j, k, l) >= -1.0_wp*mixm) .and. &
                (q_cons_vf(vp + contxb - 1)%sf(j, k, l) >= -1.0_wp*mixm)) then
 
                q_cons_vf(lp + contxb - 1)%sf(j, k, l) = 0.0_wp
 
                q_cons_vf(vp + contxb - 1)%sf(j, k, l) = 0.0_wp
 
                rm = q_cons_vf(lp + contxb - 1)%sf(j, k, l) + q_cons_vf(vp + contxb - 1)%sf(j, k, l)
 
            end if
 
        end if
 
        ! Defining the correction in terms of an absolute value might not be the best practice.
        ! Maybe a good way to do this is to partition the partial densities, giving a small percentage of the total reacting density
        mct = 2*mixm
 
        ! correcting the partial densities of the reacting fluids. What to do for the nonreacting ones?
        if (q_cons_vf(lp + contxb - 1)%sf(j, k, l) < 0.0_wp) then
 
            q_cons_vf(lp + contxb - 1)%sf(j, k, l) = mct*rm
 
            q_cons_vf(vp + contxb - 1)%sf(j, k, l) = (1.0_wp - mct)*rm
 
        elseif (q_cons_vf(vp + contxb - 1)%sf(j, k, l) < 0.0_wp) then
 
            q_cons_vf(lp + contxb - 1)%sf(j, k, l) = (1.0_wp - mct)*rm
 
            q_cons_vf(vp + contxb - 1)%sf(j, k, l) = mct*rm
 
        end if
    end subroutine s_correct_partial_densities
 
    !>  This auxiliary subroutine finds the Saturation temperature for a given
        !!      saturation pressure through a newton solver
        !!  @param pSat Saturation Pressure
        !!  @param TSat Saturation Temperature
        !!  @param TSIn equilibrium Temperature
    elemental subroutine s_tsat(pSat, TSat, TSIn)
 
# 691 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#ifdef _CRAYFTN
# 691 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!DIR$ INLINEALWAYS s_TSat
# 691 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif MFC_OpenACC
# 691 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$acc routine seq 
# 691 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#elif MFC_OpenMP
# 691 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
    
# 691 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
# 691 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
!$omp declare target device_type(any) 
# 691 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
#endif
# 693 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
 
        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
            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
 
            end do
 
        end if
 
    end subroutine s_tsat
 
    !>  This subroutine finalizes the phase change module
    impure subroutine s_finalize_relaxation_solver_module
    end subroutine s_finalize_relaxation_solver_module
 
#endif
 
end module m_phase_change