From 936c676fc6ee8ba331997969ac951d29a0196723 Mon Sep 17 00:00:00 2001
From: Spencer Bryngelson <sbryngelson@gmail.com>
Date: Sat, 8 Aug 2026 21:34:38 -0400
Subject: [PATCH 1/3] Pass an eos_state to s_compute_speed_of_sound

s_compute_speed_of_sound took ten arguments, and the contract between two of
them was invisible at the call site: H must include qv, because the routine
subtracts qv/rho internally. Nothing said so, and four call sites got it wrong
(#1707) while a fifth was right only because two separate omissions cancelled.

Introduces type(eos_state) in m_derived_types carrying the scalars the routine
needs, and two constructors:

  - s_eos_state derives H from the other members, so it cannot disagree with
    qv. The defect in #1707 is unrepresentable through this path.
  - s_eos_state_roe takes H explicitly, for the Roe-averaged Riemann paths,
    the chemistry Roe branch and the relativistic branch, all of which pass an
    H that is deliberately not the exact state enthalpy.

The call becomes (state, adv, c) instead of ten positional arguments. All 24
call sites across nine files are converted, and the state is added to the
private() clause of every GPU parallel loop that builds one.

adv stays a separate argument rather than a component. A derived-type
component cannot have a runtime extent, and num_fluids is a parameter only
under case optimization, so dimension(num_fluids) does not compile in a
general build; padding to num_fluids_max would compile but place ten reals in
a per-cell private struct on device.

The three probe sites in simulation/m_data_output that open-coded H now use
s_eos_state, which incorporates the #1707 correction as a consequence of the
interface rather than as a separate patch. That changes the reported probe
sound speed where qv /= 0; no golden observes it, because the packer keeps
only the last column of probe output (#1711, fixed separately in #1712).

Verified: builds clean; full suite 627 passed, 0 failed, no golden
regenerated. A scan confirms no GPU parallel loop uses a state without
declaring it private.
---
 src/common/m_derived_types.fpp           | 23 +++++++
 src/common/m_variables_conversion.fpp    | 77 ++++++++++++++++++------
 src/post_process/m_data_output.fpp       |  4 +-
 src/post_process/m_start_up.fpp          | 14 +++--
 src/simulation/m_cbc.fpp                 |  6 +-
 src/simulation/m_data_output.fpp         | 21 ++++---
 src/simulation/m_riemann_solver_hll.fpp  | 18 +++---
 src/simulation/m_riemann_solver_hllc.fpp | 52 ++++++++--------
 src/simulation/m_riemann_solver_hlld.fpp | 19 +++---
 src/simulation/m_riemann_solver_lf.fpp   | 11 ++--
 src/simulation/m_time_steppers.fpp       |  8 ++-
 11 files changed, 166 insertions(+), 87 deletions(-)

diff --git a/src/common/m_derived_types.fpp b/src/common/m_derived_types.fpp
index 872c04b39..3a45cd2bc 100644
--- a/src/common/m_derived_types.fpp
+++ b/src/common/m_derived_types.fpp
@@ -121,6 +121,29 @@ module m_derived_types
         type(int_bounds_info) :: x, y, z
     end type bc_xyz_info
 
+    !> Thermodynamic state handed to the equation-of-state routines.
+    !>
+    !> Carries scalars only. The volume-fraction array `adv` stays a separate argument: a derived-type
+    !> component cannot have a runtime extent, and `num_fluids` is a parameter only under case
+    !> optimization, so `dimension(num_fluids)` will not compile in a general build. Padding to
+    !> num_fluids_max would compile but put ten reals into a per-cell private struct on device.
+    !>
+    !> `H` is the specific total enthalpy and must include `qv`, because the sound-speed relation
+    !> subtracts `qv/rho`. Build states with f_eos_state so that invariant holds by construction
+    !> rather than by convention; f_eos_state_roe exists for the Roe-averaged paths, which supply an
+    !> `H` that is deliberately not the exact state enthalpy.
+    !> Contains no allocatable members - safe to use inside device routines.
+    type eos_state
+        real(wp) :: rho      !< Mixture density
+        real(wp) :: pres     !< Pressure
+        real(wp) :: gamma    !< Stiffened-gas gamma (1/(Gamma-1))
+        real(wp) :: pi_inf   !< Stiffened-gas stiffness
+        real(wp) :: qv       !< Heat of formation (volumetric)
+        real(wp) :: vel_sum  !< |u|^2
+        real(wp) :: H        !< Specific total enthalpy, including qv
+        real(wp) :: c_c      !< Roe-averaged chemistry sound-speed term (0 when unused)
+    end type eos_state
+
     !> QBMM moment index mappings - separate from bub beg/end so eqn_idx contains no allocatables.
     type qbmm_idx_info
         integer, dimension(:), allocatable     :: rs       !< R moment indices per bubble bin
diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp
index 78d59323d..ee2c611fe 100644
--- a/src/common/m_variables_conversion.fpp
+++ b/src/common/m_variables_conversion.fpp
@@ -25,7 +25,7 @@ module m_variables_conversion
         & s_convert_mixture_to_mixture_variables, s_convert_species_to_mixture_variables, &
         & s_convert_species_to_mixture_variables_kernel, s_convert_conservative_to_primitive_variables, &
         & s_convert_primitive_to_conservative_variables, s_convert_primitive_to_flux_variables, s_compute_pressure, &
-        & s_compute_species_fraction, s_compute_speed_of_sound, s_compute_fast_magnetosonic_speed, &
+        & s_compute_species_fraction, s_compute_speed_of_sound, s_compute_fast_magnetosonic_speed, s_eos_state, s_eos_state_roe, &
         & s_finalize_variables_conversion_module, gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps
 
     real(wp), allocatable, dimension(:)   :: Gs_vc
@@ -1164,54 +1164,93 @@ contains
     end subroutine s_finalize_variables_conversion_module
 
     !> 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)
+    !> Build an exact thermodynamic state. The specific total enthalpy is derived here rather than supplied, so it cannot disagree
+    !! with qv: H = (E + p)/rho with E = gamma*p + pi_inf + qv + rho*|u|^2/2. Every caller that previously open-coded the closed
+    !! form should use this, which makes the defect in #1707 unrepresentable.
+    subroutine s_eos_state(s, pres, rho, gamma, pi_inf, qv, vel_sum)
+
+        $:GPU_ROUTINE(function_name='s_eos_state', parallelism='[seq]', cray_inline=True)
+
+        type(eos_state), intent(out) :: s
+        real(wp), intent(in)         :: pres, rho, gamma, pi_inf, qv, vel_sum
+
+        s%pres = pres
+        s%rho = rho
+        s%gamma = gamma
+        s%pi_inf = pi_inf
+        s%qv = qv
+        s%vel_sum = vel_sum
+        s%c_c = 0._wp
+        s%H = ((gamma + 1._wp)*pres + pi_inf + qv)/rho + 5.e-1_wp*vel_sum
+
+    end subroutine s_eos_state
+
+    !> Build a state whose enthalpy is supplied by the caller. The Roe-averaged Riemann paths, the chemistry Roe branch and the
+    !! relativistic branch all pass an H that is deliberately not the exact enthalpy of the state, so they cannot use s_eos_state.
+    subroutine s_eos_state_roe(s, pres, rho, gamma, pi_inf, qv, vel_sum, H, c_c)
+
+        $:GPU_ROUTINE(function_name='s_eos_state_roe', parallelism='[seq]', cray_inline=True)
+
+        type(eos_state), intent(out)   :: s
+        real(wp), intent(in)           :: pres, rho, gamma, pi_inf, qv, vel_sum, H
+        real(wp), intent(in), optional :: c_c
+
+        s%pres = pres
+        s%rho = rho
+        s%gamma = gamma
+        s%pi_inf = pi_inf
+        s%qv = qv
+        s%vel_sum = vel_sum
+        s%H = H
+        s%c_c = 0._wp
+        if (present(c_c)) s%c_c = c_c
+
+    end subroutine s_eos_state_roe
+
+    subroutine s_compute_speed_of_sound(s, adv, c)
 
         $:GPU_ROUTINE(parallelism='[seq]')
 
-        real(wp), intent(in) :: pres
-        real(wp), intent(in) :: rho, gamma, pi_inf, qv
-        real(wp), intent(in) :: H
+        type(eos_state), intent(in) :: s
         #: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))
+            if (avg_state == avg_state_roe .and. abs(s%c_c) > verysmall) then
+                c = sqrt(s%c_c - (s%gamma - 1.0_wp)*(s%vel_sum - s%H))
             else
-                c = sqrt((1.0_wp + 1.0_wp/gamma)*pres/rho)
+                c = sqrt((1.0_wp + 1.0_wp/s%gamma)*s%pres/s%rho)
             end if
         else if (relativity) then  ! Relativistic sound speed
-            c = sqrt((1._wp + 1._wp/gamma)*pres/rho/H)
+            c = sqrt((1._wp + 1._wp/s%gamma)*s%pres/s%rho/s%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)))
+                blkmod1 = ((gammas(1) + 1._wp)*s%pres + pi_infs(1))/gammas(1)
+                blkmod2 = ((gammas(2) + 1._wp)*s%pres + pi_infs(2))/gammas(2)
+                c = (1._wp/(s%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))
+                    c = c + adv(q)*gs_min(q)*(s%pres + pi_infs(q)/(gammas(q) + 1._wp))
                 end do
