qenv.f

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

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

Variables on this page (36)

c4pil, cgl, cln, cmixlq, ebeta, ed, edel, efxion, eo, ep, er, es, et, evel, iovim, lpumod, qd, qdel, qdela, qdl, qemu, qfs, qo, qol, qp, qpl, qqcp, qqdp, qqdt, qqed, qqod, qqot, qrmu, qt, qtl, stotal

Purpose

The right-hand side of the envelope integration: with $\log P$ as the independent variable, it returns the derivatives of $\log m$ (relative to the total), $\log T$ and $\log r$, from hydrostatic equilibrium, mass conservation and the temperature gradient.

How it fits

Inputs: X0 $=\log P$; Y(1..3) $= \log(m/M_\ast)$, $\log T$, $\log r$; the luminosity B, the distortion factors FPL, FTL, TEFFL, $X$, $Z$; STOTAL ($\log M_\ast$) from COMP_MOD; the mixing length CMIXLQ from QENVMIX_MOD, set by ENVINT before every path that reaches QENV (the argument list must stay that of the EXTERNAL derivative routine, so CMIXL cannot be passed); LPUMOD from PULSE1_MOD.
Who calls it: ENVINT, directly and as the EXTERNAL derivative routine passed to BSSTEP and MMID.
What it calls: EQSTAT, GETOPAC, and TPGRAD (with CMIXLQ, and with IOVIM = -1 in DPMIX_MOD, so no overshoot point is flagged).
Outputs: DYDX(1..3); the counter KENV; the point's structure, gradients, $\beta$, convective velocity and ionization fractions in ENVPRT_MOD (always, since 07/02, 08/13 and 08/25), read by ENVINT, GETNEWENV and STARIN; and the quantities for pulsation output in PULSE2_MOD when saving. EFXION, QQDP, QQDT and QQCP are set always (lines 89–94 ↗) and again, to the same values, inside IF(LSAVE .OR. LPUMOD) (lines 107–109 ↗, line 120 ↗, line 125 ↗, line 127 ↗); the second assignments are redundant.

The method

With $f_P$ = FPL, $$\frac{d\log m}{d\log P} = -\frac{4\pi r^4 P}{G m^2 f_P},\qquad \frac{d\log T}{d\log P} = \nabla,\qquad \frac{d\log r}{d\log P} = -\frac{r P}{G m \rho}\,f_P,$$ with $\nabla$ from TPGRAD (lines 72–74 ↗; on $f_P$ in the last, suspected bug 1).

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. (harmless in the runs) The distortion factor in $d\log r/d\log P$.

(line 72 ↗; line 74 ↗) The interior equations in COEFFT have $d\log P/d\log m \propto f_P$ and $d\log r/d\log m$ without $f_P$ (coefft.f, lines 247–252 ↗); dividing gives $d\log r/d\log P \propto 1/f_P$. QENV divides by $f_P$ in $d\log m/d\log P$ (consistent) but multiplies by it in $d\log r/d\log P$. What follows: In a rotating model the slope $d\log r/d\log P$ is $f_P^2$ times the one COEFFT's form gives, so the fitting point's $\log r$ (RS) and the boundary plane of SURFBC are off; the effect is small where $f_P \approx 1$. Only with rotation (LROT, off by default; CRRECT passes $f_P$ of the outermost point, otherwise 1).

A possible fix:

qenv.f

@@ lines 71–77 @@
      *     LDERIV,LCONV,FPL,FTL,TEFFL,CMIXLQ)
       DYDX(1) = -DEXP(CLN*(C4PIL+4.0D0*RL+PL-CGL-SL-SL))/FPL
       DYDX(2) = DEL
-      DYDX(3) = -DEXP(CLN*(PL+RL-CGL-SL-DL))*FPL
+C 11/09/26 KC via Claude
+C code changed to divide d log r/d log P by FPL, as the interior
+C equations of COEFFT imply (FPL multiplies d log P/d log m only).
+      DYDX(3) = -DEXP(CLN*(PL+RL-CGL-SL-DL))/FPL
       KENV = KENV + 1
 C 07/02 ALWAYS STORE THE BASIC STRUCTURE VARIABLES.
       EP = PL

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 condition ($f_P \ne 1$) held in all five rotating runs 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; one of these, with FCORE = 1.0, stopped in its first second), from about 2.4 million to 11.9 million envelope steps each, with FPL from 0.924 to 3.21. With the page's fix, in the rotating 3 $M_\odot$ case (testsuite/Test_3Msolar_GS_rot_gray_OPALSCV_SF2_GS98_OP_CF10_RGp2diffrotI3cz2) the track moves within the noise ($\log L$ by up to 0.00515 dex, against 0.00758 for a null change, nearly every opacity (all with $3 \le \log T \le 9$ and $-6 \le \log\rho \le 9.75$) times $1 + 10^{-10}$; 1,844 models instead of 1,853, against 1,839). In the rotating 0.8 $M_\odot$ case (testsuite/Test_m0p8feh+0p0GS98_rotp8p0d0p1mag1e4std_TAMS), $\log T_\mathrm{eff}$ moves by up to 0.000105 and central $X$ by 0.000155, against 0.0000272 and 0.0000336 for a null change, four to five times the noise, while $\log L$ (0.000224 against 0.00014) and $\log T_c$ stay at it: not clearly above the noise. Both runs end normally. The slope is wrong in every rotating model (FPL reached 3.21), and the error (a factor $f_P^2$ in the slope) grows as $f_P$ departs from 1; the runs measured no clear effect at these rotation rates, and faster rotators were not run. Intent: none found.

