getyrec7.f

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

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

Variables on this page (2)

clsun, json

Purpose

GETYREC7 reads a stellar model stored in the older YREC7 text format ("llp 4/16/03"; its header calls it GETY7) into memory. Since the YREC7 format lacks the physics flags of the newer Model2 format, it sets those to defaults ("?" strings, false flags).

How it fits

Who calls it: only STARIN, when the input model's first record starts with NMOD; Model2 files (MOD2) go to GETMODEL2, which returns the same arguments.
Who writes the format: PUTYREC7 (its formats match these), which nothing calls; so YREC7 files now come only from older codes.
Input: the file on unit IREAD, rewound before and after; ISHORT, the unit of its error messages (lines 140–141 ↗, line 249 ↗), so that, unlike in GETMODEL2, this argument is used.
Outputs (arguments): the header (MODEL, M, SMASS, TEFFL, BL, HSTOT, DAGE, DDAGE), LROT, LEXCOM, JCORE, JENV, CMIXL, the luminosities TLUMX (converted to solar units if above $10^{20}$), the envelope triangle (TRIT, TRIL, PS, TS, RS, CFENV, FTRI), the shells (HS, HR, HL, HP, HT, HD, LC), HCOMP and OMEGA; plus the Model2-only items, set to defaults: the strings ATM, EOS, HIK, ALOK = "?," COMPMIX = "AG93," all flags false, PDISK = 7.2722×10$^{-6}$ rad s$^{-1}$ (a 10-day period), TDISK = 0, WMAX = 3×10$^{-4}$ (lines 86–89 ↗). These three do not change the run's physics: STARIN receives them in its locals PDISK0, TDISK0, WMAX0 (starin.f, lines 280–285 ↗), which it never reads (unused?), and the namelist's PDISK, TDISK, WMAX are untouched (settled from the code, Nov 9, 2026).
Modules: CONSTANTS_MOD (formerly /CONST/), for $L_\odot$.

The format, as read.

Each shell stores $X$ and $Z$ as integers in units of $10^{-6}$; $Y = 1 - X - Z - X(^3\mathrm{He})$ is recomputed. The run stops (exit status 3, an input error) if a shell's $\log m$ is negative or above the total.

Question 1 (lines 148–149 ↗). Precision of $X$ and $Z$. The format keeps $X$ and $Z$ to $10^{-6}$. Is that enough when models with diffusion are restarted from this format?
The other species (rows 4–11) are stored only from JCORE to JENV, the region between the convective core and envelope, and copied into the homogeneous core and envelope. With LEXCOM, the light elements (rows 12–15) follow the same pattern, two shells per record.
With LROT, $\Omega$ is stored as $-\log_{10}\Omega$; values of 58.9 or more mean $\Omega = 0$. If those records are missing, $\Omega$ is set to zero with a message.

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. (not reached) The flags record is read one column late.

(lines 104–105 ↗; putyrec7.f, lines 65–66 ↗) PUTYREC7 writes the record as 2(L1,1X),'F F ',2I4,F7.3, four characters between the flags and JCORE, and the YREC7 starting models in startmodels/seed/ have the same layout (F T F F 1 1 1.735); GETYREC7 reads it as 2(L1,1X),5X,2I4,F7.3, skipping five. Each field is then read one column late: JCORE from columns 10–13, which hold its last three digits and the first column of JENV; JENV from columns 14–17, its last three digits and the first column of CMIXL. With blanks ignored, the values come out right while JCORE and JENV are below 1000 (all 41 YREC7 files in the repository have at most 391 shells); with a four-digit JENV (M may reach JSON = 5000) its first digit is appended to JCORE and lost from JENV.

A possible fix:

getyrec7.f

@@ lines 102–108 @@
 C Get the flags describing the input physics and the information that
 C is stored in the model (rotating or not?, extendec compotition, etc.)
       READ(IREAD,30) LROT,LEXCOM,JCORE,JENV,CMIXL