-                c = c/rho
+                c = c/s%rho
             else if (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
+                    c = (1._wp/s%gamma + 1._wp)*(s%pres + s%pi_inf/(s%gamma + 1._wp))/s%rho
                 else
-                    c = (1._wp/gamma + 1._wp)*(pres + pi_inf/(gamma + 1._wp))/(rho*(1._wp - adv(num_fluids)))
+                    c = (1._wp/s%gamma + 1._wp)*(s%pres + s%pi_inf/(s%gamma + 1._wp))/(s%rho*(1._wp - adv(num_fluids)))
                 end if
             else
-                c = (H - 5.e-1*vel_sum - qv/rho)/gamma
+                c = (s%H - 5.e-1*s%vel_sum - s%qv/s%rho)/s%gamma
             end if
 
             if (mixture_err .and. c < 0._wp) then
diff --git a/src/post_process/m_data_output.fpp b/src/post_process/m_data_output.fpp
index 1d5ddfedb..0a36631ad 100644
--- a/src/post_process/m_data_output.fpp
+++ b/src/post_process/m_data_output.fpp
@@ -1243,6 +1243,7 @@ contains
         real(wp), dimension(num_vels) :: vel
         real(wp), dimension(num_fluids) :: adv
         integer :: i, j, k, l, s  !< looping indices
+        type(eos_state) :: eos_s
 
         Egk = 0._wp
         Elp = 0._wp
@@ -1288,7 +1289,8 @@ contains
 
                     H = ((gamma + 1._wp)*pres + pi_inf + qv)/rho
 
-                    call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv, 0._wp, 0._wp, c, qv)
+                    call s_eos_state(eos_s, pres, rho, gamma, pi_inf, qv, 0._wp)
+                    call s_compute_speed_of_sound(eos_s, adv, c)
 
                     Ma = maxvel/c
                     if (Ma > MaxMa .and. (adv(1) > (1.0_wp - 1.0e-10_wp))) then
diff --git a/src/post_process/m_start_up.fpp b/src/post_process/m_start_up.fpp
index 9c697c0ea..0c7fa0859 100644
--- a/src/post_process/m_start_up.fpp
+++ b/src/post_process/m_start_up.fpp
@@ -187,10 +187,11 @@ contains
              & -offset_z%beg:p + offset_z%end) :: liutex_mag
         real(wp), dimension(-offset_x%beg:m + offset_x%end,-offset_y%beg:n + offset_y%end,-offset_z%beg:p + offset_z%end, &
              & 3) :: liutex_axis
-        integer       :: i, j, k, l, kx, ky, kz, kf, j_glb, k_glb, l_glb
-        character(50) :: filename
-        logical       :: file_exists
-        integer       :: x_beg, x_end, y_beg, y_end, z_beg, z_end
+        integer         :: i, j, k, l, kx, ky, kz, kf, j_glb, k_glb, l_glb
+        character(50)   :: filename
+        logical         :: file_exists
+        integer         :: x_beg, x_end, y_beg, y_end, z_beg, z_end
+        type(eos_state) :: eos_s
 
         if (output_partial_domain) then
             call s_define_output_region
@@ -533,8 +534,9 @@ contains
 
                         H = ((gamma_sf(i, j, k) + 1._wp)*pres + pi_inf_sf(i, j, k) + qv_sf(i, j, k))/rho_sf(i, j, k)
 
-                        call s_compute_speed_of_sound(pres, rho_sf(i, j, k), gamma_sf(i, j, k), pi_inf_sf(i, j, k), H, adv, &
-                                                      & 0._wp, 0._wp, c, qv_sf(i, j, k))
+                        call s_eos_state_roe(eos_s, pres, rho_sf(i, j, k), gamma_sf(i, j, k), pi_inf_sf(i, j, k), qv_sf(i, j, k), &
+                                             & 0._wp, H)
+                        call s_compute_speed_of_sound(eos_s, adv, c)
 
                         out%q_sf(i, j, k) = c
                     end do
diff --git a/src/simulation/m_cbc.fpp b/src/simulation/m_cbc.fpp
index 4964113c4..357660d19 100644
--- a/src/simulation/m_cbc.fpp
+++ b/src/simulation/m_cbc.fpp
@@ -511,6 +511,7 @@ contains
         real(wp)               :: Cv, Cp, e_mix, Mw, R_gas
         real(wp)               :: vel_K_sum, vel_dv_dt_sum
         integer                :: i, j, k, r  !< Generic loop iterators
+        type(eos_state)        :: eos_s
         ! Reshaping of inputted data and association of the FD and PI coefficients, or CBC coefficients, respectively, hinging on
         ! selected CBC coordinate direction
 
