liburn.f

561 lines · 2 definitions · on GitHub (7be4133)

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

Variables on this page (28)

bbe9old, bli6old, bli7old, cln, compacc, delam, delrm, dmdt0, fmassacc, hro, ishort, jenvos, json, lbliold, linstb, lmdot, lovste, lrot, pscal0, pscal1, rbase0, rbe9, rbe90, rli6, rli60, rli7, rli70, sdel

Purpose

LIBURN burns the light elements ⁶Li, ⁷Li and ⁹Be (HCOMP(13), HCOMP(14), HCOMP(15)) over one time step: in the radiative interior shell by shell, and in the surface convection zone as one mixed reservoir, with the rates changing during the step. Its accuracy is controlled by Richardson-type extrapolation over ever finer substeps. Its comment: "lithium and beryllium burning in a non-rotating model."

How it fits

Who calls it: the main program, main.f, at the end of each step when light elements are followed (LEXCOM) and the model is not rotating with rotational mixing; and GETW for rotating models. The cheaper variant LIBURN2 (one step, no extrapolation) is used by GETW and BURSMIX.
Inputs: DELTS (the step; main.f passes it in seconds), HCOMP, HR, HS1, HS2, HT, M; JNENV0 and JNENV, the index of the base of the surface convection zone at the start and at the end of the step. The rates at the start (OLDRAT_MOD) and end (NEWRAT_MOD) of the step, from LIRATE88.

The step is in seconds. The comments call DELT "the fractional timestep in years" (lines 204–208 ↗; lines 382–385 ↗), but DT = FRAC*DELTS, main.f passes DELTS in seconds (main.f, lines 1068–1069 ↗), and the rates from LIRATE88 are per second (RLI6 = FLI6*RHOX, $\rho X N_A\langle\sigma v\rangle$: lirate88.f, lines 131–133 ↗), so the comment's "in years" is stale.
Output: the new HCOMP(13..15).
Calls: RATEXT (rational extrapolation) and SAFEDIVEXP (below).
Modules: NEWRAT_MOD, OLDRAT_MOD (the rates; with LIRATE88, LIBURN2, GETW, BURSMIX); SCRTCH_MOD (SDEL, the gradients $\nabla_\mathrm{rad}$, $\nabla$, $\nabla_\mathrm{ad}$ from PHYSIC); MDPHY_MOD (the same for rotating models); OLDMOD_MOD (HRO, the old radii); LIOV_MOD (PSCAL0, PSCAL1, the overshoot distances at the start and end of the step, set in main.f); PREVCZ_MOD (RBASE0, the end-of-step rates BLI6OLD… with the flag LBLIOLD that says they are valid, and JENVOS, kept for the next step); DEUTER_MOD (FMASSACC, for mass accretion); NMLPAR_MOD (the namelist switches and values: LOVSTE, envelope overshoot; LROT, LINSTB; for accretion LMDOT, DMDT0, COMPACC); CONSTANTS_MOD, IOUNITS_MOD, DIMENSIONS_MOD.

The method

Base of the convection zone. Locate the base between the last convective point and the one below by linear interpolation of $\nabla_\mathrm{ad} - \nabla_\mathrm{rad}$: an offset $f_x \in [-\tfrac12, \tfrac12]$ of a zone. With envelope overshoot, move the base down by the overshoot distance (PSCAL0, PSCAL1) at the start and end of the step. At the start, the base's radius is the midpoint between HRO(JNENV0) and the point below; if the start-of-step model is fully convective (JNENV0 = 1), the radius below shell 1 is taken as zero. At the end, the offset of the new base is measured in the zone width of the end-of-step model (HR).
Radiative interior (points below both bases, from the centre out, stopping at the first point where the ⁹Be rate, at the start or at the end, is at most $10^{-32}$; points where all three abundances are below $10^{-24}$ are skipped). Above $\log T = 7$ the elements are destroyed outright. Elsewhere, for $n$ = 2, 4, 6, 8, 12, … , 96 substeps of $h = \Delta t/n$: the rate at substep $k$ is interpolated linearly in its log between start and end, $$\ln\lambda_k = \ln\lambda_0 + \tfrac{k - 1/2}{n}\,(\ln\lambda_1 - \ln\lambda_0),$$ and the abundance decays as $x \leftarrow x\,e^{-\lambda_k h}$ (set to 0 when $\lambda_k h \ge 30$). The results for successive $n$ are extrapolated to $h \to 0$ in $h^2$ by RATEXT; it stops when the relative error estimates are below $10^{-6}$. Values below $10^{-24}$ are set to 0.

Question 1 (line 194 ↗; liburn2.f, line 158 ↗). The radiative loop's exit. It leaves the loop at the first point (from the center) where the ⁹Be rate is at most $10^{-32}$, assuming the rates only fall outward; LIBURN2 skips such points instead. The last loop, over the region left by a retreating zone, leaves in the same way here (line 523 ↗) and in LIBURN2 (liburn2.f, line 334 ↗), so the choice reaches all three loops. Which is meant?

Question 2 (lines 93–94 ↗; line 197 ↗; lines 232–246 ↗; lines 280–284 ↗). Hard-coded values. The destruction temperature $\log T = 7$, the cut $\lambda h \ge 30$, the floor $10^{-24}$, the substep counts and the tolerances $10^{-6}$ are fixed in the code. Are they meant to stay fixed, or should any be a namelist setting (a run with tolerances ten times smaller would show whether the results depend on them)? 3. Convection zone. The zone's starting abundance is the mass-weighted mean over the zone at the start (or at the end, if the zone deepened). Its rate is the mass-weighted mean $\bar\lambda = \sum\lambda_i\Delta m_i/\sum\Delta m_i$ at the start (or the previous step's end value, BLI6OLD…, when the flag LBLIOLD is set) and at the end, the bottom point's mass and rate adjusted for the offset $f_x$. The same substepping and extrapolation as in step 2 follow. With accretion (LMDOT, DMDT0 > 0), accreted material of composition COMPACC is added at each substep, assumed to burn for half a substep on arrival, and mixed in by mass with the factor FMASSACC. All points from the shallower of the two bases (MAXJ) to the surface get the result, and the end rates are stored for the next step with LBLIOLD set. The flag is cleared by main.f at each kind card, by LIBURN when the whole zone is below the burning threshold (the zone is skipped), and by LIBURN2, which always computes its start rates. The adjusted rates of the bottom point stay in NEWRAT_MOD (only its mass is restored). 4. Region left by a retreating zone (points JENV0..JENV−1, between the old and the new base): each is burned once with a rate weighted geometrically between its own end-of-step rate and the zone's start rate, $$x \leftarrow x\,\exp\!\left(-\Delta t\,\lambda_\mathrm{rad}^{\,f}\,\bar\lambda_0^{\,1-f}\right),\qquad f = \frac{m - m_\mathrm{beg}}{m_\mathrm{beg} - m_\mathrm{end}},$$ with $m$ = HS1 (the mass inside the point, in g), and $m_\mathrm{beg}$, $m_\mathrm{end}$ the masses at the old and the new base (HSBEG, HSEND; midpoints between the base point and the one below). The comment says that a point a fraction $f$ of the way (in mass) from the old base to the new one spent the fraction $f$ of the step in the zone, so $f$ should lie between 0 and 1; as written, $f$ is negative (bug 1 below).

Design

With accretion it prints a line to the terminal (WRITE(*,913)) at every level of the extrapolation.

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 sign of `FRAD`.

(lines 509–528 ↗) For the points left by a retreating zone (JENV > JENV0), DMASS = HSBEG - HSEND is negative, since HS1 grows outward, while HS1(I) - HSBEG is positive, so FRAD lies between −1 and 0. The rate $\lambda_\mathrm{rad}^{f}\bar\lambda_0^{1-f}$ then extrapolates beyond the zone's rate instead of interpolating between the two rates. Even with the sign corrected, the code gives the radiative rate the weight $f$, while the comment (lines 509–512) says that a point the fraction $f$ of the way from the old base to the new one spent the fraction $f$ of the step in the zone, which calls for the weight $f$ on the zone's rate: A point at the old base ($f$ = 0) left the zone at once. What follows: The light elements in the region a retreating zone leaves are burned at a rate up to $\bar\lambda_0^2/\lambda_\mathrm{rad}$, outside the range of the two rates. Only with LEXCOM (off by default), on steps where the zone retreats.

A possible fix:

liburn.f

@@ lines 516–531 @@
          HSBEG = 0.0D0
       ENDIF
       HSEND = 0.5D0*(HS1(JENV)+HS1(JENV-1))
-      DMASS = HSBEG - HSEND
+C 11/09/26 KC via Claude
+C code changed to make FRAD the fraction of the step spent in the CZ,
+C between 0 at the old base and 1 at the new one, as the comment above
+C says, and to give it the weight of the convective rate.
+      DMASS = HSEND - HSBEG
       DO 140 I = JENV0,JENV-1
 C MHP 9/91 CHANGE TO AVOID DIVISION BY ZERO.
 C SKIP IF SHELL TEMPERATURE DROPS BELOW BURNING THRESHOLD.
          IF(RBE9(I).LE.1.0D-32)GOTO 145
          FRAD = (HS1(I)-HSBEG)/DMASS
-C USE FRAD*RADIATIVE RATE AND (1-FRAD)*CONVECTIVE RATE.
-         DLI6 = DELTS*EXP(FRAD*LOG(RLI6(I))+(1.0D0-FRAD)*BLI60)
-         DLI7 = DELTS*EXP(FRAD*LOG(RLI7(I))+(1.0D0-FRAD)*BLI70)
-         DBE9 = DELTS*EXP(FRAD*LOG(RBE9(I))+(1.0D0-FRAD)*BBE90)
+C USE FRAD*CONVECTIVE RATE AND (1-FRAD)*RADIATIVE RATE.
+         DLI6 = DELTS*EXP((1.0D0-FRAD)*LOG(RLI6(I))+FRAD*BLI60)
+         DLI7 = DELTS*EXP((1.0D0-FRAD)*LOG(RLI7(I))+FRAD*BLI70)
+         DBE9 = DELTS*EXP((1.0D0-FRAD)*LOG(RBE9(I))+FRAD*BBE90)
 C***REMEMBER TO ADD FAILSAFES FOR LARGE DEPLETION***
 C KC 2025-05-31 PREVENT FLOATING POINT EXCEPTION
 C          HCOMP(13,I) = HCOMP(13,I)/EXP(DLI6)