- 30   FORMAT(2(L1,1X),5X,2I4,F7.3)
+C 11/09/26 KC via Claude
+C code changed to skip the four columns PUTYREC7 writes ('F F ')
+C between the flags and JCORE, not five.
+ 30   FORMAT(2(L1,1X),4X,2I4,F7.3)

 C get the LUMINOSITIES---PP(1-2-3)-CNO-HE-NU-GRAV
       READ(IREAD,40) (TLUMX(J),J=1,7)

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: None of the 33 runs started from a YREC7 file (the NMOD seeds in startmodels/seed/, which the cases of examples/run_from_seed_to_start/ and examples/run_rotation_to_RG/ read). The fix was not run. With those seeds the values would come out right, as above; what would reach the misreading is a YREC7 file with a four-digit JENV. No comment says the fifth column is meant; PUTYREC7 writes four.

Definitions

subroutineGETYREC7

lines 9–262 · GitHub

Reads a YREC7-format model file (header, envelope data, shells, composition in compressed form, rotation) and sets the Model2-only items to defaults.

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.

Changed since a8a5950. The stop on a bad shell mass now exits with status 3, the code for input errors (KNOWN_ISSUES.md, issue 38).

1C 2C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 3C 4C GETY7 - Read a YREC7 format stellar model file into memory 5C 6C llp 4/16/03 7C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 8 9 SUBROUTINE GETYREC7(BL,CFENV,CMIXL,DAGE,DDAGE,FTRI,HCOMP,HD,HL, 10 * * HP,HR,HS,HSTOT,HT,IREAD,ISHORT,JCORE,JENV,LC,LEXCOM,LROT,M, 11 * * MODEL,OMEGA,PS,RS,SMASS,TEFFL,TLUMX,TRIL,TRIT,TS, 12 & & ATM,EOS,HIK,LDIFY,LDIFZ,LDISK,LINSTB,LJDOT0,ALOK, 13 & & LOVSTC,LOVSTE,LOVSTM,LPUREZ,LSEMIC,COMPMIX,PDISK,TDISK,WMAX) 14C 10/15/26 MHP via Claude 15C code changed to take JSON from DIMENSIONS_MOD instead of declaring 16C PARAMETER(JSON=5000) here. 17 USE DIMENSIONS_MOD 18C 10/14/26 MHP via Claude 19C code changed to convert COMMON/CONSTANTS/ into CONSTANTS_MOD. 20 USE CONSTANTS_MOD 21C 10/20/26 MHP via Claude 22C code changed to replace implicit typing with IMPLICIT NONE and 23C explicit declarations. 24 IMPLICIT NONE 25 REAL*8 BL,CCCMAX,CFENV,CMIXL,DAGE,DDAGE,FTRI,HCOMP,HD,HL,HP,HR,HS, 26 * * HSTOT,HT,OMEGA,OMEGAL,PDISK,PS,RS,SMASS,TDISK,TEFFL,TLUMX, 27 * * TRIL,TRIT,TS,WMAX 28 INTEGER I,II,IO,IREAD,ISHORT,IX,IZ,J,JCORE,JENV,K,KT,M,MODEL 29 LOGICAL*4 LC,LDIFY,LDIFZ,LDISK,LEXCOM,LINSTB,LJDOT0,LOVSTC,LOVSTE, 30 * * LOVSTM,LPUREZ,LROT,LSEMIC 31C First three lines above are YREC7 inputs 32C Last two lines are MODEL2 add-ons 33 34 35c CHARACTER*4 ATM, LOK, HIK, COMPMIX 36c MHP 4/25 changed LOK name to make it unique, used elsewhere 37 CHARACTER*4 ATEMP,ATM,HIK,ALOK,COMPMIX 38 CHARACTER*6 EOS 39 40C 9/25/26 MHP via Claude 41C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/, 42C and COMMON/CONST3/ into a single COMMON/CONSTANTS/ 43 44 DIMENSION CFENV(9),HCOMP(15,JSON),HD(JSON),HL(JSON),HR(JSON), 45 * *HP(JSON),HS(JSON),HT(JSON),LC(JSON),OMEGA(JSON),OMEGAL(JSON), 46 * *PS(3),RS(3),TS(3),TLUMX(8),TRIL(3),TRIT(3) 47 48 save 49 50C Initialization 51 52C Zero the composition arrays (HCOMP) and the OMEGA array. 53C This is agreed upon protection in case they are inadvertently used in the program, 54C but not present in the old YREC input data. 55 DO I = 1, JSON 56 OMEGA(I) = 0D0 57 DO J = 1, 15 58 HCOMP(J,I) = 0D0 59 ENDDO 60 ENDDO 61 62C Set Model2-sspecific strings to "?" 63c ATM = " ? " 64c EOS = " ? " 65c HIK = " ? " 66c LOK = " ? " 67c MHP 4/25 padded character strings to proper length 68 ATM = " ? " 69 EOS = " ? " 70 HIK = " ? " 71 ALOK = " ? " 72 COMPMIX = "AG93" 73 74C Set Model2-specific flags false 75 LDIFY = .FALSE. 76 LDIFZ = .FALSE. 77 LDISK = .FALSE. 78 LINSTB = .FALSE. 79 LJDOT0 = .FALSE. 80 LOVSTC = .FALSE. 81 LOVSTE = .FALSE. 82 LOVSTM = .FALSE. 83 LPUREZ = .FALSE. 84 LSEMIC = .FALSE. 85 86C Set Model2-specific constants to default 87 PDISK = 7.2722D-6 88 TDISK = 0D0 89 WMAX = 3D-4 90 91C End of initialization 92 93C Read in an input model from the LU IREAD file in YREC7 Format 94 95 REWIND(IREAD) 96 97C get the header record 98 READ(IREAD,10) ATEMP,MODEL,M,SMASS,TEFFL,BL,HSTOT, 99 * * DAGE,DDAGE 100 10 10 FORMAT(A4,1X,2I5,3F8.5,1X,F14.10,F13.9,F13.0) 101 102C Get the flags describing the input physics and the information that 103C is stored in the model (rotating or not?, extendec compotition, etc.) 104 READ(IREAD,30) LROT,LEXCOM,JCORE,JENV,CMIXL 105 30 30 FORMAT(2(L1,1X),5X,2I4,F7.3) 106 107C get the LUMINOSITIES---PP(1-2-3)-CNO-HE-NU-GRAV 108 READ(IREAD,40) (TLUMX(J),J=1,7) 109 40 40 FORMAT(10X,7E10.3) 110 111C If TLUMX are in ergs, convert to solar units. Decide by 112C comparing to 10**20. If larger, divide by CLSUN. 113 CCCMAX = DMAX1(TLUMX(1),TLUMX(2),TLUMX(3),TLUMX(4),TLUMX(5), 114 * * DABS(TLUMX(6)),TLUMX(7)) 115 IF(CCCMAX.GT.1.0D20) THEN 116 DO J = 1,7 117 TLUMX(J) = TLUMX(J)/CLSUN 118 ENDDO 119 ENDIF 120 121C Get the ENVELOPE DATA 122C FTRI is 1,normally. It is set to -1 if any of the record numbers 123C for the envelope triangle records was set to -1 by WRTLST. 124 FTRI = 1D0 125 DO 80 I = 1,3 126 READ(IREAD,70) IO,TRIT(I),TRIL(I),PS(I),TS(I),RS(I), 127 * * (CFENV(I+I+I-3+J),J = 1,3) 128 70 70 FORMAT(3X,I2,F7.5,4F8.5,3E12.5) 129 IF(IO.LT.0) FTRI = -1D0 130 80 80 CONTINUE 131 132C READ IN HENYEY POINTS IN 4 PARTS 133C 134C Read FIRST PART:M,R,L,P,T,RHO,CONV(T/F),X,Y,AND Z - one line per shell 135 DO 110 I = 1,M 136 READ(IREAD,100) HS(I),HR(I),HL(I),HP(I),HT(I),HD(I),LC(I), 137 * * IX,IZ 138 100 100 FORMAT(0PF13.10,F10.7,1PE14.7,0PF10.7,2F10.7,L1,2I6) 139 IF(HS(I).LT.0D0.OR.HS(I).GT.HSTOT) THEN 140 WRITE(ISHORT,1000) 141 WRITE(ISHORT,1050) I 142 1000 1000 FORMAT(1X,39('>'),40('<')/1X,'RUN STOPPED') 143 1050 1050 FORMAT(' ERROR IN SUBROUTINE GETY7'/1X,'GLITCH IN SHELL', 144 * * I3,', SHELL MASS LESS THAN ZERO OR GREATER THAN STAR MASS') 145C 10/30/26 MHP via Claude: bare STOP -> STOP 3, the input-error code (issue 38). 146 STOP 3 147 ENDIF 148 HCOMP(1,I) = 1.0D-6*DFLOAT(IX) 149 HCOMP(3,I) = 1.0D-6*DFLOAT(IZ) 150 110 110 CONTINUE 151 152C Read SECOND PART:ELEMENT ABUNDANCES: HE3,CNO CYCLE ELEMENTS. 153C ABUNDANCES IN SURFACE AND CENTRAL CONVECTION ZONES STORED WITH 1 154C VALUE PER ZONE RATHER THAN 1 VALUE PER SHELL 155 IF(JENV.LE.JCORE) THEN 156 READ(IREAD,200)(HCOMP(I,JCORE),I=4,11) 157 READ(IREAD,200)(HCOMP(I,JENV),I=4,11) 158 ELSE 159 DO 210 J = JCORE,JENV 160 READ(IREAD,200) (HCOMP(I,J),I = 4,11) 161 200 200 FORMAT(8(1PE9.3,1X)) 162 210 210 CONTINUE 163 ENDIF 164 IF(JCORE.GT.1) THEN 165C CONVECTIVE CORE- ASSIGN FIRST COMPOSITION VALUE TO SHELLS 1-JCORE 166 DO 230 J = 1,JCORE-1 167 DO 220 I = 4,11 168 HCOMP(I,J) = HCOMP(I,JCORE) 169 220 220 CONTINUE 170 230 230 CONTINUE 171 ENDIF 172 IF(JENV.LT.M) THEN 173C CONVECTIVE SURFACE- ASSIGN LAST COMPOSITION TO SHELLS JENV-M 174 DO 250 J = JENV+1,M 175 DO 240 I = 4,11 176 HCOMP(I,J) = HCOMP(I,JENV) 177 240 240 CONTINUE 178 250 250 CONTINUE 179 ENDIF 180C DEFINE HE4 = 1 - X - Z - HE3. 181 DO 260 I = 1,M 182 HCOMP(2,I) = 1.0D0 - HCOMP(1,I) - HCOMP(3,I) - HCOMP(4,I) 183 260 260 CONTINUE 184 185C READ IN H2,LI6,LI7,BE9 (EXTENDED COMPOSITION) 186 187 IF(LEXCOM) THEN 188 IF(JENV.EQ.1)THEN 189C FULLY CONVECTIVE MODEL 190 READ(IREAD,300)(HCOMP(I,1),I=12,15) 191 300 300 FORMAT(4(1PE9.3,1X)) 192 DO 310 J = 1,M 193 DO 305 I = 12,15 194 HCOMP(I,J) = HCOMP(I,1) 195 305 305 CONTINUE 196 310 310 CONTINUE 197 ELSE 198C GENERAL CASE 199C THESE ABUNDANCES ARE READ IN WITH 2 SHELLS PER LINE 200 DO 315 K = JCORE,JENV-1,2 201 READ(IREAD,200)((HCOMP(I,J),I = 12,15),J = K,K+1) 202 315 315 CONTINUE 203C IF AN ODD NUMBER OF ABUNDANCES STORED, READ IN LAST VALUE 204 KT = JENV-1 - JCORE 205 IF(MOD(KT,2).NE.0)READ(IREAD,300)(HCOMP(I,K),I = 12,15) 206 IF(JCORE.GT.1) THEN 207C CONVECTIVE CORE- ASSIGN FIRST VALUE TO SHELLS 1-JCORE 208 DO 330 J = 1,JCORE-1 209 DO 320 I = 12,15 210 HCOMP(I,J) = HCOMP(I,JCORE) 211 320 320 CONTINUE 212 330 330 CONTINUE 213 ENDIF 214 IF(JENV.LT.M) THEN 215C CONVECTIVE SURFACE- ASSIGN LAST VALUE TO SHELLS JENV - M 216 DO 350 J = JENV+1,M 217 DO 340 I = 12,15 218 HCOMP(I,J) = HCOMP(I,JENV) 219 340 340 CONTINUE 220 350 350 CONTINUE 221 ENDIF 222 ENDIF 223 ENDIF 224 225C Read in FOURTH PART: - LOG J/M STORED , 8 SHELLS PER LINE 226 227 IF(LROT) THEN 228C READ OMEGA IN. If OMEGA records are missing go to OMEGA BYPASS below 229 READ(IREAD,500,END=9999)(OMEGAL(II),II = 1,M) 230 500 500 FORMAT(0P8F10.7) 231 DO 510 I = 1,M 232 IF(OMEGAL(I) .LT. 58.9D0) THEN 233 OMEGA(I) = 10D0**(-OMEGAL(I)) 234 ELSE 235 OMEGA(I) = 0D0 236 ENDIF 237 510 510 CONTINUE 238 ENDIF 239 240C KEEP IREAD OPEN 241 REWIND IREAD 242 243 RETURN 244 245C OMEGA BYPASS 246C Come here if OMEGA records are missing -- Can happen if 247C a new omega file is being generated (LWNEW is true in STARIN) 248 2499999 9999 WRITE(ISHORT,9998) "GETYREC7: OMEGA records are missing from ", 250 * * " input model file - OMEGA array zeroed" 2519998 9998 FORMAT(2A) 252 253 DO I = 1, JSON 254 OMEGA(I) = 0D0 255 ENDDO 256 257C KEEP IREAD OPEN 258 REWIND IREAD 259 260 RETURN 261 262 END