@@ -597,7 +598,7 @@ contains
                                     & dalpha_rho_ds, dpres_ds, dvel_dt, dadv_dt, dalpha_rho_dt, L, lambda, Ys, dYs_dt, dYs_ds, &
                                     & h_k, Cp_i, Gamma_i, Xs, drho_dt, dpres_dt, dpi_inf_dt, dqv_dt, dgamma_dt, rho, pres, E, H, &
                                     & gamma, pi_inf, qv, c, Ma, T, sum_Enthalpies, Cv, Cp, e_mix, Mw, R_gas, vel_K_sum, &
-                                    & vel_dv_dt_sum, i, j]', copyin='[dir_idx]')
+                                    & vel_dv_dt_sum, i, j, eos_s]', copyin='[dir_idx]')
                 do r = is3%beg, is3%end
                     do k = is2%beg, is2%end
                         ! Transferring the Primitive Variables
@@ -661,7 +662,8 @@ contains
                         H = (E + pres)/rho
 
                         ! Compute mixture sound speed
-                        call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv_local, vel_K_sum, 0._wp, c, qv)
+                        call s_eos_state_roe(eos_s, pres, rho, gamma, pi_inf, qv, vel_K_sum, H)
+                        call s_compute_speed_of_sound(eos_s, adv_local, c)
 
                         ! First-Order Spatial Derivatives of Primitive Variables
 
diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp
index e4c32e2c9..1e3c6a9fe 100644
--- a/src/simulation/m_data_output.fpp
+++ b/src/simulation/m_data_output.fpp
@@ -189,6 +189,7 @@ contains
         real(wp)               :: Rc_min_loc, Rc_min_glb      !< Rc stability extrema on local and global grids
         real(wp)               :: icfl, vcfl, ccfl, Rc
         integer                :: fl                          !< Fluid loop iterator
+        type(eos_state)        :: eos_s
 
         icfl_max_loc = 0._wp
         vcfl_max_loc = 0._wp
@@ -196,14 +197,15 @@ contains
         Rc_min_loc = huge(1.0_wp)
         ! Computing Stability Criteria at Current Time-step
         $:GPU_PARALLEL_LOOP(collapse=3, private='[j, k, l, vel, alpha, Re, rho, vel_sum, pres, gamma, pi_inf, c, H, qv, icfl, &
-                            & vcfl, Rc, ccfl, fl]', reduction='[[icfl_max_loc, vcfl_max_loc, ccfl_max_loc], [Rc_min_loc]]', &
-                            & reductionOp='[max, min]')
+                            & vcfl, Rc, ccfl, fl, eos_s]', reduction='[[icfl_max_loc, vcfl_max_loc, &
+                            & ccfl_max_loc], [Rc_min_loc]]', reductionOp='[max, min]')
         do l = 0, p
             do k = 0, n
                 do j = 0, m
                     call s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, qv, j, k, l)
 
-                    call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, alpha, vel_sum, 0._wp, c, qv)
+                    call s_eos_state_roe(eos_s, pres, rho, gamma, pi_inf, qv, vel_sum, H)
+                    call s_compute_speed_of_sound(eos_s, alpha, c)
 
                     if (any_non_newtonian) then
                         Re(1) = 0._wp
@@ -1158,6 +1160,7 @@ contains
         real(wp)                        :: rad, thickness    !< For integral quantities
         logical                         :: trigger           !< For integral quantities
         real(wp)                        :: rhoYks(1:num_species)
+        type(eos_state)                 :: eos_s
 
         T = dflt_T_guess
 
@@ -1292,8 +1295,8 @@ contains
                     end if
 
                     ! Compute mixture sound Speed
-                    call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, ((gamma + 1._wp)*pres + pi_inf)/rho, alpha, 0._wp, &
-                                                  & 0._wp, c, qv)
+                    call s_eos_state(eos_s, pres, rho, gamma, pi_inf, qv, 0._wp)
+                    call s_compute_speed_of_sound(eos_s, alpha, c)
 
                     accel = accel_mag(j - 2, k, l)
                 end if
@@ -1372,8 +1375,8 @@ contains
                             Rdot(:) = nRdot(:)/nbub
                         end if
                         ! Compute mixture sound speed
-                        call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, ((gamma + 1._wp)*pres + pi_inf)/rho, alpha, &
-                                                      & 0._wp, 0._wp, c, qv)
+                        call s_eos_state(eos_s, pres, rho, gamma, pi_inf, qv, 0._wp)
+                        call s_compute_speed_of_sound(eos_s, alpha, c)
                     end if
                 end if
             else
@@ -1431,8 +1434,8 @@ contains
                             end if
 
                             ! Compute mixture sound speed
-                            call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, ((gamma + 1._wp)*pres + pi_inf)/rho, alpha, &
-                                                          & 0._wp, 0._wp, c, qv)
+                            call s_eos_state(eos_s, pres, rho, gamma, pi_inf, qv, 0._wp)
+                            call s_compute_speed_of_sound(eos_s, alpha, c)
 
                             accel = accel_mag(j - 2, k - 2, l - 2)
                         end if
diff --git a/src/simulation/m_riemann_solver_hll.fpp b/src/simulation/m_riemann_solver_hll.fpp
index e222642c4..c89ce76fc 100644
--- a/src/simulation/m_riemann_solver_hll.fpp
+++ b/src/simulation/m_riemann_solver_hll.fpp
@@ -97,6 +97,7 @@ contains
         type(riemann_states_vec3) :: cm        !< Conservative momentum variables
         integer :: i, j, k, l                  !< Generic loop iterators
         integer :: Re_size_loc1, Re_size_loc2  !< host copies of Re_size; amdflang reads the declare-target original stale cross-TU
+        type(eos_state) :: eos_s_L, eos_s_R, eos_s_avg
         ! Populating the buffers of the left and right Riemann problem states variables, based on the choice of boundary conditions
 
         call s_populate_riemann_states_variables_buffers(qL_prim_rsx_vf, dqL_prim_dx_vf, dqL_prim_dy_vf, dqL_prim_dz_vf, &
@@ -120,8 +121,8 @@ contains
                                     & Y_L, Y_R, MW_L, MW_R, R_gas_L, R_gas_R, Cp_L, Cp_R, Cv_L, Cv_R, Gamm_L, Gamm_R, gamma_L, &
                                     & gamma_R, pi_inf_L, pi_inf_R, qv_L, qv_R, qv_avg, c_L, c_R, G_L, G_R, damage_L, damage_R, &
                                     & rho_avg, H_avg, c_avg, gamma_avg, ptilde_L, ptilde_R, vel_L_rms, vel_R_rms, vel_avg_rms, &
-                                    & Ms_L, Ms_R, pres_SL, pres_SR, alpha_L_sum, alpha_R_sum, flux_tau_L, flux_tau_R]', &
-                                    & copyin='[norm_dir]', firstprivate='[Re_size_loc1, Re_size_loc2]')
+                                    & Ms_L, Ms_R, pres_SL, pres_SR, alpha_L_sum, alpha_R_sum, flux_tau_L, flux_tau_R, eos_s_L, &
+                                    & eos_s_R, eos_s_avg]', copyin='[norm_dir]', firstprivate='[Re_size_loc1, Re_size_loc2]')
                 do l = ${Z_BND}$%beg, ${Z_BND}$%end
                     do k = ${Y_BND}$%beg, ${Y_BND}$%end
                         do j = ${X_BND}$%beg, ${X_BND}$%end
