seculr.f

The secular transport of angular momentum, and of the composition it carries, by rotational instabilities, over one time step. Its comment: "SECULR CALCULATES THE TRANSFER OF ANGULAR MOMENTUM (AND THE ASSOCIATED TRANSPORT OF COMPOSITION) DUE TO SECULAR ROTATIONAL INSTABILITIES." The mechanisms considered are "secular and dynamical shear, Eddington circulation, and the GSF instability." It is the central routine of rotational mixing: it sets up, iterates and checks the diffusion of $\Omega$ and of the abundances.

• Who calls it: GETW, once per time step, for a rotating model.
• Inputs (arguments): the step DT; the model (HD, HG, HI, HL, HP, HR, HS, HS1, HS2, HT, HCOMP, M); the rotation (OMEGA, specific angular momentum HJM and its saved copy HJMSAV, ETA2, QIW, R0); the composition at the start of the step, HCOMPS (filled by MIDMOD; GETW owns the array); the convection flags LCZ and the range of radiative interfaces IMIN–IMAX; and for the wind, BL, HSTOT, SMASS, TEFFL, HSBOT, HSTOP, LCZSUR (the star has a surface convection zone).
• Outputs: OMEGA, HJM, HCOMP updated; HV (the circulation velocities), COD2 (composition diffusion coefficients), HICZ, OMEGAS; the flags LOK (converged) and LREDO, IREDO (the step must be cut, set by the checks).
• What it calls, in order: SETUPV (factors for the velocities), GETQUA, VCIRC (stability tests and velocities; sets LDO), GETFC (with LVFC), MIXCZ, the wind routines MWIND and MCOWIND, SOLID (solid-body rotation), CODIFF (diffusion coefficients), ROTGRID (the equally spaced grid), DCOEFT and TRIDIA (the diffusion equation for $\Omega$ and its tridiagonal solution) or DADCOEFT (diffusion plus advection, LDIFAD), EQ2MOD (back to the model's grid), NDIFCOM (composition diffusion), CHECKJ and CHECKC (convergence and step-cutting checks).
• Modules, main ones: NMLPAR_MOD (namelist inputs: ITDIF2, the maximum number of iterations; LDISK, TDISK, disk locking; LSOLID, IMPJMOD; LVFC, LDIFAD; LJDOT0, the wind; LCODM, CODM, the constant diffusion coefficient; LNOJ; LEXCOM), DISK_MOD (SAGE, the age), INTVAR_MOD (RM, the mid-shell radii), OLDROT2_MOD (WMST), EGRID_MOD (NTOT), CONSTANTS_MOD, IOUNITS_MOD (IMODPT) and DIMENSIONS_MOD (JSON).
• Arrays it owns: FACI (factors for the velocities; formerly /INTFAC/), filled by SETUPV and read by VCIRC; VFC (formerly /VARFC/), cleared here at each call, filled by GETFC and read by CODIFF.

1. Disk locking: if LDISK and the age is below TDISK, the surface rate is held fixed (LOCK).

2. Shear. At each interface the logarithmic gradient of $\Omega$ is centered on the mid-shell radius $r_m$: $$\frac{d\Omega}{d\ln r}\bigg|_i \approx \frac{\Omega_i - \Omega_{i-1}}{4}\left[\frac{1}{\ln r_i - \ln r_m} + \frac{1}{\ln r_m - \ln r_{i-1}}\right],$$ which is $(\Omega_i - \Omega_{i-1})/(\ln r_i - \ln r_{i-1})$ when $r_m$ is the geometric mean. Here $r_m$ = RM is the arithmetic mean $\tfrac12(r_i + r_{i-1})$ (SETUPV), and since $\tfrac14(1/x + 1/y) \ge 1/(x+y)$, the gradient comes out larger than the simple difference, by a fraction of about $(\ln r_i - \ln r_{i-1})^2/16$.

3. Iteration IT = 1 … ITDIF2. From the second iteration on, $\Omega$ is replaced by the average of the start and the latest solution, $\Omega \leftarrow \tfrac12(\Omega + \Omega_\mathrm{start})$, to get, in the comment's words, "a corrected V and thus a more accurate run of diffusion coefficients"; the velocities are recomputed, $j$ and the first four abundances are reset to their starting values, and convection zones are remixed (MIXCZ; the flag LIMIX set before it, line 253↗, is no longer passed to it: unused?).

4. Wind. For an isolated surface convection zone (no diffusion across its base, HV(IMAX) = 0) the wind torque is applied to that zone (MWIND); for solid-body models (LSOLID or IMPJMOD = 1) to the whole star, followed by SOLID. Otherwise the wind enters the diffusion equation as a boundary term (MCOWIND, giving WIND1, WIND2), or not at all when locked.

5. Unstable regions. If VCIRC found no instability (LDO false), the routine ends. Otherwise CODIFF gives the coefficients, and each run of interfaces with $v > 0$ (HV) is solved on its own (or, with LCODM and CODM > 0, the whole domain IMIN-1…IMAX as one region, "MHP 8/13," lines 308–314↗; that this also runs for solid-body models, where the solution for $\Omega$, the only place CODM acts, is skipped, is KNOWN_ISSUES.md's issue 1, open, with a proposed gate): ROTGRID maps it to an equally spaced grid, the diffusion equation for $\Omega$ is solved implicitly (DCOEFT + TRIDIA, or DADCOEFT; the equation, from DCOEFT's header, is $\frac{\partial\Omega}{\partial t} = \frac{1}{4\pi\rho r^2}\frac{1}{I/M}\frac{\partial}{\partial r}\bigl(D\frac{\partial\Omega}{\partial r}\bigr)$, "the diffusion equation we are solving," $D$ being CODIFF's $D_J = 4\pi f_w v r$ multiplied in ROTGRID by "rho r^4 by d chi/dr" and by "I/MR^2" (HI/HS2, the shell's mass), for its equally spaced coordinate $\chi$; the composition obeys a diffusion equation of the same kind with $D_X$), EQ2MOD maps the change in $j$ back, and NDIFCOM diffuses species 1–4 ($X$, $Y$, $Z$, He3: those that change the mean molecular weight; the comment names only H, He3, and He4, lines 413–415↗, but the call covers $Z$ too, lines 419–424↗). A region that ROTGRID flags (LDUM2: a single radiative interface between two convection zones) is skipped. For solid-body models (LSOLID or IMPJMOD = 1) the solution for $\Omega$ is skipped ("MHP 3/09"), but NDIFCOM still diffuses the composition. The center and the surface shells are treated as convective during this.