Checked by running YREC (Nov 9, 2026, the counts of testsuite/Test_m1p0feh+0p0GS98_base_HeIgnite 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): 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), 23 had steps where the zone retreated, and on every one FRAD lay between $-0.989$ and $-0.0007$, outside 0 to 1, as read: up to 11,399 points in testsuite/Test_m1p0feh+0p0GS98_base_HeIgnite, 1,973 in testsuite/Test_solarM22M_difkur, and fewer in the solar, low-mass, rotating, and helium-burning cases. With the fix, testsuite/Test_solarM22M_difkur keeps its 1,702 models and 4.568 Gyr, and its track's $\log L$, $\log T_\mathrm{eff}$, $\log T_c$, and central $X$ and $Y$ are unchanged, where 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}$) moves $\log L$ by up to 0.000167 dex; its .track file differs within 2%, and its .short and .store files beyond 2%; testsuite/Test_m0p3feh+0p0GS98_scv_ZAMS is identical in all three files. The track comparison covers only $\log L$, $\log T_\mathrm{eff}$, $\log T_c$, and central $X$ and $Y$, so the light elements, which this fix changes, were not compared beyond the 2% test, and the giant-branch case, with the most retreating steps (11,399), was not run with the fix. Intent: The comment above (lines 509–512) describes the opposite weighting. LIBURN2 has the same lines (liburn2.f, suspected bug 1).

2. (harmless in the runs) Points newly in a deepened zone are not burned.

(lines 191–193 ↗; lines 373–379 ↗; lines 496–500 ↗) When the zone deepens (JENV < JENV0), its starting abundance is the mean over the end-of-step zone, JENV to M, which includes the points JENV to JENV0−1; but the result is written only from MAXJ = JENV0 outward, and the radiative loop stops below JENV. Those points keep their abundances from before the step, so the zone is left with an unburned bottom and a burned top until it is next mixed (at the latest by MIXCZ, which main.f calls once the next step's structure has converged, main.f, line 977 ↗), and the burning of their mass over this step is lost. Only with LEXCOM (off by default), on steps where the zone deepens.

A possible fix:

liburn.f

@@ lines 493–499 @@
       IF(FLI7.LT.1.0D-24)FLI7=0.0D0
       FBE9 = YEXT(3)
       IF(FBE9.LT.1.0D-24)FBE9=0.0D0
-      DO 130 I = MAXJ,M
+C 11/09/26 KC via Claude
+C code changed to give the burned abundances to the whole zone at the
+C end of the step, including points that joined it as it deepened.
+      DO 130 I = JENV,M
          HCOMP(13,I) = FLI6
          HCOMP(14,I) = FLI7
          HCOMP(15,I) = FBE9

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 zone deepened in 23 of the 34 runs, up to 485 steps (testsuite/Test_m1p0feh+0p0GS98_base_HeIgnite), 359 in the 1 $M_\odot$ helium-burning case run with the AAG21 table, and 154 in testsuite/Test_3Msolar_GS_rot_gray_OPALSCV_SF2_GS98_OP_CF10_RGp2diffrotI3cz2, so the points that had just joined the zone were left unburned. With the fix, that rotating 3 $M_\odot$ case keeps its 1,853 models and 0.439709 Gyr, and its track's $\log L$, $\log T_\mathrm{eff}$, $\log T_c$, and central $X$ and $Y$ are unchanged, where 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}$) moves $\log L$ by up to 0.00758 dex; its .track, .short, and .store files change within 2%; testsuite/Test_m0p3feh+0p0GS98_scv_ZAMS is identical. The light elements were not compared beyond the 2% test, and the giant-branch case, where the zone deepens most, was not run with the fix. LIBURN2 has the same lines (liburn2.f, suspected bug 2).

3. (never occurred) ⁹Be is not guarded against overflow.

(lines 529–538 ↗) In the last loop, ⁶Li and ⁷Li go through SAFEDIVEXP, but ⁹Be is still divided by EXP(DBE9), and the comment "remember to add failsafes for large depletion" remains. For DBE9 above 709.78, EXP overflows: In the default build the quotient is then 0, the intended result; in the development build, with -ffpe-trap=overflow (GNU_DEV_FFLAGS, src/Makefile line 61), the run stops. Since ⁹Be burns at higher temperatures than ⁶Li and ⁷Li, its exponent is probably the smallest of the three, so this needs a point where the lithium exponents are larger still. Only with LEXCOM (off by default).

A possible fix:

liburn.f

@@ lines 534–540 @@
 C          HCOMP(14,I) = HCOMP(14,I)/EXP(DLI7)
          CALL SAFEDIVEXP(HCOMP(14,I),DLI7)
          IF(HCOMP(14,I).LT.1.0D-24)HCOMP(14,I)=0.0D0
-         HCOMP(15,I) = HCOMP(15,I)/EXP(DBE9)
+C 11/09/26 KC via Claude
+C code changed to guard Be9 as Li6 and Li7 are.
+         CALL SAFEDIVEXP(HCOMP(15,I),DBE9)
          IF(HCOMP(15,I).LT.1.0D-24)HCOMP(15,I)=0.0D0
   140 CONTINUE
   145 CONTINUE

Checked by running YREC (Nov 9, 2026, the counts of testsuite/Test_m1p0feh+0p0GS98_base_HeIgnite 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 division ran in 23 of the 34 runs, up to 11,399 times (testsuite/Test_m1p0feh+0p0GS98_base_HeIgnite), with DBE9 from about $2\times10^{-24}$ to 201, far below 709.78, so EXP never overflowed. It needs a retreating zone that leaves points hot enough to burn ⁹Be by more than $e^{709}$ in one step, which no case gives.

4. (harmless in the runs) The bottom point's rates are left adjusted.

(lines 334–343 ↗; lines 371–372 ↗) The comment says the rate and mass of the bottom point are adjusted for the offset $f_x$; afterwards only HS2(JENV) is restored, and RLI6, RLI7, RBE9 at JENV keep the adjusted values in NEWRAT_MOD. In main.f nothing reads them before LIRATE88 recomputes both sets at the next step; but GETW copies them into the start rates OLDRAT_MOD for the next substep (getw.f, lines 409–413 ↗), where a later mean start rate (loop 68, when LBLIOLD is false) weights the adjusted rate with the unadjusted mass. Correct for main.f, fragile for GETW.

A possible fix:

liburn.f

@@ lines 53–58 @@
      *       FM,FM0,FRAC,FRAD,FX,HCOMP,HR,HS1,HS2,HS2SAV,HSAV,HSBEG,
      *       HSEND,HT,R,RATBE9,RATLI6,RATLI7,RBASE,ROS,TBE9,TLI6,TLI7,
      *       XEST,YERR,YEST,YEXT,YTOL
+      REAL*8 RBE9SV,RLI6SV,RLI7SV
       INTEGER I,J,JENV,JENV0,JNENV,JNENV0,K,K2,M,MAXJ,MINJ,NSTEPE
       LOGICAL*4 LDOACC,LDONE

@@ lines 334–340 @@
 C USE THE LOCATION OF THE TRUE EDGE OF THE CONVECTION ZONE (FX, FOUND AT
 C BEGINNING OF SR) TO ADJUST THE BURNING RATE AND MASS OF THE BOTTOM POINT
 C SUCH THAT IT INCLUDES THE ENTIRE C.Z.
+C 11/09/26 KC via Claude
+C code changed to restore the bottom point's rates with its mass.
       HS2SAV = HS2(JENV)
+      RLI6SV = RLI6(JENV)
+      RLI7SV = RLI7(JENV)
+      RBE9SV = RBE9(JENV)
       IF(JENV.GT.1.AND.JENV.LT.M)THEN
          HS2(JENV) = HS2(JENV)+FX*(HS1(JENV)-HS1(JENV-1))
          RLI6(JENV) = RLI6(JENV)+0.5D0*FX*(RLI6(JENV-1)-RLI6(JENV))
@@ lines 368–375 @@
       BLI6 = LOG(BLI6/FM)
       BLI7 = LOG(BLI7/FM)
       BBE9 = LOG(BBE9/FM)
-C RESTORE ORIGINAL MASS TO BASE OF C.Z.
+C RESTORE ORIGINAL MASS AND RATES TO BASE OF C.Z.
       HS2(JENV) = HS2SAV
+      RLI6(JENV) = RLI6SV
+      RLI7(JENV) = RLI7SV
+      RBE9(JENV) = RBE9SV
 C FOR THE STARTING ABUNDANCE, USE THE MIXED ABUNDANCE FOR THE
 C DEEPER CZ.
       IF(JENV.LT.JENV0)THEN

with RLI6SV, RLI7SV, RBE9SV added to the REAL*8 declarations.

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 adjusted rates were left at the bottom point in 24 of the 34 runs (offsets FX from $-0.4997$ to 0.5), up to 7,145 times (testsuite/Test_m1p0feh+0p0GS98_base_HeIgnite), and 2,251 times in the rotating testsuite/Test_m0p8feh+0p0GS98_rotp8p0d0p1mag1e4std_TAMS, where GETW copies them into the start rates. With the fix, that rotating case keeps an identical track and .store file (763 models, 22.3501 Gyr), where 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}$) moves $\log L$ by up to 0.00014 dex; testsuite/Test_m0p3feh+0p0GS98_scv_ZAMS is identical too. Fragile code, with no effect today.

Definitions

subroutineLIBURN