@@ -318,17 +319,18 @@ contains
 
                             @:compute_average_state()
 
-                            call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, H_L, alpha_L, vel_L_rms, 0._wp, c_L, &
-                                                          & qv_L)
+                            call s_eos_state_roe(eos_s_L, pres_L, rho_L, gamma_L, pi_inf_L, qv_L, vel_L_rms, H_L)
+                            call s_compute_speed_of_sound(eos_s_L, alpha_L, c_L)
 
-                            call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, H_R, alpha_R, vel_R_rms, 0._wp, c_R, &
-                                                          & qv_R)
+                            call s_eos_state_roe(eos_s_R, pres_R, rho_R, gamma_R, pi_inf_R, qv_R, vel_R_rms, H_R)
+                            call s_compute_speed_of_sound(eos_s_R, alpha_R, c_R)
 
                             !> The computation of c_avg does not require all the variables, and therefore the non '_avg'
                             ! variables are placeholders to call the subroutine.
 
-                            call s_compute_speed_of_sound(pres_R, rho_avg, gamma_avg, pi_inf_R, H_avg, alpha_R, vel_avg_rms, &
-                                                          & c_sum_Yi_Phi, c_avg, qv_avg)
+                            call s_eos_state_roe(eos_s_avg, pres_R, rho_avg, gamma_avg, pi_inf_R, qv_avg, vel_avg_rms, H_avg, &
+                                                 & c_sum_Yi_Phi)
+                            call s_compute_speed_of_sound(eos_s_avg, alpha_R, c_avg)
 
                             if (mhd) then
                                 call s_compute_fast_magnetosonic_speed(rho_L, c_L, B%L, norm_dir, c_fast%L, H_L)
diff --git a/src/simulation/m_riemann_solver_hllc.fpp b/src/simulation/m_riemann_solver_hllc.fpp
index 93a77ac07..561f83d48 100644
--- a/src/simulation/m_riemann_solver_hllc.fpp
+++ b/src/simulation/m_riemann_solver_hllc.fpp
@@ -111,9 +111,10 @@ contains
         real(wp) :: vel_L_tmp, vel_R_tmp
         real(wp) :: rho_Star, E_Star, p_Star, p_K_Star, vel_K_star
         real(wp) :: pres_SL, pres_SR, Ms_L, Ms_R
-        real(wp) :: zcoef, pcorr                !< low Mach number correction
-        integer  :: i, j, k, l, q               !< Generic loop iterators
-        integer  :: Re_size_loc1, Re_size_loc2  !< host copies of Re_size; amdflang reads the declare-target original stale cross-TU
+        real(wp) :: zcoef, pcorr               !< low Mach number correction
+        integer :: i, j, k, l, q               !< Generic loop iterators
+        integer :: Re_size_loc1, Re_size_loc2  !< host copies of Re_size; amdflang reads the declare-target original stale cross-TU
+        type(eos_state) :: eos_s_L, eos_s_R, eos_s_avg
         ! Populating the buffers of the left and right Riemann problem states variables, based on the choice of boundary conditions
 
         call s_populate_riemann_states_variables_buffers(qL_prim_rsx_vf, dqL_prim_dx_vf, dqL_prim_dy_vf, dqL_prim_dz_vf, &
@@ -144,7 +145,7 @@ contains
                                         & ptilde_L, ptilde_R, vel_L_rms, vel_R_rms, vel_avg_rms, vel_L_tmp, vel_R_tmp, Ms_L, &
                                         & Ms_R, pres_SL, pres_SR, alpha_L_sum, alpha_R_sum, rho_Star, E_Star, p_Star, p_K_Star, &
                                         & vel_K_star, s_L, s_R, s_M, s_P, s_S, xi_M, xi_P, xi_L, xi_R, xi_L_m1, xi_R_m1, xi_MP, &
-                                        & xi_PP]', firstprivate='[Re_size_loc1, Re_size_loc2]')
+                                        & xi_PP, eos_s_L, eos_s_R, eos_s_avg]', firstprivate='[Re_size_loc1, Re_size_loc2]')
                     do l = ${Z_BND}$%beg, ${Z_BND}$%end
                         do k = ${Y_BND}$%beg, ${Y_BND}$%end
                             do j = ${X_BND}$%beg, ${X_BND}$%end
@@ -231,16 +232,16 @@ contains
 
                                 @:compute_average_state()
 
-                                call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, H_L, alpha_L, vel_L_rms, 0._wp, &
-                                                              & c_L, qv_L)
+                                call s_eos_state_roe(eos_s_L, pres_L, rho_L, gamma_L, pi_inf_L, qv_L, vel_L_rms, H_L)
+                                call s_compute_speed_of_sound(eos_s_L, alpha_L, c_L)
 
-                                call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, H_R, alpha_R, vel_R_rms, 0._wp, &
-                                                              & c_R, qv_R)
+                                call s_eos_state_roe(eos_s_R, pres_R, rho_R, gamma_R, pi_inf_R, qv_R, vel_R_rms, H_R)
+                                call s_compute_speed_of_sound(eos_s_R, alpha_R, c_R)
 
                                 !> The computation of c_avg does not require all the variables, and therefore the non '_avg'
                                 ! variables are placeholders to call the subroutine.
-                                call s_compute_speed_of_sound(pres_R, rho_avg, gamma_avg, pi_inf_R, H_avg, alpha_R, vel_avg_rms, &
-                                                              & 0._wp, c_avg, qv_avg)
+                                call s_eos_state_roe(eos_s_avg, pres_R, rho_avg, gamma_avg, pi_inf_R, qv_avg, vel_avg_rms, H_avg)
+                                call s_compute_speed_of_sound(eos_s_avg, alpha_R, c_avg)
 
                                 if (viscous) then
                                     $:GPU_LOOP(parallelism='[seq]')
@@ -431,7 +432,7 @@ contains
                                         & alpha_R_sum, s_L, s_R, s_M, s_P, s_S, xi_M, xi_P, xi_L, xi_R, xi_L_m1, xi_R_m1, xi_MP, &
                                         & xi_PP, nbub_L, nbub_R, PbwR3Lbar, PbwR3Rbar, R3Lbar, R3Rbar, R3V2Lbar, R3V2Rbar, Ys_L, &
                                         & Ys_R, Cp_iL, Cp_iR, Xs_L, Xs_R, Gamma_iL, Gamma_iR, Yi_avg, Phi_avg, h_iL, h_iR, &
-                                        & h_avg_2]', firstprivate='[Re_size_loc1, Re_size_loc2]')
+                                        & h_avg_2, eos_s_L, eos_s_R, eos_s_avg]', firstprivate='[Re_size_loc1, Re_size_loc2]')
                     do l = ${Z_BND}$%beg, ${Z_BND}$%end
                         do k = ${Y_BND}$%beg, ${Y_BND}$%end
                             do j = ${X_BND}$%beg, ${X_BND}$%end
@@ -603,16 +604,16 @@ contains
                                     end do
                                 end if
 
