pindex.f

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

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

Variables on this page (2)

json, nprtpt

Purpose

Chooses which shells of a model are printed (ID(1..IDM)): every NPRTPT-th shell, plus every shell of the hydrogen-burning shell when LSHELL is true. The first and last shells are always included.

How it fits

Inputs: the edges of the burning shell JXBEG, JXEND, the flag LSHELL, the number of shells M; NPRTPT from NMLPAR_MOD (formerly /CCOUT1/; set in PARMIN).
Who calls it: only PUTSTORE, which passes LSHELL, JXBEG and JXEND as locals it never sets (putstore.f, line 262 ↗), so the LSHELL branch is not taken wherever the compiler zero-fills SAVEd locals (settled from the code, Nov 9, 2026).
Outputs: the list ID and its length IDM.

The method

Without a shell: $1, k, 2k, \dots, M$ with $k$ = NPRTPT. With a shell: the multiples of $k$ below JXBEG, every shell from JXBEG to JXEND, the multiples of $k$ above, and $M$.

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. (cannot be reached) `IEND` undefined.

(lines 48–56 ↗) In the LSHELL branch, if NPRTPT $\ge$ JXBEG, the first block is skipped and IF(IEND .EQ. JXBEG) reads IEND from an earlier call (SAVE), or undefined on the first. The result decides whether JXBEG is listed twice or not at all. Only if the LSHELL branch is taken (above).

A possible fix:

pindex.f

@@ lines 41–48 @@

 C  ID IS THE ARRAY OF THE SHELLS TO BE PRINTED OUT;
 C  I.E. IF SHELLS 1,3,5,..TO BE PRINTED ID(1)=1,ID(2)=3,..
+C 11/09/26 KC via Claude
+C code changed to set IEND first, so that the test on JXBEG does not
+C read a value from an earlier call.
       ID(1) = 1
       IDM = 2
+      IEND = 0
       IF(LSHELL) THEN
 C  PRINT OUT EVERY ZONE IN H-BURNING SHELL
        IF(NPRTPT.LT.JXBEG) THEN

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 34 runs: 26 of YREC's cases, from a 0.04 $M_\odot$ brown dwarf and the solar models to 9 $M_\odot$, rotation, and helium burning, and 8 made-up ones. It cannot be reached at 7be4133: PINDEX's only caller, PUTSTORE, passes its own LSHELL, which nothing sets (putstore.f, suspected bug 3); in the runs it was false at every call. With PUTSTORE's bug fixed, the branch would run at each old-format store after the core's hydrogen is exhausted, but these runs hold no counts from PINDEX with that fix. The fix was not run. Intent: none found.

2. (cannot be reached) `M` listed twice.

(lines 65–72 ↗; line 83 ↗) In the LSHELL branch, if M is a multiple of NPRTPT (always, with the default NPRTPT = 1), the loop adds M and so does the final ID(IDM) = M; the general branch guards against this (line 77 ↗), the shell branch does not. PUTSTORE would then write the surface shell twice. Only if the LSHELL branch is taken (above).

A possible fix:

pindex.f

@@ lines 64–70 @@
    20    CONTINUE
        IF(NPRTPT.LT.M) THEN
           IBEG = INT(JXEND/NPRTPT+1)*NPRTPT
+C 11/09/26 KC via Claude
+C code changed to leave M to the final ID(IDM) = M, as the general
+C case does.
           IEND = INT(M/NPRTPT)*NPRTPT
