PyREC scaffolding (from YREC development, 7be4133)/ hpoint.fGlossary

hpoint.f

878 lines · 1 definitions · on GitHub (7be4133)

Python version: hpoint.py, 731 lines, on the right (switch to Python above the code).

Variables on this page (68)

alphae, amuo, cc23, cgl, cln, clndp, cp, cpo, delao, delo, delro, drate0, epso, esumo, eta2x, hamu, hcompp, hdo, hgo, hio, hjx, hll0, hlo, hlsuro, hpo, hpp0, hpttol, hro, hrr0, hss, hto, htt0, idebug, iowr, ishort, json, lco, lczo, ldebug, lexcom, lovste, lrot, npoint, oo, qdto, qwrmax, qwrmx0, qwrst, r0x, sdel, seg, senv, sesum, so, sqdt, svel, tefflo, tenv0, tenv1, thdif, thdifo, tho, velo, visc, visco, wold, xnew, znew

Purpose

HPOINT rezones the model: at the start of each step it chooses a new set of mass points so that the jumps between neighbours in mass, luminosity, pressure, $X$ and $Z$ stay below user limits, keeps real discontinuities (convection-zone edges, composition jumps) as fixed points, and interpolates every model variable to the new points. It also adds or removes outer points when the outermost point drifts out of the allowed envelope temperature range, and recomputes the shell masses and, for rotating models, the rotation quantities.

