From 69fcc7d99d6a1ed55eeefcda8eba685a590e70af Mon Sep 17 00:00:00 2001
From: Fahad Nabid <55717325+fahnab666@users.noreply.github.com>
Date: Mon, 20 Jul 2026 19:24:59 -0400
Subject: [PATCH 1/6] PR1: add central thermodynamics interface, relocate
 pressure query

Introduce m_thermodynamics as the single equation-of-state interface and move s_compute_pressure into it unchanged. m_variables_conversion re-exports the routine, so all ten call sites are untouched. The leaf module uses only global parameters and the thermochemistry surface, so no module cycle forms. The stiffened-gas and Pyrometheus branches are preserved verbatim.

No changed answers: the branch-representative golden subset (stiffened gas, MHD, bubbles, hypoelastic, phase change, chemistry) passes with zero golden regeneration across all three targets.

First increment of PR1. Sound speed and the remaining query sites follow.
---
 docs/module_categories.json           |  1 +
 src/common/m_thermodynamics.fpp       | 95 +++++++++++++++++++++++++++
 src/common/m_variables_conversion.fpp | 69 +------------------
 3 files changed, 98 insertions(+), 67 deletions(-)
 create mode 100644 src/common/m_thermodynamics.fpp

diff --git a/docs/module_categories.json b/docs/module_categories.json
index 1391908cad..bc8550454a 100644
--- a/docs/module_categories.json
+++ b/docs/module_categories.json
@@ -13,6 +13,7 @@
             "m_riemann_solver_hllc",
             "m_muscl",
             "m_variables_conversion",
+            "m_thermodynamics",
             "m_thinc"
         ]
     },
diff --git a/src/common/m_thermodynamics.fpp b/src/common/m_thermodynamics.fpp
new file mode 100644
index 0000000000..7c6b92c223
--- /dev/null
+++ b/src/common/m_thermodynamics.fpp
@@ -0,0 +1,95 @@
+!>
+!! @file m_thermodynamics.fpp
+!! @brief Contains module m_thermodynamics
+
+#:include 'macros.fpp'
+#:include 'case.fpp'
+
+!> @brief Central equation-of-state interface. Every genuine thermodynamic query (pressure here; energy, temperature, sound speed,
+!! and derivatives in later PRs) is answered through one module so no solver path carries its own EOS algebra. Compile-time
+!! chemistry selects the adapter: the stiffened-gas analytic branch or the Pyrometheus temperature-then-query branch. This is a leaf
+!! module (it uses global parameters and the thermochemistry surface, never m_variables_conversion) so no cycle forms.
+module m_thermodynamics
+
+    use m_derived_types
+    use m_global_parameters
+    use m_constants, only: model_eqns_4eq
+    use m_thermochem, only: num_species, get_temperature, get_pressure
+
+    implicit none
+
+    private
+    public :: s_compute_pressure
+
+contains
+
+    !> Pressure from the conserved energy and mixture properties. The #:if chemistry split selects the stiffened-gas adapter
+    !! (analytic inverse of the mixture EOS) or the Pyrometheus adapter (solve temperature, then query pressure). Arithmetic is
+    !! relocated unchanged from the conversion module.
+    subroutine s_compute_pressure(energy, alf, dyn_p, pi_inf, gamma, rho, qv, rhoYks, pres, T, stress, mom, G, pres_mag)
+
+        $:GPU_ROUTINE(function_name='s_compute_pressure',parallelism='[seq]', cray_noinline=True)
+
+        real(stp), intent(in)           :: energy, alf
+        real(wp), intent(in)            :: dyn_p
+        real(wp), intent(in)            :: pi_inf, gamma, rho, qv
+        real(wp), intent(out)           :: pres
+        real(wp), intent(inout)         :: T
+        real(stp), intent(in), optional :: stress, mom
+        real(wp), intent(in), optional  :: G, pres_mag
+
+        ! Chemistry
+        real(wp), dimension(1:num_species), intent(in) :: rhoYks
+        real(wp), dimension(1:num_species)             :: Y_rs
+        real(wp)                                       :: E_e
+        real(wp)                                       :: e_Per_Kg, Pdyn_Per_Kg
+        real(wp)                                       :: T_guess
+        integer                                        :: s  !< Generic loop iterator
+        #:if not chemistry
+            ! Depending on model_eqns and bubbles_euler, the appropriate procedure for computing pressure is targeted by the
+            ! procedure pointer
+
+            if (mhd) then
+                ! MHD pressure: subtract magnetic pressure from total energy
+                pres = (energy - dyn_p - pi_inf - qv - pres_mag)/gamma
+            else if ((model_eqns /= model_eqns_4eq) .and. (bubbles_euler .neqv. .true.)) then
+                ! Gamma/pi_inf model or five-equation model (Allaire et al. JCP 2002): p from mixture EOS
+                pres = (energy - dyn_p - pi_inf - qv)/gamma
+            else if ((model_eqns /= model_eqns_4eq) .and. bubbles_euler) then
+                ! Bubble-augmented pressure with void fraction correction
+                pres = ((energy - dyn_p)/(1._wp - alf) - pi_inf - qv)/gamma
+            else
+                ! Four-equation model (Kapila et al. PoF 2001): Tait EOS inversion
+                pres = (pref + pi_inf)*(energy/(rhoref*(1 - alf)))**(1/gamma + 1) - pi_inf
+            end if
+
+            if (hypoelasticity .and. present(G)) then
+                ! Subtract elastic strain energy before computing pressure (hypoelastic model)
+                E_e = 0._wp
+                do s = eqn_idx%stress%beg, eqn_idx%stress%end
+                    if (G > 0) then
+                        E_e = E_e + ((stress/rho)**2._wp)/(4._wp*G)
+                        ! Double for shear stresses
+                        if (any(s == shear_indices)) then
+                            E_e = E_e + ((stress/rho)**2._wp)/(4._wp*G)
+                        end if
+                    end if
+                end do
+
+                pres = (energy - 0.5_wp*(mom**2._wp)/rho - pi_inf - qv - E_e)/gamma
+            end if
+        #:else
+            ! Reacting mixture pressure from temperature and species
+            Y_rs(:) = rhoYks(:)/rho
+            e_Per_Kg = energy/rho
+            Pdyn_Per_Kg = dyn_p/rho
+
+            T_guess = T
+
+            call get_temperature(e_Per_Kg - Pdyn_Per_Kg, T_guess, Y_rs, .true., T)
+            call get_pressure(rho, T, Y_rs, pres)
+        #:endif
+
+    end subroutine s_compute_pressure
+
+end module m_thermodynamics
diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp
index fe17b89c27..e13cba6818 100644
--- a/src/common/m_variables_conversion.fpp
+++ b/src/common/m_variables_conversion.fpp
@@ -18,6 +18,8 @@ module m_variables_conversion
     use m_thermochem, only: num_species, get_temperature, get_pressure, gas_constant, get_mixture_molecular_weight, &
         & get_mixture_energy_mass
 
