putyrec7.f

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

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

Variables on this page (2)

clsun, json

Purpose

Writes a stellar model in the older "YREC7" text format (LLP 4/16/03), which YREC can still read as a starting model. PUTMODEL2 writes the newer "MODEL2" format.

How it fits

Inputs: the model and its global quantities, as arguments (as for PUTMODEL2, without the physics labels). The comments "Last two lines are MODEL2 add-ons" (line 31 ↗) and "Any error messages go to LU ISHORT" (line 51 ↗) describe the argument list before ISHORT and the MODEL2 add-ons were removed from it (lines 9–10 ↗; lines 28–29 ↗).
Who calls it: no routine in src/ calls it. It is kept on purpose: KNOWN_ISSUES.md (issue 24, lines 1658–1666) records the developers' ruling to keep it, since its reader is live and the format may have legacy use, while fourteen other unreachable routines were removed (settled from YREC's public sources, Nov 9, 2026).
Who reads the format: STARIN recognizes the leading NMOD and calls GETYREC7.
Output: the file on unit IWRITE. The routine also changes two of the caller's arrays: TLUMX is converted to erg s⁻¹ in place when it is in solar units (lines 71–77 ↗), and OMEGA is floored at $10^{-59}$ (line 131 ↗).

The format. A header line NMOD (model number, shells, mass, $\log T_\mathrm{eff}$, $\log L$, $\log M$, age, step); a flags line (LROT, LEXCOM, JCORE, JENV, CMIXL); a TLUMX line (erg s⁻¹); three ENV lines; one line per shell ($\log m$, $\log r$, $L$, $\log P$, $\log T$, $\log\rho$, convection flag, and $X$ and $Z$ as integers in units of $10^{-6}$); the abundances 4–11 only from JCORE to JENV (each convection zone stores one composition); with LEXCOM, abundances 12–15 two shells per line; with LROT, $|\log_{10}\Omega|$, eight per line. For a fully convective star (JENV ≤ JCORE), two lines of abundances 4–11 are written, for JCORE and JENV (lines 102–104 ↗), although the comment says "write out 1 comp" (line 101 ↗); GETYREC7 reads two as well (getyrec7.f, lines 156–157 ↗), so only the comment is wrong (settled from the code, Nov 9, 2026).

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) The flags record does not match its reader.

(lines 65–66 ↗; getyrec7.f, lines 104–105 ↗) This routine writes the flags as 2(L1,1X),'F F ',2I4,F7.3, four characters between the flags and JCORE; GETYREC7 reads them as 2(L1,1X),5X,2I4,F7.3, skipping five, so every field it reads is shifted by one column. With Fortran's default treatment of blanks in input (ignored), the values still come back right while JENV < 1000 and CMIXL < 100; when JENV has four digits, its first digit is read as the last digit of JCORE, and JENV loses it. Only when a model written by PUTYREC7, which nothing in src/ calls, is read back.

A possible fix (in the writer, since GETYREC7 presumably matches the files YREC7 itself wrote):

putyrec7.f

@@ lines 63–69 @@

 C write PHYSICS FLAGS- ROTATION,EXTENDED COMP,CENTRAL AND SURFACE CZ'S
       WRITE(IWRITE,30)LROT,LEXCOM,JCORE,JENV,CMIXL
-   30 FORMAT(2(L1,1X),'F F ',2I4,F7.3)
+C 11/09/26 KC via Claude
+C code changed to write five characters between the flags and JCORE,
+C as GETYREC7 reads them.
+   30 FORMAT(2(L1,1X),'F F  ',2I4,F7.3)

 C write LUMINOSITIES
 C If TLUMX are in solar units, convert to ergs.  Decide by

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). Nothing in src/ calls PUTYREC7, so it cannot be reached at 7be4133; it is fragile code. Intent: KNOWN_ISSUES.md keeps the routine on purpose ("it may have legacy use even though nothing in the tree calls it") but says nothing of its format. Not run: a model written by PUTYREC7 and read back by GETYREC7, which needs a caller.

Definitions

subroutinePUTYREC7

