trapzd.f

1. (harmless in the runs) The interpolation does not move with the midpoint.

(lines 34–35↗; lines 44–47↗) Inside the loop over midpoints, $\rho$, $\Omega^2$ and $\eta_2+2$ are interpolated at RHOP+DRHO*DEL (and likewise), the same value for every midpoint $y_j$; only $m$ uses $y_j$ (SMT, line 45↗). With FPFT's JMAX = 2 (fpft.f, line 42↗) only the call with N = 2 refines, with IT = 1 and DEL = $b_2 - b_1$, so that the "midpoint" values are those of the outer shell, while $m$ and $r_0$ are taken at the midpoint. What follows: The distortion integral AINT, and through FUNC the gravity $g$ and FPFT's $f_P$, $f_T$, are computed with the wrong integrand at the midpoint; the error is of the order of the change of $\rho$, $\Omega^2$ and $\eta_2$ between neighbouring shells. Only with rotation (LROT, off by default).

A possible fix:

trapzd.f

@@ lines 41–50 @@
        DO 10 J = 1,IT
           R03T = Y**3
 C INTERPOLATE RHO,M,OMEGA,ETA2+2 BETWEEN SHELL I AND SHELL I-1
-          RHOT = RHOP+DRHO*DEL
+C 11/09/26 KC via Claude
+C code changed to interpolate RHO, OMEGA**2 and ETA2+2 at the midpoint
+C Y, as M already is.
+          RHOT = RHOP+DRHO*(Y-B1)
           SMT = SMP+DM*(Y**2 - B1**2)
-          W2T = W2P + DW2*DEL
-          ETA22T = ETA22P + DETA2*DEL
+          W2T = W2P + DW2*(Y-B1)
+          ETA22T = ETA22P + DETA2*(Y-B1)
 C CALCULATE Q BETWEEN SHELLS
           Q0 = (RHOT*W2T*R03T*(3.0D0+ETA22T)/(SMT*ETA22T))*R03T*Y
 C Q0 = RHO*W2*R07T*(3.0D0+ETA22)/(SM*ETA22)

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): The midpoint was evaluated in all five rotating test cases, from about 1.2 to 19.8 million times, and the interpolated density was off by more than 0.1% in most of them (1.1 to 19.7 million times; by up to 100%). With the fix, the 3 $M_\odot$ case (testsuite/Test_3Msolar_GS_rot_gray_OPALSCV_SF2_GS98_OP_CF10_RGp2diffrotI3cz2) stays within the noise: $\log L$ moves by up to 0.009 dex, against 0.0076 for a null change (nearly every opacity, all with $3 \le \log T \le 9$ and $-6 \le \log\rho \le 9.75$, times $1 + 10^{-10}$, on the same case), central $\log T$ by 0.0010 against 0.0008. In the 0.8 $M_\odot$ case (testsuite/Test_m0p8feh+0p0GS98_rotp8p0d0p1mag1e4std_TAMS) $\log L$ and central $\log T$ stay within the noise, and central $X$ and $Y$ moved by up to 0.00015 at 0.395 Gyr, about four and a half times a null change, the same amount at the same age as with SHAPE's fix of $\eta_2$ (SHAPE, bug 1). Rechecked on Nov 10, 2026, with a build without counting lines in SHAPE and TRAPZD, and beside it a null change of the distortion itself (the first shell's $\eta_2$ times $1 + 10^{-10}$): Model by model, every change, the null one included, stays below $10^{-7}$ in central $X$ until about 0.046 Gyr, where central $X$ stops falling, and grows to some $10^{-5}$ by 0.16 Gyr, so this case amplifies whatever perturbs it there, and the size of the change depends on the build. In that build, TRAPZD's fix again moves central $X$ by 0.00015 at 0.3955 Gyr ($1.7\times10^{-4}$ at that one model, against at most $6.1\times10^{-5}$ for the null change, and about $6\times10^{-5}$ against $4\times10^{-5}$ over the main sequence after it), while SHAPE's fix now gives $5.8\times10^{-5}$ at 0.222 Gyr: The shared number was a coincidence of the case's sensitivity, not one fix acting through the other. TRAPZD's change is about three times the null change at one model and of its size elsewhere, and $\log L$ moves by up to 0.00018 dex against 0.00019: not clearly above the noise. The 2 $M_\odot$ cases, which rotate faster, were not run with the fix. Intent: The comment (line 43↗) says the values are interpolated between the shells, as $m$ is, so this looks like a slip.