wrtmil.f

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

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

Variables on this page (14)

c4pi, cgl, cln, clsun, cmsun, imilne, json, nprtpt, sdel, seg, sesum, stotal, xenv, zenv

Purpose

WRTMIL writes one model's Milne (homology) invariants $U$ and $V$, with $W$ and the polytropic index $n + 1$, to the Milne output file, every NPRTPT shells and the last one.

How it fits

Inputs: the model (HCOMP, $\log\rho$, $L/L_\odot$, $\log P$, $\log r$, $m$ in grams HS1, M, MODEL), and through modules the energy rates SESUM and SEG and the actual gradient SDEL(2,·) (SCRTCH_MOD), the surface composition and total mass (COMP_MOD), NPRTPT (NMLPAR_MOD), the constants (CONSTANTS_MOD) and the unit IMILNE (IOUNITS_MOD).
Who calls it: PUTSTORE and WRTOUT, when LMILNE is set and other conditions hold: PUTSTORE only when LSTATM or LSTENV is set, and WRTOUT only in its pulsation branch (LPULSE on a model that is not stored). So LMILNE alone writes nothing.

Question 1 (putstore.f, lines 342–343 ↗; wrtout.f, lines 683–684 ↗). LMILNE alone writes nothing. The calls also need LSTATM or LSTENV (PUTSTORE), or LPULSE on a model that is not stored (WRTOUT). Is the Milne output meant to depend on those switches?
Output: the unit IMILNE (21, fixed in PARMIN, parmin.f, line 674 ↗), opened by PARMIN in append mode (parmin.f, lines 956–959 ↗).

The method

At each printed shell: $$U = \frac{4\pi\rho r^3}{m},\qquad V = \frac{G m \rho}{P r},\qquad W = U\,\frac{m\,(\epsilon_\mathrm{nuc} + \epsilon_7)}{L},\qquad n + 1 = \frac{1}{\nabla},$$ where $\epsilon_\mathrm{nuc}$ is SESUM and $\epsilon_7$ is SEG(7,·) (the gravitational energy rate, by the column names stitch.f gives), and $\nabla$ = SDEL(2,·), the actual gradient.

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 file is closed after each model.

(line 109 ↗; parmin.f, lines 956–959 ↗) WRTMIL ends with CLOSE(IMILNE), but the file is opened only once, in PARMIN (MHP 8/25, line 36 ↗). The next call writes to unit 21 unopened, which gfortran, like most compilers, sends to a default file (fort.21), so the Milne file holds only the first model written. Only with LMILNE (off by default).

A possible fix:

wrtmil.f

@@ lines 105–111 @@
        NP1 = 1.0D0/SDEL(2,M)
        WRITE(IMILNE,10)M,HS1(M),R,P,D,HCOMP(1,M),HL(M),
      *                   U,V,W,NP1
+C 11/09/26 KC via Claude
+C code changed to keep IMILNE open, since PARMIN opens it only once
+C for the run.
       ENDIF
-      CLOSE(IMILNE)
       RETURN
       END

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): LMILNE is off by default (nmlpar_mod.f, line 58 ↗) and set in none of YREC's cases or documents. Not run: LMILNE on, with two or more models written. Intent: The comment at line 36 ↗ says the file is now opened in PARMIN, so the CLOSE looks left over from when WRTMIL opened the file itself.

Definitions

subroutineWRTMIL

lines 10–111 · GitHub

Milne invariants $U$, $V$, with $W$ and $n+1 = 1/\nabla$, for selected shells of one model.

Uses modules COMP_MOD, CONSTANTS_MOD, DIMENSIONS_MOD, IOUNITS_MOD, NMLPAR_MOD, SCRTCH_MOD

Called by PUTSTORE, WRTOUT

Change log

The changes since rp_1

Changed since rp_1. Declarations only: COMMON blocks to modules (/SCRTCH/, /COMP/, /CCOUT1/ to NMLPAR_MOD, /LUOUT/ to IOUNITS_MOD, the constants to CONSTANTS_MOD), IMPLICIT NONE.

