tpgrad.f

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

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

Variables on this page (19)

alfmlt, cc13, cdelrl, cgl, cln, clsunl, cmixl2, cmixl3, cmxmlt, iov1, iov2, iovim, ishort, ladov, lsdepth, phmlt, spotf, spotx, tgcut

Purpose

TPGRAD decides whether a point is convective and gives its temperature gradient $\nabla = d\ln T/d\ln P$. It computes the radiative gradient, compares it with the adiabatic one (the Schwarzschild criterion), and in convective regions solves the mixing-length theory's cubic for the actual gradient and the convective velocity; on request it also gives the gradient's derivatives for the Henyey solution. Its comments assume that the equation of state and the opacity have been computed already.

How it fits

Inputs: arguments: $\log T$, $T$, $\log P$, $P$, $\rho$ (D), $\log r$, $\log m$ (SL), $L/L_\odot$ (B), the opacity $\kappa$ (O) and its derivatives, the equation-of-state quantities ($\delta = -$QDT, $c_P$ = QCP, $\nabla_\mathrm{ad}$ = DELA, and derivatives), the rotation factors FPL, FTL, $\log T_\mathrm{eff}$, and the mixing length CMIXL (set for each run by the main program, CMIXL = CMIXLA(NK), and passed down by the callers). Through modules: the constants CDELRL, CMIXL2, CMIXL3, set by SETUPS, and CC13 (CONSTANTS_MOD); the temperature above which MLT is skipped, TGCUT, the overshoot switch LADOV and the spot parameters SPOTF, SPOTX, LSDEPTH (NMLPAR_MOD); the overshoot shells IOV1, IOV2, IOVIM (DPMIX_MOD); ISHORT (IOUNITS_MOD).
Who calls it: the routines that evaluate the physics at a point: PHYSIC, COEFFT, QENV, SCONVEC and STARIN.
Outputs: LCONV, DELR ($\nabla_\mathrm{rad}$), DEL ($\nabla$), VEL (the convective velocity, cm/s), and with LDERIV the logarithmic derivatives of $\nabla$ (QACT, QACP, QACR). It sets ALFMLT, PHMLT, CMXMLT in PUALPHA_MOD to zero (lines 77–79 ↗) and never computes them, so the pulsation output's mixing-length columns are zeros: ENVINT writes the three variables, and WRTMOD writes literal zeros in their place, "which were always zero" (settled from the code, Nov 9, 2026: wrtmod.f, lines 291–296 ↗, coefft.f, lines 500–502 ↗, envint.f, lines 590–592 ↗).

The method

Radiative gradient. With CDELRL = $-\log(4\pi/3) - \log\sigma - \log 16$ (from SETUPS), $$\nabla_\mathrm{rad} = \frac{3\,\kappa\,L\,P}{64\pi\,\sigma\,G\,m\,T^4}\;\frac{f_T}{f_P},$$ with $L = $ B$\,L_\odot$, and $f_P$, $f_T$ the rotational distortion factors.
Radiative if $\nabla_\mathrm{rad} - \nabla_\mathrm{ad} \le 10^{-6}$: $\nabla = \nabla_\mathrm{rad}$, $v = 0$; but inside an adiabatic overshoot region (LADOV, shell IOVIM between IOV1 and IOV2) $\nabla = \nabla_\mathrm{ad}$. This is the Schwarzschild criterion, without the composition term $\nabla_\mu$; the YREC paper (§2.3.3) says the code "considers the Schwarzschild criterion for instability, not the Ledoux criterion," and semiconvection (LSEMIC, SCONVEC) is the option for $\mu$-gradient regions (settled from the YREC paper, Nov 9, 2026: lines 80–93 ↗).

Question 1 (line 83 ↗). The threshold $10^{-6}$. A point is radiative when $\nabla_\mathrm{rad} - \nabla_\mathrm{ad} \le 10^{-6}$, not $\le 0$. Is the margin meant to keep LC from flickering at a zone's edge from one iteration to the next, and does it shift the edges measurably? 3. Convective: above TGCUT (cores), $\nabla = \nabla_\mathrm{ad}$ and $v = 10^{-11}$, without MLT. The YREC paper (§2.3.3) gives the reason: "In convective cores the specific heat is quite high and there is an option to skip MLT calculations; in this case the adopted temperature gradient is the adiabatic one." The $v = 10^{-11}$ cm/s set there is not a convective velocity, so VEL in these regions should not be used for overturn times or mixing.

Question 2 (lines 97–107 ↗). Skipping MLT in cores. Above TGCUT the gradient is set to $\nabla_\mathrm{ad}$ without solving the cubic. Deep in a convective core, $\nabla - \nabla_\mathrm{ad}$ would be tiny (of order $10^{-7}$ or less), so the skip should change the structure negligibly while saving time: Is that the reasoning, and has the effect been measured? 4. Spots (G. Somers 9/14): the radiative flux in convective regions is reduced, $\nabla_\mathrm{rad} \to \nabla_\mathrm{rad}/(f x^4 + 1 - f)$, with $f$ = SPOTF; the contrast $x$ is SPOTX, or with LSDEPTH varies with depth as the comment derives, $x(T) = 1 - (1 - x_\mathrm{surf})\,T_\mathrm{amb}/T$. The correction comes after the point is found convective (but see bug 1).