How it fits

  • Who calls it: the main program, main.f, once per evolutionary step.
  • Inputs (arguments, changed in place): the model (HS = $\log m$, HP, HT, HR, HL, HD, HCOMP, LC, M), HSTOT ($\log M_*$), the envelope-fitting values PS, RS, ESTORE, ISTORE, the indices JCORE, JENV (edges of the convective core and envelope), JXBEG (H-shell base), LSHELL, the model number MODEL, the rotation arrays (OMEGA, ETA2, R0, HI, HJM, HKEROT, FP, FT, HG, QIW), the totals SJTOT, SKEROT, TEFFL, and CMIXL (the run's mixing-length parameter, which main.f sets from CMIXLA(NK); passed on to PHYSIC).
  • Outputs: the rezoned model; HS1 (mass at each point) and HS2 (mass of each shell); LRESET (true when outer points were added or removed); and in modules, the previous model on the new grid for the next step.
  • Calls: SPLINC and SPLINTD2 (splines to place the points), OSPLIN (an "oscillatory" spline for interpolating the variables); for rotation OVROT, MOMI, GETROT, FPFT; PHYSIC (the physics at the new points); SURFOPAC (a new surface opacity table when surface $X$ changed).
  • Modules, the main ones: NMLPAR_MOD (the namelist parameters: HPTTOL(1..12), the zoning limits, from PARMIN; TENV0, TENV1; LEXCOM, whether light elements 12–15 are followed; LROT; the overshoot settings LOVSTE, ALPHAE; the debug output LDEBUG, NPOINT); OLDMOD_MOD (after the rezoning, HSS, HPO, HTO, HRO, HLO, HDO and HCOMPP hold the start-of-step model on the new grid, and TEFFLO, HLSURO the start-of-step $T_\mathrm{eff}$ and surface luminosity, saved before the rezoning; read by KEMCOM, LIBURN, LIRATE88, CRRECT and others); ENTROP_MOD (the previous step's rates of change of $\log T$, $\log P$, $\ln L$, $\log r$, which main.f uses for its first guess of the new structure; shared with main.f only); DEUTER_MOD (DRATE0); OLDROT_MOD, OLDROT2_MOD, OLDPHY_MOD, DWMAX_MOD (rotation and the physics at the start of the step, for the angular-momentum transport in GETW and MIDMOD); SCRTCH_MOD, TEMP_MOD (the physics just computed by PHYSIC); COMP_MOD (XNEW, ZNEW, SENV); CONSTANTS_MOD; IOUNITS_MOD (ISHORT, IOWR, IDEBUG); JSON from DIMENSIONS_MOD; 14 modules in all.

The method

  1. Outer boundary. If $\log T$ at the outermost point is below TENV0, drop the outer points down to the first one hotter than TENV0 (and put the dropped mass into the envelope, SENV). If it is above TENV1, add one point from the stored envelope integration: $\Delta\log P$ from ESTORE, $\Delta\log T = 0.25\,\Delta\log P$, $\Delta\log r$ from ESTORE, the same $L$ and composition, and $\log\rho = \log P - \log T - 8$.

  2. Fixed points. Flag the edge of the convective core, the base of the convective envelope, the base of the H shell, every point where $|\Delta X|$ > HPTTOL(3), $|\Delta Z|$ > HPTTOL(4), or (rotating models) $|\Delta\log\Omega|$ > HPTTOL(12), the base of envelope overshoot and the base of a finely zoned region around the envelope's base (HPTTOL(7) in $\log P$). At most 100 flags; sort them and remove duplicates.

  3. New mass points. Build a cumulative "distance" along the star, $$\chi_j = \chi_{j-1} + \frac{\Delta\log m}{\delta_m} + \frac{\Delta L}{\delta_L} - \frac{\Delta\log P}{\delta_P},$$ with $\delta_m$ = HPTTOL(2), $\delta_L$ = HPTTOL(9) times the largest $L$, and $\delta_P$ = HPTTOL(10) in the finely zoned region, HPTTOL(11) below it, HPTTOL(8) in the convective envelope (the $L$ term only where $L$ increases outward). Between consecutive flagged points, place points evenly in $\chi$, one per unit, and find their masses from a spline of $\log m$ against $\chi$. The idea: $\chi$ adds up the changes of $\log m$, $L$ and $\log P$ along the star, each in units of its tolerance, so that between neighbouring new points their sum is one tolerance; points crowd where any of the three changes fast. The code also sets PMAX1–PMAX5 (by their comment, the largest $\Delta\log P$ in each region), HPMAX(1), IIFLG, I0, and, for rotating models, $\log\Omega$ in OMEGAL, and declares ID(7) with its DATA; nothing reads them (unused?), and the placement uses HPTTOL directly (lines 229–239↗; lines 325–333↗; line 353↗; lines 359–361↗; line 128↗; line 137↗).

  4. Refinement in X and Z. Where $X$ increases outward between new points by more than HPTTOL(5), or $Z$ decreases outward by more than HPTTOL(6), insert evenly spaced points, stopping the run (STOP 1) if they would take the model past JSON points. Then drop points closer than HPTTOL(1) in $\log m$.

  5. Convection flags. Relocate JCORE and JENV on the new grid and reset LC.

  6. Interpolation. With OSPLIN, interpolate to the new points: the entropy-term variables (ENTROP), the composition HCOMP and the starting composition HCOMPP, the deuterium burning rate DRATE0, then $\log P$, $\log r$, $L$, $\log T$, $\log\rho$, and for rotating models $\Omega$, $j$, $f_P$, $f_T$, $r_0$, $\eta_2$; and the energy generation at the start of the step when needed (see bug 1).

  7. Masses. With $m_i = 10^{\log m_i}$: $$\texttt{HS1}_i = m_i,\qquad \texttt{HS2}_i = \tfrac12\,(m_{i+1} - m_{i-1}),\qquad \texttt{HS2}_1 = \tfrac12(m_2 + m_1),\qquad \texttt{HS2}_M = M_* - \tfrac12(m_{M-1}+m_M).$$

  8. Rotation. First guess of the moment of inertia $I_i = \tfrac23\,\Delta m_i\,r_i^2$; overshoot zones (OVROT); moment of inertia (MOMI); $\Omega$ from the conserved specific angular momentum (GETROT); the distortion terms (FPFT); MOMI again; then the totals $J = \sum j_i\Delta m_i$ and $K = \sum \tfrac12\Omega_i j_i \Delta m_i$, printed against the previous ones; and the old rotation state saved in OLDROT.

  9. Physics at the start of the step (rotating models, or light elements with envelope overshoot): PHYSIC at the new points; the results stored in OLDPHY, DWMAX, and the shear term $$\theta_i = \frac{\bar r^{\,4}}{G\,\bar m}\,\bar\Omega\,\frac{\Delta\Omega}{\Delta r}$$ (THO, with geometric means of $r$ and $m$ over the two points) and $d\Omega/dr$ (QWRST) in OLDROT2.

  10. A new surface opacity table (SURFOPAC) if surface $X$ changed by more than $10^{-8}$.

Design

It prints the flag points and the old and new numbers of points to the short output every step, and a full table of old and new points to the debug unit every NPOINT models when LDEBUG is on. Several module arrays of the old model (HPO, HTO, HDO) serve as scratch space before receiving their final contents.

From YREC's documentation

Physics Options (Rarely used options) gives TENV0 and TENV1 (defaults 3.0 and 9.0, in $\log T$) as the minimum and maximum temperature of the envelope integration, and says "Code stops below this" and "Code stops above this." The mesh tolerances HPTTOL that HPOINT applies are documented under their newer names on REMAP's page.

Suspected bugs, checked by running the code

Each bug is marked by what running the code showed: (changes results, large) or (small), the fix moving a run clearly beyond a null change; (stops the run); (wrong output), the fix changing only what is written; (harmless in the runs), within the noise; (never occurred), the code ran but the fault did not; (not reached); (cannot be reached), a guard prevents it; (withdrawn), the suspicion was wrong.

1. (harmless in the runs) The start-of-step energy generation is interpolated against the new mass points.

(line 602↗; lines 733–740↗; lines 751–758↗) The block "MHP 6/00 interpolated in energy generation at start of timestep" runs after loop 1000 has copied the new points HSS into HS, so OSPLIN fits SESUM and SESUM + SEG(6,·) + SEG(7,·), which are on the old points, against HS(1..MNEW) (the new points) and, when the old model had more points, the old HS(MNEW+1..M) after them, whose masses are then not increasing. ESUMO and EPSO thus hold roughly the old values copied index by index, shifted wherever points were added or removed. What follows: MIDMOD interpolates ESUMM and EPSM from them (midmod.f, line 159↗; midmod.f, lines 166–168↗), which SETUPV (setupv.f, lines 155–156↗) and the rotational mixing grids (mixgrid.f, line 156↗; rotgrid.f, line 251↗) use, in rotating runs (LROT, off by default) and with LEXCOM and LOVSTE.

A possible fix (the block moved before loop 1000):

hpoint.f

@@ lines 729–734 @@
       ENDIF


+C 11/09/26 KC via Claude
+C code changed to interpolate the energy generation at the start of the
+C timestep before HS takes the new points, so that SESUM is fitted
+C against its own (old) mass points.
+C MHP 6/00 INTERPOLATED IN ENERGY GENERATION AT START OF TIMESTEP
+      IF(LROT .OR. (LEXCOM .AND. LOVSTE))THEN
+         CALL OSPLIN(HSS,ESUMO,HS,SESUM,NTAB,NTOT)
+         DO IM = 1,M
+            YTAB(IM) = SESUM(IM)+SEG(6,IM)+SEG(7,IM)
+         END DO
+         CALL OSPLIN(HSS,EPSO,HS,YTAB,NTAB,NTOT)
+      ENDIF
 C TRANSFER NEW POINTS.
       DO 1000 J = 1,MNEW
        HS(J) = HSS(J)
@@ lines 748–761 @@
           R0(J) = R0X(J)
  1005    CONTINUE
       ENDIF
-C MHP 6/00 INTERPOLATED IN ENERGY GENERATION AT START OF TIMESTEP
-      IF(LROT .OR. (LEXCOM .AND. LOVSTE))THEN
-         CALL OSPLIN(HSS,ESUMO,HS,SESUM,NTAB,NTOT)
-         DO IM = 1,M
-            YTAB(IM) = SESUM(IM)+SEG(6,IM)+SEG(7,IM)
-         END DO
-         CALL OSPLIN(HSS,EPSO,HS,YTAB,NTAB,NTOT)
-      ENDIF
       WRITE(ISHORT,1020) M,MNEW
  1020 FORMAT(' POINTS  OLD',I5,'   NEW',I5)
       M = MNEW

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): taken at every rezoning in 7 of 34 runs (26 of YREC's cases, from 0.04 to 9 $M_\odot$, with rotation, settling, and helium burning, and 8 made-up ones): the five rotating ones (testsuite/Test_solarGS98_difrotmixper18p4d1p0, testsuite/Test_m0p8feh+0p0GS98_rotp8p0d0p1mag1e4std_TAMS, testsuite/Test_m2p0_feh+0p0_vzams150_IES1_TAMS and _IES2_TAMS, and testsuite/Test_3Msolar_GS_rot_gray_OPALSCV_SF2_GS98_OP_CF10_RGp2diffrotI3cz2), examples/run_from_zahb_to_tahb/ZAHB_m0p8feh+0p0_TAHB (LOVSTE), and a made-up helium-burning case (examples/run_from_zahb_to_tahb/ZAHB_m1p0feh+0p0_TAHB with the AAG21 opacity table), 546 to 2,250 times per run. Every time, the old and new points differed (models of 350 to 3,455 points). With the fix, the 3 $M_\odot$ giant's track is identical, and the 0.8 $M_\odot$ track moves by up to $6.4\times10^{-5}$ dex in $\log L$, against $1.4\times10^{-4}$ 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}$): within the noise, with the same number of models (1,853 and 763). The fix was not run on the 2 $M_\odot$, solar, or helium-burning cases, where these rates may weigh more. Intent: none found.

