massloss.f

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

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

Variables on this page (21)

ap, at, c4pi, c4pil, ca3, cgl, clsunl, cmsun, creim, crsun, csecyr, csig, csigl, dmdt0, jcz, json, lmdot, lreimer, sacc, scen, taucz

Purpose

MASSLOSS is the driver for mass loss or accretion in a step (MHP 10/01). It sets the rate $\dot M$ (constant, or from a Reimers formula), computes quantities of the convective envelope and of accreted matter, and calls MDOT to change the model's mass. Sign convention: $\dot M > 0$ is a mass gain (accretion), $\dot M < 0$ a loss. Its header comment: "This is experimental code and valid for Allard atmospheres only" (llp 06/15/2009). KNOWN_ISSUES.md (issue 12) calls this section "experimental, partially implemented": The quantities it computes for physics "not wired up yet" are "not finished," not "dead." The entropies serve accretion: a commented-out line in MDOT used SACC - SCEN, the entropy of the accreted matter relative to the envelope's (whether that difference was meant to decide if accretion expands the envelope is asked on MDOT's page). As the code stands, MDOT (its step 4) scales the expansion or contraction by the envelope's entropy SCEN alone, and SCEN set to 0 at the end of this routine switches the effect off.

How it fits

Who calls it: the main program, main.f, in a step when LAGE is true (main.f:700–705), before the lithium burning of LIBURN; the main program sets LNEWFIT false just before the call.
Inputs: BL ($\log L/L_\odot$), DAGE (Gyr), DELTS (s), the model (HCOMP, HD, HP, HR, HS, HS1, HS2, HT, HJM, OMEGA, M), HSTOT (log of the total mass in g), SMASS (mass in $M_\odot$), TEFFL, JENV (the base of the surface convection zone).
Outputs: LNEWFIT, set true only for accretion; the main program then calls GETNEWENV to restore the earlier fitting point (main.f:1008). Through MDOT: the changed model, DELTS, SMASS, HSTOT, LNEW, SENVOLD. In MASSCHG2_MOD: SACC and SCEN (MDOT reads SCEN; COEFFT reads only DLOGT and DLOGP, which MDOT sets). In DEUTER_MOD: JCZ, the base of the convective envelope (or M if there is none).
Calls: ALSURFP (surface pressure from Allard atmospheres), EQSTAT (twice), MDOT.
Modules: NMLPAR_MOD (the run's parameters, formerly in COMMON MASSCHG: DMDT0, CREIM, LREIMER, LMDOT, COMPACC; set in PARMIN, defaults DMDT0 = −1e−14 $M_\odot$/yr, CREIM = −4e−13, LMDOT and LREIMER false), MASSCHG2_MOD, ATMPRT_MOD (the surface point from ALSURFP), OVRTRN_MOD (TAUCZ, the convective turnover time; the comment at line 104 says it is computed in MIXCZ, but it is now set by GETTAU, through TAUINTNEW (the new method), TAUINT or TAUCAL), DEUTER_MOD, CONSTANTS_MOD (the constants), and DIMENSIONS_MOD (JSON).

The method

If mass change is off (LMDOT false), set LNEWFIT false and return. Otherwise $\dot M$ = DMDT0 (in $M_\odot$/yr, as MDOT converts it).
Radius from $L = 4\pi R^2\sigma T_\mathrm{eff}^4$: $\log R = \tfrac12(\log L - \log 4\pi\sigma - 4\log T_\mathrm{eff})$, with $L$ in erg/s.
Reimers (LREIMER): $\dot M = \eta\,L/(g R)$ with $\eta$ = CREIM and $g = GM/R^2$, all in cgs (bug 1 on the units). The age in seconds, AGETOT, is computed and not used (unused?, line 97 ↗).
If there is a convective envelope (JENV < M), print TAUCZ (in years) and $R/R_\odot$ (labeled MASSLOSS: TAUCZ (YR), R/RSUN), and set JCZ = JENV; otherwise set TAUCZ = 0 (in the shared block) and JCZ = M. The envelope's mass $M_\mathrm{cz}$ = SUMDM $= \sum_{j=\mathrm{JENV}}^{M}$ HS2(j).
If $\dot M > 0$ (accretion), set LNEWFIT true, and with $e_\mathrm{acc} = f\,GM/R$ (EACC0; $f$ = FACC = 1), which the comments describe as the kinetic energy per gram of accreted matter (from the virial theorem):
over the convective envelope, the mass-weighted means of $e_\mathrm{acc}[1 + \tfrac12(R/r - 1)]$ (this mean overwrites EACC), of the thermal energy $P\beta/\rho$ (ETHAV), and of the entropy-like quantity $$s = \frac{P\beta}{\rho T}\left(\tfrac32\ln T - \ln\rho\right),\qquad \beta = 1 - \frac{aT^4}{3P}$$ (SCEN; $a/3$ = CA3; $P\beta/(\rho T)$ = RMU is the gas constant over the mean molecular weight);
at the photosphere from ALSURFP and the EOS: SACC, the same $s$, only if the lookup succeeded (LAlFail false); otherwise a message says that the accreted matter's entropy was not computed, and SACC keeps its earlier value;
an alternative surface state "from Stahler 1988," with $\dot M$ in g/s: $P = \frac{\dot M}{4\pi}\sqrt{2f\,g/R^3}$, $T = \big(\tfrac34\,\dot M\,e_\mathrm{acc}/(4\pi R^2\sigma)\big)^{1/4}$ (with EACC0), its $s$ in SACC2;
finally SCEN is set to 0 in both branches, so the envelope's mean entropy computed above is discarded, and SACC2 is not used (unused?: The lines that used it are commented out, lines 221–224 ↗); the code is unfinished (issue 12 above; settled from YREC's notes, Nov 9, 2026: lines 152–177 ↗, lines 209–227 ↗).
Call MDOT with $\dot M$, the mean $e_\mathrm{acc}$, the mean thermal energy, $M_\mathrm{cz}$, and RMU. MDOT recomputes RMU before using it, and uses EACC and ETHAV only for ETHNEW and FAC, which it does not use (mdot.f, lines 162–174 ↗): The three arguments have no effect (unused?), so for a loss it does not matter that they are stale values from an earlier call (SAVE), or undefined on a first call.

Design

It writes TAUCZ (in years) and $R/R_\odot$ to the terminal on every call with a convective envelope. Its locals are SAVEd.

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. (changes results, large) Checked by KC (Nov 10, 2026) The Reimers rate is in g/s, read as $M_\odot$/yr.

(lines 100–103 ↗; mdot.f, line 111 ↗; mdot.f, line 159 ↗) DMDT = CREIM*10.0D0**(BL + CLSUNL)/GSURF/RTOT is $\eta\,L/(gR) = \eta\,LR/(GM)$ with $L$, $R$, $M$ in cgs, a rate in g/s times $\eta$: for the Sun $L_\odot R_\odot/(GM_\odot) \approx 2.0\times10^{18}$ g/s. MDOT reads DMDT as $M_\odot$/yr (it multiplies by CMSUN/CSECYR). With the default CREIM = −4e−13, Reimers' coefficient for $L$, $R$, $M$ in solar units and $\dot M$ in $M_\odot$/yr, the Sun would lose about $8\times10^{5}\,M_\odot$/yr, and MDOT's step limits would cut the step accordingly. Only with LMDOT and LREIMER (both off by default).

A possible fix: Reimers' law in solar units:

massloss.f

@@ lines 98–105 @@
 C USE A REIMERS FORMULA TO COMPUTE MDOT IF DESIRED; OVERWRITES
 C CONSTANT MDOT.  IN THIS EXPRESSION MDOT=K*L/G/R.
       IF(LDOMDOT .AND.LREIMER)THEN
-         GSURF = 10.0D0**(CGL)*STOT/RTOT**2
-         DMDT = CREIM*10.0D0**(BL + CLSUNL)/GSURF/RTOT
+C 11/09/26 KC via Claude
+C code changed to evaluate the Reimers rate with L, R and M in solar
+C units, for which CREIM is defined, so that DMDT is in MSUN/YR.
+         DMDT = CREIM*10.0D0**BL*(RTOT/CRSUN)/SMASS
       ENDIF
 C 02/12 MHP TAUCZ NOW COMPUTED PRIOR TO CALL IN MIXCZ
 C CONVECTIVE OVERTURN TIMESCALE

How the code itself shows it. DMDT holds one quantity for every option: MASSLOSS starts it at the constant rate, DMDT = DMDT0 (line 80 ↗), whose default is -1.0D-14 (nmlpar_mod.f, line 101 ↗), and the Reimers line "OVERWRITES CONSTANT MDOT" (lines 98–99 ↗). MDOT then converts it under the comment "CONVERT FROM SOLAR MASSES/YEAR TO GM/SEC", RATE = ABS(DMDT)*CMSUN/CSECYR (mdot.f, lines 110–111 ↗), and its commented-out test value, DMDT = 4.0D-10 (mdot.f, line 106 ↗), is a rate in $M_\odot$/yr as well. The Reimers line instead takes $L$, $g$, and $R$ in cgs, as lines 90–95 ↗ compute them ("COMPUTE GLOBAL QUANTITIES (RADIUS,MASS,AGE) IN CGS UNITS."), so its DMDT is in g/s times CREIM. The default CREIM = -4.0D-13 (nmlpar_mod.f, line 106 ↗) is the size of Reimers' coefficient for $L$, $R$, and $M$ in solar units and $\dot M$ in $M_\odot$/yr; the fix evaluates the law so, with the RTOT and SMASS this routine already has.

Checked by running YREC (Nov 9, 2026, the 1 $M_\odot$ runs repeated on Nov 10 with a private start model; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): YREC's cases leave LMDOT off, so two made-up cases were run with LMDOT and LREIMER on, as the cases' namelist comment offers ("Reimers are currently enabled"): testsuite/Test_m3p0feh+0p0GS98_std25_TAMS (3 $M_\odot$, main sequence) and testsuite/Test_m1p0feh+0p0GS98_base_HeIgnite (1 $M_\odot$, to helium ignition). The Reimers line ran 8,341 and 920 times. In the 1 $M_\odot$ run DMDT was $-5.9\times10^{5}$ to $-1.8\times10^{6}$, where Reimers' law in solar units gives $-2.9\times10^{-13}$ to $-9.0\times10^{-13}$ $M_\odot$/yr; in the 3 $M_\odot$ run it reached $-6.7\times10^{9}$, against $-3.3\times10^{-9}$. Unfixed, the 3 $M_\odot$ star fell to 0.11 $M_\odot$ by $5.5\times10^{-5}$ Gyr (8,363 models) and failed to converge; the 1 $M_\odot$ star, starting at the end of the main sequence, fell to 0.912 $M_\odot$ in steps of a few milliseconds (htimer.f, suspected bug 1) and failed to converge, still at 9.864 Gyr (919 models; 893 models, 0.914 $M_\odot$, and the same age with 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}$). With the fix, the 3 $M_\odot$ star keeps its mass (2.99997 $M_\odot$ at model 129, 0.0017 Gyr), and the 1 $M_\odot$ star climbs the giant branch to 12.595 Gyr at 0.707 $M_\odot$; both then stop without a message at mdot.f, line 118 ↗ (mdot.f, suspected bug 2). The effect is real: It changes how the runs end, far beyond the noise. No comment says the cgs form is meant; the default CREIM is Reimers' coefficient for solar units.

Definitions

subroutineMASSLOSS

lines 9–245 · GitHub

Sets $\dot M$ (constant or Reimers), computes the accretion energetics and entropies, and calls MDOT (above).

Uses modules ATMPRT_MOD, CONSTANTS_MOD, DEUTER_MOD, DIMENSIONS_MOD, MASSCHG2_MOD, NMLPAR_MOD, OVRTRN_MOD

Called by (file level)

Flowchart of its calls, in order (4)

if (DMDT.GT.0.0D0), lines 131–228

ALSURFP 182

if (.NOT.LAlFail), lines 191–201

EQSTAT 196

if (FACC.GT.0.0D0), lines 209–224

EQSTAT 215

if (LDOMDOT): MDOT 240

Change log

The changes since rp_1

Changed since rp_1. COMMON blocks to modules, IMPLICIT NONE; the unused SCRTCH dropped; the call to MDOT no longer passes BL, so its 23 arguments line up (settles former question 6); MASSCHG2's order corrected to …, DLOGP, DLOGT, as in MDOT and COEFFT (settles former question 7).

Changed since a8a5950. LAlFail is now checked: after a failed ALSURFP lookup, the surface point (the previous lookup's AT, AP) is no longer used, and the accreted matter's entropy is skipped with a message, which settles that part of question 1; the diagnostic print of TAUCZ and $R/R_\odot$ is labeled, which question 4 now reflects.

1C MHP 10/01 2C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 3C MASS LOSS DRIVER ROUTINE 4C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 5 6C This is experimental code and valid for Allard atmospheres only. 7C llp 06/15/2009 8 9 SUBROUTINE MASSLOSS(BL,DAGE,DELTS,HCOMP,HD,HJM,HP,HR,HS, 10 * * HS1,HS2,HSTOT,HT,JENV,LNEW,M,OMEGA,SMASS, 11 * * TEFFL,SENVOLD,LNEWFIT) 12C 10/15/26 MHP via Claude 13C code changed to take JSON from DIMENSIONS_MOD instead of declaring 14C PARAMETER(JSON=5000) here. 15 USE DIMENSIONS_MOD 16C 10/14/26 MHP via Claude 17C code changed to convert COMMON/OVRTRN/ into OVRTRN_MOD. 18 USE OVRTRN_MOD 19C 10/14/26 MHP via Claude 20C code changed to convert COMMON/MASSCHG2/ into MASSCHG2_MOD. 21 USE MASSCHG2_MOD 22C 10/14/26 MHP via Claude 23C code changed to convert COMMON/DEUTER/ into DEUTER_MOD. 24 USE DEUTER_MOD 25C 10/14/26 MHP via Claude 26C code changed to convert COMMON/ATMPRT/ into ATMPRT_MOD. 27 USE ATMPRT_MOD 28C 10/14/26 MHP via Claude 29C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD. 30 USE NMLPAR_MOD 31C 10/14/26 MHP via Claude 32C code changed to convert COMMON/CONSTANTS/ into CONSTANTS_MOD. 33 USE CONSTANTS_MOD 34C 10/20/26 MHP via Claude 35C code changed to replace implicit typing with IMPLICIT NONE and 36C explicit declarations. 37 IMPLICIT NONE 38 REAL*8 AGETOT,AMU,BETA,BETA14,BETAI,BL,D,DAGE,DELA,DELACC,DELTS, 39 * * DL,DMDT,DMDT2,EACC,EACC0,EMU,ETA,ETH,ETHACC,ETHACCBAR, 40 * * ETHAV,FACC,FXION,GL,GSURF,HCOMP,HD,HJM,HP,HR,HS,HS1,HS2, 41 * * HSTOT,HT,OMEGA,P,PL,QAP,QAT,QCP,QCPP,QCPT,QDP,QDT,QDTP, 42 * * QDTT,R,RL,RMU,RTOT,SACC2,SENVOLD,SLOCAL,SMASS,STOT,SUMDM, 43 * * SUMETH,T,TEFFL,TL,X,Z 44 INTEGER J,JENV,KSAHA,M 45 LOGICAL*4 LALFAIL,LATMO,LDERIV,LDOMDOT,LNEW,LNEWFIT,LPRT 46 DIMENSION HCOMP(15,JSON),HD(JSON),HP(JSON),HS(JSON), 47 * * HS1(JSON),HS2(JSON),HT(JSON),HJM(JSON), 48 * * OMEGA(JSON),HR(JSON),FXION(3) 49C 9/25/26 MHP via Claude 50C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/, 51C and COMMON/CONST3/ into a single COMMON/CONSTANTS/ 52C 10/12/26 MHP via Claude 53C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 54C namelist variables that the physics routines read from 27 COMMON 55C blocks which also held derived or working members. Those blocks keep 56C only their non-namelist members; five of them -- NEWMX, VMULT, 57C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 58C Derivation-only namelist variables from the same blocks went to 59C COMMON/NMLBYPASS/ instead, not here. 60C 61C Declared identically in all 121 files that hold any of these members. 62C Generated from one member list, so every copy is byte-identical; the 63C linker's 'size of symbol changed' warning is the check that they are. 64C 10/11/26 MHP via Claude 65C code changed to fix the member ORDER of COMMON/MASSCHG2/: this file 66C had DLOGT,DLOGP while mdot.f (which writes them) and coefft.f (which 67C reads them) both declare DLOGP,DLOGT -- i.e. the two were 68C transposed here, so DLOGT in this file aliased mdot.f's DLOGP and 69C vice versa. Harmless so far only because this file uses neither 70C (nor SMASS0); it uses only SACC and SCEN, which are at the same 71C positions in both orderings, so this correction is behavior-neutral. 72C 10/13/26 MHP via Claude: removed the unused COMMON/SCRTCH/ declaration 73C (no member referenced here; residue of the obsolete-option removals). 74C G Somers 3/17, ADDING NEW TAUCZ COMMON BLOCK 75C MHP 5/02 EFFICIENCY FACTOR FOR THE THERMAL ENERGY CONTENT 76C OF ACCRETED MATTER. 77 DATA FACC/1.0D0/ 78 SAVE 79C INITIALIZE MASS LOSS AT DEFAULT RATE 80 DMDT = DMDT0 81 IF(LMDOT)THEN 82 LDOMDOT = .TRUE. 83 ELSE 84 LDOMDOT = .FALSE. 85 LNEWFIT = .FALSE. 86 GOTO 9999 87 ENDIF 88C IF(.NOT.LDOMDOT)RETURN 89C TEFFL IS THE BASE 10 LOG OF THE EFFECTIVE TEMPERATURE 90C COMPUTE GLOBAL QUANTITIES (RADIUS,MASS,AGE) IN CGS UNITS. 91C RADIUS 92 RL = 0.5d0*(BL + CLSUNL - C4PIL - CSIGL - 4.0D0*TEFFL) 93 RTOT = 10.0D0**RL 94C MASS 95 STOT = SMASS*CMSUN 96C AGE 97 AGETOT = DAGE*1.0D9*CSECYR 98C USE A REIMERS FORMULA TO COMPUTE MDOT IF DESIRED; OVERWRITES 99C CONSTANT MDOT. IN THIS EXPRESSION MDOT=K*L/G/R. 100 IF(LDOMDOT .AND.LREIMER)THEN 101 GSURF = 10.0D0**(CGL)*STOT/RTOT**2 102 DMDT = CREIM*10.0D0**(BL + CLSUNL)/GSURF/RTOT 103 ENDIF 104C 02/12 MHP TAUCZ NOW COMPUTED PRIOR TO CALL IN MIXCZ 105C CONVECTIVE OVERTURN TIMESCALE 106 IF(JENV.LT.M)THEN 107C TAUCZ = 0.0D0 108C DO I = JENV+1,M 109C V = 0.5D0*(SVEL(I-1)+SVEL(I)) 110C DR = 10.0D0**(HR(I))-10.0D0**(HR(I-1)) 111C IF(V.GT.0.0D0)THEN 112C TAUCZ = TAUCZ + DR/V 113C ENDIF 114C END DO 115C 10/29/26 MHP via Claude: label this diagnostic, which was unlabelled. 116 WRITE(*,'(A,1P2E12.4)')' MASSLOSS: TAUCZ (YR), R/RSUN ', 117 * * TAUCZ/CSECYR,RTOT/CRSUN 118 JCZ = JENV 119 ELSE 120 TAUCZ = 0.0D0 121 JCZ = M 122 ENDIF 123C MHP 8/10 124C THE RUNNING TOTAL OF THE MASS OF THE 125C CONVECTIVE ENVELOPE, NEEDED ELSEWHERE AS WELL 126 SUMDM = 0.0D0 127 DO J = JENV, M 128 SUMDM = SUMDM + HS2(J) 129 END DO 130C COMPUTE THE MEAN MOLECULAR WEIGHT FOR THE CZ 131 IF(DMDT.GT.0.0D0)THEN 132 LNEWFIT = .TRUE. 133 T = 10.0D0**HT(JENV) 134 P = 10.0D0**HP(JENV) 135 D = 10.0D0**HD(JENV) 136 BETA = 1.0D0- (CA3*T**4/P) 137 RMU = P*BETA/(D*T) 138C FROM THE VIRIAL THEOREM, A MAXIMUM OF 139C 1/2 OF THE GRAVITATIONAL POTENTIAL ENERGY OF 140C ACCRETED MATTER WILL BE CONVERTED INTO KINETIC 141C ENERGY. THE PARAMETER FACC IS USED TO INFER 142C THE KINETIC ENERGY PER GRAM WHICH IS COMPARED 143C WITH THE MEAN THERMAL ENERGY CONTENT OF THE 144C CONVECTION ZONE. 145 EACC = FACC*STOT*(10.0D0**CGL)/RTOT 146 EACC0 = EACC 147C DETERMINE THE MASS-WEIGHTED THERMAL ENERGY (PER GM) 148C IN EACH SHELL OF THE CONVECTIVE ENV. 149 SUMETH = 0.0D0 150C SUMDM = 0.0D0 151 ETHACCBAR = 0.0D0 152 SCEN = 0.0D0 153 DO J = JENV, M 154 T = 10.0D0**HT(J) 155 P = 10.0D0**HP(J) 156 D = 10.0D0**HD(J) 157 R = 10.0D0**HR(J) 158 DELACC = 0.5D0*EACC*(RTOT/R - 1.0D0) 159 ETHACC = EACC + DELACC 160 ETHACCBAR = ETHACCBAR+HS2(J)*ETHACC 161 BETA = 1.0D0- (CA3*T**4/P) 162 RMU = P*BETA/(D*T) 163 SLOCAL = RMU*(1.5D0*LOG(T)-LOG(D)) 164 SCEN = SCEN+SLOCAL*HS2(J) 165C THE THERMAL ENERGY PER GM IN THE JTH SHELL IS 166 ETH = P*BETA/D 167C THE RUNNING TOTAL OF THERMAL ENERGY THROUGHOUT THE 168C CONVECTIVE ENVELOPE IS 169 SUMETH = SUMETH + ETH*HS2(J) 170C THE RUNNING TOTAL OF THE MASS OF THE 171C CONVECTIVE ENVELOPE IS 172C SUMDM = SUMDM + HS2(J) 173 END DO 174C THE AVERAGE THERMAL ENERGY OF THE UNPERTURBED CE (W/O ACCN) IS 175 ETHAV = SUMETH/SUMDM 176 EACC = ETHACCBAR/SUMDM 177 SCEN = SCEN/SUMDM 178 LPRT = .FALSE. 179 GL = CGL + HSTOT - 2.0D0*RL 180C This is experimental code and valid for Allard atmospheres only. 181C llp 06/15/2009 182 CALL ALSURFP(TEFFL,GL,LPRT,LAlFail) 183 X = HCOMP(1,M) 184 Z = HCOMP(3,M) 185 LDERIV = .FALSE. 186 KSAHA = 1 187 LATMO = .TRUE. 188C 10/29/26 MHP via Claude 189C Use the Allard surface point only if the lookup succeeded; AT and AP 190C otherwise still hold the previous lookup's values. 191 IF(.NOT.LAlFail)THEN 192 T = 10.0D0**AT 193 P = 10.0D0**AP 194 TL = AT 195 PL = AP 196 CALL EQSTAT(TL,T,PL,P,DL,D,X,Z,BETA,BETAI,BETA14,FXION,RMU, 197 * * AMU,EMU,ETA,QDT,QDP,QCP,DELA,QDTT,QDTP,QAT,QAP,QCPT, 198 * * QCPP,LDERIV,LATMO,KSAHA) 199 BETA = 1.0D0- (CA3*T**4/P) 200 RMU = P*BETA/(D*T) 201 SACC = RMU*(1.5D0*LOG(T)-LOG(D)) 202 ELSE 203 WRITE(*,*)'MASSLOSS: Allard surface lookup failed; ', 204 * * 'accreted-matter entropy not computed' 205 ENDIF 206C WRITE(*,911)TL,PL,SACC,SCEN 207C 911 FORMAT(' TSUR,PSUR ',2F8.5,' SACC ',1PE12.3,' SCORE ',E12.3) 208C ALTERNATE EXPRESSION FOR SURFACE PRESSURE AND LUMINOSITY, FROM STAHLER 1988 209 IF(FACC.GT.0.0D0)THEN 210 DMDT2 = DMDT*CMSUN/CSECYR 211 P = DMDT2/C4PI*SQRT(2.0D0*FACC*10.0D0**GL/RTOT**3) 212 T = (DMDT2*EACC0*0.75D0/C4PI/RTOT**2/CSIG)**0.25D0 213 PL = LOG10(P) 214 TL = LOG10(T) 215 CALL EQSTAT(TL,T,PL,P,DL,D,X,Z,BETA,BETAI,BETA14,FXION,RMU, 216 * * AMU,EMU,ETA,QDT,QDP,QCP,DELA,QDTT,QDTP,QAT,QAP,QCPT, 217 * * QCPP,LDERIV,LATMO,KSAHA) 218 BETA = 1.0D0- (CA3*T**4/P) 219 RMU = P*BETA/(D*T) 220 SACC2 = RMU*(1.5D0*LOG(T)-LOG(D)) 221C WRITE(*,911)TL,PL,SACC2,SCEN 222C SACC = MAX(SACC,SACC2) 223 SCEN = 0.0D0 224C SCEN = SACC2 - SACC 225 ELSE 226 SCEN = 0.0D0 227 ENDIF 228 ENDIF 229C CALL MASS LOSS OR ACCRETION ROUTINE 230C 10/13/26 MHP via Claude 231C BUG FIX: this call passed BL as a 24th actual argument to MDOT, which 232C declares only 23 dummies -- BL landed in DELTS, DELTS in HCOMP and so 233C on, every argument shifted by one, with SENVOLD having no dummy at all. 234C MDOT's body never references BL, so it was pruned from the dummy list 235C (as other routines were in the KC 2025-05-31 pass) without this call 236C site being updated. Dropping BL realigns all 23. F77 does no 237C cross-file interface checking, so nothing flagged it; it has never 238C misfired only because every namelist sets LMDOT = .FALSE., which makes 239C MASSLOSS return before reaching this line. 240 IF(LDOMDOT)CALL MDOT(DELTS,HCOMP,HD,HJM,HP,HR,HS,HS1,HS2, 241 * * HSTOT,HT,JENV,LNEW,M,OMEGA,RMU,RTOT,SMASS,DMDT, 242 * * EACC,ETHAV,SUMDM,SENVOLD) 243 9999 9999 CONTINUE 244 RETURN 245 END

1"""Python version of massloss.f (YREC 7be4133): the driver for mass loss or 2accretion in a step.""" 3 4import math 5import numpy as np 6 7from pyrec import atmprt_mod 8from pyrec import constants_mod 9from pyrec import deuter_mod 10from pyrec import masschg2_mod 11from pyrec import nmlpar_mod 12from pyrec import ovrtrn_mod 13from pyrec.alsurfp import alsurfp 14from pyrec.eqstat import eqstat 15from pyrec.mdot import mdot 16 17# Efficiency factor for the thermal energy content of accreted matter (DATA 18# FACC/1.0D0/) (MHP 5/02). 19FACC = 1.0 20 21# SAVEd locals passed to MDOT that are set only for accretion: for a loss 22# MDOT gets the values of an earlier call (undefined on the first; 0.0 here). 23# MDOT does not use them (see the page). 24EACC = 0.0 25ETHAV = 0.0 26RMU = 0.0 27 28 29def massloss(bl: float, dage: float, delts: float, hcomp: np.ndarray, hd: np.ndarray, hjm: np.ndarray, hp: np.ndarray, 30 hr: np.ndarray, hs: np.ndarray, hs1: np.ndarray, hs2: np.ndarray, hstot: float, ht: np.ndarray, jenv: int, 31 lnew: bool, m: int, omega: np.ndarray, smass: float, teffl: float, senvold: float, 32 lnewfit: bool) -> tuple[float, float, bool, float, float, bool]: 33 """Set the mass-change rate (constant, or a Reimers formula), compute 34 the convective envelope's mass and, for accretion, the energetics and 35 entropies of the envelope and of accreted matter, then call MDOT to 36 change the model's mass. dmdt > 0 is accretion, < 0 a loss. 37 Experimental, and valid for Allard atmospheres only (llp 06/15/2009). 38 39 Translates: massloss.f#MASSLOSS (7be4133) 40 41 Arguments: bl (log L/Lsun), dage (Gyr), delts (s); the model hcomp, 42 hd, hjm, hp, hr, hs, hs1, hs2, ht, omega (arrays, index 0 unused, 43 changed in place by MDOT); hstot (log total mass, g); jenv (base of 44 the surface convection zone); lnew; m (number of shells); smass 45 (Msun); teffl (log Teff); senvold; lnewfit is output only. 46 47 Returns (delts, hstot, lnew, smass, senvold, lnewfit): the first four 48 and senvold as changed by MDOT (unchanged if LMDOT is false), and 49 lnewfit, True only for accretion (the caller then restores the 50 envelope fitting point). 51 52 """ 53 54 global EACC, ETHAV, RMU 55 c = constants_mod 56 nml = nmlpar_mod 57 58 # Initialize mass loss at the default rate (Fortran lines 80–87). 59 dmdt = nml.dmdt0 60 if nml.lmdot: 61 ldomdot = True 62 else: 63 ldomdot = False 64 lnewfit = False 65 return delts, hstot, lnew, smass, senvold, lnewfit # label 9999 66 67 # Global quantities in cgs: radius, mass, age (Fortran lines 92–97). 68 rl = 0.5 * (bl + c.clsunl - c.c4pil - c.csigl - 4.0 * teffl) 69 rtot = 10.0**rl 70 stot = smass * c.cmsun 71 agetot = dage * 1.0e9 * c.csecyr # noqa: F841 (unused in the Fortran too) 72 73 # Reimers formula, MDOT = K*L/g/R (overwrites the constant rate) 74 # (Fortran lines 100–103). 75 if ldomdot and nml.lreimer: 76 gsurf = 10.0**c.cgl * stot / rtot**2 77 # Suspected bug (see the page): cgs, but MDOT reads Msun/yr. 78 dmdt = nml.creim * 10.0**(bl + c.clsunl) / gsurf / rtot 79 80 # Convective overturn time (computed earlier, in MIXCZ's call) (Fortran 81 # lines 106–122). 82 if jenv < m: 83 # Question (see the page): a debug line. 84 print(f' MASSLOSS: TAUCZ (YR), R/RSUN {ovrtrn_mod.taucz / c.csecyr:12.4E}' 85 f'{rtot / c.crsun:12.4E}') 86 deuter_mod.jcz = jenv 87 else: 88 ovrtrn_mod.taucz = 0.0 89 deuter_mod.jcz = m 90 91 # The mass of the convective envelope (MHP 8/10) (Fortran lines 92 # 126–129). 93 sumdm = 0.0 94 for j in range(jenv, m + 1): 95 sumdm = sumdm + hs2[j] 96 97 # Accretion (Fortran lines 131–228). 98 if dmdt > 0.0: 99 lnewfit = True 100 t = 10.0**ht[jenv] 101 p = 10.0**hp[jenv] 102 d = 10.0**hd[jenv] 103 beta = 1.0 - (c.ca3 * t**4 / p) 104 RMU = p * beta / (d * t) 105 106 # Kinetic energy per gram of accreted matter (virial theorem) 107 # (Fortran lines 145–146). 108 EACC = FACC * stot * (10.0**c.cgl) / rtot 109 eacc0 = EACC 110 111 # Mass-weighted means over the convective envelope (Fortran lines 112 # 149–177). 113 sumeth = 0.0 114 ethaccbar = 0.0 115 masschg2_mod.scen = 0.0 116 for j in range(jenv, m + 1): 117 t = 10.0**ht[j] 118 p = 10.0**hp[j] 119 d = 10.0**hd[j] 120 r = 10.0**hr[j] 121 delacc = 0.5 * EACC * (rtot / r - 1.0) 122 ethacc = EACC + delacc 123 ethaccbar = ethaccbar + hs2[j] * ethacc 124 beta = 1.0 - (c.ca3 * t**4 / p) 125 RMU = p * beta / (d * t) 126 slocal = RMU * (1.5 * math.log(t) - math.log(d)) 127 masschg2_mod.scen = masschg2_mod.scen + slocal * hs2[j] 128 eth = p * beta / d 129 sumeth = sumeth + eth * hs2[j] 130 ETHAV = sumeth / sumdm 131 EACC = ethaccbar / sumdm 132 masschg2_mod.scen = masschg2_mod.scen / sumdm 133 134 # The photosphere from the Allard tables, and its entropy (Fortran 135 # lines 178–205). 136 lprt = False 137 gl = c.cgl + hstot - 2.0 * rl 138 # Experimental: valid for Allard atmospheres only. 139 # Question (see the page): Allard atmospheres only. 140 lalfail = alsurfp(teffl, gl, lprt, False) 141 x = hcomp[1, m] 142 z = hcomp[3, m] 143 fxion = [0.0] * 4 144 lderiv = False 145 ksaha = 1 146 latmo = True 147 # Use the Allard surface point only if the lookup succeeded; AT and 148 # AP otherwise still hold the previous lookup's values. 149 if not lalfail: 150 t = 10.0**atmprt_mod.at 151 p = 10.0**atmprt_mod.ap 152 tl = atmprt_mod.at 153 pl = atmprt_mod.ap 154 (t, p, dl, d, beta, betai, beta14, RMU, amu, emu, eta, qdt, qdp, 155 qcp, dela, qdtt, qdtp, qat, qap, qcpt, qcpp, ksaha) = eqstat( 156 tl, t, pl, p, 0.0, 0.0, x, z, 0.0, 0.0, 0.0, fxion, 0.0, 0.0, 157 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 158 lderiv, latmo, ksaha) 159 beta = 1.0 - (c.ca3 * t**4 / p) 160 RMU = p * beta / (d * t) 161 masschg2_mod.sacc = RMU * (1.5 * math.log(t) - math.log(d)) 162 else: 163 print('MASSLOSS: Allard surface lookup failed; ' 164 'accreted-matter entropy not computed') 165 # EQSTAT's outputs not yet set in this call, passed in below 166 # (any value; Fortran keeps SAVEd values from an earlier call). 167 dl = betai = beta14 = amu = emu = eta = 0.0 168 qdt = qdp = qcp = dela = qdtt = qdtp = qat = qap = 0.0 169 qcpt = qcpp = 0.0 170 171 # Alternate surface pressure and temperature, from Stahler 1988 172 # (Fortran lines 209–227). 173 if FACC > 0.0: 174 dmdt2 = dmdt * c.cmsun / c.csecyr 175 p = dmdt2 / c.c4pi * math.sqrt(2.0 * FACC * 10.0**gl / rtot**3) 176 t = (dmdt2 * eacc0 * 0.75 / c.c4pi / rtot**2 / c.csig)**0.25 177 pl = math.log10(p) 178 tl = math.log10(t) 179 (t, p, dl, d, beta, betai, beta14, RMU, amu, emu, eta, qdt, qdp, 180 qcp, dela, qdtt, qdtp, qat, qap, qcpt, qcpp, ksaha) = eqstat( 181 tl, t, pl, p, dl, d, x, z, beta, betai, beta14, fxion, RMU, 182 amu, emu, eta, qdt, qdp, qcp, dela, qdtt, qdtp, qat, qap, 183 qcpt, qcpp, lderiv, latmo, ksaha) 184 beta = 1.0 - (c.ca3 * t**4 / p) 185 RMU = p * beta / (d * t) 186 # SACC2 is not used and SCEN is zeroed: unfinished code (see the 187 # page). 188 sacc2 = RMU * (1.5 * math.log(t) - math.log(d)) # noqa: F841 (unused in the Fortran too) 189 masschg2_mod.scen = 0.0 190 else: 191 masschg2_mod.scen = 0.0 192 193 # Mass loss or accretion (MDOT, with its 23 arguments) (Fortran line 194 # 240). 195 if ldomdot: 196 (delts, hstot, lnew, RMU, smass, dmdt, senvold) = mdot( 197 delts, hcomp, hd, hjm, hp, hr, hs, hs1, hs2, hstot, ht, jenv, 198 lnew, m, omega, RMU, rtot, smass, dmdt, EACC, ETHAV, sumdm, 199 senvold) 200 return delts, hstot, lnew, smass, senvold, lnewfit # label 9999