[logan-nc]

Fixes #321.

Visual review package (figures, before/after): https://claude.ai/artifact/ETbq28nZE7HmCG7WuRa14Z

Problem

read_kin samples the kinetic table on 101 regular nodes up to ψ_N = 1 and cubic-splines ν_i, ν_e and ln Λ on them. The CI kinetic table ends at ψ_N = 0.998, so the ψ = 1 node is extrapolated to T ≈ 1e-4 eV. There ln Λ = −7.6 and ν_i = −2e14 /s, and the global cubic spline of ν rings inward from that one node. On develop, ν_i first goes negative at ψ = 0.8725 and changes sign 13 times before the edge. The torque density goes negative (anti-damping) from ψ ≈ 0.87 to the edge. A table with T(1) = 0, as an analytic equilibrium would give, is worse: ν = ∞ at the node makes the whole ν spline NaN, and DCON aborts in LSODE at ψ = 0.13.

Change (pentrc/inputs.f90, read_kin only, +22/−5)

• Inside the ν and ln Λ formulas only, floor T at 1 meV and ln Λ at 1, printing a warning when either applies. The stored n and T are not modified.
• After the final spline_fit(kin,"extrap"), replace the node slopes of the ln Λ, ν_i and ν_e columns with pchip slopes (spline_fit_pchip from EQUIL, RDCON, PENTRC - FEATURE - Monotone (pchip) fits for the tabulated Zeff and kinetic profiles #300). spline_eval uses only node values and slopes, so these columns become a local, monotone, non-overshooting interpolant, and a steep node cannot ring inward. The n, T, ω_E and Zeff columns keep their cubic spline.
• error stop on n ≤ 0 or T < 0 at the sampled nodes. These previously gave NaN and an LSODE abort.

The table sampling is unchanged. An earlier draft held the end values beyond the table instead, but that bent T inside the table's last grid interval, so it was dropped. See the design note in the review package.

Results (gfortran, CI Solovev examples)

	Torque [N·m]	δW_k	DCON Im(W)	ψ surfaces
kinetic, develop	5.3732e-3	2.903e-5	−7.27e-3	283
kinetic, this PR	5.7035e-3	2.122e-5	−1.608e-2	81
kinetic, this PR with a T(1)=0 table	5.7037e-3	2.122e-5	−1.608e-2	81
ideal, develop	2.0834e-3	6.093e-6		162
ideal, this PR	2.1478e-3	5.876e-6		112

• n, T, ω_E, ω_N and ω_T are bit-identical to develop, both in the dense spline dump and on all DCON surfaces. ν_i agrees to < 1e-5 for ψ < 0.8 and is positive and monotone through the edge (0 sign changes).
• GPEC outputs are identical, except the kinetic example's gpec_control_output, which follows DCON's kinetic W.
• The CI kinetic example now prints the two floor warnings, because its extrapolated ψ = 1 node is below 1 meV.
• Expect the kinetic-case golden numbers to move. The +6 % torque / −27 % δW_k relative to develop is the ringing being removed, and develop's totals were close only because its positive and negative edge bands partly cancelled. Evidence that the new values are right: the T(1)=0 table, with a different edge node, gives the same torque to 3e-5.

Bad-input behavior

Kinetic table	develop	this PR
CI table (ends at ψ = 0.998)	runs, ν rings	runs
T = 0 at ψ = 1	LSODE abort, ν NaN everywhere	runs (floors warn)
T < 0 at a row	LSODE abort	ERROR: read_kin - temperatures must be non-negative
n = 0 at ψ = 1	LSODE abort	ERROR: read_kin - densities must be positive

Review

Two independent red-team reviews of the slope overwrite (kin%fs1(:,6:8) = tmp_nu%fs1):

• Code paths: every consumer (torque.F90, gpout.f, DCON kin_flag via fourfit.F) evaluates through fs/fs1 only. Nothing refits kin after read_kin, no code takes a derivative of columns 6–8, and OpenMP copyin deep-copies the slopes.
• Numerics: a Fortran test linking the built spline.o reproduces scipy's PchipInterpolator to 1e-15. It shows no overshoot on 1e6 → 1e14 edge profiles (where the cubic dips to −5e12) and leaves the untouched columns bit-identical.

🤖 Generated with Claude Code

https://claude.ai/code/session_01VMDGyfQKRZQMdk3NJXsKQE