getmodel2.f

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

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

Variables on this page (2)

clsun, json

Purpose

GETMODEL2 reads a stellar model stored in the "Model2" text format ("llp 4/16/03") into memory: the header (model number, mass, $T_\mathrm{eff}$, $L$, age), the physics flags of the run that wrote it, the luminosities by source, the envelope fitting data and the shell-by-shell structure and composition.

How it fits

Who calls it: only STARIN, when the first record of the input model starts with MOD2; for the older YREC7 format (NMOD) it calls GETYREC7 instead.
Who writes the format: PUTMODEL2 (and the format is mentioned in PUTSTORE and WRTLST).
Input: the file on unit IREAD, rewound first.
Outputs (arguments): the header values (MODEL, M, SMASS, TEFFL, BL, HSTOT, DAGE, DDAGE); the flags and names of the input physics (EOS, ATM, ALOK, HIK, COMPMIX, LDIFY, LDIFZ, LROT, …, and the disk-locking TDISK, PDISK, WMAX); the luminosities TLUMX (the three pp chains, CNO, helium burning, neutrinos, gravitational; converted to solar units if they exceed $10^{20}$, that is, if they were written in erg s$^{-1}$); the envelope triangle TRIT, TRIL, PS, TS, RS, CFENV and its orientation FTRI (±1: SURFBC flips its sign each time it reflects a corner of the triangle, and uses it in its test of whether ($T_\mathrm{eff}$, $L$) lies inside; PUTMODEL2 writes the record numbers as −1, −2, −3 when it is negative, and GETMODEL2 reads the sign back); and per shell HS, HR, HL, HP, HT, HD, OMEGA, LC and 15 abundances HCOMP.

Question 1 (lines 75–83 ↗). Units guessed from size. The luminosities are taken as erg s$^{-1}$ if any exceeds $10^{20}$, else as solar units. Could the format carry the unit instead?
Modules: CONSTANTS_MOD (formerly /CONST/), for $L_\odot$ (CLSUN).

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. (never occurred) The number of shells is not checked.

(lines 57–58 ↗, lines 103–106 ↗) The header record gives M and, in the same READ, stores HS(M); the loop over the shells then fills HS(1:M) and the other arrays of JSON = 5000 elements. Nothing checks that 1 ≤ M ≤ JSON, so a damaged file, or a model of more than 5000 shells, writes outside the arrays (the default build has no bounds checks, src/Makefile line 132; make devel adds them). A read error or an end of file, which has no END= or ERR=, stops the run with the compiler's run-time message, which does not name the file.

A possible fix:

getmodel2.f

@@ lines 24–31 @@
 C explicit declarations.
       IMPLICIT NONE
       REAL*8 BL,CCCMAX,CFENV,CMIXL,DAGE,DDAGE,FTRI,HCOMP,HD,HL,HP,HR,HS,
-     *       HSTOT,HT,OMEGA,PDISK,PS,RS,SMASS,TDISK,TEFFL,TLUMX,TRIL,
-     *       TRIT,TS,WMAX
+     *       HSM,HSTOT,HT,OMEGA,PDISK,PS,RS,SMASS,TDISK,TEFFL,TLUMX,
+     *       TRIL,TRIT,TS,WMAX
       INTEGER I,IO,IREAD,J,JCORE,JENV,K,M,MODEL
       LOGICAL*4 LC,LDIFY,LDIFZ,LDISK,LEXCOM,LINSTB,LJDOT0,LOVSTC,LOVSTE,
      *       LOVSTM,LPUREZ,LROT,LSEMIC
@@ lines 54–61 @@
       REWIND(IREAD)

 C get the header record
+C 11/09/26 KC via Claude
+C code changed to check the number of shells before it is used as an
+C index, as GETYREC7 checks the shells' masses.
       READ(IREAD,10) ATEMP,MODEL,M,SMASS,TEFFL,BL,HSTOT,DAGE,