The returned DELR stays uncorrected, while DEL and VEL use the corrected gradient; the derivatives of DELDEL do not follow the correction (suspected bug 2). 5. Mixing-length theory. With $g = Gm/r^2$, $H_p = P/(\rho g)$, the mixing length $\ell = \alpha H_p$, with $\alpha$ = CMIXL, $\phi = \alpha\,\kappa\rho H_p$ (the optical thickness of an element), CMIXL2 = 1/3 and CMIXL3 = $16\sqrt2\,\sigma$: $$A_1 = \frac{16\sqrt2\,\sigma\,T^3\kappa}{c_P\,(1 + \phi^2/3)\,\sqrt{(\nabla_\mathrm{rad} - \nabla_\mathrm{ad})\,g\,\delta/H_p}},\qquad A_3 = \frac{3\phi^2}{4\,(1 + \phi^2/3)\,A_1},$$ and the cubic $$A_3 V^3 + V^2 + A_1 V - 1 = 0$$ is solved for $V$ by a Newton iteration with a second-order correction (at most 25 steps, tolerance $10^{-10}$), starting from $V = 1/A_1$, or $V = A_3^{-1/3}$ when $A_3 > 10^3$. Then $$\nabla = \nabla_\mathrm{ad} + (\nabla_\mathrm{rad} - \nabla_\mathrm{ad})\,V(V + A_1),\qquad v = V\,\sqrt{\frac{g\,\delta\,\alpha^2 H_p\,(\nabla_\mathrm{rad} - \nabla_\mathrm{ad})}{8}} .$$ This is mixing-length theory (Böhm-Vitense 1958, as the YREC paper, §2.3.3, says). From the code's constants, for $\phi \gg 1$ the factor $1/(1 + \phi^2/3)$ drops out and $A_1 \to 2U/\sqrt{\nabla_\mathrm{rad} - \nabla_\mathrm{ad}}$, with $U = \dfrac{3acT^3}{c_P\rho^2\kappa\ell^2}\sqrt{\dfrac{8H_p}{g\delta}}$ (with $\ell = \alpha H_p$ in cm, $ac = 4\sigma$); $V$ is then the convective velocity in units of $\sqrt{g\delta\ell^2(\nabla_\mathrm{rad} - \nabla_\mathrm{ad})/(8H_p)}$. If the cubic does not converge, the point is taken as adiabatic, with a message. If $A_1$'s square root would be of a negative number, the run stops (STOP 912, "RUN ENDED"). The routine also computes DELPM and RRR, which nothing reads (unused?).

Question 3 (lines 139–147 ↗; line 161 ↗). Which variant of the theory. The optically thin factor $1/(1 + \phi^2/3)$ and CMIXL3 $= 16\sqrt2\,\sigma$ fix the variant of mixing-length theory; the calibrated $\alpha$ depends on it. Which published formulation do these constants follow? The factor looks like the correction for optically thin elements of the kind Henyey, Vardya & Bodenheimer (1965) introduced, and $U$ above like the efficiency in Kippenhahn, Weigert & Weiss's textbook: Are they? 6. Derivatives (LDERIV): of $A_1$, $A_3$ and $V$ with respect to $\ln T$, $\ln P$, $\ln r$, combined into those of $\nabla$.

From YREC's documentation

Physics Options: above $\log T$ = TGCUT (default 6.9) convection is taken as adiabatic; LADOV (experimental) makes overshoot regions adiabatic. Starspots: SPOTF is the filling factor (default 0; 0.25 is the bound for active main-sequence stars), SPOTX the spots' temperature contrast (default 1; 0.85 is typical of sunspots; change it only if SPOTF > 0), and LSDEPTH lowers the contrast with depth, assuming a constant field, with no dramatic effect on the results.

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. (never occurred) Spots can stop the run on a point that is convective without them.

