mix.f

MIX updates the composition of a model over one time step: nuclear burning, mixing of convection zones, semiconvection, gravitational settling and deuterium burning. Its header gives the outline: (1) find the convection zones and the mixed regions with overshoot; (2) compute the burning rates for the given composition, density and temperature; (3) solve implicitly for the abundances in radiative regions; (4) solve for burning in each convection zone; (5) semiconvection if desired; (6) microscopic diffusion if desired. The header also says that MIX is called "prior to rezoning and on the third level of iteration in the Henyey portion." So it runs as a predictor and a corrector: the main program saves the start-of-step composition in HCOMPP and calls MIX with ITLVL = 1 on the old structure, before rezoning, to predict the composition (with the explicit equilibrium burning of step 4); then, in every Henyey iteration from level 3 on, CRRECT restores HCOMPP and calls MIX again with the current $T$ and $\rho$. The level-1 result is only a starting guess.

• Inputs: the model (HD, HL, HP, HR, HS = $\log m$, HS1 = $m$, HS2 = shell masses, HSTOT, HT, LC, M), the composition HCOMP(15,M), the step DELTS (s), the iteration level ITLVL, TEFFL and the mixing length CMIXL, passed on to SCONVEC. Parameters come from NMLPAR_MOD: DPENV, LSEMIC, the settling switches and limits (LDIFY, XMIN, YMIN, DT_GS, LNEWDIF), TCUT, the accretion settings (LMDOT, DMDT0, COMPACC) and LEXCOM, all set by PARMIN.
• Who calls it: CRRECT and the main program in main.f.
• What it calls: CONVEC (zones), RATES and DEUTRATE (rates), KEMCOM (implicit burning), EQBURN (explicit burning at ITLVL = 1), SCONVEC (semiconvection), MICRODIFF or GRSETT (settling), and DBURN (deuterium).
• Outputs: HCOMP in place; the time step in years DDAGE; the zone edges JCORE, JENV, MXZONE, MXZON0. The deuterium rate goes to DRATE in DEUTER_MOD. The mixed zones are printed to ISHORT.
• Modules: besides those above, CONSTANTS_MOD (CLN, CGL, CSECYR), IOUNITS_MOD (ISHORT) and DIMENSIONS_MOD (JSON).

1. Deep mixing. If DPENV < 1 (and ITLVL > 1), every shell with $m \le \mathrm{DPENV}\,M_\ast$ is flagged convective in a copy of the flags (LCC) before CONVEC finds the zones; a commented-out variant used a fraction of the peak luminosity instead (MHP 1/95). Otherwise LCC = LC.

   DPENV mixes from the center. The first comment says "DPENV MIXES THE STAR FROM THE SURFACE TO MASS FRACTION DPENV, E.G. DPENV = 0.7 MEANS THE OUTER .3 OF THE STAR IS MIXED" (lines 127–128↗), but the code flags the shells inside $\mathrm{DPENV}\,M_\ast$ (lines 129–141↗), as the next comment ("MIX FROM CENTER TO A FIXED MASS FRACTION," line 130↗) and the documentation ("Full mixing from center to fractional mass depth," docs/source/namelist_physics.md, line 296) say, so the first comment is stale. 2. Rates. For every shell hotter than TCUT(1), RATES gives the rates HR1–HR13, HF1, HF2, and DEUTRATE the deuterium rate if the surface (or accreted) deuterium exceeds $10^{-11}$. Cooler shells get zero rates. 3. Implicit burning. KEMCOM solves shell by shell in radiative zones (down to the TCUT(1) cutoff) and zone by zone in convective ones, which also homogenizes them. 4. Explicit burning (only ITLVL = 1): EQBURN gives hydrogen and helium burning rates assuming equilibrium ³He and CN abundances, which overwrite $X$, $Z$, ¹²C and ¹⁶O: $X \to X + \dot X\,\Delta t$, $Z \to Z - \dot Y\,\Delta t$, with $\Delta t$ = DDAGE$\times10^{-9}$ (Gyr). Here $\dot Y < 0$ (DYDT = -4.0D0*(R412 + 3.0D0*R3A), from the ¹²C + α and triple-α rates), so $Z$ rises: the helium burnt becomes carbon and oxygen, counted in $Z$. EQBURN starts from the start-of-step composition HCOMPP, not from KEMCOM's result, and overwrites only those four (and only where its rates are nonzero); what KEMCOM gave for ³He, the other CNO isotopes (¹³C, ¹⁴N, ¹⁵N, ¹⁷O, ¹⁸O) and the light elements is kept.

   DBURN's rate is per Gyr, like DT. The explicit and deuterium burning use DT = DDAGE*1.0D-9, in Gyr (line 293↗; lines 485–503↗), which matches EQBURN's rates, "per gigayear per amu" (KEMCOM's comment, kemcom.f, lines 85–88↗). DBURN's rate DRATE = X*RDEUT*C21 (deutrate.f, line 81↗) states no units, but RDEUT is per second, with a factor $N_A$/AMDEUT (deutrate.f, line 62↗), and C21 = 5.240358D-8 (deutrate.f, line 30↗) is $3.1558 \times 10^{16}$ s/Gyr divided by $N_A$, so DRATE is per Gyr per amu too.

MIX switches settling off for the rest of the run (it sets LDIFY = .FALSE. in NMLPAR_MOD) once the central hydrogen falls below XMIN. A model without a surface convection zone gets only a warning ("diffusion not meaningful") and no settling.

Physics Options: LSEMIC turns on semiconvection; DPENV (experimental) mixes fully "from center to fractional mass depth"; for settling, DT_GS (default 0.1) splits the step into several when the global step exceeds DT_GS times the settling time at the base of the convection zone, XMIN (1.0D-3) turns diffusion off once the central $X$ is below it, and YMIN (1.0D-3), in the documentation's words, "once surface helium Y < YMIN"; LNEWDIF (legacy) turns on the updated settling routines.