-     *   DDAGE,HS(1),HS(M)
+     *   DDAGE,HS(1),HSM
+      IF(M.LT.1 .OR. M.GT.JSON) THEN
+         WRITE(*,*) 'GETMODEL2: NUMBER OF SHELLS ',M,
+     *      ' OUTSIDE 1 TO ',JSON,'. RUN STOPPED'
+         STOP 3
+      ENDIF
+      HS(M) = HSM
  10   FORMAT(A5,2I8,5F16.12,1PE18.10,0P2F16.12)

 C Get the flags describing the input physics and the information that

with HSM declared REAL*8.

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 header was read in all 33 runs (YREC's cases, and the made-up variants of them used for other pages), once each, three times in the solar calibration examples/run_standard_solar_model/Test_solarcal_noGS_norot, with M from 441 to 2,118 shells and model numbers from −3 to 852. M never left 1 to JSON, so no run wrote outside the arrays, and the fix was not run. What would reach it: a damaged model file, or one of more than 5,000 shells; nothing stops such a file before the arrays are overrun. Any input model in this format takes this path. No comment says the check was left out on purpose.

Definitions

subroutineGETMODEL2

lines 9–111 · GitHub

Reads a Model2-format model file: header, physics flags, luminosities, envelope data and the shells.

Uses modules CONSTANTS_MOD, DIMENSIONS_MOD

Called by STARIN

Change log

The changes since rp_1

Changed since rp_1. Declarations only: COMMON blocks to modules, IMPLICIT NONE.

1C 2C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 3C 4C GETMODEL2 - Read a Model2 format stellar model file into memory 5C 6C llp 4/16/03 7C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 8 9 SUBROUTINE GETMODEL2(BL,CFENV,CMIXL,DAGE,DDAGE,FTRI,HCOMP,HD,HL, 10C * HP,HR,HS,HSTOT,HT,IREAD,ISHORT,JCORE,JENV,LC,LEXCOM,LROT,M, ! KC 2025-05-31 11 * * HP,HR,HS,HSTOT,HT,IREAD,JCORE,JENV,LC,LEXCOM,LROT,M, 12 * * MODEL,OMEGA,PS,RS,SMASS,TEFFL,TLUMX,TRIL,TRIT,TS, 13 & & ATM,EOS,HIK,LDIFY,LDIFZ,LDISK,LINSTB,LJDOT0,ALOK, 14 & & LOVSTC,LOVSTE,LOVSTM,LPUREZ,LSEMIC,COMPMIX,PDISK,TDISK,WMAX) 15C 10/15/26 MHP via Claude 16C code changed to take JSON from DIMENSIONS_MOD instead of declaring 17C PARAMETER(JSON=5000) here. 18 USE DIMENSIONS_MOD 19C 10/14/26 MHP via Claude 20C code changed to convert COMMON/CONSTANTS/ into CONSTANTS_MOD. 21 USE CONSTANTS_MOD 22C 10/20/26 MHP via Claude 23C code changed to replace implicit typing with IMPLICIT NONE and 24C explicit declarations. 25 IMPLICIT NONE 26 REAL*8 BL,CCCMAX,CFENV,CMIXL,DAGE,DDAGE,FTRI,HCOMP,HD,HL,HP,HR,HS, 27 * * HSTOT,HT,OMEGA,PDISK,PS,RS,SMASS,TDISK,TEFFL,TLUMX,TRIL, 28 * * TRIT,TS,WMAX 29 INTEGER I,IO,IREAD,J,JCORE,JENV,K,M,MODEL 30 LOGICAL*4 LC,LDIFY,LDIFZ,LDISK,LEXCOM,LINSTB,LJDOT0,LOVSTC,LOVSTE, 31 * * LOVSTM,LPUREZ,LROT,LSEMIC 32C First three lines above are YREC7 inputs 33C Last two lines are MODEL2 add-ons 34C MHP 4/25 chanted LOK to ALOK to avoid variable name conflicts 35 36 37 CHARACTER*6 EOS 38 CHARACTER*5 ATEMP 39c CHARACTER*4 ATM, LOK, HIK, COMPMIX 40c MHP 4/25 changed LOK name to make it unique, used elsewhere 41 CHARACTER*4 ATM, ALOK, HIK, COMPMIX 42 43C 9/25/26 MHP via Claude 44C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/, 45C and COMMON/CONST3/ into a single COMMON/CONSTANTS/ 46 47 48 DIMENSION HCOMP(15,JSON),HD(JSON),HL(JSON),HP(JSON),HR(JSON), 49 & & HS(JSON),HT(JSON),LC(JSON),TRIT(3),TRIL(3),PS(3),TS(3),RS(3), 50 & & CFENV(9),TLUMX(8),OMEGA(JSON) 51 52 save 53 54 REWIND(IREAD) 55 56C get the header record 57 READ(IREAD,10) ATEMP,MODEL,M,SMASS,TEFFL,BL,HSTOT,DAGE, 58 * * DDAGE,HS(1),HS(M) 59 10 10 FORMAT(A5,2I8,5F16.12,1PE18.10,0P2F16.12) 60 61C Get the flags describing the input physics and the information that 62C is stored in the model (rotating or not?, extendec compotition, etc.) 63 READ(IREAD,30) JCORE,JENV,CMIXL,EOS,ATM,ALOK,HIK,LPUREZ, 64 & & COMPMIX,LEXCOM,LDIFY,LDIFZ,LSEMIC,LOVSTC,LOVSTE,LOVSTM, 65 & & LROT,LINSTB,LJDOT0,LDISK,TDISK,PDISK,WMAX 66 30 30 FORMAT(2I8,F16.10,1X,A6,1X,3(A4,1X),L1,1X,A4,1X,11(L1,1X), 67 & & 3(1PE18.10)) 68 69C Get the LUMINOSITIES---PP(1-2-3)-CNO-HE-NU-GRAV 70C CONVERT TO SOLAR LUMINOSITIES 71 READ(IREAD,40) (TLUMX(J),J=1,7) 72c 40 FORMAT('TLUMX',5X,1P7E17.9) 73 40 40 FORMAT(10X,1P7E17.9) 74 75C If TLUMX are in ergs, convert to solar units. Decide by 76C comparing to 10**20. IF larger, divide by CLSUN. 77 CCCMAX = DMAX1(TLUMX(1),TLUMX(2),TLUMX(3),TLUMX(4),TLUMX(5), 78 * * DABS(TLUMX(6)),TLUMX(7)) 79 IF(CCCMAX.GT.1.0D20) THEN 80 DO J = 1,7 81 TLUMX(J) = TLUMX(J)/CLSUN 82 ENDDO 83 ENDIF 84 85C Get the ENVELOPE DATA 86C FTRI is 1,normally. It is set to -1 if any of the record numbers 87C for the envelope triangle records was set to -1 by WRTLST. 88 FTRI = 1D0 89 DO I = 1,3 90 READ(IREAD,50)IO,TRIT(I),TRIL(I),PS(I),TS(I),RS(I), 91 & & (CFENV(3*I-3+J),J=1,3) 92 IF(IO .LT. 0) FTRI = -1D0 93 ENDDO 94c 50 FORMAT('ENV',I2,5F16.12,1P3E20.12) 95 50 50 FORMAT(4X,I2,5F16.12,1P3E20.12) 96 97C Skip over column headings for per-shell data 98 READ(IREAD,*) 99 READ(IREAD,*) 100 READ(IREAD,*) 101 102C Read per-shell data (HENYEY POINGS) 103 DO I = 1,M 104 READ(IREAD,60) K,HS(I),HR(I),HL(I),HP(I),HT(I),HD(I), 105 * * OMEGA(I),LC(I),(HCOMP(J,I),J=1,15) 106 ENDDO 107 60 60 FORMAT(I6,0P2F18.14,1PE24.16,0P3F18.14,1PE24.16,1x,L1, 108 & & 3(0PF12.9),12(0PE16.8)) 109 110 RETURN 111 END

1"""Python version of getmodel2.f (YREC 7be4133): read a stellar model in 2the Model2 text format.""" 3 4from pyrec import constants_mod 5import numpy as np 6from pyrec import iounits_mod 7 8 9def _fields(line, spec): 10 """Read one record by a FORMAT, given as (kind, width[, decimals]) 11 items: 'a' text, 'i' integer, 'f' real (F or E), 'l' logical, 'x' 12 skipped. Blanks inside numbers are ignored and a blank field reads as 13 0, as Fortran does by default; a real without a decimal point has its 14 last `decimals` digits taken as decimals. (The 1P scale factor of the 15 FORMATs only matters for fields without an exponent, which PUTMODEL2 16 does not write, and is not applied.)""" 17 18 out, pos = [], 0 19 for item in spec: 20 kind, w = item[0], item[1] 21 text = line[pos:pos + w] if line is not None else '' 22 pos += w 23 if kind == 'x': 24 continue 25 if kind == 'a': 26 out.append(text.ljust(w)) 27 continue 28 t = text.replace(' ', '') 29 if kind == 'l': 30 t = t.lstrip('.') 31 out.append(t[:1].upper() == 'T') 32 elif kind == 'i': 33 out.append(int(t) if t else 0) 34 else: 35 t = t.upper().replace('D', 'E') 36 if not t: 37 out.append(0.0) 38 continue 39 mant, sep, expo = t.partition('E') 40 value = float(mant) 41 if '.' not in mant: 42 value = value / 10.0 ** item[2] 43 if sep: 44 value = value * 10.0 ** int(expo) 45 out.append(value) 46 return out 47 48 49def getmodel2(bl: float, cfenv: np.ndarray, cmixl: float, dage: float, ddage: float, ftri: float, hcomp: np.ndarray, 50 hd: np.ndarray, hl: np.ndarray, hp: np.ndarray, hr: np.ndarray, hs: np.ndarray, hstot: float, 51 ht: np.ndarray, iread: int, jcore: int, jenv: int, lc: np.ndarray, lexcom: bool, lrot: bool, m: int, 52 model: int, omega: np.ndarray, ps: np.ndarray, rs: np.ndarray, smass: float, teffl: float, 53 tlumx: np.ndarray, tril: np.ndarray, trit: np.ndarray, ts: np.ndarray, atm: str, eos: str, hik: str, 54 ldify: bool, ldifz: bool, ldisk: bool, linstb: bool, ljdot0: bool, alok: str, lovstc: bool, lovste: bool, 55 lovstm: bool, lpurez: bool, lsemic: bool, compmix: str, pdisk: float, tdisk: float, 56 wmax: float) -> tuple[float, float, float, float, float, float, int, int, bool, bool, int, int, float, 57 float, str, str, str, bool, bool, bool, bool, bool, str, bool, bool, bool, bool, bool, str, float, float, 58 float]: 59 """Read a Model2-format model from unit iread (rewound first): the 60 header, the physics flags of the run that wrote it, the luminosities 61 by source (converted to solar units if written in erg/s), the 62 envelope triangle with its orientation ftri, and the shells. 63 64 Translates: getmodel2.f#GETMODEL2 (7be4133) 65 66 Arguments: iread (the unit); every other scalar is only set. 67 Arrays changed in place: cfenv, hcomp, hd, hl, hp, hr, hs, ht, lc, 68 omega, ps, rs, tlumx, tril, trit, ts. 69 70 Returns (bl, cmixl, dage, ddage, ftri, hstot, jcore, jenv, lexcom, 71 lrot, m, model, smass, teffl, atm, eos, hik, ldify, ldifz, ldisk, 72 linstb, ljdot0, alok, lovstc, lovste, lovstm, lpurez, lsemic, 73 compmix, pdisk, tdisk, wmax). 74 75 """ 76 77 # Fortran line 54. 78 iounits_mod.rewind(iread) 79 80 # Header (FORMAT 10): A5,2I8,5F16.12,1PE18.10,0P2F16.12 (Fortran lines 81 # 56–59). 82 # Suspected bug (see the page): m is not checked against 1..json 83 # before hs[m] is set. 84 (atemp, model, m, smass, teffl, bl, hstot, dage, ddage, hs1, 85 hsm) = _fields(iounits_mod.read(iread), 86 [('a', 5), ('i', 8), ('i', 8)] + [('f', 16, 12)] * 5 87 + [('f', 18, 10), ('f', 16, 12), ('f', 16, 12)]) 88 hs[1] = hs1 89 hs[m] = hsm 90 91 # The run's physics (FORMAT 30) (Fortran lines 61–67). 92 vals = _fields(iounits_mod.read(iread), 93 [('i', 8), ('i', 8), ('f', 16, 10), ('x', 1), ('a', 6), 94 ('x', 1)] + [('a', 4), ('x', 1)] * 3 95 + [('l', 1), ('x', 1), ('a', 4), ('x', 1)] 96 + [('l', 1), ('x', 1)] * 11 + [('f', 18, 10)] * 3) 97 (jcore, jenv, cmixl, eos, atm, alok, hik, lpurez, compmix, lexcom, 98 ldify, ldifz, lsemic, lovstc, lovste, lovstm, lrot, linstb, ljdot0, 99 ldisk, tdisk, pdisk, wmax) = vals 100 101 # Luminosities by source (FORMAT 40); in erg/s if any exceeds 1e20 102 # (Fortran lines 69–83). 103 tlumx[1:8] = _fields(iounits_mod.read(iread), 104 [('x', 10)] + [('f', 17, 9)] * 7) 105 # Question (see the page): the unit is guessed from the size. 106 cccmax = max(tlumx[1], tlumx[2], tlumx[3], tlumx[4], tlumx[5], 107 abs(tlumx[6]), tlumx[7]) 108 if cccmax > 1.0e20: 109 for j in range(1, 7 + 1): 110 tlumx[j] = tlumx[j] / constants_mod.clsun 111 112 # The envelope triangle (FORMAT 50); a negative record number means 113 # its orientation is reversed (Fortran lines 85–95). 114 ftri = 1.0 115 for i in range(1, 3 + 1): 116 vals = _fields(iounits_mod.read(iread), 117 [('x', 4), ('i', 2)] + [('f', 16, 12)] * 5 118 + [('f', 20, 12)] * 3) 119 io, trit[i], tril[i], ps[i], ts[i], rs[i] = vals[:6] 120 for j in range(1, 3 + 1): 121 cfenv[3 * i - 3 + j] = vals[5 + j] 122 if io < 0: 123 ftri = -1.0 124 125 # Three header records of the shell table (Fortran lines 97–100). 126 iounits_mod.read(iread) 127 iounits_mod.read(iread) 128 iounits_mod.read(iread) 129 130 # The shells (FORMAT 60) (Fortran lines 102–108). 131 spec60 = ([('i', 6), ('f', 18, 14), ('f', 18, 14), ('f', 24, 16)] 132 + [('f', 18, 14)] * 3 + [('f', 24, 16), ('x', 1), ('l', 1)] 133 + [('f', 12, 9)] * 3 + [('f', 16, 8)] * 12) 134 for i in range(1, m + 1): 135 vals = _fields(iounits_mod.read(iread), spec60) 136 (k, hs[i], hr[i], hl[i], hp[i], ht[i], hd[i], omega[i], 137 lc[i]) = vals[:9] 138 for j in range(1, 15 + 1): 139 hcomp[j, i] = vals[8 + j] 140 141 return (bl, cmixl, dage, ddage, ftri, hstot, jcore, jenv, lexcom, lrot, 142 m, model, smass, teffl, atm, eos, hik, ldify, ldifz, ldisk, 143 linstb, ljdot0, alok, lovstc, lovste, lovstm, lpurez, lsemic, 144 compmix, pdisk, tdisk, wmax)