(lines 82–83 ↗; lines 127–135 ↗; lines 145–156 ↗) The Schwarzschild test uses $\nabla_\mathrm{rad}$ without the spot correction. For a point found convective, DELDEL is then recomputed with the spotted $\nabla_\mathrm{rad}/(f x^4 + 1 - f)$, which is smaller; if this makes DELDEL zero or negative (a point near a zone's edge, where $\nabla_\mathrm{rad} - \nabla_\mathrm{ad}$ is small), TEST is not positive and the run stops (STOP 912, "RUN ENDED"). Only with spots (SPOTF ≠ 0; the default is 0).

A possible fix (one choice: treat the point as radiative with the spotted gradient):

tpgrad.f

@@ lines 131–137 @@
          ELSE
             DEEPX = SPOTX
          ENDIF
+C 11/09/26 KC via Claude
+C code changed to treat a point as radiative when the spot correction
+C makes it stable, instead of stopping at the square root below.
          DELDEL = DELR/(SPOTF * DEEPX**4.0 + 1.0 - SPOTF) - DELA
+         IF(DELDEL.LE.1.0D-6)THEN
+            LCONV = .FALSE.
+            DEL = DELR/(SPOTF * DEEPX**4.0 + 1.0 - SPOTF)
+            VEL = 0.0D0
+            GOTO 200
+         ENDIF
       ENDIF
 C G Somers END
       G = DEXP(CLN*(CGL + SL - RL - RL))

Checked by running YREC (Nov 9, 2026; 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), one is spotted, testsuite/Test_m0p1feh+0p0GS98_spot25_ZAMS: The spot correction was applied 1,418,129 times, at $\log T$ = 3.46 to 6.66, and DELDEL stayed at 0.0788 or more, so it never turned a point stable and the run never stopped there. YREC's other spotted cases (testsuite/Test_m0p3feh+0p0GS98_spot25_ZAMS and testsuite/Test_m1p0feh+0p0GS98_spot25_ZAMS) were not run; reaching it needs a spotted point near a zone's edge where the correction makes DELDEL vanish. (KNOWN_ISSUES.md mentions the spot recomputation as a place to check for a DELDEL that is NaN, not this case.)

2. (changes results, small) Checked by KC (Nov 10, 2026) With spots, the derivatives of the gradient are those of the unspotted one.

(lines 127–135 ↗; lines 201–202 ↗) With spots, DELDEL is $\nabla_\mathrm{rad}/S - \nabla_\mathrm{ad}$, with $S = f x^4 + 1 - f$ (line 134 ↗), but its logarithmic derivatives QDDELT and QDDELP (lines 201–202 ↗) put the unspotted DELR in their numerators and divide by the spotted DELDEL, so they are too small by about the factor $S$ ($S < 1$) where DELDEL is dominated by $\nabla_\mathrm{rad}$; with LSDEPTH, $x$ depends on $T$, and $d\ln S/d\ln T = 4 f x^3 (1 - x)/S$ is missing too. What follows: The derivatives of the convective gradient (QACT, QACP), which enter the Henyey solution's matrix, are wrong in spotted convective regions; the converged model should change only within the convergence tolerance. Whether the returned DELR should carry the factor $1/S$ as well is a separate choice. Only with spots (SPOTF ≠ 0; the default is 0).

A possible fix:

tpgrad.f

@@ lines 30–36 @@
      *       QA1R,QA1T,QA3P,QA3R,QA3T,QACP,QACR,QACT,QAP,QAT,QCP,QCPP,
      *       QCPT,QDDELP,QDDELT,QDELAP,QDELAT,QDP,QDT,QDTP,QDTT,QOD,QOT,
      *       QVP,QVR,QVT,RL,RRR,SL,T,TEFFL,TEMP1,TEMP2,TEMP3,TEMPOP,
-     *       TEMPOT,TEST,TL,V,VD,VEL,VP,VTOL
+     *       TEMPOT,TEST,TL,V,VD,VEL,VP,VTOL,SPFAC,DSPFT
       INTEGER ITER
       LOGICAL*4 LCONV,LDERIV

@@ lines 133–138 @@
          ENDIF
          DELDEL = DELR/(SPOTF * DEEPX**4.0 + 1.0 - SPOTF) - DELA
       ENDIF
+C 11/10/26 KC via Claude
+C code changed to keep the spot factor S and its derivative d ln S/d ln T
+C (with LSDEPTH, X depends on T), for the derivatives of DELDEL below.
+      SPFAC = 1.0D0
+      DSPFT = 0.0D0
+      IF(SPOTF .NE. 0.00)THEN
+         SPFAC = SPOTF * DEEPX**4.0 + 1.0 - SPOTF
+         IF(LSDEPTH) DSPFT = 4.0D0*SPOTF*DEEPX**3*(1.0D0 - DEEPX)/SPFAC
+      ENDIF
 C G Somers END
       G = DEXP(CLN*(CGL + SL - RL - RL))
       PRESHT = P/(D*G)
@@ lines 198–205 @@
        QDELAP = QAP*DELA
        TEMPOT = QOD*QDT + QOT
        TEMPOP = QOD*QDP
-       QDDELT = ((-4.0D0 + TEMPOT)*DELR - QDELAT)/DELDEL
-       QDDELP = ((+1.0D0 + TEMPOP)*DELR - QDELAP)/DELDEL
+C 11/10/26 KC via Claude
+C code changed to differentiate the spotted DELR/S that DELDEL holds.
+       QDDELT = ((-4.0D0 + TEMPOT - DSPFT)*DELR/SPFAC - QDELAT)/DELDEL
+       QDDELP = ((+1.0D0 + TEMPOP)*DELR/SPFAC - QDELAP)/DELDEL
        TEMP1 = 2.0D0*CMIXL2*PHI2*PHIPHI
        QA1T = TEMPOT*(1D0-TEMP1)-QCPT-0.5D0*(QDT+QDDELT+QDTT/QDT)+3D0
        QA1P = TEMPOP*(1D0-TEMP1)-TEMP1-QCPP-0.5D0

How the code itself shows it. Without spots, DELDEL = DELR - DELA (line 82 ↗), and the numerators (-4.0D0 + TEMPOT)*DELR - QDELAT and (+1.0D0 + TEMPOP)*DELR - QDELAP (lines 201–202 ↗) are its derivatives with respect to $\ln T$ and $\ln P$ ($\nabla_\mathrm{rad} \propto \kappa P/T^4$, QDELAT = QAT*DELA), so QDDELT and QDDELP are $d\ln$DELDEL$/d\ln T$ and $d\ln$DELDEL$/d\ln P$, which QA1T and QA1P use as such (-0.5D0*(QDT+QDDELT+QDTT/QDT), lines 204–206 ↗; $A_1 \propto$ DELDEL$^{-1/2}$, line 147 ↗). With spots, DELDEL = DELR/(SPOTF * DEEPX**4.0 + 1.0 - SPOTF) - DELA (line 134 ↗), and the numerators are left unchanged.

Checked by running YREC (Nov 10, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): In testsuite/Test_m1p0feh+0p0GS98_spot25_ZAMS (1 $M_\odot$, SPOTF = 0.25, SPOTX = 0.85, LSDEPTH), the derivatives were computed 974,865 times in spotted convective points, at $\log T$ = 3.64 to 6.57, with $S$ from 0.878 to 0.99982, so the condition held at every one. With the fix, $\log L$ moves by up to 0.00055 dex and $\log T_\mathrm{eff}$ by 0.000077 (at 0.038 Gyr), against 0.000051 and 0.0000088 for 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}$, on the same case); central $\log T$ by 0.000096 against 0.0000034, and central $X$ by 0.0000095 against 0.00000036; the track has 277 models in both, and ends at 0.1598 Gyr instead of 0.1596. The effect is 9 to 28 times the null change's, small in size, as expected of a change to the matrix alone. YREC's other spotted cases (testsuite/Test_m0p1feh+0p0GS98_spot25_ZAMS and testsuite/Test_m0p3feh+0p0GS98_spot25_ZAMS) were not run.