6. Checks. CHECKJ checks the solution and updates $\Omega$; CHECKC checks the composition. On convergence (LOK), or after ITDIF2 iterations, species 5–11 (or 5–15 with LEXCOM) are diffused with NDIFCOM. This call uses neither the last region's IBEG, IEND nor its grid: With LOK true, as it always is here (line 448↗, line 452↗), NDIFCOM finds every unstable region itself and builds each one's grid (ndifcom.f, lines 64–121↗) (settled from the code, Nov 9, 2026).

Printing is switched off (LPRT = .FALSE., "G Somers 11/14"), so the block that writes the composition changes to IMODPT never runs. Large commented-out blocks remain (the old EQUAL, EQUAL2, DIFCOM, BSCOMP routines, and a coupling-factor diagnostic).

Physics Options (Angular momentum transport): LINSTB includes rotationally induced mixing and hydrodynamic angular-momentum transport; LSOLID enforces solid-body rotation everywhere (for strongly coupled magnetic models); WALPCZ makes $\Omega \propto r^{\mathrm{WALPCZ}}$ in the convection zone (0 solid body, −2 uniform specific angular momentum); IMPJMOD couples globally: 0 no effect, 1 solid body everywhere, 2 solid-body core with a decoupled convection zone, 3 solid-body core coupled at the zone's base. ITDIF2 is the number of iterations between nuclear burning and rotational mixing (1 = none), DJOK the tolerance of rotational mixing.