+    use m_thermodynamics, only: s_compute_pressure
+
     implicit none
 
     private
@@ -79,73 +81,6 @@ contains
 
     end subroutine s_convert_to_mixture_variables
 
-    !> Compute the pressure from the appropriate equation of state
-    subroutine s_compute_pressure(energy, alf, dyn_p, pi_inf, gamma, rho, qv, rhoYks, pres, T, stress, mom, G, pres_mag)
-
-        $:GPU_ROUTINE(function_name='s_compute_pressure',parallelism='[seq]', cray_noinline=True)
-
-        real(stp), intent(in)           :: energy, alf
-        real(wp), intent(in)            :: dyn_p
-        real(wp), intent(in)            :: pi_inf, gamma, rho, qv
-        real(wp), intent(out)           :: pres
-        real(wp), intent(inout)         :: T
-        real(stp), intent(in), optional :: stress, mom
-        real(wp), intent(in), optional  :: G, pres_mag
-
-        ! Chemistry
-        real(wp), dimension(1:num_species), intent(in) :: rhoYks
-        real(wp), dimension(1:num_species)             :: Y_rs
-        real(wp)                                       :: E_e
-        real(wp)                                       :: e_Per_Kg, Pdyn_Per_Kg
-        real(wp)                                       :: T_guess
-        integer                                        :: s  !< Generic loop iterator
-        #:if not chemistry
-            ! Depending on model_eqns and bubbles_euler, the appropriate procedure for computing pressure is targeted by the
-            ! procedure pointer
-
-            if (mhd) then
-                ! MHD pressure: subtract magnetic pressure from total energy
-                pres = (energy - dyn_p - pi_inf - qv - pres_mag)/gamma
-            else if ((model_eqns /= model_eqns_4eq) .and. (bubbles_euler .neqv. .true.)) then
-                ! Gamma/pi_inf model or five-equation model (Allaire et al. JCP 2002): p from mixture EOS
-                pres = (energy - dyn_p - pi_inf - qv)/gamma
-            else if ((model_eqns /= model_eqns_4eq) .and. bubbles_euler) then
-                ! Bubble-augmented pressure with void fraction correction
-                pres = ((energy - dyn_p)/(1._wp - alf) - pi_inf - qv)/gamma
-            else
-                ! Four-equation model (Kapila et al. PoF 2001): Tait EOS inversion
-                pres = (pref + pi_inf)*(energy/(rhoref*(1 - alf)))**(1/gamma + 1) - pi_inf
-            end if
-
-            if (hypoelasticity .and. present(G)) then
-                ! Subtract elastic strain energy before computing pressure (hypoelastic model)
-                E_e = 0._wp
-                do s = eqn_idx%stress%beg, eqn_idx%stress%end
-                    if (G > 0) then
-                        E_e = E_e + ((stress/rho)**2._wp)/(4._wp*G)
-                        ! Double for shear stresses
-                        if (any(s == shear_indices)) then
-                            E_e = E_e + ((stress/rho)**2._wp)/(4._wp*G)
-                        end if
-                    end if
-                end do
-
-                pres = (energy - 0.5_wp*(mom**2._wp)/rho - pi_inf - qv - E_e)/gamma
-            end if
-        #:else
-            ! Reacting mixture pressure from temperature and species
-            Y_rs(:) = rhoYks(:)/rho
-            e_Per_Kg = energy/rho
-            Pdyn_Per_Kg = dyn_p/rho
-
-            T_guess = T
-
-            call get_temperature(e_Per_Kg - Pdyn_Per_Kg, T_guess, Y_rs, .true., T)
-            call get_pressure(rho, T, Y_rs, pres)
-        #:endif
-
-    end subroutine s_compute_pressure
-
     !> Convert mixture variables to density, gamma, pi_inf, and qv for the gamma/pi_inf model. Given conservative or primitive
     !! variables, transfers the density, specific heat ratio function and the liquid stiffness function from q_vf to rho, gamma and
     !! pi_inf.

From 82214f290a56e6e693f54df6c5bcc9d120c790f4 Mon Sep 17 00:00:00 2001
From: Fahad Nabid <55717325+fahnab666@users.noreply.github.com>
Date: Mon, 20 Jul 2026 19:57:56 -0400
Subject: [PATCH 2/6] PR1: route sound speed through the thermodynamics
 interface

Move s_compute_speed_of_sound into m_thermodynamics unchanged and re-export it through m_variables_conversion, matching the pressure op. To keep the interface a leaf module, relocate the pre and post process declaration of the per-fluid property tables (gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps) from m_variables_conversion into m_global_parameters_common; the simulation target already keeps them in global parameters. Allocation, population, and deallocation are unchanged.

No changed answers: the full golden suite passes with zero regeneration (616 of 616 runnable cases) across all three targets, including the sound-speed branches (Wood alt_soundspeed, six-equation, bubble, relativistic, chemistry).

Second increment of PR1. Internal energy and temperature query sites follow.
---
 src/common/m_global_parameters_common.fpp |  9 +++
 src/common/m_thermodynamics.fpp           | 68 ++++++++++++++++++++-
 src/common/m_variables_conversion.fpp     | 73 ++---------------------
 3 files changed, 80 insertions(+), 70 deletions(-)

diff --git a/src/common/m_global_parameters_common.fpp b/src/common/m_global_parameters_common.fpp
index 13e84c8b79..5aaae06f21 100644
--- a/src/common/m_global_parameters_common.fpp
+++ b/src/common/m_global_parameters_common.fpp
@@ -123,6 +123,15 @@ module m_global_parameters_common
     integer                 :: mpi_info_int
     !> @}
 
