tauintnew.f

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

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

Variables on this page (5)

cln, cmsun, crsun, json, taucz

Purpose

TAUINTNEW computes the convective overturn time scale $\tau_\mathrm{cz}$ (TAUCZ) and the radius of the base of the convection zone (RBCZ), by the method G. Somers introduced in 3/17. Its comment: search from the centre up for the base of the surface convection zone, and from the surface down; if both agree the zone is unique and the "old top-down sinkline method" is used; if not, several convective cells exist, and the crossing time of the lowest cell in the envelope is taken.

How it fits

Inputs: the full profile from GETTAU (interior and envelope; M points, of which MINT are interior shells): HS1, LCZ, HR, HP, HD, HG, SVEL, and DEL1, DEL2: for interior shells GETTAU passes SDEL(1,·) and SDEL(3,·), the radiative and adiabatic gradients, but for the stitched envelope points EDELS(1,·) and EDELS(2,·), the radiative and actual gradients.
Who calls it: only GETTAU, when LNEWTCZ is set; otherwise GETTAU calls TAUINT.
What it calls: KSPLINE and KSPLINT, a natural cubic spline through four points.
Outputs: TAUCZ in OVRTRN_MOD (read by the mass-loss and rotation routines, by their names), and the argument RBCZ in solar radii.
Modules: OVRTRN_MOD (TAUCZ), CONSTANTS_MOD (CLN, CMSUN, CRSUN), JSON from DIMENSIONS_MOD.

The method

Find the zone. From shell IBEG (first 1) upward: pass the convective core, if any (its top JCORE, 0 at the start of each call), then the first convective shell above a radiative region is IMAX. From the surface downward, the base of the outermost zone is JMAX: the first point J below a convective point J+1 that is radiative, plus one, so that radiative points at the surface are passed over (JvS 02/2026; if no such point exists, JMAX = 1). A star with no such zone is radiative ($\tau_\mathrm{cz} = 0$); one convective throughout is fully convective.

Question 1 (lines 73–82 ↗; lines 115–119 ↗). JCORE may be stale. JCORE is set only when the upward search leaves a convective region that started at IBEG. It is 0 at the start of each call, but not reset at label 50, so when a restart begins in a radiative shell, JCORE keeps its value from the earlier pass. KC's check avoids an index of 0, but not a stale index. Should JCORE be reset at label 50?