Definitions

subroutineTPGRAD

lines 6–227 · GitHub

The Schwarzschild test and the temperature gradient: radiative, adiabatic (cores, overshoot), or from the mixing-length cubic, with the convective velocity and, on request, the gradient's derivatives.

Uses modules CONSTANTS_MOD, DPMIX_MOD, IOUNITS_MOD, NMLPAR_MOD, PUALPHA_MOD

Called by COEFFT, PHYSIC, QENV, SCONVEC ×3, STARIN

Compiler warnings (1)

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

line 127: Inequality comparison for REAL(8) [-Wcompare-reals]

Change log

The changes since rp_1

Changed since rp_1. COMMON blocks to modules, IMPLICIT NONE (the unused CTLIM dropped); CMIXL is now an argument instead of a COMMON variable; the cubic's starting value uses CC13 instead of the literal 0.333333333; the negative-square-root message now says "RUN ENDED" instead of "assumed adiabatic convection," which settles former question 2 (stopping is intended), so former questions 3–5 are now 2–4.

Changed since a8a5950. The convective velocity's second sign test is gone: Its radicand is the first test's times $H_p^2$, and a non-positive first test already stops the run (STOP 912), so its branch, which set VEL = 1.0D0-11, that is $-10$, could not be reached (settles former question 1).

