kemcom.f

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

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

Variables on this page (8)

abstol, cmin, hcompp, ishort, json, kemmax, reltol, tcut

Purpose

KEMCOM solves the nuclear burning network for one shell, or one convection zone, over a time step: it finds the new mass fractions of H, ³He, ⁴He, ¹²C, ¹³C, ¹⁴N and ¹⁶O simultaneously and implicitly, by Newton iteration, from the reaction rates computed in ENGEB (its comment: "Don VandenBerg's notes explain what the SR is doing"). ¹⁸O is updated explicitly afterwards.

How it fits

Who calls it: MIX and ROTMIX, for each radiative shell (IBEGIN = IEND) and each convection zone (IBEGIN < IEND).
Inputs: HT ($\log T$), the rates HR1–HR8, HR10–HR12 and the branching ratio HF1 (per Gyr per amu, without the composition factors, as its comment says), HS2 (shell masses), DDAGE (the step in years); the abundances at the start of the step HCOMPP from the module OLDMOD_MOD.
Output: the new HCOMP (elements 1–7, 9 and 11; element 3, $Z$, is either reset or set to its zone mean, see step 4) for shells IBEGIN..IEND.
Calls: SIMEQC, the linear solver for the 7×8 Newton system.
Modules (COMMON blocks until rp_1): NMLPAR_MOD for the settings CMIN, ABSTOL, RELTOL, KEMMAX (formerly /BURTOL/, set by PARMIN) and TCUT(1), the lowest burning temperature (formerly /CTLIM/); OLDMOD_MOD (the previous model: HCOMPP; written by HPOINT and others); IOUNITS_MOD (the output unit).

The method

Mixing. For a convection zone, the starting abundances are the mass-weighted means $\bar x_j = \sum_i x_{j,i}\,\Delta m_i / \sum_i \Delta m_i$ over the zone (elements 1–11), and so are the rates. If the zone's (or shell's) first point is below TCUT(1), it only homogenizes and returns.

Question 1 (lines 196–222 ↗). Mass-averaged rates in convective cores. In a convection zone the rates themselves are mass-averaged and the zone burns as one, with the mean abundances. Averaging the rates rather than the products rate × abundances is exact only for instantaneous mixing. Is that adequate for ³He and the CN isotopes in convective cores, where the burning time of ¹²C can be shorter than the convective turnover time, and is a diffusive treatment of convective mixing planned? 2. Equations. Only the seven abundant species are followed, as the YREC paper's "full non-equilibrium burning" (§2.3.4, with ¹⁸O, step 4). The coefficients below show that a rate can stand for more than one reaction: $-3r_1X^2$ takes three protons per pp reaction (p+p→d and d+p→³He), and $-2X\,r_6\,{}^{14}\mathrm{N}$ two per ¹⁴N+p; ¹⁵N and ¹⁷O, "not currently used" (the comment), do not appear. Which chains each rate stands for, and why the nuclei between are skipped, is the question at the end of this step. Primordial deuterium and the light elements are burned separately, by DBURN and LIBURN. The rates HR$_j$ are, in the comment's words, "per gigayear per amu," to be "multiplied by the mass fractions of the reactants and the atomic weight (in amu) of the product to get the rate of change of the mass fraction," which gives the form below (how $N_A$ and the $1/A$ factors enter is rates.f, question 3). With $\Delta t$ in Gyr and $A_k$ the species' mass numbers (1, 3, 4, 12, 13, 14, 16), the residuals are $$G_k(w) = w_k - w_k^0 - A_k\,\Delta t\,R_k(w),$$ where $w = (X, {}^3\mathrm{He}, Y, {}^{12}\mathrm{C}, {}^{13}\mathrm{C}, {}^{14}\mathrm{N}, {}^{16}\mathrm{O})$ are the end-of-step abundances (so the scheme is implicit) and, with $r_n$ the rate GRn and $f_3$ = HF1, $f_4 = 1-f_3$ (the branching of ¹⁵N + p to ¹²C + α or to ¹⁶O): $$R_X = -3r_1X^2 + 2r_2\,{}^3\mathrm{He}^2 - r_3\,{}^3\mathrm{He}\,Y - r_4X\,{}^{12}\mathrm{C} - r_5X\,{}^{13}\mathrm{C} - 2X\,(r_6\,{}^{14}\mathrm{N} + r_7\,{}^{16}\mathrm{O}),$$ $$R_{^3\mathrm{He}} = r_1X^2 - 2r_2\,{}^3\mathrm{He}^2 - r_3\,{}^3\mathrm{He}\,Y,$$ $$R_Y = r_2\,{}^3\mathrm{He}^2 + r_3\,{}^3\mathrm{He}\,Y + f_3r_6X\,{}^{14}\mathrm{N} + r_7X\,{}^{16}\mathrm{O} - r_{10}Y\,{}^{12}\mathrm{C} - r_8Y\,{}^{13}\mathrm{C} - r_{11}Y\,{}^{14}\mathrm{N} - 3r_{12}Y^3,$$ $$R_{^{12}\mathrm{C}} = r_{12}Y^3 - r_4X\,{}^{12}\mathrm{C} - r_{10}Y\,{}^{12}\mathrm{C} + f_3r_6X\,{}^{14}\mathrm{N},\qquad R_{^{13}\mathrm{C}} = r_4X\,{}^{12}\mathrm{C} - r_5X\,{}^{13}\mathrm{C} - r_8Y\,{}^{13}\mathrm{C},$$ $$R_{^{14}\mathrm{N}} = r_5X\,{}^{13}\mathrm{C} + r_7X\,{}^{16}\mathrm{O} - r_6X\,{}^{14}\mathrm{N} - r_{11}Y\,{}^{14}\mathrm{N},\qquad R_{^{16}\mathrm{O}} = f_4r_6X\,{}^{14}\mathrm{N} + r_{10}Y\,{}^{12}\mathrm{C} - r_7X\,{}^{16}\mathrm{O} + r_8Y\,{}^{13}\mathrm{C}.$$ The reactions by the comments: $r_1$ pp, $r_2$ ³He+³He, $r_3$ ³He+⁴He, $r_4$ ¹²C+p, $r_5$ ¹³C+p, $r_6$ ¹⁴N+p (with ¹⁵N+p), $r_7$ ¹⁶O+p (with ¹⁷O+p), $r_8$ ¹³C+α, $r_{10}$ ¹²C+α, $r_{11}$ ¹⁴N+α, $r_{12}$ triple α; $r_9$ (¹⁶O+α) and $r_{13}$ (¹²C+¹²C) are zeroed.