+    !> @name Per-fluid property tables. In the simulation target these live in m_global_parameters; for pre and post process they
+    !! live here so the thermodynamics interface can read them from global parameters without a module cycle.
+    !> @{
+#ifndef MFC_SIMULATION
+    real(wp), allocatable, public, dimension(:) :: gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps
+    $:GPU_DECLARE(create='[gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps]')
+#endif
+    !> @}
+
 contains
 
     !> Initialize equation-index state (eqn_idx, sys_size, b_size, tensor_size) from the namelist parameters. This is the shared
diff --git a/src/common/m_thermodynamics.fpp b/src/common/m_thermodynamics.fpp
index 7c6b92c223..6882d371ec 100644
--- a/src/common/m_thermodynamics.fpp
+++ b/src/common/m_thermodynamics.fpp
@@ -13,13 +13,16 @@ module m_thermodynamics
 
     use m_derived_types
     use m_global_parameters
-    use m_constants, only: model_eqns_4eq
+    use m_constants, only: model_eqns_4eq, model_eqns_5eq, model_eqns_6eq, avg_state_roe, verysmall, sgm_eps
     use m_thermochem, only: num_species, get_temperature, get_pressure
 
     implicit none
 
     private
     public :: s_compute_pressure
+#ifndef MFC_PRE_PROCESS
+    public :: s_compute_speed_of_sound
+#endif
 
 contains
 
@@ -92,4 +95,67 @@ contains
 
     end subroutine s_compute_pressure
 
+#ifndef MFC_PRE_PROCESS
+    !> Speed of sound from thermodynamic state, supporting the several equation-of-state models. Arithmetic is relocated unchanged
+    !! from the conversion module: the non-chemistry branches hold c squared until a single final sqrt, the chemistry branch returns
+    !! c already square-rooted, and the mixture_err floor is preserved.
+    subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv, vel_sum, c_c, c, qv)
+
+        $:GPU_ROUTINE(parallelism='[seq]')
+
+        real(wp), intent(in) :: pres
+        real(wp), intent(in) :: rho, gamma, pi_inf, qv
+        real(wp), intent(in) :: H
+        #:if not MFC_CASE_OPTIMIZATION and USING_AMD
+            real(wp), dimension(3), intent(in) :: adv
+        #:else
+            real(wp), dimension(num_fluids), intent(in) :: adv
+        #:endif
+        real(wp), intent(in)  :: vel_sum
+        real(wp), intent(in)  :: c_c
+        real(wp), intent(out) :: c
+        real(wp)              :: blkmod1, blkmod2
+        integer               :: q
+
+        if (chemistry) then  ! Reacting mixture sound speed
+            if (avg_state == avg_state_roe .and. abs(c_c) > verysmall) then
+                c = sqrt(c_c - (gamma - 1.0_wp)*(vel_sum - H))
+            else
+                c = sqrt((1.0_wp + 1.0_wp/gamma)*pres/rho)
+            end if
+        else if (relativity) then  ! Relativistic sound speed
+            c = sqrt((1._wp + 1._wp/gamma)*pres/rho/H)
+        else
+            if (alt_soundspeed) then  ! Wood's mixture sound speed via bulk moduli
+                blkmod1 = ((gammas(1) + 1._wp)*pres + pi_infs(1))/gammas(1)
+                blkmod2 = ((gammas(2) + 1._wp)*pres + pi_infs(2))/gammas(2)
+                c = (1._wp/(rho*(adv(1)/blkmod1 + adv(2)/blkmod2)))
+            else if (model_eqns == model_eqns_6eq) then  ! Six-equation model sound speed
+                c = 0._wp
+                $:GPU_LOOP(parallelism='[seq]')
+                do q = 1, num_fluids
+                    c = c + adv(q)*gs_min(q)*(pres + pi_infs(q)/(gammas(q) + 1._wp))
+                end do
+                c = c/rho
+            else if (((model_eqns == model_eqns_4eq) .or. (model_eqns == model_eqns_5eq .and. bubbles_euler))) then
+                ! Sound speed for bubble mixture to order O(\alpha)
+
+                if (mpp_lim .and. (num_fluids > 1)) then
+                    c = (1._wp/gamma + 1._wp)*(pres + pi_inf/(gamma + 1._wp))/rho
+                else
+                    c = (1._wp/gamma + 1._wp)*(pres + pi_inf/(gamma + 1._wp))/(rho*(1._wp - adv(num_fluids)))
+                end if
+            else
+                c = (H - 5.e-1*vel_sum - qv/rho)/gamma
+            end if
+
+            if (mixture_err .and. c < 0._wp) then
+                c = 100._wp*sgm_eps
+            else
+                c = sqrt(c)
+            end if
+        end if
+
+    end subroutine s_compute_speed_of_sound
+#endif
 end module m_thermodynamics
diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp
index e13cba6818..4a3f2e0c6b 100644
--- a/src/common/m_variables_conversion.fpp
+++ b/src/common/m_variables_conversion.fpp
@@ -14,11 +14,14 @@ module m_variables_conversion
     use m_helper_basic
     use m_helper
     use m_constants, only: riemann_solver_hll, riemann_solver_hlld, model_eqns_gamma_law, model_eqns_5eq, model_eqns_6eq, &
-        & model_eqns_4eq, avg_state_roe
+        & model_eqns_4eq
     use m_thermochem, only: num_species, get_temperature, get_pressure, gas_constant, get_mixture_molecular_weight, &
         & get_mixture_energy_mass
 
     use m_thermodynamics, only: s_compute_pressure
+#ifndef MFC_PRE_PROCESS
+    use m_thermodynamics, only: s_compute_speed_of_sound
+#endif
 
     implicit none
 
@@ -41,12 +44,6 @@ module m_variables_conversion
 #endif
     s_finalize_variables_conversion_module
 
-    ! In simulation, gammas, pi_infs, and qvs are already declared in m_global_variables
-#ifndef MFC_SIMULATION
-    real(wp), allocatable, public, dimension(:) :: gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps
-    $:GPU_DECLARE(create='[gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps]')
-#endif
-
     real(wp), allocatable, dimension(:)   :: Gs_vc
     integer, allocatable, dimension(:)    :: bubrs_vc
     real(wp), allocatable, dimension(:,:) :: Res_vc