1C 2C 3C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 4C TPGRAD 5C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 6 SUBROUTINE TPGRAD(TL,T,PL,P,D,RL,SL,B,O,QDT,QDP,QOT,QOD,QCP, 7 * *DEL,DELR,DELA,QDTT,QDTP,QAT,QAP,QACT,QACP,QACR,QCPT,QCPP, 8 * *VEL,LDERIV,LCONV,FPL,FTL,TEFFL,CMIXL) 9C 10/14/26 MHP via Claude 10C code changed to convert COMMON/PUALPHA/ into PUALPHA_MOD. 11 USE PUALPHA_MOD 12C 10/14/26 MHP via Claude 13C code changed to convert COMMON/DPMIX/ into DPMIX_MOD. 14 USE DPMIX_MOD 15C 10/14/26 MHP via Claude 16C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD. 17 USE NMLPAR_MOD 18C 10/14/26 MHP via Claude 19C code changed to convert COMMON/CONSTANTS/ into CONSTANTS_MOD. 20 USE CONSTANTS_MOD 21C 10/14/26 MHP via Claude 22C code changed to convert COMMON/IOUNITS/ into IOUNITS_MOD. 23 USE IOUNITS_MOD 24C 10/20/26 MHP via Claude 25C code changed to replace implicit typing with IMPLICIT NONE and 26C explicit declarations. 27 IMPLICIT NONE 28 REAL*8 A1,A3,A3P,ATEFFL,B,CMIXL,D,DDEL,DEEPX,DEL,DELA,DELDEL,DELI, 29 * * DELPM,DELR,FPL,FTL,G,O,P,PHI,PHI2,PHIPHI,PL,PRESHT,QA1P, 30 * * QA1R,QA1T,QA3P,QA3R,QA3T,QACP,QACR,QACT,QAP,QAT,QCP,QCPP, 31 * * QCPT,QDDELP,QDDELT,QDELAP,QDELAT,QDP,QDT,QDTP,QDTT,QOD,QOT, 32 * * QVP,QVR,QVT,RL,RRR,SL,T,TEFFL,TEMP1,TEMP2,TEMP3,TEMPOP, 33 * * TEMPOT,TEST,TL,V,VD,VEL,VP,VTOL 34 INTEGER ITER 35 LOGICAL*4 LCONV,LDERIV 36 37C DL,OL,X,Z,LOCOND USED BY OPACTY 38C COMPUTES RADIATIVE GRADIENT AND COMPARES WITH ADIABATIC GRADIENT 39C COMPUTES CONVECTIVE GRADIENT VIA MIXING LENGTH THEORY IF APPLICABLE 40C ASSUMES EQSTAT AND OPACTY HAVE BEEN CALLED 41C RETURNS LCONV = .TRUE. IF CONVECTIVE 42C DELR = RADIATIVE GRADIENT 43C DEL = ACTUAL GRADIENT 44C QACT,QACP = NAT-LOG DERIVATIVES OF THE CONVECTIVE GRADIENT 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 9/25/26 MHP via Claude 52C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/, 53C and COMMON/CONST3/ into a single COMMON/CONSTANTS/ 54C 9/24/26 MHP via Claude 55C code changed to pass CMIXL as an explicit argument instead of via 56C common block, to reduce implicit global state ahead of the F90 57C module conversion 58C 10/12/26 MHP via Claude 59C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 60C namelist variables that the physics routines read from 27 COMMON 61C blocks which also held derived or working members. Those blocks keep 62C only their non-namelist members; five of them -- NEWMX, VMULT, 63C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 64C Derivation-only namelist variables from the same blocks went to 65C COMMON/NMLBYPASS/ instead, not here. 66C 67C Declared identically in all 121 files that hold any of these members. 68C Generated from one member list, so every copy is byte-identical; the 69C linker's 'size of symbol changed' warning is the check that they are. 70C 10/13/26 MHP via Claude 71C code changed to drop COMMON/CTLIM/, declared here but never used. 72C DBG 7/95 To store variables for pulse output 73C G Somers END 74 DATA VTOL/1.0D-10/ 75 SAVE 76 77 ALFMLT=0.0D0 78 PHMLT=0.0D0 79 CMXMLT=0.0D0 80 DELR = O*B*DEXP(CLN*(PL - SL - 4D0*TL + CLSUNL - CGL + CDELRL))* 81 * *FTL/FPL 82 DELDEL = DELR - DELA 83 IF(DELDEL.LE.1.0D-6) THEN 84C ZONE IS RADIATIVE 85 LCONV = .FALSE. 86 DEL = DELR 87 VEL=0.0D0 88 IF (LADOV .AND. IOVIM.GE.IOV1 .AND. IOVIM.LE.IOV2 89 * * .AND. IOVIM.NE.-1) THEN 90 DEL = DELA 91 END IF 92 GOTO 200 93 ENDIF 94C ZONE IS CONVECTIVE 95 LCONV = .TRUE. 96 DEL = DELA 97C SKIP MIXING LENGTH THEORY FOR CORES 98 IF(TL.GT.TGCUT) THEN 99 VEL = 1.0D-11 100 IF(LDERIV) THEN 101C DERIVATIVES OF CONVECTIVE GRADIENT 102 QACT = QAT 103 QACP = QAP 104 QACR = 0.0D0 105 ENDIF 106 GOTO 200 107 ENDIF 108C G Somers 9/14, Add the ability to include spots, which alter 109C the radiative flux in the convective reigons. This is done by 110C setting the flux to Lspotted = L/[f*x^4 + (1-f)], where f is the 111C spot filling factor and x is the temperature contrast between 112C the spotted surface and the normal surface (0 <= x <= 1). We also 113C alter the surface boundary condition used in envint.f, by looking 114C up the pressure at the un-spotted T (ATEFFL) instead of Teff. 115C 116C We have also included the ability to have the temperature contrast 117C SPOTX vary with depth, as is the case in stars. In this case, SPOTX 118C is locally determined by X(T) = 1 - (1-Xsurf)(ATeff/T), where ATeff 119C is the surface temperature of the unspotted regions. This is the 120C analytic expression resulting from a user specified Xsurf, the 121C assumption of magnetic equipartition at the surface (Pmag = Psurf/2), 122C and the assumption of constant Pmag with depth. 123C 124C This flux alters DELR, so recalculate DELDEL with the correction. Only 125C do this if the spot filling factor is non-zero, and the envelope is 126C convective. 127 IF(SPOTF .NE. 0.00)THEN 128 IF(LSDEPTH)THEN 129 ATEFFL = TEFFL - 0.25*LOG10(SPOTF * SPOTX**4.0 + 1.0 - SPOTF) 130 DEEPX = 1.0 - (1.0 - SPOTX)*(10.**ATEFFL)/(10.**TL) 131 ELSE 132 DEEPX = SPOTX 133 ENDIF 134 DELDEL = DELR/(SPOTF * DEEPX**4.0 + 1.0 - SPOTF) - DELA 135 ENDIF 136C G Somers END 137 G = DEXP(CLN*(CGL + SL - RL - RL)) 138 PRESHT = P/(D*G) 139 PHI = CMIXL*D*O*PRESHT 140 PHI2 = PHI*PHI 141 PHIPHI = 1.0D0/(1.0D0 + CMIXL2*PHI2) 142C SOLVE CUBIC A3*V**3 + V**2 + A1*V - 1.0 = 0.0 143C ENSURE THAT THE SQUARE ROOT BEING TAKEN IS THAT OF A 144C POSITIVE NUMBER 145 TEST = DELDEL*G*(-QDT)/PRESHT 146 IF(TEST.GT.0.0D0)THEN 147 A1 = CMIXL3*PHIPHI*T**3*O/(QCP*DSQRT(TEST)) 148 ELSE 149C PATHOLOGICAL FAILURE IN THE CUBIC SOLVER. 150 WRITE(*,201)DELDEL,QDT,TEST 151 201 201 FORMAT('***WARNING***'/'TPGRAD TRIED TO TAKE SQUARE ROOT', 152 * * ' OF NEGATIVE NUMBER - RUN ENDED'/ 153 * * ' DELDEL ',1PE12.3,' QDT ',E12.3,' DISC. ',E12.3) 154 WRITE(*,911)O,B,PL,SL,TL,FTL,FPL,DELR,DELA,RL 155 911 911 FORMAT(1P10E12.3) 156 STOP 912 157C GOTO 15 158 ENDIF 159C A1 = CMIXL3*PHIPHI*T**3*O/(QCP*DSQRT(DELDEL*G*(-QDT)/PRESHT)) 160 V = 1.0D0/A1 161 A3 = 0.75D0*PHI2*PHIPHI/A1 162 A3P = 3.0D0*A3 163 IF(A3.GT.1.0D+3) V = A3**(-CC13) 164 DO 10 ITER = 1,25 165 V = DMIN1(V,1.0D0) 166 VP = A1 + V*(2.0D0 + V*A3P) 167 VD = (-1.0D0 + V*(A1 + V*(1.0D0 + V*A3)))/VP 168 VD = VD*(1.0D0 + VD*(1.0D0 + V*A3P)/VP) 169 V = V - VD 170 IF(DABS(VD).LT.VTOL) GOTO 30 171 10 10 CONTINUE 172C 15 CONTINUE 173 WRITE(ISHORT,20) PL,TL,O,QCP,QDT 174 20 20 FORMAT(' -----CUBIC NON-CONVERGENCE(PL,TL,CAPPA,CP,QDT)=' 175 * * ,2F10.6/1P3E12.6) 176C ASSUME ADIABATIC CONVECTION 177 VEL = 1.0D-11 178 IF(LDERIV) THEN 179C DERIVATIVES OF CONVECTIVE GRADIENT 180 QACT = QAT 181 QACP = QAP 182 QACR = 0.0D0 183 ENDIF 184 GOTO 200 185 30 30 DDEL = DELDEL*V*(V+A1) 186 DEL = DELA + DDEL 187C CALCULATE CONVECTIVE VELOCITY 188C 10/30/26 MHP via Claude 189C code changed to drop a dead sign test. The radicand is the TEST above 190C times PRESHT**2, and a non-positive TEST already stopped (STOP 912). 191C Its else branch set VEL = 1.0D0-11, which is -10, not 1.0D-11. 192 VEL = V*DSQRT(G*(-QDT)*(CMIXL**2)*PRESHT*DELDEL/8.0D0) 193 DELPM=DEL-V*V*DELDEL 194 RRR = 10.0D0**RL 195 IF(LDERIV) THEN 196C DERIVATIVES OF CONVECTIVE GRADIENT 197 QDELAT = QAT*DELA 198 QDELAP = QAP*DELA 199 TEMPOT = QOD*QDT + QOT 200 TEMPOP = QOD*QDP 201 QDDELT = ((-4.0D0 + TEMPOT)*DELR - QDELAT)/DELDEL 202 QDDELP = ((+1.0D0 + TEMPOP)*DELR - QDELAP)/DELDEL 203 TEMP1 = 2.0D0*CMIXL2*PHI2*PHIPHI 204 QA1T = TEMPOT*(1D0-TEMP1)-QCPT-0.5D0*(QDT+QDDELT+QDTT/QDT)+3D0 205 QA1P = TEMPOP*(1D0-TEMP1)-TEMP1-QCPP-0.5D0 206 * * *(QDP+QDDELP+QDTP/QDT) + 0.5D0 207 QA1R = -2.0D0*TEMP1 + 2.0D0 208 QA3T = 2.0D0*PHIPHI*TEMPOT - QA1T 209 QA3P = 2.0D0*PHIPHI*(1.0D0+TEMPOP) - QA1P 210 QA3R = 4.0D0*PHIPHI - QA1R 211 TEMP1 = A1*V 212 TEMP2 = A3*V*V*V 213 TEMP3 = 1.0D0/VP 214 QVT = -(TEMP1*QA1T + TEMP2*QA3T)*TEMP3 215 QVP = -(TEMP1*QA1P + TEMP2*QA3P)*TEMP3 216 QVR = -(TEMP1*QA1R + TEMP2*QA3R)*TEMP3 217 TEMP1 = DELDEL*(A1+V+V) 218 TEMP2 = DELDEL*A1*V 219 DELI = 1.0D0/DEL 220 QACT = (QDELAT+QDDELT*DDEL+TEMP1*QVT+TEMP2*QA1T)*DELI 221 QACP = (QDELAP+QDDELP*DDEL+TEMP1*QVP+TEMP2*QA1P)*DELI 222 QACR = (TEMP1*QVR+TEMP2*QA1R)*DELI 223 ENDIF 224 200 200 CONTINUE 225 226 RETURN 227 END

