deutrate.f

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

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

Variables on this page (7)

cavog, cc13, cln, drate, drate0, jcz, taucz

Purpose

This file computes, at one point, the rate of deuterium burning by $^2$H(p,$\gamma$)$^3$He out of equilibrium ("MHP 5/02"), for the burning done in DBURN and DBURNM. Screening is left out, "reasonable given the low temperatures involved."

How it fits

Inputs: the arguments DL ($\log\rho$), TL ($\log T$), X (hydrogen), I (the point), ITLVL (the level of iteration); JCZ (the base of the surface convection zone) from DEUTER_MOD; TAUCZ (convective overturn time) from OVRTRN_MOD; CLN, CC13 and CAVOG (Avogadro's number, $N_A = 1/m_u$) from CONSTANTS_MOD; the deuteron mass in amu, AMDEUT $= 2.013553$, set here (the same value as ANUC(3) in ENGEB; the comment asks to keep the two in step).
Who calls it: MIX, for every point from the centre out until $T$ falls below TCUT(1), when the surface deuterium (HCOMP(12,M), or the accreted COMPACC(12) if larger, during accretion) exceeds $10^{-11}$.
Outputs: DRATE(I) in DEUTER_MOD, and at the first level of iteration also DRATE0(I) (the rate at the start of the step). DBURN reads both; MIDMOD interpolates between them for DBURNM.

The method

With $T_9 = T/10^9$ K, the reaction rate per mole (the fit of Caughlan and Fowler type, by its form): $$N_A\langle\sigma v\rangle = 2.240\times10^{3}\;T_9^{-2/3}\;e^{-3.72\,T_9^{-1/3}}\left(1 + 0.112\,T_9^{1/3} + 3.38\,T_9^{2/3} + 2.65\,T_9\right),$$ then, with $A_D$ = AMDEUT $= 2.013553$, $$\mathtt{DRATE} = X\,\rho\,N_A\langle\sigma v\rangle \cdot \frac{N_A}{A_D}\cdot C_{21}, \qquad C_{21} = 5.240358\times10^{-8} \approx \frac{3.156\times10^{16}\ \mathrm{s\,Gyr^{-1}}}{N_A},$$ so DRATE $\approx X\rho N_A\langle\sigma v\rangle/A_D$ per Gyr (the value of $C_{21}$ is read as seconds per Gyr over $N_A$; this is an inference). The burning routines then use $X_D \propto e^{-2\,\mathtt{DRATE}\,t}$, so the destruction rate per deuteron is $\lambda = 2\,\mathtt{DRATE} = (2/A_D)\,X\rho N_A\langle\sigma v\rangle \approx 0.993\,X\rho N_A\langle\sigma v\rangle$ per Gyr. The comment calls CAVOG/AMDEUT "Avogadro's number divided by the mass of the deuteron in amu."

The convection-zone cap. In the surface convection zone ($I \ge$ JCZ, with $\tau_\mathrm{cz} > 1$), $X\cdot$RDEUT is capped at $(N_A/A_D)/\tau_\mathrm{cz}$ (CAVOG/AMDEUT/TAUCZ), so that deuterium burns "on a time scale no shorter than the convective overturn timescale." TAUCAL forms TAUCZ as a length in cm over a convective velocity (CZWIDE/CVEAVG or PSCA/CVEL, taucal.f, lines 159–165 ↗, taucal.f, line 179 ↗), so in seconds, with YREC's cgs velocities; then DRATE $= C_{21}\,(N_A/A_D)/\tau_\mathrm{cz}$ is $1/(A_D\,\tau_\mathrm{cz})$ in Gyr$^{-1}$, and the destruction rate per deuteron, $2\,$DRATE, is $0.993/\tau_\mathrm{cz}$: the cap is consistent with the comment (settled from the code, Nov 9, 2026: lines 65–66 ↗, line 81 ↗).

Definitions

subroutineDEUTRATE

lines 5–86 · GitHub

The rate of $^2$H(p,$\gamma$)$^3$He at one point, without screening, in units the burning routines expect (per Gyr), capped in the surface convection zone by the overturn time; sets DRATE (and DRATE0 at the first level).

Uses modules CONSTANTS_MOD, DEUTER_MOD, OVRTRN_MOD

Called by MIX

Change log

The changes since rp_1

Changed since rp_1. C21 is now set in double precision (5.240358D-8), which settles the former question on its single-precision constant; COMMON blocks to modules, IMPLICIT NONE.

Changed since 09c3b9b. The rate's leading power of $T_9$ is now $T_9^{-2/3}$ (T9M23), as in the usual fit, not $T_9^{+2/3}$, which settles former question 1 (the rate rises by $T_9^{-4/3}$, about $10^{4}$ at $10^6$ K); the factor $3.0115\times10^{23}$ ($6.023\times10^{23}/2$) is now CAVOG/AMDEUT, Avogadro's number from CONSTANTS_MOD over the deuteron's mass in amu (a new DATA value), in the rate and in the cap (both 0.7% smaller); former question 2 is now 1.

1C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 2C DEUTRATE - COMPUTE THE RATE OF NONEQUILIBRIUM DEUTERIUM BURNING 3C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 4C SUBROUTINE DEUTRATE(DL,TL,X,XH2,I,ITLVL) ! KC 2025-05-31 5 SUBROUTINE DEUTRATE(DL,TL,X,I,ITLVL) 6C 10/14/26 MHP via Claude 7C code changed to convert COMMON/OVRTRN/ into OVRTRN_MOD. 8 USE OVRTRN_MOD 9C 10/14/26 MHP via Claude 10C code changed to convert COMMON/DEUTER/ into DEUTER_MOD. 11 USE DEUTER_MOD 12C 10/14/26 MHP via Claude 13C code changed to convert COMMON/CONSTANTS/ into CONSTANTS_MOD. 14 USE CONSTANTS_MOD 15C 10/20/26 MHP via Claude 16C code changed to replace implicit typing with IMPLICIT NONE and 17C explicit declarations. 18 IMPLICIT NONE 19 REAL*8 AMDEUT,C21,DL,RDEUT,RDEUT2,RDEUTMAX,RHO,T9,T9M1,T9M13,T9M23,T9P13, 20 * * T9P23,TFACDEUT,TL,X,Z 21 INTEGER I,ITLVL 22C 9/25/26 MHP via Claude 23C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/, 24C and COMMON/CONST3/ into a single COMMON/CONSTANTS/ 25C 9/24/26 MHP via Claude 26C code changed to remove dead CMIXL from COMMON/CONST3/; it is now 27C passed as an explicit argument where actually used 28C 10/13/26 MHP via Claude: removed the unused COMMON/NMLPAR/ declaration 29C (no member referenced here; residue of the obsolete-option removals). 30 DATA C21/5.240358D-8/ 31C AMDEUT IS THE BARE DEUTERON MASS IN AMU, THE SAME VALUE AS ANUC(3) 32C IN ENGEB.F. KEEP THE TWO IN STEP. 33 DATA AMDEUT/2.013553D0/ 34 SAVE 35C T9P13 IS THE TEMPERATURE IN UNITS OF 10^9 DEGREES K TO THE PLUS 1/3 36C POWER. MINUS IS DENOTED BY M. HERE T9 IS THE TEMPERATURE IN UNITS 37C OF 10^9 K, CONVERTED FROM THE LOG_10 (T) AND RHO IS THE DENSITY IN 38C CGS UNITS. 39 RHO=EXP(CLN*DL) 40 T9 = EXP(CLN*(TL - 9.0D0)) 41 T9P13 = T9**CC13 42 T9P23 = T9P13**2 43 T9M13=1.0D0/T9P13 44 T9M23=T9M13**2 45 T9M1=1.0D0/T9 46C MHP 5/02 ADD DEUTERIUM BURNING TERM TO THE CODE 47C IF DEUTERIUM IS ABOVE A MINIMUM THRESHOLD VALUE. 48C RDEUT IS THE RATE (EXCLUDING FACTORS OF THE 49C ABUNDANCES) AND QRTDEUT IS THE DERIVATIVE W/R/T T. 50C NOTE THAT SCREENING IS EXCLUDED - REASONABLE GIVEN 51C THE LOW TEMPERATURES INVOLVED. 52 Z = -3.72D0*T9M13 53 TFACDEUT = 1.0D0+0.112D0*T9P13+3.38D0*T9P23+2.65D0*T9 54C FACTOR OF CAVOG/AMDEUT REFLECTS AVAGADROS NUMBER DIVIDED BY THE 55C MASS OF THE DEUTERON IN AMU 56C 10/27/26 MHP via Claude 57C code changed to correct the same two errors as in engeb.f: the 58C leading power of T9 is now T9M23, not T9P23, and 3.0115D23 (that is, 59C 6.023D23/2) is replaced by CAVOG/AMDEUT (CONSTANTS_MOD), here and in 60C RDEUTMAX below. Dividing by AMDEUT rather than 2 matches the 61C deuteron mass ENGEB uses, ANUC(3). 62 RDEUT = RHO*2.240D3*T9M23*EXP(Z)*TFACDEUT*CAVOG/AMDEUT 63C NOW LIMIT DEUTERIUM BURNING IN A SURFACE CZ TO BE ON A TIME SCALE 64C NO SHORTER THAN THE CONVECTIVE OVERTURN TIMESCALE. 65 IF(I.GE.JCZ .AND. TAUCZ.GT.1.0D0)THEN 66 RDEUTMAX = CAVOG/AMDEUT/TAUCZ 67 RDEUT2 = RDEUT*X 68 IF(I.EQ.JCZ)THEN 69C JVS 0712 Commented out write command 70c WRITE(*,911)RDEUT2,RDEUTMAX,TAUCZ 71C 911 FORMAT(1P3E15.8) 72 ENDIF 73 IF(RDEUT2.GT.RDEUTMAX)THEN 74C JVS 0712 Commented out write command 75c WRITE(*,*)RDEUT2,RDEUTMAX 76 IF(X.GT.1.0D-6)THEN 77 RDEUT = RDEUTMAX/X 78 ENDIF 79 ENDIF 80 ENDIF 81 DRATE(I) = X*RDEUT*C21 82 IF(ITLVL.EQ.1)THEN 83 DRATE0(I) = DRATE(I) 84 ENDIF 85 RETURN 86 END

1"""Python version of deutrate.f (YREC 7be4133): the rate of deuterium 2burning out of equilibrium at one point.""" 3 4import math 5 6from pyrec import constants_mod 7from pyrec import deuter_mod 8from pyrec import ovrtrn_mod 9 10# DATA C21/5.240358D-8/: seconds per Gyr over Avogadro's number, probably 11# (Fortran line 30). 12C21 = 5.240358e-8 13# DATA AMDEUT/2.013553D0/: the bare deuteron mass in amu, the same value as 14# ANUC(3) in engeb.f (keep the two in step) (Fortran lines 31–33). 15AMDEUT = 2.013553 16 17 18def deutrate(dl: float, tl: float, x: float, i: int, itlvl: int) -> None: 19 """Rate of 2H(p,gamma)3He at point i, without screening, per Gyr, 20 capped in the surface convection zone so that deuterium burns no 21 faster than the convective overturn time. 22 23 Translates: deutrate.f#DEUTRATE (7be4133) 24 25 Arguments: dl, tl (log rho, log T); x, the hydrogen mass fraction; 26 i, the point; itlvl, the level of iteration. 27 28 Returns nothing: sets deuter_mod.drate[i], and at the first level of 29 iteration also deuter_mod.drate0[i]. 30 31 """ 32 33 # Density and powers of T9 (Fortran lines 39–45). 34 c = constants_mod 35 rho = math.exp(c.cln * dl) 36 t9 = math.exp(c.cln * (tl - 9.0)) 37 t9p13 = t9**c.cc13 38 t9p23 = t9p13**2 39 t9m13 = 1.0 / t9p13 40 t9m23 = t9m13**2 41 t9m1 = 1.0 / t9 # noqa: F841 (unused in the Fortran too) 42 43 # The rate without screening (Fortran lines 52–62). 44 z = -3.72 * t9m13 45 tfacdeut = 1.0 + 0.112 * t9p13 + 3.38 * t9p23 + 2.65 * t9 46 # Avogadro's number over the deuteron's mass in amu. 47 rdeut = (rho * 2.240e3 * t9m23 * math.exp(z) * tfacdeut * c.cavog 48 / AMDEUT) 49 50 # Cap in the surface convection zone by the overturn time, in s 51 # (Fortran lines 63–80). 52 # Question (see the page): what the test taucz > 1 is for. 53 if i >= deuter_mod.jcz and ovrtrn_mod.taucz > 1.0: 54 rdeutmax = c.cavog / AMDEUT / ovrtrn_mod.taucz 55 rdeut2 = rdeut * x 56 if rdeut2 > rdeutmax: 57 if x > 1.0e-6: 58 rdeut = rdeutmax / x 59 60 # The rate per Gyr, and at the first level the start-of-step rate 61 # (Fortran lines 81–84). 62 deuter_mod.drate[i] = x * rdeut * C21 63 if itlvl == 1: 64 deuter_mod.drate0[i] = deuter_mod.drate[i]