rotgrid.f

For the transport of angular momentum and the mixing it drives, maps one radiative region (shells IBEG to IEND) onto an auxiliary grid of points equally spaced in a mesh coordinate $\chi$, and computes there what the diffusion solver needs: each zone's mass, moment of inertia, angular momentum and angular velocity at the zone centers, and the diffusion coefficients, multiplied by the Jacobian of the change of variables, at the zone edges.

• Who calls it: SECULR (the secular transport of angular momentum and composition), for each region it treats.
• What it calls: GETGRID, which defines the coordinate $\chi$ (CHI at the model's shells, ECHI at the new points, spacing DCHI, NTOT points, and their masses ES1; all in EGRID_MOD), and OSPLIN for every interpolation, through the work arrays of SPLIN_MOD.
• Inputs (arguments): the model's log density HD, log radius HR and radius HRU, log pressure HP, luminosity HL, log mass HS, masses HS1 (enclosed) and HS2 (of each shell), moments of inertia HI, specific angular momenta HJMSAV, angular velocities WSAV, the convection flags LCZ, and the diffusion coefficients COD (angular momentum) and COD2 (composition). Also, from modules: the mesh tolerances HPTTOL, the switch LDIFAD and the factor FW (all in NMLPAR_MOD, formerly /CTOL/, /VARFC/ and /VMULT/), the energy generation EPSM (MDPHY_MOD), with LDIFAD the advective and diffusive velocities VESA, VESD (DIFAD2_MOD), and the constants (CONSTANTS_MOD).
• Outputs: on the new grid, EM (zone masses), EI (moments of inertia), EJ (angular momenta), EW (angular velocities), ECOD, ECOD2 (the two diffusion coefficients with their factors), and with LDIFAD ECOD3, ECOD4 in DIFAD_MOD; DR is set to DCHI. SECULR solves the diffusion on this grid and maps the result back. LDUM2 returns true, and nothing else is done, when the region has only one interface (IEND - IBEG <= 1).

1. Zone centers. $\log r$, $I/(m r^2)$ (EI0), $j = J/m$ and $\Omega$ are interpolated from the shells onto the new points. Then, for the interior zones, $$\Delta m_k = \tfrac12\,(m_{k+1} - m_{k-1}),\qquad I_k = \frac{I}{m r^2}\Big|_k\,\Delta m_k\,r_k^2,\qquad J_k = j_k\,\Delta m_k .$$
2. The two end zones run from the half-way point to the next equally spaced point out to the half-way point to the next original shell (or the center, or the surface, where the mass is $10^{\texttt{HSTOT}}$). If the region borders a convection zone, all of that zone's shells are added to the end zone, so a convection zone moves as one block. The end zones take $r^2$ from the shells IBEG and IEND (HRU, line 138↗, line 165↗) rather than from the interpolated $\log r$; it is the same value, since the first and last grid points lie at those shells (ECHI(1) = CHI(1), and ECHI(NTOT) = CHI(NTAB) up to rounding, since here NTOT = NTAB; getgrid.f, line 62↗, getgrid.f, lines 72–81↗), where the spline passes through the shells' values (settled from the code, Nov 9, 2026). The loops that add a convection zone also set I0 = I + 1 and I1 = I - 1 (line 143↗, line 170↗), with I left over from the zone-center loop rather than the loop variable II; I0 and I1 are read nowhere (unused?).
3. Zone edges. The diffusion coefficients are interpolated at $\chi_k - \tfrac12\Delta\chi$, assumed constant beyond the region's ends. With LDIFAD, the advective and diffusive terms $0.2\cdot4\pi f_W\,v$ are interpolated too, and the composition coefficient is added to the diffusive one (following Zahn 1992, as the comment says). These ECOD3, ECOD4 (lines 213–238↗, lines 275–280↗) are unused?: with LDIFAD, SECULR calls DADCOEFT, which sets ECOD3 and ECOD4 afresh at every iteration (dadcoeft.f, lines 268–270↗) before it reads them, and DCOEFT, which also reads them, runs only without LDIFAD (CHECKJ prints them after DADCOEFT); so the composition coefficient added here does not reach the solver. The factor 0.2 in FAC is the same as DADCOEFT's (a question on that page).

4. The Jacobian. The mesh coordinate is defined by the comment as $$\chi = \frac{\log m}{\delta_m} + \frac{L}{\delta_L} - \frac{\log P}{\delta_P},$$ with $\delta_m$ = HPTTOL(2), $\delta_L$ = HPTTOL(9) times the surface luminosity, $\delta_P$ = HPTTOL(11), so that $$\frac{d\chi}{dr} = \frac{4\pi\rho r^2}{\ln 10\;\delta_m\,m} + \frac{4\pi\rho r^2\,\epsilon}{\delta_L} + \frac{G m \rho}{\ln 10\;\delta_P\,P\,r^2}.$$ The composition coefficient is multiplied by $\rho r^2\,d\chi/dr$, and the angular-momentum coefficients by $\rho r^2\,(I/m)\,d\chi/dr$ (for a thin spherical shell, $I/m = \tfrac23 r^2$, the comment notes), each interpolated in $\log_{10}$.

   The comment's $d\chi/dr$ omits $P$. The comment writes the third term as RHO*GM/(LN10*DP*R**2) (line 248↗); the code computes $G m\rho/(P r^2)$, with -HP(II) in the exponent (lines 251–252↗), which is the correct $d\log P/dr$ term, so the comment should read RHO*GM/(LN10*DP*P*R**2).