lines 8–146 · GitHub

Uses modules CONSTANTS_MOD, DIMENSIONS_MOD

Change log

The changes since rp_1

Changed since rp_1. Declarations only: /CONST/ to CONSTANTS_MOD, JSON from DIMENSIONS_MOD, IMPLICIT NONE.

1C 2C 3C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 4C PUTYREC7 - Write out a stallar model in YREC7 format 5C 6C llp 4/16/03 7C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 8 SUBROUTINE PUTYREC7(BL,CFENV,CMIXL,DAGE,DDAGE,FTRI,HCOMP,HD,HL, 9C * HP,HR,HS,HSTOT,HT,IWRITE,ISHORT,JCORE,JENV,LC,LEXCOM,LROT,M, 10C * MODEL,OMEGA,PS,RS,SMASS,TEFFL,TLUMX,TRIL,TRIT,TS, ! KC 2025-05-31 11 * * HP,HR,HS,HSTOT,HT,IWRITE,JCORE,JENV,LC,LEXCOM,LROT,M, 12 * * MODEL,OMEGA,PS,RS,SMASS,TEFFL,TLUMX,TRIL,TRIT,TS) 13C 10/15/26 MHP via Claude 14C code changed to take JSON from DIMENSIONS_MOD instead of declaring 15C PARAMETER(JSON=5000) here. 16 USE DIMENSIONS_MOD 17C 10/14/26 MHP via Claude 18C code changed to convert COMMON/CONSTANTS/ into CONSTANTS_MOD. 19 USE CONSTANTS_MOD 20C 10/20/26 MHP via Claude 21C code changed to replace implicit typing with IMPLICIT NONE and 22C explicit declarations. 23 IMPLICIT NONE 24 REAL*8 BL,CCCMAX,CFENV,CMIXL,DAGE,DDAGE,FTRI,HCOMP,HD,HL,HP,HR,HS, 25 * * HSTOT,HT,OMEGA,OMEGAL,PS,RS,SMASS,TEFFL,TLUMX,TRIL,TRIT,TS 26 INTEGER I,I0,IWRITE,IX,IZ,J,J1,JCORE,JENV,K,KTEST,M,MODEL,N1,N2 27 LOGICAL*4 LC,LEXCOM,LROT 28C & ATM,EOS,HIK,LDIFY,LDIFZ,LDISK,LINSTB,LJDOT0,ALOK, 29C & LOVSTC,LOVSTE,LOVSTM,LPUREZ,LSEMIC,COMPMIX,PDISK,TDISK,WMAX) ! KC 2025-05-31 30C First three lines above are YREC7 inputs 31C Last two lines are MODEL2 add-ons 32 33C WRITE MODEL OUT IN ASCII FORMAT 34 35C 9/25/26 MHP via Claude 36C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/, 37C and COMMON/CONST3/ into a single COMMON/CONSTANTS/ 38 39c CHARACTER*4 ATM, LOK, HIK, COMPMIX 40c MHP 4/25 changed LOK name to make it unique, used elsewhere 41C CHARACTER*4 ATM,HIK,ALOK,COMPMIX 42C CHARACTER*6 EOS 43 44 DIMENSION HCOMP(15,JSON),HD(JSON),HL(JSON),HP(JSON),HR(JSON), 45 * * HS(JSON),HT(JSON),LC(JSON),TRIT(3),TRIL(3),PS(3),TS(3),RS(3), 46 * * CFENV(9),TLUMX(8),OMEGA(JSON),OMEGAL(JSON) 47 48 save 49 50C PUTYREC7 writes the model provided to it out to the LU IWRITE file. 51C Any error messages go to LU ISHORT. 52 53 54C THE 6 HEADER RECORDS WRITTEN FIRST 55 56 IF (DDAGE .GE. 1D6) THEN 57 WRITE(IWRITE,10) MODEL,M,SMASS,TEFFL,BL,HSTOT,DAGE,DDAGE 58 10 10 FORMAT('NMOD ',2I5,3F8.5,1X,F14.10,F13.9,F13.0) 59 ELSE 60 WRITE(IWRITE,20) MODEL,M,SMASS,TEFFL,BL,HSTOT,DAGE,DDAGE 61 20 20 FORMAT('NMOD ',2I5,3F8.5,1X,F14.10,F13.9,1PE13.6) 62 ENDIF 63 64C write PHYSICS FLAGS- ROTATION,EXTENDED COMP,CENTRAL AND SURFACE CZ'S 65 WRITE(IWRITE,30)LROT,LEXCOM,JCORE,JENV,CMIXL 66 30 30 FORMAT(2(L1,1X),'F F ',2I4,F7.3) 67 68C write LUMINOSITIES 69C If TLUMX are in solar units, convert to ergs. Decide by 70C comparing to 10**20. If smaller, multiply by CLSUN. 71 CCCMAX = DMAX1(TLUMX(1),TLUMX(2),TLUMX(3),TLUMX(4),TLUMX(5), 72 * * DABS(TLUMX(6)),TLUMX(7)) 73 IF(CCCMAX.LE.1.0D20) THEN 74 DO J = 1,7 75 TLUMX(J) = TLUMX(J) * CLSUN 76 ENDDO 77 ENDIF 78 WRITE(IWRITE,40) (TLUMX(J),J=1,7) 79 40 40 FORMAT('TLUMX',5X,1P7E10.3) 80 81C write ENVELOPE DATA 82 DO 60 I = 1,3 83 I0 = I 84 IF(FTRI.LT.0D0) I0 = -I 85 WRITE(IWRITE,50)I0,TRIT(I),TRIL(I),PS(I),TS(I),RS(I), 86 * *(CFENV(3*I-3+J),J=1,3) 87 50 50 FORMAT('ENV',I2,F7.5,4F8.5,1P3E12.5) 88 60 60 CONTINUE 89 90C write HENYEY POINTS, one line per shell 91 DO 110 I = 1,M 92 IX = IDINT(1.0D6*HCOMP(1,I) + 0.50D0) 93 IZ = IDINT(1.0D6*HCOMP(3,I) + 0.50D0) 94 WRITE(IWRITE,100) HS(I),HR(I),HL(I),HP(I),HT(I),HD(I),LC(I), 95 * * IX,IZ 96 100 100 FORMAT(0PF13.10,F10.7,1PE14.7,0PF10.7,2F10.7,L1,2I6) 97 110 110 CONTINUE 98 99C WRITE OUT COMPOSITION ARRAY - CENTRAL AND SURFACE 100C CONVECTION ZONES HAVE ONLY 1 COMPOSITION STORED PER ZONE 101C IF STAR IS FULLY CONVECTIVE, WRITE OUT 1 COMP. 102 IF(JENV.LE.JCORE) THEN 103 WRITE(IWRITE,200)(HCOMP(I,JCORE),I=4,11) 104 WRITE(IWRITE,200)(HCOMP(I,JENV),I=4,11) 105 ELSE 106 DO 210 J = JCORE,JENV 107 WRITE(IWRITE,200) (HCOMP(I,J),I = 4,11) 108 200 200 FORMAT(8(1PE9.3,1X)) 109 210 210 CONTINUE 110 ENDIF 111C EXTENDED COMP - WRITE OUT 1 ABUND IF FULLY CONVECTIVE, 112C OTHERWISE WRITE OUT 2 POINTS PER LINE. 113 IF(LEXCOM) THEN 114 IF(JENV.LE.JCORE) THEN 115 WRITE(IWRITE,220)(HCOMP(I,JENV),I=12,15) 116 220 220 FORMAT(4(1PE9.3,1X)) 117 ELSE 118 DO 230 K = JCORE,JENV-1,2 119 WRITE(IWRITE,200)((HCOMP(I,J),I = 12,15),J = K,K+1) 120 230 230 CONTINUE 121C IF AN ODD NUMBER OF ABUNDANCES EXISTS, WRITE IN LAST VALUE 122 KTEST = JENV-1 - JCORE 123 IF(MOD(KTEST,2).NE.0) 124 * * WRITE(IWRITE,220)(HCOMP(I,K),I = 12,15) 125 ENDIF 126 ENDIF 127C WRITE OUT RUN OF OMEGA,STORED 8 ELEMENTS PER LINE 128C (- THE LOG OF OMEGA IS STORED) 129 IF(LROT) THEN 130 DO 300 I = 1,M 131 IF (OMEGA(I).LE.1.D-59) OMEGA(I)=1.D-59 132 OMEGAL(I) = DABS(DLOG10(OMEGA(I))) 133 300 300 CONTINUE 134 N1 = INT(M/8) 135 N2 = M - N1*8 136 J1 = 1 137 DO 320 I = 1,N1 138 WRITE(IWRITE,310)(OMEGAL(J),J = J1,J1 + 7) 139 310 310 FORMAT(0P8F10.7) 140 J1 = J1 + 8 141 320 320 CONTINUE 142 IF(N2.GT.0) WRITE(IWRITE,310)(OMEGAL(J),J=J1,J1+N2-1) 143 ENDIF 144 145 RETURN 146 END

