putstore.f

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

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

Variables on this page (48)

cc13, cgl, cln, clsun, derot, gam1, ishort, istor, json, kttau, ldh, ldify, ldifz, ldisk, lexcom, linstb, ljdot0, lmilne, lopal95, lopale06, lovstc, lovste, lovstm, lphhd, lpurez, lrot, lscv, lsemic, lstatm, lstch, lstenv, lstmod, lstore, lstphys, lstrot, pdisk, sbeta, sdel, seg, sesum, seta, so, svel, tdisk, ves, vgsf, vss, wmax

Purpose

Writes a verbose stellar model to the .store file (ISTOR), "either at specified ages, every NPRTMOD models, or at the end of runs." Its header (G. Somers 11/14) calls it "mostly a clone of wrtlst.f and putmodel2.f, combined." Since JvS 09/25 it has two formats: the older one (a MODEL2-style model followed by optional per-shell physics and rotation columns) and, with LSTCH, the "stitched" format, in which STITCH writes interior, envelope and atmosphere points together.

How it fits

Inputs: the model as arguments (structure, composition, envelope data, luminosities TLUMX, rotation: OMEGA, ETA2, R0, FP, FT, HJM, HI), the mixing length CMIXL, XMSOL (passed on to WRTMOD) and the initial-mixture label COMPMIX (from STARIN, through the main program); the physics flags and the output switches (LSTORE, LSTATM, LSTENV, LSTMOD, LSTPHYS, LSTROT, LSTCH, LMILNE) from NMLPAR_MOD, and LPHHD from CCOUT_MOD; the per-shell physics in SCRTCH_MOD, SOUND_MOD (GAM1), TEMP2_MOD (the rotational mixing velocities) and ROTEN_MOD; the constants from CONSTANTS_MOD and the units from IOUNITS_MOD.
Who calls it: WRTOUT and the main program in main.f.
What it calls: STITCH (new format), PINDEX (which shells to print), WRTMIL (if also LMILNE) and WRTMOD (envelope and atmosphere, if LSTATM or LSTENV).
Output: the .store file, opened by PARMIN (which also writes the old format's header key). Its MOD2 blocks can serve as starting models (STARIN reads the MOD2 format).

The method

Labels for the physics: atmosphere (EDD, KS, HRA, KUR, ALL from KTTAU; KTTAU = 5 gets none, bug 2), equation of state (SAHA, OPAL06, SCV, with +DH and OPAL 2006; SCV+O6 is never written, bug 1), low-temperature opacity (always NONE: the labels ALEX and KURZ went with the Alexander 1995 and Kurucz 1990 tables) and high-temperature opacity (OP95, or NONE). The message for the label ' ? ' (lines 176–180 ↗) is never written, since no branch sets that label (settled from the code, Nov 9, 2026).

Question 1 (lines 164–169 ↗). No low-temperature opacity label. With the Alexander 1995 and Kurucz 1990 tables removed, ALOK is always NONE: the Ferguson et al. (2005) tables (LALEX06), the only low-temperature tables left, were never labelled. Should LALEX06 set ALOK (for example to AL06)?
Stitched format (LSTCH): on the first call (LPHHD) a header key, the flags line and the column titles; then STITCH (with CMIXL).
Old format: the MOD2 header line, the flags, the TLUMX line (converted to erg s⁻¹ in place if the largest is $\le 10^{20}$, as in PUTMODEL2, whose in-place conversion is putmodel2.f, suspected bug 1), three ENV lines; then, with LSTMOD, one line per shell selected by PINDEX (whose burning-shell arguments are never set, bug 3): the basic columns, led by the model number; with LSTPHYS opacity, gravity $g = Gm/r^2$, $\nabla_\mathrm{rad}$, $\nabla$, $\nabla_\mathrm{ad}$, convective velocity, $\Gamma_1$, SBETA, SETA and the energy generation terms; with LSTROT (and LROT) the distortion parameter and the ratio of polar to equatorial radius, $$A = \frac{5}{3}\,\frac{\Omega^2 r_0^3}{G m\,(2+\eta_2)},\qquad \frac{R_p}{R_e} = \frac{1 - A}{1 + A/2},$$ where $A$ is the amplitude of the $P_2$ distortion of the level surface, $r(\theta) = r_0\,[1 - A\,P_2(\cos\theta)]$, in the Clairaut–Radau treatment of FPFT and SHAPE (Kippenhahn & Thomas 1970): $P_2 = 1$ at the pole and $-\tfrac12$ at the equator, hence $r_0(1 - A)$ and $r_0(1 + A/2)$; then FP, FT, $j$, $I$, the rotational kinetic energy change, and the mixing velocities $v_\mathrm{ES}$, $v_\mathrm{GSF}$, $v_\mathrm{SS}$ and their sum. Missing groups are written as zeros. The fifth energy column, SEG(5,·), is the triple-α term: COEFFT stores in it EG(5), ENGEB's E3AL (coefft.f, lines 286–296 ↗; coefft.f, lines 433–435 ↗; engeb.f, lines 1125–1127 ↗), as the stitched header's "TRIPLE_ALPHA" says; the old header's "3HE" is wrong (bug 4) (settled from the code, Nov 9, 2026).

Question 2 (lines 307–312 ↗). Ionization fractions written as zeros. The columns HII, HEII, HEIII are written as 0.0 (the line with SFXION is commented out), although COEFFT fills SFXION for every shell (coefft.f, line 450 ↗) and STITCH writes it (stitch.f, lines 117–120 ↗). Is that intended?

Question 3 (lines 133–174 ↗; lines 213–251 ↗). Duplicated code. The labels, the MOD2 header, and the TLUMX conversion repeat PUTMODEL2 and WRTLST, and the copies have already drifted (KTTAU = 5 is labelled in WRTLST, not here). Should one routine be shared?

For users (from YREC's Quick Start). This is the .store file (FSTOR): snapshots of the model's structure during the run. Which snapshots, and what they contain, the control namelist decides: LSTORE (NPRTMOD, every how many models; NPRTPT, every how many shells), LSTPCH (the final model), and LSTATM, LSTENV, LSTPHYS, LSTROT (the atmosphere, the envelope, the microphysics, the rotation), with LSTCH for the stitched format.

From YREC's documentation

Control Options (Output):

variable	meaning (documentation)
`LSTORE`	store model structures to the `.store` file at the chosen points (`LSTMOD` should be set the same; it may be redundant)
`LSTCH`	write interior, envelope and atmosphere in one consistent format (recommended)
`LSTATM`, `LSTENV`	add the atmosphere, the envelope; with both off nothing beyond the fitting point is written
`LSTPHYS`	add the microphysics ($\epsilon$, $\kappa$, …)
`LSTROT`	add the rotation data
`LSTPCH`	write the final model's details

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 SCV + OPAL 2006 label is lost.

(lines 149–157 ↗) In the LSCV branch, IF (LOPALE06) EOS='SCV+O6' sits inside IF (LDH), where it is overwritten at once; without LDH, SCV with OPAL 2006 is labelled SCV in the flags line of both formats. Only with LSCV and LOPALE06 on and LDH off. The same lines are in WRTLST (wrtlst.f, suspected bug 1).

A possible fix:

putstore.f

@@ lines 146–157 @@
 C Determine equation of state flag, EOS
       EOS='SAHA  '
       IF (LDH) EOS='SAH+DH'
+C 11/09/26 KC via Claude
+C code changed to label SCV with OPAL 2006 and without Debye-Huckel
+C SCV+O6; the test sat inside IF (LDH), where it was overwritten.
       IF (LSCV) THEN
          EOS='SCV   '
 C 9/24/26 MHP via Claude
 C code changed to remove the OPAL 1995/2001 equations of state
-         IF (LDH) THEN
          IF (LOPALE06) EOS='SCV+O6'
+         IF (LDH) THEN
             EOS='SCV+DH'
             IF (LOPALE06) EOS='SCDHO6'
          ENDIF

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 LSCV branch ran in 30 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, but the setting that goes wrong, SCV with OPAL 2006 and without LDH, held in none: YREC's cases that use SCV with OPAL 2006 also set LDH. Any of them with LDH off would reach it. The fix was not run. Intent: none found; the misplaced line is already in the first commit of the history (Apr 10, 2025), indented at the level of IF (LDH), as if pasted one line low.

2. (wrong output) Checked by KC (Nov 10, 2026) Atmosphere 5 has no label.

(lines 135–145 ↗) KTTAU = 5 (Kurucz/Castelli, JNT 6/14) has no branch, so ATM keeps its value from the previous call (it is SAVEd) or, at the first call, is undefined, and that is written to the .store file's flags line. WRTLST has the branch (wrtlst.f, lines 103–105 ↗). Only with KTTAU = 5 (the default is 0) and LSTORE.

A possible fix:

putstore.f

@@ lines 142–147 @@
          ATM='KUR '
       ELSEIF (KTTAU .EQ. 4) THEN
          ATM='ALL '
+C 11/09/26 KC via Claude
+C code changed to label the Kurucz/Castelli atmosphere, as WRTLST does.
+      ELSEIF (KTTAU .EQ. 5) THEN
+         ATM='K/C '
       ENDIF
 C Determine equation of state flag, EOS
       EOS='SAHA  '

How the outputs show it. In the made-up run (testsuite/Test_solarGS98_kurucz with KTTAU = 5 and FATM = atmos/CastelliKurucz/atmk2004p00.tab), the .store file's flags line, written once, at the first store (model 23; lines 189–192 ↗), holds four NUL bytes in the ATM column, between "SCDHO6" and "NONE": ATM was never set. The .last file, written at the last model (1224) by WRTLST (wrtlst.f, lines 103–105 ↗) through PUTMODEL2 (putmodel2.f, line 80 ↗), gives "K/C" in the same column.

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): YREC's testsuite/Test_solarGS98_kurucz, run with KTTAU = 5 and FATM = atmos/CastelliKurucz/atmk2004p00.tab, as the README of those tables asks ("You NEED TO CHANGE KTTAU to 5," input/atmos/CastelliKurucz/README, line 15), reaches it at both of its stores (models 23 and 1000), so ATM was written without a value of its own. With the fix, the track is identical and the harness finds the .store file within 2%; the null change moves $\log L$ by up to 0.000144 dex. The harm is the label, which the fix corrects. Intent: none found; WRTLST got the branch ("JNT 06/14"), and PUTSTORE, written later (G Somers, 11/14), lacks it.

3. (never occurred) Checked by KC (Nov 10, 2026) `PINDEX` gets arguments that are never set.

(lines 51–52 ↗; line 262 ↗) JXBEG, JXEND, and LSHELL are neither arguments nor in a module here, so they are undefined locals (zero and .FALSE. with a compiler that zero-initializes SAVEd locals). The burning-shell option of PINDEX, which prints every shell of the hydrogen-burning shell (pindex.f, lines 46–70 ↗), therefore never applies, and the old format prints only every NPRTPT-th shell there. Only in the old format (LSTCH off) with LSTMOD, and with NPRTPT above 1 (its default is 1, nmlpar_mod.f, line 169 ↗, kept by every case of YREC's that writes the old format).

A possible fix: pass LSHELL, JXBEG, and JXEND to PUTSTORE as arguments, from WRTOUT (which has them as arguments) and from the main program (which sets them through FINDSH, main.f, lines 578–579 ↗), as they were passed to WRTMOD in the commented-out call (lines 344–346 ↗). The branch this reaches has two fragile points of its own (pindex.f, suspected bug 1; pindex.f, suspected bug 2).

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): PINDEX got JXBEG = JXEND = 0 and LSHELL false at every call, in the 14 runs that write the old format. In 10 of them, the callers' LSHELL was true at a store, with the shell between points 26 and 2830. But no shell went unlisted: Every case of YREC's that writes the old format keeps NPRTPT at 1 (its default, nmlpar_mod.f, line 169 ↗), and with 1, PINDEX's general branch lists every shell, 2 to M−1 and then M (pindex.f, lines 73–84 ↗); testsuite/Test_m2p0_feh+0p0_vzams150_IES2_TAMS lists all 748 shells of its model 921 once, and testsuite/Test_m0p95feh+0p0_GN93_TAHB all 3,268 and 3,425 shells of its models 0 and 943 (rechecked on Nov 10). The run with the fix differs, but not because of this bug: Its track is identical, and only the .store file changes beyond 2% (the null change moves $\log L$ by up to 1.67e-7 dex in the first case). With LSHELL passed, PINDEX takes its shell branch, which with NPRTPT = 1 also lists every shell once (loops 10, 20, and 30, pindex.f, lines 46–72 ↗) and then lists M a second time (pindex.f, line 83 ↗): That is PINDEX's own suspected bug, M listed twice (pindex.f, suspected bug 2), which nothing reaches today and this fix would. So this fix should go with that one. The bug itself shows only with NPRTPT above 1, which no case of YREC's sets. Intent: The commented-out call to WRTMOD (lines 344–346 ↗) passed these three; nothing says they were dropped on purpose.

4. (wrong output) Checked by KC (Nov 10, 2026) The old header calls the triple-α column "3HE."

(line 275 ↗; line 309 ↗) The fifth energy column, SEG(5,I), is the triple-α term (Explanation, step 3), as the stitched format's header says (line 208 ↗). What follows: a misleading column name in the old .store format. Only with LSTORE, LSTMOD, and LSTPHYS, and LSTCH off.

A possible fix:

putstore.f

@@ lines 272–278 @@
      1'    O17             O18             H2              Li6             Li7',
      1'             Be9           OPAC       GRAV        DELR        DEL      ',
      1'   DELA       V_CONV     GAM1      HII     HEII     HEIII    BETA      ',
-     1'ETA       PPI         PPII       PPIII        CNO         3HE         ',
+C 11/09/26 KC via Claude
+C code changed to name the fifth energy column for the triple-alpha
+C term it holds, as FORMAT 1014 does.
+     1'ETA       PPI         PPII       PPIII        CNO         3ALPHA      ',
      1'E_NUC        E_NEU       E_GRAV          A           RP/RE           FP',
      1'            FT           J/M          MOMENT        DEL_KE       V_ES ',
      1'      V_GSF      V_SS       VTOT   ',/)