Definitions

subroutineQENV

lines 6–136 · GitHub

Uses modules COMP_MOD, CONSTANTS_MOD, DPMIX_MOD, ENVPRT_MOD, PULSE1_MOD, PULSE2_MOD, QENVMIX_MOD

Called by ENVINT ×2

Flowchart of its calls, in order (3)

EQSTAT 64→GETOPAC 67→TPGRAD 69

Change log

The changes since rp_1

Changed since rp_1. The MHD equation of state removed: EQSTAT is always called (no LMHD, no MEQOS); TPGRAD takes the mixing length as a new last argument, CMIXLQ from QENVMIX_MOD; COMMON blocks to modules (the unused /PULSE/, /FLAG/ and /MHD/ dropped), IMPLICIT NONE.

1C 2C 3C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 4C QENV 5C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 6 SUBROUTINE QENV(X0,Y,DYDX,B,FPL,FTL,GL,LATMO,LDERIV,LOCOND, 7 * * LSAVE,RL,TEFFL,X,Z,KENV,KSAHA) 8C 10/14/26 MHP via Claude 9C code changed to convert COMMON/QENVMIX/ into QENVMIX_MOD. 10 USE QENVMIX_MOD 11C 10/14/26 MHP via Claude 12C code changed to convert COMMON/PULSE2/ into PULSE2_MOD. 13 USE PULSE2_MOD 14C 10/14/26 MHP via Claude 15C code changed to convert COMMON/PULSE1/ into PULSE1_MOD. 16 USE PULSE1_MOD 17C 10/14/26 MHP via Claude 18C code changed to convert COMMON/ENVPRT/ into ENVPRT_MOD. 19 USE ENVPRT_MOD 20C 10/14/26 MHP via Claude 21C code changed to convert COMMON/DPMIX/ into DPMIX_MOD. 22 USE DPMIX_MOD 23C 10/14/26 MHP via Claude 24C code changed to convert COMMON/COMP/ into COMP_MOD. 25 USE COMP_MOD 26C 10/14/26 MHP via Claude 27C code changed to convert COMMON/CONSTANTS/ into CONSTANTS_MOD. 28 USE CONSTANTS_MOD 29C 10/20/26 MHP via Claude 30C code changed to replace implicit typing with IMPLICIT NONE and 31C explicit declarations. 32 IMPLICIT NONE 33 REAL*8 AMU,B,BETA,BETA1,BETA14,D,DEL,DELA,DELR,DL,DYDX,EMU,ETA, 34 * * FPL,FTL,FXION,GL,O,OL,P,PL,QACP,QACR,QACT,QAP,QAT,QCP,QCPP, 35 * * QCPT,QDP,QDT,QDTP,QDTT,QOD,QOT,RL,RMU,SL,T,TEFFL,TL,VEL,X, 36 * * X0,Y,Z 37 INTEGER KENV,KSAHA 38 LOGICAL*4 LATMO,LCONV,LDERIV,LOCOND,LSAVE 39C 9/24/26 MHP via Claude 40C QENV's argument list must stay fixed: it is invoked polymorphically as 41C the DERIV callback by MMID/BSSTEP (see EXTERNAL QENV in envint.f), 42C which always calls it with this exact original argument list. CMIXL 43C cannot be added here as a normal argument; it is relayed instead via 44C COMMON/QENVMIX/, set by ENVINT immediately before every code path that 45C may reach QENV (direct calls and the BSSTEP/MMID callback alike). 46 47 48 DIMENSION Y(3),DYDX(3),FXION(3) 49C 10/13/26 MHP via Claude: removed the unused COMMON/NMLPAR/ declaration 50C (no member referenced here; residue of the obsolete-option removals). 51C 10/13/26 MHP via Claude 52C code changed to drop COMMON/PULSE/, declared here but never used. 53C 9/25/26 MHP via Claude 54C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/, 55C and COMMON/CONST3/ into a single COMMON/CONSTANTS/ 56 SAVE 57 58 PL = X0 59 SL = Y(1) + STOTAL 60 TL = Y(2) 61 RL = Y(3) 62C 9/22/26 MHP via Claude 63C code changed to remove the MHD equation of state 64 CALL EQSTAT(TL,T,PL,P,DL,D,X,Z,BETA,BETA1,BETA14,FXION,RMU, 65 * * AMU,EMU,ETA,QDT,QDP,QCP,DELA,QDTT,QDTP,QAT,QAP,QCPT,QCPP, 66 * * LDERIV,LATMO,KSAHA) 67 CALL GETOPAC(DL, TL, X, Z, O, OL, QOD, QOT, FXION) 68 IOVIM = -1 69 CALL TPGRAD(TL,T,PL,P,D,RL,SL,B,O,QDT,QDP,QOT,QOD,QCP,DEL, 70 * * DELR,DELA,QDTT,QDTP,QAT,QAP,QACT,QACP,QACR,QCPT,QCPP,VEL, 71 * * LDERIV,LCONV,FPL,FTL,TEFFL,CMIXLQ) 72 DYDX(1) = -DEXP(CLN*(C4PIL+4.0D0*RL+PL-CGL-SL-SL))/FPL 73 DYDX(2) = DEL 74 DYDX(3) = -DEXP(CLN*(PL+RL-CGL-SL-DL))*FPL 75 KENV = KENV + 1 76C 07/02 ALWAYS STORE THE BASIC STRUCTURE VARIABLES. 77 EP = PL 78 ET = TL 79 ES = SL - STOTAL 80 ER = RL 81 ED = DL 82 EVEL = VEL 83C JVS 08/13 ALWAYS STORE GRADIENTS (FOR TRACKING CZ) 84 EDEL(1) = DELR 85 EDEL(2) = DELA 86 EDEL(3) = DEL 87 EBETA = BETA ! added 03/14 88C JVS 08/25 ALSO ALWAYS SAVE ADDITIONAL INFO FOR PROFILE 89 EFXION(1) = FXION(1) 90 EFXION(2) = FXION(2) 91 EFXION(3) = FXION(3) 92 QQDP = QDP 93 QQDT = QDT 94 QQCP = QCP 95 96 IF(LSAVE .OR. LPUMOD) THEN 97C EP = PL 98C ET = TL 99C ES = SL - STOTAL 100C ER = RL 101C ED = DL 102 EO = O 103C EBETA = BETA 104C EDEL(1) = DELR 105C EDEL(2) = DELA 106C EDEL(3) = DEL 107 EFXION(1) = FXION(1) 108 EFXION(2) = FXION(2) 109 EFXION(3) = FXION(3) 110C EVEL = VEL 111 QTL = TL 112 QT = DEXP(CLN*TL) 113 QPL = PL 114 QP = DEXP(CLN*PL) 115 QDL = DL 116 QD = DEXP(CLN*DL) 117 QO = O 118 QOL = OL 119 QFS = DEXP(CLN*(SL-STOTAL)) 120 QQDP = QDP 121 QQED = 0.0D0 122 QQOD = QOD 123 QQOT = QOT 124 QDEL = DEL 125 QQDT = QDT 126 QDELA = DELA 127 QQCP = QCP 128 QRMU = RMU 129 QEMU = EMU 130 ENDIF 131 132C KC 2025-05-31 THESE MUST BE RETAINED FOR EXTERNAL PROCEDURE COMPATIBILITY. 133 IF (.FALSE.) PRINT *, GL, LOCOND 134 135 RETURN 136 END

