pdist.f

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

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

Variables on this page (10)

iopatm, iopenv, iopmod, iowr, ishort, lpulse, poa, pob, poc, pomax

Purpose

Decides when to write a model for pulsation calculations: it "calculates the distance the star has travelled in the HR diagram; if far enough, output pulsation model." When the accumulated distance passes a threshold, it opens three new output files named after the run and the model number, and sets LPULSE so that the output routines write the model to them.

How it fits

Inputs: the previous point (POL1 $=\log L$, POT1 $=\log T_\mathrm{eff}$, POA1 = age) and the running total POLEN; the current BL ($\log L/L_\odot$), TEFFL, MODELN, the run index NK and the age DAGE (arguments); the weights POA, POB, POC and the threshold POMAX in NMLPAR_MOD (namelist CONTROL, read by PARMIN); the base file names FPMOD, FPENV, FPATM (arguments, since MHP 8/25).
Who calls it: the main program in main.f, after each model when LPOUT is set.
Outputs: POLEN and the previous point updated; LPULSE in NMLPAR_MOD; the units IOPMOD, IOPENV, IOPATM (from IOUNITS_MOD) opened as <FPMOD>NN_MMMM and so on, with NN = NK (the main program's loop over runs, DO 200 NK = 1, NUMRUN) and MMMM the model number; a line to IOWR and ISHORT. While LPULSE is true, WRTMOD writes the model to IOPMOD and passes LPULPT on to ENVINT, which writes the atmosphere and envelope to IOPATM and IOPENV from the points QATM and QENV store.

Question 1 (line 94 ↗; line 103 ↗; line 112 ↗). STATUS='NEW'. The files are opened as new, so rerunning a model sequence in the same directory stops at the first existing pulsation file. Is the stop meant to protect earlier pulsation files, or should they be opened with STATUS='UNKNOWN' (and rewound), as PARMIN opens the other outputs?

The method

The squared step in the HR diagram and in age, with weights $a, b, c$, is accumulated: $$\ell^2 \mathrel{+}= (a\,\Delta\log L)^2 + (b\,\Delta\log T_\mathrm{eff})^2 + (c\,\Delta t)^2 .$$ When $\ell^2 \ge$ POMAX, the files are opened and $\ell^2$ is reset to zero; otherwise the three units are closed and LPULSE is set false.

POMAX bounds a squared distance. POLEN adds the squares of the weighted steps, $(\mathtt{POA}\,\Delta\log L)^2 + (\mathtt{POB}\,\Delta\log T_\mathrm{eff})^2 + (\mathtt{POC}\,\Delta\mathrm{age})^2$ (lines 76–79 ↗; printed as "LEN SQ"), and is compared with POMAX (line 83 ↗), so the default POMAX = 0.1 means a distance of about 0.32. YREC's documentation describes none of POMAX, POA, POB, and POC; it could say that the threshold is a squared distance.

From YREC's documentation

Control Options (Output): LPULSE writes pulsation output, to FPMOD (the interior), FPENV (the envelope) and FPATM (the atmosphere), and needs LSTATM and LSTENV for the outer two; IPVER chooses its format.

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) Two different lines under the same labels.

(lines 119–122 ↗) The same format writes the step's differences (TM1–TM3) to IOWR but their weighted squares (TMP1–TMP3) to ISHORT, under the same labels "D L," "D Te," "D AGE," so the two files disagree and the ISHORT line is mislabeled. Only with LPOUT (off by default).

A possible fix: a choice for the developers: write the same quantities to both units, or give the weighted squares their own labels.

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. PDIST runs only with LPOUT on; its default is false, no case sets it, and no document offers it. LPOUT is in the namelist group that PARMIN reads (parmin.f, line 355 ↗), so any case with it on reaches this line at every model. The fix (measurement only: the differences written to ISHORT as well) was not run. Intent: none found.

Definitions

subroutinePDIST

lines 14–125 · GitHub