Question 2 (line 74 ↗; lines 85–98 ↗). A fully convective star. The code reaches "fully convective" only when every point from IBEG to M is convective, and then sets IMAX = 1, JCORE = 1. A star convective throughout its interior is not stitched (GETTAU stitches the envelope only when its CHKPRS = HP(IMAX)-HP(M) < 1, and here IMAX = 1: gettau.f, lines 88–103 ↗), so M = MINT and the near-core test is skipped; but after a restart that begins inside a zone reaching M, with the envelope stitched (M ≠ MINT), the near-core test HP(1) - HP(1) < 1 holds, IBEG = 2, and the search starts again, possibly in the same state. And if the outermost points are radiative (point M, which GETTAU's own search skips), the whole convective interior is taken as a convective core, and $\tau_\mathrm{cz} = 0$. How should a fully convective star be handled?
Reject false zones and search again from above: a zone less than one decade in pressure above the core, $\log P(\mathrm{JCORE}) - \log P(\mathrm{IMAX}) < 1$ (unless the profile is the interior only); a zone fewer than 3 shells thick; in the interior-only case with two zones (IJMAX ≠ 0, M = MINT), IMAX is set to JMAX; and, at the end, a result $\tau_\mathrm{cz} > 10^{20}$. For a zone, HSTOP > HSBOT, so LCZSUR is always true (lines 137–150 ↗); HSTOP and HSBOT serve only that test.

Nothing bounds the restarts or keeps their tests inside the profile (suspected bug 2).
The base: $r_b$ where $\nabla_1 - \nabla_2$ changes sign, by linear interpolation in $\log r$ (for a fully convective star, $\log r_b = 1$: The comment's "very small number" is RBCZ = 1., then EXP(CLN*RBCZ), that is 10 cm, harmless next to a star's radius; settled from the code, Nov 9, 2026: lines 157–168 ↗).

Question 3 (lines 163–166 ↗; gettau.f, lines 155–157 ↗; gettau.f, lines 177–179 ↗). Two kinds of gradient. DEL2 is the adiabatic gradient for interior shells but the actual gradient for the stitched envelope points, so the sign change of $\nabla_1 - \nabla_2$ means different things on the two sides of the fitting point. Intended?
One zone (IMAX = JMAX): go up from the base to where $r - r_b = H_p$, with $H_p = P/(\rho g)$. Around that point, four shells give $\delta_k = (r_k - r_b) - H_{p,k}$, and natural cubic splines in $\delta$ evaluated at $\delta = 0$ give $v$ and $H_p$ (the four shells KK-2 to KK+1, with KK kept between 3 and M-1, so that they lie inside the profile); then $$\tau_\mathrm{cz} = \frac{H_p}{v}\bigg|_{\,r - r_b = H_p} .$$ If already at IMAX $r - r_b > H_p$, $\tau_\mathrm{cz} = H_p/v$ there; but the code then goes on to step 5, which replaces it (bug 1). A third spline gives the radius at that point, ENVRCL, which nothing reads (unused?, lines 218–223 ↗).
Several cells, or $r_b + H_p$ outside the star: the lowest zone's top $r_t$ (again where $\nabla_1 - \nabla_2$ changes sign), its width $\Delta r = r_t - r_b$, and $v$ at its middle, $r_b + \Delta r/2$, from a four-point spline in $r$ (the same limits on KK): $$\tau_\mathrm{cz} = \frac{\Delta r}{v(r_b + \Delta r/2)} .$$

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 overturn time at `IMAX` is overwritten.

(lines 176–181 ↗; lines 229–231 ↗; lines 235–272 ↗) When the first convective shell is already more than $H_p$ above the base (VALTEST > 0 at IMAX), the code sets TAUCZ = PSCA/CVEL, but the ENDIF moved on 2025-05-31 to avoid a block mismatch (lines 229–230 ↗) now closes the IF before label 110, so this branch falls through to the averaging method, which sets TAUCZ = DELR/CVEL instead. Before the move, label 110 was inside the ELSE block and the branch went on past it. On the default path (LNEWTCZ true), for one zone (IJMAX = 0) whose first shell lies more than one scale height above the base, as on a coarse grid.

A possible fix:

tauintnew.f

@@ lines 179–184 @@
             CVEL = SVEL(IMAX)
             PSCA = PSCA2
             TAUCZ = PSCA/CVEL
+C 11/09/26 KC via Claude
+C code changed to skip the averaging method, as the branch did before
+C the ENDIF was moved (2025-05-31).
+            GOTO 140
          ELSE
             DO K = IMAX+1,M-1,1
                PSCA1 = PSCA2

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 the 33 runs no model took the branch at lines 176–181 ↗, that is, none had its first convective shell more than one pressure scale height above the base of a single zone. The fix was not run. What would reach it: a single surface zone with few shells near its base. Intent: The move is marked "KC 2025-05-31 MOVED ENDIF HERE TO AVOID BLOCK MISMATCH" (lines 229–230 ↗); it made the jump to label 110 legal, and the fall-through looks like a side effect of the move, not a design.

2. (never occurred) The restarts can read past the profile, and can loop forever.

(line 75 ↗; lines 115–131 ↗; lines 284–287 ↗; lines 92–98 ↗) The searches restart with GOTO 50 for a zone near the core (lines 115–118 ↗), a zone thinner than 3 shells (lines 124–131 ↗), and $\tau_\mathrm{cz} > 10^{20}$ (lines 284–287 ↗). Each restart raises IBEG, but nothing stops it from passing M, after which LCZ(IBEG) is read beyond the profile (line 75 ↗), and the 3-shell test reads LCZ(IMAX+2) beyond M when IMAX ≥ M-1, and LCZ(IMAX+1) too when IMAX = M (line 124 ↗). Past its MM points, GETTAU's array LCF holds values from earlier calls (gettau.f, line 59 ↗; gettau.f, line 158 ↗; gettau.f, line 180 ↗). If such a value is true, the fully convective branch sets I = 1 and JCORE = 1 (lines 92–98 ↗) and the tests that follow send the search back to IBEG = 2, where it can come round again; and a profile convective at every point, the envelope's surface included, with the envelope stitched (M ≠ MINT) restarts at IBEG = 2 without end, since the near-core test then holds at once, $\log P_1 - \log P_1 = 0 < 1$ (lines 115–118 ↗). Only on the default path (LNEWTCZ on), for a profile whose only convective shells above the core lie within two points of the surface, or a stitched profile convective at every point, the envelope's surface included.

A possible fix:

tauintnew.f

@@ lines 20–26 @@
       REAL*8 CVEL,DD1,DD2,DEL1,DEL2,DELR,ENVRCL,FX,HD,HG,HP,HR,HS1,
      *       HSBOT,HSTOP,PSCA,PSCA1,PSCA2,RBCZ,RTCZ,RTEST,RTEST1,RTEST2,
      *       SVEL,VALTEST,XDEL,XRAD,YPSC,YRAD,YSPL,YVEL
-      INTEGER I,IBEG,IJMAX,IMAX,ITOP,J,JCORE,JMAX,K,KK,M,MINT
+      INTEGER I,IBEG,IJMAX,IMAX,ITOP,J,JCORE,JMAX,K,KK,M,MINT,NREST
       LOGICAL*4 LALLCZ,LCZ,LCZSUR,LRADZONE
 C 10/13/26 MHP via Claude: removed the unused COMMON/COMP/ declaration
 C (no member referenced here; residue of the obsolete-option removals).
@@ lines 71–78 @@
 C is radiative never set JCORE, so the test below read the value left
 C from an earlier call (SAVE).
       JCORE = 0
+C 11/10/26 KC via Claude
+C code changed to bound the restarts: once IBEG passes M, or after M
+C restarts (a search that comes round again), the search ends as for
+C a radiative envelope (TAUCZ = 0), and LCZ is not read past M.
+      NREST = 0
    50 LRADZONE = .TRUE.
-      IF (LCZ(IBEG)) LRADZONE = .FALSE.
+      NREST = NREST + 1
+      IF(IBEG.GT.M .OR. NREST.GT.M) IBEG = M+1
+      IF(IBEG.LE.M) THEN
+         IF(LCZ(IBEG)) LRADZONE = .FALSE.
+      ENDIF
       DO 60 I = IBEG,M,1
          IF (.NOT.LRADZONE.AND..NOT.LCZ(I)) THEN
             JCORE = I-1
@@ lines 112–118 @@
 C ENSURE THIS CZ IS SUFFICIENTLY FAR FROM THE CONVECTIVE CORE. IN SOME CASES
 C THERE ARE A FEW SHELLS RIGHT ABOVE THE CORE THAT ARE CONVECTIVE. SKIP THESE.
 C KC 2025-05-30 ADDED IF CHECK TO AVOID RUNTIME ERROR.
-      IF (IMAX .GT. 0 .AND. JCORE .GT. 0) THEN
+C 11/10/26 KC via Claude
+C code changed to skip this test for a fully convective model, which
+C has no core: there HP(JCORE)-HP(IMAX) is 0, and with the envelope
+C stitched (M.NE.MINT) the search went back to IBEG = 2 without end.
+      IF (IMAX .GT. 0 .AND. JCORE .GT. 0 .AND. .NOT.LALLCZ) THEN
          IF(M.NE.MINT .AND. HP(JCORE)-HP(IMAX).LT.1)THEN
             IBEG = IMAX+1
             GOTO 50
@@ lines 121–126 @@
 C ENSURE THIS CZ IS AT LEAST 3 SHELLS THICK. IF NOT, WE ARE LIKELY SNAGGED
 C ON A SINGLE SHELL CONV ZONE DEEP IN THE INTERIOR. RECOMMENSE SEARCH FROM
 C ABOVE THIS SHELL.
+C 11/10/26 KC via Claude
+C code changed to treat a zone with fewer than 3 shells below M as thin
+C without reading LCZ past M.
+      IF(IMAX+2.GT.M) THEN
+         IBEG = M+1
+         GOTO 50
+      ENDIF
       IF(.NOT.(LCZ(IMAX).AND.LCZ(IMAX+1).AND.LCZ(IMAX+2)))THEN
          IF(.NOT.LCZ(IMAX+1))THEN
             IBEG = IMAX+1

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_m1p4feh+0p0GS98_dif_TAMS the search restarted 120 times, at IBEG = 298 to 352 in profiles of 407 to 427 points, never more than 20 times in a call; in the other four cases run (testsuite/Test_m1p0feh+0p0GS98_base_TAMS, testsuite/Test_m0p3feh+0p0GS98_base_ZAMS, testsuite/Test_m0p8feh+0p0GS98_rotp8p0d0p1mag1e4std_TAMS, and examples/run_standard_solar_model/Test_solarcal_noGS_norot) it never restarted. IBEG never passed M, and the 3-shell test never had IMAX + 2 > M, so nothing was read past the profile and every search ended. The fix was not run: Its guards never act in these runs, so it changes nothing there. The fully convective testsuite/Test_m0p3feh+0p0GS98_base_ZAMS (LNEWTCZ and LCALCENV on by default) never restarted, presumably because its profile's outermost points are radiative, so that it is not convective at every point; this was not checked. What would reach it: a profile whose only convective shells above the core lie within two points of the surface, or a stitched profile convective at every point, the envelope's surface included. Intent: none found.

Definitions

subroutineTAUINTNEW

lines 4–289 · GitHub

The newer overturn time (G. Somers 3/17): finds the base of the main surface convection zone from both directions, then $H_p/v$ one pressure scale height above the base (by four-point splines), or the lowest cell's width over its central velocity; also returns the base's radius in solar units.

Uses modules CONSTANTS_MOD, DIMENSIONS_MOD, OVRTRN_MOD

Called by GETTAU

Flowchart of its calls, in order (8)

if (LCZSUR), lines 155–277

else, lines 182–230

for K = IMAX+1,M-1,1, lines 183–228

if (VALTEST.GT.0.0), lines 190–227

KSPLINE 210→KSPLINT 211→KSPLINE 216→KSPLINT 217→KSPLINE 222→KSPLINT 223

KSPLINE 269→KSPLINT 270

Change log

The changes since rp_1

Changed since rp_1. COMMON blocks to modules, IMPLICIT NONE (the unused blocks dropped); the downward search for the base of the outermost zone now requires a convective point above the radiative one (JvS 02/2026), so a radiative layer at the surface no longer makes JMAX = M. No question is settled: question 6 (now the question on a fully convective star) concerns the upward search.

Changed since a8a5950. JCORE starts at 0 in each call, so a model with a radiative centre no longer reads a value left from an earlier call (narrows question 1); the four-point stencils keep KK between 3 and M-1, so they no longer read point M+1 (settles former question 3).

1C$$$$$$ 2C RECURSIVE SUBROUTINE TAUINTNEW(HCOMP,HS2,HS1,LCZ,HR,HP,HD,HG,M,MINT,SVEL, 3C * DEL1,DEL2,HSTOT,RBCZ) ! KC 2025-05-31 4 SUBROUTINE TAUINTNEW(HS1,LCZ,HR,HP,HD,HG,M,MINT,SVEL, 5 * * DEL1,DEL2,RBCZ) 6C 10/15/26 MHP via Claude 7C code changed to take JSON from DIMENSIONS_MOD instead of declaring 8C PARAMETER(JSON=5000) here. 9 USE DIMENSIONS_MOD 10C 10/14/26 MHP via Claude 11C code changed to convert COMMON/OVRTRN/ into OVRTRN_MOD. 12 USE OVRTRN_MOD 13C 10/14/26 MHP via Claude 14C code changed to convert COMMON/CONSTANTS/ into CONSTANTS_MOD. 15 USE CONSTANTS_MOD 16C 10/20/26 MHP via Claude 17C code changed to replace implicit typing with IMPLICIT NONE and 18C explicit declarations. 19 IMPLICIT NONE 20 REAL*8 CVEL,DD1,DD2,DEL1,DEL2,DELR,ENVRCL,FX,HD,HG,HP,HR,HS1, 21 * * HSBOT,HSTOP,PSCA,PSCA1,PSCA2,RBCZ,RTCZ,RTEST,RTEST1,RTEST2, 22 * * SVEL,VALTEST,XDEL,XRAD,YPSC,YRAD,YSPL,YVEL 23 INTEGER I,IBEG,IJMAX,IMAX,ITOP,J,JCORE,JMAX,K,KK,M,MINT 24 LOGICAL*4 LALLCZ,LCZ,LCZSUR,LRADZONE 25C 10/13/26 MHP via Claude: removed the unused COMMON/COMP/ declaration 26C (no member referenced here; residue of the obsolete-option removals). 27C 10/13/26 MHP via Claude: removed the unused COMMON/DEUTER/ declaration 28C (no member referenced here; residue of the obsolete-option removals). 29C 10/13/26 MHP via Claude: removed the unused COMMON/NMLPAR/ declaration 30C (no member referenced here; residue of the obsolete-option removals). 31C 9/25/26 MHP via Claude 32C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/, 33C and COMMON/CONST3/ into a single COMMON/CONSTANTS/ 34C 10/13/26 MHP via Claude 35C code changed to remove COMMON/JTEST/: IMAX1, IMAX2 and LJVS are 36C referenced nowhere in executable code. The only code that used LJVS is 37C a commented-out block in wrtout.f, removed with it. 38C Needs access to this common block: grants knowledge of envellope 39C KC 2025-05-30 reordered common block elements 40C JvS 08/25 Updated with new elements 41C 10/13/26 MHP via Claude: removed the unused COMMON/ENVSTRUCT/ declaration 42C (no member referenced here; residue of the obsolete-option removals). 43C DIMENSION HCOMP(15,JSON),HS2(JSON),LCZ(JSON), ! KC 2025-05-31 44 DIMENSION LCZ(JSON), 45 * * HR(JSON), HP(JSON), HD(JSON), HG(JSON), HS1(JSON), ! XSUM(15) 46 * * DEL1(JSON), DEL2(JSON), SVEL(JSON) 47C G Somers; Adding vectors for cubic spline int. 48 DIMENSION XDEL(4),XRAD(4),YRAD(4),YVEL(4),YPSC(4),YSPL(4) 49 SAVE 50C 51C G Somers 3/17, IMPLEMENTING NEW TAUCZ METHOD. 52C THE METHOD WORKS AS FOLLOWS. SEARCH FROM THE CENTER OF THE STAR UP UNTIL YOU 53C FIND THE BASE OF THE SURFACE CONVECTION ZONE. THEN SEARCH FROM THE SURFACE 54c DOWN UNTIL YOU FIND THE SAME. IF BOTH DIRECTIONS GIVE SAME ANSWER, THEN 55c SURFACE CONVECTION ZONE IS UNIQUE, AND THE OLD TOP-DOWN SINKLINE METHOD CAN 56c BE EMPLOYED. IF THEY ARE DIFFERENT, THEN MULTIPLE SURFACE CONVECTION CELLS 57C EXIST. IN THIS CASE, DETERMINE THE CROSSING TIME OF THE LOWEST CONVECTION 58C CELL IN THE ENVELOPE. 59C 60C IMAX: BASE OF THE MAIN INTERNAL CONVECTION ZONE. 61C JMAX: BASE OF THE UPPER-MOST CONVECTION ZONE. 62C RBCZ: RADIUS AT WHICH THE BCZ EXISTS 63C TAUCZ: OVERTURN TIMESCALE 64C 65C DETERMINE EXTENT OF SURFACE CONVECTION ZONE. 66 LALLCZ = .FALSE. 67C FIND THE TOP OF THE MAIN INTERNAL CONVECTION ZONE. 68 IBEG = 1 69C 10/30/26 MHP via Claude 70C code changed to start with no core (JCORE = 0). A model whose centre 71C is radiative never set JCORE, so the test below read the value left 72C from an earlier call (SAVE). 73 JCORE = 0 74 50 50 LRADZONE = .TRUE. 75 IF (LCZ(IBEG)) LRADZONE = .FALSE. 76 DO 60 I = IBEG,M,1 77 IF (.NOT.LRADZONE.AND..NOT.LCZ(I)) THEN 78 JCORE = I-1 79 LRADZONE = .TRUE. 80 ENDIF 81 IF (LRADZONE.AND.LCZ(I)) GOTO 70 82 60 60 CONTINUE 83C IF THE CODE GETS HERE, EITHER THE STAR IS FULLY CONVECTIVE, OR 84C HAS A COMPLETELY RADIATIVE ENVLOPE, MAYBE WITH A CONVECTIVE CORE. 85 IF(LRADZONE)THEN ! RADIATIVE ENVELOPE 86 LALLCZ = .FALSE. 87 LCZSUR = .FALSE. 88 IMAX = M 89 JMAX = M 90 JCORE = M 91 GOTO 100 92 ELSE ! FULLY CONVECTIVE 93 LALLCZ = .TRUE. 94 LCZSUR = .TRUE. 95 I = 1 96 J = 0 97 JCORE = 1 98 GOTO 90 99 ENDIF 100 70 70 CONTINUE 101C FIND THE BASE OF THE UPPER-MOST CONVECTION ZONE 102C JvS 02/2026: ADDED LCZ(J+1) TO TRAP OUT CASES WHERE THE LAST (FEW) POINT(S) IN 103C THE MODEL IS(ARE) ALREADY RADIATIVE 104 DO 80 J = M-1,1,-1 105 IF(LCZ(J+1) .AND. .NOT.LCZ(J)) GOTO 90 106 80 80 CONTINUE 107 90 90 IMAX = I 108 JMAX = J + 1 109 IJMAX = JMAX-IMAX 110C LOCATING THE CORRECT CZ INCURES VARIOUS PATHOLOGIES. TRAP THEM OUT. 111C 112C ENSURE THIS CZ IS SUFFICIENTLY FAR FROM THE CONVECTIVE CORE. IN SOME CASES 113C THERE ARE A FEW SHELLS RIGHT ABOVE THE CORE THAT ARE CONVECTIVE. SKIP THESE. 114C KC 2025-05-30 ADDED IF CHECK TO AVOID RUNTIME ERROR. 115 IF (IMAX .GT. 0 .AND. JCORE .GT. 0) THEN 116 IF(M.NE.MINT .AND. HP(JCORE)-HP(IMAX).LT.1)THEN 117 IBEG = IMAX+1 118 GOTO 50 119 ENDIF 120 ENDIF 121C ENSURE THIS CZ IS AT LEAST 3 SHELLS THICK. IF NOT, WE ARE LIKELY SNAGGED 122C ON A SINGLE SHELL CONV ZONE DEEP IN THE INTERIOR. RECOMMENSE SEARCH FROM 123C ABOVE THIS SHELL. 124 IF(.NOT.(LCZ(IMAX).AND.LCZ(IMAX+1).AND.LCZ(IMAX+2)))THEN 125 IF(.NOT.LCZ(IMAX+1))THEN 126 IBEG = IMAX+1 127 ELSE 128 IBEG = IMAX+2 129 ENDIF 130 GOTO 50 131 ENDIF 132C IN THE CASE WHERE ONLY THE INTERIOR MODEL IS CONSIDERED, USE JMAX. 133 IF (IJMAX.NE.0.AND.M.EQ.MINT) THEN 134 IMAX = JMAX 135 IJMAX = 0 136 ENDIF 137C HSTOP IS THE MASS AT THE TOP OF THE C.Z. 138C HSBOT IS THE MASS AT THE BOTTOM OF THE C.Z. 139 HSTOP = HS1(M) 140 IF(JMAX.GT.1)THEN 141 HSBOT = 0.5D0*(HS1(JMAX)+HS1(JMAX-1)) 142 ELSE 143 HSBOT = 0.0D0 144 ENDIF 145C LCZSUR=T IF A SURFACE C.Z.DEEP ENOUGH FOR ANGULAR MOMENTUM LOSS EXISTS 146 IF((HSTOP-HSBOT)/CMSUN.GT.0.0D0)THEN 147 LCZSUR = .TRUE. 148 ELSE 149 LCZSUR= .FALSE. 150 ENDIF 151 100 100 CONTINUE 152 153C CALCULATE TAUCZ TOP DOWN FOR FULLY AND PARTIALLY CONVECTIVE. FOR A FULLY 154C RADIATIVE ENVELOPE, SET TAUCZ = 0.0. 155 IF(LCZSUR)THEN 156C INFER RBCZ, THE RADIUS AT THE BASE OF THE SURFACE CZ. 157 IF(LALLCZ)THEN 158C SET TO CENTER OF STAR, I.E. A VERY SMALL NUMBER 159C IN CGS UNITS, RELATIVE TO THE RADII OF STARS. 160 RBCZ = 1. 161 ELSE 162C PINPOINT RCZ LINEARLY 163 DD2 = DEL1(IMAX-1)-DEL2(IMAX-1) 164 DD1 = DEL1(IMAX)-DEL2(IMAX) 165 FX = DD2/(DD2-DD1) 166 RBCZ = HR(IMAX-1)+FX*(HR(IMAX)-HR(IMAX-1)) 167 ENDIF 168 RBCZ = EXP(CLN*RBCZ) 169C IF THERE ARE MULTIPLE CONVECTION CELLS IN THE PREDOMINATELY 170C RADIATIVE ENVELOPE, USE AVERAGING METHOD. 171 IF (IJMAX.NE.0) GOTO 110 172C CALCULATE HP 173 PSCA2 = EXP(CLN*(HP(IMAX)-HD(IMAX)))/HG(IMAX) 174 RTEST2 = EXP(CLN*HR(IMAX))-RBCZ 175 VALTEST = RTEST2 - PSCA2 176 IF(VALTEST.GT.0.0)THEN 177C INNER MOST POINT PSCA IS BELOW BCZ 178 K = IMAX 179 CVEL = SVEL(IMAX) 180 PSCA = PSCA2 181 TAUCZ = PSCA/CVEL 182 ELSE 183 DO K = IMAX+1,M-1,1 184 PSCA1 = PSCA2 185 RTEST1 = RTEST2 186 PSCA2 = EXP(CLN*(HP(K)-HD(K)))/HG(K) 187 RTEST2 = EXP(CLN*HR(K))-RBCZ 188 VALTEST = RTEST2 - PSCA2 189C FIND LOCATION WHERE HP = R 190 IF(VALTEST.GT.0.0)THEN 191C IF K = IMAX+1, MOVE AWAY FROM EDGE CASE. 192 KK = K 193 IF (K.EQ.IMAX+1) KK = K+1 194C 10/30/26 MHP via Claude: keep the 4-point stencil KK-2..KK+1 inside 195C 1..M; next to the surface it read point M+1. 196 KK = MAX(3,MIN(KK,M-1)) 197C FIND V WITH CUBIC SPLINE 198 XDEL(1) = EXP(CLN*HR(KK-2))-RBCZ- 199 * * EXP(CLN*(HP(KK-2)-HD(KK-2)))/HG(KK-2) 200 XDEL(2) = EXP(CLN*HR(KK-1))-RBCZ- 201 * * EXP(CLN*(HP(KK-1)-HD(KK-1)))/HG(KK-1) 202 XDEL(3) = EXP(CLN*HR(KK))-RBCZ- 203 * * EXP(CLN*(HP(KK)-HD(KK)))/HG(KK) 204 XDEL(4) = EXP(CLN*HR(KK+1))-RBCZ- 205 * * EXP(CLN*(HP(KK+1)-HD(KK+1)))/HG(KK+1) 206 YVEL(1) = SVEL(KK-2) 207 YVEL(2) = SVEL(KK-1) 208 YVEL(3) = SVEL(KK) 209 YVEL(4) = SVEL(KK+1) 210 CALL KSPLINE(XDEL,YVEL,YSPL) 211 CALL KSPLINT(XDEL,YVEL,YSPL,0.0D0,CVEL) 212 YPSC(1) = EXP(CLN*(HP(KK-2)-HD(KK-2)))/HG(KK-2) 213 YPSC(2) = EXP(CLN*(HP(KK-1)-HD(KK-1)))/HG(KK-1) 214 YPSC(3) = EXP(CLN*(HP(KK)-HD(KK)))/HG(KK) 215 YPSC(4) = EXP(CLN*(HP(KK+1)-HD(KK+1)))/HG(KK+1) 216 CALL KSPLINE(XDEL,YPSC,YSPL) 217 CALL KSPLINT(XDEL,YPSC,YSPL,0.0D0,PSCA) 218 YRAD(1) = HR(KK-2) 219 YRAD(2) = HR(KK-1) 220 YRAD(3) = HR(KK) 221 YRAD(4) = HR(KK+1) 222 CALL KSPLINE(XDEL,YRAD,YSPL) 223 CALL KSPLINT(XDEL,YRAD,YSPL,0.0D0,ENVRCL) 224C DEFINE TAUCZ 225 TAUCZ = PSCA/CVEL 226 GOTO 140 227 ENDIF 228 END DO 229C KC 2025-05-31 MOVED ENDIF HERE TO AVOID BLOCK MISMATCH. 230 ENDIF 231 110 110 CONTINUE 232c IF CODE GETS HERE, THEN 1 PSCA ABOVE THE BCZ IS OUTSIDE OF THE STAR, 233c OR MULTIPLE CONVECTION CELLS RESIDE WITHIN THE RADIATIVE ENVELOPE. 234C IN THIS CASE, EVALUATE CVEL AT CENTER OF BOTTOM CZ. FIRST, FIND UPPER EDGE. 235 DO K = IMAX,M-1,1 236 IF(.NOT.LCZ(K))GOTO 120 237 ENDDO 238 120 120 ITOP = K-1 239c PINPOINT TOP OF CZ LINEARLY 240 DD2 = DEL1(ITOP)-DEL2(ITOP) 241 DD1 = DEL1(ITOP+1)-DEL2(ITOP+1) 242 FX = DD2/(DD2-DD1) 243 RTCZ = HR(ITOP)+FX*(HR(ITOP+1)-HR(ITOP)) 244 RTCZ = EXP(CLN*RTCZ) 245C CALCULATE WIDTH OF THE CONVECTION CELL 246 DELR = RTCZ - RBCZ 247C FIND INDICES SURROUNDING LINEAR CENTER OF CZ 248 RTEST = LOG10(RBCZ + 0.5*DELR) 249 DO K = IMAX,ITOP+1,1 250 IF(HR(K).GT.RTEST) GOTO 130 251 ENDDO 252 130 130 CONTINUE 253 RTEST = EXP(CLN*RTEST) 254C IF K = IMAX+1, MOVE AWAY FROM EDGE CASE. 255 KK = K 256 IF (K.EQ.IMAX+1) KK = K+1 257C 10/30/26 MHP via Claude: keep the 4-point stencil KK-2..KK+1 inside 258C 1..M; next to the surface it read point M+1. 259 KK = MAX(3,MIN(KK,M-1)) 260C CALCULATE CVEL AT THAT CENTRAL LOCATION 261 XRAD(1) = EXP(CLN*HR(KK-2)) 262 XRAD(2) = EXP(CLN*HR(KK-1)) 263 XRAD(3) = EXP(CLN*HR(KK)) 264 XRAD(4) = EXP(CLN*HR(KK+1)) 265 YVEL(1) = SVEL(KK-2) 266 YVEL(2) = SVEL(KK-1) 267 YVEL(3) = SVEL(KK) 268 YVEL(4) = SVEL(KK+1) 269 CALL KSPLINE(XRAD,YVEL,YSPL) 270 CALL KSPLINT(XRAD,YVEL,YSPL,RTEST,CVEL) 271C DEFINE TAUCZ 272 TAUCZ = DELR/CVEL 273 140 140 CONTINUE 274C KC 2025-05-31 MOVED ENDIF ABOVE TO AVOID BLOCK MISMATCH. 275C ENDIF 276C CONVERT CORE RADIUS INTO SOLAR UNITS 277 RBCZ = RBCZ/CRSUN 278 ELSE 279 TAUCZ = 0.0D0 280 RBCZ = 0.0 281 ENDIF 282C ENSURE THAT TAUCZ WAS NOT ACCIDENTALLY CALCULATED 283C DEEP IN THE STELLAR INTERIOR. IF YES, REDO CALCULATION. 284 IF (TAUCZ.GT.1.0E20) THEN 285 IBEG = IMAX + 1 286 GOTO 50 287 ENDIF 288 RETURN 289 END

1"""Python version of tauintnew.f (YREC 7be4133): the convective overturn 2time by G. Somers' method (3/17).""" 3 4import math 5 6import numpy as np 7 8from pyrec import constants_mod 9from pyrec import ovrtrn_mod 10from pyrec.kspline import kspline 11from pyrec.ksplint import ksplint 12 13# TAUINTNEW has SAVE, but since 7be4133 every local is set in each call 14# before it is read (JCORE = 0 at the start), so none is kept here. 15 16 17def _psca(hp: np.ndarray, hd: np.ndarray, hg: np.ndarray, k: int) -> float: 18 # The pressure scale height P/(rho g) at shell k. 19 cln = constants_mod.cln 20 return math.exp(cln * (hp[k] - hd[k])) / hg[k] 21 22 23def tauintnew(hs1: np.ndarray, lcz: np.ndarray, hr: np.ndarray, hp: np.ndarray, hd: np.ndarray, hg: np.ndarray, m: int, 24 mint: int, svel: np.ndarray, del1: np.ndarray, del2: np.ndarray, rbcz: float) -> float: 25 """Compute the convective overturn time TAUCZ (in ovrtrn_mod) and the 26 radius of the base of the convection zone: find the zone from the 27 centre up and from the surface down; for a single zone, H_P/v one 28 pressure scale height above its base; otherwise the lowest cell's width 29 over the velocity at its middle. 30 31 Translates: tauintnew.f#TAUINTNEW (7be4133) 32 33 Arguments: the profile of m points (mint of them interior shells), index 34 0 unused: hs1 (mass, g), lcz (convective or not), hr, hp, hd (log r, P, 35 rho), hg (gravity), svel (convective velocity), del1, del2 (the two 36 gradients compared to find the zone's edges); rbcz is output only. 37 38 Returns rbcz, the base's radius in solar radii (0 for a radiative 39 envelope). 40 41 """ 42 43 cln = constants_mod.cln 44 xdel = np.zeros(5) 45 xrad = np.zeros(5) 46 yrad = np.zeros(5) 47 yvel = np.zeros(5) 48 ypsc = np.zeros(5) 49 yspl = np.zeros(5) 50 51 lallcz = False 52 ibeg = 1 53 jcore = 0 # no core to start with (MHP 10/30/26) 54 while True: # label 50 55 # Find the top of the main internal convection zone (Fortran lines 56 # 74–99). 57 # Question (see the page): jcore keeps its value from an earlier 58 # pass when a restart begins in a radiative shell. 59 lradzone = not lcz[ibeg] 60 found = False 61 for i in range(ibeg, m + 1): # loop 60 62 if not lradzone and not lcz[i]: 63 jcore = i - 1 64 lradzone = True 65 if lradzone and lcz[i]: 66 found = True 67 break # to label 70 68 if not found: 69 # Fully convective, or a radiative envelope (maybe with a 70 # convective core). 71 if lradzone: 72 lallcz = False 73 lczsur = False 74 imax = m 75 jmax = m 76 jcore = m 77 to_100 = True 78 else: 79 lallcz = True 80 lczsur = True 81 i = 1 82 j = 0 83 jcore = 1 84 to_100 = False 85 else: 86 # label 70: find the base of the uppermost convection zone 87 # (Fortran lines 100–106). 88 # LCZ(J+1) traps the case where the last points are already 89 # radiative (JvS 02/2026). 90 j = 0 # the Fortran's J after the loop runs out 91 for jj in range(m - 1, 0, -1): # loop 80 92 if lcz[jj + 1] and not lcz[jj]: 93 j = jj 94 break # to label 90 95 to_100 = False 96 97 if not to_100: 98 # label 90 (Fortran lines 107–109). 99 imax = i 100 jmax = j + 1 101 ijmax = jmax - imax 102 103 # Skip convective shells just above the convective core. 104 # The check on imax and jcore (KC 2025-05-30) (Fortran lines 105 # 115–120). 106 # Question (see the page): a fully convective star. 107 if imax > 0 and jcore > 0: 108 if m != mint and hp[jcore] - hp[imax] < 1: 109 ibeg = imax + 1 110 continue # GOTO 50 111 112 # At least 3 shells thick, or search again from above (Fortran 113 # lines 124–131). 114 if not (lcz[imax] and lcz[imax + 1] and lcz[imax + 2]): 115 if not lcz[imax + 1]: 116 ibeg = imax + 1 117 else: 118 ibeg = imax + 2 119 continue # GOTO 50 120 121 # Interior model only: use jmax (Fortran lines 133–136). 122 if ijmax != 0 and m == mint: 123 imax = jmax 124 ijmax = 0 125 126 # Masses at the top and bottom of the zone, only for LCZSUR, 127 # which is true for any zone (Fortran lines 139–150). 128 hstop = hs1[m] 129 if jmax > 1: 130 hsbot = 0.5 * (hs1[jmax] + hs1[jmax - 1]) 131 else: 132 hsbot = 0.0 133 lczsur = (hstop - hsbot) / constants_mod.cmsun > 0.0 134 135 # label 100: taucz from the top down; 0 for a radiative envelope. 136 # The base's radius (Fortran lines 155–168). 137 if lczsur: 138 if lallcz: 139 rbcz = 1.0 # exp(cln * 1) below: 10 cm 140 else: 141 # The base where del1 - del2 changes sign, linear in log r. 142 dd2 = del1[imax - 1] - del2[imax - 1] 143 dd1 = del1[imax] - del2[imax] 144 fx = dd2 / (dd2 - dd1) 145 rbcz = hr[imax - 1] + fx * (hr[imax] - hr[imax - 1]) 146 rbcz = math.exp(cln * rbcz) 147 148 # One zone: H_P / v where r - r_b = H_P (Fortran lines 171–230). 149 to_110 = ijmax != 0 # several cells: the averaging method 150 if not to_110: 151 psca2 = _psca(hp, hd, hg, imax) 152 rtest2 = math.exp(cln * hr[imax]) - rbcz 153 valtest = rtest2 - psca2 154 if valtest > 0.0: 155 # The innermost point's scale height reaches below. 156 k = imax 157 cvel = svel[imax] 158 psca = psca2 159 ovrtrn_mod.taucz = psca / cvel 160 # Suspected bug (see the page): with the ENDIF where it 161 # now is, this branch falls through to label 110, which 162 # overwrites taucz (the ENDIF moved, KC 2025-05-31). 163 to_110 = True 164 else: 165 to_110 = True 166 for k in range(imax + 1, m): 167 psca1 = psca2 # noqa: F841 (unused in the Fortran too) 168 rtest1 = rtest2 # noqa: F841 (unused in the Fortran too) 169 psca2 = _psca(hp, hd, hg, k) 170 rtest2 = math.exp(cln * hr[k]) - rbcz 171 valtest = rtest2 - psca2 172 # Where r - r_b = H_P. 173 if valtest > 0.0: 174 kk = k 175 if k == imax + 1: 176 kk = k + 1 177 # Keep the stencil inside 1..m (MHP 10/30/26). 178 kk = max(3, min(kk, m - 1)) 179 for n in range(1, 5): 180 kn = kk - 3 + n 181 xdel[n] = (math.exp(cln * hr[kn]) - rbcz 182 - _psca(hp, hd, hg, kn)) 183 for n in range(1, 5): 184 yvel[n] = svel[kk - 3 + n] 185 kspline(xdel, yvel, yspl) 186 cvel = ksplint(xdel, yvel, yspl, 0.0, 0.0) 187 for n in range(1, 5): 188 ypsc[n] = _psca(hp, hd, hg, kk - 3 + n) 189 kspline(xdel, ypsc, yspl) 190 psca = ksplint(xdel, ypsc, yspl, 0.0, 0.0) 191 for n in range(1, 5): 192 yrad[n] = hr[kk - 3 + n] 193 kspline(xdel, yrad, yspl) 194 envrcl = ksplint(xdel, yrad, yspl, 0.0, 0.0) # noqa: F841 (unused in the Fortran too) 195 ovrtrn_mod.taucz = psca / cvel 196 to_110 = False 197 break # to label 140 198 199 if to_110: 200 # label 110: r_b + H_P is outside the star, or several cells 201 # lie in the radiative envelope; use the lowest cell's 202 # middle. First find its upper edge (Fortran lines 231–272). 203 k = m # the Fortran's K after the loop runs out 204 for kk2 in range(imax, m): 205 if not lcz[kk2]: 206 k = kk2 207 break # to label 120 208 itop = k - 1 # label 120 209 dd2 = del1[itop] - del2[itop] 210 dd1 = del1[itop + 1] - del2[itop + 1] 211 fx = dd2 / (dd2 - dd1) 212 rtcz = hr[itop] + fx * (hr[itop + 1] - hr[itop]) 213 rtcz = math.exp(cln * rtcz) 214 delr = rtcz - rbcz 215 # The shells around the cell's middle. 216 rtest = math.log10(rbcz + 0.5 * delr) 217 k = itop + 2 # the Fortran's K after the loop runs out 218 for kk2 in range(imax, itop + 2): 219 if hr[kk2] > rtest: 220 k = kk2 221 break # to label 130 222 rtest = math.exp(cln * rtest) # label 130 223 kk = k 224 if k == imax + 1: 225 kk = k + 1 226 kk = max(3, min(kk, m - 1)) # MHP 10/30/26 227 for n in range(1, 5): 228 xrad[n] = math.exp(cln * hr[kk - 3 + n]) 229 yvel[n] = svel[kk - 3 + n] 230 kspline(xrad, yvel, yspl) 231 cvel = ksplint(xrad, yvel, yspl, rtest, 0.0) 232 ovrtrn_mod.taucz = delr / cvel 233 234 # label 140: the base's radius in solar units (Fortran line 277). 235 rbcz = rbcz / constants_mod.crsun 236 else: 237 ovrtrn_mod.taucz = 0.0 238 rbcz = 0.0 239 240 # A taucz found deep in the interior: search again above (Fortran 241 # lines 284–287). 242 # Question (see the page): restarts without a limit. 243 if ovrtrn_mod.taucz > 1.0e20: 244 ibeg = imax + 1 245 continue # GOTO 50 246 return rbcz

tauintnew.f

Purpose

How it fits

The method

Suspected bugs, checked by running the code

1. (not reached) The overturn time at IMAX is overwritten.

2. (never occurred) The restarts can read past the profile, and can loop forever.

Definitions

subroutineTAUINTNEW

Change log

1. (not reached) The overturn time at `IMAX` is overwritten.