Replace the enclosing-scope fypp macros in the simulation solvers with explicit interfaces

[sbryngelson]

Nine fypp macros in the simulation solvers expand into the enclosing scope. Most take no arguments at
all: they read and assign local variables of whatever routine expands them, so a call site says nothing
about what is read or written, and neither the compiler nor a reader can check the coupling.

Why this is worth fixing

roe_avg computes vel_avg_rms over every velocity component and then throws it away six lines later:

vel_avg_rms = 0._wp
$:GPU_LOOP(parallelism='[seq]')
do i = 1, num_vels
    vel_avg_rms = vel_avg_rms + (sqrt(rho_L)*vel_L(i) + sqrt(rho_R)*vel_R(i))**2._wp/(sqrt(rho_L) + sqrt(rho_R))**2._wp
end do

H_avg = ...
gamma_avg = ...

vel_avg_rms = (sqrt(rho_L)*vel_L(1) + sqrt(rho_R)*vel_R(1))**2._wp/(sqrt(rho_L) + sqrt(rho_R))**2._wp

The loop is dead and the Roe-averaged vel_avg_rms keeps only the first component. In a subroutine with
intent(out) this is a visible double assignment; inside a 45-line zero-argument macro it is not. It
predates the current work and is not reachable by default (only the pressure-based wave-speed estimate
reads the averaged state), which is why it has gone unnoticed.

compute_low_Mach_correction is the same hazard in the other direction: under low_Mach == 2 it
overwrites vel_L(dir_idx(1)) and vel_R(dir_idx(1)) in place. A call site reading
@:compute_low_Mach_correction() gives no indication that the velocity it is about to use has changed.

compute_hypo_elastic_energy(EL, ER, shear_cond) takes three arguments but still reads i, tau_e_L,
tau_e_R, G_L and G_R implicitly, so the arguments suggest a contract the macro does not keep.

Inventory

macro	sites	file
compute_low_Mach_correction	8	include/inline_riemann.fpp
compute_axis_inv_re	4	m_viscous.fpp
compute_average_state	3	include/inline_riemann.fpp
compute_capillary_stress_tensor	3	include/inline_capillary.fpp
hll_flux_component(LHS, I)	3	m_riemann_solver_hypo_hlld.fpp
compute_elastic_wave_speeds_lr	2	include/inline_riemann.fpp
arithmetic_avg	1	include/inline_riemann.fpp
roe_avg	1	include/inline_riemann.fpp
compute_hypo_elastic_energy	1	include/inline_riemann.fpp

26 call sites.

Proposed replacements

Each becomes a $:GPU_ROUTINE(parallelism='[seq]', cray_inline=True) procedure with explicit intent,
which is the pattern the equation-of-state operators already use (s_compute_mixture_coefficients,
f_bulk_modulus, f_pressure).

! one-line formula, already argument-taking - a function, not a macro
f_hll_flux(S_L, S_R, F_L_i, F_R_i, U_L_i, U_R_i) result(flux)

! the increment, guarded and doubled for shear, rather than two mutated accumulators
f_elastic_energy(tau, G, is_shear) result(dE)

f_elastic_wave_speed(vel, c, G, tau, rho, sgn) result(s)

s_compute_axis_inv_re(grad_x_vf, grad_y_vf, grad_z_vf, alpha_visc, Res_viscous, j, k, l, Re_visc)
s_compute_capillary_stress_tensor(sigma, w1, w2, w3, normW, Omega)

! the mutation becomes visible in the signature
s_compute_low_Mach_correction(c_L, c_R, rho_L, rho_R, vel_L_rms, vel_R_rms, s_L, s_R, s_M, s_P, &
                              vel_L, vel_R, zcoef, pcorr)   ! vel_L, vel_R intent(inout)

s_compute_average_state(rho_L, rho_R, vel_L, vel_R, H_L, H_R, gamma_L, gamma_R, qv_L, qv_R, &
                        rho_avg, vel_avg_rms, H_avg, gamma_avg, qv_avg)

compute_average_state is the widest and should go last. Its chemistry path additionally produces
Cp_avg, Cv_avg, T_avg, Yi_avg, Phi_avg and c_sum_Yi_Phi; splitting the chemistry branch into
its own routine keeps both signatures reasonable.

Whether the dead vel_avg_rms loop should be deleted or the overwrite removed is a physics question,
not a refactoring one, and should be settled separately - the Roe average of |u|^2 is normally the sum
over components, which suggests the overwrite is the error rather than the loop.

Constraint

These are the hottest loops in the code, and #1714 showed the failure mode: a derived-type dummy on a
[seq] device routine forces a memory ABI and collapsed register promotion across the whole enclosing
loop body, costing 20% on amdflang. s_compute_axis_inv_re would take grad_x_vf and friends, which is
exactly that shape.

Scalar-only [seq] operators have since been measured as free: converting the bulk modulus and pressure
inversion to f_bulk_modulus/f_pressure gave 0 regression(s) across 471 kernels in static GPU
resources and no wall-clock change on 5eq_rk3_weno3_hllc (MI210, amdflang, OpenMP offload).

So each conversion should be landed with a static resource diff of the offload image, not assumed
neutral - and the derived-type ones should be measured before the scalar ones are used to justify them.