Uses modules IOUNITS_MOD, NMLPAR_MOD

Called by (file level)

Change log

The changes since rp_1

Changed since rp_1. NK and DAGE are now arguments (from /ZRAMP/ and /THEAGE/ before; main.f is the only caller); the other COMMON blocks to modules, IMPLICIT NONE; the unused /PULSE/ members and dead /LUNUM/ units dropped.

1C 2C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 3C PDIST 4C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 5 6C calculates the distance the star has travelled in the HR diagram 7C if far enough, output pulsation model 8C MHP 8/25 passed file names directly and removed common block with file names 9C 10/9/26 MHP via Claude 10C code changed to pass NK as an explicit argument instead of via 11C COMMON/ZRAMP/, which this file otherwise never used (NK was its only 12C member ever referenced here). See the CALL PDIST site in main.f, 13C which already has NK as its own local run-index variable. 14 SUBROUTINE PDIST(POL1,POT1,POA1,POLEN,BL,TEFFL,MODELN,NK,FPATM, 15 * * FPENV,FPMOD,DAGE) 16C 10/14/26 MHP via Claude 17C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD. 18 USE NMLPAR_MOD 19C 10/14/26 MHP via Claude 20C code changed to convert COMMON/IOUNITS/ into IOUNITS_MOD. 21 USE IOUNITS_MOD 22C 10/20/26 MHP via Claude 23C code changed to replace implicit typing with IMPLICIT NONE and 24C explicit declarations. 25 IMPLICIT NONE 26 REAL*8 BL,DAGE,POA1,POL1,POLEN,POT1,TEFFL,TM1,TM2,TM3,TMP1,TMP2, 27 * * TMP3 28 INTEGER IOCCOL,MODELN,NK 29 30 CHARACTER*256 COUT,CTEMP 31C Removed variables not used in the routine 32 CHARACTER*256 FPMOD, FPENV, FPATM 33C MHP 10/02 added proper dimensions to last 2 variables 34C 10/13/26 MHP via Claude 35C code changed to receive DAGE, the stellar age, as an argument instead of 36C through COMMON/THEAGE/. main.f is PDIST's only caller. 37C 10/12/26 MHP via Claude 38C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 39C namelist variables that the physics routines read from 27 COMMON 40C blocks which also held derived or working members. Those blocks keep 41C only their non-namelist members; five of them -- NEWMX, VMULT, 42C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 43C Derivation-only namelist variables from the same blocks went to 44C COMMON/NMLBYPASS/ instead, not here. 45C 46C Declared identically in all 121 files that hold any of these members. 47C Generated from one member list, so every copy is byte-identical; the 48C linker's 'size of symbol changed' warning is the check that they are. 49C 10/13/26 MHP via Claude 50C code changed to drop COMMON/PULSE/, declared here but never used. 51C MHP 8/25 Removed character file names from common block 52C 9/24/26 MHP via Claude 53C code changed to remove dead ILLDAT, IKUR from COMMON/LUNUM/ 54C 9/24/26 MHP via Claude 55C code changed to remove dead IDYN from COMMON/LUNUM/ 56C 9/25/26 MHP via Claude 57C code changed to consolidate all I/O logical unit numbers, previously 58C scattered across COMMON/LUOUT/, COMMON/LUNUM/, COMMON/ALEX06/, 59C COMMON/ACDPTH/, COMMON/NEWOPAC/, COMMON/LOPAL95/, COMMON/ATMOS2/, 60C COMMON/ALATM03/, COMMON/OPALEOS/, and COMMON/SCV2/, into a single 61C COMMON/IOUNITS/ 62C 9/25/26 MHP via Claude 63C code changed to remove dead COMMON/LUNUM/ISCOMP (declared, 64C assigned a unit number, but never opened or read anywhere) 65C MHP 8/25 Removed common block with file names 66C COMMON/LUFNM/ FLAST, FFIRST, FRUN, FSTAND, FFERMI, 67C 1 FDEBUG, FTRACK, FSHORT, FMILNE, FMODPT, 68C 2 FSTOR, FPMOD, FPENV, FPATM, FDYN, 69C 3 FLLDAT, FSNU, FSCOMP, FKUR, 70C 4 FMHD1, FMHD2, FMHD3, FMHD4, FMHD5, FMHD6, FMHD7, FMHD8 71 SAVE 72 73 TM1 = BL-POL1 74 TM2 = TEFFL-POT1 75 TM3 = DAGE-POA1 76 TMP1 = POA*TM1*POA*TM1 77 TMP2 = POB*TM2*POB*TM2 78 TMP3 = POC*TM3*POC*TM3 79 POLEN = POLEN+TMP1+TMP2+TMP3 80 POL1 = BL 81 POT1 = TEFFL 82 POA1 = DAGE 83 IF (POLEN .GE. POMAX) THEN 84 LPULSE = .TRUE. 85 POLEN = 0.0D0 86 IOCCOL = INDEX(FPMOD,' ')-1 87 IF (MODELN.LT.10000)THEN 88 WRITE(CTEMP,'(I2.2,''_'',I4.4)')NK, MODELN 89 COUT = FPMOD(1:IOCCOL) // CTEMP(1:7) 90 ELSE 91 WRITE(CTEMP,'(I2.2,''_'',I5.5)')NK, MODELN 92 COUT = FPMOD(1:IOCCOL) // CTEMP(1:8) 93 ENDIF 94 OPEN(IOPMOD,FILE=COUT,STATUS='NEW',FORM='FORMATTED') 95 IOCCOL = INDEX(FPENV,' ')-1 96 IF (MODELN.LT.10000)THEN 97 WRITE(CTEMP,'(I2.2,''_'',I4.4)')NK, MODELN 98 COUT = FPENV(1:IOCCOL) // CTEMP(1:7) 99 ELSE 100 WRITE(CTEMP,'(I2.2,''_'',I5.5)')NK, MODELN 101 COUT = FPENV(1:IOCCOL) // CTEMP(1:8) 102 ENDIF 103 OPEN(IOPENV,FILE=COUT,STATUS='NEW',FORM='FORMATTED') 104 IOCCOL = INDEX(FPATM,' ')-1 105 IF (MODELN.LT.10000)THEN 106 WRITE(CTEMP,'(I2.2,''_'',I4.4)')NK, MODELN 107 COUT = FPATM(1:IOCCOL) // CTEMP(1:7) 108 ELSE 109 WRITE(CTEMP,'(I2.2,''_'',I5.5)')NK, MODELN 110 COUT = FPATM(1:IOCCOL) // CTEMP(1:8) 111 ENDIF 112 OPEN(IOPATM,FILE=COUT,STATUS='NEW',FORM='FORMATTED') 113 ELSE 114 CLOSE(IOPMOD) 115 CLOSE(IOPENV) 116 CLOSE(IOPATM) 117 LPULSE = .FALSE. 118 END IF 119 WRITE(IOWR,276)POLEN,TM1,TM2,TM3 120 WRITE(ISHORT,276)POLEN,TMP1,TMP2,TMP3 121 276 276 FORMAT(5X,' LEN SQ=',1PE10.3, ' D L=',1PE10.3, 122 * * ' D Te=',1PE10.3,' D AGE=',1PE10.3) 123 124 RETURN 125 END