lines 8–543 · GitHub

Light-element burning with substeps and extrapolation, radiative and convective parts (above).

Uses modules CONSTANTS_MOD, DEUTER_MOD, DIMENSIONS_MOD, IOUNITS_MOD, LIOV_MOD, MDPHY_MOD, NEWRAT_MOD, NMLPAR_MOD, OLDMOD_MOD, OLDRAT_MOD, PREVCZ_MOD, SCRTCH_MOD

Called by (file level), GETW

Flowchart of its calls, in order (4)

for I = 1,MINJ-1 (DO 50), lines 193–285

for J = 1,11 (DO 40), lines 203–273

RATEXT 255

for J = 1,11 (DO 120), lines 380–486

RATEXT 470

for I = JENV0,JENV-1 (DO 140), lines 520–539

SAFEDIVEXP 532→SAFEDIVEXP 535

Compiler warnings (1)

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

line 138: Equality comparison for REAL(8) [-Wcompare-reals]

subroutineSAFEDIVEXP

lines 548–561 · GitHub

KC's guard (2025-05-31): replaces $x$ by $x/e^{a}$ when $a < 709.78$ (the largest argument for which $e^a$ fits in double precision), and by 0 otherwise.

Called by LIBURN ×2

Change log

The changes since rp_1

Changed since rp_1. Declarations only: COMMON blocks to modules (the namelist members of /DPMIX/, /ROT/ and /MASSCHG/ to NMLPAR_MOD, /CONST1/ to CONSTANTS_MOD, /LUOUT/ to IOUNITS_MOD), JSON from DIMENSIONS_MOD, IMPLICIT NONE (also in SAFEDIVEXP).

Changed since a8a5950. The previous step's end-of-step zone rates are now reused, through the new flag LBLIOLD in PREVCZ_MOD; the old test BLI6OLD.LE.0 was always true for the log of a small rate, so the start rates were always recomputed. A fully convective start-of-step model no longer reads HRO(0). The overshoot zone width is taken from the end-of-step model alone (HR(J+1), not HRO(J+1)), which settles former question 2; the later questions are renumbered.