2. (never occurred) The flag array can overflow.

(line 73↗; lines 223–227↗; lines 264–267↗; lines 287–290↗; line 294↗) NFLAG has 100 elements, and the gradient tests stop at IFLAG = 100; but when they end with IFLAG = 99, the base of the overshoot region and the base of the finely zoned region are flagged after them, and label 130 writes NFLAG(101) = M, past the array, which the sort then reads. This needs LOVSTE (off by default) and a nonzero HPTTOL(7) (zero by default, but every case of YREC sets it through LNEWVARS, with TOL_CZBASE_FINE_WIDTH from 0.1 to 0.4: remap.f, line 249↗), a convective envelope, and exactly 98 flags set before them.

A possible fix:

hpoint.f

@@ lines 261–267 @@
           ELSE
              PMAX3 = HPTTOL(10)
           ENDIF
-            IF(JOS.GT.1)THEN
+C 11/09/26 KC via Claude
+C code changed to add the overshoot and fine-zoning flags only while
+C NFLAG has room for them and for M.
+            IF(JOS.GT.1 .AND. IFLAG.LT.100)THEN
                NFLAG(IFLAG) = JOS
                IFLAG = IFLAG + 1
             ENDIF
@@ lines 284–290 @@
           IF(MOD(DELP,HPTTOL(10)).NE.0D0)NP=NP+1
           PMAX2 = DELP/DFLOAT(NP)
           IF(PMAX2.EQ.0.0D0)PMAX2 = HPTTOL(10)