1"""Python version of pdist.f (YREC 7be4133): when to write a model for 2pulsation calculations.""" 3 4from pyrec import iounits_mod 5from pyrec import nmlpar_mod 6 7 8def _fmt276(polen: float, d1: float, d2: float, d3: float) -> str: 9 """FORMAT 276: 5X,' LEN SQ=',1PE10.3,' D L=',1PE10.3,' D Te=',1PE10.3, 10 ' D AGE=',1PE10.3. 11 12 """ 13 14 return (" " * 5 + f" LEN SQ={polen:10.3E} D L={d1:10.3E}" 15 f" D Te={d2:10.3E} D AGE={d3:10.3E}") 16 17 18def _pulse_name(base: str, nk: int, modeln: int) -> str: 19 """The file name <base>NN_MMMM (or NN_MMMMM from model 10000 on). 20 21 base is a CHARACTER*256 in Fortran, padded with blanks: the name ends 22 at its first blank, INDEX(base, ' ') - 1. A Python string has no 23 padding, so a name without blanks is taken whole (Fortran lines 24 86–93, 95–102, 104–111). 25 26 """ 27 28 ioccol = base.find(" ") if " " in base else len(base) 29 if modeln < 10000: 30 ctemp = f"{nk:02d}_{modeln:04d}" 31 else: 32 ctemp = f"{nk:02d}_{modeln:05d}" 33 return base[:ioccol] + ctemp 34 35 36def pdist(pol1: float, pot1: float, poa1: float, polen: float, bl: float, teffl: float, modeln: int, nk: int, 37 fpatm: str, fpenv: str, fpmod: str, dage: float) -> tuple[float, float, float, float]: 38 """Accumulate the squared distance the star has moved in the HR diagram 39 and in age since the last pulsation model; when it reaches POMAX, open 40 three new pulsation output files and set LPULSE, otherwise close them 41 and clear LPULSE. 42 43 Translates: pdist.f#PDIST (7be4133) 44 45 Arguments: pol1, pot1, poa1 (log L, log Teff and age at the previous 46 model), polen (the accumulated squared distance); bl, teffl, the 47 model number modeln, the run index nk and the age dage now; the base 48 file names fpatm, fpenv, fpmod. 49 50 Returns (pol1, pot1, poa1, polen), updated. Sets nmlpar_mod.lpulse. 51 52 """ 53 54 # The step and its weighted squares, accumulated (Fortran lines 73–82). 55 # Question (see the page): POMAX is compared with a squared distance. 56 tm1 = bl - pol1 57 tm2 = teffl - pot1 58 tm3 = dage - poa1 59 tmp1 = nmlpar_mod.poa * tm1 * nmlpar_mod.poa * tm1 60 tmp2 = nmlpar_mod.pob * tm2 * nmlpar_mod.pob * tm2 61 tmp3 = nmlpar_mod.poc * tm3 * nmlpar_mod.poc * tm3 62 polen = polen + tmp1 + tmp2 + tmp3 63 pol1 = bl 64 pot1 = teffl 65 poa1 = dage 66 67 # Far enough: open the three pulsation files and set LPULSE; otherwise 68 # close them (Fortran lines 83–118). Question (see the page): 69 # STATUS='NEW' stops if the file exists. 70 if polen >= nmlpar_mod.pomax: 71 nmlpar_mod.lpulse = True 72 polen = 0.0 73 iounits_mod.open(iounits_mod.iopmod, _pulse_name(fpmod, nk, modeln), 74 status="new") 75 iounits_mod.open(iounits_mod.iopenv, _pulse_name(fpenv, nk, modeln), 76 status="new") 77 iounits_mod.open(iounits_mod.iopatm, _pulse_name(fpatm, nk, modeln), 78 status="new") 79 else: 80 iounits_mod.close(iounits_mod.iopmod) 81 iounits_mod.close(iounits_mod.iopenv) 82 iounits_mod.close(iounits_mod.iopatm) 83 nmlpar_mod.lpulse = False 84 85 # The step, to IOWR and ISHORT (Fortran lines 119–122). Suspected bug 86 # (see the page): the differences go to IOWR, their weighted squares 87 # to ISHORT, under the same labels. 88 iounits_mod.write(iounits_mod.iowr, _fmt276(polen, tm1, tm2, tm3)) 89 iounits_mod.write(iounits_mod.ishort, _fmt276(polen, tmp1, tmp2, tmp3)) 90 91 return pol1, pot1, poa1, polen