1"""Python version of tpgrad.f (YREC 7be4133): the Schwarzschild test, the 2temperature gradient and the convective velocity at a point.""" 3 4import math 5 6from pyrec import constants_mod 7from pyrec import dpmix_mod 8from pyrec import iounits_mod 9from pyrec import nmlpar_mod 10from pyrec import pualpha_mod 11 12VTOL = 1.0e-10 # DATA VTOL: the cubic's tolerance 13 14 15def tpgrad(tl: float, t: float, pl: float, p: float, d: float, rl: float, sl: float, b: float, o: float, qdt: float, 16 qdp: float, qot: float, qod: float, qcp: float, dele: float, delr: float, dela: float, qdtt: float, 17 qdtp: float, qat: float, qap: float, qact: float, qacp: float, qacr: float, qcpt: float, qcpp: float, 18 vel: float, lderiv: bool, lconv: bool, fpl: float, ftl: float, teffl: float, cmixl: float) -> tuple[float, 19 float, float, float, float, float, bool]: 20 """Decide whether a point is convective and give its temperature 21 gradient: the radiative gradient, the Schwarzschild test against the 22 adiabatic one, and in convective regions the mixing-length cubic for 23 the actual gradient and the convective velocity; with lderiv, the 24 gradient's logarithmic derivatives for the Henyey solution. 25 26 Translates: tpgrad.f#TPGRAD (7be4133) 27 28 Arguments: log T, T, log P, P, rho (d), log r, log m (sl), L/Lsun (b), 29 the opacity o and the equation-of-state quantities (qdt = -delta, qdp, 30 qot, qod, qcp = c_P, dela = nabla_ad, and their derivatives qdtt, qdtp, 31 qat, qap, qcpt, qcpp), lderiv, the rotation factors fpl, ftl, log Teff 32 and cmixl (the mixing length in pressure scale heights). dele (the 33 Fortran's DEL), delr, qact, qacp, qacr, vel and lconv are outputs; the 34 derivatives keep their values where the Fortran does not set them. 35 36 Returns (dele, delr, qact, qacp, qacr, vel, lconv). 37 38 """ 39 40 c = constants_mod 41 nml = nmlpar_mod 42 43 # The pulsation output's MLT values, always zero (Fortran lines 77–79). 44 pualpha_mod.alfmlt = 0.0 45 pualpha_mod.phmlt = 0.0 46 pualpha_mod.cmxmlt = 0.0 47 48 # The radiative gradient and the Schwarzschild test (Fortran lines 49 # 80–83). 50 # Question (see the page): the margin 1e-6. 51 delr = o * b * math.exp(c.cln * (pl - sl - 4.0 * tl + c.clsunl - c.cgl 52 + c.cdelrl)) * ftl / fpl 53 deldel = delr - dela 54 if deldel <= 1.0e-6: 55 # Radiative (adiabatic inside an adiabatic overshoot region) 56 # (Fortran lines 84–92). 57 lconv = False 58 dele = delr 59 vel = 0.0 60 if (nml.ladov and dpmix_mod.iovim >= dpmix_mod.iov1 61 and dpmix_mod.iovim <= dpmix_mod.iov2 62 and dpmix_mod.iovim != -1): 63 dele = dela 64 return dele, delr, qact, qacp, qacr, vel, lconv # label 200 65 66 # Convective (Fortran lines 94–96). 67 lconv = True 68 dele = dela 69 70 # Cores: no mixing-length theory (Fortran lines 97–107). 71 # Question (see the page): skipping MLT in cores. 72 if tl > nml.tgcut: 73 vel = 1.0e-11 74 if lderiv: 75 qact = qat 76 qacp = qap 77 qacr = 0.0 78 return dele, delr, qact, qacp, qacr, vel, lconv 79 80 # Spots reduce the radiative flux in convective regions (G Somers 9/14) 81 # (Fortran lines 108–136). 82 # Question (see the page): delr itself and its derivatives are left 83 # uncorrected. 84 # Suspected bug (see the page): the radiative/convective test is made 85 # without the correction, and a deldel the correction makes negative 86 # stops the run below. 87 if nml.spotf != 0.00: 88 if nml.lsdepth: 89 ateffl = teffl - 0.25 * math.log10( 90 nml.spotf * nml.spotx**4.0 + 1.0 - nml.spotf) 91 deepx = 1.0 - (1.0 - nml.spotx) * (10.0**ateffl) / (10.0**tl) 92 else: 93 deepx = nml.spotx 94 deldel = delr / (nml.spotf * deepx**4.0 + 1.0 - nml.spotf) - dela 95 96 # Gravity, the scale height and the element's optical thickness (Fortran 97 # lines 137–141). 98 # Question (see the page): which variant of the theory. 99 g = math.exp(c.cln * (c.cgl + sl - rl - rl)) 100 presht = p / (d * g) 101 phi = cmixl * d * o * presht 102 phi2 = phi * phi 103 phiphi = 1.0 / (1.0 + c.cmixl2 * phi2) 104 105 # Solve the cubic a3 v^3 + v^2 + a1 v - 1 = 0: a1, or a stop on a 106 # negative square root (Fortran lines 142–159). 107 test = deldel * g * (-qdt) / presht 108 if test > 0.0: 109 a1 = c.cmixl3 * phiphi * t**3 * o / (qcp * math.sqrt(test)) 110 else: 111 print(f"***WARNING***\nTPGRAD TRIED TO TAKE SQUARE ROOT" 112 f" OF NEGATIVE NUMBER - RUN ENDED\n" 113 f" DELDEL {deldel:12.3E} QDT {qdt:12.3E} DISC. {test:12.3E}") 114 print("".join(f"{v:12.3E}" for v in 115 (o, b, pl, sl, tl, ftl, fpl, delr, dela, rl))) 116 raise SystemExit(912) 117 118 # The Newton iteration (Fortran lines 160–171). 119 v = 1.0 / a1 120 a3 = 0.75 * phi2 * phiphi / a1 121 a3p = 3.0 * a3 122 if a3 > 1.0e+3: 123 v = a3**(-c.cc13) 124 converged = False 125 for it in range(1, 26): # loop 10 126 v = min(v, 1.0) 127 vp = a1 + v * (2.0 + v * a3p) 128 vd = (-1.0 + v * (a1 + v * (1.0 + v * a3))) / vp 129 vd = vd * (1.0 + vd * (1.0 + v * a3p) / vp) 130 v = v - vd 131 if abs(vd) < VTOL: 132 converged = True 133 break # to label 30 134 135 # No convergence: assume adiabatic convection (Fortran lines 172–184). 136 if not converged: 137 iounits_mod.write(iounits_mod.ishort, 138 f" -----CUBIC NON-CONVERGENCE(PL,TL,CAPPA,CP,QDT)=" 139 f"{pl:10.6f}{tl:10.6f}\n" 140 f"{o:12.6E}{qcp:12.6E}{qdt:12.6E}") # format 20 141 # Assume adiabatic convection. 142 vel = 1.0e-11 143 if lderiv: 144 qact = qat 145 qacp = qap 146 qacr = 0.0 147 return dele, delr, qact, qacp, qacr, vel, lconv 148 149 # label 30: the gradient (Fortran lines 185–186). 150 ddel = deldel * v * (v + a1) 151 dele = dela + ddel 152 153 # The convective velocity; its radicand is positive, since a 154 # non-positive test above stopped (MHP 10/30/26) (Fortran lines 155 # 187–194). 156 vel = v * math.sqrt(g * (-qdt) * (cmixl**2) * presht * deldel / 8.0) 157 delpm = dele - v * v * deldel # noqa: F841 (unused in the Fortran too) 158 rrr = 10.0**rl # noqa: F841 (unused in the Fortran too) 159 160 # Derivatives of the convective gradient (Fortran lines 195–223). 161 if lderiv: 162 qdelat = qat * dela 163 qdelap = qap * dela 164 tempot = qod * qdt + qot 165 tempop = qod * qdp 166 qddelt = ((-4.0 + tempot) * delr - qdelat) / deldel 167 qddelp = ((+1.0 + tempop) * delr - qdelap) / deldel 168 temp1 = 2.0 * c.cmixl2 * phi2 * phiphi 169 qa1t = (tempot * (1.0 - temp1) - qcpt 170 - 0.5 * (qdt + qddelt + qdtt / qdt) + 3.0) 171 qa1p = (tempop * (1.0 - temp1) - temp1 - qcpp 172 - 0.5 * (qdp + qddelp + qdtp / qdt) + 0.5) 173 qa1r = -2.0 * temp1 + 2.0 174 qa3t = 2.0 * phiphi * tempot - qa1t 175 qa3p = 2.0 * phiphi * (1.0 + tempop) - qa1p 176 qa3r = 4.0 * phiphi - qa1r 177 temp1 = a1 * v 178 temp2 = a3 * v * v * v 179 temp3 = 1.0 / vp 180 qvt = -(temp1 * qa1t + temp2 * qa3t) * temp3 181 qvp = -(temp1 * qa1p + temp2 * qa3p) * temp3 182 qvr = -(temp1 * qa1r + temp2 * qa3r) * temp3 183 temp1 = deldel * (a1 + v + v) 184 temp2 = deldel * a1 * v 185 deli = 1.0 / dele 186 qact = (qdelat + qddelt * ddel + temp1 * qvt + temp2 * qa1t) * deli 187 qacp = (qdelap + qddelp * ddel + temp1 * qvp + temp2 * qa1p) * deli 188 qacr = (temp1 * qvr + temp2 * qa1r) * deli 189 190 return dele, delr, qact, qacp, qacr, vel, lconv # label 200