-                                call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, H_L, alpha_L, vel_L_rms, 0._wp, &
-                                                              & c_L, qv_L)
+                                call s_eos_state_roe(eos_s_L, pres_L, rho_L, gamma_L, pi_inf_L, qv_L, vel_L_rms, H_L)
+                                call s_compute_speed_of_sound(eos_s_L, alpha_L, c_L)
 
-                                call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, H_R, alpha_R, vel_R_rms, 0._wp, &
-                                                              & c_R, qv_R)
+                                call s_eos_state_roe(eos_s_R, pres_R, rho_R, gamma_R, pi_inf_R, qv_R, vel_R_rms, H_R)
+                                call s_compute_speed_of_sound(eos_s_R, alpha_R, c_R)
 
                                 !> The computation of c_avg does not require all the variables, and therefore the non '_avg'
                                 ! variables are placeholders to call the subroutine.
-                                call s_compute_speed_of_sound(pres_R, rho_avg, gamma_avg, pi_inf_R, H_avg, alpha_R, vel_avg_rms, &
-                                                              & 0._wp, c_avg, qv_avg)
+                                call s_eos_state_roe(eos_s_avg, pres_R, rho_avg, gamma_avg, pi_inf_R, qv_avg, vel_avg_rms, H_avg)
+                                call s_compute_speed_of_sound(eos_s_avg, alpha_R, c_avg)
 
                                 if (viscous) then
                                     $:GPU_LOOP(parallelism='[seq]')
@@ -808,8 +809,8 @@ contains
                                         & xi_R, xi_L_m1, xi_R_m1, Ms_L, Ms_R, pres_SL, pres_SR, vel_L, vel_R, Re_L, Re_R, &
                                         & alpha_L, alpha_R, alpha_rho_L, alpha_rho_R, alpha_lim_L, alpha_lim_R, s_L, s_R, s_S, &
                                         & vel_avg_rms, pcorr, zcoef, vel_L_tmp, vel_R_tmp, Ys_L, Ys_R, Xs_L, Xs_R, Gamma_iL, &
-                                        & Gamma_iR, Cp_iL, Cp_iR, Yi_avg, Phi_avg, h_iL, h_iR, h_avg_2]', copyin='[is1, is2, &
-                                        & is3]', firstprivate='[Re_size_loc1, Re_size_loc2]')
+                                        & Gamma_iR, Cp_iL, Cp_iR, Yi_avg, Phi_avg, h_iL, h_iR, h_avg_2, eos_s_L, eos_s_R, &
+                                        & eos_s_avg]', copyin='[is1, is2, is3]', firstprivate='[Re_size_loc1, Re_size_loc2]')
                     do l = ${Z_BND}$%beg, ${Z_BND}$%end
                         do k = ${Y_BND}$%beg, ${Y_BND}$%end
                             do j = ${X_BND}$%beg, ${X_BND}$%end
@@ -941,16 +942,17 @@ contains
 
                                 @:compute_average_state()
 
-                                call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, H_L, alpha_L, vel_L_rms, 0._wp, &
-                                                              & c_L, qv_L)
+                                call s_eos_state_roe(eos_s_L, pres_L, rho_L, gamma_L, pi_inf_L, qv_L, vel_L_rms, H_L)
+                                call s_compute_speed_of_sound(eos_s_L, alpha_L, c_L)
 
-                                call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, H_R, alpha_R, vel_R_rms, 0._wp, &
-                                                              & c_R, qv_R)
+                                call s_eos_state_roe(eos_s_R, pres_R, rho_R, gamma_R, pi_inf_R, qv_R, vel_R_rms, H_R)
+                                call s_compute_speed_of_sound(eos_s_R, alpha_R, c_R)
 
                                 !> The computation of c_avg does not require all the variables, and therefore the non '_avg'
                                 !  variables are placeholders to call the subroutine.
-                                call s_compute_speed_of_sound(pres_R, rho_avg, gamma_avg, pi_inf_R, H_avg, alpha_R, vel_avg_rms, &
-                                                              & c_sum_Yi_Phi, c_avg, qv_avg)
+                                call s_eos_state_roe(eos_s_avg, pres_R, rho_avg, gamma_avg, pi_inf_R, qv_avg, vel_avg_rms, H_avg, &
+                                                     & c_sum_Yi_Phi)
+                                call s_compute_speed_of_sound(eos_s_avg, alpha_R, c_avg)
 
                                 if (viscous) then
                                     if (chemistry) then
diff --git a/src/simulation/m_riemann_solver_hlld.fpp b/src/simulation/m_riemann_solver_hlld.fpp
index c49992bed..c60398c22 100644
--- a/src/simulation/m_riemann_solver_hlld.fpp
+++ b/src/simulation/m_riemann_solver_hlld.fpp
@@ -56,10 +56,11 @@ contains
         ! normal velocity, and x is the normal direction Note: Bx is omitted as the magnetic flux is always zero in the normal
         ! direction
 
-        real(wp) :: sqrt_rhoL_star, sqrt_rhoR_star, denom_ds, sign_Bx
-        real(wp) :: vL_star, vR_star, wL_star, wR_star
-        real(wp) :: v_double, w_double, By_double, Bz_double, E_doubleL, E_doubleR, E_double
-        integer  :: i, j, k, l
+        real(wp)        :: sqrt_rhoL_star, sqrt_rhoR_star, denom_ds, sign_Bx
+        real(wp)        :: vL_star, vR_star, wL_star, wR_star
+        real(wp)        :: v_double, w_double, By_double, Bz_double, E_doubleL, E_doubleR, E_double
+        integer         :: i, j, k, l
+        type(eos_state) :: eos_s_L, eos_s_R
 
         call s_populate_riemann_states_variables_buffers(qL_prim_rsx_vf, dqL_prim_dx_vf, dqL_prim_dy_vf, dqL_prim_dz_vf, &
             & qR_prim_rsx_vf, dqR_prim_dx_vf, dqR_prim_dy_vf, dqR_prim_dz_vf, norm_dir, ix, iy, iz)
@@ -78,7 +79,7 @@ contains
                                     & U_doubleL, U_doubleR, F_L, F_R, F_starL, F_starR, F_hlld, s_L, s_R, s_M, s_starL, s_starR, &
                                     & pTot_L, pTot_R, p_star, rhoL_star, rhoR_star, E_starL, E_starR, sqrt_rhoL_star, &
                                     & sqrt_rhoR_star, denom_ds, sign_Bx, vL_star, vR_star, wL_star, wR_star, v_double, w_double, &
-                                    & By_double, Bz_double, E_doubleL, E_doubleR, E_double]', copyin='[norm_dir]')
+                                    & By_double, Bz_double, E_doubleL, E_doubleR, E_double, eos_s_L, eos_s_R]', copyin='[norm_dir]')
                 do l = ${Z_BND}$%beg, ${Z_BND}$%end
                     do k = ${Y_BND}$%beg, ${Y_BND}$%end
                         do j = ${X_BND}$%beg, ${X_BND}$%end
@@ -147,10 +148,10 @@ contains
                             H_no_mag%R = (E%R + pres%R - pres_mag%R)/rho%R
 
                             ! (2) Compute fast wave speeds