@@ -1184,68 +1181,6 @@ contains
 
     end subroutine s_finalize_variables_conversion_module
 
-#ifndef MFC_PRE_PROCESS
-    !> Compute the speed of sound from thermodynamic state variables, supporting multiple equation-of-state models.
-    subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv, vel_sum, c_c, c, qv)
-
-        $:GPU_ROUTINE(parallelism='[seq]')
-
-        real(wp), intent(in) :: pres
-        real(wp), intent(in) :: rho, gamma, pi_inf, qv
-        real(wp), intent(in) :: H
-        #:if not MFC_CASE_OPTIMIZATION and USING_AMD
-            real(wp), dimension(3), intent(in) :: adv
-        #:else
-            real(wp), dimension(num_fluids), intent(in) :: adv
-        #:endif
-        real(wp), intent(in)  :: vel_sum
-        real(wp), intent(in)  :: c_c
-        real(wp), intent(out) :: c
-        real(wp)              :: blkmod1, blkmod2
-        integer               :: q
-
-        if (chemistry) then  ! Reacting mixture sound speed
-            if (avg_state == avg_state_roe .and. abs(c_c) > verysmall) then
-                c = sqrt(c_c - (gamma - 1.0_wp)*(vel_sum - H))
-            else
-                c = sqrt((1.0_wp + 1.0_wp/gamma)*pres/rho)
-            end if
-        else if (relativity) then  ! Relativistic sound speed
-            c = sqrt((1._wp + 1._wp/gamma)*pres/rho/H)
-        else
-            if (alt_soundspeed) then  ! Wood's mixture sound speed via bulk moduli
-                blkmod1 = ((gammas(1) + 1._wp)*pres + pi_infs(1))/gammas(1)
-                blkmod2 = ((gammas(2) + 1._wp)*pres + pi_infs(2))/gammas(2)
-                c = (1._wp/(rho*(adv(1)/blkmod1 + adv(2)/blkmod2)))
-            else if (model_eqns == model_eqns_6eq) then  ! Six-equation model sound speed
-                c = 0._wp
-                $:GPU_LOOP(parallelism='[seq]')
-                do q = 1, num_fluids
-                    c = c + adv(q)*gs_min(q)*(pres + pi_infs(q)/(gammas(q) + 1._wp))
-                end do
-                c = c/rho
-            else if (((model_eqns == model_eqns_4eq) .or. (model_eqns == model_eqns_5eq .and. bubbles_euler))) then
-                ! Sound speed for bubble mixture to order O(\alpha)
-
-                if (mpp_lim .and. (num_fluids > 1)) then
-                    c = (1._wp/gamma + 1._wp)*(pres + pi_inf/(gamma + 1._wp))/rho
-                else
-                    c = (1._wp/gamma + 1._wp)*(pres + pi_inf/(gamma + 1._wp))/(rho*(1._wp - adv(num_fluids)))
-                end if
-            else
-                c = (H - 5.e-1*vel_sum - qv/rho)/gamma
-            end if
-
-            if (mixture_err .and. c < 0._wp) then
-                c = 100._wp*sgm_eps
-            else
-                c = sqrt(c)
-            end if
-        end if
-
-    end subroutine s_compute_speed_of_sound
-#endif
-
 #ifndef MFC_PRE_PROCESS
     !> Compute the fast magnetosonic wave speed from the sound speed, density, and magnetic field components.
     subroutine s_compute_fast_magnetosonic_speed(rho, c, B, norm, c_fast, h)

From fb226c50cc02127b3f43282b7c90c9055901cafd Mon Sep 17 00:00:00 2001
From: Fahad Nabid <55717325+fahnab666@users.noreply.github.com>
Date: Mon, 20 Jul 2026 20:30:40 -0400
Subject: [PATCH 3/6] Route primitive-to-conserved energy through the
 thermodynamics interface

Add s_compute_internal_energy to m_thermodynamics as the inverse of
s_compute_pressure over the same model branches (MHD, five-equation,
bubble, four-equation, chemistry). The primitive-to-conservative energy
build in m_variables_conversion now delegates to it, so the conversion
routine no longer carries its own EOS energy algebra. Arithmetic and the
original branch comments are relocated unchanged; the full golden suite
is bit-identical (616 passed, 0 failed).
---
 src/common/m_thermodynamics.fpp       | 42 ++++++++++++++++++++++++-
 src/common/m_variables_conversion.fpp | 44 +++++++++------------------
 2 files changed, 55 insertions(+), 31 deletions(-)

diff --git a/src/common/m_thermodynamics.fpp b/src/common/m_thermodynamics.fpp
index 6882d371ec..2dc133df12 100644
--- a/src/common/m_thermodynamics.fpp
+++ b/src/common/m_thermodynamics.fpp
@@ -14,12 +14,14 @@ module m_thermodynamics
     use m_derived_types
     use m_global_parameters
     use m_constants, only: model_eqns_4eq, model_eqns_5eq, model_eqns_6eq, avg_state_roe, verysmall, sgm_eps
-    use m_thermochem, only: num_species, get_temperature, get_pressure
+    use m_thermochem, only: num_species, get_temperature, get_pressure, gas_constant, get_mixture_molecular_weight, &
+        & get_mixture_energy_mass
 
     implicit none
 
     private
     public :: s_compute_pressure
+    public :: s_compute_internal_energy
 #ifndef MFC_PRE_PROCESS
     public :: s_compute_speed_of_sound
 #endif
@@ -95,6 +97,44 @@ contains
 
     end subroutine s_compute_pressure
 
