checkc.f

This file checks the composition after a step of rotational mixing for anomalous abundances (below 0 or above 1), cuts the diffusion time step if it finds any, and recomputes the mean molecular weight, which the mixing changes.

• Inputs: the arguments HCOMP (mass fractions), IT (the iteration on the diffusion coefficients), LPRT, M; ITDIF2 (the last iteration, for the print) and LEXCOM from NMLPAR_MOD; the envelope abundances XENV, ZENV, AMUENV from COMP_MOD and YENV, Y3ENV from COMP2_MOD; the start-of-step composition HCOMPP from OLDMOD_MOD.
• Who calls it: SECULR, the secular (rotational) mixing step.
• Outputs: DT (halved on an error), IREDO, LOK, LREDO; small anomalies in HCOMP are repaired in place; the mean molecular weight AMUM in MDPHY_MOD, read by VCIRC, SETUPV (the $\mu$ gradient), GETW and MIDMOD.

1. Species checked: on every iteration only $X$, $Y$ and $^3$He (1, 2, 4); species 5 to 11 (5 to 15 with LEXCOM) are mixed by SECULR once, after its iteration loop, and the code's comment names the clip to $[0,1]$ at the start of the next step in MIX and ROTMIX as their safeguard. $Z$ (3) is skipped, since it "is not diffused as a unit." Only interior points $2..M-1$ are checked.

2. A value outside $[0, 1]$ that is tiny in size, $|X_{j,i}| < 10^{-5} X_{j,M}$ (relative to the surface value), is replaced by $\max(0, \tfrac12(X_{j,i+1} + X_{j,i-1}))$ (lines 109–117↗): the comment explains that mixing into a burning region can give small negative abundances. The branch inside it that sets the centre or surface value to 0 (I.EQ.1.OR.I.EQ.M, lines 110–111↗) is never taken, since the loop runs over $2..M-1$ (line 103↗): unused? Otherwise the time step is halved (LREDO), up to 3 times; then the run stops with exit status 1. On the first anomaly that is not repaired, the routine jumps to its end, so $\mu$ is not recomputed in that call.

3. At the last iteration, with LPRT, print the largest fractional change of an abundance over the step, and the surface lithium if LEXCOM.

4. For IT > 1, the mean molecular weight of a fully ionized gas: per ion, from the surface value and the local differences; per electron, directly: $$\frac{1}{\mu_\mathrm{ion}} = \frac{1}{\mu_\mathrm{ion,env}} + \sum_k \frac{X_k - X_{k,\mathrm{env}}}{A_k}, \qquad \frac{1}{\mu_e} = \frac{X}{A_\mathrm{H}} + 2\left(\frac{Y}{A_4} + \frac{X_3}{A_3}\right) + \frac{Z}{2},$$ for $k$ = H, $^4$He, $Z$ (with $A = 12$), $^3$He; then $\mu = \mu_\mathrm{ion}\mu_e/(\mu_\mathrm{ion} + \mu_e)$, that is $1/\mu = 1/\mu_\mathrm{ion} + 1/\mu_e$. The comment notes that partial ionization would need a correction in massive stars.