-            IF(JFINE.GT.1)THEN
+            IF(JFINE.GT.1 .AND. IFLAG.LT.100)THEN
                NFLAG(IFLAG) = JFINE
                IFLAG = IFLAG + 1
             ENDIF

Checked by running YREC (Nov 9, 2026, the counts of testsuite/Test_m1p0feh+0p0GS98_base_HeIgnite repeated on Nov 10 with a private start model; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): The base of the finely zoned region was flagged in 26 of the 34 runs, up to 7,196 times per run (testsuite/Test_m1p0feh+0p0GS98_base_HeIgnite), but with at most 21 flags set before it, far from the 98 needed; NFLAG(101) was never written. Not reached: It would take a run with LOVSTE and about a hundred jumps in $X$, $Z$, or $\Omega$ to flag. Intent: The message "MORE THAN 100 FLAG POINTS-FIRST 100 RETAINED" (line 225↗) shows that 100 flags were meant to be the limit.

3. (changes results, large) Checked by KC (Nov 10, 2026) ALPHAE taken as a depth in $\log_{10}P$.

(lines 252–255↗; main.f, lines 803–804↗) The base of envelope overshoot is the first point where $\log_{10}P - \log_{10}P(\texttt{JENV}) >$ ALPHAE. The documentation gives ALPHAE as the "envelope overshoot in pressure scale heights," and main.f uses it so (PSCAL0 = ALPHAE $\cdot H_p$); an overshoot of $\alpha H_p$ is about $\Delta\ln P = \alpha$, that is $\Delta\log_{10}P \approx 0.43\,\alpha$. The point flagged as the base of the overshoot region is thus about 2.3 times deeper in $\log P$ than the overshoot, where the mesh keeps a fixed point (PMAX3, computed from it, is unused); the run below shows that this changes the track. With LOVSTE (off by default).

A possible fix:

hpoint.f

@@ lines 61–66 @@
      *       HT,OMEGA,OMEGAL,PMAX1,PMAX2,PMAX3,PMAX4,PMAX5,PS,QIW,R0,RS,
      *       SJTOT,SKEROT,SUMJ,SUMKE,TEFFL,TESTP,WLBOT,WLTOP,WMID,X,
      *       XTAB,Y,Y2,YTAB
