getw.f

GETW evolves the internal rotation over one model time step. In its header's words (MHP 11/91, replacing the old FINDW): it "enforces local conservation of angular momentum in radiative regions, and enforces the user-specified rotation law in convective regions"; it solves for the transport of angular momentum and the associated mixing if desired (LINSTB), accounts for angular momentum loss by a magnetic wind (LJDOT0), and, with instabilities, burns lithium and beryllium (LEXCOM). The method is that of Pinsonneault, Kawaler, Sofia and Demarque (1989, ApJ 338, 424), as the header cites. It is the central routine of YREC's rotation.

• Who calls it: only the main program (main.f), once per model step.
• Inputs: the model at the end of the step (HCOMP, HD, HL, HP, HR, HS, HS1, HS2, HT, LC, M, HSTOT, BL, TEFFL, SMASS), the form factors FP, FT, the step DELTS, the limit DWMAX, HJM ($j$ per shell), the mixing-length parameter CMIXL (passed to PHYSIC and GETTAU) and the print flag LPRT0 (set by the main program); both are arguments since rp_1. The model at the start of the step comes from the modules OLDMOD_MOD, OLDPHY_MOD, OLDROT_MOD, OLDROT2_MOD.
• What it calls (in order of use): OVROT (convective and overshoot zones for angular momentum: the flags LCZ), GETROT ($\Omega$ from $j$), MWIND (wind loss) and WCZIMP (rotation law in the surface convection zone after the loss), PHYSIC (gradients, viscosities, diffusivities), GETTAU (convective overturn time for the wind; its argument RBCZ, the radius of the convection zone's base, goes to ENVR, a local scalar never read afterwards, lines 243–244↗: unused?), MIDMOD (the structure interpolated in time within the step), MIXCZ (homogenize convection zones), SECULR (the secular instabilities, solved as diffusion equations), CONVEC, LIRATE88, LIBURN or LIBURN2 (light-element burning), and BURSMIX (Bulirsch–Stoer mixing plus burning).
• Outputs: HJM, OMEGA, HI, ETA2, QIW, R0 (new rotation), HCOMP (mixed and burned composition), LJDOT, JENV0; and through modules: HCOMPP in OLDMOD_MOD (beginning-of-step composition for the next step), the lithium rates in OLDRAT_MOD, PSCAL0/PSCAL1 in LIOV_MOD, FRACSTEP and TAUCZ in OVRTRN_MOD. HCOMPS, the composition at the start of the step (the former /OLDAB/), is now GETW's own array, filled by MIDMOD and passed to SECULR.
• Modules (COMMON blocks until rp_1, about thirty then, several of them unused): NMLPAR_MOD for the run's settings (LINSTB, LSOLID, IMPJMOD, LJDOT0, LDISK, TDISK, DTDIF, DTWIND, LEXCOM, and the overshoot's LOVSTE, ALPHAE), DISK_MOD (SAGE), OLDMOD_MOD, OLDROT_MOD, OLDROT2_MOD (previous model), MDPHY_MOD, TEMP_MOD, OLDPHY_MOD (physics at mid-step and at the start), NEWRAT_MOD, OLDRAT_MOD (lithium rates), QUADRU_MOD (quadrupole moment, from GETQUA), LIOV_MOD (overshoot), OVRTRN_MOD, DIMENSIONS_MOD, CONSTANTS_MOD.

1. OVROT marks the shells that rotate as convective (convection zones plus overshoot). Disk locking holds the surface rotation while the age is below TDISK.

2. Without instabilities (LINSTB false): $\Omega$ from $j$ (GETROT; with LSOLID or IMPJMOD = 1 the whole star is treated as convective, so it rotates as a solid body). If the wind is on and the surface is convective, MWIND removes angular momentum over the whole step from the surface convection zone (or the whole star), with FRACSTEP = 0.5 ("timestep average loss rate": the point in the step at which MWIND interpolates $\tau_\mathrm{cz}$, so 0.5 uses its mid-step value), and WCZIMP resets its rotation. Done.