1 2 3 4 5 6C 7C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 8C WRTMIL 9C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 10 SUBROUTINE WRTMIL(HCOMP,HD,HL,HP,HR,HS1,M,MODEL) 11C 10/15/26 MHP via Claude 12C code changed to take JSON from DIMENSIONS_MOD instead of declaring 13C PARAMETER(JSON=5000) here. 14 USE DIMENSIONS_MOD 15C 10/14/26 MHP via Claude 16C code changed to convert COMMON/SCRTCH/ into SCRTCH_MOD. 17 USE SCRTCH_MOD 18C 10/14/26 MHP via Claude 19C code changed to convert COMMON/COMP/ into COMP_MOD. 20 USE COMP_MOD 21C 10/14/26 MHP via Claude 22C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD. 23 USE NMLPAR_MOD 24C 10/14/26 MHP via Claude 25C code changed to convert COMMON/CONSTANTS/ into CONSTANTS_MOD. 26 USE CONSTANTS_MOD 27C 10/14/26 MHP via Claude 28C code changed to convert COMMON/IOUNITS/ into IOUNITS_MOD. 29 USE IOUNITS_MOD 30C 10/20/26 MHP via Claude 31C code changed to replace implicit typing with IMPLICIT NONE and 32C explicit declarations. 33 IMPLICIT NONE 34 REAL*8 D,HCOMP,HD,HL,HP,HR,HS1,P,R,SMTOT,U,V,W 35 INTEGER I,IBEG,IEND,M,MODEL 36C MHP 8/25 Milne output file is now opened in parmin,unused variables removed 37C 9/25/26 MHP via Claude 38C code changed to consolidate all I/O logical unit numbers, previously 39C scattered across COMMON/LUOUT/, COMMON/LUNUM/, COMMON/ALEX06/, 40C COMMON/ACDPTH/, COMMON/NEWOPAC/, COMMON/LOPAL95/, COMMON/ATMOS2/, 41C COMMON/ALATM03/, COMMON/OPALEOS/, and COMMON/SCV2/, into a single 42C COMMON/IOUNITS/ 43C 10/12/26 MHP via Claude 44C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 45C namelist variables that the physics routines read from 27 COMMON 46C blocks which also held derived or working members. Those blocks keep 47C only their non-namelist members; five of them -- NEWMX, VMULT, 48C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 49C Derivation-only namelist variables from the same blocks went to 50C COMMON/NMLBYPASS/ instead, not here. 51C 52C Declared identically in all 121 files that hold any of these members. 53C Generated from one member list, so every copy is byte-identical; the 54C linker's 'size of symbol changed' warning is the check that they are. 55C 9/25/26 MHP via Claude 56C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/, 57C and COMMON/CONST3/ into a single COMMON/CONSTANTS/ 58 DIMENSION HCOMP(15,JSON),HD(JSON),HL(JSON),HP(JSON),HR(JSON), 59 * * HS1(JSON) 60 REAL*8 NP1 61 SAVE 62 63C FIND THE MILNE INVARIANTS U AND V, ALONG WITH THE INDEX N+1. 64C WRITE THEM OUT TO LOGICAL UNIT IMILNE. 65C HEADER 66 SMTOT = DEXP(CLN*STOTAL)/CMSUN 67 WRITE(IMILNE,5) MODEL,XENV,ZENV,SMTOT 68 5 5 FORMAT(10X,'MODEL',I5,' XENV =',1PD10.3,' ZENV =',D10.3, 69 * * ' MASS(SOLAR UNITS) =',D10.3) 70C CALCULATE FOR FIRST POINT(ALWAYS DONE) 71 D = DEXP(CLN*HD(1)) 72 P = DEXP(CLN*HP(1)) 73 R = DEXP(CLN*HR(1)) 74 U = C4PI*D*R**3/HS1(1) 75 V = DEXP(CLN*CGL)*HS1(1)*D/(P*R) 76 W = U*HS1(1)*(SESUM(1)+SEG(7,1))/(HL(1)*CLSUN) 77 NP1 = 1.0D0/SDEL(2,1) 78 WRITE(IMILNE,10)1,HS1(1),R,P,D,HCOMP(1,1),HL(1),U,V,W,NP1 79 10 10 FORMAT(1X,I4,10(1PE11.3)) 80C PRINT OUT EVERY NPRTPT POINTS;LAST POINT ALWAYS PRINTED. 81 IEND = 1 82 IF(NPRTPT.LE.M) THEN 83 IBEG = MAX(2,NPRTPT) 84 IEND = M - MOD(M,NPRTPT) 85 DO 20 I = IBEG,IEND,NPRTPT 86 D = DEXP(CLN*HD(I)) 87 P = DEXP(CLN*HP(I)) 88 R = DEXP(CLN*HR(I)) 89 U = C4PI*D*R**3/HS1(I) 90 V = DEXP(CLN*CGL)*HS1(I)*D/(P*R) 91 W = U*HS1(I)*(SESUM(I)+SEG(7,I))/(HL(I)*CLSUN) 92 NP1 = 1.0D0/SDEL(2,I) 93 WRITE(IMILNE,10)I,HS1(I),R,P,D,HCOMP(1,I),HL(I), 94 * * U,V,W,NP1 95 20 20 CONTINUE 96 ENDIF 97 IF(IEND.LT.M) THEN 98C PRINT OUT LAST POINT IF NPRTPT DOESNT DIVIDE EVENLY INTO M. 99 D = DEXP(CLN*HD(M)) 100 P = DEXP(CLN*HP(M)) 101 R = DEXP(CLN*HR(M)) 102 U = C4PI*D*R**3/HS1(M) 103 V = DEXP(CLN*CGL)*HS1(M)*D/(P*R) 104 W = U*HS1(M)*(SESUM(M)+SEG(7,M))/(HL(M)*CLSUN) 105 NP1 = 1.0D0/SDEL(2,M) 106 WRITE(IMILNE,10)M,HS1(M),R,P,D,HCOMP(1,M),HL(M), 107 * * U,V,W,NP1 108 ENDIF 109 CLOSE(IMILNE) 110 RETURN 111 END