1"""Python version of getyrec7.f (YREC 7be4133): read a stellar model in 2the older YREC7 text format.""" 3 4from pyrec import constants_mod 5import numpy as np 6from pyrec import dimensions_mod 7from pyrec import iounits_mod 8 9 10def _fields(line, spec): 11 """Read one record by a FORMAT, given as (kind, width[, decimals]) 12 items: 'a' text, 'i' integer, 'f' real (F or E), 'l' logical, 'x' 13 skipped. Blanks inside numbers are ignored and a blank field reads as 14 0, as Fortran does by default; a real without a decimal point has its 15 last `decimals` digits taken as decimals. (The 1P scale factor of the 16 FORMATs only matters for fields without an exponent, and is not 17 applied.)""" 18 19 out, pos = [], 0 20 for item in spec: 21 kind, w = item[0], item[1] 22 text = line[pos:pos + w] 23 pos += w 24 if kind == 'x': 25 continue 26 if kind == 'a': 27 out.append(text.ljust(w)) 28 continue 29 t = text.replace(' ', '') 30 if kind == 'l': 31 t = t.lstrip('.') 32 out.append(t[:1].upper() == 'T') 33 elif kind == 'i': 34 out.append(int(t) if t else 0) 35 else: 36 t = t.upper().replace('D', 'E') 37 if not t: 38 out.append(0.0) 39 continue 40 mant, sep, expo = t.partition('E') 41 value = float(mant) 42 if '.' not in mant: 43 value = value / 10.0 ** item[2] 44 if sep: 45 value = value * 10.0 ** int(expo) 46 out.append(value) 47 return out 48 49 50def _read_values(unit, n, spec): 51 """Read n values with a FORMAT that is reused for further records 52 (Fortran's format reversion); None at the end of the file.""" 53 54 values = [] 55 while len(values) < n: 56 line = iounits_mod.read(unit) 57 if line is None: 58 return None 59 values += _fields(line, spec) 60 return values[:n] 61 62 63def getyrec7(bl: float, cfenv: np.ndarray, cmixl: float, dage: float, ddage: float, ftri: float, hcomp: np.ndarray, 64 hd: np.ndarray, hl: np.ndarray, hp: np.ndarray, hr: np.ndarray, hs: np.ndarray, hstot: float, 65 ht: np.ndarray, iread: int, ishort: int, jcore: int, jenv: int, lc: np.ndarray, lexcom: bool, lrot: bool, 66 m: int, model: int, omega: np.ndarray, ps: np.ndarray, rs: np.ndarray, smass: float, teffl: float, 67 tlumx: np.ndarray, tril: np.ndarray, trit: np.ndarray, ts: np.ndarray, atm: str, eos: str, hik: str, 68 ldify: bool, ldifz: bool, ldisk: bool, linstb: bool, ljdot0: bool, alok: str, lovstc: bool, lovste: bool, 69 lovstm: bool, lpurez: bool, lsemic: bool, compmix: str, pdisk: float, tdisk: float, 70 wmax: float) -> tuple[float, float, float, float, float, float, int, int, bool, bool, int, int, float, 71 float, str, str, str, bool, bool, bool, bool, bool, str, bool, bool, bool, bool, bool, str, float, float, 72 float]: 73 """Read a YREC7-format model from unit iread (rewound before and 74 after): the header, the luminosities, the envelope triangle, the 75 shells with X and Z in units of 1e-6, the other species stored only 76 between the convective core and envelope, and Omega (as -log Omega) 77 with rotation. The Model2-only items get defaults. 78 79 Translates: getyrec7.f#GETYREC7 (7be4133) 80 81 Arguments: iread (the unit), ishort (for the error messages); every 82 other scalar is only set. 83 Arrays changed in place: cfenv, hcomp, hd, hl, hp, hr, hs, ht, lc, 84 omega, ps, rs, tlumx, tril, trit, ts. 85 86 Returns (bl, cmixl, dage, ddage, ftri, hstot, jcore, jenv, lexcom, 87 lrot, m, model, smass, teffl, atm, eos, hik, ldify, ldifz, ldisk, 88 linstb, ljdot0, alok, lovstc, lovste, lovstm, lpurez, lsemic, 89 compmix, pdisk, tdisk, wmax). 90 91 """ 92 93 # Zero the composition and Omega, in case they are used though the 94 # file lacks them (Fortran lines 50–60). 95 json = dimensions_mod.json 96 for i in range(1, json + 1): 97 omega[i] = 0.0 98 for j in range(1, 15 + 1): 99 hcomp[j, i] = 0.0 100 101 # The Model2-only items, which this format lacks (Fortran lines 62–89). 102 atm = " ? " 103 eos = " ? " 104 hik = " ? " 105 alok = " ? " 106 compmix = "AG93" 107 ldify = False 108 ldifz = False 109 ldisk = False 110 linstb = False 111 ljdot0 = False 112 lovstc = False 113 lovste = False 114 lovstm = False 115 lpurez = False 116 lsemic = False 117 pdisk = 7.2722e-6 118 tdisk = 0.0 119 wmax = 3e-4 120 121 # Fortran line 95. 122 iounits_mod.rewind(iread) 123 124 # Header (FORMAT 10) (Fortran lines 97–100). 125 (atemp, model, m, smass, teffl, bl, hstot, dage, ddage) = _fields( 126 iounits_mod.read(iread), 127 [('a', 4), ('x', 1), ('i', 5), ('i', 5)] + [('f', 8, 5)] * 3 128 + [('x', 1), ('f', 14, 10), ('f', 13, 9), ('f', 13, 0)]) 129 130 # Flags (FORMAT 30) (Fortran lines 102–105). 131 # Suspected bug (see the page): 5X skips one column more than PUTYREC7 132 # writes; kept bug for bug. 133 lrot, lexcom, jcore, jenv, cmixl = _fields( 134 iounits_mod.read(iread), 135 [('l', 1), ('x', 1), ('l', 1), ('x', 1), ('x', 5), ('i', 4), 136 ('i', 4), ('f', 7, 3)]) 137 138 # Luminosities by source (FORMAT 40); in erg/s if any exceeds 1e20 139 # (Fortran lines 107–119). 140 tlumx[1:8] = _fields(iounits_mod.read(iread), 141 [('x', 10)] + [('f', 10, 3)] * 7) 142 cccmax = max(tlumx[1], tlumx[2], tlumx[3], tlumx[4], tlumx[5], 143 abs(tlumx[6]), tlumx[7]) 144 if cccmax > 1.0e20: 145 for j in range(1, 7 + 1): 146 tlumx[j] = tlumx[j] / constants_mod.clsun 147 148 # The envelope triangle (FORMAT 70) (Fortran lines 121–130). 149 ftri = 1.0 150 for i in range(1, 3 + 1): # label 80 151 vals = _fields(iounits_mod.read(iread), 152 [('x', 3), ('i', 2), ('f', 7, 5)] + [('f', 8, 5)] * 4 153 + [('f', 12, 5)] * 3) 154 io, trit[i], tril[i], ps[i], ts[i], rs[i] = vals[:6] 155 for j in range(1, 3 + 1): 156 cfenv[i + i + i - 3 + j] = vals[5 + j] 157 if io < 0: 158 ftri = -1.0 159 160 # The shells (FORMAT 100), X and Z in units of 1e-6 (Fortran lines 161 # 132–150). 162 spec100 = ([('f', 13, 10), ('f', 10, 7), ('f', 14, 7), ('f', 10, 7), 163 ('f', 10, 7), ('f', 10, 7), ('l', 1), ('i', 6), ('i', 6)]) 164 for i in range(1, m + 1): # label 110 165 (hs[i], hr[i], hl[i], hp[i], ht[i], hd[i], lc[i], ix, 166 iz) = _fields(iounits_mod.read(iread), spec100) 167 if hs[i] < 0.0 or hs[i] > hstot: 168 iounits_mod.write(ishort, ' ' + '>' * 39 + '<' * 40) 169 iounits_mod.write(ishort, ' RUN STOPPED') 170 iounits_mod.write(ishort, ' ERROR IN SUBROUTINE GETY7') 171 iounits_mod.write(ishort, f' GLITCH IN SHELL{i:3d}, SHELL MASS ' 172 'LESS THAN ZERO OR GREATER THAN STAR MASS') 173 raise SystemExit(3) 174 hcomp[1, i] = 1.0e-6 * float(ix) 175 hcomp[3, i] = 1.0e-6 * float(iz) 176 177 # Species 4-11, stored from jcore to jenv only (FORMAT 200) (Fortran 178 # lines 152–183). 179 spec200 = [('f', 9, 3), ('x', 1)] * 8 180 if jenv <= jcore: 181 hcomp[4:12, jcore] = _fields(iounits_mod.read(iread), spec200) 182 hcomp[4:12, jenv] = _fields(iounits_mod.read(iread), spec200) 183 else: 184 for j in range(jcore, jenv + 1): # label 210 185 hcomp[4:12, j] = _fields(iounits_mod.read(iread), spec200) 186 if jcore > 1: 187 # The homogeneous core. 188 for j in range(1, jcore): # label 230 189 for i in range(4, 11 + 1): # label 220 190 hcomp[i, j] = hcomp[i, jcore] 191 if jenv < m: 192 # The homogeneous envelope. 193 for j in range(jenv + 1, m + 1): # label 250 194 for i in range(4, 11 + 1): # label 240 195 hcomp[i, j] = hcomp[i, jenv] 196 for i in range(1, m + 1): # label 260 197 hcomp[2, i] = 1.0 - hcomp[1, i] - hcomp[3, i] - hcomp[4, i] 198 199 # The light elements (species 12-15), two shells per record (Fortran 200 # lines 185–223). 201 if lexcom: 202 spec300 = [('f', 9, 3), ('x', 1)] * 4 203 if jenv == 1: 204 # Fully convective: one record for the whole star. 205 hcomp[12:16, 1] = _fields(iounits_mod.read(iread), spec300) 206 for j in range(1, m + 1): # label 310 207 for i in range(12, 15 + 1): # label 305 208 hcomp[i, j] = hcomp[i, 1] 209 else: 210 for k in range(jcore, jenv, 2): # label 315 211 vals = _fields(iounits_mod.read(iread), spec200) 212 hcomp[12:16, k] = vals[0:4] 213 hcomp[12:16, k + 1] = vals[4:8] 214 # Fortran's K after the loop: one step past its last value. 215 kt = jenv - 1 - jcore 216 k = jcore + 2 * max(0, int((kt + 2) / 2)) 217 if kt % 2 != 0: 218 # An odd number of shells: the last one alone. 219 hcomp[12:16, k] = _fields(iounits_mod.read(iread), spec300) 220 if jcore > 1: 221 for j in range(1, jcore): # label 330 222 for i in range(12, 15 + 1): # label 320 223 hcomp[i, j] = hcomp[i, jcore] 224 if jenv < m: 225 for j in range(jenv + 1, m + 1): # label 350 226 for i in range(12, 15 + 1): # label 340 227 hcomp[i, j] = hcomp[i, jenv] 228 229 # Rotation (FORMAT 500): -log Omega, 58.9 or more for Omega = 0 230 # (Fortran lines 225–238); missing records go to the bypass (lines 231 # 245–260). 232 if lrot: 233 omegal = _read_values(iread, m, [('f', 10, 7)] * 8) 234 if omegal is None: # END=9999 235 iounits_mod.write(ishort, 236 "GETYREC7: OMEGA records are missing from " 237 " input model file - OMEGA array zeroed") 238 for i in range(1, json + 1): 239 omega[i] = 0.0 240 iounits_mod.rewind(iread) 241 return (bl, cmixl, dage, ddage, ftri, hstot, jcore, jenv, 242 lexcom, lrot, m, model, smass, teffl, atm, eos, hik, 243 ldify, ldifz, ldisk, linstb, ljdot0, alok, lovstc, 244 lovste, lovstm, lpurez, lsemic, compmix, pdisk, tdisk, 245 wmax) 246 for i in range(1, m + 1): # label 510 247 if omegal[i - 1] < 58.9: 248 omega[i] = 10.0 ** (-omegal[i - 1]) 249 else: 250 omega[i] = 0.0 251 252 # Keep iread open, rewound (Fortran lines 240–243). 253 iounits_mod.rewind(iread) 254 return (bl, cmixl, dage, ddage, ftri, hstot, jcore, jenv, lexcom, lrot, 255 m, model, smass, teffl, atm, eos, hik, ldify, ldifz, ldisk, 256 linstb, ljdot0, alok, lovstc, lovste, lovstm, lpurez, lsemic, 257 compmix, pdisk, tdisk, wmax)