1"""Python version of qenv.f (YREC 7be4133): the right-hand side of the 2envelope integration in log P.""" 3 4import math 5 6import numpy as np 7 8from pyrec import comp_mod 9from pyrec import constants_mod 10from pyrec import dpmix_mod 11from pyrec import envprt_mod 12from pyrec import pulse1_mod 13from pyrec import pulse2_mod 14from pyrec import qenvmix_mod 15from pyrec.eqstat import eqstat 16from pyrec.getopac import getopac 17from pyrec.tpgrad import tpgrad 18 19 20def qenv(x0: float, y: np.ndarray, dydx: np.ndarray, b: float, fpl: float, ftl: float, gl: float, latmo: bool, 21 lderiv: bool, locond: bool, lsave: bool, rl: float, teffl: float, x: float, z: float, kenv: int, 22 ksaha: int) -> tuple[float, int, int]: 23 """Derivatives of log(m/M), log T and log r with respect to log P in 24 the envelope, from hydrostatic equilibrium, mass conservation and the 25 temperature gradient of tpgrad. Also stores the point in envprt_mod 26 (always) and pulse2_mod (when saving). 27 28 Translates: qenv.f#QENV (7be4133) 29 30 Arguments: x0 = log P; y (array of 3, index 0 unused) = log(m/M), 31 log T, log r; the luminosity b, distortion factors fpl, ftl, teffl, 32 x, z; lsave (store the point for output). The mixing length comes from 33 qenvmix_mod.cmixlq, since the argument list is that of the derivative 34 routine called by bsstep and mmid (gl and locond are not used). 35 Arrays changed in place: dydx[1..3]. 36 37 Returns (rl, kenv, ksaha): rl, set to y[3]; the envelope-call 38 counter, incremented; the counter of Saha calls, from eqstat. 39 40 """ 41 42 c = constants_mod 43 fxion = np.zeros(4) # index 0 unused 44 45 # The state, its equation of state, opacity and gradient (Fortran lines 46 # 58–71). 47 pl = x0 48 sl = y[1] + comp_mod.stotal 49 tl = y[2] 50 rl = y[3] 51 (t, p, dl, d, beta, beta1, beta14, rmu, amu, emu, eta, qdt, qdp, qcp, 52 dela, qdtt, qdtp, qat, qap, qcpt, qcpp, ksaha) = eqstat( 53 tl, 0.0, pl, 0.0, 0.0, 0.0, x, z, 0.0, 0.0, 0.0, fxion, 0.0, 0.0, 54 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 55 lderiv, latmo, ksaha) 56 o, ol, qod, qot = getopac(dl, tl, x, z, 0.0, 0.0, 0.0, 0.0, fxion) 57 dpmix_mod.iovim = -1 58 (dele, delr, qact, qacp, qacr, vel, lconv) = tpgrad( 59 tl, t, pl, p, d, rl, sl, b, o, qdt, qdp, qot, qod, qcp, 0.0, 60 0.0, dela, qdtt, qdtp, qat, qap, 0.0, 0.0, 0.0, qcpt, qcpp, 61 0.0, lderiv, False, fpl, ftl, teffl, qenvmix_mod.cmixlq) 62 63 # The derivatives (Fortran lines 72–75). 64 dydx[1] = -math.exp(c.cln * (c.c4pil + 4.0 * rl + pl - c.cgl - sl - sl)) / fpl 65 dydx[2] = dele 66 # Suspected bug (see the page): FPL multiplies here, where COEFFT's form 67 # would divide. 68 dydx[3] = -math.exp(c.cln * (pl + rl - c.cgl - sl - dl)) * fpl 69 kenv = kenv + 1 70 71 # Always store the basic structure variables (07/02) (Fortran lines 72 # 76–94). 73 e = envprt_mod 74 e.ep = pl 75 e.et = tl 76 e.es = sl - comp_mod.stotal 77 e.er = rl 78 e.ed = dl 79 e.evel = vel 80 # Always store the gradients (for tracking the convection zone) (JVS 08/13). 81 e.edel[1] = delr 82 e.edel[2] = dela 83 e.edel[3] = dele 84 e.ebeta = beta 85 # Also always save these, for the profile (JVS 08/25). 86 e.efxion[1] = fxion[1] 87 e.efxion[2] = fxion[2] 88 e.efxion[3] = fxion[3] 89 q = pulse2_mod 90 q.qqdp = qdp 91 q.qqdt = qdt 92 q.qqcp = qcp 93 94 # Store more for output (Fortran lines 96–130). 95 if lsave or pulse1_mod.lpumod: 96 e.eo = o 97 # These repeat the assignments above. 98 e.efxion[1] = fxion[1] 99 e.efxion[2] = fxion[2] 100 e.efxion[3] = fxion[3] 101 q.qtl = tl 102 q.qt = math.exp(c.cln * tl) 103 q.qpl = pl 104 q.qp = math.exp(c.cln * pl) 105 q.qdl = dl 106 q.qd = math.exp(c.cln * dl) 107 q.qo = o 108 q.qol = ol 109 q.qfs = math.exp(c.cln * (sl - comp_mod.stotal)) 110 q.qqdp = qdp 111 q.qqed = 0.0 112 q.qqod = qod 113 q.qqot = qot 114 q.qdel = dele 115 q.qqdt = qdt 116 q.qdela = dela 117 q.qqcp = qcp 118 q.qrmu = rmu 119 q.qemu = emu 120 121 return rl, kenv, ksaha