+    !> Conserved energy from the primitive pressure and mixture properties: the inverse of s_compute_pressure over the same model
+    !! branches. The chemistry adapter builds the mixture internal energy from temperature and species. Arithmetic is relocated
+    !! unchanged from the conversion module.
+    subroutine s_compute_internal_energy(pres, alf, dyn_p, pi_inf, gamma, rho, qv, Ys, energy, pres_mag)
+
+        $:GPU_ROUTINE(function_name='s_compute_internal_energy',parallelism='[seq]', cray_noinline=True)
+
+        real(stp), intent(in)                          :: pres, alf
+        real(wp), intent(in)                           :: dyn_p, pi_inf, gamma, rho, qv
+        real(stp), intent(out)                         :: energy
+        real(wp), intent(in), optional                 :: pres_mag
+        real(wp), dimension(1:num_species), intent(in) :: Ys
+        real(wp)                                       :: e_mix, mix_mol_weight, T
+        #:if not chemistry
+            ! Computing the energy from the pressure
+            if (mhd) then
+                ! MHD energy includes magnetic pressure contribution
+                energy = gamma*pres + dyn_p + pres_mag + pi_inf + qv
+            else if ((model_eqns /= model_eqns_4eq) .and. (bubbles_euler .neqv. .true.)) then
+                ! Five-equation model (Allaire et al. JCP 2002): E = Gamma*p + 0.5*rho*|u|^2 + pi_inf + qv
+                energy = gamma*pres + dyn_p + pi_inf + qv
+            else if ((model_eqns /= model_eqns_4eq) .and. bubbles_euler) then
+                ! Bubble-augmented energy with void fraction correction
+                energy = dyn_p + (1._wp - alf)*(gamma*pres + pi_inf)
+            else
+                ! Four-equation model (Kapila et al. PoF 2001): Tait EOS, no conserved energy variable
+                energy = 0._wp
+            end if
+        #:else
+            ! Reacting mixture: compute conserved energy from species mass fractions and temperature
+            call get_mixture_molecular_weight(Ys, mix_mol_weight)
+            T = pres*mix_mol_weight/(gas_constant*rho)
+            call get_mixture_energy_mass(T, Ys, e_mix)
+            energy = dyn_p + rho*e_mix
+        #:endif
+
+    end subroutine s_compute_internal_energy
+
 #ifndef MFC_PRE_PROCESS
     !> Speed of sound from thermodynamic state, supporting the several equation-of-state models. Arithmetic is relocated unchanged
     !! from the conversion module: the non-chemistry branches hold c squared until a single final sqrt, the chemistry branch returns
diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp
index 4a3f2e0c6b..695f19e194 100644
--- a/src/common/m_variables_conversion.fpp
+++ b/src/common/m_variables_conversion.fpp
@@ -18,7 +18,7 @@ module m_variables_conversion
     use m_thermochem, only: num_species, get_temperature, get_pressure, gas_constant, get_mixture_molecular_weight, &
         & get_mixture_energy_mass
 
-    use m_thermodynamics, only: s_compute_pressure
+    use m_thermodynamics, only: s_compute_pressure, s_compute_internal_energy
 #ifndef MFC_PRE_PROCESS
     use m_thermodynamics, only: s_compute_speed_of_sound
 #endif
@@ -731,7 +731,6 @@ contains
         real(wp), dimension(2)           :: Re_K
         integer                          :: i, j, k, l  !< Generic loop iterators
         real(wp), dimension(num_species) :: Ys
-        real(wp)                         :: e_mix, mix_mol_weight, T
         real(wp)                         :: pres_mag
         real(wp)                         :: Ga          !< Lorentz factor (gamma in relativity)
         real(wp)                         :: h           !< relativistic enthalpy
@@ -829,43 +828,28 @@ contains
                     end do
 
                     if (chemistry) then
-                        ! Reacting mixture: compute conserved energy from species mass fractions and temperature
                         do i = eqn_idx%species%beg, eqn_idx%species%end
                             Ys(i - eqn_idx%species%beg + 1) = q_prim_vf(i)%sf(j, k, l)
                             q_cons_vf(i)%sf(j, k, l) = rho*q_prim_vf(i)%sf(j, k, l)
                         end do
+                    end if
 
-                        call get_mixture_molecular_weight(Ys, mix_mol_weight)
-                        T = q_prim_vf(eqn_idx%E)%sf(j, k, l)*mix_mol_weight/(gas_constant*rho)
-                        call get_mixture_energy_mass(T, Ys, e_mix)
-
-                        q_cons_vf(eqn_idx%E)%sf(j, k, l) = dyn_pres + rho*e_mix
-                    else
-                        ! Computing the energy from the pressure
-                        if (mhd) then
-                            if (n == 0) then
-                                pres_mag = 0.5_wp*(Bx0**2 + q_prim_vf(eqn_idx%B%beg)%sf(j, k, &
-                                                   & l)**2 + q_prim_vf(eqn_idx%B%beg + 1)%sf(j, k, l)**2)
-                            else
-                                pres_mag = 0.5_wp*(q_prim_vf(eqn_idx%B%beg)%sf(j, k, l)**2 + q_prim_vf(eqn_idx%B%beg + 1)%sf(j, &
-                                                   & k, l)**2 + q_prim_vf(eqn_idx%B%beg + 2)%sf(j, k, l)**2)
-                            end if
-                            ! MHD energy includes magnetic pressure contribution
-                            q_cons_vf(eqn_idx%E)%sf(j, k, l) = gamma*q_prim_vf(eqn_idx%E)%sf(j, k, &
-                                      & l) + dyn_pres + pres_mag + pi_inf + qv
-                        else if ((model_eqns /= model_eqns_4eq) .and. (bubbles_euler .neqv. .true.)) then
-                            ! Five-equation model (Allaire et al. JCP 2002): E = Gamma*p + 0.5*rho*|u|^2 + pi_inf + qv
-                            q_cons_vf(eqn_idx%E)%sf(j, k, l) = gamma*q_prim_vf(eqn_idx%E)%sf(j, k, l) + dyn_pres + pi_inf + qv
-                        else if ((model_eqns /= model_eqns_4eq) .and. (bubbles_euler)) then
-                            ! Bubble-augmented energy with void fraction correction
-                            q_cons_vf(eqn_idx%E)%sf(j, k, l) = dyn_pres + (1._wp - q_prim_vf(eqn_idx%alf)%sf(j, k, &
-                                      & l))*(gamma*q_prim_vf(eqn_idx%E)%sf(j, k, l) + pi_inf)
+                    if (mhd) then
+                        if (n == 0) then
+                            pres_mag = 0.5_wp*(Bx0**2 + q_prim_vf(eqn_idx%B%beg)%sf(j, k, &
+                                               & l)**2 + q_prim_vf(eqn_idx%B%beg + 1)%sf(j, k, l)**2)
                         else