-                            call s_compute_speed_of_sound(pres%L, rho%L, gamma%L, pi_inf%L, H_no_mag%L, alpha_L, vel_rms%L, &
-                                                          & 0._wp, c%L, qv%L)
-                            call s_compute_speed_of_sound(pres%R, rho%R, gamma%R, pi_inf%R, H_no_mag%R, alpha_R, vel_rms%R, &
-                                                          & 0._wp, c%R, qv%R)
+                            call s_eos_state_roe(eos_s_L, pres%L, rho%L, gamma%L, pi_inf%L, qv%L, vel_rms%L, H_no_mag%L)
+                            call s_compute_speed_of_sound(eos_s_L, alpha_L, c%L)
+                            call s_eos_state_roe(eos_s_R, pres%R, rho%R, gamma%R, pi_inf%R, qv%R, vel_rms%R, H_no_mag%R)
+                            call s_compute_speed_of_sound(eos_s_R, alpha_R, c%R)
                             call s_compute_fast_magnetosonic_speed(rho%L, c%L, B%L, norm_dir, c_fast%L, H_no_mag%L)
                             call s_compute_fast_magnetosonic_speed(rho%R, c%R, B%R, norm_dir, c_fast%R, H_no_mag%R)
 
diff --git a/src/simulation/m_riemann_solver_lf.fpp b/src/simulation/m_riemann_solver_lf.fpp
index 512a576e7..ed70e05ae 100644
--- a/src/simulation/m_riemann_solver_lf.fpp
+++ b/src/simulation/m_riemann_solver_lf.fpp
@@ -93,6 +93,7 @@ contains
         integer :: i, j, k, l  !< Generic loop iterators
         integer :: Re_size_loc1, Re_size_loc2  !< host copies of Re_size; amdflang reads the declare-target original stale cross-TU
         integer, dimension(3) :: idx_right_phys  !< Physical (j,k,l) indices for right state.
+        type(eos_state) :: eos_s_L, eos_s_R
         ! Populating the buffers of the left and right Riemann problem states variables, based on the choice of boundary conditions
 
         call s_populate_riemann_states_variables_buffers(qL_prim_rsx_vf, dqL_prim_dx_vf, dqL_prim_dy_vf, dqL_prim_dz_vf, &
@@ -116,7 +117,7 @@ contains
                                     & c_avg, pres_L, pres_R, rho_L, rho_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, qv_L, qv_R, c_L, &
                                     & c_R, E_L, E_R, H_L, H_R, ptilde_L, ptilde_R, s_M, s_P, xi_M, xi_P, Cp_avg, Cv_avg, T_avg, &
                                     & eps, c_sum_Yi_Phi, Cp_L, Cp_R, Cv_L, Cv_R, R_gas_L, R_gas_R, MW_L, MW_R, T_L, T_R, Y_L, &
-                                    & Y_R]', firstprivate='[Re_size_loc1, Re_size_loc2]')
+                                    & Y_R, eos_s_L, eos_s_R]', firstprivate='[Re_size_loc1, Re_size_loc2]')
                 do l = ${Z_BND}$%beg, ${Z_BND}$%end
                     do k = ${Y_BND}$%beg, ${Y_BND}$%end
                         do j = ${X_BND}$%beg, ${X_BND}$%end
@@ -291,11 +292,11 @@ contains
                                 H_R = (E_R + pres_R)/rho_R
                             end if
 
-                            call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, H_L, alpha_L, vel_L_rms, 0._wp, c_L, &
-                                                          & qv_L)
+                            call s_eos_state_roe(eos_s_L, pres_L, rho_L, gamma_L, pi_inf_L, qv_L, vel_L_rms, H_L)
+                            call s_compute_speed_of_sound(eos_s_L, alpha_L, c_L)
 
-                            call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, H_R, alpha_R, vel_R_rms, 0._wp, c_R, &
-                                                          & qv_R)
+                            call s_eos_state_roe(eos_s_R, pres_R, rho_R, gamma_R, pi_inf_R, qv_R, vel_R_rms, H_R)
+                            call s_compute_speed_of_sound(eos_s_R, alpha_R, c_R)
 
                             if (mhd) then
                                 call s_compute_fast_magnetosonic_speed(rho_L, c_L, B%L, norm_dir, c_fast%L, H_L)
diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp
index daf1c2a73..aa685807d 100644
--- a/src/simulation/m_time_steppers.fpp
+++ b/src/simulation/m_time_steppers.fpp
@@ -659,14 +659,15 @@ contains
         real(wp)               :: dt_local
         integer                :: j, k, l  !< Generic loop iterators
         integer                :: fl       !< Fluid loop iterator
+        type(eos_state)        :: eos_s
 
         if (.not. igr) then
             call s_convert_conservative_to_primitive_variables(q_cons_ts(1)%vf, q_T_sf, q_prim_vf, idwint)
         end if
 
         dt_local = huge(1.0_wp)
-        $:GPU_PARALLEL_LOOP(collapse=3, private='[vel, alpha, Re, rho, vel_sum, pres, gamma, pi_inf, c, H, qv, fl, max_dt]', &
-                            & reduction='[[dt_local]]', reductionOp='[min]')
+        $:GPU_PARALLEL_LOOP(collapse=3, private='[vel, alpha, Re, rho, vel_sum, pres, gamma, pi_inf, c, H, qv, fl, max_dt, &
+                            & eos_s]', reduction='[[dt_local]]', reductionOp='[min]')
         do l = 0, p
             do k = 0, n
                 do j = 0, m
@@ -677,7 +678,8 @@ contains
                     end if
 
                     ! Compute mixture sound speed
-                    call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, alpha, vel_sum, 0._wp, c, qv)
+                    call s_eos_state_roe(eos_s, pres, rho, gamma, pi_inf, qv, vel_sum, H)
+                    call s_compute_speed_of_sound(eos_s, alpha, c)
 
                     if (any_non_newtonian) then
                         Re(1) = 0._wp

From 102c54c7b13ec42f83577aeec9be99d47332f833 Mon Sep 17 00:00:00 2001
From: Spencer Bryngelson <sbryngelson@gmail.com>
Date: Sat, 8 Aug 2026 22:27:22 -0400
Subject: [PATCH 2/3] Build hot-path eos_states with the intrinsic structure
 constructor

Profiling the +5-8% regression on the Riemann benchmarks showed it was not
inlining loss: both before and after, s_compute_speed_of_sound is a real
non-inlined call. The cost was the constructor call added on top of it.

Object-code counts for m_riemann_solver_hllc (arm64, GNU 15.2):

    baseline                bl 179, relocs speed_of_sound 27, eos_state 0
    s_eos_state_roe call    bl 206, relocs speed_of_sound 27, eos_state 27
    structure constructor   bl 179, relocs speed_of_sound 27, eos_state 0

Nine call sites x 3 relocations each = the 27 added calls, exactly one per
site. Replacing the constructor call at the 16 Riemann sites with Fortran's
intrinsic structure constructor builds the state in place and restores the
baseline call count.

The derived-type constructors remain for the diagnostic sites, where deriving
H from the other members is the point and the call cost is irrelevant.

