Refactor `m_riemann_solvers`: reduce duplication, improve structure, add schemes

[sbryngelson]

Consolidates #63, #529, #612, and #818.

m_riemann_solvers has grown into one of the largest and most tangled modules in MFC. This issue tracks a phased refactor with several independent, actionable items — some are easy wins suitable as first contributions, others are larger architectural changes.

---

1. Easy wins (good first issues)

1a. Merge the bubbles/no-bubbles HLLC branches (#63)

HLLC currently has two near-identical branches for model_eqns == 2 .and. bubbles (L1421–1856) and model_eqns == 2 (L1858–2164). These can be merged into a single branch with if (bubbles_euler) then guards for the bubbles-specific terms, then the first branch deleted.

1b. Extract s_compute_enthalpy helper (#612)

Enthalpy at left/right states is computed identically in at least six places:

• m_riemann_solvers.fpp L361–473, L913–1028, L1270–1328, L1519–1636
• m_cbc.fpp L769–808
• m_data_output.fpp L252–276
• m_start_up.f90 (post-process) L451–452

A single s_compute_enthalpy subroutine eliminates all of this repetition.

---

2. Structural refactor within HLLC (#818)

HLLC has 4 large physics branches (4-Eqn / 5-Eqn / 5-Eqn bubbles / 6-Eqn) with shared scaffolding (populate buffers, initialize, finalize) duplicated across each. The shared structure should be factored out, with each branch reduced to its distinct core logic:

subroutine s_hllc_riemann_solver(...)
    call s_populate_riemann_states_variables_buffers(...)
    call s_initialize_riemann_solver(flux_src_vf, norm_dir)

    if (model_eqns == 3) then
        call s_hllc_6eqns(...)
    elseif (model_eqns == 4) then
        call s_hllc_4eqns(...)
    elseif (model_eqns == 2) then
        if (bubbles_euler) then
            call s_hllc_5eqns_bubbles(...)
        else
            call s_hllc_5eqns_no_bubbles(...)
        end if
    end if

    if (viscous) call s_compute_viscous_source_flux(...)
    if (surface_tension) call s_compute_capillary_source_flux(...)
    call s_finalize_riemann_solver(...)
end subroutine

---

3. Group scalars into derived types (#818)

The solvers currently pass around a large number of independent scalars (rho_L/R, pres_L/R, E_L/R, H_L/R, gamma_L/R, vel(:), alpha(:), …). Grouping these into a simple derived type (no allocatables, so GPU-safe) would make interfaces much cleaner with no expected performance cost:

type riemann_side_state
    real(wp) :: rho, pres, E, H, gamma, pi_inf, qv, vel_rms
    real(wp), dimension(num_vels)   :: vel
    real(wp), dimension(num_fluids) :: alpha, alpha_rho
end type riemann_side_state

type(riemann_side_state) :: L, R
call build_mixture_state(L, qL_prim_rs_vf, j, k, l)
call build_mixture_state(R, qR_prim_rs_vf, j+1, k, l)

Profiling should confirm no regression before merging.

---

4. Separate physical fluxes from numerical fluxes (#529)

The module currently conflates physical flux expressions (e.g., $\rho u^2 + p$) with numerical flux assembly (HLL/HLLC wave-speed weighting). Separating these makes it far easier to add new solvers and reason about correctness.

---

5. Add simpler Riemann solvers (#529)

Once physical/numerical fluxes are separated, adding Lax-Friedrichs (global and local), Rusanov, and potentially Roe becomes straightforward. These existed in earlier archived MFC versions and are useful for numerical verification and testing.

---

6. Cross-solver deduplication (#818)

HLL, HLLC, HLLD, and LF share: left/right state construction, EOS calls, wave speed estimates, and geometric flux treatment. Only the flux assembly kernel differs. Factoring out the shared pieces reduces the total module size substantially. In practice, physics coverage differs per solver (elasticity only in HLL, bubbles only in HLLC, MHD only in HLLD), so cross-solver refactoring should target only the clearly shared scaffolding.

Note on inlining and GPU: Most code is currently inlined for performance. Subroutine extraction should be paired with fypp inlining or verified to be performance-neutral. The xi-based branch-free wave-region selection should be preserved as-is for GPU.

---

Suggested order of attack

1. s_compute_enthalpy extraction (1b) — isolated, testable, easy
2. Merge bubbles branches in HLLC (1a) — well-scoped
3. Scalar → derived type (3) — needs profiling
4. HLLC internal refactor (2)
5. Physical/numerical flux separation (4)
6. Additional solvers (5)
7. Cross-solver deduplication (6)

Closes #63, #529, #612, #818.