1"""Python version of wrtmil.f (YREC 7be4133): the Milne (homology) 2invariants of one model, to the Milne output file.""" 3 4import math 5 6import numpy as np 7 8from pyrec import comp_mod 9from pyrec import constants_mod 10from pyrec import iounits_mod 11from pyrec import nmlpar_mod 12from pyrec import scrtch_mod 13 14 15def _row(i: int, hcomp: np.ndarray, hd: np.ndarray, hl: np.ndarray, hp: np.ndarray, hr: np.ndarray, 16 hs1: np.ndarray) -> str: 17 """One shell's line (FORMAT 10): U, V, W and n + 1 = 1/nabla (the 18 Fortran repeats these lines at each shell: lines 71–79, 86–94, and 19 99–107). 20 21 """ 22 23 c = constants_mod 24 s = scrtch_mod 25 d = math.exp(c.cln * hd[i]) 26 p = math.exp(c.cln * hp[i]) 27 r = math.exp(c.cln * hr[i]) 28 u = c.c4pi * d * r**3 / hs1[i] 29 v = math.exp(c.cln * c.cgl) * hs1[i] * d / (p * r) 30 # Question (see the page): W includes the gravitational term seg[7]. 31 w = u * hs1[i] * (s.sesum[i] + s.seg[7, i]) / (hl[i] * c.clsun) 32 np1 = 1.0 / s.sdel[2, i] 33 vals = (hs1[i], r, p, d, hcomp[1, i], hl[i], u, v, w, np1) 34 return f' {i:4d}' + ''.join(f'{x:11.3E}' for x in vals) 35 36 37def _d(x: float) -> str: 38 """Fortran's 1PD10.3: like 1PE10.3, with a D exponent. 39 40 """ 41 42 return f'{x:10.3E}'.replace('E', 'D') 43 44 45def wrtmil(hcomp: np.ndarray, hd: np.ndarray, hl: np.ndarray, hp: np.ndarray, hr: np.ndarray, hs1: np.ndarray, m: int, 46 model: int) -> None: 47 """Write the model's header and U, V, W and n + 1 at the first shell, 48 every nprtpt-th shell and the last one, then close the unit. 49 50 Translates: wrtmil.f#WRTMIL (7be4133) 51 52 Arguments: the model hcomp[species, shell], log rho hd, L/Lsun hl, 53 log P hp, log r hr, m in grams hs1 (index 0 unused); m, the number of 54 shells; model, the model's number. None is returned. 55 56 """ 57 58 c = constants_mod 59 io = iounits_mod 60 nprtpt = nmlpar_mod.nprtpt 61 62 # The header; Fortran's 1PD10.3 writes a D exponent (Fortran lines 63 # 66–69). 64 smtot = math.exp(c.cln * comp_mod.stotal) / c.cmsun 65 io.write(io.imilne, 66 f'{"":10}MODEL{model:5d} XENV ={_d(comp_mod.xenv)}' 67 f' ZENV ={_d(comp_mod.zenv)} MASS(SOLAR UNITS) ={_d(smtot)}') 68 69 # The first shell, always written (Fortran lines 70–79). 70 io.write(io.imilne, _row(1, hcomp, hd, hl, hp, hr, hs1)) 71 72 # Every nprtpt-th shell (Fortran lines 81–96). 73 iend = 1 74 if nprtpt <= m: 75 ibeg = max(2, nprtpt) 76 iend = m - m % nprtpt 77 for i in range(ibeg, iend + 1, nprtpt): # loop 20 78 io.write(io.imilne, _row(i, hcomp, hd, hl, hp, hr, hs1)) 79 80 # The last shell, if nprtpt does not divide m (Fortran lines 97–108). 81 if iend < m: 82 io.write(io.imilne, _row(m, hcomp, hd, hl, hp, hr, hs1)) 83 84 # Suspected bug (see the page): closed after each model, opened only 85 # once (Fortran line 109). 86 io.close(io.imilne)