1"""Python version of putyrec7.f (YREC 7be4133): write a model in the older 2YREC7 text format.""" 3 4import math 5 6import numpy as np 7 8from pyrec import constants_mod 9from pyrec import dimensions_mod 10from pyrec import iounits_mod 11 12# The local array OMEGAL (SAVE). 13omegal = np.zeros(dimensions_mod.json + 1) 14 15 16def _l1(flag: bool) -> str: 17 """Fortran's L1.""" 18 19 return "T" if flag else "F" 20 21 22def _e1(x: float, w: int, d: int) -> str: 23 """Fortran's 1PEw.d: one digit before the point, e.g. 1.234E+05.""" 24 25 return f"{x:{w}.{d}E}" 26 27 28def putyrec7(bl: float, cfenv: np.ndarray, cmixl: float, dage: float, ddage: float, ftri: float, hcomp: np.ndarray, 29 hd: np.ndarray, hl: np.ndarray, hp: np.ndarray, hr: np.ndarray, hs: np.ndarray, hstot: float, 30 ht: np.ndarray, iwrite: int, jcore: int, jenv: int, lc: np.ndarray, lexcom: bool, lrot: bool, m: int, 31 model: int, omega: np.ndarray, ps: np.ndarray, rs: np.ndarray, smass: float, teffl: float, 32 tlumx: np.ndarray, tril: np.ndarray, trit: np.ndarray, ts: np.ndarray) -> None: 33 """Write the model to unit iwrite in the YREC7 format: header, flags, 34 TLUMX (erg/s), three ENV lines, one line per shell (X and Z as 35 integers in units of 1e-6), the abundances 4-11 from jcore to jenv, 36 with lexcom the abundances 12-15, with lrot |log10 Omega|. 37 38 Translates: putyrec7.f#PUTYREC7 (7be4133) 39 40 Arguments: as in the Fortran (arrays with index 0 unused; trit, tril, 41 ps, ts, rs [1..3], cfenv [1..9], tlumx [1..8]). tlumx (converted from 42 solar units) and omega (floored at 1e-59) are changed in place. 43 Returns nothing. 44 45 """ 46 47 w = iounits_mod.write 48 49 # The header record (Fortran lines 56–62). 50 head = (f"NMOD {model:5d}{m:5d}{smass:8.5f}{teffl:8.5f}{bl:8.5f} " 51 f"{hstot:14.10f}{dage:13.9f}") 52 if ddage >= 1e6: 53 w(iwrite, head + f"{ddage:#13.0f}") # F13.0 keeps the point 54 else: 55 w(iwrite, head + _e1(ddage, 13, 6)) 56 57 # The flags (Fortran lines 65–66). Suspected bug (see the page): 58 # GETYREC7 reads this record with 5X, not 4. 59 w(iwrite, f"{_l1(lrot)} {_l1(lexcom)} F F {jcore:4d}{jenv:4d}" 60 f"{cmixl:7.3f}") 61 62 # The luminosities, converted to erg/s in the caller's array if in 63 # solar units (Fortran lines 71–79). 64 cccmax = max(tlumx[1], tlumx[2], tlumx[3], tlumx[4], tlumx[5], 65 abs(tlumx[6]), tlumx[7]) 66 if cccmax <= 1.0e20: 67 for j in range(1, 8): 68 tlumx[j] = tlumx[j] * constants_mod.clsun 69 w(iwrite, "TLUMX" + " " * 5 70 + "".join(_e1(tlumx[j], 10, 3) for j in range(1, 8))) 71 72 # The envelope data (Fortran lines 82–88). 73 for i in range(1, 4): 74 i0 = i 75 if ftri < 0.0: 76 i0 = -i 77 w(iwrite, f"ENV{i0:2d}{trit[i]:7.5f}{tril[i]:8.5f}{ps[i]:8.5f}" 78 f"{ts[i]:8.5f}{rs[i]:8.5f}" 79 + "".join(_e1(cfenv[3 * i - 3 + j], 12, 5) for j in range(1, 4))) 80 81 # One line per shell; X and Z as integers in units of 1e-6 82 # (Fortran lines 91–97). Question (see the page): the precision. 83 for i in range(1, m + 1): 84 ix = int(1.0e6 * hcomp[1, i] + 0.50) 85 iz = int(1.0e6 * hcomp[3, i] + 0.50) 86 w(iwrite, f"{hs[i]:13.10f}{hr[i]:10.7f}" + _e1(hl[i], 14, 7) 87 + f"{hp[i]:10.7f}{ht[i]:10.7f}{hd[i]:10.7f}{_l1(lc[i])}" 88 f"{ix:6d}{iz:6d}") 89 90 def fmt200(values: list) -> str: 91 # FORMAT(8(1PE9.3,1X)); more than 8 values continue on a new line. 92 vals = [_e1(v, 9, 3) + " " for v in values] 93 return "\n".join("".join(vals[n:n + 8]) for n in range(0, len(vals), 8)) 94 95 # The abundances 4-11 from jcore to jenv; two lines for a fully 96 # convective star, as GETYREC7 reads (Fortran lines 102–110). 97 if jenv <= jcore: 98 w(iwrite, fmt200([hcomp[i, jcore] for i in range(4, 12)])) 99 w(iwrite, fmt200([hcomp[i, jenv] for i in range(4, 12)])) 100 else: 101 for j in range(jcore, jenv + 1): 102 w(iwrite, fmt200([hcomp[i, j] for i in range(4, 12)])) 103 104 # With lexcom, the abundances 12-15, two shells per line (Fortran 105 # lines 113–126). 106 if lexcom: 107 if jenv <= jcore: 108 w(iwrite, fmt200([hcomp[i, jenv] for i in range(12, 16)])) 109 else: 110 for k in range(jcore, jenv, 2): 111 w(iwrite, fmt200([hcomp[i, j] for j in (k, k + 1) 112 for i in range(12, 16)])) 113 # Fortran's K after the loop: one step past the last value. 114 ntrip = max(0, (jenv - 1 - jcore) // 2 + 1) 115 k = jcore + 2 * ntrip 116 ktest = jenv - 1 - jcore 117 if ktest % 2 != 0: 118 w(iwrite, fmt200([hcomp[i, k] for i in range(12, 16)])) 119 120 # With lrot, |log10 Omega|, eight per line; omega is floored in the 121 # caller's array (Fortran lines 129–143). Question (see the page): the 122 # sign of log Omega is lost. 123 if lrot: 124 for i in range(1, m + 1): 125 if omega[i] <= 1.0e-59: 126 omega[i] = 1.0e-59 127 omegal[i] = abs(math.log10(omega[i])) 128 n1 = m // 8 129 n2 = m - n1 * 8 130 j1 = 1 131 for i in range(1, n1 + 1): 132 w(iwrite, "".join(f"{omegal[j]:10.7f}" for j in range(j1, j1 + 8))) 133 j1 = j1 + 8 134 if n2 > 0: 135 w(iwrite, "".join(f"{omegal[j]:10.7f}" 136 for j in range(j1, j1 + n2)))