+C 11/10/26 KC via Claude
+C code changed to declare the scale height and radius used to find
+C the base of the envelope overshoot.
+      REAL*8 PSCALD,ROS
       INTEGER I,I0,ID,IEND,IFLAG,IIFLG,IM,ISTORE,IXFLAG,IXX,J,JBEG,
      *       JCORE,JEND,JENV,JFINE,JJ,JOS,JPREV,JXBEG,K,KHI,KLO,M,
      *       MINMIN,MM,MNEW,MODEL,MRZONE,MXZONE,NFLAG,NN,NP,NPT,NRZONE,
@@ lines 249–256 @@
           I = JENV
           JOS = JENV
        ELSE
+C 11/10/26 KC via Claude
+C code changed to find the base of the overshoot region as OVERSH sets
+C it: the first point more than ALPHAE pressure scale heights (from
+C HSUBP, at the base of the convection zone) below JENV in radius,
+C not ALPHAE in log P.
+          CALL HSUBP(HCOMP,HD,HP,HR,HS,HT,JENV,PSCALD)
+          ROS = DEXP(CLN*HR(JENV)) - ALPHAE*PSCALD
           DO 191 JOS = JENV-1,1,-1
-             IF(HP(JOS)-HP(JENV).GT.ALPHAE)GOTO 193
+             IF(DEXP(CLN*HR(JOS)).LT.ROS)GOTO 193
   191       CONTINUE
   193       JOS = JOS + 1
           DELP = HP(JOS)-HP(JENV)

How the code itself shows it. Every other line that reads ALPHAE takes it as a length in pressure scale heights. OVERSH, which sets the mixed region whose base this point marks, says so, "ALPHAE IS THE DESIRED OVERSHOOT (IN SCALE HEIGHTS).", and multiplies the scale height from HSUBP by it, PSCALD = PSCALD*ALPHAE (oversh.f, lines 102–104↗); it then moves the edge of the mixed region down by that radial distance (oversh.f, lines 124–130↗). main.f forms the same length for light-element burning, PSCAL0 = ALPHAE*EXP(CLNDP*(HP(JENV)+2.0D0*HR(JENV)-HD(JENV)-CGL-HS(JENV))), that is $\alpha_e P r^2/(\rho G m)$ (main.f, lines 803–804↗, also main.f, lines 1025–1026↗), and GETW the same with its midpoint arrays (getw.f, lines 390–391↗), and the documentation reads "Envelope overshoot in pressure scale heights" (docs/source/namelist_physics.md, line 69). Over a distance $\alpha_e H_p$, $\ln P$ changes by about $\alpha_e$; HP holds $\log_{10}P$, so the matching difference is about ALPHAE/CLN, with CLN = CLNDP = DLOG(10.0D0) (setups.f, lines 134–135↗), the version the run below measured; it is exact only if $H_p$ were constant over the region, and $H_p$ shrinks inward. The fix above therefore uses OVERSH's own criterion, the radius $\alpha_e H_p$ below the base of the zone with $H_p$ from HSUBP there (hsubp.f, line 75↗), so that HPOINT's point falls at the edge OVERSH mixes to. That is the only overshoot scheme YREC has for the envelope; core overshoot can also be limited to a fraction BETAC of the core's radius (LOVMAX, oversh.f, lines 86–89↗), which does not apply here.

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): taken at every model of the two runs with envelope overshoot (ALPHAE = 0.5), examples/run_from_zahb_to_tahb/ZAHB_m0p8feh+0p0_TAHB (1,803 times) and the made-up examples/run_from_zahb_to_tahb/ZAHB_m1p0feh+0p0_TAHB with the AAG21 opacity table, also run for suspected bug 1 (2,250), and every time the base found with $\alpha_e/\ln 10$ was another point (the points flagged lay between 1,654 and 3,155, the corrected ones between 1,684 and 3,170). With the fix as first written (ALPHAE/CLN in place of ALPHAE, an approximate version of the fix above), ZAHB_m0p8feh+0p0_TAHB changes well above the noise: $\log L$ by up to 0.334 dex, against $2.0\times10^{-7}$ for a null change; $\log T_\mathrm{eff}$ by 0.0356, against $3.3\times10^{-7}$; $\log T_c$ by 0.0752, against $1.6\times10^{-7}$; and central $Y$ by up to 0.279, against $6.9\times10^{-8}$. The run ends after 1,855 models at 0.2356 Gyr instead of 1,804 at 0.2297. The fix acts through the mesh alone, so this also shows how sensitive helium burning is to the zoning. Intent: The documentation gives ALPHAE in pressure scale heights, and main.f uses it so (main.f, lines 803–804↗).