1C 2C 3C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 4C LIBURN 5C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 6C 11/91 HR ADDED TO CALL. 7C SUBROUTINE LIBURN(DELTS,HCOMP,HD,HR,HS1,HS2,HT,JNENV,JNENV0,M) ! KC 2025-05-31 8 SUBROUTINE LIBURN(DELTS,HCOMP,HR,HS1,HS2,HT,JNENV,JNENV0,M) 9C 10/15/26 MHP via Claude 10C code changed to take JSON from DIMENSIONS_MOD instead of declaring 11C PARAMETER(JSON=5000) here. 12 USE DIMENSIONS_MOD 13C 10/14/26 MHP via Claude 14C code changed to convert COMMON/SCRTCH/ into SCRTCH_MOD. 15 USE SCRTCH_MOD 16C 10/14/26 MHP via Claude 17C code changed to convert COMMON/PREVCZ/ into PREVCZ_MOD. 18 USE PREVCZ_MOD 19C 10/14/26 MHP via Claude 20C code changed to convert COMMON/OLDRAT/ into OLDRAT_MOD. 21 USE OLDRAT_MOD 22C 10/14/26 MHP via Claude 23C code changed to convert COMMON/OLDMOD/ into OLDMOD_MOD. 24 USE OLDMOD_MOD 25C 10/14/26 MHP via Claude 26C code changed to convert COMMON/NEWRAT/ into NEWRAT_MOD. 27 USE NEWRAT_MOD 28C 10/14/26 MHP via Claude 29C code changed to convert COMMON/MDPHY/ into MDPHY_MOD. 30 USE MDPHY_MOD 31C 10/14/26 MHP via Claude 32C code changed to convert COMMON/LIOV/ into LIOV_MOD. 33 USE LIOV_MOD 34C 10/14/26 MHP via Claude 35C code changed to convert COMMON/DEUTER/ into DEUTER_MOD. 36 USE DEUTER_MOD 37C 10/14/26 MHP via Claude 38C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD. 39 USE NMLPAR_MOD 40C 10/14/26 MHP via Claude 41C code changed to convert COMMON/CONSTANTS/ into CONSTANTS_MOD. 42 USE CONSTANTS_MOD 43C 10/14/26 MHP via Claude 44C code changed to convert COMMON/IOUNITS/ into IOUNITS_MOD. 45 USE IOUNITS_MOD 46C 10/20/26 MHP via Claude 47C code changed to replace implicit typing with IMPLICIT NONE and 48C explicit declarations. 49 IMPLICIT NONE 50 REAL*8 BBE9,BBE90,BLI6,BLI60,BLI7,BLI70,CZBE9,CZLI6,CZLI7,DBE9, 51 * * DDJ,DDJM1,DELTS,DLI6,DLI7,DMASS,DR,DT,FBE9,FBE90,FBE9A, 52 * * FBE9ADD,FLI6,FLI60,FLI6A,FLI6ADD,FLI7,FLI70,FLI7A,FLI7ADD, 53 * * FM,FM0,FRAC,FRAD,FX,HCOMP,HR,HS1,HS2,HS2SAV,HSAV,HSBEG, 54 * * HSEND,HT,R,RATBE9,RATLI6,RATLI7,RBASE,ROS,TBE9,TLI6,TLI7, 55 * * XEST,YERR,YEST,YEXT,YTOL 56 INTEGER I,J,JENV,JENV0,JNENV,JNENV0,K,K2,M,MAXJ,MINJ,NSTEPE 57 LOGICAL*4 LDOACC,LDONE 58 59C 9/25/26 MHP via Claude 60C code changed to consolidate all I/O logical unit numbers, previously 61C scattered across COMMON/LUOUT/, COMMON/LUNUM/, COMMON/ALEX06/, 62C COMMON/ACDPTH/, COMMON/NEWOPAC/, COMMON/LOPAL95/, COMMON/ATMOS2/, 63C COMMON/ALATM03/, COMMON/OPALEOS/, and COMMON/SCV2/, into a single 64C COMMON/IOUNITS/ 65C 9/25/26 MHP via Claude 66C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/, 67C and COMMON/CONST3/ into a single COMMON/CONSTANTS/ 68C MHP 10/91 COMMON BLOCK ADDED FOR OVERSHOOT. 69C 10/12/26 MHP via Claude 70C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 71C namelist variables that the physics routines read from 27 COMMON 72C blocks which also held derived or working members. Those blocks keep 73C only their non-namelist members; five of them -- NEWMX, VMULT, 74C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 75C Derivation-only namelist variables from the same blocks went to 76C COMMON/NMLBYPASS/ instead, not here. 77C 78C Declared identically in all 121 files that hold any of these members. 79C Generated from one member list, so every copy is byte-identical; the 80C linker's 'size of symbol changed' warning is the check that they are. 81C 10/13/26 MHP via Claude: removed the unused COMMON/DPMIX/ declaration 82C (no member referenced here; residue of the obsolete-option removals). 83C MHP 11/91 COMMON BLOCK ADDED FOR OVERSHOOT. 84C 11/91 COMMON BLOCK ADDED FOR THE EXACT LOCATION OF THE BASE OF THE 85C SURFACE CZ AND THE INTEGRATED BURNING RATES THROUGHOUT IT. USED TO 86C REMEMBER THE PREVIOUS "END OF TIMESTEP VALUES" FOR THE NEW 87C "BEGINNING OF TIMESTEP VALUES". 88C MHP 05/02 INCLUDE MASS ACCRETION AND ITS EFFECT ON ABUNDANCES 89 DIMENSION HCOMP(15,JSON),HS2(JSON),HT(JSON),HS1(JSON) ! HD(JSON) 90 DIMENSION RATLI6(JSON),RATLI7(JSON),RATBE9(JSON) 91 DIMENSION HSAV(3,JSON),HR(JSON) 92 DIMENSION NSTEPE(11),YTOL(3),YEST(3),YERR(3),YEXT(3) 93 DATA NSTEPE/2,4,6,8,12,16,24,32,48,64,96/ 94 DATA YTOL/1.0D-6,1.0D-6,1.0D-6/ 95C DATA TLIM/6.0D0/ 96 SAVE 97C DETERMINE LITHIUM AND BERYLLIUM BURNING IN A NON-ROTATING MODEL. 98C THE BURNING RATES DEPEND ON THE LOCAL T AND RHO, THE ABUNDANCE, AND 99C IN THE CASE OF CONVECTION ZONES THE MIXING. 100C 101C IN RADIATIVE REGIONS BURNING IS DONE IMPLICITLY, CORRECTING THE RATES FOR THE 102C CHANGE IN ABUNDANCE AS A FUNCTION OF TIME. 103C 104C IN CONVECTIVE REGIONS THE PROGRAM ITERATES BETWEEN BURNING AND MIXING 105C TO CORRECT THE RATES FOR THE CHANGE IN ABUNDANCE AS A FUNCTION OF TIME. 106C 107C SAVE ORIGINAL ABUNDANCES. 108 DO 10 I = 1,M 109 HSAV(1,I) = HCOMP(13,I) 110 HSAV(2,I) = HCOMP(14,I) 111 HSAV(3,I) = HCOMP(15,I) 112 10 10 CONTINUE 113C THE DEGREE OF LITHIUM BURNING IN A SURFACE CZ DEPENDS SENSITIVELY 114C ON THE TEMPERATURE AT ITS BASE - SO ACCURATELY LOCATING IS IMPORTANT. 115C DETERMINE THE TRUE LOCATION (FX) OF THE BASE OF THE CZ AT THE END 116C OF THE TIMESTEP, AND THE LOCATION OF THE EDGE OF OVERSHOOT REGIONS 117C IF APPLICABLE. 118 IF(JNENV.GT.1.AND.JNENV.LT.M)THEN 119 IF(LROT.AND.LINSTB)THEN 120 DDJ = DELAM(JNENV) - DELRM(JNENV) 121 DDJM1 = DELAM(JNENV-1) - DELRM(JNENV-1) 122 ELSE 123C EVALUATE DEL(AD) - DEL(RAD) AT THE LAST CONVECTIVE POINT AND THE ONE 124C BELOW IT. 125 DDJ = SDEL(3,JNENV)-SDEL(1,JNENV) 126 DDJM1 = SDEL(3,JNENV-1)-SDEL(1,JNENV-1) 127 ENDIF 128C USE LINEAR INTERPOLATION TO FIND THE DISTANCE OF THE TRUE LOCATION 129C OF THE BASE FROM THE ZONE MIDPOINT. IF FX IS NEGATIVE,THEN THE TRUE 130C BASE IS HIGHER; IF IT IS POSITIVE, THE TRUE BASE IS LOWER. 131 FX = MAX(-0.5D0,0.5D0-DDJM1/(DDJM1-DDJ)) 132 FX = MIN(0.5D0,FX) 133 IF(.NOT.LOVSTE)THEN 134 JENV = JNENV 135 JENV0 = JNENV0 136 ELSE 137C STARTING CZ DEPTH 138 IF(RBASE0.EQ.0.0D0)THEN 139C 10/30/26 MHP via Claude 140C code changed to take the radius below shell 1 as zero (the centre). 141C With a fully convective start-of-step model JNENV0 is 1 and this read 142C HRO(0). 143 IF(JNENV0.GT.1)THEN 144 RBASE0 = 0.5D0*(EXP(CLN*HRO(JNENV0)) 145 * * +EXP(CLN*HRO(JNENV0-1))) 146 ELSE 147 RBASE0 = 0.5D0*EXP(CLN*HRO(1)) 148 ENDIF 149 ROS = RBASE0 - PSCAL0 150 DO J = JNENV0-1,1,-1 151 R = EXP(CLN*HRO(J)) 152 IF(R.LT.ROS)THEN 153 JENV0 = J + 1 154 GOTO 11 155 ENDIF 156 END DO 157 JENV0 = 1 158 11 11 CONTINUE 159 ELSE 160 JENV0 = JENVOS 161 ENDIF 162C ENDING CZ DEPTH : DETERMINE OVERSHOOT FROM TRUE CZ BASE. 163 DR = EXP(CLN*HR(JNENV))-EXP(CLN*HR(JNENV-1)) 164 RBASE = 0.5D0*(EXP(CLN*HR(JNENV))+EXP(CLN*HR(JNENV-1))) 165 * * -FX*DR 166 RBASE0 = RBASE 167 ROS = RBASE - PSCAL1 168 DO J = JNENV-1,1,-1 169 R = EXP(CLN*HR(J)) 170 IF(R.LT.ROS)THEN 171 JENV = J + 1 172C 10/30/26 MHP via Claude 173C code changed to take the zone width from the end-of-step model alone. 174C The outer edge was the start-of-step HRO, the inner one the current HR. 175 DR = EXP(CLN*HR(J+1))-R 176 FX = 0.5D0-((ROS-R)/DR) 177 FX = MAX(-0.5D0,FX) 178 FX = MIN(0.5D0,FX) 179 GOTO 12 180 ENDIF 181 END DO 182 JENV = 1 183 12 12 CONTINUE 184 ENDIF 185 ELSE 186 JENV = JNENV 187 JENV0 = JNENV0 188 ENDIF 189 JENVOS = JENV 190C RADIATIVE INTERIOR. 191 MINJ = MIN(JENV,JENV0) 192 MAXJ = MAX(JENV,JENV0) 193 DO 50 I = 1,MINJ-1 194 IF(RBE9(I).LE.1.0D-32 .OR. RBE90(I).LE.1.0D-32)GOTO 60 195 IF(HCOMP(13,I).LT.1.0D-24.AND.HCOMP(14,I).LT.1.0D-24 196 * * .AND.HCOMP(15,I).LT.1.0D-24)GOTO 50 197 IF(HT(I).GT.7.0D0)THEN 198 HCOMP(13,I) = 0.0D0 199 HCOMP(14,I) = 0.0D0 200 HCOMP(15,I) = 0.0D0 201 GOTO 50 202 ENDIF 203 DO 40 J = 1,11 204C PERFORM BURNING IN NSTEPE SMALLER TIMESTEPS. 205C DEFINE FRAC - THE FRACTION OF THE TOTAL TIMESTEP IN EACH SMALL ONE - 206C AND DELT - THE FRACTIONAL TIMESTEP IN YEARS. 207 FRAC = 1.0D0/DFLOAT(NSTEPE(J)) 208 DT = FRAC*DELTS 209 LDONE = .FALSE. 210C RESTORE INITIAL ABUNDANCES. 211 HCOMP(13,I) = HSAV(1,I) 212 HCOMP(14,I) = HSAV(2,I) 213 HCOMP(15,I) = HSAV(3,I) 214C STORE STARTING REACTION RATES. 215 TLI6 = LOG(RLI60(I))-0.5D0*FRAC* 216 * * (LOG(RLI6(I))-LOG(RLI60(I))) 217 TLI7 = LOG(RLI70(I))-0.5D0*FRAC* 218 * * (LOG(RLI7(I))-LOG(RLI70(I))) 219 TBE9 = LOG(RBE90(I))-0.5D0*FRAC* 220 * * (LOG(RBE9(I))-LOG(RBE90(I))) 221 DO 20 K = 1,NSTEPE(J) 222C INCREMENT THE REACTION RATES. 223 TLI6 = TLI6+FRAC*(LOG(RLI6(I))-LOG(RLI60(I))) 224 TLI7 = TLI7+FRAC*(LOG(RLI7(I))-LOG(RLI70(I))) 225 TBE9 = TBE9+FRAC*(LOG(RBE9(I))-LOG(RBE90(I))) 226 RATLI6(I) = DT*EXP(TLI6) 227 RATLI7(I) = DT*EXP(TLI7) 228 RATBE9(I) = DT*EXP(TBE9) 229C THE REACTION RATES ARE OF THE FORM 230C DLI/DT = F(RHO,T,X)*LI =>DLN LI = DT*F(RHO,T,X) 231C SOLVE FOR D LN (SPECIES)/DT AND ZERO OUT IF THE DEPLETION IS TOO HIGH. 232 IF(RATLI6(I).LT.3.0D1)THEN 233 HCOMP(13,I) = HCOMP(13,I)/EXP(RATLI6(I)) 234 ELSE 235 HCOMP(13,I) = 0.0D0 236 ENDIF 237 IF(RATLI7(I).LT.3.0D1)THEN 238 HCOMP(14,I) = HCOMP(14,I)/EXP(RATLI7(I)) 239 ELSE 240 HCOMP(14,I) = 0.0D0 241 ENDIF 242 IF(RATBE9(I).LT.3.0D1)THEN 243 HCOMP(15,I) = HCOMP(15,I)/EXP(RATBE9(I)) 244 ELSE 245 HCOMP(15,I) = 0.0D0 246 ENDIF 247 20 20 CONTINUE 248C STORE ABUNDANCES AS A FUNCTION OF TIMESTEP IN VECTOR YEST. 249 YEST(1) = HCOMP(13,I) 250 YEST(2) = HCOMP(14,I) 251 YEST(3) = HCOMP(15,I) 252 XEST = FRAC**2 253C USE A RATIONAL POLYNOMIAL EXTRAPOLATOR TO EXTRAPOLATE TO THE SOLUTION FOR 254C ZERO TIMESTEP. 255 CALL RATEXT(J,XEST,YEST,YEXT,YERR,3,7) 256C CHECK IF THE SOLUTION HAS CONVERGED. 257 IF(J.GT.1)THEN 258 LDONE = .TRUE. 259 DO 30 K2=1,3 260 IF(YEST(K2).LT.1.0D-24)GOTO 30 261 YERR(K2) = ABS(YERR(K2)/YEST(K2)) 262 IF(YERR(K2).GT.YTOL(K2))LDONE=.FALSE. 263 30 30 CONTINUE 264 IF(LDONE)THEN 265C WRITE(ISHORT,912)I,J,(YEXT(K2)K2=1,3) 266C 912 FORMAT(1X,'CONVERGED',I5,' LEVEL ',I2,' LI6 ',1P, 267C * E12.5,' LI7 ',E12.5,' BE9 ',E12.5) 268 GOTO 45 269 ENDIF 270 ENDIF 271 40 40 CONTINUE 272C IF THE PROGRAM GETS HERE THEN IT FAILED TO CONVERGE TO WITHIN 273C THE SPECIFIED TOLERANCE IN THE MAXIMUM NUMBER OF ITERATIONS. 274 WRITE(ISHORT,911)I,(YERR(K2),K2=1,3) 275 911 911 FORMAT(1X,'***LIBURN CONVERGENCE FAILURE IN SHELL ',I4, 276 * * 'ERRORS '/1P3E10.3) 277 45 45 CONTINUE 278C WRITE NEW ABUNDANCES AND EXIT. 279 HCOMP(13,I)=YEXT(1) 280 IF(HCOMP(13,I).LT.1.0D-24)HCOMP(13,I)=0.0D0 281 HCOMP(14,I)=YEXT(2) 282 IF(HCOMP(14,I).LT.1.0D-24)HCOMP(14,I)=0.0D0 283 HCOMP(15,I)=YEXT(3) 284 IF(HCOMP(15,I).LT.1.0D-24)HCOMP(15,I)=0.0D0 285 50 50 CONTINUE 286 60 60 CONTINUE 287C CONVECTION ZONE. 288C 289C SKIP IF WHOLE CZ IS BELOW THE BURNING THRESHOLD. 290C 11/02/26 MHP via Claude 291C code changed to clear LBLIOLD when the CZ is skipped, so a later step 292C does not reuse a stale end-of-step rate. 293 IF(RBE90(JENV0).LE.1.0D-32.OR.RBE9(JENV).LE.1.0D-32)THEN 294 LBLIOLD = .FALSE. 295 GOTO 200 296 ENDIF 297C FIND RATES AT THE BEGINNING OF THE TIMESTEP (USING THE DEPTH AT THE START). 298 FLI60 = 0.0D0 299 FLI70 = 0.0D0 300 FBE90 = 0.0D0 301 FM0 = 0.0D0 302 DO 65 I = JENV0,M 303 FLI60 = FLI60+HCOMP(13,I)*HS2(I) 304 FLI70 = FLI70+HCOMP(14,I)*HS2(I) 305 FBE90 = FBE90+HCOMP(15,I)*HS2(I) 306 FM0 = FM0 + HS2(I) 307 65 65 CONTINUE 308C 67 CONTINUE 309 FLI60 = FLI60/FM0 310 FLI70 = FLI70/FM0 311 FBE90 = FBE90/FM0 312C 11/02/26 MHP via Claude 313C code changed to test LBLIOLD. BLI6OLD is the log of a small rate, so 314C BLI6OLD.LE.0 was always true and the ELSE branch never ran. 315 IF(.NOT.LBLIOLD)THEN 316C COMPUTE MASS-WEIGHTED AVERAGE RATES AT THE START OF THE STEP. 317 BLI60 = 0.0D0 318 BLI70 = 0.0D0 319 BBE90 = 0.0D0 320 DO 68 I = JENV0,M 321 BLI60 = BLI60 + RLI60(I)*HS2(I) 322 BLI70 = BLI70 + RLI70(I)*HS2(I) 323 BBE90 = BBE90 + RBE90(I)*HS2(I) 324 68 68 CONTINUE 325 BLI60 = LOG(BLI60/FM0) 326 BLI70 = LOG(BLI70/FM0) 327 BBE90 = LOG(BBE90/FM0) 328 ELSE 329C USE THE RATE FROM THE END OF THE PREVIOUS TIMESTEP. 330 BLI60 = BLI6OLD 331 BLI70 = BLI7OLD 332 BBE90 = BBE9OLD 333 ENDIF 334C USE THE LOCATION OF THE TRUE EDGE OF THE CONVECTION ZONE (FX, FOUND AT 335C BEGINNING OF SR) TO ADJUST THE BURNING RATE AND MASS OF THE BOTTOM POINT 336C SUCH THAT IT INCLUDES THE ENTIRE C.Z. 337 HS2SAV = HS2(JENV) 338 IF(JENV.GT.1.AND.JENV.LT.M)THEN 339 HS2(JENV) = HS2(JENV)+FX*(HS1(JENV)-HS1(JENV-1)) 340 RLI6(JENV) = RLI6(JENV)+0.5D0*FX*(RLI6(JENV-1)-RLI6(JENV)) 341 RLI7(JENV) = RLI7(JENV)+0.5D0*FX*(RLI7(JENV-1)-RLI7(JENV)) 342 RBE9(JENV) = RBE9(JENV)+0.5D0*FX*(RBE9(JENV-1)-RBE9(JENV)) 343 ENDIF 344C FIND RATES AT THE END OF THE TIMESTEP (USING THE DEPTH AT THE END). 345C ALSO STORE INITIAL ABUNDANCES(FLI60,FLI70,FBE90). 346C FM IS THE TOTAL MASS IN THE CZ. 347 FLI6 = 0.0D0 348 FLI7 = 0.0D0 349 FBE9 = 0.0D0 350 FM = 0.0D0 351 BLI6 = 0.0D0 352 BLI7 = 0.0D0 353 BBE9 = 0.0D0 354 DO 70 I = JENV,M 355 BLI6 = BLI6 + RLI6(I)*HS2(I) 356 BLI7 = BLI7 + RLI7(I)*HS2(I) 357 BBE9 = BBE9 + RBE9(I)*HS2(I) 358 FLI6 = FLI6+HCOMP(13,I)*HS2(I) 359 FLI7 = FLI7+HCOMP(14,I)*HS2(I) 360 FBE9 = FBE9+HCOMP(15,I)*HS2(I) 361 FM = FM + HS2(I) 362 70 70 CONTINUE 363C 75 CONTINUE 364 FLI6 = FLI6/FM 365 FLI7 = FLI7/FM 366 FBE9 = FBE9/FM 367C MASS-WEIGHTED AVERAGE RATES AT THE END OF THE STEP. 368 BLI6 = LOG(BLI6/FM) 369 BLI7 = LOG(BLI7/FM) 370 BBE9 = LOG(BBE9/FM) 371C RESTORE ORIGINAL MASS TO BASE OF C.Z. 372 HS2(JENV) = HS2SAV 373C FOR THE STARTING ABUNDANCE, USE THE MIXED ABUNDANCE FOR THE 374C DEEPER CZ. 375 IF(JENV.LT.JENV0)THEN 376 FLI60 = FLI6 377 FLI70 = FLI7 378 FBE90 = FBE9 379 ENDIF 380 DO 120 J = 1,11 381C PERFORM BURNING IN NSTEPE SMALLER TIMESTEPS. 382C DEFINE FRAC - THE FRACTION OF THE TOTAL TIMESTEP IN EACH SMALL ONE - 383C AND DELT - THE FRACTIONAL TIMESTEP IN YEARS. 384 FRAC = 1.0D0/DFLOAT(NSTEPE(J)) 385 DT = FRAC*DELTS 386 LDONE = .FALSE. 387C INITIALIZE ABUNDANCES. 388 FLI6 = FLI60 389 FLI7 = FLI70 390 FBE9 = FBE90 391C INITIALIZE RATES. 392 CZLI6 = BLI60-0.5D0*FRAC*(BLI6-BLI60) 393 CZLI7 = BLI70-0.5D0*FRAC*(BLI7-BLI70) 394 CZBE9 = BBE90-0.5D0*FRAC*(BBE9-BBE90) 395C BURN IN NSTEPE SMALL TIMESTEPS. 396C MHP 05/02 ADDED THE IMPACT OF MASS ACCRETION ON LIGHT 397C ELEMENT ABUNDANCES. INITIALIZE ARRAYS; SPECIES ARE 398C ASSUMED TO EXPERIENCE BURNING FOR ON AVERAGE 1/2 OF 399C THE STEP WHERE THEY ARE INITIALLY ACCRETED AND ARE 400C THEN FULLY BURNED. 401 IF(LMDOT.AND.DMDT0.GT.0.0D0)THEN 402 LDOACC = .TRUE. 403 FLI6A = 0.0D0 404 FLI7A = 0.0D0 405 FBE9A = 0.0D0 406 ELSE 407 LDOACC = .FALSE. 408 ENDIF 409 DO 100 K = 1,NSTEPE(J) 410C INCREMENT RATES BY LINEAR INTERPOLATION IN THE LOG. 411 CZLI6 = CZLI6 + FRAC*(BLI6-BLI60) 412 CZLI7 = CZLI7 + FRAC*(BLI7-BLI70) 413 CZBE9 = CZBE9 + FRAC*(BBE9-BBE90) 414 DLI6 = DT*EXP(CZLI6) 415 DLI7 = DT*EXP(CZLI7) 416 DBE9 = DT*EXP(CZBE9) 417 IF(DLI6.LT.3.0D1 .AND. FLI6.GT.1.0D-24)THEN 418 FLI6 = FLI6/EXP(DLI6) 419 ELSE 420 FLI6 = 0.0D0 421 ENDIF 422 IF(FLI6.LT.1.0D-24)FLI6 = 0.0D0 423 IF(DLI7.LT.3.0D1 .AND. FLI7.GT.1.0D-24)THEN 424 FLI7 = FLI7/EXP(DLI7) 425 ELSE 426 FLI7 = 0.0D0 427 ENDIF 428 IF(FLI7.LT.1.0D-24)FLI7 = 0.0D0 429 IF(DBE9.LT.3.0D1 .AND. FBE9.GT.1.0D-24)THEN 430 FBE9 = FBE9/EXP(DBE9) 431 ELSE 432 FBE9 = 0.0D0 433 ENDIF 434 IF(FBE9.LT.1.0D-24)FBE9 = 0.0D0 435C MHP 05/02 MASS ACCRETION. THE NET LITHIUM ADDED IN 436C A GIVEN STEP IS XXXADD WHILE THE LITHIUM ADDED IN PRIOR 437C STEPS IS XXXA; THE TWO ARE ADDED AFTER BURNING. 438C N.B. SINCE THIS IS A HALF-LIFE PROBLEM IT IS OK TO BURN 439C THE COMPONENTS SEPARATELY. 440 IF(LDOACC)THEN 441 FLI6ADD = COMPACC(13)*FRAC/EXP(0.5D0*DLI6) 442 FLI6A = FLI6A/EXP(DLI6) + FLI6ADD 443 FLI7ADD = COMPACC(14)*FRAC/EXP(0.5D0*DLI7) 444 FLI7A = FLI7A/EXP(DLI7) + FLI7ADD 445 FBE9ADD = COMPACC(15)*FRAC/EXP(0.5D0*DBE9) 446 FBE9A = FBE9A/EXP(DBE9) + FBE9ADD 447 ENDIF 448 100 100 CONTINUE 449C DO A MASS-WEIGHTED AVERAGE OF THE ORIGINAL LIGHT 450C ELEMENT CONTENT AND THE NET AMOUNT ADDED. 451C FMASSACC = DMDT*DT/ORIGINAL CZ MASS 452C NOTE: THIS FORMULATION ASSUMES THAT ALL ACCRETED MATTER 453C LANDS IN A CZ WITH THE END-OF-TIMESTEP DEPTH. 454C NOT STRICTLY TRUE, BUT NOT A BAD APPROXIMATION EITHER. 455 IF(LDOACC)THEN 456 WRITE(*,913)FLI60,FLI6,FLI6A,FLI70,FLI7,FLI7A,FMASSACC 457 913 913 FORMAT(1P7E12.3) 458 FLI6 = (FLI6*FM+FLI6A*FMASSACC)/(FM+FMASSACC) 459 FLI7 = (FLI7*FM+FLI7A*FMASSACC)/(FM+FMASSACC) 460 FBE9 = (FBE9*FM+FBE9A*FMASSACC)/(FM+FMASSACC) 461 ENDIF 462C STORE THE ABUNDANCE AS A FUNCTION OF THE TIME STEP. 463C STORE ABUNDANCES AS A FUNCTION OF TIMESTEP IN VECTOR YEST. 464 YEST(1) = FLI6 465 YEST(2) = FLI7 466 YEST(3) = FBE9 467 XEST = FRAC**2 468C USE A RATIONAL POLYNOMIAL EXTRAPOLATOR TO EXTRAPOLATE TO THE SOLUTION FOR 469C ZERO TIMESTEP. 470 CALL RATEXT(J,XEST,YEST,YEXT,YERR,3,7) 471C CHECK IF THE SOLUTION HAS CONVERGED. 472 IF(J.GT.1)THEN 473 LDONE = .TRUE. 474 DO 110 K2=1,3 475 IF(YEST(K2).LT.1.0D-24)GOTO 110 476 YERR(K2) = ABS(YERR(K2)/YEST(K2)) 477 IF(YERR(K2).GT.YTOL(K2))LDONE=.FALSE. 478 110 110 CONTINUE 479 IF(LDONE)THEN 480C WRITE(ISHORT,912)JENV,J,(YEXT(K2),YERR(K2),K2=1,3) 481 GOTO 125 482 ENDIF 483 ENDIF 484 120 120 CONTINUE 485C IF THE PROGRAM GETS HERE THEN IT FAILED TO CONVERGE TO WITHIN 486C THE SPECIFIED TOLERANCE IN THE MAXIMUM NUMBER OF ITERATIONS. 487 WRITE(ISHORT,911)JENV,(YERR(K2),K2=1,3) 488 125 125 CONTINUE 489C WRITE NEW ABUNDANCES AND EXIT. 490 FLI6 = YEXT(1) 491 IF(FLI6.LT.1.0D-24)FLI6=0.0D0 492 FLI7 = YEXT(2) 493 IF(FLI7.LT.1.0D-24)FLI7=0.0D0 494 FBE9 = YEXT(3) 495 IF(FBE9.LT.1.0D-24)FBE9=0.0D0 496 DO 130 I = MAXJ,M 497 HCOMP(13,I) = FLI6 498 HCOMP(14,I) = FLI7 499 HCOMP(15,I) = FBE9 500 130 130 CONTINUE 501C STORE ENDING RATE FOR USE AT THE BEGINNING OF THE NEXT STEP. 502 BLI6OLD = BLI6 503 BLI7OLD = BLI7 504 BBE9OLD = BBE9 505 LBLIOLD = .TRUE. 506C NOW SOLVE FOR ABUNDANCES IN THE REGION WHICH BEGAN CONVECTIVE AND 507C ENDED RADIATIVE. 508 IF(JENV.LE.JENV0)GOTO 200 509C FIND STARTING AND ENDING LOCATION IN MASS OF THE CZ BASE, AND 510C ASSUME THAT THE FRACTION OF TIME SPENT IN THE CZ IS PROPORTIONAL TO 511C THE LOCATION IN MASS (I.E. A POINT 1/3 FROM THE OLD TO THE NEW CZ 512C BASE IN MASS SPENT 1/3 OF THE TIME IN THE CZ AND 2/3 OUT OF IT). 513 IF(JENV0.GT.1)THEN 514 HSBEG = 0.5D0*(HS1(JENV0)+HS1(JENV0-1)) 515 ELSE 516 HSBEG = 0.0D0 517 ENDIF 518 HSEND = 0.5D0*(HS1(JENV)+HS1(JENV-1)) 519 DMASS = HSBEG - HSEND 520 DO 140 I = JENV0,JENV-1 521C MHP 9/91 CHANGE TO AVOID DIVISION BY ZERO. 522C SKIP IF SHELL TEMPERATURE DROPS BELOW BURNING THRESHOLD. 523 IF(RBE9(I).LE.1.0D-32)GOTO 145 524 FRAD = (HS1(I)-HSBEG)/DMASS 525C USE FRAD*RADIATIVE RATE AND (1-FRAD)*CONVECTIVE RATE. 526 DLI6 = DELTS*EXP(FRAD*LOG(RLI6(I))+(1.0D0-FRAD)*BLI60) 527 DLI7 = DELTS*EXP(FRAD*LOG(RLI7(I))+(1.0D0-FRAD)*BLI70) 528 DBE9 = DELTS*EXP(FRAD*LOG(RBE9(I))+(1.0D0-FRAD)*BBE90) 529C***REMEMBER TO ADD FAILSAFES FOR LARGE DEPLETION*** 530C KC 2025-05-31 PREVENT FLOATING POINT EXCEPTION 531C HCOMP(13,I) = HCOMP(13,I)/EXP(DLI6) 532 CALL SAFEDIVEXP(HCOMP(13,I),DLI6) 533 IF(HCOMP(13,I).LT.1.0D-24)HCOMP(13,I)=0.0D0 534C HCOMP(14,I) = HCOMP(14,I)/EXP(DLI7) 535 CALL SAFEDIVEXP(HCOMP(14,I),DLI7) 536 IF(HCOMP(14,I).LT.1.0D-24)HCOMP(14,I)=0.0D0 537 HCOMP(15,I) = HCOMP(15,I)/EXP(DBE9) 538 IF(HCOMP(15,I).LT.1.0D-24)HCOMP(15,I)=0.0D0 539 140 140 CONTINUE 540 145 145 CONTINUE 541 200 200 CONTINUE 542 RETURN 543 END 544 545C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 546C KC 2025-05-31 SAFEDIVEXP 547C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 548 SUBROUTINE SAFEDIVEXP(NUMERATOR, EXPONENT) 549C 10/20/26 MHP via Claude 550C code changed to replace implicit typing with IMPLICIT NONE. 551 IMPLICIT NONE 552 REAL*8 NUMERATOR, EXPONENT 553 554 IF (EXPONENT .LT. 709.7827D0) THEN 555 NUMERATOR = NUMERATOR/EXP(EXPONENT) 556 ELSE 557 NUMERATOR = 0.0D0 558 ENDIF 559 560 RETURN 561 END SUBROUTINE SAFEDIVEXP

