PyREC scaffolding (from YREC development, 7be4133)/ sconvec.fGlossary

sconvec.f

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

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

Variables on this page (10)

c4pi, cln, clsun, ishort, json, ldh, xxdh, yydh, zdh, zzdh

Purpose

Tests the edges of the first convection zone (every zone at rp_1) for semiconvection and, where the region next to it is unstable, extends the zone, following Castellani et al. (1971, Ap&SS 10, 340), as the header comment says. It then merges zones that have come to overlap.

How it fits

  • Who calls it: MIX, before the composition is mixed in the convection zones, when semiconvection is on (LSEMIC) and from the second iteration (ITLVL > 1).
  • Inputs: the step DELTS (s), the model's $\log\rho$ (HD), $L$ (HL, in solar units, by the factor CLSUN), $\log P$, $\log r$, $\log m$, $\log T$ (HP, HR, HS, HT), the composition HCOMP, TEFFL, the mixing length CMIXL (passed on to TPGRAD), and the zones MXZONE(12,2) (bottom and top shell of each), NZONE.
  • What it calls: for each test point, EQSTAT (equation of state: $\nabla_\mathrm{ad}$, the mean molecular weights), GETOPAC (opacity) and TPGRAD (the radiative gradient $\nabla_\mathrm{rad}$).
  • Outputs: MXZONE and NZONE, changed; a line in the short output (ISHORT) for each edge it moves and each merger. Also, HD at the shells it tests is overwritten with the density returned by EQSTAT.
  • Modules: DEBHU_MOD (it loads XXDH, YYDH, ZZDH, ZDH(1..3) with each test point's composition for the Debye–Hückel correction in EQSTAT, when LDH is set), CONSTANTS_MOD, NMLPAR_MOD, IOUNITS_MOD (ISHORT), DIMENSIONS_MOD.

The method

For the first zone only (a local switch, L_ONLY_CORE, set to true; with it false, every zone) and each edge that is not the center or the surface (LFLAGD, LFLAGU); a fully convective star is skipped:

  1. At the first radiative shell outside the edge (JMR), compute $\nabla_\mathrm{rad}$ with the zone's composition. If $\nabla_\mathrm{rad} < \nabla_\mathrm{ad}$, the edge is stable. Otherwise keep $\nabla_\mathrm{rad,1}$, the mean molecular weight $\mu_\mathrm{CZ}$ (XINTMU = AMU + EMU), and $1 - \nabla_\mathrm{ad}/\nabla_\mathrm{rad}$.
  2. Recompute at the same shell with its own composition: $\nabla_\mathrm{rad,0}$, and the ratio $f = \nabla_\mathrm{rad,1}/\nabla_\mathrm{rad,0}$ (FDEL).

  3. The farthest the instability can reach in one step is the propagation speed (eq. 5′, p. 347 of the paper) times the step: $$R_\mathrm{max} = \Bigl(1 - \frac{\nabla_\mathrm{ad}}{\nabla_\mathrm{rad,1}}\Bigr)\,\frac{L}{10\cdot4\pi P\,r^2}\,\Delta t,$$ with $L$, $P$, $r$ at the zone's edge shell.

    "Overshoot" in the comments means semiconvection. The comments at line 198↗, line 205↗, and line 219↗ speak of "OVERSHOOT," but the routine is the semiconvection of Castellani et al. (1971) (lines 63–65↗), and its own merger message says "MERGED DUE TO SEMICONVECTION" (line 282↗), so the word is a leftover. 5. Mergers. If the top of zone $j$ is now above the bottom of zone $j+1$, they become one zone, and the list is shifted down (but see bug 1).

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) A merger keeps the upper zone's top.

(lines 278–291↗) When zones $j$ and $j+1$ overlap, the merged zone takes the bottom of zone $j$ and the top of zone $j+1$ (MXZONE(J+1,1) = MXZONE(J,1), then the list is shifted). The walk upward from zone 1's top (step 4) passes through any zone above it and can end beyond that zone's top; the merged zone then ends at the upper zone's top, below the edge just found, and the semiconvective region above it is not mixed. Only with LSEMIC (off by default) and a second zone within $R_\mathrm{max}$ of the first.