Benchmarks (idle machine, mean of two baselines): hll -0.12%, lf -2.12%,
hypo_hll -1.00%, ibm -0.57%, igr -0.09%, viscous -1.26%, all within the
measured +/-3.5% baseline noise band. hllc remains +5.63%, consistent across
two independent runs and above both baseline samples; it has nine call sites,
the most of any solver, so the residual is attributed to materialising the
state rather than to call overhead.

Full suite 627 passed, 0 failed.
---
 src/simulation/m_riemann_solver_hll.fpp  |  7 +++----
 src/simulation/m_riemann_solver_hllc.fpp | 20 ++++++++++----------
 src/simulation/m_riemann_solver_hlld.fpp |  4 ++--
 src/simulation/m_riemann_solver_lf.fpp   |  4 ++--
 4 files changed, 17 insertions(+), 18 deletions(-)

diff --git a/src/simulation/m_riemann_solver_hll.fpp b/src/simulation/m_riemann_solver_hll.fpp
index c89ce76fc..6ad0ebb56 100644
--- a/src/simulation/m_riemann_solver_hll.fpp
+++ b/src/simulation/m_riemann_solver_hll.fpp
@@ -319,17 +319,16 @@ contains
 
                             @:compute_average_state()
 
-                            call s_eos_state_roe(eos_s_L, pres_L, rho_L, gamma_L, pi_inf_L, qv_L, vel_L_rms, H_L)
+                            eos_s_L = eos_state(rho_L, pres_L, gamma_L, pi_inf_L, qv_L, vel_L_rms, H_L, 0._wp)
                             call s_compute_speed_of_sound(eos_s_L, alpha_L, c_L)
 
-                            call s_eos_state_roe(eos_s_R, pres_R, rho_R, gamma_R, pi_inf_R, qv_R, vel_R_rms, H_R)
+                            eos_s_R = eos_state(rho_R, pres_R, gamma_R, pi_inf_R, qv_R, vel_R_rms, H_R, 0._wp)
                             call s_compute_speed_of_sound(eos_s_R, alpha_R, c_R)
 
                             !> The computation of c_avg does not require all the variables, and therefore the non '_avg'
                             ! variables are placeholders to call the subroutine.
 
-                            call s_eos_state_roe(eos_s_avg, pres_R, rho_avg, gamma_avg, pi_inf_R, qv_avg, vel_avg_rms, H_avg, &
-                                                 & c_sum_Yi_Phi)
+                            eos_s_avg = eos_state(rho_avg, pres_R, gamma_avg, pi_inf_R, qv_avg, vel_avg_rms, H_avg, c_sum_Yi_Phi)
                             call s_compute_speed_of_sound(eos_s_avg, alpha_R, c_avg)
 
                             if (mhd) then
diff --git a/src/simulation/m_riemann_solver_hllc.fpp b/src/simulation/m_riemann_solver_hllc.fpp
index 561f83d48..2f2d6be3b 100644
--- a/src/simulation/m_riemann_solver_hllc.fpp
+++ b/src/simulation/m_riemann_solver_hllc.fpp
@@ -232,15 +232,15 @@ contains
 
                                 @:compute_average_state()
 
-                                call s_eos_state_roe(eos_s_L, pres_L, rho_L, gamma_L, pi_inf_L, qv_L, vel_L_rms, H_L)
+                                eos_s_L = eos_state(rho_L, pres_L, gamma_L, pi_inf_L, qv_L, vel_L_rms, H_L, 0._wp)
                                 call s_compute_speed_of_sound(eos_s_L, alpha_L, c_L)
 
-                                call s_eos_state_roe(eos_s_R, pres_R, rho_R, gamma_R, pi_inf_R, qv_R, vel_R_rms, H_R)
+                                eos_s_R = eos_state(rho_R, pres_R, gamma_R, pi_inf_R, qv_R, vel_R_rms, H_R, 0._wp)
                                 call s_compute_speed_of_sound(eos_s_R, alpha_R, c_R)
 
                                 !> The computation of c_avg does not require all the variables, and therefore the non '_avg'
                                 ! variables are placeholders to call the subroutine.
-                                call s_eos_state_roe(eos_s_avg, pres_R, rho_avg, gamma_avg, pi_inf_R, qv_avg, vel_avg_rms, H_avg)
+                                eos_s_avg = eos_state(rho_avg, pres_R, gamma_avg, pi_inf_R, qv_avg, vel_avg_rms, H_avg, 0._wp)
                                 call s_compute_speed_of_sound(eos_s_avg, alpha_R, c_avg)
 
                                 if (viscous) then
@@ -604,15 +604,15 @@ contains
                                     end do
                                 end if
 
-                                call s_eos_state_roe(eos_s_L, pres_L, rho_L, gamma_L, pi_inf_L, qv_L, vel_L_rms, H_L)
+                                eos_s_L = eos_state(rho_L, pres_L, gamma_L, pi_inf_L, qv_L, vel_L_rms, H_L, 0._wp)
                                 call s_compute_speed_of_sound(eos_s_L, alpha_L, c_L)
 
-                                call s_eos_state_roe(eos_s_R, pres_R, rho_R, gamma_R, pi_inf_R, qv_R, vel_R_rms, H_R)
+                                eos_s_R = eos_state(rho_R, pres_R, gamma_R, pi_inf_R, qv_R, vel_R_rms, H_R, 0._wp)
                                 call s_compute_speed_of_sound(eos_s_R, alpha_R, c_R)
 
                                 !> The computation of c_avg does not require all the variables, and therefore the non '_avg'
                                 ! variables are placeholders to call the subroutine.
-                                call s_eos_state_roe(eos_s_avg, pres_R, rho_avg, gamma_avg, pi_inf_R, qv_avg, vel_avg_rms, H_avg)
+                                eos_s_avg = eos_state(rho_avg, pres_R, gamma_avg, pi_inf_R, qv_avg, vel_avg_rms, H_avg, 0._wp)
                                 call s_compute_speed_of_sound(eos_s_avg, alpha_R, c_avg)
 
                                 if (viscous) then
@@ -942,16 +942,16 @@ contains
 
                                 @:compute_average_state()
 
-                                call s_eos_state_roe(eos_s_L, pres_L, rho_L, gamma_L, pi_inf_L, qv_L, vel_L_rms, H_L)
+                                eos_s_L = eos_state(rho_L, pres_L, gamma_L, pi_inf_L, qv_L, vel_L_rms, H_L, 0._wp)
                                 call s_compute_speed_of_sound(eos_s_L, alpha_L, c_L)
 
-                                call s_eos_state_roe(eos_s_R, pres_R, rho_R, gamma_R, pi_inf_R, qv_R, vel_R_rms, H_R)
+                                eos_s_R = eos_state(rho_R, pres_R, gamma_R, pi_inf_R, qv_R, vel_R_rms, H_R, 0._wp)
                                 call s_compute_speed_of_sound(eos_s_R, alpha_R, c_R)
 
                                 !> The computation of c_avg does not require all the variables, and therefore the non '_avg'
                                 !  variables are placeholders to call the subroutine.