Question 2 (lines 264–284 ↗). Which chains the rates stand for. The code follows seven species and its coefficients count three protons per pp reaction. That reads as if short-lived nuclei (¹³N, ¹⁵O, ¹⁷F, which β-decay within minutes) and those burned fast at these temperatures (²H made by the pp reaction, ⁷Be, ⁷Li, ¹⁵N, ¹⁷O) were taken to be processed at once, so that $r_1$ stands for p+p→d and d+p→³He, $r_3$ for ³He+⁴He through ⁷Be and ⁷Li to two ⁴He, $r_6$ for ¹⁴N+p through ¹⁵O and ¹⁵N+p, and $r_7$ for ¹⁶O+p through ¹⁷F and ¹⁷O+p, which would be why the CN cycle here has no ¹⁵N. Is that the design, and where does it stop holding (for ¹⁷O, say, at low temperature)? 3. Newton iteration. The Jacobian $\partial G_k/\partial w_l$ fills A(1..49) column by column and the residuals $G_k$ fill A(50..56) (the array B of rp_1, by EQUIVALENCE, now addressed as A(49+I)); entries below $10^{-8}$ off the diagonal are zeroed, and SIMEQC (which solves $J\,\delta = G$, leaving $\delta$ in A(50..56)) is skipped when the matrix is diagonal. Then $w \leftarrow w - \delta$; abundances below CMIN are set to 0. It stops when $\max|\delta_k| <$ ABSTOL and $\max|\delta_k/w_k| <$ RELTOL, or stops the run after KEMMAX iterations (STOP 1, the run-time error code). The first guess is the result for the shell below (or the shell's current HCOMP for the first shell).

Question 3 (lines 349–360 ↗). The diagonal shortcut. When all off-diagonal entries are below $10^{-8}$, SIMEQC is skipped and the correction is the residual $G_k$ itself, not $G_k/J_{kk}$. The diagonal entries are $1 + A_k\Delta t\,(\dots)$, close to 1 only for small rates; is the undivided step intended?

The comment on the first guess is out of date. It still speaks of calls from CRRECT with ITLVL > 1 ("USE THE PREVIOUS RESULT," lines 234–237 ↗), an argument now removed; the code takes the result for the shell below, or the shell's own HCOMP for shell 1, with no ITLVL branch (lines 238–254 ↗), as the comment's next lines say (settled from the code, Nov 10, 2026).

Question 4 (line 354 ↗; line 395 ↗). Hard-coded thresholds. $10^{-8}$ for zeroing matrix entries and $5\times10^{-7}$ for resetting $Z$. Are these still appropriate, or should they be settings beside CMIN, ABSTOL and RELTOL? 4. After convergence. ¹⁸O explicitly: $x_{18} = x_{18}^0 + 18\,\Delta t\,r_{11}\,Y\,{}^{14}\mathrm{N}$. If the previous model's $X$ at the zone's top (HCOMPP(1,IEND)) is below $5\times10^{-7}$, $Z$ is reset to $1 - X - {}^3\mathrm{He} - Y$ from the new values; otherwise $Z$ is the zone's mean at the start of the step. All shells of the zone get the same new values. Mixing is thus taken to be instantaneous: the YREC paper (§2.3.5) calls this usually an excellent approximation, which breaks down in such cases as AGB shell flashes.

Question 5 (line 393 ↗). ¹⁸O is only produced. The explicit update adds ¹⁴N(α,γ) production but no destruction. Intended?

Question 6 (lines 395–399 ↗). The switch for $Z$. Until the previous model's $X$ at the zone's top falls below $5\times10^{-7}$, $Z$ is held at its starting mean, although CNO processing and ¹²C(α,γ) change its parts; after, $Z = 1 - X - {}^3\mathrm{He} - Y$, so the C and O made by helium burning count as metals (as the opacity's $Z > 0.1$ branch in GETOPAC assumes). Why the switch at that threshold, and does $X + {}^3\mathrm{He} + Y + Z$ stay exactly 1 before it?

In a burning zone, elements 8 and 10 (¹⁵N and ¹⁷O) are left as they were (suspected bug 2).

From YREC's documentation

Physics Options (Nuclear reaction & burning options) gives the solver's settings: CMIN (1.0D-20), the smallest abundance followed; ABSTOL and RELTOL (1.0D-5, 1.0D-4), the absolute and relative convergence criteria for the changes in abundance; KEMMAX (50), the largest number of substeps; and TCUT(1) (6.5), the $\log T$ above which $\epsilon_\mathrm{nuc}$ is computed, the other TCUT entries being thresholds for later burning stages. Its advice: only TCUT is adjusted from the defaults (6.5 is too high for deuterium burning on the pre-main sequence), and the others should not be changed without good reason.

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) The error message gives a fixed limit of 50.

(lines 377–385 ↗) When the iteration has not converged after KEMMAX iterations, the message says "RUN STOPPED AFTER 50 ATTEMPTS," while the limit is KEMMAX (50 by default, a namelist setting). With another KEMMAX, the message gives the wrong count.

A possible fix:

kemcom.f

@@ lines 377–386 @@
          IF(ICON.GE.KEMMAX) THEN
 C MHP 10/02 IU not defined
 C            WRITE (ISHORT,1000) IU
-            WRITE (ISHORT,1000) IBEGIN
+C 11/09/26 KC via Claude
+C code changed to print the iteration limit KEMMAX, not a fixed 50.
+            WRITE (ISHORT,1000) IBEGIN,KEMMAX
  1000       FORMAT(1X,39('>'),40('<')/1X,'ERROR IN SUBROUTINE KEMCOM'/
      *      1X,'UNABLE TO SOLVE FOR NEW ABUNDANCES IN SHELL',I4/1X,
-     *      'RUN STOPPED AFTER 50 ATTEMPTS')
+     *      'RUN STOPPED AFTER',I4,' ATTEMPTS')
 C 10/30/26 MHP via Claude: bare STOP -> STOP 1, the run-time error code (issue 38).
             STOP 1
          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 in any of 34 runs (26 of YREC's cases, from a 0.04 $M_\odot$ brown dwarf and the solar models to 9 $M_\odot$, rotation, and helium burning, and 8 made-up ones): The network always converged within KEMMAX = 50 iterations, the value in every case, so the message was never printed. It is wrong only with another KEMMAX, a namelist setting that no case changes, and then only when the network fails.

2. (never occurred) ¹⁵N and ¹⁷O are not mixed in burning zones.

(lines 400–410 ↗; lines 121–127 ↗) When a convection zone burns, loop 40 gives every shell of the zone the new values of nine species but leaves elements 8 and 10 of HCOMP (¹⁵N and ¹⁷O) as they were, whereas the non-burning branch homogenizes all 11 (lines 121–127 ↗). Within the step nothing repairs it: CRRECT resets HCOMP to HCOMPP before calling MIX at ITLVL > 2 (crrect.f, lines 150–158 ↗), and MIX's own homogenization (mix.f, lines 360–397 ↗) runs only with LSEMIC or at ITLVL = 1; after convergence, MIXCZ mixes all 11 species again, over the shells that LC flags convective (main.f, line 977 ↗; mixcz.f, lines 63–110 ↗). What follows: ¹⁵N and ¹⁷O can stay unmixed in a burning convection zone. They enter only the Debye–Hückel composition ZDH (with LDH, off by default: crrect.f, lines 141–142 ↗; coefft.f, lines 212–213 ↗) and the output.

A possible fix:

kemcom.f

@@ lines 407–412 @@
          HCOMP(9,IU) = W(7)
          HCOMP(3,IU) = Z
          HCOMP(11,IU) = XO18
+C 11/10/26 KC via Claude
+C code changed to give N15 and O17 their zone means, as the
+C non-burning branch does.
+         HCOMP(8,IU) = XSUM(8)
+         HCOMP(10,IU) = XSUM(10)
    40 CONTINUE
   200 RETURN
       END

Checked by running YREC (Nov 10, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): A burning convection zone reached loop 40 11,559 times in the five cases run: 807 times in testsuite/Test_m1p0feh+0p0GS98_base_TAMS, 266 in testsuite/Test_m0p3feh+0p0GS98_base_ZAMS (which sets LDH, so ZDH is read), 1,489 in testsuite/Test_m1p4feh+0p0GS98_dif_TAMS (with diffusion), 4,567 in testsuite/Test_m0p8feh+0p0GS98_rotp8p0d0p1mag1e4std_TAMS (with rotation), and 4,430 in examples/run_standard_solar_model/Test_solarcal_noGS_norot. Not once did a shell's ¹⁵N or ¹⁷O differ from the zone's mean by more than a relative $10^{-10}$: The zones were already uniform in both, so leaving them changed nothing. The fix was not run, since it would move these values by less than that. What would reach it: a burning zone spreading over layers of different ¹⁵N or ¹⁷O, which KEMCOM's network, not following them, cannot make. Intent: The comment calls the two "not currently used" (line 79 ↗).

Definitions

subroutineKEMCOM

lines 4–412 · GitHub

Implicit Newton solution of the 7-species network over one step (equations above), with explicit ¹⁸O.

Uses modules DIMENSIONS_MOD, IOUNITS_MOD, NMLPAR_MOD, OLDMOD_MOD

Called by MIX ×2, ROTMIX ×2

Flowchart of its calls, in order (1)

if (J1J.NE.0): SIMEQC 360

Change log

The changes since rp_1

Changed since rp_1. The EQUIVALENCE of B with A(50..56) removed, B(I) written as A(49+I) (the same storage, so the same results); COMMON blocks to modules, IMPLICIT NONE. No question settled.

Changed since a8a5950. The stop after KEMMAX iterations now exits with status 1 (STOP 1), the run-time error code, instead of a bare STOP, which exited with 0 (KNOWN_ISSUES.md, issue 38).

1C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 2C KEMCOM 3C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 4 SUBROUTINE KEMCOM(HT,IBEGIN,IEND,HR1,HR2,HR3,HR4,HR5,HR6,HR7,HR8, 5C *HR9,HR10,HR11,HR12,HR13,HF1,HS2,HCOMP,DDAGE,ITLVL) ! KC 2025-05-31 6 * *HR10,HR11,HR12,HF1,HS2,HCOMP,DDAGE) 7C 10/15/26 MHP via Claude 8C code changed to take JSON from DIMENSIONS_MOD instead of declaring 9C PARAMETER(JSON=5000) here. 10 USE DIMENSIONS_MOD 11C 10/14/26 MHP via Claude 12C code changed to convert COMMON/OLDMOD/ into OLDMOD_MOD. 13 USE OLDMOD_MOD 14C 10/14/26 MHP via Claude 15C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD. 16 USE NMLPAR_MOD 17C 10/14/26 MHP via Claude 18C code changed to convert COMMON/IOUNITS/ into IOUNITS_MOD. 19 USE IOUNITS_MOD 20C 10/20/26 MHP via Claude 21C code changed to replace implicit typing with IMPLICIT NONE and 22C explicit declarations. 23 IMPLICIT NONE 24 REAL*8 A,BIGB,BIGD,DCMIN,DDAGE,DELT,DELT12,DELT13,DELT14,DELT16, 25 * * DELT3,DELT4,F3,F4,GR1,GR10,GR11,GR12,GR13,GR2,GR3,GR4,GR5, 26 * * GR6,GR7,GR8,GR9,HCOMP,HF1,HR1,HR10,HR11,HR12,HR2,HR3,HR4, 27 * * HR5,HR6,HR7,HR8,HS2,HT,W,WTSUM,X2R,X3R,X4R,X6R,XB,XO18,XPR, 28 * * XSUM,Y3R,YPR,Z 29 INTEGER I,IBEGIN,ICON,IEND,IU,J,J1J,J2J 30 31C 10/12/26 MHP via Claude 32C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 33C namelist variables that the physics routines read from 27 COMMON 34C blocks which also held derived or working members. Those blocks keep 35C only their non-namelist members; five of them -- NEWMX, VMULT, 36C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 37C Derivation-only namelist variables from the same blocks went to 38C COMMON/NMLBYPASS/ instead, not here. 39C 40C Declared identically in all 121 files that hold any of these members. 41C Generated from one member list, so every copy is byte-identical; the 42C linker's 'size of symbol changed' warning is the check that they are. 43C 10/13/26 MHP via Claude 44C code changed to drop COMMON/CTLIM/, declared here but never used. 45C 9/25/26 MHP via Claude 46C code changed to consolidate all I/O logical unit numbers, previously 47C scattered across COMMON/LUOUT/, COMMON/LUNUM/, COMMON/ALEX06/, 48C COMMON/ACDPTH/, COMMON/NEWOPAC/, COMMON/LOPAL95/, COMMON/ATMOS2/, 49C COMMON/ALATM03/, COMMON/OPALEOS/, and COMMON/SCV2/, into a single 50C COMMON/IOUNITS/ 51C 10/13/26 MHP via Claude 52C code changed to remove EQUIVALENCE (A(50),B(1)), the only EQUIVALENCE 53C left in the code. A holds the 7x7 Jacobian in A(1:49) and the 54C right-hand side in A(50:56); SIMEQC, which declares DIMENSION A(56), 55C is handed the whole array and solves in place, leaving the solution in 56C A(50:56). B existed only as an alias for reading that solution back, 57C in three lines of the convergence check below, so those now index A 58C directly as A(49+I) -- identical by the definition of the EQUIVALENCE, 59C since the loop runs I = 1,7. A itself is untouched, so the contiguity 60C SIMEQC relies on is preserved. 61 DIMENSION A(56),W(7),XSUM(11) 62 DIMENSION HT(JSON),HR1(JSON),HR2(JSON),HR3(JSON),HR4(JSON), 63C *HR5(JSON),HR6(JSON),HR7(JSON),HR8(JSON),HR9(JSON),HR10(JSON), 64C *HR11(JSON),HR12(JSON),HR13(JSON),HF1(JSON),HCOMP(15,JSON), ! KC 2025-05-31 65 * *HR5(JSON),HR6(JSON),HR7(JSON),HR8(JSON),HR10(JSON), 66 * *HR11(JSON),HR12(JSON),HF1(JSON),HCOMP(15,JSON), 67 * * HS2(JSON) 68 SAVE 69 70C KEMCOM USES THE REACTION RATES COMPUTED IN ENGEB TO IMPLICITLY SOLVE 71C FOR THE ABUNDANCES OF H,HE3,HE4,C12,C13,N14,O16, AND O18 SIMULTANEUOSLY. 72C DON VANDENBERGS NOTES EXPLAIN WHAT THE SR IS DOING. 73C 74C INPUT VARIABLES : 75C 76C DDAGE - MODEL TIME STEP IN YEARS. 77C HCOMP - ARRAY CONTAINING RUN OF MASS FRACTIONS OF DIFFERENT SPECIES. 78C HCOMP(1,..)=H; 2=HE4; 3=Z; 4=HE3; 5=C12; 6=C13; 7=N14; 9=O16; 11=O18. 79C ELEMENTS 8(N15) AND 10(O17) NOT CURRENTLY USED;12-15 ARE LIGHT ELEMENTS 80C WHOSE BURNING IS TREATED ELSEWHERE. 81C HCOMPP - THE ARRAY OF ABUNDANCES AT THE START OF THE TIMESTEP. 82C HS2 - RUN OF MASS CONTAINED IN EACH SHELL. 83C HT - RUN OF MODEL TEMPERATURE. 84C IBEGIN, IEND - STARTING AND ENDING SHELLS; DIFFERENT FOR A CONVECTION ZONE. 85C HR1-HR13,HF1 - REACTION RATES (EXCLUDING TERMS THAT DEPEND ONLY ON 86C THE COMPOSITION) PER GIGAYEAR PER AMU. NEEDS TO BE MULTIPLIED BY THE 87C MASS FRACTIONS OF THE REACTANTS AND THE ATOMIC WEIGHT (IN AMU) OF THE 88C PRODUCT TO GET THE RATE OF CHANGE OF THE MASS FRACTION. 89C 90C NUMERICAL PARAMETERS IN COMMON BLOCK /BURTOL/: 91C ABUNDANCES BELOW CMIN ARE ZEROED OUT. 92C ABSTOL AND RELTOL ARE THE ABSOLUTE AND RELATIVE TOLERANCES FOR CONVERGENCE. 93C KEMMAX IS THE MAXIMUM NUMBER OF ITERATIONS BEFORE THE PROGRAM WILL HALT. 94C 95C OUTPUT VARIABLES : THE NEW RUN OF ABUNDANCE HCOMP. 96C 97 IF(IBEGIN.NE.IEND)THEN 98C HOMOGENIZE CONVECTION ZONES. 99C XSUM IS THE MASS-WEIGHTED AVERAGE ABUNDANCE FOR THE CZ. 100C WTSUM IS THE TOTAL MASS OF THE CZ. 101C INITIALIZE SUMS. 102 WTSUM = 0.0D0 103 DO 1 J = 1,11 104 XSUM(J) = 0.0D0 105 1 1 CONTINUE 106 DO 5 I = IBEGIN,IEND 107 WTSUM = WTSUM + HS2(I) 108 DO 3 J = 1,11 109 XSUM(J) = XSUM(J)+HCOMPP(J,I)*HS2(I) 110 3 3 CONTINUE 111 5 5 CONTINUE 112 DO 7 J = 1,11 113 XSUM(J) = XSUM(J)/WTSUM 114 7 7 CONTINUE 115 ELSE 116 DO 9 J = 1,11 117 XSUM(J) = HCOMPP(J,IBEGIN) 118 9 9 CONTINUE 119 ENDIF 120C SKIP BURNING CALCULATIONS IF STARTING SHELL BELOW T CUTOFF FOR REACTIONS. 121 IF(HT(IBEGIN).LT.TCUT(1))THEN 122 DO 13 I = IBEGIN,IEND 123 DO 11 J = 1,11 124 HCOMP(J,I) = XSUM(J) 125 11 11 CONTINUE 126 13 13 CONTINUE 127 GOTO 200 128 ENDIF 129C 130C SET UP NUMERICAL PARAMETERS. 131C 132C TIMESTEP IN GIGAYEARS. THE OTHER DELTS ARE THIS TIMESTEP MULTIPLIED 133C BY THE ATOMIC MASS OF THE DIFFERENT SPECIES (IN AMU). 134 DELT=DDAGE*1.0D-9 135 DELT3 = 3.0D0*DELT 136 DELT4 = 4.0D0*DELT 137 DELT12 = 12.0D0*DELT 138 DELT13 = 13.0D0*DELT 139 DELT14 = 14.0D0*DELT 140 DELT16 = 16.0D0*DELT 141 DCMIN=CMIN 142C COUNTER FOR THE NUMBER OF ITERATIONS. 143 ICON=0 144C 145C NUCLEAR REACTION RATES. 146C 147C THESE REACTIONS ARE NO LONGER INCLUDED AND ARE ZEROED OUT. 148 GR9 = 0.0D0 149 GR13 = 0.0D0 150 IF(IBEGIN.EQ.IEND)THEN 151C PP 152 GR1 = HR1(IBEGIN) 153C HE3,HE3 154 GR2 = HR2(IBEGIN) 155C HE3,HE4 156 GR3 = HR3(IBEGIN) 157C C12,P 158 GR4 = HR4(IBEGIN) 159C C13,P 160 GR5 = HR5(IBEGIN) 161C N14,P + N15,P 162 GR6 = HR6(IBEGIN) 163C O16,P + O17,P. 164 GR7 = HR7(IBEGIN) 165C C13,ALPHA 166 GR8 = HR8(IBEGIN) 167C O16,ALPHA (NOT USED) 168C GR9 = HR9(IBEGIN) 169C C12,ALPHA 170 GR10 = HR10(IBEGIN) 171C N14,ALPHA 172 GR11 = HR11(IBEGIN) 173C TRIPLE ALPHA 174 GR12 = HR12(IBEGIN) 175C C12,C12 (NOT USED) 176C GR13 = HR13(IBEGIN) 177C BRANCHING RATIO FOR N15,P : 178C F3 = FRACTION GOING TO C12+ALPHA 179C F4 = FRACTION GOING TO O16 180 F3 = HF1(IBEGIN) 181 F4 = 1.0D0 - F3 182 ELSE 183C USE THE MASS-WEIGHTED AVERAGE RATES FOR THE CZ 184 GR1 = 0.0D0 185 GR2 = 0.0D0 186 GR3 = 0.0D0 187 GR4 = 0.0D0 188 GR5 = 0.0D0 189 GR6 = 0.0D0 190 GR7 = 0.0D0 191 GR8 = 0.0D0 192 GR10 = 0.0D0 193 GR11 = 0.0D0 194 GR12 = 0.0D0 195 F3 = 0.0D0 196 DO 15 I = IBEGIN,IEND 197 GR1 = GR1 + HS2(I)*HR1(I) 198 GR2 = GR2 + HS2(I)*HR2(I) 199 GR3 = GR3 + HS2(I)*HR3(I) 200 GR4 = GR4 + HS2(I)*HR4(I) 201 GR5 = GR5 + HS2(I)*HR5(I) 202 GR6 = GR6 + HS2(I)*HR6(I) 203 GR7 = GR7 + HS2(I)*HR7(I) 204 GR8 = GR8 + HS2(I)*HR8(I) 205 GR10 = GR10 + HS2(I)*HR10(I) 206 GR11 = GR11 + HS2(I)*HR11(I) 207 GR12 = GR12 + HS2(I)*HR12(I) 208 F3 = F3 + HS2(I)*HF1(I) 209 15 15 CONTINUE 210 GR1 = GR1/WTSUM 211 GR2 = GR2/WTSUM 212 GR3 = GR3/WTSUM 213 GR4 = GR4/WTSUM 214 GR5 = GR5/WTSUM 215 GR6 = GR6/WTSUM 216 GR7 = GR7/WTSUM 217 GR8 = GR8/WTSUM 218 GR10 = GR10/WTSUM 219 GR11 = GR11/WTSUM 220 GR12 = GR12/WTSUM 221 F3 = F3/WTSUM 222 F4 = 1.0D0 - F3 223 ENDIF 224C 225C INITIAL ABUNDANCES OF SPECIES. 226C 227 XPR = XSUM(1) 228 Y3R = XSUM(4) 229 YPR = XSUM(2) 230 X2R = XSUM(5) 231 X3R = XSUM(6) 232 X4R = XSUM(7) 233 X6R = XSUM(9) 234C STARTING GUESS FOR FINAL ABUNDANCES OF SPECIES : 235C FOR CALLS IN CRRECT (ITLVL>1) USE THE PREVIOUS RESULT. 236C OTHERWISE USE THE FINAL ABUNDANCES FOUND FOR THE PREVIOUS SHELL 237C UNLESS I = 1. 238 IF(IBEGIN.NE.1) THEN 239 W(1) = HCOMP(1,IBEGIN-1) 240 W(2) = HCOMP(4,IBEGIN-1) 241 W(3) = HCOMP(2,IBEGIN-1) 242 W(4) = HCOMP(5,IBEGIN-1) 243 W(5) = HCOMP(6,IBEGIN-1) 244 W(6) = HCOMP(7,IBEGIN-1) 245 W(7) = HCOMP(9,IBEGIN-1) 246 ELSE 247 W(1) = HCOMP(1,IBEGIN) 248 W(2) = HCOMP(4,IBEGIN) 249 W(3) = HCOMP(2,IBEGIN) 250 W(4) = HCOMP(5,IBEGIN) 251 W(5) = HCOMP(6,IBEGIN) 252 W(6) = HCOMP(7,IBEGIN) 253 W(7) = HCOMP(9,IBEGIN) 254 ENDIF 255C 256C ITERATION FOR NEW ABUNDANCES. 257C 258 10 10 CONTINUE 259C FUNCTIONS TO BE MINIMIZED. THESE ARE OF THE FORM 260C A(#)=ABUNDANCE(END STEP)-ABUNDANCE (START STEP)-TIMESTEP*(D(SPECIES)/DT) 261C BECAUSE THE BURNING RATES D(SPECIES)/DT USE THE ABUNDANCES AT THE 262C END OF THE TIMESTEP W(1)-W(7), THIS IS AN IMPLICIT SCHEME. 263C EQUATION FOR H (X) 264 A(50) = W(1)-XPR-DELT*(-3.D0*GR1*W(1)**2+2.D0*GR2*W(2)**2-GR3 265 * * *W(2)*W(3)-GR4*W(1)*W(4)-GR5*W(1)*W(5)-2.D0*W(1)*(GR6*W(6)+ 266 * * GR7*W(7))) 267C HE3 268 A(51)= W(2)-Y3R-DELT3*(GR1*W(1)**2-2.D0*GR2*W(2)**2-GR3*W(2) 269 * * *W(3)) 270C HE4 (Y) 271 A(52) = W(3)-YPR-DELT4*(GR2*W(2)**2+GR3*W(2)*W(3)+F3*GR6*W(6) 272 * * *W(1)+GR7*W(1)*W(7)-GR10*W(3)*W(4)-GR8*W(3)*W(5)-GR11*W(3) 273 * * *W(6)-3.D0*GR12*W(3)**3) 274C C12 275 A(53) = W(4)-X2R-DELT12*(GR12*W(3)**3-GR4*W(1)*W(4)-GR10*W(3) 276 * * *W(4)+F3*GR6*W(1)*W(6)) 277C C13 278 A(54)=W(5)-X3R-DELT13*(GR4*W(1)*W(4)-GR5*W(1)*W(5)-GR8*W(3)*W(5)) 279C N14 280 A(55) = W(6)-X4R-DELT14*(GR5*W(1)*W(5)+GR7*W(1)*W(7)-GR6*W(1)* 281 * * W(6) - GR11*W(3)*W(6)) 282C O16*W(4) 283 A(56) = W(7)-X6R-DELT16*(F4*GR6*W(1)*W(6)+GR10*W(3)*W(4)-GR7*W(1) 284 * * *W(7)+GR8*W(3)*W(5)) 285C NOTE THAT THE O18 ABUNDANCE IS COMPUTED EXPLICITLY AT THE END OF THE 286C SUBROUTINE. 287C THE REMAINING VALUES OF A CORRESPOND TO THE DERIVATIVES OF THE ABOVE 7 288C FUNCTIONS WITH RESPECT TO THE 7 SPECIES. 289C D(RATES)/DX 290 A(1) = 1.D0+DELT*(6.D0*GR1*W(1)+GR4*W(4)+GR5*W(5)+2.D0*(GR6*W(6)+ 291 * * GR7*W(7))) 292 A(2) = -2.D0*DELT3*GR1*W(1) 293 A(3) = -DELT4*(F3*GR6*W(6)+GR7*W(7)) 294 A(4) = -DELT12*(F3*GR6*W(6)-GR4*W(4)) 295 A(5) = -DELT13*(GR4*W(4)-GR5*W(5)) 296 A(6) = -DELT14*(GR5*W(5)+GR7*W(7)-GR6*W(6)) 297 A(7) = -DELT16*(F4*GR6*W(6)-GR7*W(7)) 298C D(RATES)/DHE3 299 A(8) = -DELT*(4.D0*GR2*W(2)-GR3*W(3)) 300 A(9) = 1.D0+DELT3*(4.D0*GR2*W(2)+GR3*W(3)) 301 A(10) = -DELT4*(2.D0*GR2*W(2)+GR3*W(3)) 302 A(11) = 0.D0 303 A(12) = 0.D0 304 A(13) = 0.D0 305 A(14) = 0.D0 306C D(RATES)/DHE4 307 A(15) = DELT*GR3*W(2) 308 A(16) = DELT3*GR3*W(2) 309 A(17) = 1.D0-DELT4*(GR3*W(2)-GR10*W(4)-GR8*W(5)-GR11*W(6) 310 * * -9.D0*GR12*W(3)**2) 311 A(18) = -DELT12*(3.D0*GR12*W(3)**2-GR10*W(4)) 312 A(19) = DELT13*GR8*W(5) 313 A(20) = DELT14*GR11*W(6) 314 A(21) = -DELT16*(GR10*W(4)+GR8*W(5)) 315C D(RATES)/DC12 316 A(22) = DELT*GR4*W(1) 317 A(23) = 0.D0 318 A(24) = DELT4*GR10*W(3) 319 A(25) = 1.D0+DELT12*(GR4*W(1)+GR10*W(3)) 320 A(26) = -DELT13*GR4*W(1) 321 A(27) = 0.D0 322 A(28) = -DELT16*GR10*W(3) 323C D(RATES)/DC13 324 A(29) = DELT*GR5*W(1) 325 A(30) = 0.D0 326 A(31) = DELT4*GR8*W(3) 327 A(32) = 0.D0 328 A(33) = 1.D0+DELT13*(GR5*W(1)+GR8*W(3)) 329 A(34) = -DELT14*GR5*W(1) 330 A(35) = -DELT16*GR8*W(3) 331C D(RATES)/DN14 332 A(36) = 2.D0*DELT*GR6*W(1) 333 A(37) = 0.D0 334 A(38) = -DELT4*(F3*GR6*W(1)-GR11*W(3)) 335 A(39) = -DELT12*F3*GR6*W(1) 336 A(40) = 0.D0 337 A(41) = 1.D0+DELT14*(GR6*W(1)+GR11*W(3)) 338 A(42) = -F4*DELT16*GR6*W(1) 339C D(RATES)/DO16 340 A(43) = 2.D0*DELT*GR7*W(1) 341 A(44) = 0.D0 342 A(45) = -DELT4*GR7*W(1) 343 A(46) = 0.D0 344 A(47) = 0.D0 345 A(48) = -DELT14*GR7*W(1) 346 A(49) = 1.D0+DELT16*GR7*W(1) 347C REDUCE THE 7X8 SYSTEM TO FIND THE CORRECTIONS TO THE RELATIVE 348C ABUNDANCES. 349 J1J = 0 350 J2J = 1 351 DO 20 I = 1,49 352 IF(I.EQ.J2J) THEN 353 J2J = J2J+8 354 ELSE IF(DABS(A(I)).LT.1.D-8) THEN 355 A(I) = 0.D0 356 ELSE 357 J1J = 1 358 ENDIF 359 20 20 CONTINUE 360 IF(J1J.NE.0) CALL SIMEQC(A,8,7) 361C CHECK TO SEE IF THE SYSTEM HAS CONVERGED WITHIN THE DESIRED TOLERANCES. 362 BIGB = 0.D0 363 BIGD = 0.D0 364 DO 30 I = 1,7 365 W(I) = W(I)-A(49+I) 366 IF(DABS(A(49+I)).GE.BIGD) BIGD = DABS(A(49+I)) 367 IF(W(I).LT.DCMIN) THEN 368 W(I) = 0.D0 369 ELSE 370 XB = DABS(A(49+I)/W(I)) 371 IF(XB.GE.BIGB) BIGB = XB 372 ENDIF 373 30 30 CONTINUE 374 IF(BIGD.GE.ABSTOL.OR.BIGB.GE.RELTOL) THEN 375C SYSTEM NOT CONVERGED - SEE IF MAXIMUM NUMBER OF ITERATIONS EXCEEDED. 376 ICON = ICON+1 377 IF(ICON.GE.KEMMAX) THEN 378C MHP 10/02 IU not defined 379C WRITE (ISHORT,1000) IU 380 WRITE (ISHORT,1000) IBEGIN 381 1000 1000 FORMAT(1X,39('>'),40('<')/1X,'ERROR IN SUBROUTINE KEMCOM'/ 382 * * 1X,'UNABLE TO SOLVE FOR NEW ABUNDANCES IN SHELL',I4/1X, 383 * * 'RUN STOPPED AFTER 50 ATTEMPTS') 384C 10/30/26 MHP via Claude: bare STOP -> STOP 1, the run-time error code (issue 38). 385 STOP 1 386 ELSE 387 GOTO 10 388 ENDIF 389 ENDIF 390C SYSTEM HAS CONVERGED. 391C UPDATE COMPOSITION MATRIX. 392C UPDATE O18. 393 XO18 = XSUM(11)+18.0D0*DELT*GR11*W(3)*W(6) 394C CHANGE METAL ABUNDANCE IF X<5.0D-7. 395 IF(HCOMPP(1,IEND).LT.5.0D-7)THEN 396 Z = 1.0D0-W(1)-W(2)-W(3) 397 ELSE 398 Z = XSUM(3) 399 ENDIF 400 DO 40 IU = IBEGIN,IEND 401 HCOMP(1,IU) = W(1) 402 HCOMP(4,IU) = W(2) 403 HCOMP(2,IU) = W(3) 404 HCOMP(5,IU) = W(4) 405 HCOMP(6,IU) = W(5) 406 HCOMP(7,IU) = W(6) 407 HCOMP(9,IU) = W(7) 408 HCOMP(3,IU) = Z 409 HCOMP(11,IU) = XO18 410 40 40 CONTINUE 411 200 200 RETURN 412 END 413