A possible fix:

sconvec.f

@@ lines 282–288 @@
    93    FORMAT(2X,'MIXED ZONES MERGED DUE TO SEMICONVECTION'
      *          /2X,'OLD',2('[',I3,'-',I3,']'),
      *          ' NEW','[',I3,'-',I3,']')
+C 11/09/26 KC via Claude
+C code changed to keep the higher of the two tops when two zones merge.
          MXZONE(J+1,1) = MXZONE(J,1)
+         MXZONE(J+1,2) = MAX(MXZONE(J,2),MXZONE(J+1,2))
          DO 90 K = J,NZONE-1
             DO 95 J2 = 1,2
                MXZONE(K,J2) = MXZONE(K+1,J2)

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): In the one run with semiconvection, sample_cases/UpperMS/Sample_m9p0feh+0p0GS98_yaleeos_sc_TAMS, no two zones merged. It needs LSEMIC (which parmin.f refuses with any overshoot) and a second zone within the walk's reach; YREC's other semiconvective case, sample_cases/UpperMS/Sample_m9p0feh+0p0GS98_yaleeos_sc_TAHB (core helium burning), was not run.

2. (harmless in the runs) The radius is not updated in the walk.

(line 144↗; lines 227–251↗) In the loop over shells II, every quantity passed to TPGRAD is reset for the shell, except RL, which stays at HR(JMR). TPGRAD's $\nabla_\mathrm{rad}$ does not use the radius (tpgrad.f, lines 80–81↗), so the stability test is not affected today; only the gravity $G m/r^2$ that TPGRAD computes for the convective case is, and its results are not used here. A change that made $\nabla_\mathrm{rad}$ depend on $r$ would break the test.

A possible fix:

sconvec.f

@@ lines 225–231 @@
                IF(RSUM.GT.RMAX) GO TO 33
 C TEST ON RESCALED RAD. GRADIENT FOR CONVECTION
                B = HL(II)
+C 11/09/26 KC via Claude
+C code changed to set RL for each shell of the walk, as the other
+C structure variables are.
                SL = HS(II)
+               RL = HR(II)
                PL = HP(II)
                TL = HT(II)
                DL = HD(II)

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): The walk ran once in the 34 runs, at shell 380 of sample_cases/UpperMS/Sample_m9p0feh+0p0GS98_yaleeos_sc_TAMS, the only run with semiconvection, and there RL was not that shell's radius. With the fix, the track and .store file are identical (706 models, 0.0253726 Gyr), where a null change (nearly every opacity, all with $3 \le \log T \le 9$ and $-6 \le \log\rho \le 9.75$, times $1 + 10^{-10}$) moves $\log L$ by up to $7.56\times10^{-6}$ dex, as expected, since $\nabla_\mathrm{rad}$ does not use the radius. Fragile code.

Definitions

subroutineSCONVEC

lines 4–302 · GitHub

Uses modules CONSTANTS_MOD, DEBHU_MOD, DIMENSIONS_MOD, IOUNITS_MOD, NMLPAR_MOD

Called by MIX

Flowchart of its calls, in order (9)
for I = 1,ILOOP (DO 210), lines 98–273
for K = 1,2 (DO 200), lines 121–271
EQSTAT159→GETOPAC163→TPGRAD164→EQSTAT186→GETOPAC190→TPGRAD191
for II = JBEG,JEND,JSTEP (DO 32), lines 218–255
EQSTAT243→GETOPAC247→TPGRAD248

Change log

The changes since rp_1

Changed since rp_1. Only the first zone in MXZONE is tested (L_ONLY_CORE; see the question under the method); the mixing length CMIXL is a new argument, passed to TPGRAD; the unused /NWLAOL/ and /NEWOPAC/ are dropped; COMMON blocks to modules, IMPLICIT NONE.