-                                call s_eos_state_roe(eos_s_avg, pres_R, rho_avg, gamma_avg, pi_inf_R, qv_avg, vel_avg_rms, H_avg, &
-                                                     & c_sum_Yi_Phi)
+                                eos_s_avg = eos_state(rho_avg, pres_R, gamma_avg, pi_inf_R, qv_avg, vel_avg_rms, H_avg, &
+                                                      & c_sum_Yi_Phi)
                                 call s_compute_speed_of_sound(eos_s_avg, alpha_R, c_avg)
 
                                 if (viscous) then
diff --git a/src/simulation/m_riemann_solver_hlld.fpp b/src/simulation/m_riemann_solver_hlld.fpp
index c60398c22..37d1b893e 100644
--- a/src/simulation/m_riemann_solver_hlld.fpp
+++ b/src/simulation/m_riemann_solver_hlld.fpp
@@ -148,9 +148,9 @@ contains
                             H_no_mag%R = (E%R + pres%R - pres_mag%R)/rho%R
 
                             ! (2) Compute fast wave speeds
-                            call s_eos_state_roe(eos_s_L, pres%L, rho%L, gamma%L, pi_inf%L, qv%L, vel_rms%L, H_no_mag%L)
+                            eos_s_L = eos_state(rho%L, pres%L, gamma%L, pi_inf%L, qv%L, vel_rms%L, H_no_mag%L, 0._wp)
                             call s_compute_speed_of_sound(eos_s_L, alpha_L, c%L)
-                            call s_eos_state_roe(eos_s_R, pres%R, rho%R, gamma%R, pi_inf%R, qv%R, vel_rms%R, H_no_mag%R)
+                            eos_s_R = eos_state(rho%R, pres%R, gamma%R, pi_inf%R, qv%R, vel_rms%R, H_no_mag%R, 0._wp)
                             call s_compute_speed_of_sound(eos_s_R, alpha_R, c%R)
                             call s_compute_fast_magnetosonic_speed(rho%L, c%L, B%L, norm_dir, c_fast%L, H_no_mag%L)
                             call s_compute_fast_magnetosonic_speed(rho%R, c%R, B%R, norm_dir, c_fast%R, H_no_mag%R)
diff --git a/src/simulation/m_riemann_solver_lf.fpp b/src/simulation/m_riemann_solver_lf.fpp
index ed70e05ae..9501a8100 100644
--- a/src/simulation/m_riemann_solver_lf.fpp
+++ b/src/simulation/m_riemann_solver_lf.fpp
@@ -292,10 +292,10 @@ contains
                                 H_R = (E_R + pres_R)/rho_R
                             end if
 
-                            call s_eos_state_roe(eos_s_L, pres_L, rho_L, gamma_L, pi_inf_L, qv_L, vel_L_rms, H_L)
+                            eos_s_L = eos_state(rho_L, pres_L, gamma_L, pi_inf_L, qv_L, vel_L_rms, H_L, 0._wp)
                             call s_compute_speed_of_sound(eos_s_L, alpha_L, c_L)
 
-                            call s_eos_state_roe(eos_s_R, pres_R, rho_R, gamma_R, pi_inf_R, qv_R, vel_R_rms, H_R)
+                            eos_s_R = eos_state(rho_R, pres_R, gamma_R, pi_inf_R, qv_R, vel_R_rms, H_R, 0._wp)
                             call s_compute_speed_of_sound(eos_s_R, alpha_R, c_R)
 
                             if (mhd) then

From c3199034245377eaaae92c3c806c046379facbe2 Mon Sep 17 00:00:00 2001
From: Spencer Bryngelson <sbryngelson@gmail.com>
Date: Sun, 9 Aug 2026 20:34:11 -0500
Subject: [PATCH 3/3] Address review: drop the never-passed optional c_c and
 fix eos_state docs

---
 src/common/m_derived_types.fpp        |  4 ++--
 src/common/m_variables_conversion.fpp | 10 ++++------
 2 files changed, 6 insertions(+), 8 deletions(-)

diff --git a/src/common/m_derived_types.fpp b/src/common/m_derived_types.fpp
index 3a45cd2bc..bb01c8a89 100644
--- a/src/common/m_derived_types.fpp
+++ b/src/common/m_derived_types.fpp
@@ -129,8 +129,8 @@ module m_derived_types
     !> num_fluids_max would compile but put ten reals into a per-cell private struct on device.
     !>
     !> `H` is the specific total enthalpy and must include `qv`, because the sound-speed relation
-    !> subtracts `qv/rho`. Build states with f_eos_state so that invariant holds by construction
-    !> rather than by convention; f_eos_state_roe exists for the Roe-averaged paths, which supply an
+    !> subtracts `qv/rho`. Build states with s_eos_state so that invariant holds by construction
+    !> rather than by convention; s_eos_state_roe exists for the Roe-averaged paths, which supply an
     !> `H` that is deliberately not the exact state enthalpy.
     !> Contains no allocatable members - safe to use inside device routines.
     type eos_state
diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp
index ee2c611fe..57f6b2029 100644
--- a/src/common/m_variables_conversion.fpp
+++ b/src/common/m_variables_conversion.fpp
@@ -1163,7 +1163,6 @@ contains
 
     end subroutine s_finalize_variables_conversion_module
 
-    !> Compute the speed of sound from thermodynamic state variables, supporting multiple equation-of-state models.
     !> Build an exact thermodynamic state. The specific total enthalpy is derived here rather than supplied, so it cannot disagree
     !! with qv: H = (E + p)/rho with E = gamma*p + pi_inf + qv + rho*|u|^2/2. Every caller that previously open-coded the closed
     !! form should use this, which makes the defect in #1707 unrepresentable.
@@ -1187,13 +1186,12 @@ contains
 
     !> Build a state whose enthalpy is supplied by the caller. The Roe-averaged Riemann paths, the chemistry Roe branch and the
     !! relativistic branch all pass an H that is deliberately not the exact enthalpy of the state, so they cannot use s_eos_state.
-    subroutine s_eos_state_roe(s, pres, rho, gamma, pi_inf, qv, vel_sum, H, c_c)
+    subroutine s_eos_state_roe(s, pres, rho, gamma, pi_inf, qv, vel_sum, H)
 
         $:GPU_ROUTINE(function_name='s_eos_state_roe', parallelism='[seq]', cray_inline=True)
 
-        type(eos_state), intent(out)   :: s
-        real(wp), intent(in)           :: pres, rho, gamma, pi_inf, qv, vel_sum, H
-        real(wp), intent(in), optional :: c_c
+        type(eos_state), intent(out) :: s
+        real(wp), intent(in)         :: pres, rho, gamma, pi_inf, qv, vel_sum, H
 
         s%pres = pres
         s%rho = rho
@@ -1203,10 +1201,10 @@ contains
         s%vel_sum = vel_sum
         s%H = H
         s%c_c = 0._wp
-        if (present(c_c)) s%c_c = c_c
 
     end subroutine s_eos_state_roe
 
+    !> Compute the speed of sound from thermodynamic state variables, supporting multiple equation-of-state models.
     subroutine s_compute_speed_of_sound(s, adv, c)
 
         $:GPU_ROUTINE(parallelism='[seq]')