1"""Python version of kemcom.f (YREC development, 7be4133): the implicit 2solution of the CNO and pp burning network for one shell or one convection 3zone.""" 4 5import numpy as np 6 7from pyrec import iounits_mod 8from pyrec import nmlpar_mod 9from pyrec import oldmod_mod 10from pyrec.simeqc import simeqc 11 12# KEMCOM has SAVE: its work arrays are module-level. A(1..49) holds the 13# Jacobian column by column, A(50..56) the residuals (once a separate 14# array B, EQUIVALENCEd to it), which SIMEQC replaces by the Newton correction. 15a = np.zeros(57) 16w = np.zeros(8) 17xsum = np.zeros(12) 18 19 20def kemcom(ht: np.ndarray, ibegin: int, iend: int, hr1: np.ndarray, hr2: np.ndarray, hr3: np.ndarray, hr4: np.ndarray, 21 hr5: np.ndarray, hr6: np.ndarray, hr7: np.ndarray, hr8: np.ndarray, hr10: np.ndarray, hr11: np.ndarray, 22 hr12: np.ndarray, hf1: np.ndarray, hs2: np.ndarray, hcomp: np.ndarray, ddage: float) -> None: 23 """New abundances of H, 3He, 4He, 12C, 13C, 14N and 16O after one step. 24 25 Solves the network implicitly by Newton iteration (VandenBerg's 26 notes), for one radiative shell (IBEGIN = IEND) or one convection zone 27 (mass-weighted means of abundances and rates); 18O is updated 28 explicitly afterwards. 29 30 Translates: kemcom.f#KEMCOM (7be4133) 31 32 Parameters 33 ---------- 34 ht : array 35 log T. 36 ibegin, iend : int 37 First and last shell of the zone. 38 hr1 ... hr12, hf1 : arrays 39 Reaction rates (per Gyr per amu, without composition factors) and 40 the 15N + p branching ratio. 41 hs2 : array 42 Shell masses. 43 hcomp : array 44 Abundances; elements 1-7, 9 and 11 set for shells IBEGIN..IEND. 45 ddage : float 46 Time step (years). 47 48 """ 49 50 nml = nmlpar_mod 51 hcompp = oldmod_mod.hcompp 52 53 # Abundances at the start of the step, mixed over a zone (Fortran lines 54 # 97–119). 55 if ibegin != iend: 56 wtsum = 0.0 57 for j in range(1, 12): # label 1 58 xsum[j] = 0.0 59 for i in range(ibegin, iend + 1): # labels 5, 3 60 wtsum = wtsum + hs2[i] 61 for j in range(1, 12): 62 xsum[j] = xsum[j] + hcompp[j, i] * hs2[i] 63 for j in range(1, 12): # label 7 64 xsum[j] = xsum[j] / wtsum 65 else: 66 for j in range(1, 12): # label 9 67 xsum[j] = hcompp[j, ibegin] 68 69 # Too cool to burn: homogenize only (Fortran lines 121–128). 70 if ht[ibegin] < nml.tcut[1]: 71 for i in range(ibegin, iend + 1): # labels 13, 11 72 for j in range(1, 12): 73 hcomp[j, i] = xsum[j] 74 return # GOTO 200 75 76 # The step in Gyr, times the species' masses (Fortran lines 134–143). 77 delt = ddage * 1.0e-9 78 delt3 = 3.0 * delt 79 delt4 = 4.0 * delt 80 delt12 = 12.0 * delt 81 delt13 = 13.0 * delt 82 delt14 = 14.0 * delt 83 delt16 = 16.0 * delt 84 dcmin = nml.cmin 85 icon = 0 86 87 # Fortran lines 148–149. 88 gr9 = 0.0 # noqa: F841 (unused in the Fortran too) 89 gr13 = 0.0 # noqa: F841 (unused in the Fortran too) 90 91 # The rates: the shell's, or the zone's mass-weighted means, summed in 92 # the Fortran's order (Fortran lines 150–223). 93 if ibegin == iend: 94 k = ibegin 95 gr1, gr2, gr3, gr4 = hr1[k], hr2[k], hr3[k], hr4[k] 96 gr5, gr6, gr7, gr8 = hr5[k], hr6[k], hr7[k], hr8[k] 97 gr10, gr11, gr12 = hr10[k], hr11[k], hr12[k] 98 f3 = hf1[ibegin] 99 f4 = 1.0 - f3 100 else: 101 gr1 = gr2 = gr3 = gr4 = gr5 = gr6 = gr7 = gr8 = 0.0 102 gr10 = gr11 = gr12 = 0.0 103 f3 = 0.0 104 for i in range(ibegin, iend + 1): # label 15 105 gr1 = gr1 + hs2[i] * hr1[i] 106 gr2 = gr2 + hs2[i] * hr2[i] 107 gr3 = gr3 + hs2[i] * hr3[i] 108 gr4 = gr4 + hs2[i] * hr4[i] 109 gr5 = gr5 + hs2[i] * hr5[i] 110 gr6 = gr6 + hs2[i] * hr6[i] 111 gr7 = gr7 + hs2[i] * hr7[i] 112 gr8 = gr8 + hs2[i] * hr8[i] 113 gr10 = gr10 + hs2[i] * hr10[i] 114 gr11 = gr11 + hs2[i] * hr11[i] 115 gr12 = gr12 + hs2[i] * hr12[i] 116 f3 = f3 + hs2[i] * hf1[i] 117 gr1 = gr1 / wtsum 118 gr2 = gr2 / wtsum 119 gr3 = gr3 / wtsum 120 gr4 = gr4 / wtsum 121 gr5 = gr5 / wtsum 122 gr6 = gr6 / wtsum 123 gr7 = gr7 / wtsum 124 gr8 = gr8 / wtsum 125 gr10 = gr10 / wtsum 126 gr11 = gr11 / wtsum 127 gr12 = gr12 / wtsum 128 f3 = f3 / wtsum 129 f4 = 1.0 - f3 130 131 # Starting abundances (note the order of w: X, 3He, Y, 12C, 13C, 14N, 132 # 16O = HCOMP elements 1, 4, 2, 5, 6, 7, 9) (Fortran lines 227–233). 133 xpr = xsum[1] 134 y3r = xsum[4] 135 ypr = xsum[2] 136 x2r = xsum[5] 137 x3r = xsum[6] 138 x4r = xsum[7] 139 x6r = xsum[9] 140 141 # First guess: the shell below, or the shell itself for the first one 142 # (Fortran lines 238–254). 143 # Question (see the page): the comment's ITLVL no longer exists. 144 k = ibegin - 1 if ibegin != 1 else ibegin 145 w[1] = hcomp[1, k] 146 w[2] = hcomp[4, k] 147 w[3] = hcomp[2, k] 148 w[4] = hcomp[5, k] 149 w[5] = hcomp[6, k] 150 w[6] = hcomp[7, k] 151 w[7] = hcomp[9, k] 152 153 while True: # label 10 154 # Residuals G_k (Fortran lines 264–284). 155 a[50] = w[1] - xpr - delt * (-3.0 * gr1 * w[1] ** 2 + 2.0 * gr2 * w[2] ** 2 156 - gr3 * w[2] * w[3] - gr4 * w[1] * w[4] 157 - gr5 * w[1] * w[5] 158 - 2.0 * w[1] * (gr6 * w[6] + gr7 * w[7])) 159 a[51] = w[2] - y3r - delt3 * (gr1 * w[1] ** 2 - 2.0 * gr2 * w[2] ** 2 160 - gr3 * w[2] * w[3]) 161 a[52] = w[3] - ypr - delt4 * (gr2 * w[2] ** 2 + gr3 * w[2] * w[3] 162 + f3 * gr6 * w[6] * w[1] + gr7 * w[1] * w[7] 163 - gr10 * w[3] * w[4] - gr8 * w[3] * w[5] 164 - gr11 * w[3] * w[6] 165 - 3.0 * gr12 * w[3] ** 3) 166 a[53] = w[4] - x2r - delt12 * (gr12 * w[3] ** 3 - gr4 * w[1] * w[4] 167 - gr10 * w[3] * w[4] 168 + f3 * gr6 * w[1] * w[6]) 169 a[54] = w[5] - x3r - delt13 * (gr4 * w[1] * w[4] - gr5 * w[1] * w[5] 170 - gr8 * w[3] * w[5]) 171 a[55] = w[6] - x4r - delt14 * (gr5 * w[1] * w[5] + gr7 * w[1] * w[7] 172 - gr6 * w[1] * w[6] - gr11 * w[3] * w[6]) 173 a[56] = w[7] - x6r - delt16 * (f4 * gr6 * w[1] * w[6] + gr10 * w[3] * w[4] 174 - gr7 * w[1] * w[7] + gr8 * w[3] * w[5]) 175 176 # The Jacobian, column by column (A(1..7) is d G / d X, ...) 177 # (Fortran lines 290–346). 178 a[1] = 1.0 + delt * (6.0 * gr1 * w[1] + gr4 * w[4] + gr5 * w[5] 179 + 2.0 * (gr6 * w[6] + gr7 * w[7])) 180 a[2] = -2.0 * delt3 * gr1 * w[1] 181 a[3] = -delt4 * (f3 * gr6 * w[6] + gr7 * w[7]) 182 a[4] = -delt12 * (f3 * gr6 * w[6] - gr4 * w[4]) 183 a[5] = -delt13 * (gr4 * w[4] - gr5 * w[5]) 184 a[6] = -delt14 * (gr5 * w[5] + gr7 * w[7] - gr6 * w[6]) 185 a[7] = -delt16 * (f4 * gr6 * w[6] - gr7 * w[7]) 186 187 a[8] = -delt * (4.0 * gr2 * w[2] - gr3 * w[3]) 188 a[9] = 1.0 + delt3 * (4.0 * gr2 * w[2] + gr3 * w[3]) 189 a[10] = -delt4 * (2.0 * gr2 * w[2] + gr3 * w[3]) 190 a[11] = 0.0 191 a[12] = 0.0 192 a[13] = 0.0 193 a[14] = 0.0 194 195 a[15] = delt * gr3 * w[2] 196 a[16] = delt3 * gr3 * w[2] 197 a[17] = 1.0 - delt4 * (gr3 * w[2] - gr10 * w[4] - gr8 * w[5] 198 - gr11 * w[6] - 9.0 * gr12 * w[3] ** 2) 199 a[18] = -delt12 * (3.0 * gr12 * w[3] ** 2 - gr10 * w[4]) 200 a[19] = delt13 * gr8 * w[5] 201 a[20] = delt14 * gr11 * w[6] 202 a[21] = -delt16 * (gr10 * w[4] + gr8 * w[5]) 203 204 a[22] = delt * gr4 * w[1] 205 a[23] = 0.0 206 a[24] = delt4 * gr10 * w[3] 207 a[25] = 1.0 + delt12 * (gr4 * w[1] + gr10 * w[3]) 208 a[26] = -delt13 * gr4 * w[1] 209 a[27] = 0.0 210 a[28] = -delt16 * gr10 * w[3] 211 212 a[29] = delt * gr5 * w[1] 213 a[30] = 0.0 214 a[31] = delt4 * gr8 * w[3] 215 a[32] = 0.0 216 a[33] = 1.0 + delt13 * (gr5 * w[1] + gr8 * w[3]) 217 a[34] = -delt14 * gr5 * w[1] 218 a[35] = -delt16 * gr8 * w[3] 219 220 a[36] = 2.0 * delt * gr6 * w[1] 221 a[37] = 0.0 222 a[38] = -delt4 * (f3 * gr6 * w[1] - gr11 * w[3]) 223 a[39] = -delt12 * f3 * gr6 * w[1] 224 a[40] = 0.0 225 a[41] = 1.0 + delt14 * (gr6 * w[1] + gr11 * w[3]) 226 a[42] = -f4 * delt16 * gr6 * w[1] 227 228 a[43] = 2.0 * delt * gr7 * w[1] 229 a[44] = 0.0 230 a[45] = -delt4 * gr7 * w[1] 231 a[46] = 0.0 232 a[47] = 0.0 233 a[48] = -delt14 * gr7 * w[1] 234 a[49] = 1.0 + delt16 * gr7 * w[1] 235 236 # Zero small off-diagonal entries (the diagonal is every 8th, from 237 # 1); solve only if any off-diagonal entry is left (Fortran lines 238 # 349–360). 239 j1j = 0 240 j2j = 1 241 for i in range(1, 50): # label 20 242 if i == j2j: 243 j2j = j2j + 8 244 elif abs(a[i]) < 1.0e-8: 245 a[i] = 0.0 246 else: 247 j1j = 1 248 # Question (see the page): without SIMEQC the correction is the 249 # residual itself, not divided by the diagonal. 250 if j1j != 0: 251 simeqc(a, 8, 7) 252 253 # Newton update (Fortran lines 362–373). 254 bigb = 0.0 255 bigd = 0.0 256 for i in range(1, 8): # label 30 257 w[i] = w[i] - a[49 + i] 258 if abs(a[49 + i]) >= bigd: 259 bigd = abs(a[49 + i]) 260 if w[i] < dcmin: 261 w[i] = 0.0 262 else: 263 xb = abs(a[49 + i] / w[i]) 264 if xb >= bigb: 265 bigb = xb 266 267 # The convergence test (Fortran lines 374–389). 268 if bigd >= nml.abstol or bigb >= nml.reltol: 269 icon = icon + 1 270 if icon >= nml.kemmax: 271 # Suspected bug (see the page): the message says 50, the 272 # limit is KEMMAX. 273 iounits_mod.write(iounits_mod.ishort, 274 " " + ">" * 39 + "<" * 40 + "\n" 275 " ERROR IN SUBROUTINE KEMCOM\n" 276 f" UNABLE TO SOLVE FOR NEW ABUNDANCES IN SHELL{ibegin:4d}\n" 277 " RUN STOPPED AFTER 50 ATTEMPTS") # format 1000 278 raise SystemExit(1) # STOP 1, the run-time error code 279 continue # GOTO 10 280 break 281 282 # 18O explicitly (Fortran line 393). Question (see the page): production 283 # only. 284 xo18 = xsum[11] + 18.0 * delt * gr11 * w[3] * w[6] 285 286 # Z: reset, or the zone's mean (Fortran lines 395–399). 287 # Question (see the page): the threshold 5e-7 is hard-coded. 288 if hcompp[1, iend] < 5.0e-7: 289 z = 1.0 - w[1] - w[2] - w[3] 290 else: 291 z = xsum[3] 292 293 # The new abundances for every shell of the zone (Fortran lines 294 # 400–410). Question (see the page): elements 8 and 10 are not mixed. 295 for iu in range(ibegin, iend + 1): # label 40 296 hcomp[1, iu] = w[1] 297 hcomp[4, iu] = w[2] 298 hcomp[2, iu] = w[3] 299 hcomp[5, iu] = w[4] 300 hcomp[6, iu] = w[5] 301 hcomp[7, iu] = w[6] 302 hcomp[9, iu] = w[7] 303 hcomp[3, iu] = z 304 hcomp[11, iu] = xo18 305 # label 200