1"""Python version of liburn.f (YREC development, 7be4133).""" 2 3import math 4 5import numpy as np 6 7from pyrec import constants_mod 8from pyrec import deuter_mod 9from pyrec import dimensions_mod 10from pyrec import iounits_mod 11from pyrec import liov_mod 12from pyrec import mdphy_mod 13from pyrec import newrat_mod 14from pyrec import nmlpar_mod 15from pyrec import oldmod_mod 16from pyrec import oldrat_mod 17from pyrec import prevcz_mod 18from pyrec import scrtch_mod 19from pyrec.ratext import ratext 20 21# DATA: the numbers of substeps tried, and the relative tolerances 22# (Fortran lines 93–94). 23NSTEPE = [0, 2, 4, 6, 8, 12, 16, 24, 32, 48, 64, 96] 24YTOL = [0.0, 1.0e-6, 1.0e-6, 1.0e-6] 25 26 27def _converged(yest: np.ndarray, yerr: np.ndarray) -> bool: 28 """Loops 30 and 110: relative errors of the extrapolation within YTOL 29 (Fortran lines 257–263 and 472–478). 30 31 Abundances below 1e-24 are not tested. ``yerr`` is made relative in 32 place, as in the Fortran. 33 34 """ 35 36 ldone = True 37 for k2 in range(1, 4): 38 if yest[k2] < 1.0e-24: 39 continue 40 yerr[k2] = abs(yerr[k2] / yest[k2]) 41 if yerr[k2] > YTOL[k2]: 42 ldone = False 43 return ldone 44 45 46def _write_911(i: int, yerr: np.ndarray) -> None: 47 """Format 911: convergence failure of the extrapolation (Fortran lines 48 274–276 and 487). 49 50 """ 51 52 iounits_mod.write(iounits_mod.ishort, 53 f" ***LIBURN CONVERGENCE FAILURE IN SHELL {i:4d}ERRORS ") 54 iounits_mod.write(iounits_mod.ishort, 55 "".join(f"{yerr[k2]:10.3E}" for k2 in range(1, 4))) 56 57 58def liburn(delts: float, hcomp: np.ndarray, hr: np.ndarray, hs1: np.ndarray, hs2: np.ndarray, ht: np.ndarray, 59 jnenv: int, jnenv0: int, m: int) -> None: 60 """Burn Li6, Li7 and Be9 over one time step (non-rotating mixing). 61 62 The radiative interior is burned shell by shell and the surface 63 convection zone as one mixed reservoir, the rates interpolated in their 64 logs between the start (``oldrat_mod``) and end (``newrat_mod``) of the 65 step, with substeps and rational extrapolation (RATEXT) to zero step 66 size. Points left behind by a retreating zone are burned once with a 67 blended rate. 68 69 Translates: liburn.f#LIBURN (7be4133) 70 71 Parameters 72 ---------- 73 delts : float 74 The time step (main.f passes seconds; see the page's question on 75 the comments' units). 76 hcomp : ndarray, shape (16, JSON + 1) 77 Composition; rows 13, 14, 15 (Li6, Li7, Be9) updated in place. 78 hr : ndarray 79 log r at each point. 80 hs1, hs2 : ndarray 81 Mass inside each point, and each point's mass (g). ``hs2[jenv]`` 82 is changed and restored. 83 ht : ndarray 84 log T at each point. 85 jnenv, jnenv0 : int 86 The base of the surface convection zone at the end and the start 87 of the step. 88 m : int 89 Number of points. 90 91 Also changes ``newrat_mod`` (the bottom point's rates, a suspected 92 bug on the page) and 93 ``prevcz_mod`` (``rbase0``, ``bli6old``, ``bli7old``, ``bbe9old``, 94 ``lbliold``, ``jenvos``). Nothing is returned. 95 96 """ 97 98 # Module data (Fortran lines 12–45). 99 json = dimensions_mod.json 100 cln = constants_mod.cln 101 rli6 = newrat_mod.rli6 102 rli7 = newrat_mod.rli7 103 rbe9 = newrat_mod.rbe9 104 rli60 = oldrat_mod.rli60 105 rli70 = oldrat_mod.rli70 106 rbe90 = oldrat_mod.rbe90 107 hro = oldmod_mod.hro 108 sdel = scrtch_mod.sdel 109 110 # Local arrays (Fortran lines 89–96). SAVE in the Fortran: no local is 111 # read before it is set in a call. 112 hsav = np.zeros((4, json + 1)) 113 ratli6 = np.zeros(json + 1) 114 ratli7 = np.zeros(json + 1) 115 ratbe9 = np.zeros(json + 1) 116 yest = np.zeros(4) 117 yerr = np.zeros(4) 118 yext = np.zeros(4) 119 fx = 0.0 120 121 # Save the original abundances (loop 10, Fortran lines 107–112). 122 for i in range(1, m + 1): # loop 10 123 hsav[1, i] = hcomp[13, i] 124 hsav[2, i] = hcomp[14, i] 125 hsav[3, i] = hcomp[15, i] 126 127 # The base of the convection zone: the offset fx (in zones) where 128 # del_ad - del_rad crosses zero, between jnenv - 1 and jnenv (Fortran 129 # lines 113–189). 130 if 1 < jnenv < m: 131 if nmlpar_mod.lrot and nmlpar_mod.linstb: 132 ddj = mdphy_mod.delam[jnenv] - mdphy_mod.delrm[jnenv] 133 ddjm1 = mdphy_mod.delam[jnenv - 1] - mdphy_mod.delrm[jnenv - 1] 134 else: 135 ddj = sdel[3, jnenv] - sdel[1, jnenv] 136 ddjm1 = sdel[3, jnenv - 1] - sdel[1, jnenv - 1] 137 fx = max(-0.5, 0.5 - ddjm1 / (ddjm1 - ddj)) 138 fx = min(0.5, fx) 139 if not nmlpar_mod.lovste: 140 jenv = jnenv 141 jenv0 = jnenv0 142 else: 143 # Envelope overshoot: move both bases down by PSCAL0, PSCAL1; 144 # the start of the step (Fortran lines 137–161). 145 if prevcz_mod.rbase0 == 0.0: 146 # The radius below shell 1 is zero (the centre). 147 if jnenv0 > 1: 148 prevcz_mod.rbase0 = 0.5 * (math.exp(cln * hro[jnenv0]) 149 + math.exp(cln * hro[jnenv0 - 1])) 150 else: 151 prevcz_mod.rbase0 = 0.5 * math.exp(cln * hro[1]) 152 ros = prevcz_mod.rbase0 - liov_mod.pscal0 153 jenv0 = 1 154 for j in range(jnenv0 - 1, 0, -1): 155 r = math.exp(cln * hro[j]) 156 if r < ros: 157 jenv0 = j + 1 158 break # to label 11 159 else: 160 jenv0 = prevcz_mod.jenvos 161 162 # The end of the step, from the true base (Fortran lines 163 # 162–183). 164 dr = math.exp(cln * hr[jnenv]) - math.exp(cln * hr[jnenv - 1]) 165 rbase = (0.5 * (math.exp(cln * hr[jnenv]) 166 + math.exp(cln * hr[jnenv - 1])) - fx * dr) 167 prevcz_mod.rbase0 = rbase 168 ros = rbase - liov_mod.pscal1 169 jenv = 1 170 for j in range(jnenv - 1, 0, -1): 171 r = math.exp(cln * hr[j]) 172 if r < ros: 173 jenv = j + 1 174 dr = math.exp(cln * hr[j + 1]) - r 175 fx = 0.5 - ((ros - r) / dr) 176 fx = max(-0.5, fx) 177 fx = min(0.5, fx) 178 break # to label 12 179 else: 180 jenv = jnenv 181 jenv0 = jnenv0 182 prevcz_mod.jenvos = jenv 183 184 # The radiative interior, below both bases (loop 50, Fortran lines 185 # 190–286). 186 minj = min(jenv, jenv0) 187 maxj = max(jenv, jenv0) 188 for i in range(1, minj): 189 # Question (see the page): leaves the loop, not the point. 190 if rbe9[i] <= 1.0e-32 or rbe90[i] <= 1.0e-32: 191 break # to label 60 192 if (hcomp[13, i] < 1.0e-24 and hcomp[14, i] < 1.0e-24 193 and hcomp[15, i] < 1.0e-24): 194 continue 195 if ht[i] > 7.0: 196 hcomp[13, i] = 0.0 197 hcomp[14, i] = 0.0 198 hcomp[15, i] = 0.0 199 continue 200 201 # Substeps of 1/2, 1/4, ... 1/96 of the step, extrapolated to zero 202 # (loop 40, Fortran lines 203–271). 203 ldone = False 204 for j in range(1, 12): # loop 40 205 frac = 1.0 / float(NSTEPE[j]) 206 dt = frac * delts 207 hcomp[13, i] = hsav[1, i] 208 hcomp[14, i] = hsav[2, i] 209 hcomp[15, i] = hsav[3, i] 210 211 # log rates, started half a substep before the step (Fortran 212 # lines 214–220). 213 tli6 = (math.log(rli60[i]) 214 - 0.5 * frac * (math.log(rli6[i]) - math.log(rli60[i]))) 215 tli7 = (math.log(rli70[i]) 216 - 0.5 * frac * (math.log(rli7[i]) - math.log(rli70[i]))) 217 tbe9 = (math.log(rbe90[i]) 218 - 0.5 * frac * (math.log(rbe9[i]) - math.log(rbe90[i]))) 219 220 # Burn, the rates rising linearly in their logs (loop 20, 221 # Fortran lines 221–247). 222 for k in range(1, NSTEPE[j] + 1): # loop 20 223 tli6 = tli6 + frac * (math.log(rli6[i]) - math.log(rli60[i])) 224 tli7 = tli7 + frac * (math.log(rli7[i]) - math.log(rli70[i])) 225 tbe9 = tbe9 + frac * (math.log(rbe9[i]) - math.log(rbe90[i])) 226 ratli6[i] = dt * math.exp(tli6) 227 ratli7[i] = dt * math.exp(tli7) 228 ratbe9[i] = dt * math.exp(tbe9) 229 if ratli6[i] < 3.0e1: 230 hcomp[13, i] = hcomp[13, i] / math.exp(ratli6[i]) 231 else: 232 hcomp[13, i] = 0.0 233 if ratli7[i] < 3.0e1: 234 hcomp[14, i] = hcomp[14, i] / math.exp(ratli7[i]) 235 else: 236 hcomp[14, i] = 0.0 237 if ratbe9[i] < 3.0e1: 238 hcomp[15, i] = hcomp[15, i] / math.exp(ratbe9[i]) 239 else: 240 hcomp[15, i] = 0.0 241 242 # Extrapolate to zero step and test (Fortran lines 248–270). 243 yest[1] = hcomp[13, i] 244 yest[2] = hcomp[14, i] 245 yest[3] = hcomp[15, i] 246 xest = frac ** 2 247 ratext(j, xest, yest, yext, yerr, 3, 7) 248 if j > 1: 249 ldone = _converged(yest, yerr) 250 if ldone: 251 break # to label 45 252 253 # Convergence failure (Fortran lines 272–276); label 45, the new 254 # abundances (Fortran lines 277–284). 255 if not ldone: 256 _write_911(i, yerr) 257 hcomp[13, i] = yext[1] 258 if hcomp[13, i] < 1.0e-24: 259 hcomp[13, i] = 0.0 260 hcomp[14, i] = yext[2] 261 if hcomp[14, i] < 1.0e-24: 262 hcomp[14, i] = 0.0 263 hcomp[15, i] = yext[3] 264 if hcomp[15, i] < 1.0e-24: 265 hcomp[15, i] = 0.0 266 267 # label 60: the convection zone; skipped if it does not burn (Fortran 268 # lines 286–296). 269 if rbe90[jenv0] <= 1.0e-32 or rbe9[jenv] <= 1.0e-32: 270 prevcz_mod.lbliold = False 271 return # label 200 272 273 # Mean abundances over the zone at the start of the step (loop 65, 274 # Fortran lines 297–311). 275 fli60 = 0.0 276 fli70 = 0.0 277 fbe90 = 0.0 278 fm0 = 0.0 279 for i in range(jenv0, m + 1): 280 fli60 = fli60 + hcomp[13, i] * hs2[i] 281 fli70 = fli70 + hcomp[14, i] * hs2[i] 282 fbe90 = fbe90 + hcomp[15, i] * hs2[i] 283 fm0 = fm0 + hs2[i] 284 fli60 = fli60 / fm0 285 fli70 = fli70 / fm0 286 fbe90 = fbe90 / fm0 287 288 # log of the zone's mean start rates, or the previous step's end ones 289 # (Fortran lines 312–333). 290 if not prevcz_mod.lbliold: 291 bli60 = 0.0 292 bli70 = 0.0 293 bbe90 = 0.0 294 for i in range(jenv0, m + 1): # loop 68 295 bli60 = bli60 + rli60[i] * hs2[i] 296 bli70 = bli70 + rli70[i] * hs2[i] 297 bbe90 = bbe90 + rbe90[i] * hs2[i] 298 bli60 = math.log(bli60 / fm0) 299 bli70 = math.log(bli70 / fm0) 300 bbe90 = math.log(bbe90 / fm0) 301 else: 302 bli60 = prevcz_mod.bli6old 303 bli70 = prevcz_mod.bli7old 304 bbe90 = prevcz_mod.bbe9old 305 306 # The bottom point's mass and rates adjusted for the offset fx 307 # (Fortran lines 334–343). 308 # Suspected bug (see the page): the rates are not restored afterwards. 309 hs2sav = hs2[jenv] 310 if 1 < jenv < m: 311 hs2[jenv] = hs2[jenv] + fx * (hs1[jenv] - hs1[jenv - 1]) 312 rli6[jenv] = rli6[jenv] + 0.5 * fx * (rli6[jenv - 1] - rli6[jenv]) 313 rli7[jenv] = rli7[jenv] + 0.5 * fx * (rli7[jenv - 1] - rli7[jenv]) 314 rbe9[jenv] = rbe9[jenv] + 0.5 * fx * (rbe9[jenv - 1] - rbe9[jenv]) 315 316 # Mean abundances and rates over the zone at the end, then the mass 317 # restored (loop 70, Fortran lines 344–372). 318 fli6 = 0.0 319 fli7 = 0.0 320 fbe9 = 0.0 321 fm = 0.0 322 bli6 = 0.0 323 bli7 = 0.0 324 bbe9 = 0.0 325 for i in range(jenv, m + 1): 326 bli6 = bli6 + rli6[i] * hs2[i] 327 bli7 = bli7 + rli7[i] * hs2[i] 328 bbe9 = bbe9 + rbe9[i] * hs2[i] 329 fli6 = fli6 + hcomp[13, i] * hs2[i] 330 fli7 = fli7 + hcomp[14, i] * hs2[i] 331 fbe9 = fbe9 + hcomp[15, i] * hs2[i] 332 fm = fm + hs2[i] 333 fli6 = fli6 / fm 334 fli7 = fli7 / fm 335 fbe9 = fbe9 / fm 336 bli6 = math.log(bli6 / fm) 337 bli7 = math.log(bli7 / fm) 338 bbe9 = math.log(bbe9 / fm) 339 hs2[jenv] = hs2sav 340 341 # A deepened zone starts from the end-of-step mean (Fortran lines 342 # 373–379). 343 if jenv < jenv0: 344 fli60 = fli6 345 fli70 = fli7 346 fbe90 = fbe9 347 348 # Substeps and extrapolation for the zone (loop 120, Fortran lines 349 # 380–484). 350 ldone = False 351 for j in range(1, 12): # loop 120 352 frac = 1.0 / float(NSTEPE[j]) 353 dt = frac * delts 354 fli6 = fli60 355 fli7 = fli70 356 fbe9 = fbe90 357 czli6 = bli60 - 0.5 * frac * (bli6 - bli60) 358 czli7 = bli70 - 0.5 * frac * (bli7 - bli70) 359 czbe9 = bbe90 - 0.5 * frac * (bbe9 - bbe90) 360 361 # Accretion: material of composition COMPACC added each substep 362 # (Fortran lines 396–408). 363 if nmlpar_mod.lmdot and nmlpar_mod.dmdt0 > 0.0: 364 ldoacc = True 365 fli6a = 0.0 366 fli7a = 0.0 367 fbe9a = 0.0 368 else: 369 ldoacc = False 370 371 # Burn in substeps (loop 100, Fortran lines 409–448). 372 for k in range(1, NSTEPE[j] + 1): # loop 100 373 czli6 = czli6 + frac * (bli6 - bli60) 374 czli7 = czli7 + frac * (bli7 - bli70) 375 czbe9 = czbe9 + frac * (bbe9 - bbe90) 376 dli6 = dt * math.exp(czli6) 377 dli7 = dt * math.exp(czli7) 378 dbe9 = dt * math.exp(czbe9) 379 if dli6 < 3.0e1 and fli6 > 1.0e-24: 380 fli6 = fli6 / math.exp(dli6) 381 else: 382 fli6 = 0.0 383 if fli6 < 1.0e-24: 384 fli6 = 0.0 385 if dli7 < 3.0e1 and fli7 > 1.0e-24: 386 fli7 = fli7 / math.exp(dli7) 387 else: 388 fli7 = 0.0 389 if fli7 < 1.0e-24: 390 fli7 = 0.0 391 if dbe9 < 3.0e1 and fbe9 > 1.0e-24: 392 fbe9 = fbe9 / math.exp(dbe9) 393 else: 394 fbe9 = 0.0 395 if fbe9 < 1.0e-24: 396 fbe9 = 0.0 397 if ldoacc: 398 # Accreted material burns half a substep on arrival (Fortran 399 # lines 435–447). 400 fli6add = nmlpar_mod.compacc[13] * frac / math.exp(0.5 * dli6) 401 fli6a = fli6a / math.exp(dli6) + fli6add 402 fli7add = nmlpar_mod.compacc[14] * frac / math.exp(0.5 * dli7) 403 fli7a = fli7a / math.exp(dli7) + fli7add 404 fbe9add = nmlpar_mod.compacc[15] * frac / math.exp(0.5 * dbe9) 405 fbe9a = fbe9a / math.exp(dbe9) + fbe9add 406 407 # Mass-weighted mean with the accreted matter (Fortran lines 408 # 449–461). 409 if ldoacc: 410 fmassacc = deuter_mod.fmassacc 411 print("".join(f"{v:12.3E}" for v in (fli60, fli6, fli6a, fli70, 412 fli7, fli7a, fmassacc))) 413 fli6 = (fli6 * fm + fli6a * fmassacc) / (fm + fmassacc) 414 fli7 = (fli7 * fm + fli7a * fmassacc) / (fm + fmassacc) 415 fbe9 = (fbe9 * fm + fbe9a * fmassacc) / (fm + fmassacc) 416 417 # Extrapolate to zero step and test (Fortran lines 462–483). 418 yest[1] = fli6 419 yest[2] = fli7 420 yest[3] = fbe9 421 xest = frac ** 2 422 ratext(j, xest, yest, yext, yerr, 3, 7) 423 if j > 1: 424 ldone = _converged(yest, yerr) 425 if ldone: 426 break # to label 125 427 428 # Convergence failure (Fortran lines 485–487); label 125, the new 429 # abundances (Fortran lines 488–495). 430 if not ldone: 431 _write_911(jenv, yerr) 432 fli6 = yext[1] 433 if fli6 < 1.0e-24: 434 fli6 = 0.0 435 fli7 = yext[2] 436 if fli7 < 1.0e-24: 437 fli7 = 0.0 438 fbe9 = yext[3] 439 if fbe9 < 1.0e-24: 440 fbe9 = 0.0 441 442 # The zone's new abundances (loop 130, Fortran lines 496–500). 443 # Suspected bug (see the page): from maxj, so points newly in a 444 # deepened zone keep their old abundances. 445 for i in range(maxj, m + 1): # loop 130 446 hcomp[13, i] = fli6 447 hcomp[14, i] = fli7 448 hcomp[15, i] = fbe9 449 450 # The end rates, for the start of the next step (Fortran lines 451 # 501–505). 452 prevcz_mod.bli6old = bli6 453 prevcz_mod.bli7old = bli7 454 prevcz_mod.bbe9old = bbe9 455 prevcz_mod.lbliold = True 456 457 # The region left by a retreating zone: burned once, blended rate 458 # (Fortran lines 506–540). 459 if jenv <= jenv0: 460 return # label 200 461 if jenv0 > 1: 462 hsbeg = 0.5 * (hs1[jenv0] + hs1[jenv0 - 1]) 463 else: 464 hsbeg = 0.0 465 hsend = 0.5 * (hs1[jenv] + hs1[jenv - 1]) 466 # Suspected bug (see the page): dmass < 0 here, so frad lies in (-1, 0), 467 # and the weights are the other way round from the comment. 468 dmass = hsbeg - hsend 469 for i in range(jenv0, jenv): # loop 140 470 if rbe9[i] <= 1.0e-32: 471 break # to label 145 472 frad = (hs1[i] - hsbeg) / dmass 473 dli6 = delts * math.exp(frad * math.log(rli6[i]) + (1.0 - frad) * bli60) 474 dli7 = delts * math.exp(frad * math.log(rli7[i]) + (1.0 - frad) * bli70) 475 dbe9 = delts * math.exp(frad * math.log(rbe9[i]) + (1.0 - frad) * bbe90) 476 hcomp[13, i] = safedivexp(hcomp[13, i], dli6) 477 if hcomp[13, i] < 1.0e-24: 478 hcomp[13, i] = 0.0 479 hcomp[14, i] = safedivexp(hcomp[14, i], dli7) 480 if hcomp[14, i] < 1.0e-24: 481 hcomp[14, i] = 0.0 482 # Suspected bug (see the page): Be9 is not guarded by SAFEDIVEXP 483 # (math.exp raises OverflowError where the Fortran's EXP overflows). 484 hcomp[15, i] = hcomp[15, i] / math.exp(dbe9) 485 if hcomp[15, i] < 1.0e-24: 486 hcomp[15, i] = 0.0 487 488 489def safedivexp(numerator: float, exponent: float) -> float: 490 """KC's guard (2025-05-31): numerator / exp(exponent), or 0 on overflow. 491 492 Translates: liburn.f#SAFEDIVEXP (7be4133) 493 494 Returns ``numerator``, the Fortran's in-out argument. 495 496 """ 497 498 # Fortran lines 554–558. 499 if exponent < 709.7827: 500 numerator = numerator / math.exp(exponent) 501 else: 502 numerator = 0.0 503 return numerator

liburn.f

Purpose

How it fits

The method

Design

Suspected bugs, checked by running the code

1. (harmless in the runs) The sign of FRAD.

2. (harmless in the runs) Points newly in a deepened zone are not burned.

3. (never occurred) ⁹Be is not guarded against overflow.

4. (harmless in the runs) The bottom point's rates are left adjusted.

Definitions

subroutineLIBURN

subroutineSAFEDIVEXP

Change log

1. (harmless in the runs) The sign of `FRAD`.