4. (not reached) The debug table prints new values under "old points."

(lines 436–441↗; lines 493–499↗; lines 648–655↗; lines 713–720↗) Under the heading "OLD POINTS," the columns X, Z, and O16 print HKO, HFO, and HCOMP(9,·): HKO and HFO hold $X$ and $Z$ interpolated at the new points (before the $X$ and $Z$ insertions of step 4), and HCOMP has already been moved to the new points (loop 850). The comments also say these values go into HIO and HGO. Only with LDEBUG (off by default).

A possible fix:

hpoint.f

@@ lines 645–650 @@
   914 CONTINUE


+C 11/09/26 KC via Claude
+C code changed to keep the old X, Z and O16 for the debug table, which
+C printed values at the new points under OLD POINTS.
+      DO J = 1,M
+         HKO(J) = HCOMP(1,J)
+         HFO(J) = HCOMP(3,J)
+         YTAB(J) = HCOMP(9,J)
+      END DO
       DO 850 I = 1,IEND
        DO 833 J = 1,M
           HPO(J) = HCOMP(I,J)
@@ lines 711–723 @@
   910    FORMAT('1',20X,'OLD POINTS',54X,'NEW POINTS'/2(3X,'N',5X,'S',
      *   8X,'P',7X,'T',7X,'R',8X,'L',7X,'X',4X,'Z',3X,'O16',1X) )
          WRITE(IDEBUG,920) (I,HS(I),HP(I),HT(I),HR(I),HL(I),
-     *   HKO(I),HFO(I),HCOMP(9,I),I,HSS(I),HPO(I),HTO(I),
+     *   HKO(I),HFO(I),YTAB(I),I,HSS(I),HPO(I),HTO(I),
      *   HRO(I),HLO(I),HCOMP(1,I),HCOMP(3,I),HCOMP(9,I), I = 1,MINMIN)
   920    FORMAT( 2(1X,I3,F11.7,F8.4,F8.5,F8.4,1PE9.2,0PF6.3,2F5.3) )
          IF(M.GT.MINMIN) THEN
             MINMIN = MINMIN + 1
             WRITE(IDEBUG,930) (I,HS(I),HP(I),HT(I),HR(I),HL(I),
-     *      HKO(I),HFO(I),HCOMP(9,I),I=MINMIN,M)
+     *      HKO(I),HFO(I),YTAB(I),I=MINMIN,M)
   930       FORMAT( 1X,I3,F11.7,F8.4,F8.5,F8.4,1PE9.2,0PF6.3,2F5.3)
          ELSE IF(MNEW.GT.MINMIN) THEN
             MINMIN = MINMIN + 1

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): not reached in any of the 34 runs; it needs LDEBUG with NPOINT below 9999, which no case sets. Intent: HKO and HFO are reused as scratch arrays ("USING HIO AS DUMMY ARRAY," line 448↗), so the columns look like a side effect of that reuse.

5. (never occurred) A third equal flag survives.

(lines 309–320↗) After removing a duplicate, the loop still advances I, so the entry moved into place is not compared with the one before it, and of three equal flags in a row two remain. The duplicate makes the placement add a point at almost the same $\chi$, hence the same mass (lines 398–419↗), which the removal of points closer than HPTTOL(1) ($10^{-8}$ in $\log m$ by default) then deletes (lines 542–556↗); so this is correct today, but a zero HPTTOL(1) would leave two points at the same mass.

A possible fix:

hpoint.f

@@ lines 306–311 @@
   150 CONTINUE
       IF(.NOT.LDONE) GOTO 140
 C ENSURE THAT POINTS ARENT FLAGGED MORE THAN ONCE.
+C 11/09/26 KC via Claude
+C code changed to compare the entry moved into place with the one
+C before it, so that no flag is kept twice.
       I = 2
   160 CONTINUE
       IF(NFLAG(I).EQ.NFLAG(I-1))THEN
@@ lines 315–322 @@
   170       CONTINUE
        ENDIF
        IFLAG = IFLAG - 1
+      ELSE
+       I = I + 1
       ENDIF
