microdiff.f

MICRODIFF is the driver of YREC's newer treatment of gravitational settling (microscopic diffusion) of helium, metals and light elements. Its header (Garrett Somers, 05/2015) says it is adapted from GRSETT, "which is retained for backwards compatability," and "follows the formulism of Bahcall and Loeb 1989." It moves the model onto an equally spaced grid in radius, diffuses each species there, and maps the changes back onto the model's mass grid.

• Inputs (arguments): the time step DT (s), the composition HCOMP(15,M), $d\ln P/dr$ (HQPR), $\log r$ (HR, cm), $\log\rho$ (HD), the unlogged mass HS1 (g), $\log T$ (HT), the convection flags LC, the number of points M and the total mass STOT. STOT is in grams; MICRODIFF_SETUP converts copies of DT, STOT, HS1, and HQPR (DTU, STOTU, HS1U, and HQPRU, which MICRODIFF hands to all four parts of the method) in place to Bahcall and Loeb's units (radius in $R_\odot$ = 6.9598×10¹⁰ cm, temperature in $10^7$ K, mass in $100\ \mathrm{g\,cm^{-3}}\times R_\odot^3$, time in $2.7\times10^{13}$ yr, a unit that "includes factor of 2.2 from ln lambda"; the same as SETUP_GRSETT), and MICRODIFF_ETM converts them back; the caller's values are not changed. The temperature gradient $\nabla = d\ln T/d\ln P$ (SDEL(2,·) in SCRTCH_MOD), named in the header, is read not here but by MICRODIFF_MTE.
• Who calls it: MIX (when the settling flags are on and LNEWDIF selects the new routine) and ROTMIX (the same inside the rotational-mixing step).
• What it calls: the four parts of the method, in order: MICRODIFF_SETUP (units and boundaries), MICRODIFF_MTE (model to equal grid), MICRODIFF_RUN (once per diffused species) and MICRODIFF_ETM (equal grid back to model).
• Output: the updated HCOMP, in place.
• Modules: NMLPAR_MOD, for the namelist switches LDIFY, LDIFZ and LDIFLI (read by PARMIN and shared with GRSETT, MIX, ROTMIX and the other MICRODIFF routines), and DIMENSIONS_MOD for JSON. The other members of rp_1's /GRAVST/, /GRAVS3/ and /GRAVS4/ (GRTOL, NITER_GS, FGRY, FGRZ, LTHOUL, LNEWDIF) are in NMLPAR_MOD too, but not used here.

After the set-up (which returns at once for a fully convective model, or one with no hydrogen or helium left), up to three kinds of species are diffused, each by a call to MICRODIFF_RUN that returns the change in abundance on the equal grid:

1. Hydrogen if LDIFY: column $j=1$ of the three-species vector ESPEC (X, Y, and a third species), with the metals as the third species, taken as iron ($A = 55.86$, $Z = 26$). Helium then follows from $Y = 1 - X - Z - X_4$ in MICRODIFF_ETM. LDIFY is always on when MICRODIFF runs: MIX and ROTMIX call it only under IF(LDIFY), and PARMIN turns LDIFY on with LDIFZ or LDIFLI; so EDELX is always set before MICRODIFF_ETM applies it (settled from the code, Nov 9, 2026: mix.f, line 412↗, rotmix.f, line 151↗, parmin.f, lines 1239–1252↗; lines 157–182↗).

2. Metals if LDIFZ: column $j=3$, again as iron.

3. Light elements if LDIFLI: NLIGHT = 3 species in HCOMP rows CSLIGHT = 13, 14, 15, with atomic weights 6.015, 7.016, 9.012 and charges 3, 3, 4 (⁶Li, ⁷Li and ⁹Be: STARIN fills these rows from XLI6_INI, XLI7_INI and XBE9_INI), each in turn as the third species. The last point is set as ESPEC(3,I) = ELI(II,NPT) with I = NPT, the value the loop variable has after DO I=1,NPT-1 completes (the Fortran standard), so it is correct as written (lines 220–224↗).

The changes EDELX, EDELZ and EDELLI are passed to MICRODIFF_ETM, which interpolates them back and applies them.

The physics included. Gravitational settling, thermal diffusion (the $A_T\nabla$ term) and concentration diffusion, from Thoul et al.'s solution of Burgers' equations (THDIFF, called by MICRODIFF_COD), with every species fully ionized (fixed charges) and no radiative accelerations (levitation). The results are therefore least trustworthy where radiative accelerations or partial ionization matter, as in the thin envelopes of F and hotter stars.

Grids. Variables prefixed E are on the equally spaced grid in radius; those also suffixed _H are at the midpoints between its points. Both are needed by the two-step Lax–Wendroff scheme.