-                            ! Four-equation model (Kapila et al. PoF 2001): Tait EOS, no conserved energy variable
-                            q_cons_vf(eqn_idx%E)%sf(j, k, l) = 0._wp
+                            pres_mag = 0.5_wp*(q_prim_vf(eqn_idx%B%beg)%sf(j, k, l)**2 + q_prim_vf(eqn_idx%B%beg + 1)%sf(j, k, &
+                                               & l)**2 + q_prim_vf(eqn_idx%B%beg + 2)%sf(j, k, l)**2)
                         end if
+                    else
+                        pres_mag = 0._wp
                     end if
 
+                    call s_compute_internal_energy(q_prim_vf(eqn_idx%E)%sf(j, k, l), q_prim_vf(eqn_idx%alf)%sf(j, k, l), &
+                                                   & dyn_pres, pi_inf, gamma, rho, qv, Ys, q_cons_vf(eqn_idx%E)%sf(j, k, l), &
+                                                   & pres_mag=pres_mag)
+
                     ! Six-equation model (Saurel et al. JCP 2009): compute per-phase internal energies
                     if (model_eqns == model_eqns_6eq) then
                         do i = 1, num_fluids

From e1ac11f5fb8799fe5c83e77315a70263031abc11 Mon Sep 17 00:00:00 2001
From: Fahad Nabid <55717325+fahnab666@users.noreply.github.com>
Date: Mon, 20 Jul 2026 21:09:54 -0400
Subject: [PATCH 4/6] PR1 increment 4: route rho*c^2 bulk-modulus derivative
 through m_thermodynamics

Add f_bulk_modulus, the stiffened-gas thermodynamic derivative
((gamma+1)p+pi_inf)/gamma = rho*c^2, as an inlinable pure function in the
central thermodynamics interface. Route the six identical inline sites to it:
the two alt_soundspeed Wood-mixture branches in s_compute_speed_of_sound, the
two five-equation K-div-u blkmod sites in m_rhs, and the two post-process
derived sound-speed blkmod sites. The fused /(gamma*rho) post-process branch is
left inline (different arithmetic, not this derivative). Bit-identical: 616
passed, 0 failed.
---
 src/common/m_thermodynamics.fpp          | 19 +++++++++++++++++--
 src/post_process/m_derived_variables.fpp |  5 +++--
 src/simulation/m_rhs.fpp                 |  9 +++++----
 3 files changed, 25 insertions(+), 8 deletions(-)

diff --git a/src/common/m_thermodynamics.fpp b/src/common/m_thermodynamics.fpp
index 2dc133df12..91b9a0f6fa 100644
--- a/src/common/m_thermodynamics.fpp
+++ b/src/common/m_thermodynamics.fpp
@@ -24,6 +24,7 @@ module m_thermodynamics
     public :: s_compute_internal_energy
 #ifndef MFC_PRE_PROCESS
     public :: s_compute_speed_of_sound
+    public :: f_bulk_modulus
 #endif
 
 contains
@@ -136,6 +137,20 @@ contains
     end subroutine s_compute_internal_energy
 
 #ifndef MFC_PRE_PROCESS
+    !> Stiffened-gas bulk modulus, the thermodynamic derivative $\rho c^2 = ((\gamma+1)p+\pi_\infty)/\gamma$. Shared by the
+    !! alt_soundspeed Wood mixture sound speed here, the five-equation K-div-u source in m_rhs, and the post-process derived sound
+    !! speed, so the one derivative lives in one place.
+    pure function f_bulk_modulus(gamma, pres, pi_inf) result(blkmod)
+
+        $:GPU_ROUTINE(parallelism='[seq]')
+
+        real(wp), intent(in) :: gamma, pres, pi_inf
+        real(wp)             :: blkmod
+
+        blkmod = ((gamma + 1._wp)*pres + pi_inf)/gamma
+
+    end function f_bulk_modulus
+
     !> Speed of sound from thermodynamic state, supporting the several equation-of-state models. Arithmetic is relocated unchanged
     !! from the conversion module: the non-chemistry branches hold c squared until a single final sqrt, the chemistry branch returns
     !! c already square-rooted, and the mixture_err floor is preserved.
@@ -167,8 +182,8 @@ contains
             c = sqrt((1._wp + 1._wp/gamma)*pres/rho/H)
         else
             if (alt_soundspeed) then  ! Wood's mixture sound speed via bulk moduli
-                blkmod1 = ((gammas(1) + 1._wp)*pres + pi_infs(1))/gammas(1)
-                blkmod2 = ((gammas(2) + 1._wp)*pres + pi_infs(2))/gammas(2)
+                blkmod1 = f_bulk_modulus(gammas(1), pres, pi_infs(1))
+                blkmod2 = f_bulk_modulus(gammas(2), pres, pi_infs(2))
                 c = (1._wp/(rho*(adv(1)/blkmod1 + adv(2)/blkmod2)))
             else if (model_eqns == model_eqns_6eq) then  ! Six-equation model sound speed
                 c = 0._wp
diff --git a/src/post_process/m_derived_variables.fpp b/src/post_process/m_derived_variables.fpp
index 05f1fc1805..26f6372d5f 100644
--- a/src/post_process/m_derived_variables.fpp
+++ b/src/post_process/m_derived_variables.fpp
@@ -11,6 +11,7 @@ module m_derived_variables
     use m_mpi_proxy
     use m_helper_basic
     use m_variables_conversion
+    use m_thermodynamics, only: f_bulk_modulus
     use m_constants, only: model_eqns_gamma_law
 
     implicit none
@@ -105,8 +106,8 @@ contains
                         q_sf(i, j, k) = (((gamma_sf(i, j, k) + 1._wp)*q_prim_vf(eqn_idx%E)%sf(i, j, k) + pi_inf_sf(i, j, &
                              & k))/(gamma_sf(i, j, k)*rho_sf(i, j, k)))
                     else
-                        blkmod1 = ((gammas(1) + 1._wp)*q_prim_vf(eqn_idx%E)%sf(i, j, k) + pi_infs(1))/gammas(1)
-                        blkmod2 = ((gammas(2) + 1._wp)*q_prim_vf(eqn_idx%E)%sf(i, j, k) + pi_infs(2))/gammas(2)
+                        blkmod1 = f_bulk_modulus(gammas(1), q_prim_vf(eqn_idx%E)%sf(i, j, k), pi_infs(1))
+                        blkmod2 = f_bulk_modulus(gammas(2), q_prim_vf(eqn_idx%E)%sf(i, j, k), pi_infs(2))
                         q_sf(i, j, k) = (1._wp/(rho_sf(i, j, k)*(q_prim_vf(eqn_idx%adv%beg)%sf(i, j, &
                              & k)/blkmod1 + (1._wp - q_prim_vf(eqn_idx%adv%beg)%sf(i, j, k))/blkmod2)))
                     end if
diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp
index 8176f7b838..0f316c7315 100644
--- a/src/simulation/m_rhs.fpp
+++ b/src/simulation/m_rhs.fpp
@@ -13,6 +13,7 @@ module m_rhs
     use m_global_parameters
     use m_mpi_proxy
     use m_variables_conversion
+    use m_thermodynamics, only: f_bulk_modulus
     use m_weno
     use m_constants, only: riemann_solver_hll, riemann_solver_hlld, model_eqns_6eq, int_comp_mthinc, recon_type_weno, &
         & recon_type_muscl
@@ -870,10 +871,10 @@ contains
             do q_loop = 0, p
                 do l_loop = 0, n
                     do k_loop = 0, m
-                        blkmod1(k_loop, l_loop, q_loop) = ((gammas(1) + 1._wp)*q_prim_vf%vf(eqn_idx%E)%sf(k_loop, l_loop, &
-                                & q_loop) + pi_infs(1))/gammas(1)
-                        blkmod2(k_loop, l_loop, q_loop) = ((gammas(2) + 1._wp)*q_prim_vf%vf(eqn_idx%E)%sf(k_loop, l_loop, &
-                                & q_loop) + pi_infs(2))/gammas(2)
+                        blkmod1(k_loop, l_loop, q_loop) = f_bulk_modulus(gammas(1), q_prim_vf%vf(eqn_idx%E)%sf(k_loop, l_loop, &
+                                & q_loop), pi_infs(1))
+                        blkmod2(k_loop, l_loop, q_loop) = f_bulk_modulus(gammas(2), q_prim_vf%vf(eqn_idx%E)%sf(k_loop, l_loop, &
+                                & q_loop), pi_infs(2))
                         alpha1(k_loop, l_loop, q_loop) = q_cons_vf%vf(eqn_idx%adv%beg)%sf(k_loop, l_loop, q_loop)
 
                         if (bubbles_euler) then

From d456e8fab8717203b936a0da29d3bb3bc8130af2 Mon Sep 17 00:00:00 2001
From: Fahad Nabid <55717325+fahnab666@users.noreply.github.com>
Date: Tue, 21 Jul 2026 00:49:58 -0400
Subject: [PATCH 5/6] Guard void-fraction index in prim->cons energy call for
 non-bubble configs

eqn_idx%alf is 0 outside euler-bubble models (model_eqns=3 in pre/post_process, where the assignment is MFC_SIMULATION-only). PR1 hoisted q_prim_vf(eqn_idx%alf) into an unconditional actual argument, so bounds-checked reldebug CI traps on q_prim_vf(0). The value is read only in the routine's bubble branch, so clamp the index with max(1, ...); results are bit-identical.
---
 src/common/m_variables_conversion.fpp | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp
index 695f19e194..6cefde3148 100644
--- a/src/common/m_variables_conversion.fpp
+++ b/src/common/m_variables_conversion.fpp
@@ -846,7 +846,8 @@ contains
                         pres_mag = 0._wp
                     end if
 