3. With instabilities, the step is divided into $N$ diffusion substeps, $\Delta t = \Delta t_\mathrm{model}/N$, $N = \min(\lceil \mathrm{DWMAX}/\mathrm{DTDIF}\rceil, \lceil \mathrm{DTWIND}/\mathrm{DTDIF}\rceil)$. All three are fractional changes of $\Omega$, not times: DWMAX is the largest $2|\Omega-\Omega_\mathrm{old}|/(\Omega+\Omega_\mathrm{old})$ of the last step (WTIME), DTWIND the largest $\Delta\Omega/\Omega$ allowed per model step (WTIME's comment), DTDIF the largest per substep (REMAP sets them from TIME_DW_GLOBAL and TIME_DW_MIX); DWMAX = 0 is replaced by DTWIND. PHYSIC and GETTAU are called once. Then, for each substep:

4. MIDMOD interpolates the structure linearly in time to the substep's fraction $f$ = SUMDT/DELTS, and the time derivative of $\theta$ (Maeder & Zahn 1998, per the comment) is updated from $\Omega$: $\theta \propto \bar\Omega\,\Delta\Omega/\Delta r$. GETW stores, at each interface, THP $= (r^2/g)\,\bar\Omega\,d\Omega/dr$, with THM $= r^2/g$ from SETUPV, and with LDIFAD VCIRC turns its change over the substep into the $d\Theta/dt$ term (THN) of the meridional circulation. At the first substep THP is the previous model's value THO (lines 271–278↗); at later ones it is recomputed from OMEGA, which is not updated until the end of the step (bug 1);

5. for a fully convective star (LALLCZ), MWIND and WCZIMP remove angular momentum over the substep here; otherwise the wind's quantities (HICZ, HSBOT, HSTOP, OMEGAS, LCZSUR) go to SECULR, which calls MWIND itself;

6. unless the star is fully convective, MIXCZ homogenizes convection zones and SECULR solves the diffusion equations for $j$ and composition. If it asks for a shorter step (LREDO), CHECKJ or CHECKC has already halved DT; $j$, the composition and $\mu$ are restored, SUMDT goes back by $2\,$DT (the old substep), and the substep is redone;
7. the composition at the start of the next substep is stored (HCOMPP);
8. with the extended composition (LEXCOM): the depth of the convection zone (CONVEC), the overshoot depth $\ell = \alpha_e H_p$ with $H_p = Pr^2/(\rho Gm)$, the light-element rates (LIRATE88: its end-of-step rates RLI6, RLI7, RBE9, on the model HDM, HTM from MIDMOD) and burning (LIBURN, or LIBURN2 with the Bulirsch–Stoer option), mixing, and abundances below $10^{-24}$ set to zero. The Bulirsch–Stoer option is unused?: LBURS is set false at line 140↗ and never changed, so the LIBURN2 branch (lines 414–418↗) and the call to BURSMIX (lines 455–461↗) are never reached. So BURSMIX and BSROTMIX are not reached at 7be4133; their suspected bugs (bursmix.f, suspected bug 1, bsrotmix.f, suspected bug 1) are marked so.
9. At the end, $\Omega$ is recomputed from the final $j$ with GETROT.
10. If LPRT0 is set, it prints the surface quadrupole term: with $R$ from $L$ and $T_\mathrm{eff}$, $Q_\mathrm{surf} = (r_M/R)^3\,\phi_M$, $\Phi = GM/R$, and $-\tfrac32 Q_\mathrm{surf}/\Phi$. The argument BL is first overwritten with $\log_{10}$ HL(M) (line 472↗), the value and units (solar) that CRRECT gives it (crrect.f, line 112↗); the units match, as line 473 adds CLSUNL to it (settled from the code, Nov 9, 2026).