+          IF (IEND .EQ. M) IEND = IEND - NPRTPT
           DO 30 J = IBEG,IEND,NPRTPT
              ID(IDM) = J
              IDM = IDM + 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 the same 34 runs. It cannot be reached at 7be4133: PINDEX's only caller, PUTSTORE, passes its own LSHELL, which nothing sets (putstore.f, suspected bug 3); in the runs it was false at every call. With PUTSTORE's bug fixed, the branch would run at each old-format store after the core's hydrogen is exhausted, but these runs hold no counts from PINDEX with that fix. With NPRTPT = 1, the default and the value in YREC's cases, it would then hold at every such store. The run of PUTSTORE's fix agrees: In testsuite/Test_m0p95feh+0p0_GN93_TAHB and testsuite/Test_m2p0_feh+0p0_vzams150_IES2_TAMS the .store file changed beyond 2% while the track stayed identical, as M written twice would make it (an inference from the code; the stored models were not compared line by line). The fix was not run. Intent: none found; the general branch guards against it (line 77 ↗).

Definitions

subroutinePINDEX

lines 6–85 · GitHub

Uses modules DIMENSIONS_MOD, NMLPAR_MOD

Called by PUTSTORE

Change log

The changes since rp_1

Changed since rp_1. Declarations only: COMMON blocks to modules (/CCOUT1/, now gone, to NMLPAR_MOD; the unused /CCOUT/ dropped), JSON from DIMENSIONS_MOD, IMPLICIT NONE.

1C 2C 3C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 4C PINDEX 5C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 6 SUBROUTINE PINDEX(JXBEG,JXEND,LSHELL,M,ID,IDM) 7C 10/15/26 MHP via Claude 8C code changed to take JSON from DIMENSIONS_MOD instead of declaring 9C PARAMETER(JSON=5000) here. 10 USE DIMENSIONS_MOD 11C 10/14/26 MHP via Claude 12C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD. 13 USE NMLPAR_MOD 14C 10/20/26 MHP via Claude 15C code changed to replace implicit typing with IMPLICIT NONE and 16C explicit declarations. 17 IMPLICIT NONE 18 INTEGER IBEG,ID,IDM,IEND,IXBEG,J,JXBEG,JXEND,M 19 LOGICAL*4 LSHELL 20 21C 10/11/26 MHP via Claude 22C code changed to remove LSCRIB from COMMON/CCOUT/: it was declared 23C positionally in all 17 files sharing this block but never read or 24C written in executable code anywhere, and no namelist file sets it. 25C 10/12/26 MHP via Claude 26C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 27C namelist variables that the physics routines read from 27 COMMON 28C blocks which also held derived or working members. Those blocks keep 29C only their non-namelist members; five of them -- NEWMX, VMULT, 30C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 31C Derivation-only namelist variables from the same blocks went to 32C COMMON/NMLBYPASS/ instead, not here. 33C 34C Declared identically in all 121 files that hold any of these members. 35C Generated from one member list, so every copy is byte-identical; the 36C linker's 'size of symbol changed' warning is the check that they are. 37C 10/13/26 MHP via Claude: removed the unused COMMON/CCOUT/ declaration 38C (no member referenced here; residue of the obsolete-option removals). 39 DIMENSION ID(JSON) 40 SAVE 41 42C ID IS THE ARRAY OF THE SHELLS TO BE PRINTED OUT; 43C I.E. IF SHELLS 1,3,5,..TO BE PRINTED ID(1)=1,ID(2)=3,.. 44 ID(1) = 1 45 IDM = 2 46 IF(LSHELL) THEN 47C PRINT OUT EVERY ZONE IN H-BURNING SHELL 48 IF(NPRTPT.LT.JXBEG) THEN 49 IBEG = MAX(2,NPRTPT) 50 IEND = INT(JXBEG/NPRTPT)*NPRTPT 51 DO 10 J = IBEG,IEND,NPRTPT 52 ID(IDM) = J 53 IDM = IDM + 1 54 10 10 CONTINUE 55 ENDIF 56 IF(IEND .EQ. JXBEG) THEN 57 IXBEG = JXBEG + 1 58 ELSE 59 IXBEG = JXBEG 60 END IF 61 DO 20 J = IXBEG,JXEND 62 ID(IDM) = J 63 IDM = IDM + 1 64 20 20 CONTINUE 65 IF(NPRTPT.LT.M) THEN 66 IBEG = INT(JXEND/NPRTPT+1)*NPRTPT 67 IEND = INT(M/NPRTPT)*NPRTPT 68 DO 30 J = IBEG,IEND,NPRTPT 69 ID(IDM) = J 70 IDM = IDM + 1 71 30 30 CONTINUE 72 ENDIF 73 ELSE IF(NPRTPT.LT.M) THEN 74C GENERAL CASE; PRINT OUT EVERY NPRTPT POINTS. 75 IBEG = MAX(2,NPRTPT) 76 IEND = INT(M/NPRTPT)*NPRTPT 77 IF (IEND .EQ. M) IEND = IEND - NPRTPT 78 DO 40 J = IBEG,IEND,NPRTPT 79 ID(IDM) = J 80 IDM = IDM + 1 81 40 40 CONTINUE 82 ENDIF 83 ID(IDM) = M 84 RETURN 85 END