How the outputs show it. In testsuite/Test_m0p95feh+0p0_GN93_TAHB, at model 943, the .store file's column headed "3HE" (line 275 ↗) is 0.1000E+01 at the center, where the core burns helium. This column's fractions of E_NUC, summed over the masses of all 3,425 shells, give 3.083 $L_\odot$ (with the case's CLSUN of 3.828×10³³ erg s⁻¹), and the .track gives the same model's triple-α luminosity, 3a_lsun (wrtout.f, line 558 ↗), as 3.08298493E+000. At model 0 the two are 16.16 and 1.61768489E+001. The column is the triple-α term, as the stitched header says (line 208 ↗).

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 written with the physics columns in the same 14 runs. With the fix, examples/run_from_dbl_to_zams/Test_m0040_feh0_noGS_ZAMS gives an identical track and a .store file within 2%, and examples/run_standard_solar_model/Test_solarcal_noGS_norot differs exactly as for the null change ($\log L$ by up to 0.000322 dex): within the noise, as expected for a label. The harm is the label, which the fix corrects. Intent: none found; the stitched format's header names the column TRIPLE_ALPHA (line 208 ↗).

Definitions

subroutinePUTSTORE

lines 13–360 · GitHub

Uses modules CCOUT_MOD, CONSTANTS_MOD, DIMENSIONS_MOD, IOUNITS_MOD, NMLPAR_MOD, ROTEN_MOD, SCRTCH_MOD, SOUND_MOD, TEMP2_MOD

Called by (file level), WRTOUT

Flowchart of its calls, in order (4)

if (LSTCH), lines 184–210

STITCH 198

else, lines 211–354

PINDEX 262

if (LSTATM.OR.LSTENV), lines 342–349

if (LMILNE): WRTMIL 343→WRTMOD 347

Change log

The changes since rp_1

Changed since rp_1. CMIXL, XMSOL and COMPMIX are new arguments (COMPMIX was /I2O/, CMIXL was in /CONST3/), passed on: CMIXL to STITCH, both to WRTMOD; the labels of the removed options (the OPAL 1995 and 2001 equations of state, the OPAL92, LAOL89, Alexander 1995 and Kurucz 1990 opacities) are gone; each shell's line now starts with the model number, which fixes a crash (format 62 already expected it), and formats 62 and 63 put spaces between some columns; LMILNE now comes from NMLPAR_MOD (bug fix), so WRTMIL is called when it is set, which settles former question 2; COMMON blocks to modules (eleven unused blocks dropped), IMPLICIT NONE.

1C 2C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 3C PUTSTORE - Write out a verbose stellar model to the .store file. This program 4C is mostly a clone of wrtlst.f and putmodel2.f, combined. 5C 6C G Somers 11/14 7C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 8C 9C 10C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 11C 06 12C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 13 SUBROUTINE PUTSTORE(HCOMP,HD,HL,HP,HR,HS,HT,LC,TRIT,TRIL,PS, 14 * *TS,RS,CFENV,FTRI,TLUMX,JCORE,JENV,MODEL,M,SMASS,TEFFL,BL,HSTOT, 15 * *DAGE,DDAGE,OMEGA,HS1,ETA2,R0,FP,FT,HJM,HI,CMIXL,XMSOL,COMPMIX) 16C 10/15/26 MHP via Claude 17C code changed to take JSON from DIMENSIONS_MOD instead of declaring 18C PARAMETER(JSON=5000) here. 19 USE DIMENSIONS_MOD 20C 10/14/26 MHP via Claude 21C code changed to convert COMMON/TEMP2/ into TEMP2_MOD. 22 USE TEMP2_MOD 23C 10/14/26 MHP via Claude 24C code changed to convert COMMON/SOUND/ into SOUND_MOD. 25 USE SOUND_MOD 26C 10/14/26 MHP via Claude 27C code changed to convert COMMON/SCRTCH/ into SCRTCH_MOD. 28 USE SCRTCH_MOD 29C 10/14/26 MHP via Claude 30C code changed to convert COMMON/ROTEN/ into ROTEN_MOD. 31 USE ROTEN_MOD 32C 10/14/26 MHP via Claude 33C code changed to convert COMMON/CCOUT/ into CCOUT_MOD. 34 USE CCOUT_MOD 35C 10/14/26 MHP via Claude 36C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD. 37 USE NMLPAR_MOD 38C 10/14/26 MHP via Claude 39C code changed to convert COMMON/CONSTANTS/ into CONSTANTS_MOD. 40 USE CONSTANTS_MOD 41C 10/14/26 MHP via Claude 42C code changed to convert COMMON/IOUNITS/ into IOUNITS_MOD. 43 USE IOUNITS_MOD 44C 10/20/26 MHP via Claude 45C code changed to replace implicit typing with IMPLICIT NONE and 46C explicit declarations. 47 IMPLICIT NONE 48 REAL*8 A,BL,CCCMAX,CFENV,CG,CMIXL,DAGE,DDAGE,DUMA,ETA2,FM,FP,FT, 49 * * FTRI,HCOMP,HD,HI,HJM,HL,HP,HR,HS,HS1,HSTOT,HT,OMEGA,PS,R0, 50 * * RPOLEQ,RS,SG,SMASS,TEFFL,TLUMX,TRIL,TRIT,TS,VTOT,XMSOL 51 INTEGER I,I0,ID,IDM,II,J,JCORE,JENV,JXBEG,JXEND,K,M,MODEL 52 LOGICAL*4 LC,LSHELL 53 54C PUTSTORE PUTS THE MOST RECENT VERBOSE OUTPUT FILE INTO THE .STORE FILE, 55C EITHER AT SPECIFIED AGES, EVERY NPRTMOD MODELS, OR AT THE END OF RUNS. 56 57C WRITE MODEL OUT IN ASCII FORMAT 58 CHARACTER*6 EOS 59c CHARACTER*4 ATM, LOK, HIK, COMPMIX 60c MHP 4/25 changed LOK name to make it unique, used elsewhere 61 CHARACTER*4 ATM, ALOK, HIK, COMPMIX 62C MHP 8/25 Removed unused variables 63C CHARACTER*256 FLAOL, FPUREZ 64C CHARACTER*256 FOPALE,FOPALE01,FOPALE06 ! FcondOpacP 65 66C 9/25/26 MHP via Claude 67C code changed to consolidate all I/O logical unit numbers, previously 68C scattered across COMMON/LUOUT/, COMMON/LUNUM/, COMMON/ALEX06/, 69C COMMON/ACDPTH/, COMMON/NEWOPAC/, COMMON/LOPAL95/, COMMON/ATMOS2/, 70C COMMON/ALATM03/, COMMON/OPALEOS/, and COMMON/SCV2/, into a single 71C COMMON/IOUNITS/ 72C 10/12/26 MHP via Claude 73C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 74C namelist variables that the physics routines read from 27 COMMON 75C blocks which also held derived or working members. Those blocks keep 76C only their non-namelist members; five of them -- NEWMX, VMULT, 77C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 78C Derivation-only namelist variables from the same blocks went to 79C COMMON/NMLBYPASS/ instead, not here. 80C 81C Declared identically in all 121 files that hold any of these members. 82C Generated from one member list, so every copy is byte-identical; the 83C linker's 'size of symbol changed' warning is the check that they are. 84 DIMENSION HCOMP(15,JSON),HD(JSON),HL(JSON),HP(JSON),HR(JSON), 85 * * HS(JSON),HT(JSON),LC(JSON),TRIT(3),TRIL(3),PS(3),TS(3),RS(3), 86 * * CFENV(9),TLUMX(8),OMEGA(JSON),FP(JSON),FT(JSON),ETA2(JSON),R0(JSON), 87 * * HJM(JSON),HI(JSON),ID(JSON),HS1(JSON) 88C llp 3/19/03 Add COMMON block /I2O/ for info directly transferred from 89C input to output model - starting with a code for th initial model 90C compostion (COMPMIX) 91C 10/13/26 MHP via Claude 92C code changed to replace COMMON/I2O/ with an explicit argument. COMPMIX 93C is CHARACTER, which does not belong in a COMMON block. GETMODEL2, 94C GETYREC7 and PUTMODEL2 already took it as an argument; only these three 95C files shared it through the block. 96 97C llp 3/19/03 Add required COMMON blocks such that header flags 98C ATM, EOS, HIK and LOK can be determined. 99C 10/13/26 MHP via Claude: removed the unused COMMON/ATMOS/ declaration 100C (no member referenced here; residue of the obsolete-option removals). 101C MHP 8/17 ADDED EXCEN, C_2 TO COMMON BLOCK FOR MATT ET AL. 2012 CENT. TERM 102C 10/13/26 MHP via Claude: removed the unused COMMON/CWIND/ declaration 103C (no member referenced here; residue of the obsolete-option removals). 104C COMMON/CWIND/WMAX,EXMD,EXW,EXTAU,EXR,EXM,CONSTFACTOR,STRUCTFACTOR,LJDOT0 105C 10/13/26 MHP via Claude: removed the unused COMMON/DEBHU/ declaration 106C (no member referenced here; residue of the obsolete-option removals). 107C 10/13/26 MHP via Claude: removed the unused COMMON/DISK/ declaration 108C (no member referenced here; residue of the obsolete-option removals). 109C 10/13/26 MHP via Claude: removed the unused COMMON/DPMIX/ declaration 110C (no member referenced here; residue of the obsolete-option removals). 111C 10/13/26 MHP via Claude: removed the unused COMMON/NEWOPAC/ declaration 112C (no member referenced here; residue of the obsolete-option removals). 113C 10/13/26 MHP via Claude: removed the unused COMMON/SCVEOS/ declaration 114C (no member referenced here; residue of the obsolete-option removals). 115 116C INCLUDE COMMON BLOCKS WITH REQUIRED PHYSICS AND OUTPUT FLAGS 117C 10/11/26 MHP via Claude 118C code changed to remove LSCRIB from COMMON/CCOUT/: it was declared 119C positionally in all 17 files sharing this block but never read or 120C written in executable code anywhere, and no namelist file sets it. 121C 10/11/26 MHP via Claude 122C BUG FIX: this file reads LMILNE below (gating CALL WRTMIL) but never 123C declared COMMON/CCOUT2/, so it was an uninitialized local, always 124C .FALSE. -- the Milne output call here was dead code. No namelist 125C file currently sets LMILNE, so this has no effect on any existing 126C run; it only restores the intended behavior if someone enables it. 127C 9/25/26 MHP via Claude 128C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/, 129C and COMMON/CONST3/ into a single COMMON/CONSTANTS/ 130 131 SAVE 132 133C physics flags: 134C Determine atmosphere flag, ATM 135 IF (KTTAU .EQ. 0) THEN 136 ATM='EDD ' 137 ELSEIF (KTTAU .EQ. 1) THEN 138 ATM='KS ' 139 ELSEIF (KTTAU .EQ. 2) THEN 140 ATM='HRA ' 141 ELSEIF (KTTAU .EQ. 3) THEN 142 ATM='KUR ' 143 ELSEIF (KTTAU .EQ. 4) THEN 144 ATM='ALL ' 145 ENDIF 146C Determine equation of state flag, EOS 147 EOS='SAHA ' 148 IF (LDH) EOS='SAH+DH' 149 IF (LSCV) THEN 150 EOS='SCV ' 151C 9/24/26 MHP via Claude 152C code changed to remove the OPAL 1995/2001 equations of state 153 IF (LDH) THEN 154 IF (LOPALE06) EOS='SCV+O6' 155 EOS='SCV+DH' 156 IF (LOPALE06) EOS='SCDHO6' 157 ENDIF 158 ELSE 159 IF (LOPALE06) THEN 160 EOS='OPAL06' 161 IF (LDH) EOS='OP6+DH' 162 ENDIF 163 ENDIF 164C Determine low temperature opacities flag, LOK 165 ALOK='NONE' 166C 9/23/26 MHP via Claude 167C code changed to remove the Kurucz 1990 opacity tables 168C 9/24/26 MHP via Claude 169C code changed to remove the Alexander 1995 opacity tables 170C Determine high temperature opacities flag, HIK 171 HIK='NONE' 172 IF (LOPAL95) HIK='OP95' 173C 9/23/26 MHP via Claude 174C code changed to remove the LAOL89 opacity tables 175 176 IF(ATM .EQ. ' ? ') THEN 177 WRITE(ISHORT,7) 178 7 7 FORMAT('*** YREC7 input file, flags, etc., have been ', 179 * * 'defaulted. ***') 180 ENDIF 181 182C 09/25 JvS: Add secondary format option that prints stitched interior and envelope 183C points. Output is either in the old format or new format, not both. 184 IF(LSTCH)THEN 185 IF(LPHHD)THEN 186 WRITE(ISTOR,1013) ! header key 187 ! write model physics header. Should only happen upon first model output. 188 ! write physics flags: 189 WRITE(ISTOR,29) JCORE,JENV,CMIXL,EOS,ATM,ALOK,HIK,LPUREZ, 190 & & COMPMIX,LEXCOM,LDIFY,LDIFZ,LSEMIC,LOVSTC,LOVSTE,LOVSTM, 191 & & LROT,LINSTB,LJDOT0,LDISK,TDISK,PDISK,WMAX,LSTORE,LSTATM,LSTENV, 192 & & LSTMOD,LSTPHYS,LSTROT 193 29 29 FORMAT('#',2I8,F16.10,1X,A6,1X,3(A4,1X),L1,1X,A4,1X,11(L1,1X), 194 & & 3(1PE18.10),1X,6(L1,1X)) 195 WRITE(ISTOR,1014) ! profile header 196 LPHHD = .FALSE. ! Turn off the physics header fo the rest of the run. 197 ENDIF 198 CALL STITCH(HCOMP,HR,HP,HD,HS,HT,HL,HS1,OMEGA,ETA2, 199 * * HI,R0,HJM,FP,FT,TEFFL,HSTOT,BL,M,LC,MODEL,CMIXL) 200 1013 1013 FORMAT('# JCORE JENV CMIXL EOS ATM ALOK HIK LPUREZ COMPMIX', 201 1 1 ' LEXCOM LDIFY LDIFZ LSEMIC LOVSTC LOVSTE LOVSTM', 202 1 1 ' LROT LINSTB LJDOT0 LDISK TDISK PDISK WMAX LSTORE', 203 1 1 ' LSTATM LSTENV LSTMOD LSTPHYS LSTROT') 204 1014 1014 FORMAT( 205 1 1'MODEL SHELL MASS RADIUS LUMINOSITY PRESSURE TEMPERATURE DENSITY OMEGA ', 206 1 1'CONVECTIVE INTERIOR_PT ENV_PT ATM_PT H1 He4 METALS He3 C12 C13 N14 N15 O16 ', 207 1 1'O17 O18 H2 Li6 Li7 Be9 OPACITY GRAV DELR DEL DELAD V_CONV GAM1 HII HEII HEIII ', 208 1 1'BETA ETA PPI PPII PPIII CNO TRIPLE_ALPHA E_NUC E_NEU E_GRAV CP DLNRHODLNT A ', 209 1 1'RP/RE FP FT J/M MOMENT DEL_KE V_ES V_GSF V_SS VTOT ') 210 211 ELSE 212C write header records 213 IF(DAGE .lt. 1d3) THEN 214 WRITE(ISTOR,10) 'MOD2 ',MODEL,M,SMASS,TEFFL,BL,HSTOT,DAGE, 215 * * DDAGE,HS(1),HS(M) 216 10 10 FORMAT(A5,2I8,5F16.11,1PE18.10,0P2F16.12) 217 ELSE IF (DAGE .LT. 1D4) THEN 218 WRITE(ISTOR,11) 'MOD2 ',MODEL,M,SMASS,TEFFL,BL,HSTOT,DAGE, 219 * * DDAGE,HS(1),HS(M) 220 11 11 FORMAT(A5,2I8,4F16.12,F16.10,1PE18.10,0P2F16.12) 221 ELSE IF (DAGE .LT. 1D5) THEN 222 WRITE(ISTOR,12) 'MOD2 ',MODEL,M,SMASS,TEFFL,BL,HSTOT,DAGE, 223 * * DDAGE,HS(1),HS(M) 224 12 12 FORMAT(A5,2I8,4F16.12,F16.9,1PE18.10,0P2F16.12) 225 ELSE 226 WRITE(ISTOR,13) 'MOD2 ',MODEL,M,SMASS,TEFFL,BL,HSTOT,DAGE, 227 * * DDAGE,HS(1),HS(M) 228 13 13 FORMAT(A5,2I8,4F16.12,F16.8,1PE18.10,0P2F16.12) 229 ENDIF 230 231C write physics flags: 232 WRITE(ISTOR,30) JCORE,JENV,CMIXL,EOS,ATM,ALOK,HIK,LPUREZ, 233 & & COMPMIX,LEXCOM,LDIFY,LDIFZ,LSEMIC,LOVSTC,LOVSTE,LOVSTM, 234 & & LROT,LINSTB,LJDOT0,LDISK,TDISK,PDISK,WMAX,LSTORE,LSTATM,LSTENV, 235 & & LSTMOD,LSTPHYS,LSTROT 236 30 30 FORMAT(2I8,F16.10,1X,A6,1X,3(A4,1X),L1,1X,A4,1X,11(L1,1X), 237 & & 3(1PE18.10),1X,6(L1,1X)) 238 239C write luminosities and output flags 240C If TLUMX are in solar units, convert to ergs. Decide by 241C comparing to 10**20, if smaller, multiply by CLSUN 242 243 CCCMAX = DMAX1(TLUMX(1),TLUMX(2),TLUMX(3),TLUMX(4),TLUMX(5), 244 * * DABS(TLUMX(6)),TLUMX(7)) 245 IF(CCCMAX.LE.1.0D20) THEN 246 DO J = 1,7 247 TLUMX(J) = TLUMX(J) * CLSUN 248 ENDDO 249 ENDIF 250 WRITE(ISTOR,40) (TLUMX(J),J=1,7) 251 40 40 FORMAT('TLUMX',5X,1P7E17.9) 252 253C write ENVELOPE DATA 254 DO I = 1,3 255 I0 = I 256 IF(FTRI.LT.0D0) I0 = -I 257 WRITE(ISTOR,50)I0,TRIT(I),TRIL(I),PS(I),TS(I),RS(I), 258 & & (CFENV(3*I-3+J),J=1,3) 259 ENDDO 260 50 50 FORMAT('ENV',I2,5F16.12,1P3E20.12) 261 262 CALL PINDEX(JXBEG,JXEND,LSHELL,M,ID,IDM) 263 264 IF(LSTMOD)THEN 265C write column headings for all per shell information 266 WRITE(ISTOR,55) 267 55 55 FORMAT(/, 268 1 1' SHELL MASS RADIUS LUMINOSITY ', 269 1 1'PRESSURE TEMPERATURE DENSITY OMEGA ', 270 1 1' C H1 He4 METALS He3 C12 ', 271 1 1' C13 N14 N15 O16 ', 272 1 1' O17 O18 H2 Li6 Li7', 273 1 1' Be9 OPAC GRAV DELR DEL ', 274 1 1' DELA V_CONV GAM1 HII HEII HEIII BETA ', 275 1 1'ETA PPI PPII PPIII CNO 3HE ', 276 1 1'E_NUC E_NEU E_GRAV A RP/RE FP', 277 1 1' FT J/M MOMENT DEL_KE V_ES ', 278 1 1' V_GSF V_SS VTOT ',/) 279 280C write out the requested information. 281 CG=DEXP(CLN*CGL) 282 DO II = 1,IDM 283 I = ID(II) 284C write out the basic info 285C 9/28/26 MHP via Claude 286C code changed to add MODEL to this WRITE, matching FORMAT 62 (which 287C already expected it) and fixing a crash ("Expected INTEGER for item 288C 2 ... got REAL") present on main itself: FORMAT 62 was widened there 289C but this WRITE was never updated to match. 290C 10/8/26 MHP via Claude 291C code changed to remove the interior/envelope/atmosphere region-flag 292C triple added above: after consultation with other developers, this 293C doesn't fit the .store header (FORMAT 55 documents no such columns) 294C and isn't desired here. This file only ever writes interior points, 295C so the flags were always the same hard-coded, uninformative 296C .TRUE.,.FALSE.,.FALSE. -- unlike stitch.f, which genuinely needs the 297C triple since it writes interior, envelope, and atmosphere rows with 298C the same format. MODEL is kept; FORMAT 62 below now supplies only 299C LC(I)'s own logical, as it did before the region-flag triple was 300C ever added. 301 WRITE(ISTOR,62,ADVANCE='no') MODEL,I,HS(I),HR(I),HL(I), 302 * * HP(I),HT(I),HD(I),OMEGA(I),LC(I), 303 * * (HCOMP(J,I),J=1,15) 304C write out additional physics if desired 305 IF(LSTPHYS)THEN 306 SG = DEXP(CLN*(CGL - 2.0D0*HR(I)))*HS1(I) 307 WRITE(ISTOR,63,ADVANCE='no') SO(I),SG,SDEL(1,I),SDEL(2,I), 308 * * SDEL(3,I),SVEL(I),GAM1(I),0.0,0.0,0.0,SBETA(I),SETA(I), 309 * * (SEG(K,I),K=1,5),SESUM(I),SEG(6,I),SEG(7,I) 310c WRITE(ISTOR,63,ADVANCE='no') SO(I),SG,SDEL(1,I),SDEL(2,I), 311c * SDEL(3,I),SVEL(I),GAM1(I),SFXION(1,I),SFXION(2,I),SFXION(3,I), 312c * SBETA(I),SETA(I),(SEG(K,I),K=1,5),SESUM(I),SEG(6,I),SEG(7,I) 313 ELSE 314 WRITE(ISTOR,63,ADVANCE='no') 0.0,0.0,0.0,0.0,0.0,0.0,0.0, 315 * * 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0 316 ENDIF 317C write out additional rotation info if desired 318 IF(LSTROT.AND.LROT)THEN 319 FM = DEXP(CLN*HS(I)) 320 DUMA = CC13*OMEGA(I)**2/(CG*FM)*5.D0/(2.D0+ETA2(I)) 321 A = DUMA * R0(I)**3 322 RPOLEQ = (1.0D0 - A)/(1.0D0 + 0.5D0*A) 323 VTOT = VES(I)+VGSF(I)+VSS(I) 324 WRITE(ISTOR,64) A,RPOLEQ,FP(I),FT(I),HJM(I),HI(I),DEROT(I), 325 * * VES(I),VGSF(I),VSS(I),VTOT 326 ELSE 327 WRITE(ISTOR,64) 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0 328 ENDIF 329 ENDDO 330 62 62 FORMAT(I6,1X,I6,0P2F18.14,1PE24.16,0PF18.14,1X,F18.14,1X,F18.14,1PE24.16, 331 & & 1X,L1,3(0PF12.9),12(0PE16.8),2X) 332 333 63 63 FORMAT(1PE10.4,1PE11.3,E12.4,E12.4,E12.4,1PE12.4,0PF9.5,1X,F9.5,1X,F9.5, 334 & & 1X,F9.5,1X,F9.5,1X,F9.5,E12.4,E12.4,E12.4,E12.4,E12.4,E13.5,E13.5,E13.5) 335 336 64 64 FORMAT(E14.6,E14.6,E14.6,E14.6,E14.6,E14.6,E14.6,E11.3,E11.3,E11.3,E11.3, 337 & & E11.3) 338 339 ENDIF 340C now call wrtmod, with the goal of outputting the envelope and atmosphere, or 341C if required by LPULSE. 342 IF(LSTATM.OR.LSTENV)THEN 343 IF(LMILNE) CALL WRTMIL(HCOMP,HD,HL,HP,HR,HS1,M,MODEL) 344C CALL WRTMOD(M,LSHELL,JXBEG,JXEND,JCORE,JENV,HCOMP,HS1,HD,HL, 345C * HP,HR,HT,LC,MODEL,BL,TEFFL,OMEGA,FP,FT,ETA2,R0,HJM,HI,HS, 346C * DAGE) ! KC 2025-05-31 347 CALL WRTMOD(M,JENV,HCOMP,HS1,HD,HL, 348 * * HP,HR,HT,MODEL,BL,TEFFL,FP,FT,HS,DAGE,CMIXL,XMSOL) 349 ENDIF 350 351 WRITE(ISTOR,65) 352 65 65 FORMAT(/,/) 353 354 ENDIF !closes logic block for new vs old store file format 355 356C G Somers END 357CCCCCCCCCCCCCCCCCCCCCCCCC 358 359 RETURN 360 END 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385

1"""Python version of putstore.f (YREC 7be4133): writes a verbose model to 2the .store file.""" 3 4import math 5 6import numpy as np 7 8from pyrec import ccout_mod 9from pyrec import constants_mod 10from pyrec import iounits_mod 11from pyrec import nmlpar_mod 12from pyrec import roten_mod 13from pyrec import scrtch_mod 14from pyrec import sound_mod 15from pyrec import temp2_mod 16from pyrec.pindex import pindex 17from pyrec.stitch import stitch 18from pyrec.wrtmil import wrtmil 19from pyrec.wrtmod import wrtmod 20 21# Locals kept between calls (SAVE). ATM keeps its last value when KTTAU is 22# not 0-4; JXBEG, JXEND, LSHELL are never set (see the page, bug 3) 23# (Fortran lines 51–52, 131). 24atm = None 25jxbeg = 0 26jxend = 0 27lshell = False 28 29_ATM_LABELS = {0: 'EDD ', 1: 'KS ', 2: 'HRA ', 3: 'KUR ', 4: 'ALL '} 30 31# Output formats are approximate: Fortran's 0PE (0.ddd) is written as 32# Python's E (d.ddd). FORMATs 1013 and 1014 (Fortran lines 200–209). 33_FMT1013 = ('# JCORE JENV CMIXL EOS ATM ALOK HIK LPUREZ COMPMIX' 34 ' LEXCOM LDIFY LDIFZ LSEMIC LOVSTC LOVSTE LOVSTM' 35 ' LROT LINSTB LJDOT0 LDISK TDISK PDISK WMAX LSTORE' 36 ' LSTATM LSTENV LSTMOD LSTPHYS LSTROT') 37_FMT1014 = ( 38 'MODEL SHELL MASS RADIUS LUMINOSITY PRESSURE TEMPERATURE DENSITY OMEGA ' 39 'CONVECTIVE INTERIOR_PT ENV_PT ATM_PT H1 He4 METALS He3 C12 C13 N14 N15 O16 ' 40 'O17 O18 H2 Li6 Li7 Be9 OPACITY GRAV DELR DEL DELAD V_CONV GAM1 HII HEII HEIII ' 41 'BETA ETA PPI PPII PPIII CNO TRIPLE_ALPHA E_NUC E_NEU E_GRAV CP DLNRHODLNT A ' 42 'RP/RE FP FT J/M MOMENT DEL_KE V_ES V_GSF V_SS VTOT ') 43# FORMAT 55 (Fortran lines 267–278). Suspected bug (see the page): '3HE' 44# names the triple-alpha column. 45_FMT55 = ( 46 '\n' 47 ' SHELL MASS RADIUS LUMINOSITY ' 48 'PRESSURE TEMPERATURE DENSITY OMEGA ' 49 ' C H1 He4 METALS He3 C12 ' 50 ' C13 N14 N15 O16 ' 51 ' O17 O18 H2 Li6 Li7' 52 ' Be9 OPAC GRAV DELR DEL ' 53 ' DELA V_CONV GAM1 HII HEII HEIII BETA ' 54 'ETA PPI PPII PPIII CNO 3HE ' 55 'E_NUC E_NEU E_GRAV A RP/RE FP' 56 ' FT J/M MOMENT DEL_KE V_ES ' 57 ' V_GSF V_SS VTOT \n') 58 59 60def _l(flag): 61 """Fortran's L1.""" 62 63 return 'T' if flag else 'F' 64 65 66def _flags_line(jcore, jenv, cmixl, eos, alok, hik, compmix): 67 """FORMATs 29 and 30 (29 adds a leading '#'): the physics flags.""" 68 69 n = nmlpar_mod 70 return (f"{jcore:8d}{jenv:8d}{cmixl:16.10f} {eos:6s} {atm:4s} " 71 f"{alok:4s} {hik:4s} {_l(n.lpurez)} {compmix:4s} " 72 + "".join(_l(v) + " " for v in ( 73 n.lexcom, n.ldify, n.ldifz, n.lsemic, n.lovstc, n.lovste, 74 n.lovstm, n.lrot, n.linstb, n.ljdot0, n.ldisk)) 75 + f"{n.tdisk:18.10E}{n.pdisk:18.10E}{n.wmax:18.10E} " 76 + "".join(_l(v) + " " for v in ( 77 n.lstore, n.lstatm, n.lstenv, n.lstmod, n.lstphys, 78 n.lstrot))) 79 80 81def putstore(hcomp: np.ndarray, hd: np.ndarray, hl: np.ndarray, hp: np.ndarray, hr: np.ndarray, hs: np.ndarray, 82 ht: np.ndarray, lc: np.ndarray, trit: np.ndarray, tril: np.ndarray, ps: np.ndarray, ts: np.ndarray, 83 rs: np.ndarray, cfenv: np.ndarray, ftri: float, tlumx: np.ndarray, jcore: int, jenv: int, model: int, 84 m: int, smass: float, teffl: float, bl: float, hstot: float, dage: float, ddage: float, omega: np.ndarray, 85 hs1: np.ndarray, eta2: np.ndarray, r0: np.ndarray, fp: np.ndarray, ft: np.ndarray, hjm: np.ndarray, 86 hi: np.ndarray, cmixl: float, xmsol: float, compmix: str) -> None: 87 """Write the model to the .store file (unit ISTOR): either the 88 "stitched" format through stitch (LSTCH), or the older one: a MOD2 89 header, the flags, TLUMX, the envelope lines and, with LSTMOD, a line 90 per selected shell with optional physics (LSTPHYS) and rotation 91 (LSTROT) columns; then, with LSTATM or LSTENV, the envelope and 92 atmosphere through wrtmil (LMILNE) and wrtmod. 93 94 Translates: putstore.f#PUTSTORE (7be4133) 95 96 Arguments: the model's arrays (hcomp ... lc, omega, hs1, eta2, r0, fp, 97 ft, hjm, hi; index 0 unused), the envelope triangle (trit, tril, ps, 98 ts, rs, cfenv, ftri), the luminosities tlumx, and the scalars; the 99 mixing length cmixl, xmsol (passed to wrtmod) and compmix, the initial 100 mixture's code. Arrays changed in place: tlumx (converted to erg/s if 101 given in solar units). Nothing is returned. 102 103 """ 104 105 global atm 106 n = nmlpar_mod 107 c = constants_mod 108 istor = iounits_mod.istor 109 sc = scrtch_mod 110 t2 = temp2_mod 111 112 # Physics flags. Atmosphere flag ATM (Fortran lines 133–145). 113 # Suspected bug (see the page): KTTAU = 5 has no label. 114 if n.kttau in _ATM_LABELS: 115 atm = _ATM_LABELS[n.kttau] 116 117 # Equation of state flag EOS (Fortran lines 146–163). 118 eos = 'SAHA ' 119 if n.ldh: 120 eos = 'SAH+DH' 121 if n.lscv: 122 eos = 'SCV ' 123 if n.ldh: 124 # Suspected bug (see the page): overwritten on the next line. 125 if n.lopale06: 126 eos = 'SCV+O6' 127 eos = 'SCV+DH' 128 if n.lopale06: 129 eos = 'SCDHO6' 130 else: 131 if n.lopale06: 132 eos = 'OPAL06' 133 if n.ldh: 134 eos = 'OP6+DH' 135 136 # Low-temperature opacity flag, always NONE (Fortran lines 164–169). 137 # Question (see the page): no label for LALEX06. 138 alok = 'NONE' 139 140 # High-temperature opacity flag (Fortran lines 170–174). 141 hik = 'NONE' 142 if n.lopal95: 143 hik = 'OP95' 144 145 # Never true: no branch sets ' ? ' (Fortran lines 176–180). 146 if atm == ' ? ': 147 iounits_mod.write(iounits_mod.ishort, 148 '*** YREC7 input file, flags, etc., have been ' 149 'defaulted. ***') 150 151 # The stitched format, or the old one, not both (JvS 09/25) (Fortran 152 # lines 182–210). 153 if n.lstch: 154 if ccout_mod.lphhd: 155 iounits_mod.write(istor, _FMT1013) 156 iounits_mod.write(istor, '#' + _flags_line( 157 jcore, jenv, cmixl, eos, alok, hik, compmix)) 158 iounits_mod.write(istor, _FMT1014) 159 ccout_mod.lphhd = False # the physics header only once 160 stitch(hcomp, hr, hp, hd, hs, ht, hl, hs1, omega, eta2, hi, r0, hjm, 161 fp, ft, teffl, hstot, bl, m, lc, model, cmixl) 162 return 163 164 # Header record (FORMATs 10-13): the decimals depend on DAGE's size 165 # (Fortran lines 212–229). 166 if dage < 1.0e3: 167 head = "".join(f"{v:16.11f}" for v in (smass, teffl, bl, hstot, dage)) 168 else: 169 if dage < 1.0e4: 170 nd = 10 171 elif dage < 1.0e5: 172 nd = 9 173 else: 174 nd = 8 175 head = ("".join(f"{v:16.12f}" for v in (smass, teffl, bl, hstot)) 176 + f"{dage:16.{nd}f}") 177 iounits_mod.write(istor, f"MOD2 {model:8d}{m:8d}" + head 178 + f"{ddage:18.10E}{hs[1]:16.12f}{hs[m]:16.12f}") 179 180 # Physics flags (FORMAT 30) (Fortran lines 231–237). 181 iounits_mod.write(istor, _flags_line(jcore, jenv, cmixl, eos, alok, 182 hik, compmix)) 183 184 # Luminosities: if TLUMX are in solar units, convert to erg/s (decided 185 # by comparing with 1e20), in place (Fortran lines 239–251). 186 # Question (see the page): duplicated code. 187 cccmax = max(tlumx[1], tlumx[2], tlumx[3], tlumx[4], tlumx[5], 188 abs(tlumx[6]), tlumx[7]) 189 if cccmax <= 1.0e20: 190 for j in range(1, 8): 191 tlumx[j] = tlumx[j] * c.clsun 192 iounits_mod.write(istor, 'TLUMX ' 193 + "".join(f"{tlumx[j]:17.9E}" for j in range(1, 8))) 194 195 # Envelope data (FORMAT 50) (Fortran lines 253–260). 196 for i in range(1, 4): 197 i0 = i 198 if ftri < 0.0: 199 i0 = -i 200 iounits_mod.write( 201 istor, f"ENV{i0:2d}{trit[i]:16.12f}{tril[i]:16.12f}" 202 f"{ps[i]:16.12f}{ts[i]:16.12f}{rs[i]:16.12f}" 203 + "".join(f"{cfenv[3 * i - 3 + j]:20.12E}" for j in range(1, 4))) 204 205 # Which shells to print (Fortran line 262). Suspected bug (see the 206 # page): jxbeg, jxend and lshell are never set. 207 id_ = np.zeros(len(hs), dtype=int) 208 idm = pindex(jxbeg, jxend, lshell, m, id_, 0) 209 210 # One line per selected shell (Fortran lines 264–339). 211 if n.lstmod: 212 iounits_mod.write(istor, _FMT55) 213 cg = math.exp(c.cln * c.cgl) 214 for ii in range(1, idm + 1): 215 i = id_[ii] 216 217 # The basic columns, led by MODEL (FORMAT 62) (Fortran lines 218 # 284–303). 219 line = (f"{model:6d} {i:6d}{hs[i]:18.14f}{hr[i]:18.14f}" 220 f"{hl[i]:24.16E}{hp[i]:18.14f} {ht[i]:18.14f} " 221 f"{hd[i]:18.14f}{omega[i]:24.16E} {_l(lc[i])}" 222 + "".join(f"{hcomp[j, i]:12.9f}" for j in range(1, 4)) 223 + "".join(f"{hcomp[j, i]:16.8E}" for j in range(4, 16)) 224 + " ") 225 226 # Additional physics (FORMAT 63) (Fortran lines 304–316). 227 if n.lstphys: 228 sg = math.exp(c.cln * (c.cgl - 2.0 * hr[i])) * hs1[i] 229 # Question (see the page): the ionization fractions are 230 # written as zeros. 231 vals = [sc.so[i], sg, sc.sdel[1, i], sc.sdel[2, i], 232 sc.sdel[3, i], sc.svel[i], sound_mod.gam1[i], 233 0.0, 0.0, 0.0, sc.sbeta[i], sc.seta[i]] 234 vals += [sc.seg[k, i] for k in range(1, 6)] 235 vals += [sc.sesum[i], sc.seg[6, i], sc.seg[7, i]] 236 else: 237 vals = [0.0] * 20 238 line += (f"{vals[0]:10.4E}{vals[1]:11.3E}{vals[2]:12.4E}" 239 f"{vals[3]:12.4E}{vals[4]:12.4E}{vals[5]:12.4E}" 240 f"{vals[6]:9.5f}" 241 + "".join(f" {v:9.5f}" for v in vals[7:12]) 242 + "".join(f"{v:12.4E}" for v in vals[12:17]) 243 + "".join(f"{v:13.5E}" for v in vals[17:20])) 244 245 # Additional rotation info (FORMAT 64) (Fortran lines 317–328). 246 if n.lstrot and n.lrot: 247 fm = math.exp(c.cln * hs[i]) 248 duma = c.cc13 * omega[i]**2 / (cg * fm) * 5.0 / (2.0 + eta2[i]) 249 a = duma * r0[i]**3 250 rpoleq = (1.0 - a) / (1.0 + 0.5 * a) 251 vtot = t2.ves[i] + t2.vgsf[i] + t2.vss[i] 252 rot = [a, rpoleq, fp[i], ft[i], hjm[i], hi[i], 253 roten_mod.derot[i], t2.ves[i], t2.vgsf[i], t2.vss[i], 254 vtot] 255 else: 256 rot = [0.0] * 11 257 line += ("".join(f"{v:14.6E}" for v in rot[:7]) 258 + "".join(f"{v:11.3E}" for v in rot[7:])) 259 iounits_mod.write(istor, line) 260 261 # The envelope and atmosphere (Fortran lines 340–349). 262 if n.lstatm or n.lstenv: 263 if n.lmilne: 264 wrtmil(hcomp, hd, hl, hp, hr, hs1, m, model) 265 wrtmod(m, jenv, hcomp, hs1, hd, hl, hp, hr, ht, model, bl, teffl, 266 fp, ft, hs, dage, cmixl, xmsol) 267 268 # FORMAT 65, (/,/) (Fortran lines 351–352). 269 iounits_mod.write(istor, '\n\n')

putstore.f

Purpose

How it fits

The method

From YREC's documentation

Suspected bugs, checked by running the code

1. (never occurred) The SCV + OPAL 2006 label is lost.

2. (wrong output) Checked by KC (Nov 10, 2026) Atmosphere 5 has no label.

3. (never occurred) Checked by KC (Nov 10, 2026) PINDEX gets arguments that are never set.

4. (wrong output) Checked by KC (Nov 10, 2026) The old header calls the triple-α column "3HE."

Definitions

subroutinePUTSTORE

Change log

3. (never occurred) Checked by KC (Nov 10, 2026) `PINDEX` gets arguments that are never set.