-      I = I + 1
       IF(I.LE.IFLAG) GOTO 160
   180 CONTINUE
       WRITE(ISHORT,185) (NFLAG(J),J=1,IFLAG)

Checked by running YREC (Nov 9, 2026, the counts of testsuite/Test_m1p0feh+0p0GS98_base_HeIgnite repeated on Nov 10 with a private start model; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): A duplicate flag was removed in 6 of the 34 runs, up to 193 times per run (testsuite/Test_m1p0feh+0p0GS98_base_HeIgnite; also the 1.4, 3, and 9 $M_\odot$ cases), and a third equal flag never followed. The fix was not run. Every case sets HPTTOL(1) (as TOL_DM_MIN, $5\times10^{-11}$) above zero, so the removal of close points would still catch a duplicate. Intent: "ENSURE THAT POINTS ARENT FLAGGED MORE THAN ONCE." (line 308↗).

Definitions

subroutineHPOINT

lines 6–878 · GitHub

Rezoning: fixed points, even spacing in the cumulative variable $\chi$, refinement in $X$ and $Z$, interpolation of all variables, new shell masses, and the rotation and start-of-step physics on the new grid (above).

Uses modules COMP_MOD, CONSTANTS_MOD, DEUTER_MOD, DIMENSIONS_MOD, DWMAX_MOD, ENTROP_MOD, IOUNITS_MOD, NMLPAR_MOD, OLDMOD_MOD, OLDPHY_MOD, OLDROT2_MOD, OLDROT_MOD, SCRTCH_MOD, TEMP_MOD

Called by (file level)

Flowchart of its calls, in order (33)
SPLINC394
for I = 1,IFLAG, lines 398–421
for J = JPREV+1,JPREV+NPT, lines 410–417
SPLINTD2412
SPLINC435
for I = 2,MNEW, lines 437–441
SPLINTD2439
SPLINC492
for I = 2,MNEW, lines 494–499
SPLINTD2496
OSPLIN617→OSPLIN628→OSPLIN635→OSPLIN642
for I = 1,IEND (DO 850), lines 648–678
OSPLIN652→OSPLIN661
OSPLIN682→OSPLIN688→OSPLIN689→OSPLIN690→OSPLIN691→OSPLIN692
if (LROT), lines 695–702
OSPLIN696→OSPLIN697→OSPLIN698→OSPLIN699→OSPLIN700→OSPLIN701
if (LROT .OR. (LEXCOM .AND. LOVSTE)), lines 752–758
OSPLIN753→OSPLIN757
if (LROT), lines 774–827
OVROT780→MOMI786→GETROT791→FPFT793→MOMI796
if (LROT .OR. (LEXCOM .AND. LOVSTE)), lines 828–868
PHYSIC834
if (DABS(XNEW-HCOMP(1,M)).GT.1.0D-8), lines 870–875
SURFOPAC873
Compiler warnings (4)

gfortran -Wall -Wextra, compiling only; nothing was run.

  • line 259: Inequality comparison for REAL(8) [-Wcompare-reals]
  • line 269: Equality comparison for REAL(8) [-Wcompare-reals]
  • line 284: Inequality comparison for REAL(8) [-Wcompare-reals]
  • line 286: Equality comparison for REAL(8) [-Wcompare-reals]

Change log

The changes since rp_1

Changed since rp_1. CMIXL is a new argument, passed on to PHYSIC (it was in /CONST3/); COMMON blocks to modules, IMPLICIT NONE; the namelist parameters it reads (from /CTOL/, /CTLIM/, /DPMIX/, /FLAG/, /ROT/, /CCOUT1/, /CCOUT2/ before) now come from NMLPAR_MOD, and unused declarations are gone. No question is settled.

Changed since a8a5950. The $X$ and $Z$ insertions of step 4 stop (STOP 1) before running past JSON (settles former question 1); the previous model's surface luminosity is saved as HLSURO in place of the point count MO, for ENTIME (settles former question 6); DRATE0 is moved onto the new mesh unconditionally, since the old gate (LEXCOM and surface deuterium $\ge 10^{-14}$) read HCOMP(12,M) on the new mesh with the old M and, when it failed, left DRATE0 in old-mesh order; the regrid of HCOMPM, which nothing read, is gone with BURN_MOD (removed since a8a5950); the store of DDELM, which SETUPV overwrites before use, is dropped; and the stop when the whole model is cooler than TENV0 is STOP 1 (KNOWN_ISSUES.md, issue 38).