-                    call s_compute_internal_energy(q_prim_vf(eqn_idx%E)%sf(j, k, l), q_prim_vf(eqn_idx%alf)%sf(j, k, l), &
+                    ! eqn_idx%alf = 0 for non-bubble configs; clamp the index (value read only in the bubble branch)
+                    call s_compute_internal_energy(q_prim_vf(eqn_idx%E)%sf(j, k, l), q_prim_vf(max(1, eqn_idx%alf))%sf(j, k, l), &
                                                    & dyn_pres, pi_inf, gamma, rho, qv, Ys, q_cons_vf(eqn_idx%E)%sf(j, k, l), &
                                                    & pres_mag=pres_mag)
 

From c861a5c17e8daceb0ae907249cc2197761cefb83 Mon Sep 17 00:00:00 2001
From: Fahad Nabid <55717325+fahnab666@users.noreply.github.com>
Date: Tue, 21 Jul 2026 01:45:35 -0400
Subject: [PATCH 6/6] Complete the six-operation thermodynamics interface

Add the two remaining operations so all six live in m_thermodynamics:
f_compute_temperature (ideal-gas T from pressure, routed from the chemistry
cons->prim path) and f_validate_state (sound-speed admissibility predicate).
The two scattered mixture_err sound-speed guards, in the interface and in
post-process, now consult f_validate_state instead of testing c<0 inline, and
the post-process derived sound speed routes through f_bulk_modulus like the
alt_soundspeed branch. Bit-identical: floor values are unchanged and the
temperature routing keeps operand order.
---
 src/common/m_thermodynamics.fpp          | 42 +++++++++++++++++++++---
 src/common/m_variables_conversion.fpp    |  4 +--
 src/post_process/m_derived_variables.fpp |  8 ++---
 3 files changed, 43 insertions(+), 11 deletions(-)

diff --git a/src/common/m_thermodynamics.fpp b/src/common/m_thermodynamics.fpp
index 91b9a0f6fa..e032b750c5 100644
--- a/src/common/m_thermodynamics.fpp
+++ b/src/common/m_thermodynamics.fpp
@@ -5,10 +5,11 @@
 #:include 'macros.fpp'
 #:include 'case.fpp'
 
-!> @brief Central equation-of-state interface. Every genuine thermodynamic query (pressure here; energy, temperature, sound speed,
-!! and derivatives in later PRs) is answered through one module so no solver path carries its own EOS algebra. Compile-time
-!! chemistry selects the adapter: the stiffened-gas analytic branch or the Pyrometheus temperature-then-query branch. This is a leaf
-!! module (it uses global parameters and the thermochemistry surface, never m_variables_conversion) so no cycle forms.
+!> @brief Central equation-of-state interface. Every genuine thermodynamic query (pressure, internal energy, temperature, sound
+!! speed, and the rho c^2 derivative) is answered through one module so no solver path carries its own EOS algebra. Compile-time
+!! chemistry selects the adapter: the stiffened_gas adapter (analytic mixture EOS) or the pyrometheus adapter
+!! (temperature-then-query against the thermochemistry surface). This is a leaf module (it uses global parameters and the
+!! thermochemistry surface, never m_variables_conversion) so no cycle forms.
 module m_thermodynamics
 
     use m_derived_types
@@ -22,6 +23,8 @@ module m_thermodynamics
     private
     public :: s_compute_pressure
     public :: s_compute_internal_energy
+    public :: f_compute_temperature
+    public :: f_validate_state
 #ifndef MFC_PRE_PROCESS
     public :: s_compute_speed_of_sound
     public :: f_bulk_modulus
@@ -136,6 +139,35 @@ contains
 
     end subroutine s_compute_internal_energy
 
+    !> Mixture temperature from pressure via the ideal-gas relation T = p/(rho R_gas), where R_gas is the specific gas constant of
+    !! the local composition. The one temperature query shared by the conversion routines, so the chemistry cons->prim path no
+    !! longer carries its own copy. Arithmetic is relocated unchanged from the conversion module.
+    pure function f_compute_temperature(pres, rho, R_gas) result(T)
+
+        $:GPU_ROUTINE(parallelism='[seq]')
+
+        real(wp), intent(in) :: pres, rho, R_gas
+        real(wp)             :: T
+
+        T = pres/rho/R_gas
+
+    end function f_compute_temperature
+
+    !> Thermodynamic-state admissibility: a state is physical only when its squared sound speed is non-negative, i.e. the wave speed
+    !! is real. The mixture_err sound-speed guards in the interface and in post-processing both consult this one predicate instead
+    !! of testing the sign inline, so the validity criterion lives in one place. Pure and device-safe: it returns a flag and never
+    !! aborts, because GPU device code cannot. Later adapters (Mie-Gruneisen, JWL) tighten the admissible domain here.
+    pure function f_validate_state(sound_speed_sq) result(is_admissible)
+
+        $:GPU_ROUTINE(parallelism='[seq]')
+
+        real(wp), intent(in) :: sound_speed_sq
+        logical              :: is_admissible
+
+        is_admissible = (sound_speed_sq >= 0._wp)
+
+    end function f_validate_state
+
 #ifndef MFC_PRE_PROCESS
     !> Stiffened-gas bulk modulus, the thermodynamic derivative $\rho c^2 = ((\gamma+1)p+\pi_\infty)/\gamma$. Shared by the
     !! alt_soundspeed Wood mixture sound speed here, the five-equation K-div-u source in m_rhs, and the post-process derived sound
@@ -204,7 +236,7 @@ contains
                 c = (H - 5.e-1*vel_sum - qv/rho)/gamma
             end if
 
-            if (mixture_err .and. c < 0._wp) then
+            if (mixture_err .and. .not. f_validate_state(c)) then
                 c = 100._wp*sgm_eps
             else
                 c = sqrt(c)
diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp
index 6cefde3148..4f15e2b067 100644
--- a/src/common/m_variables_conversion.fpp
+++ b/src/common/m_variables_conversion.fpp
@@ -18,7 +18,7 @@ module m_variables_conversion
     use m_thermochem, only: num_species, get_temperature, get_pressure, gas_constant, get_mixture_molecular_weight, &
         & get_mixture_energy_mass
 
-    use m_thermodynamics, only: s_compute_pressure, s_compute_internal_energy
+    use m_thermodynamics, only: s_compute_pressure, s_compute_internal_energy, f_compute_temperature
 #ifndef MFC_PRE_PROCESS
     use m_thermodynamics, only: s_compute_speed_of_sound
 #endif
@@ -1032,7 +1032,7 @@ contains
                         ! Computing the energy from the internal energy of the mixture
                         call get_mixture_molecular_weight(Y_k, mix_mol_weight)
                         R_gas = gas_constant/mix_mol_weight
-                        T_K = pres_K/rho_K/R_gas
+                        T_K = f_compute_temperature(pres_K, rho_K, R_gas)
                         call get_mixture_energy_mass(T_K, Y_K, E_K)
                         E_K = rho_K*E_K + 5.e-1_wp*rho_K*vel_K_sum
                     else
diff --git a/src/post_process/m_derived_variables.fpp b/src/post_process/m_derived_variables.fpp
index 26f6372d5f..82e91b99c9 100644
--- a/src/post_process/m_derived_variables.fpp
+++ b/src/post_process/m_derived_variables.fpp
@@ -11,7 +11,7 @@ module m_derived_variables
     use m_mpi_proxy
     use m_helper_basic
     use m_variables_conversion
-    use m_thermodynamics, only: f_bulk_modulus
+    use m_thermodynamics, only: f_bulk_modulus, f_validate_state
     use m_constants, only: model_eqns_gamma_law
 
     implicit none
@@ -103,8 +103,8 @@ contains
             do j = -offset_y%beg, n + offset_y%end
                 do i = -offset_x%beg, m + offset_x%end
                     if (alt_soundspeed .neqv. .true.) then
-                        q_sf(i, j, k) = (((gamma_sf(i, j, k) + 1._wp)*q_prim_vf(eqn_idx%E)%sf(i, j, k) + pi_inf_sf(i, j, &
-                             & k))/(gamma_sf(i, j, k)*rho_sf(i, j, k)))
+                        q_sf(i, j, k) = f_bulk_modulus(gamma_sf(i, j, k), q_prim_vf(eqn_idx%E)%sf(i, j, k), pi_inf_sf(i, j, &
+                             & k))/rho_sf(i, j, k)
                     else
                         blkmod1 = f_bulk_modulus(gammas(1), q_prim_vf(eqn_idx%E)%sf(i, j, k), pi_infs(1))
                         blkmod2 = f_bulk_modulus(gammas(2), q_prim_vf(eqn_idx%E)%sf(i, j, k), pi_infs(2))
@@ -112,7 +112,7 @@ contains
                              & k)/blkmod1 + (1._wp - q_prim_vf(eqn_idx%adv%beg)%sf(i, j, k))/blkmod2)))
                     end if
 
-                    if (mixture_err .and. q_sf(i, j, k) < 0._wp) then
+                    if (mixture_err .and. .not. f_validate_state(q_sf(i, j, k))) then
                         q_sf(i, j, k) = 1.e-16_wp
                     else
                         q_sf(i, j, k) = sqrt(q_sf(i, j, k))