wrtmod.f

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

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

Variables on this page (50)

atmbeg, atmmax, atmmin, atmstp, c4pil, cgas, cgl, cln, clsunl, cmsun, crsun, crsunl, csigl, envbeg, envmax, envmin, envstp, gam1, imodpt, iopmod, ipver, isound, json, ldh, lpulse, lsound, lstatm, lstenv, lstore, pemu, pqcp, pqdp, pqdt, pqed, pqet, pqod, pqot, prmu, sdel, sesum, so, spotf, spotx, stotal, xenv0, xxdh, yydh, zdh, zenv0, zzdh

Purpose

WRTMOD writes extra profiles of one model: with LSOUND, a sound-speed and structure table (to the units IMODPT and ISOUND, the file Csound.dat); with LSTATM, LSTENV or LPULSE, an envelope integrated again for printing; and with LPULSE, a model for pulsation codes on IOPMOD, in one of three versions (IPVER).

How it fits

Inputs: the model (arguments: M, JENV, HCOMP, HS1, $\log\rho$, $L$, $\log P$, $\log r$, $\log T$, MODEL, $\log L$, $\log T_\mathrm{eff}$, FP, FT, $\log m$ HS, the age, the mixing length CMIXL and the star's mass XMSOL, a snapshot of SMASS taken by the main program), and through modules the per-shell physics: $\Gamma_1$ and the switch LSOUND (SOUND_MOD), the pulsation derivatives PQDP, PQDT, … (PULSE1_MOD), energy rates, opacity and gradients (SCRTCH_MOD), STOTAL (COMP_MOD), and the namelist switches and integration steps (LSTORE, LSTATM, LSTENV, LPULSE, IPVER, SPOTF, SPOTX, ATMSTP, ENVSTP, … in NMLPAR_MOD; formerly CCOUT, PULSE, SPOTS, ENVGEN, INTATM, INTENV).
Who calls it: PUTSTORE (when the atmosphere or envelope is stored) and WRTOUT (for pulsation output between stored models).
What it calls: ENVINT, which itself prints the envelope and atmosphere (and the pulsation envelope with LPULPT).
Outputs: IMODPT, ISOUND (the file FSOUND, Csound.dat by default, opened once by PARMIN), IOPMOD, and the envelope printed by ENVINT.

The method

Sound speed (LSOUND). For each shell, with $V = G m \rho/(P r)$ and $U = 4\pi\rho r^3/m$: $r/R_\odot$, $m/M_\odot$, $m/r^3$ (solar units), $V/\Gamma_1$, $\Gamma_1$, $$A = -V/\Gamma_1 + V\left(\frac{\partial\ln\rho}{\partial\ln P} + \nabla\,\frac{\partial\ln\rho}{\partial\ln T}\right)$$ (the code's XX4; probably the buoyancy term, by its form), $U$, and $c_s = \sqrt{\Gamma_1 P/\rho}$ in km/s. Between shells it also writes $dP/dr$-like and $d\rho/dr$-like terms (QPR, QDR) and their differences scaled by $\Gamma_1 P\,\Delta r$ and $\rho\,\Delta r$.

Question 1 (line 138 ↗; line 140 ↗). Csound.dat's header. The header line is written again at every call to the file PARMIN opened once, so a run that writes several models repeats it between their tables, with no model number or age. Is the repeated header meant to separate the models, and should it say which model follows?
Envelope. The integration steps are set to the output steps, ENVINT integrates from the surface to the last shell's pressure (for a spotted star at the ambient temperature, as in starin.f), and the steps are restored afterwards.
Pulsation model (LPULSE). A header (LPULSE alone), then, only when both LPULSE and LSTORE are set, for each shell in the list ID(1..IDM): $\log r$, $m/M$, $L$, $\log T$, $\log\rho$, $\log P$, $\epsilon$, $\kappa$ and the derivatives, in format 1, 2 or 3 (version 3 adds three columns for the MLT quantities, written as zeros).

For users (from YREC's Quick Start). With LPULSE, the pulsation output: the interior model on IOPMOD (.pmod, written here) and the envelope and atmosphere on IOPENV and IOPATM (.penv, .patm, written by ENVINT when called from here); the Quick Start names the last .atm, the test suite's namelists .patm. IMODPT is the .full file (FMODPT).

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. (wrong output) Checked by KC (Nov 10, 2026) A typo in $d\rho/dr$.

(lines 165–168 ↗) QDR2, the $d\rho/dr$-like term at shell I+1, multiplies by QPR1, the $dP/dr$-like term at shell I, instead of QPR2, as QDR1 does with QPR1. What follows: The written QDR2-based difference QQDR (lines 170–171 ↗) is wrong in the .full file (IMODPT); nothing reads it back. Only with LSOUND, which main.f sets for the last model of a run.

A possible fix:

wrtmod.f

@@ lines 165–171 @@
               QPR1 = EXP(CLN*(CGL+HS(I)+HD(I)-2.0D0*HR(I)))
               QDR1 = EXP(CLN*(HD(I)-HP(I)))*PQDP(I)*QPR1
               QPR2 = EXP(CLN*(CGL+HS(I+1)+HD(I+1)-2.0D0*HR(I+1)))
-              QDR2 = EXP(CLN*(HD(I+1)-HP(I+1)))*PQDP(I+1)*QPR1
+C 11/09/26 KC via Claude
+C code changed to use QPR2, the term at shell I+1, in QDR2.
+              QDR2 = EXP(CLN*(HD(I+1)-HP(I+1)))*PQDP(I+1)*QPR2
               QQPR = (QPR1-QPR2)/DIVP
               QQDR = (QDR1-QDR2)/DIVR
               WRITE(IMODPT,124)RMID/CRSUN,QPR1,QDR1,QQPR,QQDR

How the outputs show it. The right QDR2 is the next shell's QDR1, which the .full file writes on the next line of format 124 (line 166 ↗; line 171 ↗; the file alternates them with the shells' format-123 lines, line 152 ↗), so the right QQDR follows from the file's own columns: the written value times $(Q_{D,i}-Q_{D,i+1})/(Q_{D,i}-Q_{D,i+1}Q_{P,i}/Q_{P,i+1})$, with $Q_P$ and $Q_D$ the second and third columns. In testsuite/Test_solarGS98_dif, at the last model (1208), the 301st line of this kind (at 0.1598544 $R_\odot$) writes QQDR as -6.19613054E-21, where the columns give 1.54E-20; the two differ in sign on 766 of the 795 lines, and the median ratio is -2.5.

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 sound-speed block ran once per run, at the last model (models 748 to 1,723), in 10 of 34 runs: seven of YREC's solar cases with LPULSE (for example, testsuite/Test_solarGS98_dif), two made-up ones built on them, and a made-up run with LSTCH off (YREC's testsuite/Test_m1p0feh+0p0GS98_base_TAMS, with the default). On every shell, QPR1 differed from QPR2 by −5.6% to +10.3%, the factor by which QDR2 is off. With the fix, testsuite/Test_solarGS98_allard and testsuite/Test_solarM22M_difkur have identical tracks and .store files and .short files within 2%; the .full file that holds QQDR is not among the files the test suite compares. Intent: none found; QDR1 uses QPR1, at its own shell.

2. (wrong output) Checked by KC (Nov 10, 2026) The pulsation model's shell list is never set.

(line 51 ↗; lines 255–256 ↗; lines 265–266 ↗) IDM and ID(J) are locals, never assigned (the "long block" removed by G. Somers 11/14, line 307 ↗, or the arguments LSHELL, JXBEG, JXEND, removed since, may have set them). The header prints IDM, and the loop DO 220 J = 1,IDM uses undefined values. What follows: With the usual static, zero-filled storage of SAVEd locals, IDM = 0, so the pulsation model (LPULSE and LSTORE) has a header with "NUMBER OF SHELLS IN MODEL= 0" and no shells; with another value, ID(J) could index out of bounds. Only with LPULSE (off by default).

A possible fix: a choice for the developers: what shells the pulsation model should list (every shell, IDM = M and ID(J) = J, is the simplest), or whether the LPULSE output is still in use.

How the outputs show it. In testsuite/Test_solarGS98_dif, all 1,186 model headers of the .pmod file read "NUMBER OF SHELLS IN MODEL= 0" (lines 255–256 ↗), and the file holds nothing but these headers (7,116 lines, six a model). The .track gives the same models' shell counts in its Shls column (wrtout.f, line 558 ↗): 709 at model 22, the first header, and 797 at model 1208, the last. The envelope and atmosphere parts of the same pulsation output, the .penv and .patm files, are written in full (34 and 57 MB).

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): WRTMOD reached the pulsation header in the nine runs with LPULSE (seven of YREC's solar cases and two made-up ones built on them), at every model not stored (1,186 to 1,701 times a run), and IDM was 0 every time: The header says the model has 0 shells, and no shell is written, so the pulsation model is empty in every run that asks for it. IDM never held another value, so ID(J) never indexed out of bounds. The fix is the developers' choice, and none was run. Intent: The comment at line 307 ↗ records the removal of a block, not what was to replace it.

3. (wrong output) Checked by KC (Nov 10, 2026) Debugging prints.

(line 127 ↗; line 131 ↗) print*,'wrtmod LSOUND 1: ' writes to the screen at every call, and print*,'wrtmod LSOUND 2: ' at every call with LSOUND.

A possible fix:

wrtmod.f

@@ lines 122–135 @@
 C (no member referenced here; residue of the obsolete-option removals).
 C G Somers END

+C 11/09/26 KC via Claude
+C code changed to remove the two debugging prints of LSOUND.
       SAVE

-         print*,'wrtmod LSOUND 1: ',LSOUND
-
       IF(LSOUND)THEN

-         print*,'wrtmod LSOUND 2: ',LSOUND
-
 CFD 10/09 Add an extra output to plot the sound speed profile easyly
 C 9/25/26 MHP via Claude
 C code changed to use the managed ISOUND unit (opened once in PARMIN

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 first print ran at every call of WRTMOD: 1,186 to 1,701 times a run in the nine runs with LPULSE, and once in the run with LSTCH off; the second once per run, at the last model. With the prints removed, testsuite/Test_solarAAG21_dif and testsuite/Test_solarGS98_allard have identical tracks and .store files and .short files within 2%; the prints go to the screen, which the test suite does not compare. Intent: The matching print in the main program is commented out (main.f, lines 670–671 ↗).

4. (never occurred) `LPRT` may be unset.

(line 186 ↗; envint.f, line 206 ↗; envint.f, line 222 ↗) LPRT is a local, set to true only if LSTORE, and never set false; otherwise ENVINT gets its SAVEd or undefined value, which decides whether it prints the envelope and atmosphere. With the usual zero-filled storage it is false, which matches the design (the pulsation model, too, needs LSTORE, line 273 ↗), so this is correct today but rests on the compiler.

A possible fix:

wrtmod.f

@@ lines 183–189 @@
 C  PRINTING OUT THE RESULTS.
 C  SET UP FLAGS AND COUNTERS.
        LSBC0 = .FALSE.
-       IF(LSTORE)LPRT = .TRUE.
+C 11/09/26 KC via Claude
+C code changed to set LPRT on every call, not only when LSTORE.
+       LPRT = LSTORE
        KATM = 0
        KENV = 0
        KSAHA = 0

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. The line ran in 10 runs (1,186 to 1,701 times a run in the nine with LPULSE), always with LSTORE on, so LPRT was set every time. Not run: LPULSE with LSTORE off, which none of YREC's cases uses. Intent: none found.

Definitions

subroutineWRTMOD

lines 9–310 · GitHub

Extra profiles of one model: sound-speed table, a reprinted envelope, and the pulsation model (whose shell list is currently undefined).

Uses modules COMP_MOD, CONSTANTS_MOD, DEBHU_MOD, DIMENSIONS_MOD, IOUNITS_MOD, LABEL_MOD, NMLPAR_MOD, PULSE1_MOD, SCRTCH_MOD, SOUND_MOD

Called by PUTSTORE, WRTOUT

Flowchart of its calls, in order (1)

if (LSTATM .OR. LSTENV .OR. LPULSE), lines 181–243

ENVINT 233

Compiler warnings (2)

gfortran -Wall -Wextra, compiling only; nothing was run.

line 228: Inequality comparison for REAL(8) [-Wcompare-reals]
line 228: Inequality comparison for REAL(8) [-Wcompare-reals]

Change log

The changes since rp_1

Changed since rp_1. CMIXL and XMSOL are now arguments (CMIXL is passed on to ENVINT); the sound-speed file is written to the unit ISOUND, opened once by PARMIN, instead of being opened on unit 500 at each call, which settles part of question 1; declarations: COMMON blocks to modules, IMPLICIT NONE.

Changed since a8a5950. Version 3 of the pulsation model writes literal zeros in its three mixing-length columns, which held PUALPHA_MOD's always-zero VALFMLT, VPHMLT, and VCMXMLT (now removed), so that readers of IPVER = 3 see the same layout; PUALPHA_MOD is no longer used here.

1C 2C 3C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 4C WRTMOD 5C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 6C SUBROUTINE WRTMOD(M,LSHELL,JXBEG,JXEND,JCORE,JENV,HCOMP,HS1,HD, 7C *HL,HP,HR,HT,LC,MODEL,BL,TEFFL,OMEGA,FP,FT,ETA2,R0,HJM,HI,HS, 8C * DAGE) ! KC 2025-05-31 9 SUBROUTINE WRTMOD(M,JENV,HCOMP,HS1,HD, 10 * *HL,HP,HR,HT,MODEL,BL,TEFFL,FP,FT,HS,DAGE,CMIXL,XMSOL) 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/LABEL/ into LABEL_MOD. 17 USE LABEL_MOD 18C 10/14/26 MHP via Claude 19C code changed to convert COMMON/SOUND/ into SOUND_MOD. 20 USE SOUND_MOD 21C 10/14/26 MHP via Claude 22C code changed to convert COMMON/SCRTCH/ into SCRTCH_MOD. 23 USE SCRTCH_MOD 24C 10/14/26 MHP via Claude 25C code changed to convert COMMON/PULSE1/ into PULSE1_MOD. 26 USE PULSE1_MOD 27C 10/14/26 MHP via Claude 28C code changed to convert COMMON/DEBHU/ into DEBHU_MOD. 29 USE DEBHU_MOD 30C 10/14/26 MHP via Claude 31C code changed to convert COMMON/COMP/ into COMP_MOD. 32 USE COMP_MOD 33C 10/14/26 MHP via Claude 34C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD. 35 USE NMLPAR_MOD 36C 10/14/26 MHP via Claude 37C code changed to convert COMMON/CONSTANTS/ into CONSTANTS_MOD. 38 USE CONSTANTS_MOD 39C 10/14/26 MHP via Claude 40C code changed to convert COMMON/IOUNITS/ into IOUNITS_MOD. 41 USE IOUNITS_MOD 42C 10/20/26 MHP via Claude 43C code changed to replace implicit typing with IMPLICIT NONE and 44C explicit declarations. 45 IMPLICIT NONE 46 REAL*8 ABEG0,AMAX0,AMIN0,ATEFFL,B,BL,CMIXL,DAGE,DIVP,DIVR,DR,DUM1, 47 * * DUM2,DUM3,DUM4,EBEG0,EMAX0,EMIN0,FM,FP,FPL,FR,FS,FSI,FT, 48 * * FTL,GL,HCOMP,HD,HL,HP,HR,HS,HS1,HSTOT,HT,PELPF,PLIM,QDR1, 49 * * QDR2,QPR1,QPR2,QQDR,QQPR,QSMASS,RL,RMID,RSURFL,SV,TEFFL, 50 * * TEMPR,X,XMSOL,XX1,XX2,XX3,XX4,XX5,XXX,Z 51 INTEGER I,ID,IDM,IDUM,IXX,J,JENV,KATM,KENV,KSAHA,M,MODEL 52 LOGICAL*4 LPRT,LPULPT,LSBC0 53 54C 9/25/26 MHP via Claude 55C code changed to consolidate all I/O logical unit numbers, previously 56C scattered across COMMON/LUOUT/, COMMON/LUNUM/, COMMON/ALEX06/, 57C COMMON/ACDPTH/, COMMON/NEWOPAC/, COMMON/LOPAL95/, COMMON/ATMOS2/, 58C COMMON/ALATM03/, COMMON/OPALEOS/, and COMMON/SCV2/, into a single 59C COMMON/IOUNITS/ 60C 9/25/26 MHP via Claude 61C code changed to add ISOUND (Csound.dat unit) to COMMON/IOUNITS/ 62C 9/24/26 MHP via Claude 63C code changed to remove dead ILLDAT, IKUR from COMMON/LUNUM/ 64C 9/24/26 MHP via Claude 65C code changed to remove dead IDYN from COMMON/LUNUM/ 66C 9/25/26 MHP via Claude 67C code changed to remove dead COMMON/LUNUM/ISCOMP (declared, 68C assigned a unit number, but never opened or read anywhere) 69C 10/11/26 MHP via Claude 70C code changed to remove LSCRIB from COMMON/CCOUT/: it was declared 71C positionally in all 17 files sharing this block but never read or 72C written in executable code anywhere, and no namelist file sets it. 73C 10/12/26 MHP via Claude 74C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 75C namelist variables that the physics routines read from 27 COMMON 76C blocks which also held derived or working members. Those blocks keep 77C only their non-namelist members; five of them -- NEWMX, VMULT, 78C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 79C Derivation-only namelist variables from the same blocks went to 80C COMMON/NMLBYPASS/ instead, not here. 81C 82C Declared identically in all 121 files that hold any of these members. 83C Generated from one member list, so every copy is byte-identical; the 84C linker's 'size of symbol changed' warning is the check that they are. 85C 10/13/26 MHP via Claude: removed the unused COMMON/CCOUT/ declaration 86C (no member referenced here; residue of the obsolete-option removals). 87C 9/25/26 MHP via Claude 88C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/, 89C and COMMON/CONST3/ into a single COMMON/CONSTANTS/ 90C 9/24/26 MHP via Claude 91C code changed to pass CMIXL as an explicit argument instead of via 92C common block, to reduce implicit global state ahead of the F90 93C module conversion 94C DBG PULSE 95C 10/13/26 MHP via Claude 96C code changed to replace COMMON/PULSE/ with an explicit argument. main.f 97C snapshots XMSOL=SMASS and wrtmod.f is the only reader; it is a snapshot, 98C not the running SMASS, so it cannot be replaced by the mass already in 99C scope there. 100C 10/13/26 MHP via Claude: removed the unused COMMON/PULSE2/ declaration 101C (no member referenced here; residue of the obsolete-option removals). 102C MHP 7/96 COMMON BLOCK INSERTED FOR SOUND SPEED 103C DBG 104C DBG 7/92 COMMON BLOCK ADDED TO COMPUTE DEBYE-HUCKEL CORRECTION. 105C DBG 7/95 To store variables for pulse output 106C 10/13/26 MHP via Claude: removed the unused COMMON/ROTEN/ declaration 107C (no member referenced here; residue of the obsolete-option removals). 108 DIMENSION HCOMP(15,JSON),DUM1(4),DUM2(3),DUM3(3),DUM4(3) 109C *HS1(JSON),HD(JSON),HL(JSON),HP(JSON),HR(JSON),HT(JSON),LC(JSON) 110C DIMENSION OMEGA(JSON),FP(JSON),FT(JSON),ETA2(JSON),ID(JSON) 111C DIMENSION R0(JSON),HJM(JSON),HI(JSON),HS(JSON) ! KC 2025-05-31 112 DIMENSION HS(JSON),HS1(JSON),HD(JSON),HL(JSON),HP(JSON),HR(JSON), 113 * * HT(JSON) 114 DIMENSION FP(JSON),FT(JSON),ID(JSON) 115C KC 2025-05-30 addressed warning messages from Makefile.legacy 116C DATA IHEADR/4H****/ 117 118C G Somers 10/14, Add spot common block, and store common block. 119C 10/13/26 MHP via Claude: removed the unused COMMON/TEMP2/ declaration 120C (no member referenced here; residue of the obsolete-option removals). 121C 10/13/26 MHP via Claude: removed the unused COMMON/QUADD/ declaration 122C (no member referenced here; residue of the obsolete-option removals). 123C G Somers END 124 125 SAVE 126 127 print*,'wrtmod LSOUND 1: ',LSOUND 128 129 IF(LSOUND)THEN 130 131 print*,'wrtmod LSOUND 2: ',LSOUND 132 133CFD 10/09 Add an extra output to plot the sound speed profile easyly 134C 9/25/26 MHP via Claude 135C code changed to use the managed ISOUND unit (opened once in PARMIN 136C against the user-configurable FSOUND filename) instead of a 137C hardcoded unit 500/'Csound.dat' 138 WRITE(ISOUND,9124) 139 9123 9123 FORMAT(1X,I6,1X,2F12.8,1P,3E16.8) 140 9124 9124 FORMAT(1X,'Shell#',5x,'R/Rsun',7x,'M/Msun',7x,'Cs',15x,'Rho') 141CFD end 142 DO I = 1,M 143 FR = EXP(CLN*HR(I))/CRSUN 144 FM = EXP(CLN*HS(I))/CMSUN 145 XX1 = FM/FR**3 146 XXX = -EXP(CLN*(CGL+HS(I)+HD(I)-HP(I)-HR(I))) 147 XX2 = -XXX/GAM1(I) 148 XX3 = GAM1(I) 149 XX4 = -XX2-XXX*(PQDP(I)+SDEL(2,I)*PQDT(I)) 150 XX5 = EXP(CLN*(C4PIL+HD(I)+3.0D0*HR(I)-HS(I))) 151 SV = 1.0D-5*SQRT(GAM1(I)*EXP(CLN*(HP(I)-HD(I)))) 152 WRITE(IMODPT,123)FR,FM,XX1,XX2,XX3,XX4,XX5,SV 153 123 123 FORMAT(1X,2F12.8,1P6E16.8) 154CFD 10/09 Add an extra output to plot the sound speed profile easyly 155 WRITE(ISOUND,9123)I,FR,FM,SV,HD(I),HT(I) 156CFD end 157 158C DERIVATIVES OF DP/DR, DRHO/DR 159 IF(I.LT.M) THEN 160 RMID = 0.5D0*(EXP(CLN*HR(I))+EXP(CLN*HR(I+1))) 161 DR = EXP(CLN*HR(I+1)) - EXP(CLN*HR(I)) 162 DIVP = 0.5D0*(GAM1(I)*EXP(CLN*HP(I))+ 163 * * GAM1(I+1)*EXP(CLN*HP(I+1)))*DR 164 DIVR = 0.5D0*(EXP(CLN*HD(I))+EXP(CLN*HD(I+1)))*DR 165 QPR1 = EXP(CLN*(CGL+HS(I)+HD(I)-2.0D0*HR(I))) 166 QDR1 = EXP(CLN*(HD(I)-HP(I)))*PQDP(I)*QPR1 167 QPR2 = EXP(CLN*(CGL+HS(I+1)+HD(I+1)-2.0D0*HR(I+1))) 168 QDR2 = EXP(CLN*(HD(I+1)-HP(I+1)))*PQDP(I+1)*QPR1 169 QQPR = (QPR1-QPR2)/DIVP 170 QQDR = (QDR1-QDR2)/DIVR 171 WRITE(IMODPT,124)RMID/CRSUN,QPR1,QDR1,QQPR,QQDR 172 124 124 FORMAT(1X,F11.7,1P4E16.8) 173 ENDIF 174 END DO 175 ENDIF 176C 177C G Somers 11/14 LCHEMO BLOCK REMOVED, AS THIS INFO IS ALREADY IN .STORE. 178C 179C DBG PULSE: PRINT OUT PULSATION ENV AND ATM IN ENVINT 180C G Somers 11/14 CHANGE TO NEW I/O FLAGS. 181 IF(LSTATM .OR. LSTENV .OR. LPULSE) THEN 182C INTEGRATE AN ENVELOPE FROM THE SURFACE TO THE CONVERGED MODEL, 183C PRINTING OUT THE RESULTS. 184C SET UP FLAGS AND COUNTERS. 185 LSBC0 = .FALSE. 186 IF(LSTORE)LPRT = .TRUE. 187 KATM = 0 188 KENV = 0 189 KSAHA = 0 190C SAVE THE INTEGRATION STEP PARAMETERS AND ENFORCE THE SPACING 191C REQUESTED FOR PRINTOUT PURPOSES. 192 ABEG0 = ATMBEG 193 AMIN0 = ATMMIN 194 AMAX0 = ATMMAX 195 EBEG0 = ENVBEG 196 EMIN0 = ENVMIN 197 EMAX0 = ENVMAX 198 ATMBEG = ATMSTP 199 ATMMIN = ATMSTP 200 ATMMAX = ATMSTP 201 ENVBEG = ENVSTP 202 ENVMIN = ENVSTP 203 ENVMAX = ENVSTP 204 B = DEXP(CLN*BL) 205 RL = 0.5D0*(BL + CLSUNL - 4.0D0*TEFFL - C4PIL - CSIGL) 206 GL = CGL + STOTAL - RL - RL 207 X = HCOMP(1,M) 208 Z = HCOMP(3,M) 209 FPL = FP(M) 210 FTL = FT(M) 211 IXX=0 212 HSTOT = STOTAL 213 PLIM = HP(M) 214C DBG PULSE: ADDED ARGUEMENT TO ENVINT TO TURN ON/OFF PULSE OUTPUT 215 LPULPT = LPULSE 216 IF (LDH) THEN 217 XXDH = HCOMP(1,M) 218 YYDH = HCOMP(2,M)+HCOMP(4,M) 219 ZZDH = HCOMP(3,M) 220 ZDH(1) = HCOMP(5,M)+HCOMP(6,M) 221 ZDH(2) = HCOMP(7,M)+HCOMP(8,M) 222 ZDH(3) = HCOMP(9,M)+HCOMP(10,M)+HCOMP(11,M) 223 END IF 224C MHP 10/02 define ISTORE - used in ENVINT 225 IDUM = 0 226C G Somers 10/14, FOR SPOTTED RUNS, FIND THE 227C PRESSURE AT THE AMBIENT TEMPERATURE ATEFFL 228 IF(JENV.EQ.M.AND.SPOTF.NE.0.0.AND.SPOTX.NE.1.0)THEN 229 ATEFFL = TEFFL - 0.25*LOG10(SPOTF * SPOTX**4.0 + 1.0 - SPOTF) 230 ELSE 231 ATEFFL = TEFFL 232 ENDIF 233 CALL ENVINT(B,FPL,FTL,GL,HSTOT,IXX,LPRT,LSBC0,PLIM,RL, 234 * * ATEFFL,X,Z,DUM1,IDUM,KATM,KENV,KSAHA,DUM2, 235 * * DUM3,DUM4,LPULPT,CMIXL) 236C G Somers END 237 ATMBEG = ABEG0 238 ATMMIN = AMIN0 239 ATMMAX = AMAX0 240 ENVBEG = EBEG0 241 ENVMIN = EMIN0 242 ENVMAX = EMAX0 243 ENDIF 244C 245C G Somers 11/14 LCONZO (convection zone info) block deleted. 246C 247C G Somers 11/14 LJOUT (rotation info) block deleted. 248C 249 FSI = DEXP(-CLN*STOTAL) 250C DBG PULSE: WRITE HEADER INFORMATION FOR PULSE MODEL 251 IF(LPULSE) THEN 252 RSURFL = 0.5D0*(BL - C4PIL - CSIGL - 4.0D0*TEFFL + CLSUNL) 253 TEMPR = RSURFL - CRSUNL 254 QSMASS = XMSOL 255 WRITE (IOPMOD, 5001) MODEL,IDM,IPVER,QSMASS, 256 1 1 TEFFL,BL,TEMPR, DAGE, CMIXL, XENV0, ZENV0 257 5001 5001 FORMAT(' MODEL#=', I5, ' NUMBER OF SHELLS IN MODEL=',I5, 258 * * ' VER=',I2,/, 259 * * ' MASS=',F8.5, ' LOG(TEFF)=',F8.5,/, ' LOG(L/LSUN)=', 260 * * F16.10, ' LOG(R/RSUN)=',F16.10, /, 261 * * ' AGE=', 1PE12.5,' GYR',/, 262 * * ' MIXING LENGTH PARAMETER=', 0PF16.10,/, 263 * * ' ZAMS (X,Z)=', 2F18.10) 264 END IF 265 DO 220 J = 1,IDM 266 I = ID(J) 267 FS = FSI*HS1(I) 268C 269C G Somers 11/14 FINAL LSCRIB BLOCK DELETED 270C 271C DBG WRITE PULSE MODEL 272C PRINT*, 'LPULSE=',LPULSE 273 IF (LPULSE.AND.LSTORE) THEN 274C MHP 10/02 uncommented pelpf statement, used later in i/o 275C PELPF = CGAS * DEXP(CLN*(HT(I) + HD(I)))* PEMU(I) 276C ADDED X AND Z TO OUTPUT 277 IF ((J.EQ.2).AND.(I.EQ.1)) GOTO 5003 278 IF(IPVER.EQ.1) THEN 279 PELPF = CGAS * DEXP(CLN*(HT(I) + HD(I)))* PEMU(I) 280 WRITE(IOPMOD, 5052)HR(I),FS,HL(I),HT(I),HD(I), 281 * * HP(I), SESUM(I),SO(I), PQDP(I), PQED(I), 282 * * PQET(I), PQOD(I), PQOT(I), SDEL(2,I),SDEL(3,I), 283 * * PQCP(I), PRMU(I), PQDT(I), PELPF 284 ELSE IF (IPVER.EQ.2) THEN 285 WRITE(IOPMOD, 6052)HR(I),FS,HL(I),HT(I),HD(I), 286 * * HP(I), SESUM(I),SO(I), PQDP(I), PQED(I), 287 * * PQET(I), PQOD(I), PQOT(I), SDEL(2,I),SDEL(3,I), 288 * * PQCP(I), PRMU(I), PQDT(I), HCOMP(1,I),HCOMP(3,I) 289 ELSE IF (IPVER.EQ.3) THEN 290C DBG 7/95 Modifed to include mixing length variables 291C 11/02/26 MHP via Claude 292C code changed to write literal zeros for the mixing-length columns, 293C which were always zero; kept so IPVER=3 readers see the same layout. 294 WRITE(IOPMOD, 6053)HR(I),FS,HL(I),HT(I),HD(I),0.0D0, 295 * * HP(I), SESUM(I),SO(I), PQDP(I), PQED(I),0.0D0, 296 * * PQET(I), PQOD(I), PQOT(I), SDEL(2,I),SDEL(3,I),0.0D0, 297 * * PQCP(I), PRMU(I), PQDT(I), HCOMP(1,I),HCOMP(3,I) 298 END IF 299 5003 5003 CONTINUE 300 5052 5052 FORMAT(5E16.9,/,5E16.9,/,5E16.9,/,5E16.9) 301 6052 6052 FORMAT(5E23.16,/,5E23.16,/,5E23.16,/,5E23.16) 302 6053 6053 FORMAT(6E23.16,/,6E23.16,/,6E23.16,/,5E23.16) 303 END IF 304C DBG END 305 220 220 CONTINUE 306C 307C G Somers 11/14 REMOVED LONG BLOCK 308C 309 RETURN 310 END

1"""Python version of wrtmod.f (YREC 7be4133): extra profiles of one model 2(the sound-speed table, an envelope integrated again for printing, and the 3model for pulsation codes).""" 4 5import math 6 7import numpy as np 8 9from pyrec import comp_mod 10from pyrec import constants_mod 11from pyrec import debhu_mod 12from pyrec import dimensions_mod 13from pyrec import iounits_mod 14from pyrec import label_mod 15from pyrec import nmlpar_mod 16from pyrec import pulse1_mod 17from pyrec import scrtch_mod 18from pyrec import sound_mod 19from pyrec.envint import envint 20 21JSON = dimensions_mod.json 22 23# Locals kept between calls (SAVE). 24# Suspected bug (see the page): IDM and ID are never assigned, so the 25# pulsation model's shell list is undefined in the Fortran; here it starts 26# empty. 27id_ = np.zeros(JSON + 1, dtype=int) 28_saved = {'idm': 0, 29 # Suspected bug (see the page): LPRT is set only if LSTORE; 30 # otherwise ENVINT gets its saved (or, on the first call, 31 # undefined) value. 32 'lprt': False} 33 34 35def _l1(flag: bool) -> str: 36 """Fortran's L1.""" 37 38 return "T" if flag else "F" 39 40 41def _e0(x: float, w: int, d: int) -> str: 42 """Fortran's 0PEw.d: a mantissa in [0.1, 1), e.g. 0.12345678E+01.""" 43 44 if x == 0.0: 45 return f"{'0.' + '0' * d + 'E+00':>{w}}" 46 e = math.floor(math.log10(abs(x))) + 1 47 m = x / 10.0**e 48 if round(abs(m), d) >= 1.0: # rounding carried into the next decade 49 m /= 10.0 50 e += 1 51 s = f"{abs(m):.{d}f}"[1:] # drop the leading zero, as Fortran does 52 sign = "-" if x < 0 else "" 53 return f"{sign + '0' + s + f'E{e:+03d}':>{w}}" 54 55 56def _records(values: list, per_line: tuple, w: int, d: int) -> str: 57 """Values in 0PEw.d, per_line[k] of them on record k (FORMATs 5052, 58 6052, 6053); the last record stops when the values run out.""" 59 60 lines = [] 61 k = 0 62 for count in per_line: 63 chunk = values[k:k + count] 64 if not chunk: 65 break 66 lines.append("".join(_e0(v, w, d) for v in chunk)) 67 k += count 68 return "\n".join(lines) 69 70 71def wrtmod(m: int, jenv: int, hcomp: np.ndarray, hs1: np.ndarray, hd: np.ndarray, hl: np.ndarray, hp: np.ndarray, 72 hr: np.ndarray, ht: np.ndarray, model: int, bl: float, teffl: float, fp: np.ndarray, ft: np.ndarray, 73 hs: np.ndarray, dage: float, cmixl: float, xmsol: float) -> None: 74 """Write extra profiles of one model: with LSOUND, a sound-speed and 75 structure table (units IMODPT and ISOUND); with LSTATM, LSTENV or 76 LPULSE, an envelope integrated again by ENVINT for printing; with 77 LPULSE, a model for pulsation codes on IOPMOD (version IPVER). 78 79 Translates: wrtmod.f#WRTMOD (7be4133) 80 81 Arguments: the number of shells m, the envelope's base jenv, the model 82 (hcomp, hs1, hd, hl, hp, hr, ht, fp, ft, hs, all logarithmic as in the 83 model arrays), the model number, log L bl, log Teff teffl, the age 84 dage, the mixing length cmixl, and xmsol, the star's mass in solar 85 masses as main.f snapshots it. Nothing is returned. 86 87 """ 88 89 n = nmlpar_mod 90 c = constants_mod 91 cln = c.cln 92 io = iounits_mod 93 gam1 = sound_mod.gam1 94 pqdp = pulse1_mod.pqdp 95 pqdt = pulse1_mod.pqdt 96 sdel = scrtch_mod.sdel 97 98 # Suspected bug (see the page): debugging prints at every call 99 # (Fortran line 127). 100 print(' wrtmod LSOUND 1: ', _l1(sound_mod.lsound)) 101 102 # The sound-speed and structure table (Fortran lines 129–175). 103 if sound_mod.lsound: 104 print(' wrtmod LSOUND 2: ', _l1(sound_mod.lsound)) 105 106 # The sound-speed profile's header, on ISOUND (opened in PARMIN) 107 # (CFD 10/09) (Fortran line 138). Question (see the page): the 108 # header is written at every call. 109 io.write(io.isound, " Shell# R/Rsun M/Msun Cs" 110 " Rho") # FORMAT 9124 111 112 # One line per shell on IMODPT and ISOUND (Fortran lines 142–156). 113 for i in range(1, m + 1): 114 fr = math.exp(cln * hr[i]) / c.crsun 115 fm = math.exp(cln * hs[i]) / c.cmsun 116 xx1 = fm / fr**3 117 xxx = -math.exp(cln * (c.cgl + hs[i] + hd[i] - hp[i] - hr[i])) 118 xx2 = -xxx / gam1[i] 119 xx3 = gam1[i] 120 xx4 = -xx2 - xxx * (pqdp[i] + sdel[2, i] * pqdt[i]) 121 xx5 = math.exp(cln * (c.c4pil + hd[i] + 3.0 * hr[i] - hs[i])) 122 sv = 1.0e-5 * math.sqrt(gam1[i] * math.exp(cln * (hp[i] - hd[i]))) 123 io.write(io.imodpt, # FORMAT 123 124 f" {fr:12.8f}{fm:12.8f}{xx1:16.8E}{xx2:16.8E}" 125 f"{xx3:16.8E}{xx4:16.8E}{xx5:16.8E}{sv:16.8E}") 126 io.write(io.isound, # FORMAT 9123 127 f" {i:6d} {fr:12.8f}{fm:12.8f}{sv:16.8E}{hd[i]:16.8E}" 128 f"{ht[i]:16.8E}") 129 130 # Derivatives of dP/dr, drho/dr (Fortran lines 158–173). 131 if i < m: 132 rmid = 0.5 * (math.exp(cln * hr[i]) + math.exp(cln * hr[i + 1])) 133 dr = math.exp(cln * hr[i + 1]) - math.exp(cln * hr[i]) 134 divp = 0.5 * (gam1[i] * math.exp(cln * hp[i]) 135 + gam1[i + 1] * math.exp(cln * hp[i + 1])) * dr 136 divr = 0.5 * (math.exp(cln * hd[i]) 137 + math.exp(cln * hd[i + 1])) * dr 138 qpr1 = math.exp(cln * (c.cgl + hs[i] + hd[i] - 2.0 * hr[i])) 139 qdr1 = math.exp(cln * (hd[i] - hp[i])) * pqdp[i] * qpr1 140 qpr2 = math.exp(cln * (c.cgl + hs[i + 1] + hd[i + 1] 141 - 2.0 * hr[i + 1])) 142 # Suspected bug (see the page): QPR1 where QPR2 is expected. 143 qdr2 = math.exp(cln * (hd[i + 1] - hp[i + 1])) * pqdp[i + 1] \ 144 * qpr1 145 qqpr = (qpr1 - qpr2) / divp 146 qqdr = (qdr1 - qdr2) / divr 147 io.write(io.imodpt, # FORMAT 124 148 f" {rmid / c.crsun:11.7f}{qpr1:16.8E}{qdr1:16.8E}" 149 f"{qqpr:16.8E}{qqdr:16.8E}") 150 151 # Integrate an envelope from the surface to the converged model, 152 # printing the results (ENVINT prints them) (Fortran lines 181–243). 153 if n.lstatm or n.lstenv or n.lpulse: 154 155 # Flags and counters (Fortran lines 185–189). 156 lsbc0 = False 157 if n.lstore: 158 _saved['lprt'] = True 159 lprt = _saved['lprt'] 160 katm = 0 161 kenv = 0 162 ksaha = 0 163 164 # Save the integration step parameters and enforce the spacing 165 # requested for printout (Fortran lines 192–203). 166 abeg0 = n.atmbeg 167 amin0 = n.atmmin 168 amax0 = n.atmmax 169 ebeg0 = n.envbeg 170 emin0 = n.envmin 171 emax0 = n.envmax 172 n.atmbeg = n.atmstp 173 n.atmmin = n.atmstp 174 n.atmmax = n.atmstp 175 n.envbeg = n.envstp 176 n.envmin = n.envstp 177 n.envmax = n.envstp 178 179 # The surface quantities for ENVINT (Fortran lines 204–215). 180 b = math.exp(cln * bl) 181 rl = 0.5 * (bl + c.clsunl - 4.0 * teffl - c.c4pil - c.csigl) 182 gl = c.cgl + comp_mod.stotal - rl - rl 183 x = hcomp[1, m] 184 z = hcomp[3, m] 185 fpl = fp[m] 186 ftl = ft[m] 187 ixx = 0 188 hstot = comp_mod.stotal 189 plim = hp[m] 190 lpulpt = n.lpulse 191 192 # The Debye-Hueckel composition (Fortran lines 216–223). 193 if n.ldh: 194 debhu_mod.xxdh = hcomp[1, m] 195 debhu_mod.yydh = hcomp[2, m] + hcomp[4, m] 196 debhu_mod.zzdh = hcomp[3, m] 197 debhu_mod.zdh[1] = hcomp[5, m] + hcomp[6, m] 198 debhu_mod.zdh[2] = hcomp[7, m] + hcomp[8, m] 199 debhu_mod.zdh[3] = hcomp[9, m] + hcomp[10, m] + hcomp[11, m] 200 idum = 0 201 202 # For a spotted star, the ambient temperature (G. Somers 10/14) 203 # (Fortran lines 225–232). 204 if jenv == m and n.spotf != 0.0 and n.spotx != 1.0: 205 ateffl = teffl - 0.25 * math.log10( 206 n.spotf * n.spotx**4.0 + 1.0 - n.spotf) 207 else: 208 ateffl = teffl 209 210 # The envelope, integrated and printed (Fortran lines 233–235). 211 dum1 = np.zeros(5) 212 dum2 = np.zeros(4) 213 dum3 = np.zeros(4) 214 dum4 = np.zeros(4) 215 rl, idum, katm, kenv, ksaha = envint( 216 b, fpl, ftl, gl, hstot, ixx, lprt, lsbc0, plim, rl, ateffl, x, 217 z, dum1, idum, katm, kenv, ksaha, dum2, dum3, dum4, lpulpt, 218 cmixl) 219 220 # The integration steps restored (Fortran lines 237–242). 221 n.atmbeg = abeg0 222 n.atmmin = amin0 223 n.atmmax = amax0 224 n.envbeg = ebeg0 225 n.envmin = emin0 226 n.envmax = emax0 227 228 # Fortran line 249. 229 fsi = math.exp(-cln * comp_mod.stotal) 230 idm = _saved['idm'] 231 232 # The pulsation model's header (Fortran lines 251–264). 233 if n.lpulse: 234 rsurfl = 0.5 * (bl - c.c4pil - c.csigl - 4.0 * teffl + c.clsunl) 235 tempr = rsurfl - c.crsunl 236 qsmass = xmsol 237 io.write(io.iopmod, # FORMAT 5001 238 f" MODEL#={model:5d} NUMBER OF SHELLS IN MODEL={idm:5d}" 239 f" VER={n.ipver:2d}\n" 240 f" MASS={qsmass:8.5f} LOG(TEFF)={teffl:8.5f}\n" 241 f" LOG(L/LSUN)={bl:16.10f} LOG(R/RSUN)={tempr:16.10f}\n" 242 f" AGE={dage:12.5E} GYR\n" 243 f" MIXING LENGTH PARAMETER={cmixl:16.10f}\n" 244 f" ZAMS (X,Z)={label_mod.xenv0:18.10f}" 245 f"{label_mod.zenv0:18.10f}") 246 247 # The pulsation model, shell by shell, on the list ID(1..IDM) (Fortran 248 # lines 265–305). 249 p = pulse1_mod 250 s = scrtch_mod 251 for j in range(1, idm + 1): # loop 220 252 i = id_[j] 253 fs = fsi * hs1[i] 254 if not (n.lpulse and n.lstore): 255 continue 256 if j == 2 and i == 1: # GOTO 5003 257 continue 258 if n.ipver == 1: 259 pelpf = c.cgas * math.exp(cln * (ht[i] + hd[i])) * p.pemu[i] 260 vals = [hr[i], fs, hl[i], ht[i], hd[i], 261 hp[i], s.sesum[i], s.so[i], p.pqdp[i], p.pqed[i], 262 p.pqet[i], p.pqod[i], p.pqot[i], s.sdel[2, i], 263 s.sdel[3, i], 264 p.pqcp[i], p.prmu[i], p.pqdt[i], pelpf] 265 io.write(io.iopmod, _records(vals, (5, 5, 5, 5), 16, 9)) 266 elif n.ipver == 2: 267 vals = [hr[i], fs, hl[i], ht[i], hd[i], 268 hp[i], s.sesum[i], s.so[i], p.pqdp[i], p.pqed[i], 269 p.pqet[i], p.pqod[i], p.pqot[i], s.sdel[2, i], 270 s.sdel[3, i], 271 p.pqcp[i], p.prmu[i], p.pqdt[i], hcomp[1, i], 272 hcomp[3, i]] 273 io.write(io.iopmod, _records(vals, (5, 5, 5, 5), 23, 16)) 274 elif n.ipver == 3: 275 # The mixing-length columns (DBG 7/95), always zero, written 276 # as literal zeros to keep the layout. 277 vals = [hr[i], fs, hl[i], ht[i], hd[i], 0.0, 278 hp[i], s.sesum[i], s.so[i], p.pqdp[i], p.pqed[i], 279 0.0, 280 p.pqet[i], p.pqod[i], p.pqot[i], s.sdel[2, i], 281 s.sdel[3, i], 0.0, 282 p.pqcp[i], p.prmu[i], p.pqdt[i], hcomp[1, i], 283 hcomp[3, i]] 284 io.write(io.iopmod, _records(vals, (6, 6, 6, 5), 23, 16))

wrtmod.f

Purpose

How it fits

The method

Suspected bugs, checked by running the code

1. (wrong output) Checked by KC (Nov 10, 2026) A typo in $d\rho/dr$.

2. (wrong output) Checked by KC (Nov 10, 2026) The pulsation model's shell list is never set.

3. (wrong output) Checked by KC (Nov 10, 2026) Debugging prints.

4. (never occurred) LPRT may be unset.

Definitions

subroutineWRTMOD

Change log

4. (never occurred) `LPRT` may be unset.