mdot.f

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

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

Variables on this page (24)

ca3, cc13, cc23, cln, cmsun, compacc, csecyr, dlogp, dlogt, fczdmdt, fmassacc, ftotdmdt, json, ldisk, lmdot, lrot, lsolwind, sage, scen, senv, smass0, stotal, tdisk, walpcz

Purpose

MDOT applies a mass change (loss or accretion) to the model over a step: it limits the step so that not too much mass changes, removes mass from the surface convection zone or adds it at the surface, rescales the mass coordinates, the total mass and the shell masses, removes angular momentum with lost mass, and for accretion mixes the accreted composition into the convection zone. Sign convention: DMDT > 0 is a mass gain.

How it fits

Who calls it: only MASSLOSS, when LMDOT is true.
Inputs and outputs (arguments): DELTS (s, may be shortened), the model (HCOMP, HD, HJM, HP, HR, HS, HS1, HS2, HT, OMEGA, M, JENV), HSTOT and SMASS (updated), DMDT ($M_\odot$/yr), RTOT, SUMDM (the convection zone's mass, from MASSLOSS); EACC, ETHAV and RMU have no effect (unused?; see step 4); outputs LNEW (= true: recompute the surface boundary condition) and SENVOLD (the old HS(M) − HSTOT, for GETNEWENV).
Modules: NMLPAR_MOD, the namelist variables (FCZDMDT, FTOTDMDT, COMPACC, LMDOT, TDISK, LDISK, LROT, WALPCZ, LSOLWIND); MASSCHG2_MOD (SCEN, SMASS0, DLOGP, DLOGT; with MASSLOSS and COEFFT, which adds DLOGT, DLOGP to $\log T$, $\log P$ at and above JENV in its entropy terms when LMDOT and DMDT0 > 0); DISK_MOD (the age SAGE, for the disk lifetime); COMP_MOD (STOTAL, SENV); DEUTER_MOD (FMASSACC, read by LIBURN, DBURN and DBURNM); CONSTANTS_MOD; JSON from DIMENSIONS_MOD. The solar-wind parameters DMSUN, DMWSUN, DMWMAX (formerly COMMON MASSCHG3) are now locals set to zero by DATA.

The method

Two masses of the convection zone appear: $M_\mathrm{cz}$ = DMCZ = HS1(M) − HS1(JENV) (or HS1(M) if JENV = 1), and $M'_\mathrm{cz}$ = DMCZT = SUMDM.

Optionally (LSOLWIND with rotation, and only if DMWSUN > 0) $\dot M = \dot M_\odot\,\min\big((\Omega/\Omega_\odot)^2, (\Omega_\mathrm{max}/\Omega_\odot)^2\big)$, with $\dot M_\odot$ = DMSUN, $\Omega_\odot$ = DMWSUN, $\Omega_\mathrm{max}$ = DMWMAX. Since these three are zero, the branch is inert; the guard on DMWSUN keeps it from dividing by zero now that LSOLWIND reaches this routine from the namelist.
Rate in g/s: $|\dot M|\,M_\odot/\mathrm{yr}$. Step limits: $\Delta t \le f_\mathrm{tot}M/|\dot M|$, $\Delta t \le f_\mathrm{cz}\,M_\mathrm{cz}/|\dot M|$ ($f_\mathrm{tot}$ = FTOTDMDT, $f_\mathrm{cz}$ = FCZDMDT; a zero rate sets neither), and at most half the mass outside the fitting point, $\tfrac12(10^{\texttt{HSTOT}} - 10^{\texttt{HS(M)}})/|\dot M|$. The comment speaks of adding mass, but this limit applies to loss as well. With a disk (LDISK) and accretion, the step ends at the disk's lifetime TDISK (Gyr), and LMDOT is then turned off (bug 1). A model without a surface convection zone (JENV = M) stops the run (bug 2).
$\Delta M = \dot M\,\Delta t$ (g).
Accretion ($\Delta M > 0$):
the new mean thermal energy $\bar e = (e_\mathrm{acc}\Delta M + \bar e_\mathrm{th}M'_\mathrm{cz})/(\Delta M + M'_\mathrm{cz})$ and FAC = $\log_{10}(\bar e/\bar e_\mathrm{th})$ are computed but not used (their use is commented out);
RMU = $P\beta/(\rho T)$ (the gas constant over $\mu$) is recomputed at JENV, and a scale factor $\Delta\ln E = (\Delta M/M_\mathrm{cz})\,s_\mathrm{cen}/(\mathcal{R}/\mu)$, with $s_\mathrm{cen}$ = SCEN, gives $\Delta\log r = \tfrac23\Delta\ln E/\ln 10$, $\Delta\log T = -\Delta\log r$, $\Delta\log P = -4\Delta\log r$ (homologous expansion of the envelope: $P \propto r^{-4}$, $T \propto r^{-1}$ at fixed $m$; the mass term DLNM is set to zero);
if JENV = 1 (fully convective), $\log r$, $\log P$, $\log T$ of every point change by these amounts; otherwise only $P$ and $T$ at JENV, and the radii above it are stretched about $r_\mathrm{b} = r_\mathrm{JENV}$: $r \leftarrow r_\mathrm{b} + 10^{\Delta\log r}(r - r_\mathrm{b})$.

MASSLOSS sets SCEN = 0, so these changes are zero in practice: The code is kept for later work, as part of the experimental mass-loss section, "not finished," not "dead" (settled from KNOWN_ISSUES.md, issue 12, Nov 9, 2026).

> **Question** ([lines 173–174](#L173-L174)<a class='gh' href='https://github.com/yreclab/yrec/blob/7be4133/src/mdot.f#L173-L174' title='these lines on GitHub'>↗</a>). **The form of the expansion.** A commented-out line scaled the expansion by `SACC - SCEN`, the entropy of the accreted matter relative to the envelope's; the active one by `SCEN` alone. An expansion scaled by `SCEN` alone depends on the units, since the entropy is defined only up to an additive constant (it contains $\ln T$ and $\ln\rho$ in cgs). When the code is revived, is `SACC - SCEN` the intended form?

Angular momentum (rotation, loss): the lost matter carries the surface's specific angular momentum, $\Delta J = \tfrac23 R^2\,\Omega_s\,\Delta M$ (solid-body zone, WALPCZ ≥ 0), or $j_s\Delta M$ (WALPCZ ≤ −2), or the same with $\Omega_s = \Omega_M (R/r_M)^{\texttt{WALPCZ}}$, the law $\Omega \propto r^{\texttt{WALPCZ}}$ that WCZIMP imposes; the zone's new angular momentum is $J_\mathrm{cz}+\Delta J$, with $J_\mathrm{cz} = \sum$ HJM·HS2 from JENV to M. Its specific angular momentum HJM is scaled after step 7, once the new shell masses are known, by $(J_\mathrm{cz}+\Delta J)/\sum$ HJM·HS2, so that the zone holds exactly $J_\mathrm{cz}+\Delta J$.

Question 1 (line 206 ↗; line 280 ↗). Angular momentum of accreted matter. For accretion HJM is not changed while HS2(M) grows, so the accreted matter implicitly takes the outer point's specific angular momentum. Is that intended?
Mass coordinates (loss): points in the zone are compressed toward its base, $m_i \leftarrow m_\mathrm{b} + \frac{M_\mathrm{cz}+\Delta M}{M_\mathrm{cz}}(m_i - m_\mathrm{b})$; in a fully convective star (JENV = 1), toward the center, $m_i \leftarrow \frac{M_\mathrm{cz}+\Delta M}{M_\mathrm{cz}}\,m_i$ for every point, consistent with $M_\mathrm{cz}$ = HS1(M) there. For accretion HS1 is not changed: the gain goes to the outermost shell.

Question 2 (lines 240–258 ↗; line 280 ↗). All accreted mass in the last shell. For accretion HS1 is unchanged and the gain goes to HS2(M), which over many steps grows unless rezoning redistributes it. The main program calls HPOINT after every step (except in helium-flash runs, LKUTHE). Does its criterion split a growing outermost shell, so that the surface resolution survives long accretion?
New total mass (SMASS, SMASS0, HSTOT, STOTAL), shell masses HS2 from JENV to M − 1 as half the distance between neighbours (for JENV = 1, HS2(1) = (HS1(1) + HS1(2))/2, the innermost shell running from the center), HS2(M) = $10^{\texttt{HSTOT}}$ − (HS1(M) + HS1(M−1))/2, and SENV.
Accretion: elements 1–11 of the zone are mixed with the accreted composition COMPACC, $\bar x = \big(\sum_j x_j\,\Delta m_j + x_\mathrm{acc}\Delta M\big)/\big(\sum_j \Delta m_j+\Delta M\big)$, summed over the zone's shells JENV to M with their masses before the accretion ($\Delta m_j$ = HS2(J), less $\Delta M$ for HS2(M)), so that every species is conserved; FMASSACC = $\Delta M$.

Design

Several unconditional WRITE(*,…) lines report the step reduction, the angular momentum and the mass change.

From YREC's documentation

Physics Options lists LMDOT, "experimental mass-loss module," and LSOLWIND, mass loss scaled to the solar wind, "not implemented," among the experimental options.

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) Accretion past the disk's lifetime gives a negative step.

(lines 146–157 ↗; main.f, line 981 ↗) With LDISK and accretion, TEST = SAGE + 1.0d-9*DELTS/CSECYR is the age at the end of the step, and when it exceeds TDISK the step is cut to DELTS = (TDISK-SAGE)*1.0D9*CSECYR. If SAGE already exceeds TDISK, this is negative: DELM = DMDT*CMSUN*DELTS/CSECYR is then negative, so the accretion is applied as a loss, and the main program goes on with a negative step (what it then does was not traced). The main program sets SAGE = DAGE only after the step's call to MASSLOSS (main.f, line 981 ↗), so in the first step SAGE is unset (0 with zero-initialized storage), and from the second step on it is the model's age: a run that starts from a model older than TDISK reaches this. Only with LMDOT, LDISK (both off by default), and DMDT0 > 0.

A possible fix:

mdot.f

@@ lines 146–152 @@
 C mhp 8/10 turn mass loss off when disk exhausted only when dm /dt > 0, e.g. accretion
       IF(LDISK .AND. DMDT.GT.0.0D0)THEN
          TEST = SAGE + 1.0d-9*DELTS/CSECYR
-         IF(TEST.GT.TDISK)THEN
+C 11/09/26 KC via Claude
+C code changed to stop accreting, without a negative step, when the
+C disk is already gone at the start of the step.
+         IF(SAGE.GE.TDISK)THEN
+            DMDT = 0.0D0
+            LDONE = .TRUE.
+         ELSE IF(TEST.GT.TDISK)THEN
             DELTS = (TDISK-SAGE)*1.0D9*CSECYR
             LDONE = .TRUE.
          ELSE

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: No run had LMDOT with LDISK and accretion; the two made-up Reimers runs (massloss.f, suspected bug 1) lose mass, so they skip this test. The fix was not run. What would reach it: LMDOT, LDISK, and DMDT0 > 0, from a starting model older than TDISK. Intent: The comment above the test, "turn mass loss off when disk exhausted only when dm /dt > 0, e.g. accretion," means to end the accretion, not to reverse it.

2. (stops the run) Checked by KC (Nov 10, 2026) A stop without a message when there is no surface convection zone.

(line 118 ↗; massloss.f, lines 106–122 ↗) IF(JENV.EQ.M) STOP911 (that is, STOP 911) ends the run with only the stop code 911, written nowhere else, when the model has no surface convection zone. MASSLOSS handles that case (it sets TAUCZ = 0 and JCZ = M) and then calls MDOT anyway, so a run with mass change on reaches the stop as soon as its envelope turns radiative. Only with LMDOT (off by default).

A possible fix: say why the run stops (whether mass change should instead go on without a convection zone is the developers' choice):

mdot.f

@@ lines 115–121 @@
 C      DMCZT = 10.0D0**HSTOT - HS1(JENV)
       DMCZT = SUMDM
 C COMPUTE MASS OF SURFACE CONVECTION ZONE BELOW FITTING POINT
-      IF(JENV.EQ.M) STOP911
+C 11/09/26 KC via Claude
+C code changed to say why the run stops when there is no surface
+C convection zone.
+      IF(JENV.EQ.M) THEN
+         WRITE(*,*)'MDOT: NO SURFACE CONVECTION ZONE; STOPPED'
+         STOP 911
+      ENDIF
       IF(JENV.GT.1) THEN
          DMCZ = HS1(M) - HS1(JENV)
       ELSE

How the code itself shows it. MASSLOSS expects the case: When JENV equals M, it sets TAUCZ = 0.0D0 and JCZ = M (massloss.f, lines 106–122 ↗) rather than stopping, and MDOT's next lines handle the other end, a fully convective star, with DMCZ = HS1(M) (lines 119–123 ↗), but not this one. Elsewhere YREC meets a missing surface convection zone with a message: MIX writes "NO SURFACE CZ - DIFFUSION NOT MEANINGFUL" and goes on ("DBG 2/92 CHANGED STOP TO JUST A WARNING MESSAGE, EXECUTION CONTINUES", mix.f, lines 428–433 ↗), and ROTMIX writes the same message before it stops with STOP 1 (rotmix.f, lines 163–170 ↗). Other run-time stops, such as CHECKJ's and WIND's, also write why first and use STOP 1, "the run-time error code (issue 38)" (checkj.f, lines 152–159 ↗; wind.f, lines 72–76 ↗); 911 is used as a stop code nowhere else in active code. The fix above keeps STOP 911 and leaves two questions for the developers: Should this stop, too, use STOP 1, or should mass change go on without a convection zone, as diffusion does in MIX?

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): The line ran at every MDOT call of the two made-up Reimers runs (massloss.f, suspected bug 1), 8,341 and 920 times (JENV from 1 to 685, M from 368 to 957), and JENV = M never held: The wrong Reimers rate ruins those runs first. With the Reimers rate fixed, both end here, at STOP 911 with no message (exit status 143, that is 911 mod 256): the 3 $M_\odot$ star at model 129 (0.0017 Gyr), the 1 $M_\odot$ star at model 6,599 (12.595 Gyr, 0.707 $M_\odot$). So a run with mass change that reaches this test ends with only the stop code. No fix was run, since the page's fix adds only a message. Intent: MDOT's header says the routine has never run, since every namelist sets LMDOT = .FALSE.

Definitions

subroutineMDOT

lines 5–326 · GitHub

Step limit, mass removal or addition, angular-momentum loss, and mixing of accreted matter (above).

Uses modules COMP_MOD, CONSTANTS_MOD, DEUTER_MOD, DIMENSIONS_MOD, DISK_MOD, MASSCHG2_MOD, NMLPAR_MOD

Called by MASSLOSS

Change log

The changes since rp_1

Changed since rp_1. COMMON blocks to modules, IMPLICIT NONE; MASSCHG3 dropped (its three reals now zeroed locals, LSOLWIND from the namelist); the solar-wind branch guarded by DMWSUN > 0 (a bug fix: it would have given $\dot M$ = NaN). The call from MASSLOSS now matches the interface (fixed in massloss.f), which settles former question 1; former question 2 (MASSCHG3 set nowhere) is settled by the authors' fix, which confirms it and keeps the option inert until the three parameters are given a source.

Changed since a8a5950. Five fixes in the experimental mass-change path: a zero rate no longer divides by zero in the step limits; the wind's surface rate now follows $\Omega \propto r^{\texttt{WALPCZ}}$, the law WCZIMP imposes (the old form had the ratio inverted), which settles former question 5; HJM is now scaled after the new shell masses are known, so that the zone loses exactly $\Delta J$, which settles former question 6; in a fully convective star the mass loss compresses the points about the center, and HS2(1) is (HS1(1) + HS1(2))/2 instead of reading HS1(0), which settles former question 4; accretion mixes the zone's mass-weighted abundances instead of the base shell's. Former question 7 is now 4.

1C 2C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 3C MASS LOSS ROUTINE 4C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 5 SUBROUTINE MDOT(DELTS,HCOMP,HD,HJM,HP,HR,HS,HS1,HS2,HSTOT, 6 * * HT,JENV,LNEW,M,OMEGA,RMU,RTOT,SMASS,DMDT, 7 * * EACC,ETHAV,SUMDM,SENVOLD) 8C 10/15/26 MHP via Claude 9C code changed to take JSON from DIMENSIONS_MOD instead of declaring 10C PARAMETER(JSON=5000) here. 11 USE DIMENSIONS_MOD 12C 10/14/26 MHP via Claude 13C code changed to convert COMMON/MASSCHG2/ into MASSCHG2_MOD. 14 USE MASSCHG2_MOD 15C 10/14/26 MHP via Claude 16C code changed to convert COMMON/DISK/ into DISK_MOD. 17 USE DISK_MOD 18C 10/14/26 MHP via Claude 19C code changed to convert COMMON/DEUTER/ into DEUTER_MOD. 20 USE DEUTER_MOD 21C 10/14/26 MHP via Claude 22C code changed to convert COMMON/COMP/ into COMP_MOD. 23 USE COMP_MOD 24C 10/14/26 MHP via Claude 25C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD. 26 USE NMLPAR_MOD 27C 10/14/26 MHP via Claude 28C code changed to convert COMMON/CONSTANTS/ into CONSTANTS_MOD. 29 USE CONSTANTS_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 ABUND,BETA,CJNEW,CJTOT,D,DELJ,DELM,DELTS,DLNM,DLOGEN,DLOGR, 35 * * DMCZ,DMCZT,DMDT,DMSUN,DMSUN0,DMWMAX,DMWSUN,DTCZ,DTENV, 36 * * DTMAX,DTTOT,EACC,ETHAV,ETHNEW,F1,F2,FAC,FACR,FX,HCOMP,HD, 37 * * HJM,HP,HR,HS,HS1,HS2,HSTOT,HSTOT1,HSTOT2,HT,OMEGA,P,RATE, 38 * * RATIO,RB,RINIT,RMU,RNEW,RTOT,SENVOLD,SMASS,SUMDM,SURFI,T, 39 * * CJSUM,DMJ,SUMM,SUMX, 40 * * TEST,WSURF 41 INTEGER I,J,JENV,M 42 LOGICAL*4 LDONE,LNEW 43C 9/25/26 MHP via Claude 44C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/, 45C and COMMON/CONST3/ into a single COMMON/CONSTANTS/ 46C 9/24/26 MHP via Claude 47C code changed to remove dead CMIXL from COMMON/CONST3/; it is now 48C passed as an explicit argument where actually used 49C 10/12/26 MHP via Claude 50C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 51C namelist variables that the physics routines read from 27 COMMON 52C blocks which also held derived or working members. Those blocks keep 53C only their non-namelist members; five of them -- NEWMX, VMULT, 54C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 55C Derivation-only namelist variables from the same blocks went to 56C COMMON/NMLBYPASS/ instead, not here. 57C 58C Declared identically in all 121 files that hold any of these members. 59C Generated from one member list, so every copy is byte-identical; the 60C linker's 'size of symbol changed' warning is the check that they are. 61C mhp 8/10 added option to scale mass loss rate by rotation; requires rotating model 62C KC 2025-05-30 reordered common block elements 63C COMMON/MASSCHG3/LSOLWIND,DMSUN,DMWSUN,DMWMAX 64C 10/12/26 MHP via Claude 65C code changed to drop COMMON/MASSCHG3/ entirely. mdot.f was its only 66C declaring file once LSOLWIND moved to the consolidated block, so 67C DMSUN, DMWSUN and DMWMAX are now locals. They are given an explicit 68C zero below rather than relying on static storage happening to be 69C zeroed, which makes the guarded branch provably inert until these 70C three are given a real source. 71 DATA DMSUN,DMWSUN,DMWMAX/3*0.0D0/ 72 DIMENSION HCOMP(15,JSON),HS(JSON),HS1(JSON),HS2(JSON), 73 * * HT(JSON),HJM(JSON),OMEGA(JSON),HR(JSON), 74 * * HP(JSON),HD(JSON) 75 SAVE 76 SENVOLD = HS(M) - HSTOT 77C MHP 8/10- CHECK FOR SCALED SOLAR WIND MASS LOSS 78C 10/12/26 MHP via Claude 79C BUG FIX: added the DMWSUN guard below. DMSUN, DMWSUN and DMWMAX are 80C set NOWHERE in the code -- no assignment, no DATA statement, and no 81C namelist entry in any of the 315 input files -- and COMMON/MASSCHG3/ is 82C declared only here, so all three are permanently zero. The branch 83C therefore evaluated (OMEGA(M)/0)**2, and DMDT = 0*Infinity = NaN. 84C 85C It has never fired, for two separate reasons: parmin.f's 86C COMMON/MASSCHG3/ declaration is commented out (two variants, both 87C inactive), so the PHYSICS namelist's LSOLWIND is a local there and its 88C value never reaches this routine, where LSOLWIND is itself never 89C written and so always .FALSE.; and every one of the 315 namelists sets 90C LMDOT = .FALSE., which makes MASSLOSS return before calling MDOT. 91C 92C That matters because the namelist consolidation will put LSOLWIND in 93C the shared block and so repair the propagation. Nine namelists set 94C LSOLWIND = .TRUE. and all nine also set LROT = .TRUE., so without this 95C guard they would start entering the branch and taking a NaN mass-loss 96C rate. The option stays inert until DMSUN/DMWSUN/DMWMAX are given a 97C source, which is a separate decision for the experimental mass-loss 98C code. 99 IF(LSOLWIND .AND. LROT .AND. DMWSUN.GT.0.0D0) THEN 100 F1 = (OMEGA(M)/DMWSUN)**2 101 F2 = (DMWMAX/DMWSUN)**2 102 DMDT = DMSUN*MIN(F1,F2) 103 WRITE(*,*)OMEGA(M),DMDT 104 ENDIF 105C LMDOT = .TRUE. 106C DMDT = 4.0D-10 107C MAXIMUM FRACTION OF CZ MASS REMOVED PER TIMESTEP 108C FRAC = 0.1D0 109C IF(.NOT.LMDOT)GOTO 9999 110C CONVERT FROM SOLAR MASSES/YEAR TO GM/SEC 111 RATE = ABS(DMDT)*CMSUN/CSECYR 112C THE SUM OF THE MASSES OF ALL SHELLS (E.G. TO JENV - 1) 113C SHOULD BE USED RATHER THAN THE SUM OF THE MASSES DOWN 114C TO THE MIDPOINT OF THE BOTTOM SHELL 115C DMCZT = 10.0D0**HSTOT - HS1(JENV) 116 DMCZT = SUMDM 117C COMPUTE MASS OF SURFACE CONVECTION ZONE BELOW FITTING POINT 118 IF(JENV.EQ.M) STOP911 119 IF(JENV.GT.1) THEN 120 DMCZ = HS1(M) - HS1(JENV) 121 ELSE 122 DMCZ = HS1(M) 123 ENDIF 124C COMPUTE MAXIMUM TIMESTEP BASED ON NOT REMOVING TOO MUCH MASS FROM THE 125C SURFACE CZ (FCZDMDT) OR AS A FRACTION OF THE TOTAL MASS (FTOTDMDT) 126C 10/29/26 MHP via Claude: no mass-change limit for a zero rate (these 127C divided by zero). 128 IF(RATE.GT.0.0D0)THEN 129 DTTOT = FTOTDMDT*SMASS*CMSUN/RATE 130 DTCZ = FCZDMDT*DMCZ/RATE 131 DTMAX = MIN(DTTOT,DTCZ) 132 ELSE 133 DTMAX = DELTS 134 ENDIF 135C RESTRICT TIMESTEP TO ADD NO MORE THAN 1/2 OF THE CURRENT MASS 136C BEYOND THE FITTING POINT TO THE STAR. 137 IF(RATE.GT.0.0D0)THEN 138 DTENV = 0.5D0*(10.0D0**HSTOT-10.0D0**HS(M))/RATE 139 DTMAX = MIN(DTMAX,DTENV) 140 ENDIF 141 IF(DTMAX.LT.DELTS)THEN 142 WRITE(*,10)DELTS,DTMAX 143 10 10 FORMAT(' TIMESTEP REDUCED FOR MASS LOSS - OLD,NEW',1P2E12.3) 144 DELTS = DTMAX 145 ENDIF 146C mhp 8/10 turn mass loss off when disk exhausted only when dm /dt > 0, e.g. accretion 147 IF(LDISK .AND. DMDT.GT.0.0D0)THEN 148 TEST = SAGE + 1.0d-9*DELTS/CSECYR 149 IF(TEST.GT.TDISK)THEN 150 DELTS = (TDISK-SAGE)*1.0D9*CSECYR 151 LDONE = .TRUE. 152 ELSE 153 LDONE = .FALSE. 154 ENDIF 155 ELSE 156 LDONE = .FALSE. 157 ENDIF 158C FINAL AMOUNT OF MASS LOSS INFERRED IN CGS UNITS. 159 DELM = DMDT*CMSUN*DELTS/CSECYR 160C COMPUTE THE MEAN THERMAL ENERGY 161C CONTENT OF THE CONVECTION ZONE FOR MODELS WITH ACCRETION. 162 IF(DELM.GT.0.0D0)THEN 163 ETHNEW = (EACC*DELM + ETHAV*SUMDM)/(DELM+SUMDM) 164 FAC = LOG10(ETHNEW/ETHAV) 165C WRITE(*,912) JENV,EACC,ETHAV,ETHNEW,FAC,SUMDM,DELM 166C 912 FORMAT(I5,1P6E12.3) 167C OVERALL SCALE FACTOR IN R 168 T = 10.0D0**HT(JENV) 169 P = 10.0D0**HP(JENV) 170 D = 10.0D0**HD(JENV) 171 BETA = 1.0D0- (CA3*T**4/P) 172 RMU = P*BETA/(D*T) 173C DLOGEN = (DELM/DMCZ)*(SACC-SCEN)/RMU 174 DLOGEN = (DELM/DMCZ)*SCEN/RMU 175C DLNM = LOG(DMCZ+DELM)-LOG(DMCZ) 176 DLNM = 0.0D0 177 DLOGR = (CC23*DLOGEN - CC13*DLNM)/CLN 178 DLOGT = DLNM/CLN - DLOGR 179 DLOGP = 2.0D0*DLNM/CLN -4.0D0*DLOGR 180C WRITE(*,*)DLOGR,DLOGP,DLOGT 181 IF(JENV.EQ.1)THEN 182 DO J = JENV,M 183 HR(J) = HR(J)+DLOGR 184 HP(J) = HP(J)+DLOGP 185 HT(J) = HT(J)+DLOGT 186 END DO 187 ELSE 188 HP(JENV) = HP(JENV)+DLOGP 189 HT(JENV) = HT(JENV)+DLOGT 190 RB = 10.0D0**HR(JENV) 191 FACR = 10.0D0**DLOGR 192 DO J = JENV+1,M 193 RINIT = 10.0D0**HR(J) 194 RNEW = RB+FACR*(RINIT-RB) 195 HR(J) = LOG10(RNEW) 196 END DO 197 ENDIF 198C DO J = JENV,M 199C HT(J) = HT(J) + FAC 200C END DO 201C WRITE(*,*)HT(15),HT(16),HT(17) 202 ENDIF 203C FOR MODELS WITH ROTATION AND MASS LOSS REMOVE ANGULAR MOMENTUM. 204C IN THIS CASE WE ASSUME A THERMAL WIND WHERE THE SURFACE MATTER 205C CARRIES AWAY ONLY ITS LOCAL ANGULAR MOMENTUM PER UNIT MASS. 206 IF(LROT .AND. DELM.LT.0.0D0)THEN 207C MOMENT OF INERTIA PER UNIT MASS AT THE SURFACE. 208 SURFI = 2.0D0/3.0D0 * RTOT**2 209 IF(WALPCZ.GE.0.0D0)THEN 210C SOLID BODY CZ ROTATION 211 DELJ = OMEGA(M)*SURFI*DELM 212 ELSE IF(WALPCZ.LE.-2.0D0)THEN 213C CONSTANT J/M 214 DELJ = HJM(M)*DELM 215 ELSE 216C 10/29/26 MHP via Claude 217C Extrapolate with the law WCZ imposes, OMEGA = C*R**WALPCZ, so 218C WSURF = OMEGA(M)*(RTOT/R(M))**WALPCZ. The old form had the ratio 219C inverted, so for WALPCZ = -1 the surface spun faster than R(M). 220 WSURF = OMEGA(M)*(RTOT/10.0D0**HR(M))**WALPCZ 221 DELJ = WSURF*SURFI*DELM 222 ENDIF 223C FIND TOTAL STARTING ANGULAR MOMENTUM OF C.Z. 224 CJTOT = HJM(JENV)*HS2(JENV) 225 DO J = JENV+1,M 226 CJTOT = CJTOT + HJM(J)*HS2(J) 227 END DO 228 CJNEW = (CJTOT+DELJ) 229 WRITE(*,22)DELM,DELJ,OMEGA(M),CJNEW,CJTOT,DMCZT 230 22 22 FORMAT(1P6E13.4) 231C 10/29/26 MHP via Claude 232C J/M used to be scaled here by CJNEW/CJTOT, before the shell masses 233C shrank below. Shrinking them at that J/M removed a further 234C <J/M>*|DELM|, so the zone lost DELJ plus that. J/M is now scaled after 235C the new shell masses are known, so the zone holds exactly CJNEW. 236 ENDIF 237C REMOVE OR ADD MASS IN THE CONVECTION ZONE. 238C HS1 IS THE LOCATION OF THE ZONE CENTERS IN GM; 239C HS IS THE BASE-10 LOG OF HS1 240 IF(DELM.LT.0.0D0)THEN 241 FX = (DMCZ + DELM)/DMCZ 242C 10/29/26 MHP via Claude 243C For a fully convective star DMCZ is HS1(M), the mass from the centre, 244C so compress about the centre; compressing about HS1(1) moved HS1(M) by 245C less than DELM. 246C WRITE(*,*)JENV,M,DMCZ,DELM,FX,HS1(M),HS(M) 247 IF(JENV.EQ.1)THEN 248 DO I = 1,M 249 HS1(I) = FX*HS1(I) 250 HS(I) = LOG10(HS1(I)) 251 END DO 252 ELSE 253 DO I = JENV+1,M 254 HS1(I) = HS1(JENV)+FX*(HS1(I)-HS1(JENV)) 255 HS(I) = LOG10(HS1(I)) 256 END DO 257 ENDIF 258 ENDIF 259C WRITE(*,*)HS1(M),HS(M) 260C CORRECT TOTAL MASS IN SOLAR UNITS (SMASS) AND 261C LOG OF TOTAL MASS IN GRAMS (HSTOT,STOTAL) 262 SMASS = SMASS + DELM/CMSUN 263 SMASS0 = SMASS 264 DMSUN0 = DELM/CMSUN 265 WRITE(*,20)SMASS,DMSUN0 266 20 20 FORMAT('MASS LOSS APPLIED - NEW M,DEL M',1P2E19.10) 267 HSTOT1 = 10.0D0**HSTOT 268 HSTOT2 = HSTOT1 + DELM 269 HSTOT = LOG10(HSTOT2) 270 STOTAL = HSTOT 271C CORRECT MASS CONTENTS OF INDIVIDUAL SHELLS (HS2, IN GM) 272C 10/29/26 MHP via Claude 273C With JENV = 1 this read HS1(0). The innermost shell runs from the 274C centre to the midpoint of HS1(1) and HS1(2) (the HS1(0) = -HS1(1) 275C convention HPOINT uses), so HS2(1) = (HS1(1)+HS1(2))/2. 276 IF(JENV.EQ.1)HS2(1) = 0.5D0*(HS1(1)+HS1(2)) 277 DO I = MAX(JENV,2),M-1 278 HS2(I) = 0.5D0*(HS1(I+1)-HS1(I-1)) 279 END DO 280 HS2(M) = HSTOT2 - 0.5D0*(HS1(M)+HS1(M-1)) 281C 10/29/26 MHP via Claude: scale J/M on the new shell masses (see above). 282 IF(LROT .AND. DELM.LT.0.0D0)THEN 283 CJSUM = 0.0D0 284 DO J = JENV,M 285 CJSUM = CJSUM + HJM(J)*HS2(J) 286 END DO 287 RATIO = CJNEW/CJSUM 288 DO J = JENV,M 289 HJM(J) = HJM(J)*RATIO 290 END DO 291 ENDIF 292C 07/02RESET SENV 293 SENV = HS(M) - HSTOT 294C RECOMPUTE SURFACE BOUNDRY CONDITION 295 LNEW = .TRUE. 296C REMIX THE SURFACE CONVECTION ZONE IF MDOT IS POSITIVE. 297C WRITE(*,911)DELM,DMCZT,COMPACC(12),HCOMP(12,M), 298C * COMPACC(4),HCOMP(4,M) 299C 911 FORMAT(1P6E12.3) 300 IF(DELM.GT.0.0D0)THEN 301C 10/29/26 MHP via Claude 302C Mix by the zone's own mass-weighted abundances, sum(X dm), over the 303C shells' masses before the accreted DELM was added to HS2(M). The old 304C form used the base shell's X for the whole zone with the mass DMCZT, 305C which conserved species only if the zone was already uniform. 306 DO I = 1,11 307 SUMX = 0.0D0 308 SUMM = 0.0D0 309 DO J = JENV,M 310 DMJ = HS2(J) 311 IF(J.EQ.M)DMJ = DMJ - DELM 312 SUMX = SUMX + HCOMP(I,J)*DMJ 313 SUMM = SUMM + DMJ 314 END DO 315 ABUND = (SUMX+COMPACC(I)*DELM)/(SUMM+DELM) 316 DO J = JENV,M 317 HCOMP(I,J) = ABUND 318 END DO 319 END DO 320 FMASSACC = DELM 321 ENDIF 322C 9999 CONTINUE 323 IF(LDONE)LMDOT = .FALSE. 324C WRITE(*,*)M,HS(M),HSTOT 325 RETURN 326 END 327

1"""Python version of mdot.f (YREC 7be4133): mass loss or accretion over a 2step.""" 3 4import math 5import numpy as np 6 7from pyrec import comp_mod 8from pyrec import constants_mod 9from pyrec import deuter_mod 10from pyrec import disk_mod 11from pyrec import masschg2_mod 12from pyrec import nmlpar_mod 13 14# DATA DMSUN,DMWSUN,DMWMAX/3*0.0D0/: the scaled-solar-wind parameters, set 15# nowhere else (formerly COMMON/MASSCHG3/), so the branch below is inert. 16DMSUN = 0.0 17DMWSUN = 0.0 18DMWMAX = 0.0 19 20 21def mdot(delts: float, hcomp: np.ndarray, hd: np.ndarray, hjm: np.ndarray, hp: np.ndarray, hr: np.ndarray, 22 hs: np.ndarray, hs1: np.ndarray, hs2: np.ndarray, hstot: float, ht: np.ndarray, jenv: int, lnew: bool, m: int, 23 omega: np.ndarray, rmu: float, rtot: float, smass: float, dmdt: float, eacc: float, ethav: float, sumdm: float, 24 senvold: float) -> tuple[float, float, bool, float, float, float, float]: 25 """Apply a mass change (loss, or accretion if dmdt > 0) to the model 26 over a step: limit the step, rescale the mass coordinates, total mass 27 and shell masses, remove angular momentum with lost mass, and mix 28 accreted matter into the surface convection zone. 29 30 Translates: mdot.f#MDOT (7be4133) 31 32 Arguments: the step delts (s); the model hcomp[species, shell], hd, hjm, 33 hp, hr, hs, hs1, hs2, ht, omega (arrays, index 0 unused, changed in 34 place); log total mass hstot; jenv, the base of the surface convection 35 zone; m, the number of shells; rmu; rtot, the radius; smass (Msun); 36 dmdt (Msun/yr); eacc, ethav (no effect); sumdm, the zone's mass. 37 38 Returns (delts, hstot, lnew, rmu, smass, dmdt, senvold): the possibly 39 shortened step, the new log total mass, lnew = True (recompute the 40 surface boundary condition), rmu (recomputed for accretion), the new 41 mass, dmdt (changed only by the inert solar-wind branch) and senvold, 42 the old hs[m] - hstot. 43 44 """ 45 46 c = constants_mod 47 senvold = hs[m] - hstot # Fortran line 76. 48 49 # Scaled solar wind mass loss (inert: DMWSUN = 0) (MHP 8/10) (Fortran 50 # lines 99–104). 51 if nmlpar_mod.lsolwind and nmlpar_mod.lrot and DMWSUN > 0.0: 52 f1 = (omega[m] / DMWSUN)**2 53 f2 = (DMWMAX / DMWSUN)**2 54 dmdt = DMSUN * min(f1, f2) 55 print(omega[m], dmdt) 56 57 # Convert from solar masses per year to g/s (Fortran lines 111–116). 58 rate = abs(dmdt) * c.cmsun / c.csecyr 59 dmczt = sumdm 60 61 # Mass of the surface convection zone below the fitting point (Fortran 62 # lines 118–123). 63 # Suspected bug (see the page): a hard stop with code 911 only. 64 if jenv == m: 65 raise SystemExit(911) 66 if jenv > 1: 67 dmcz = hs1[m] - hs1[jenv] 68 else: 69 dmcz = hs1[m] 70 71 # Maximum step: not too much mass from the zone (FCZDMDT) or from the 72 # star (FTOTDMDT) (Fortran lines 128–145). 73 # A zero rate sets no limit. 74 if rate > 0.0: 75 dttot = nmlpar_mod.ftotdmdt * smass * c.cmsun / rate 76 dtcz = nmlpar_mod.fczdmdt * dmcz / rate 77 dtmax = min(dttot, dtcz) 78 else: 79 dtmax = delts 80 # At most half the mass beyond the fitting point. 81 if rate > 0.0: 82 dtenv = 0.5 * (10.0**hstot - 10.0**hs[m]) / rate 83 dtmax = min(dtmax, dtenv) 84 if dtmax < delts: 85 print(f" TIMESTEP REDUCED FOR MASS LOSS - OLD,NEW" 86 f"{delts:12.3E}{dtmax:12.3E}") # format 10 87 delts = dtmax 88 89 # Turn accretion off once the disk is exhausted (MHP 8/10) (Fortran 90 # lines 147–157). 91 # Suspected bug (see the page): a negative step if sage > tdisk. 92 if nmlpar_mod.ldisk and dmdt > 0.0: 93 test = disk_mod.sage + 1.0e-9 * delts / c.csecyr 94 if test > nmlpar_mod.tdisk: 95 delts = (nmlpar_mod.tdisk - disk_mod.sage) * 1.0e9 * c.csecyr 96 ldone = True 97 else: 98 ldone = False 99 else: 100 ldone = False 101 102 # Mass change in g (Fortran line 159). 103 delm = dmdt * c.cmsun * delts / c.csecyr 104 105 # Accretion: mean thermal energy and structure changes of the zone 106 # (Fortran lines 162–202). 107 if delm > 0.0: 108 # ethnew and fac are not used (see the page). 109 ethnew = (eacc * delm + ethav * sumdm) / (delm + sumdm) 110 fac = math.log10(ethnew / ethav) # noqa: F841 (unused in the Fortran too) 111 t = 10.0**ht[jenv] 112 p = 10.0**hp[jenv] 113 d = 10.0**hd[jenv] 114 beta = 1.0 - (c.ca3 * t**4 / p) 115 rmu = p * beta / (d * t) 116 # SCEN = 0, so these changes vanish (unfinished code, see the page). 117 # Question (see the page): SCEN alone, or SACC - SCEN? 118 dlogen = (delm / dmcz) * masschg2_mod.scen / rmu 119 dlnm = 0.0 120 dlogr = (c.cc23 * dlogen - c.cc13 * dlnm) / c.cln 121 masschg2_mod.dlogt = dlnm / c.cln - dlogr 122 masschg2_mod.dlogp = 2.0 * dlnm / c.cln - 4.0 * dlogr 123 if jenv == 1: 124 for j in range(jenv, m + 1): 125 hr[j] = hr[j] + dlogr 126 hp[j] = hp[j] + masschg2_mod.dlogp 127 ht[j] = ht[j] + masschg2_mod.dlogt 128 else: 129 hp[jenv] = hp[jenv] + masschg2_mod.dlogp 130 ht[jenv] = ht[jenv] + masschg2_mod.dlogt 131 rb = 10.0**hr[jenv] 132 facr = 10.0**dlogr 133 for j in range(jenv + 1, m + 1): 134 rinit = 10.0**hr[j] 135 rnew = rb + facr * (rinit - rb) 136 hr[j] = math.log10(rnew) 137 138 # Rotation and mass loss: the lost matter carries its local specific 139 # angular momentum (a thermal wind) (Fortran lines 206–236). 140 if nmlpar_mod.lrot and delm < 0.0: 141 surfi = 2.0 / 3.0 * rtot**2 142 walpcz = nmlpar_mod.walpcz 143 if walpcz >= 0.0: 144 # Solid-body zone. 145 delj = omega[m] * surfi * delm 146 elif walpcz <= -2.0: 147 # Constant J/M. 148 delj = hjm[m] * delm 149 else: 150 # Extrapolate with the law WCZ imposes, Omega = C*r**WALPCZ. 151 wsurf = omega[m] * (rtot / 10.0**hr[m])**walpcz 152 delj = wsurf * surfi * delm 153 # Total starting angular momentum of the zone. 154 cjtot = hjm[jenv] * hs2[jenv] 155 for j in range(jenv + 1, m + 1): 156 cjtot = cjtot + hjm[j] * hs2[j] 157 cjnew = cjtot + delj 158 print(f"{delm:13.4E}{delj:13.4E}{omega[m]:13.4E}{cjnew:13.4E}" 159 f"{cjtot:13.4E}{dmczt:13.4E}") # format 22 160 # J/M is scaled below, once the new shell masses are known. 161 162 # Remove mass in the convection zone: compress toward its base 163 # (hs1: shell centres in g; hs: log10 of hs1) (Fortran lines 240–258). 164 # A fully convective star (dmcz = hs1[m]) is compressed about the center. 165 if delm < 0.0: 166 fx = (dmcz + delm) / dmcz 167 if jenv == 1: 168 for i in range(1, m + 1): 169 hs1[i] = fx * hs1[i] 170 hs[i] = math.log10(hs1[i]) 171 else: 172 for i in range(jenv + 1, m + 1): 173 hs1[i] = hs1[jenv] + fx * (hs1[i] - hs1[jenv]) 174 hs[i] = math.log10(hs1[i]) 175 176 # New total mass in solar units (smass) and log g (hstot, stotal) 177 # (Fortran lines 262–270). 178 smass = smass + delm / c.cmsun 179 masschg2_mod.smass0 = smass 180 dmsun0 = delm / c.cmsun 181 print(f"MASS LOSS APPLIED - NEW M,DEL M{smass:19.10E}{dmsun0:19.10E}") 182 hstot1 = 10.0**hstot 183 hstot2 = hstot1 + delm 184 hstot = math.log10(hstot2) 185 comp_mod.stotal = hstot 186 187 # Shell masses hs2 (g) (Fortran lines 276–280). 188 # With jenv = 1 the innermost shell runs from the center to the 189 # midpoint of hs1[1] and hs1[2] (the hs1[0] = -hs1[1] convention). 190 if jenv == 1: 191 hs2[1] = 0.5 * (hs1[1] + hs1[2]) 192 for i in range(max(jenv, 2), m): 193 hs2[i] = 0.5 * (hs1[i + 1] - hs1[i - 1]) 194 hs2[m] = hstot2 - 0.5 * (hs1[m] + hs1[m - 1]) 195 196 # Scale J/M on the new shell masses, so the zone holds exactly cjnew 197 # (Fortran lines 282–291). 198 if nmlpar_mod.lrot and delm < 0.0: 199 cjsum = 0.0 200 for j in range(jenv, m + 1): 201 cjsum = cjsum + hjm[j] * hs2[j] 202 ratio = cjnew / cjsum 203 for j in range(jenv, m + 1): 204 hjm[j] = hjm[j] * ratio 205 206 # Fortran lines 293–295. 207 comp_mod.senv = hs[m] - hstot 208 lnew = True 209 210 # Accretion: remix the zone with the accreted composition (Fortran 211 # lines 300–321). 212 if delm > 0.0: 213 for i in range(1, 12): 214 # The zone's own sum(X dm), over the shell masses before the 215 # accreted delm was added to hs2[m]. 216 sumx = 0.0 217 summ = 0.0 218 for j in range(jenv, m + 1): 219 dmj = hs2[j] 220 if j == m: 221 dmj = dmj - delm 222 sumx = sumx + hcomp[i, j] * dmj 223 summ = summ + dmj 224 abund = (sumx + nmlpar_mod.compacc[i] * delm) / (summ + delm) 225 for j in range(jenv, m + 1): 226 hcomp[i, j] = abund 227 deuter_mod.fmassacc = delm 228 229 # Fortran line 323. 230 if ldone: 231 nmlpar_mod.lmdot = False 232 return delts, hstot, lnew, rmu, smass, dmdt, senvold