1"""Python version of pindex.f (YREC 7be4133): which shells of a model are 2printed.""" 3 4import numpy as np 5 6from pyrec import nmlpar_mod 7 8# IEND is kept between calls (SAVE). Suspected bug (see the page): the 9# LSHELL branch can read it before setting it. 10iend = 0 11 12 13def pindex(jxbeg: int, jxend: int, lshell: bool, m: int, id_: np.ndarray, idm: int) -> int: 14 """List the shells to print in id_[1..idm]: 1, every NPRTPT-th shell 15 and m; with lshell, also every shell of the hydrogen-burning shell 16 jxbeg..jxend. 17 18 Translates: pindex.f#PINDEX (7be4133) 19 20 Arguments: jxbeg, jxend, the burning shell's edges; lshell; m, the 21 number of shells; id_, set in place (ID in the Fortran); idm (set). 22 Returns idm, the length of the list. 23 24 """ 25 26 global iend 27 nprtpt = nmlpar_mod.nprtpt 28 29 # The first shell (Fortran lines 44–45). 30 id_[1] = 1 31 idm = 2 32 33 # Every zone of the H-burning shell (Fortran lines 46–72). 34 if lshell: 35 36 # Every nprtpt-th shell below it (Fortran lines 48–55). 37 if nprtpt < jxbeg: 38 ibeg = max(2, nprtpt) 39 iend = (jxbeg // nprtpt) * nprtpt 40 for j in range(ibeg, iend + 1, nprtpt): 41 id_[idm] = j 42 idm = idm + 1 43 44 # The shell itself (Fortran lines 56–64). Suspected bug (see the 45 # page): IEND may be left from an earlier call. 46 if iend == jxbeg: 47 ixbeg = jxbeg + 1 48 else: 49 ixbeg = jxbeg 50 for j in range(ixbeg, jxend + 1): 51 id_[idm] = j 52 idm = idm + 1 53 54 # Every nprtpt-th shell above it (Fortran lines 65–72). 55 if nprtpt < m: 56 ibeg = (jxend // nprtpt + 1) * nprtpt 57 iend = (m // nprtpt) * nprtpt 58 # Suspected bug (see the page): m may be listed twice here. 59 for j in range(ibeg, iend + 1, nprtpt): 60 id_[idm] = j 61 idm = idm + 1 62 63 # The general case: every nprtpt-th shell (Fortran lines 73–82). 64 elif nprtpt < m: 65 ibeg = max(2, nprtpt) 66 iend = (m // nprtpt) * nprtpt 67 if iend == m: 68 iend = iend - nprtpt 69 for j in range(ibeg, iend + 1, nprtpt): 70 id_[idm] = j 71 idm = idm + 1 72 73 # The last shell (Fortran line 83). 74 id_[idm] = m 75 return idm

pindex.f

Purpose

How it fits

The method

Suspected bugs, checked by running the code

1. (cannot be reached) IEND undefined.

2. (cannot be reached) M listed twice.

Definitions

subroutinePINDEX

Change log

1. (cannot be reached) `IEND` undefined.

2. (cannot be reached) `M` listed twice.