vcirc.f

690 lines · 1 definitions · on GitHub (7be4133)

Python version: vcirc.py, 431 lines, on the right (switch to Python above the code).

Variables on this page (62)

alm, amum, cc13, cc23, cln, ddelm, dgsf, dss, eqqcod, facd2, facd3, fact1, fact2, fact4, fact5, fact6, fact7, fact8, fc, fes, fes1, fes2, fes3, fesc, festime, fgsf, fgsfc, fgsfk, fmu, fss, fv0, gg, hgm, ies, igsf, imu, json, ldifad, lnoj, qqcod, quad, rat, rcrit, rm, thdifm, thm, thn, thp, ves, ves0, vesa, vesa0, vesd, vesd0, vesd2, vesd3, vgsf, vgsf0, viscm, vmu, vss, vss0

Purpose

VCIRC estimates the velocities of the rotational mixing processes at each interface between shells of a radiative region: meridional (Eddington–Sweet) circulation, the Goldreich–Schubert–Fricke (GSF) instability, and the secular and dynamical shear instabilities. Its header comment: "skipping convective regions, check for instability against each secular mechanism; if unstable, compute a circulation velocity." The sum of the velocities, HV, sets the diffusion coefficient for angular momentum and mixing: each "velocity" becomes a diffusion coefficient over a length, $D = 4\pi f_w\,v\,r$ (CODIFF), so that, for example, $v_\mathrm{SS} \propto K/r$ gives $D_\mathrm{SS} \propto K$.

How it fits

Inputs: arguments from SECULR: $\log r$ and $r$ (HR, HRU), the range IMIN–IMAX, the iteration IT, LCZ, OMEGA, the shear QWLNR (in SECULR, $\tfrac14(\Omega_i - \Omega_{i-1})$ times a radial factor), its limit for dynamical shear QWRMAX (from SETUPV), the total luminosity CLTOT, the time step DT, $\log P$, and the interpolation weights FACI (filled by SETUPV; formerly COMMON INTFAC, now owned by SECULR). Through modules, mostly filled by SETUPV: the velocity factors FES1–FES3, FGSFK, FACT6 (TAUKH_MOD), FACT1–FACT8 (VFACT_MOD), the shell-midpoint quantities (INTVAR_MOD: RM, HGM, ALM, …), the physics at each shell (MDPHY_MOD: $\mu$ AMUM, the viscosity VISCM, the thermal diffusivity THDIFM, …), the quadrupole moment and gravity (QUADRU_MOD: QUAD, GG), the $\theta$ terms (OLDROT2_MOD), and the namelist parameters in NMLPAR_MOD (formerly VMULT, VMULT2, NOTRAN, VARFC): the scale factors FES, FGSF, FSS, FMU, FC, FESC, FGSFC, RCRIT and the switches IES, IGSF, IMU, LNOJ, LDIFAD.
Who calls it: only SECULR, before and inside its iterations; SECULR may then call GETFC with HV.
Outputs: HV (total velocity per interface) and LDO (true if any is positive); in TEMP2_MOD, VES, VGSF, VSS, VMU and their previous values, read by CHECKJ, CODIFF, PUTSTORE and STITCH; with LDIFAD, the advective and diffusive parts VESA, VESD (DIFAD2_MOD, shared with CHECKJ and ROTGRID) and VESD2, VESD3, DSS, DGSF, FV0 (DIFAD3_MOD, shared with DADCOEFT and SETUPV); RAT (QUADD_MOD), read by CHECKJ and FPFT; and FESTIME, FV0 rescaled in place. The midpoint values VISCMI, THDIFMI, AMUMI, WM and the previous VMU, VMUP, are now local arrays kept between calls (formerly COMMON INTVR2 and PREVMU).

The method

At the interface between shells $i-1$ and $i$, with $q = d\Omega/d\ln r$ there (QWLNR, from SECULR); $q_\mathrm{max}$ (QWRMAX, from SETUPV) the largest shear stable against the dynamical shear instability, from the Richardson criterion $\mathrm{Ri} = N^2/q^2 \ge \tfrac14$; $F_1$, $F_2$, $F_3$ = FES1–FES3 and $F_K$ = FGSFK, from SETUPV:

Midpoint values. The viscosity $\nu$ and thermal diffusivity $K$ by 4-point interpolation in their logarithms (weights FACI); $\mu$ and $\bar\Omega$ (WM) as the means of the two shells.
Meridional circulation (Zahn 1991, by the comments): $v_\mathrm{ES} = \bar\Omega^2 (F_1 + \bar\Omega^2 F_2)$, multiplied by the quadrupole correction of Zahn (1992), $$R = \frac{\Phi_s + \Phi_2}{\Phi_\mathrm{std}},\quad \Phi_\mathrm{std} = \frac{4}{3}\,\frac{\bar\Omega^2 r}{\bar g},\quad \Phi_s = \frac13\,\bar\Omega^2\,\frac{d}{dr}\!\left(\frac{r^2}{g}\right),\quad \Phi_2 = -\frac{d}{dr}\!\left(\frac{Q}{g}\right),$$ with a warning printed if $R < 0$. Without LDIFAD, $v_\mathrm{ES} = |R\,v_\mathrm{ES} - F_3\,\bar\Omega\,d\Omega/dr|$. With LDIFAD (MHP 05/02, 06/02), it is split into an advective part, $v_A = \tfrac15(R\,v_\mathrm{ES} + \dot\theta)$ with $\dot\theta$ = THN, and diffusive parts of first to third order.

IES = 2 (Kippenhahn & Möllenhoff 1974, restored MHP 4/26) replaces this: $v_\mathrm{ES} = |F_1'(F_2' R\,\bar\Omega^2 + F_6\,\Delta(\Omega r)^2/\Delta r)|$, with $F_1'$, $F_2'$, $F_6$ = FACT1, FACT2, FACT6 and $R$ = PHIS2/PHISTD (the quadrupole term is computed, then dropped); QW, THN and the LDIFAD parts are not set (bug 4).

> **Question** ([lines 222–223](#L222-L223)<a class='gh' href='https://github.com/yreclab/yrec/blob/7be4133/src/vcirc.f#L222-L223' title='these lines on GitHub'>↗</a>). **`RAT` with `IES = 2`.** `RAT(I) = (PHIS2+PHI2)/PHISTD` is overwritten at once by `RAT(I) = PHIS2/PHISTD`, so the quadrupole term is dropped, unlike the default branch ([line 241](#L241)<a class='gh' href='https://github.com/yreclab/yrec/blob/7be4133/src/vcirc.f#L241' title='these lines on GitHub'>↗</a>). Intended?

Inhibition by $\mu$ gradients (Zahn & Maeder 1998): every velocity is multiplied by $$f_\mu = \frac{\max(\Delta\nabla, 10^{-3})}{\max(\Delta\nabla + \nabla_\mu, 10^{-3})},\qquad \nabla_\mu = \frac{\Delta\ln\mu}{\Delta\ln P},$$ with $\Delta\nabla$ = DDELM (in SETUPV, $\nabla_\mathrm{ad} - \nabla + \nabla_\mu$, at least a minimum). The loop also multiplies FESTIME and FV0 by $f_\mu$ in place (bug 1; FV0 is read by no other routine: unused?), and computes V2 (line 314 ↗), which is never used: unused?.

Question 1 (lines 297–320 ↗). The $\mu$ loop does not skip convective interfaces. The meridional, GSF, and shear loops skip interfaces inside convection zones (line 230 ↗, line 352 ↗); this loop does not, and rescales there the LDIFAD parts VESA, VESD, VESD2, VESD3, which loop 31 did not set at those interfaces in this call, and FESTIME. Harmless, or should it skip them?

IMU = 3 (Kippenhahn & Möllenhoff 1974, restored MHP 4/26) replaces this with "fictitious $\mu$ currents" at each radiative interface, $v_\mu = |f_\mu\,F_1' F_4\,\bar\mu\,(d\mu/dr)/\sqrt{\nu c/K}|$ (VMU, with $f_\mu$ = FMU, $F_4$ = FACT4, $c$ = CC) and $v_{\mu2} = f_\mu\,\bar g\,|d\mu/dr|/(\bar\mu\,\bar\Omega^2)$ (VMU2), which are subtracted in step 6. CC, the speed of light by the comment, is never set (bug 2). At a convective interface VMU2 keeps the value of an earlier call, but it only divides velocities that are zero there (lines 152–153 ↗, line 230 ↗, line 352 ↗), and its $\bar\Omega^2$ needs $\Omega > 0$, which the logarithms of step 4 already need (settled from the code, Nov 9, 2026). 4. GSF (Kippenhahn 1980): if specific angular momentum decreases outward ($\Omega_i r_i^2 < \Omega_{i-1} r_{i-1}^2$), $$v_\mathrm{GSF} = F_K\,K\,\frac{\bar\Omega^2}{r}\,\max\!\left(2\left|\frac{d\ln j}{dr}\right|, \tfrac14\left|\frac{d\ln\Omega}{dr}\right|\right),$$ otherwise $v_\mathrm{GSF} = \tfrac14 F_K K\,\bar\Omega^2\,|d\ln\Omega/dr|/r$ (lines 396–405 ↗). With IGSF = 0 or 2 it is damped by viscosity: zero if $|q| < q_\nu = 2\sqrt{\nu/K}\,q_\mathrm{max}$, else multiplied by $\sqrt{(|q| - q_\nu)/q_\nu}$; with IGSF = 0 also by the $\mu$ gradient, likewise.

> **Question** ([lines 396–402](#L396-L402)<a class='gh' href='https://github.com/yreclab/yrec/blob/7be4133/src/vcirc.f#L396-L402' title='these lines on GitHub'>↗</a>). **$j$ decreasing outward.** This step treats $\Omega_i r_i^2 < \Omega_{i-1}r_{i-1}^2$ as a stronger GSF case, still at a thermal-diffusion rate ($\propto K$). A specific angular momentum that decreases outward is dynamically unstable (the Rayleigh, or Solberg–Høiland, criterion). Should this case be handled as dynamical (immediate redistribution), as for the dynamical shear, rather than at the GSF's thermal rate?

Shear (Endal & Sofia, paper II; Zahn 1991). If $|q| > q_\mathrm{max}$ (dynamical shear), $v_\mathrm{SS} = \tfrac{8}{45} K (|q|/q_\mathrm{max})^2/r$, and $v_\mathrm{GSF}$ is rescaled to the gradient at the limit (a message is written to unit 6, the screen, for every such interface at every call, lines 494–496 ↗): With $\Omega_\mathrm{min} = \Omega_i + (\Omega_{i-1} - \Omega_i)\,q_\mathrm{max}/|q|$ (WMIN), the value of $\Omega_{i-1}$ at which $|q|$ would equal $q_\mathrm{max}$, the gradients of $\ln\Omega$ and $\ln j$ from $\Omega_i$ to $\Omega_\mathrm{min}$ replace those from $\Omega_i$ to $\Omega_{i-1}$ in the formulae of step 4, the branch being chosen from the actual $\Omega_i r_i^2$ and $\Omega_{i-1} r_{i-1}^2$ (lines 445–497 ↗). The rescaling divides by no zero: The branch runs only when $|q| > q_\mathrm{max}$, and $q$ is proportional to $\Omega_i - \Omega_{i-1}$ (seculr.f, line 177 ↗), so the two differ and DLNWDR0 and FX2 are positive (settled from the code, Nov 9, 2026); the logarithms need $\Omega > 0$. WMID, computed at lines 451–455 ↗, is never used: unused?. Otherwise the secular shear is unstable if $$|q| > q_\mathrm{crit} = \sqrt{\frac{\nu}{K}\,\frac{R_\mathrm{crit}}{8}}\;q_\mathrm{max},$$ and then $v_\mathrm{SS} = \tfrac{8}{45} K\,\left[(|q| - q_\mathrm{crit})/q_\mathrm{max}\right]^2/r$, with $q_\mathrm{crit}$ raised by a $\mu$ gradient if that gives the larger limit.

Question 2 (lines 445–447 ↗). Dynamical shear at a thermal rate. For $|q| > q_\mathrm{max}$, $v_\mathrm{SS} = \tfrac{8}{45}K(|q|/q_\mathrm{max})^2/r$ is scaled by the thermal diffusivity $K$, although the dynamical shear instability acts on a dynamical timescale and would mix far faster. Why is its rate tied to $K$; is the intent only to make it large enough to remove the excess shear within a step?
Scale and average. The velocities are multiplied by the user's factors FES, FGSF, FSS. With IMU = 3, after VMU is averaged with its previous value (IT > 1), the $\mu$ currents are subtracted instead: $v_\mathrm{ES} \to \max(0, F v_\mathrm{ES} - \sqrt{F v_\mathrm{ES}\,r}\,\sqrt{f_c f_{c,\mathrm{ES}}}\,v_\mu)$ with $F$ = FES, $f_c$ = FC, $f_{c,\mathrm{ES}}$ = FESC, or $|F v_\mathrm{ES}|/(1 + v_{\mu2})$ where that is not positive; likewise $v_\mathrm{GSF}$ with FGSF and FGSFC (but with FES where the result is not positive: bug 3); $v_\mathrm{SS} \to \max(0, $ FSS $v_\mathrm{SS})$. with LDIFAD, the GSF and shear velocities times $r$ join the diffusive part. After the first iteration (IT > 1) each is averaged with the previous one. HV = $v_\mathrm{ES} + v_\mathrm{GSF} + v_\mathrm{SS}$ (zero below $10^{-20}$). Without IMU = 3, VMU is only zeroed (line 155 ↗) and copied to VMUP (line 138 ↗): unused? there.

Question 3 (line 595 ↗; lines 269–271 ↗; line 605 ↗). Sign of the meridional velocity. With LDIFAD, the circulation is signed, but VES is made positive by ABS(FES*VES) for every interface before HV is formed, so HV never sees a negative VES; only VESA (the advective part) keeps its sign. So HV takes $|v_\mathrm{ES}|$ of a signed circulation: intended?
No transport in the core (LNOJ, 9/93): from IMIN outward, the velocities are zeroed while $L(r)/L_\mathrm{tot} < 0.99$, "to avoid numerical problems in the H-burning shell."

Question 4 (lines 671–687 ↗). LNOJ's threshold. The core is defined as $L(r)/L_\mathrm{tot} < 0.99$, hard-coded. Should it be a parameter?

Design

Many older prescriptions are kept as comments (James & Kahn 1971 GSF, older $\mu$ treatments); IES = 2 and IMU = 3 select the restored Kippenhahn & Möllenhoff (1974) treatments, and any other value (the defaults are 1) the Zahn (1991) circulation and the Zahn & Maeder (1998) $\mu$ factor.

From YREC's documentation

Physics Options: IGSF chooses the stability condition for the GSF and ABCD instabilities; 1 is recommended. In the code, 1 applies none, 2 the condition on the angular-velocity gradient against viscosity and thermal diffusion, and 0 adds a mean-molecular-weight condition. LDIFAD, experimental, treats meridional circulation as diffusion plus advection; FC scales the mixing coefficients relative to angular momentum.

Suspected bugs, checked by running the code

Each bug is marked by what running the code showed: (changes results, large) or (small), the fix moving a run clearly beyond a null change; (stops the run); (wrong output), the fix changing only what is written; (harmless in the runs), within the noise; (never occurred), the code ran but the fault did not; (not reached); (cannot be reached), a guard prevents it; (withdrawn), the suspicion was wrong.

1. (cannot be reached) The $\mu$ factor compounds on `FESTIME`.

(line 256 ↗; lines 304–317 ↗) The $\mu$ loop multiplies FESTIME (and FV0) by $f_\mu$ in place at every call. SETUPV sets them once per call of SECULR (setupv.f, line 337 ↗), before its iterations, and SECULR calls VCIRC once per iteration (seculr.f, lines 190–191 ↗, seculr.f, lines 230–231 ↗), so at iteration $n$, THN = FESTIME*(THM*QW-THP)/DT carries the $n-1$ earlier factors ($f_\mu^{\,n-1}$ when $\mu$ does not change); and THN enters VESA and VES (line 257 ↗, line 269 ↗), which the same loop then multiplies by $f_\mu$ once more. FESTIME is read nowhere else, and FV0 nowhere at all. What follows: THN wrong from the second iteration on, only with LDIFAD (off by default) and ITDIF2 > 1 (the default is 1).

A possible fix:

vcirc.f

@@ lines 313–320 @@
             DR = HRU(I) - HRU(I-1)
             V2 = 0.2D0*(Q0+Q1)*(OMEGA(I)-OMEGA(I-1))/DR/WM(I)
 C ADD MU GRADIENTS TO VELOCITY ESTIMATES
-            FESTIME(I) = FESTIME(I)*QQQ
-            FV0(I) = FV0(I)*QQQ
+C 11/09/26 KC via Claude
+C code changed to stop rescaling FESTIME and FV0 in place: THN takes
+C the mu factor through VES and VESA, and the rescaling compounded
+C over SECULR's iterations.
+C           FESTIME(I) = FESTIME(I)*QQQ
+C           FV0(I) = FV0(I)*QQQ
 C         WRITE(*,911)I,WM(I),VESA(I),V2,Q0,Q1,THN(I),VES(I),QQQ
 C 911           FORMAT(I5,1P8E12.3)
          END DO

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): The rescaling ran in the four rotating test cases with IMU ≠ 3 (from about 0.65 to 2.2 million times, at SECULR's iterations 1 to 10, with $f_\mu$ from 0.13 to 1.28), so FESTIME does compound. But it is read only for THN in the LDIFAD branch (line 256 ↗), and PARMIN stops any run that sets LDIFAD with rotation (parmin.f, lines 999–1006 ↗), so the wrong THN cannot be reached at 7be4133; the condition line never fired. What remains is a factor that compounds unread. Intent: none found.

2. (changes results, small) Checked by KC (Nov 10, 2026) `CC` is never set.

(lines 121–122 ↗; lines 286–288 ↗) With IMU = 3, TKHS = SQRT(VISCMI(I)*CC/THDIFMI(I)) uses the local CC, whose DATA statement (the speed of light) is commented out, so CC is undefined: zero if the compiler zeroes saved storage, which makes TKHS zero and VMU infinite (or NaN where $\mu$ is uniform). Then the reduction of step 6 (lines 550–556 ↗) never leaves a positive VEST, and VES always falls back to ABS(FES*VES(I)/FM). Only with IMU = 3 (the default is 1).

A possible fix:

vcirc.f

@@ lines 119–125 @@
      *          QWLNR(JSON),QWRMAX(JSON),HV(JSON),VMU2(JSON),
      *          HP(JSON)
 C SPEED OF LIGHT.
-C       DATA CC/2.99792458D10/
+C 11/09/26 KC via Claude
+C code changed to restore the speed of light, CC, used by TKHS.
+      DATA CC/2.99792458D10/
       SAVE

 C   SKIPPING CONVECTIVE REGIONS, CHECK FOR INSTABILITY AGAINST EACH

How the code itself shows it. CC gets a value nowhere: Its only DATA line is commented out (line 122 ↗, C DATA CC/2.99792458D10/), and no other line assigns it, so the speed of light that this line comments out is the only value the code ever gave it, the value it had before the cleanup, which is what the fix above restores. The units, however, do not fit the comment above the formula, which calls TKHS the square root of Kippenhahn's timescale $t_\mathrm{KH}^*$ (lines 284–285 ↗): VISCOS computes the kinematic viscosity (viscos.f, line 47 ↗) and the radiative diffusivity THDIF $= 16\sigma T^3/(3\kappa\rho^2 c_P)$ (viscos.f, lines 104–107 ↗), both in cm$^2$ s$^{-1}$, so VISCMI*CC/THDIFMI with CC $= c$ has units of cm s$^{-1}$, and its square root is not a time. Was CC meant to be the speed of light here, or another constant (a length, or $c$ with a factor that the comment's equation (12) supplies)? The run below measures what restoring $c$ does; the right value is for the developers.

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): The case run with IMU = 3 (testsuite/Test_m2p0_feh+0p0_vzams150_IES2_TAMS) took the path 472,022 times, at shells 2 to 748, and CC was 0 every time, so TKHS was zero; the run, in the release build, went on without a message. With the speed of light restored, the track changes far beyond the noise: $\log L$ moves by up to 0.0056 dex, $\log T_\mathrm{eff}$ by 0.0020, and central $\log T$ by 0.0055 (at its end, 1.2 Gyr), against about $2\times10^{-8}$ for a null change (nearly every opacity, all with $3 \le \log T \le 9$ and $-6 \le \log\rho \le 9.75$, times $1 + 10^{-10}$, on the same case), and central $X$ by 0.00016 (at 0.076 Gyr), with the same 546 models. The other two cases with IMU = 3 (testsuite/Test_m2p0_feh+0p0_vzams150_NoJMixIES2_TAMS and testsuite/Test_m2p0_feh+0p0_vzams150_DRIES2_TAMS) were not run. Intent: The speed of light was commented out as an unused variable (commit c6a51f3, "commented out unused variables detected by ifx"), when the IMU = 3 branch was itself commented out; when that branch was switched back on (commit 49dfbec), the DATA line (lines 121–122 ↗) was not restored: a slip in restoring the branch, though whether $c$ is the right constant there is asked below.

3. (never occurred) `FES` for the GSF velocity.

(lines 564–568 ↗) With IMU = 3, when the reduced GSF velocity is not positive, VGSF(I) = ABS(FES*VGSF(I)/FM) uses FES, the circulation's factor, where the rest of the GSF branch uses FGSF (line 559 ↗). Only with IMU = 3 and FES different from FGSF (both default to 1).

A possible fix:

vcirc.f

@@ lines 562–568 @@
                VGSFT = 0.0D0
             ENDIF
             IF(VGSFT.LE.0.0D0)THEN
-               VGSF(I) = ABS(FES*VGSF(I)/FM)
+C 11/09/26 KC via Claude
+C code changed to scale the GSF velocity by FGSF, not FES.
+               VGSF(I) = ABS(FGSF*VGSF(I)/FM)
             ELSE
                VGSF(I) = VGSFT
             ENDIF

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): In testsuite/Test_m2p0_feh+0p0_vzams150_IES2_TAMS the fallback ran 473,160 times, but with VGSF = 0 every time, so the factor multiplied zero; and FES = FGSF = 1 in every namelist, so the condition line never fired. A run with IMU = 3, FES ≠ FGSF, and a positive GSF velocity would reach it; none was run. Intent: none found.

4. (cannot be reached) `IES = 2` with `LDIFAD` reuses the previous call's parts.

(lines 202–226 ↗; lines 602–607 ↗; lines 629–632 ↗) The IES = 2 loop sets only VES and RAT; THN, VESA, VESD, VESD2, VESD3 and EQQCOD keep the values of the previous call (they are not cleared at lines 150–156 ↗). With LDIFAD, they are then multiplied by $f_\mu$ and by FES again, and VESD gains RM*(VGSF+VSS) again, at every call, so it grows from call to call. Only with IES = 2 and LDIFAD (both off by default).

A possible fix: a choice for the developers: stop with a message when IES = 2 and LDIFAD are both set, or compute the LDIFAD parts in the IES = 2 loop as the default loop does.

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): testsuite/Test_m2p0_feh+0p0_vzams150_IES2_TAMS took the IES = 2 path 1,142 times, but never with LDIFAD: PARMIN stops any run that sets LDIFAD with rotation (parmin.f, lines 999–1006 ↗), so the bug cannot be reached at 7be4133. What remains is code that no run takes.

Definitions

subroutineVCIRC

lines 3–690 · GitHub

Velocities of meridional circulation (Zahn 1991, with Zahn 1992's quadrupole term, or Kippenhahn & Möllenhoff 1974), the GSF instability (Kippenhahn 1980) and secular and dynamical shear at each radiative interface, with $\mu$-gradient inhibition (Zahn & Maeder 1998, or Kippenhahn & Möllenhoff's $\mu$ currents), the user's scale factors and averaging between iterations.

Uses modules CONSTANTS_MOD, DIFAD2_MOD, DIFAD3_MOD, DIMENSIONS_MOD, EGRIDDER_MOD, INTVAR_MOD, MDPHY_MOD, NMLPAR_MOD, OLDROT2_MOD, QUADD_MOD, QUADRU_MOD, TAUKH_MOD, TEMP2_MOD, VFACT_MOD

Called by SECULR ×2

Change log

The changes since rp_1

Changed since rp_1. The options IES = 2 (Kippenhahn & Möllenhoff circulation) and IMU = 3 (their $\mu$ currents, VMU and the restored VMU2, with the scale factors applied in the same step) are restored (MHP 4/26); FACI is now an argument; the midpoint arrays and VMUP are local; declarations: COMMON blocks to modules, IMPLICIT NONE. VMU is computed again with IMU = 3.

Changed since a8a5950. In the dynamical-shear branch, WMIN is taken from the centred gradient QWLNR that triggers the branch, so that it lies between $\Omega_i$ and $\Omega_{i-1}$ and is positive; the plain difference quotient used before could put it beyond $\Omega_{i-1}$, or at or below zero (a NaN in its logarithm). The GSF branch is chosen from the actual profile, as its normalization FX2 is, not from the limited one.

1C $$$$$$ 2C SUBROUTINE VCIRC(HJM,HR,HRU,IMIN,IMAX,IT,LCZ,M,OMEGA,LDO, ! KC 2025-05-31 3 SUBROUTINE VCIRC(HR,HRU,IMIN,IMAX,IT,LCZ,M,OMEGA,LDO, 4 * * QWLNR,QWRMAX,HV,CLTOT,DT,HP,FACI) 5C 10/15/26 MHP via Claude 6C code changed to take JSON from DIMENSIONS_MOD instead of declaring 7C PARAMETER(JSON=5000) here. 8 USE DIMENSIONS_MOD 9C 10/14/26 MHP via Claude 10C code changed to convert COMMON/VFACT/ into VFACT_MOD. 11 USE VFACT_MOD 12C 10/14/26 MHP via Claude 13C code changed to convert COMMON/TEMP2/ into TEMP2_MOD. 14 USE TEMP2_MOD 15C 10/14/26 MHP via Claude 16C code changed to convert COMMON/TAUKH/ into TAUKH_MOD. 17 USE TAUKH_MOD 18C 10/14/26 MHP via Claude 19C code changed to convert COMMON/QUADRU/ into QUADRU_MOD. 20 USE QUADRU_MOD 21C 10/14/26 MHP via Claude 22C code changed to convert COMMON/QUADD/ into QUADD_MOD. 23 USE QUADD_MOD 24C 10/14/26 MHP via Claude 25C code changed to convert COMMON/OLDROT2/ into OLDROT2_MOD. 26 USE OLDROT2_MOD 27C 10/14/26 MHP via Claude 28C code changed to convert COMMON/MDPHY/ into MDPHY_MOD. 29 USE MDPHY_MOD 30C 10/14/26 MHP via Claude 31C code changed to convert COMMON/INTVAR/ into INTVAR_MOD. 32 USE INTVAR_MOD 33C 10/14/26 MHP via Claude 34C code changed to convert COMMON/EGRIDDER/ into EGRIDDER_MOD. 35 USE EGRIDDER_MOD 36C 10/14/26 MHP via Claude 37C code changed to convert COMMON/DIFAD3/ into DIFAD3_MOD. 38 USE DIFAD3_MOD 39C 10/14/26 MHP via Claude 40C code changed to convert COMMON/DIFAD2/ into DIFAD2_MOD. 41 USE DIFAD2_MOD 42C 10/14/26 MHP via Claude 43C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD. 44 USE NMLPAR_MOD 45C 10/14/26 MHP via Claude 46C code changed to convert COMMON/CONSTANTS/ into CONSTANTS_MOD. 47 USE CONSTANTS_MOD 48C 10/20/26 MHP via Claude 49C code changed to replace implicit typing with IMPLICIT NONE and 50C explicit declarations. 51 IMPLICIT NONE 52 REAL*8 AMUMI,CC,CLTOT,DDEL,DDTEST,DDTEST2,DELL,DLNJMDR,DLNJMDR0, 53 * * DLNWDR,DLNWDR0,DR,DT,DV,FACI,FCC,FM,FX,FX2,FXX,FXX2,GMID, 54 * * HP,HR,HRU,HV,OMEGA,PHI2,PHIS2,PHISTD,Q0,Q1,Q2,QMU,QP,QQ, 55 * * QQQ,QW,QWLNR,QWR,QWRMAX,QWRMX,QWRMX2,RMID,THDIFMI,TKHS,V2, 56 * * VEST,VGSFT,VISCMI,VMU2,VMUP,WM,WMID,WMIN 57 INTEGER I,I0,I1,IMAX,IMIN,IT,M 58 LOGICAL*4 LCZ,LDO 59 60C 10/13/26 MHP via Claude 61C FACI is now a dummy argument, formerly COMMON/INTFAC/. 62 DIMENSION FACI(4,JSON) 63C 9/25/26 MHP via Claude 64C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/, 65C and COMMON/CONST3/ into a single COMMON/CONSTANTS/ 66C 10/11/26 MHP via Claude 67C code changed to remove LSCRIB from COMMON/CCOUT/: it was declared 68C positionally in all 17 files sharing this block but never read or 69C written in executable code anywhere, and no namelist file sets it. 70C 10/12/26 MHP via Claude 71C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 72C namelist variables that the physics routines read from 27 COMMON 73C blocks which also held derived or working members. Those blocks keep 74C only their non-namelist members; five of them -- NEWMX, VMULT, 75C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 76C Derivation-only namelist variables from the same blocks went to 77C COMMON/NMLBYPASS/ instead, not here. 78C 79C Declared identically in all 121 files that hold any of these members. 80C Generated from one member list, so every copy is byte-identical; the 81C linker's 'size of symbol changed' warning is the check that they are. 82C 10/13/26 MHP via Claude: removed the unused COMMON/CCOUT/ declaration 83C (no member referenced here; residue of the obsolete-option removals). 84C 10/13/26 MHP via Claude 85C code changed to replace COMMON/INTFAC/ with an explicit argument: 86C setupv.f fills FACI and vcirc.f reads it, both called only from 87C seculr.f, which now owns the storage. 88C MHP 3/92 COMMON BLOCKS ADDED FOR SHEAR VELOCITY AND MU INHIBITION. 89C 10/13/26 MHP via Claude: removed the unused COMMON/IOUNITS/ declaration 90C (no member referenced here; residue of the obsolete-option removals). 91C JvS 09/25 CHANGED VARIABLE NAMES TO BE CONSISTENT ACROSS ALL SUBROUTINES 92C COMMON/MDPHY/HAMU(JSON),CP(JSON),DELM(JSON),DELAM(JSON), 93C * DELRM(JSON),SESUM(JSON),OM(JSON),SQDT(JSON), 94C * THDIF(JSON),SVEL(JSON),VISC(JSON),EPSM(JSON) 95C MHP 06/02 ADDED FACT7 AND FACT8 FOR DIF/AD TREATMENT 96C 10/9/26 MHP via Claude 97C code changed to convert COMMON/PREVMU/ and COMMON/INTVR2/ to plain 98C local arrays: both are declared only in this file (confirmed via a 99C codebase-wide scan), so there's nothing to share via COMMON. The 100C existing blanket SAVE statement below already preserves their 101C values across calls, same as COMMON did (VMUP specifically needs 102C this: it holds the prior iteration's VMU for a damping blend when 103C IT.GT.1). 104 DIMENSION VMUP(JSON) 105 DIMENSION AMUMI(JSON),THDIFMI(JSON),VISCMI(JSON),WM(JSON) 106C 3/93 QUADRUPOLE DRIVING TERM ADDED. 107C MHP 9/93 108C 10/13/26 MHP via Claude 109C code changed to drop COMMON/VARFC/, declared here but never used. 110C MHP 11/94 ADVECTIVE VELOCITY AND DIFFUSIVE COMPONENT. 111C MHP 06/02 112C HIGHER ORDER TERMS IN M.C. VELOCITY 113C Time change of theta 114C DEFINITION TERMS FOR THE SECOND AND THIRD DERIVATIVE 115C TERMS 116C MHP 4/26 RESTORED VMU2, NOW USED 117C DIMENSION HJM(JSON),HR(JSON),HRU(JSON),LCZ(JSON),OMEGA(JSON), ! KC 2025-05-31 118 DIMENSION HR(JSON),HRU(JSON),LCZ(JSON),OMEGA(JSON), 119 * * QWLNR(JSON),QWRMAX(JSON),HV(JSON),VMU2(JSON), 120 * * HP(JSON) 121C SPEED OF LIGHT. 122C DATA CC/2.99792458D10/ 123 SAVE 124 125C SKIPPING CONVECTIVE REGIONS, CHECK FOR INSTABILITY AGAINST EACH 126C SECULAR MECHANISM.IF UNSTABLE,COMPUTE A CIRCULATION VELOCITY. 127C THE PROCEDURE FOR CALCULATING STABILITY IS AS FOLLOWS: 128C STABILITY IS CALCULATED FOR THE I/I-1 INTERFACE. 129C EACH STABILITY CRITERION IS EXPRESSED AS FOLLOWS: 130C A COMBINATION OF VARIABLES(RHO,P,ETC.)STORED AT EACH SHELL 131C MULTIPLIED BY A GRADIENT BETWEEN THE SHELLS MUST BE LESS THAN 132C SOME CRITICAL NUMBER. 133 134 LDO = .FALSE. 135C STORE THE OLD VELOCITY ESTIMATES FOR LATER USE. 136 IF(IT.GT.1)THEN 137 DO 5 I = IMIN,IMAX 138 VMUP(I) = VMU(I) 139 VGSF0(I) = VGSF(I) 140 VES0(I) = VES(I) 141 VSS0(I) = VSS(I) 142 5 5 CONTINUE 143 IF(LDIFAD)THEN 144 DO I = IMIN,IMAX 145 VESA0(I) = VESA(I) 146 VESD0(I) = VESD(I) 147 END DO 148 ENDIF 149 ENDIF 150 DO 10 I = 1,M 151 HV(I) = 0.0D0 152 VES(I) = 0.0D0 153 VGSF(I) = 0.0D0 154 VSS(I) = 0.0D0 155 VMU(I) = 0.0D0 156 10 10 CONTINUE 157C MEAN MOLECULAR WEIGHT (AS WELL AS QUANTITIES WHICH DEPEND ON IT) AND 158C THE ANGULAR VELOCITY DISTRIBUTION CHANGE DURING A DIFFUSION TIMESTEP. 159C FIND THE NEW VALUES AT THE MIDPOINT IN RADIUS BETWEEN MASS SHELLS 160C WHERE STABILITY IS EVALUATED. 4 POINT LAGRANGIAN INTERPOLATION IS USED 161C FOR THE KINEMATIC VISCOSITY AND THERMAL DIFFUSIVITY. 162C TREATMENT OF FIRST INTERFACE. 163 IF(IMIN.LT.3)THEN 164 VISCMI(2)=EXP(LOG(VISCM(1))*FACI(1,2)+LOG(VISCM(2))*FACI(2,2)+ 165 * * LOG(VISCM(3))*FACI(3,2)+LOG(VISCM(4))*FACI(4,2)) 166 THDIFMI(2)=EXP(LOG(THDIFM(1))*FACI(1,2)+LOG(THDIFM(2))*FACI(2,2) 167 * * +LOG(THDIFM(3))*FACI(3,2)+LOG(THDIFM(4))*FACI(4,2)) 168 I0=3 169 ELSE 170 I0=IMIN 171 ENDIF 172C TREATMENT OF LAST INTERFACE. 173 IF(IMAX.EQ.M)THEN 174 VISCMI(M)=EXP(LOG(VISCM(M-3))*FACI(1,M)+LOG(VISCM(M-2))*FACI(2,M) 175 * * +LOG(VISCM(M-1))*FACI(3,M)+LOG(VISCM(M))*FACI(4,M)) 176 THDIFMI(M)=EXP(LOG(THDIFM(M-3))*FACI(1,M)+LOG(THDIFM(M-2))* 177 * * FACI(2,M)+LOG(THDIFM(M-1))*FACI(3,M)+ 178 * * LOG(THDIFM(M))*FACI(4,M)) 179 I1=M-1 180 ELSE 181 I1=IMAX 182 ENDIF 183C GENERAL CASE. 184 DO 20 I = I0,I1 185 VISCMI(I)=EXP(LOG(VISCM(I-2))*FACI(1,I)+LOG(VISCM(I-1)) 186 * * *FACI(2,I)+LOG(VISCM(I))*FACI(3,I)+LOG(VISCM(I+1)) 187 * * *FACI(4,I)) 188 THDIFMI(I)=EXP(LOG(THDIFM(I-2))*FACI(1,I)+LOG(THDIFM(I-1)) 189 * * *FACI(2,I)+LOG(THDIFM(I))*FACI(3,I)+LOG(THDIFM(I+1)) 190 * * *FACI(4,I)) 191 20 20 CONTINUE 192C USE LINEAR INTERPOLATION FOR OMEGA AND MU. 193 DO I = 2,M 194 AMUMI(I) = 0.5D0*(AMUM(I)+AMUM(I-1)) 195 WM(I) = 0.5D0*(OMEGA(I)+OMEGA(I-1)) 196 END DO 197C MHP 4/26 RESTORED PRIOR KT OPTION, DEFAULT USES ZAHN & MAEDER 198C MHP 8/03 OMITTED OLD KM1974 MERIDIONAL CIRCULATION VELOCITY ESTIMATE. 199C THE IES FLAG IS THEREFORE NO LONGER IN USE. 200C 201C KIPPENHAHN AND MOLLENHOF(1974) MERIDIONAL CIRCULATION VELOCITY. 202 IF(IES.EQ.2)THEN 203 DO 30 I = IMIN,IMAX 204C SKIP CONVECTIVE INTERFACES. 205 IF(LCZ(I).AND.LCZ(I-1))GOTO 30 206C EDDINGTON CIRCULATION. VELOCITY DEFINED IN SR SETUPV. 207 DR = HRU(I)-HRU(I-1) 208 DV = ((OMEGA(I)*HRU(I))**2-(OMEGA(I-1)*HRU(I-1))**2)/DR 209C QUADRUPOLE TERM, AS PER ZAHN (1992), ADDED. 210 GMID = 0.5D0*(GG(I)+GG(I-1)) 211C THE VARIATION IN GRAVITY ON A LEVEL SURFACE CAN BE EXPRESSED IN 212C GENERAL AS 1/3 OMEGA**2 *D/DR((R**2 - PHI)/G), WHERE G IS THE 213C LOCAL AVERAGE GRAVITY (~GM/R**2) AND PHI IS THE QUADRUPOLE MOMENT. 214C THE CLASSICAL EXPRESSION NEGLECTS THE QUADRUPOLE AND THE SPATIAL 215C DERIVATIVE OF THE MASS (PHISTD). THE VELOCITY ESTIMATE FOR 216C MERIDIONAL CIRCULATION IS THEREFORE MULTIPLED BY THE RATIO OF 217C THE CORRECTED EXPRESSION TO THE CLASSICAL ONE (RAT). 218 PHISTD = 2.0D0*CC23*WM(I)**2*RM(I)/GMID 219 PHIS2 = CC13*WM(I)**2*(HRU(I)**2/GG(I)- 220 * * HRU(I-1)**2/GG(I-1))/DR 221 PHI2 = (QUAD(I-1)/GG(I-1)-QUAD(I)/GG(I))/DR 222 RAT(I) = (PHIS2+PHI2)/PHISTD 223 RAT(I) = PHIS2/PHISTD 224 VES(I) = ABS(FACT1(I)*(FACT2(I)*RAT(I)*WM(I)**2 225 * * +FACT6(I)*DV)) 226 30 30 CONTINUE 227 ELSE 228C ZAHN (1991) VELOCITY ESTIMATE. 229 DO 31 I = IMIN,IMAX 230 IF(LCZ(I).AND.LCZ(I-1))GOTO 31 231C ORIGINAL ESTIMATE,USING DG/G = W**2 R**3 / GM. 232 VES(I) = WM(I)**2*(FES1(I)+WM(I)**2*FES2(I)) 233C QUADRUPOLE TERM ADDED, DG/G COMPUTED AS PER ZAHN 1992. 234 DR = HRU(I) - HRU(I-1) 235 GMID = 0.5D0*(GG(I)+GG(I-1)) 236 PHISTD = 2.0D0*CC23*WM(I)**2*RM(I)/GMID 237 PHIS2 = CC13*WM(I)**2*(HRU(I)**2/GG(I)- 238 * * HRU(I-1)**2/GG(I-1))/DR 239 PHI2 = (QUAD(I-1)/GG(I-1)-QUAD(I)/GG(I))/DR 240 QW = WM(I)*(OMEGA(I)-OMEGA(I-1))/DR 241 RAT(I) = (PHIS2+PHI2)/PHISTD 242 IF(RAT(I).LT.0.0D0)THEN 243 WRITE(*,303)I,RAT(I),PHISTD,PHIS2,PHI2, 244 * * QUAD(I),QUAD(I-1),GG(I),GG(I-1),HRU(I), 245 * * HRU(I-1) 246 303 303 FORMAT(I5,' RAT ',1PE10.3,' POT SPH,CYL,QUA', 247 * * 3E10.3/' QUA,G,R ',6E12.3) 248 ENDIF 249 IF(.NOT.LDIFAD)THEN 250 VES(I) = ABS(RAT(I)*VES(I)-FES3(I)*QW) 251 ELSE 252C MHP 05/02 ADD FACTOR OF 1/5 HERE 253C VESA(I) = RAT(I)*VES(I) 254C VESD(I) = ABS(FES3(I)*WM(I)**2) 255C MHP 06/02 ADDED D THETA/DT TERM 256 THN(I) = FESTIME(I)*(THM(I)*QW-THP(I))/DT 257 VESA(I) = 0.2D0*(RAT(I)*VES(I)+THN(I)) 258 Q1 = FACT7(I)*WM(I)**2 259 Q2 = FACT8(I)*QW 260 Q0 = FES3(I)*WM(I)**2 261C VESD(I) = 0.2D0*ABS(FES3(I)*WM(I)**2) 262C VESD(I) = 0.2D0*ABS(Q0+Q1+Q2) 263 VESD(I) = 0.2D0*(Q0+Q1+Q2) 264C SECOND ORDER TERM 265 VESD2(I) = 0.2D0*FACD2(I)*WM(I) 266C THIRD ORDER TERM 267 VESD3(I) = 0.2D0*FACD3(I)*WM(I) 268C VES(I) = ABS(RAT(I)*VES(I)+THN(I)-FES3(I)*QW) 269 VES(I) = RAT(I)*VES(I)+THN(I)+ 270 * * (FES3(I)+FACT7(I))*QW+ 271 * * FACT8(I)*((OMEGA(I)-OMEGA(I-1))/DR)**2 272 EQQCOD(I) = QQCOD(I)*WM(I) 273 ENDIF 274 31 31 CONTINUE 275 ENDIF 276C MHP 4/26 RESTORED PRIOR KT OPTION, DEFAULT USES ZAHN & MAEDER 277C INHIBITION CAUSED BY GRADIENTS IN MEAN MOLECULAR WEIGHT. 278C INCLUDE ONLY THE ZAHN & MAEDER 1998 TREATMENT, WHERE (DEL - DEL AD) IS 279C REPLACED BY (DEL + DEL MU - DEL AD). RETAIN IMU FLAG FOR LATER USE. 280 IF(IMU.EQ.3)THEN 281C FICTITIOUS MU CURRENTS THAT OPPOSE ES CIRCULATION CALCULATED HERE. 282 DO I = IMIN,IMAX 283 IF(LCZ(I).AND.LCZ(I-1))GOTO 32 284C SQUARE ROOT OF TKH*,KIPPENHAHN,IAU#66,P.23,USING EQ.(12)FOR D ON P.25. 285C NOTE FACTOR OF G IS SUBSUMED IN FACT4. 286 TKHS = SQRT(VISCMI(I)*CC/THDIFMI(I)) 287 QMU = (AMUM(I)-AMUM(I-1))/(HRU(I)-HRU(I-1)) 288 VMU(I) = ABS(FMU*FACT1(I)*FACT4(I)*AMUMI(I)*QMU/TKHS) 289 VMU2(I)=FMU*HGM(I)*ABS(QMU)/AMUMI(I)/WM(I)**2 290 32 32 CONTINUE 291 END DO 292C LOCAL TIMESCALE ESTIMATE FOR MU INHIBITION. 293C ASSUMES V/L = OMEGA. 294 ELSE 295C ELSE IF(IMU.EQ.2)THEN 296C MHP 06/02 REPLACE WITH THE ZAHN&MAEDER 1998 PRESCRIPTION 297 DO I = IMIN,IMAX 298 QMU = LOG(AMUM(I))-LOG(AMUM(I-1)) 299 QP = CLN*(HP(I)-HP(I-1)) 300 DDEL=DDELM(I)+ QMU/QP 301 DDTEST = MAX(DDELM(I),1.0D-3) 302 DDTEST2 = MAX(DDEL,1.0D-3) 303 QQQ = DDTEST/DDTEST2 304 VES(I) = VES(I)*QQQ 305 IF(LDIFAD)THEN 306 VESA(I) = VESA(I)*QQQ 307 VESD(I) = VESD(I)*QQQ 308 VESD2(I) = VESD2(I)*QQQ 309 VESD3(I) = VESD3(I)*QQQ 310 ENDIF 311 Q1 = FACT7(I)*WM(I)**2*QQQ 312 Q0 = FES3(I)*WM(I)**2*QQQ 313 DR = HRU(I) - HRU(I-1) 314 V2 = 0.2D0*(Q0+Q1)*(OMEGA(I)-OMEGA(I-1))/DR/WM(I) 315C ADD MU GRADIENTS TO VELOCITY ESTIMATES 316 FESTIME(I) = FESTIME(I)*QQQ 317 FV0(I) = FV0(I)*QQQ 318C WRITE(*,911)I,WM(I),VESA(I),V2,Q0,Q1,THN(I),VES(I),QQQ 319C 911 FORMAT(I5,1P8E12.3) 320 END DO 321C DO 33 I = IMIN,IMAX 322C IF(LCZ(I).AND.LCZ(I-1))GOTO 33 323C DR = HRU(I) - HRU(I-1) 324C VMU(I)=FMU*HGM(I)*ABS(AMUM(I)-AMUM(I-1))/DR/AMUMI(I) 325C * /WM(I)**2 326C 33 CONTINUE 327C ALTERNATE EXPRESSION : D = D0/(1+R*DEL MU/MU) 328C THIS IS ACTUALLY SIMILAR TO THE ZM98 FORMULA, THEREFORE OBSOLETE. 329C SHOULD HAVE PUBLISHED IT - OH WELL. 330C ELSE 331C DO 34 I = IMIN,IMAX 332C IF(LCZ(I).AND.LCZ(I-1))GOTO 34 333C DR = HRU(I) - HRU(I-1) 334C DDEL = DELAMI(I)/MAX(1.0D-6,DELAMI(I)-DELMI(I)) 335C VMU(I)=FMU*RM(I)*DDEL*ABS(AMUM(I)-AMUM(I-1))/DR/AMUMI(I) 336C 34 CONTINUE 337 ENDIF 338C GSF INSTABILITY. VELOCITY DEFINED IN SR SETUPV. 339C OMIT ALL BUT THE KIPPENHAHN 1980 ESTIMATE. 340C IGSF IS RETAINED, AND THE VALUE DETERMINES WHAT IS PERMITTED TO 341C INHIBIT IT. 342C IGSF = 0: INHIBITED BY VISCOSITY AND MU GRADIENTS 343C IGSF = 1: NEITHER VISCOSITY NOR MU GRADIENTS INHIBIT 344C IGSF = 2: INHIBITED BY VISCOSITY ONLY 345C NOTE THAT THE ABCD INSTABILITY - RELATED TO THE GSF INSTABILITY - 346C WOULD NOT BE INHIBITED BY VISCOSITY/MU GRADIENTS, SO INCLUDING IT 347C IS EQUIVALENT TO SETTING IGSF=1 (THE ORIGINAL K1980 ESTIMATE). 348C THE INHIBITION FACTORS COME FROM SUBSEQUENT WORK BY SPRUIT. 349C IF(IGSF.EQ.0 .OR. IGSF.EQ.2 .OR. IGSF.EQ.1)THEN 350c KIPPENHAHN (1980) ESTIMATE 351 DO 40 I = IMIN,IMAX 352 IF(LCZ(I).AND.LCZ(I-1))GOTO 40 353C MHP 8/93 STABILITY CONDITION ADDED, NEGLECTING THE EFFECTS OF 354C MU GRADIENTS. 355 IF(IGSF.EQ.2 .OR. IGSF.EQ.0)THEN 356 QWRMX = 2.0D0*SQRT(VISCMI(I)/THDIFMI(I))*QWRMAX(I) 357 IF(ABS(QWLNR(I)).LT.QWRMX)THEN 358 VGSF(I) = 0.0D0 359 GOTO 40 360 ELSE 361 FXX = SQRT((ABS(QWLNR(I))-QWRMX)/QWRMX) 362 ENDIF 363 ELSE 364 FXX = 1.0D0 365 ENDIF 366 DR = HRU(I)-HRU(I-1) 367 IF(IGSF.EQ.0)THEN 368 QWRMX=2.0D0*SQRT(HGM(I)*ABS(AMUM(I)-AMUM(I-1)) 369 * * /DR/AMUMI(I)) 370 IF(ABS(QWLNR(I)).LT.QWRMX)THEN 371 VGSF(I) = 0.0D0 372 GOTO 40 373 ELSE 374C MHP 05/02 ADDED TESTS TO AVOID DIVIDE BY ZERO 375C FXX2 = SQRT((ABS(QWLNR(I))-QWRMX)/QWRMX) 376 QQ = 1.0D-2*QWLNR(I) 377 IF(ABS(QWRMX).LT.ABS(QQ))THEN 378 FXX2 = 1.0D0 379 ELSE 380 IF(ABS(QWRMX).LT.1.0D-20)THEN 381 FXX2 = 1.0D0 382 ELSE 383 FXX2 = SQRT((ABS(QWLNR(I))-QWRMX)/QWRMX) 384 ENDIF 385 ENDIF 386 FXX = MIN(FXX,FXX2) 387 ENDIF 388 ENDIF 389 RMID = 0.5D0*(HRU(I)+HRU(I-1)) 390 DLNWDR = ABS(LOG(OMEGA(I))-LOG(OMEGA(I-1)))/DR 391C GSF IS TRIGGERED BY D OMEGA/DZ NONZERO (I.E. ROTATION NOT 392C ON CYLINDERS,WHICH IS TRUE IN GENERAL IN OUR MODELS), OR 393C BY D LN(J/M)/DR < 0, WHICH SOMETIMES OCCURS. 394C IF D LN(J/M)/DR < 0, CHECK SECOND CRITERIA AND USE THE LARGEST 395C VELOCITY; OTHERWISE,USE ONLY D OMEGA/DZ NON-ZERO CRITERION. 396 IF(OMEGA(I)*HRU(I)**2.LT.OMEGA(I-1)*HRU(I-1)**2)THEN 397 DLNJMDR = ABS(2.0D0/RMID+(LOG(OMEGA(I))- 398 * * LOG(OMEGA(I-1)))/DR) 399 FX = MAX(2.0D0*DLNJMDR,0.25D0*DLNWDR) 400C FACTOR OF R IN THE DENOMINATOR OCCURS BECAUSE ALL THE CIRCULATION 401C VELOCITIES ARE LATER MULTIPLIED BY R (THE "LENGTH SCALE"). 402 VGSF(I) = FGSFK(I)*FX*THDIFMI(I)*WM(I)**2/RMID 403 ELSE 404 VGSF(I)=0.25D0*FGSFK(I)*THDIFMI(I)*DLNWDR*WM(I)**2/RMID 405 ENDIF 406 VGSF(I) = ABS(FXX*VGSF(I)) 407 40 40 CONTINUE 408C OMIT JAMES AND KAHN ESTIMATE 409C ELSE 410C JAMES AND KAHN (1971) ESTIMATE. 411C NOTE FACTOR OF OMEGA**4, WHICH IS PART OF THE MERIDIONAL 412C CIRCULATION VELOCITY ESTIMATE. 413C DO 41 I = IMIN,IMAX 414C IF(LCZ(I).AND.LCZ(I-1))GOTO 41 415C MHP 8/93 STABILITY CONDITION ADDED, NEGLECTING THE EFFECTS OF 416C MU GRADIENTS. 417C QWRMX = 2.0D0*SQRT(VISCMI(I)/THDIFMI(I))*QWRMAX(I) 418C IF(ABS(QWLNR(I)).LT.QWRMX)THEN 419C VGSF(I) = 0.0D0 420C GOTO 41 421C ELSE 422C FXX = SQRT((ABS(QWLNR(I))-QWRMX)/QWRMX) 423C ENDIF 424C DR = HRU(I)-HRU(I-1) 425C RMID = 0.5D0*(HRU(I)+HRU(I-1)) 426C DLNJMDR = ABS(2.0D0/RMID+(LOG(OMEGA(I))- 427C * LOG(OMEGA(I-1)))/DR) 428C IF(IGSF.EQ.3)THEN 429C DWDR = abs(OMEGA(I)-OMEGA(I-1))/DR 430C VGSF(I) = 2.0D0*FGSFJ(I)*DWDR**2 431C VGSF(I) = ABS(VGSF(I)) 432C ELSE IF(IES.EQ.2)THEN 433C VGSF(I) = 2.0D0*VES(I)*FACT3(I)*DLNJMDR**2 434C ELSE 435C VGSF(I) = 2.0D0*FGSFJ(I)*(DLNJMDR*WM(I))**2 436C VGSF(I) = ABS(FXX*VGSF(I)) 437C ENDIF 438C 41 CONTINUE 439C ENDIF 440C DIFFUSIVE AND DYNAMICAL SHEAR INSTABILITIES - REF. ENDAL&SOFIA PAPER II. 441 DO 60 I = IMIN,IMAX 442C CHECK FOR OPERATION OF DYNAMICAL SHEAR. 443C IF DYNAMICAL SHEAR IS OPERATING,SET SECULAR SHEAR VELOCITY TO MAXIMUM 444C VALUE AND COMPUTE (LARGE) DYNAMICAL SHEAR VELOCITY. 445 IF(ABS(QWLNR(I)).GT.QWRMAX(I))THEN 446 QWR = ABS(QWLNR(I)) 447 VSS(I)=8.0D0/4.5D1*THDIFMI(I)*(QWR/QWRMAX(I))**2/RM(I) 448C CORRECT GSF VELOCITY AS WELL. 449 DR = HRU(I)-HRU(I-1) 450 RMID = 0.5D0*(HRU(I)+HRU(I-1)) 451 IF(OMEGA(I).LT.OMEGA(I-1))THEN 452 WMID = OMEGA(I) + QWRMAX(I)*CLN*(HR(I)-LOG10(RMID)) 453 ELSE 454 WMID = OMEGA(I) - QWRMAX(I)*CLN*(HR(I)-LOG10(RMID)) 455 ENDIF 456C 10/30/26 MHP via Claude 457C code changed to take WMIN from the trigger's own centred gradient. 458C QWLNR is proportional to OMEGA(I)-OMEGA(I-1), so the OMEGA(I-1) at 459C which it equals QWRMAX lies this fraction of the way to OMEGA(I-1). 460C The fraction is below 1 whenever the trigger fires, so WMIN lies 461C between the two and is positive. The plain difference quotient used 462C before could put WMIN beyond OMEGA(I-1), or at or below zero. 463 WMIN = OMEGA(I) + (OMEGA(I-1)-OMEGA(I))*QWRMAX(I)/QWR 464 DLNWDR0 = abs(LOG(OMEGA(I))-LOG(OMEGA(I-1)))/DR 465 DLNJMDR0 = ABS(2.0D0/RMID+(LOG(OMEGA(I))- 466 * * LOG(OMEGA(I-1)))/DR) 467 DLNJMDR = abs(2.0D0/RMID+(LOG(OMEGA(I))-log(wmin))/DR) 468 DLNWDR = abs(LOG(OMEGA(I))-LOG(WMIN))/DR 469C MHP 05/02 ADDED FIRST BRANCH FOR IGSF=0 - SOMEHOW 470C MISSED IN CODE CHANGES EARLIER! 471c kippenhahn (1980) estimate 472C OMIT BRANCH FOR JAMES AND KAHN. 473C IF(IGSF.EQ.2 .OR. IGSF.EQ.1 .OR. IGSF.EQ.0)THEN 474C 10/30/26 MHP via Claude 475C code changed to choose the branch from the actual profile, which the 476C normalization FX2 below also uses, not from the limited one. 477 IF(OMEGA(I)*HRU(I)**2.LT.OMEGA(I-1)*HRU(I-1)**2)THEN 478 FX = MAX(2.0D0*DLNJMDR,0.25D0*DLNWDR) 479 FX2 = MAX(2.0D0*DLNJMDR0,0.25D0*DLNWDR0) 480 VGSF(I) = VGSF(I)*FX/FX2 481 ELSE 482 VGSF(I)=VGSF(I)*DLNWDR/DLNWDR0 483 ENDIF 484C ELSE IF(IGSF.EQ.3)THEN 485C DWDR = abs(OMEGA(I)-WMIN)/DR 486C DWDR0 = ABS(OMEGA(I)-OMEGA(I-1))/DR 487C VGSF(I) = VGSF(I)*(DWDR/DWDR0)**2 488C ELSE IF(IES.EQ.2)THEN 489C VGSF(I) = VGSF(I)*(DLNJMDR/DLNJMDR0)**2 490C ELSE 491C WMID0 = 0.5D0*(OMEGA(I)+OMEGA(I-1)) 492C VGSF(I) = VGSF(I)*(DLNJMDR*WMID/DLNJMDR0/WMID0)**2 493C ENDIF 494 WRITE(6,9911)I,OMEGA(I),OMEGA(I-1),WMIN 495 9911 9911 FORMAT(1X,'DYNAMICAL SHEAR-SHELL',I5,1P,' WTOP',E11.3, 496 * * ' WBOT',E11.3,' LIMIT',E11.3) 497 GOTO 60 498C *** END OF CHANGED SECTION 499 ENDIF 500C FIND MAXIMUM GRADIENT IN OMEGA ALLOWED BY SECULAR SHEAR. 501C THE RUN OF QWRMAX INPUT IS THAT ALLOWED BY THE DYNAMICAL SHEAR; 502C THE SECULAR SHEAR RICHARDSON # IS RELATED BY 503C RICHNO(SECULAR) = PRANDTL# * CRITICAL REYNOLDS#/8 * RICHNO(DYNAMICAL) 504C PRANDTL # = KINEMATIC VISCOSITY/THERMOMETRIC DIFFUSIVITY. 505C MHP 3/92 SQUARE ROOT OF PR# NEEDED, NOT PR # - ERROR CORRECTED! 506C THE VELOCITY ESTIMATE HERE IS FROM ZAHN 1991. 507 QWRMX = SQRT(VISCMI(I)/THDIFMI(I)*1.25D-1*RCRIT)*QWRMAX(I) 508 IF(ABS(QWLNR(I)).GT.QWRMX) THEN 509C UNSTABLE; CHECK FOR MU GRADIENTS. 510 IF(ABS((AMUM(I)-AMUM(I-1))/AMUMI(I)).LT.1.0D-10) THEN 511 QWRMX2 = 0.0D0 512 QWR = ABS(QWLNR(I)) - QWRMX 513 VSS(I)=8.0D0/4.5D1*THDIFMI(I)*(QWR/QWRMAX(I))**2/RM(I) 514 ELSE 515C CHECK FOR EFFECTS OF MU GRADIENT. 516C RICHNO = RHO/P*(-d lnT/d lnMU)*(del MU)/(del P)*(G/QWLNR)**2 517C WHERE -d lnT/d lnMU = (CON-1)/(1+3CON),CON=(a/3)T**4/P 518C GIVEN (1 - CON)P = CGAS*RHO*T/MU 519C FACT = (RHOM/PM)*QTMU*DMU/AMUMI/DP*HGM**2 520 QWRMX2 = 2.0D0*SQRT(MAX(1.0D-20,FACT5(I)* 521 * * ABS((AMUM(I)-AMUM(I-1))/AMUMI(I)))) 522 IF(ABS(QWLNR(I)).GT.QWRMX2) THEN 523C INTERFACE UNSTABLE WITH RESPECT TO BOTH CONDITIONS; CHOOSE THE 524C MAXIMUM GRADIENT IMPLIED BY THE SECOND CONDITION IF IT'S LARGER 525C THAN THE FIRST(I.E. IF A MU GRADIENT IS SLOWING J TRANSPORT). 526 QWRMX = MAX(QWRMX2,QWRMX) 527 QWR = ABS(QWLNR(I)) - QWRMX 528 VSS(I)=8.0D0/4.5D1*THDIFMI(I)*(QWR/QWRMAX(I))**2/RM(I) 529 ENDIF 530 ENDIF 531 ENDIF 532 60 60 CONTINUE 533C MHP 4/26 RESTORED PRIOR KT OPTION, DEFAULT USES ZAHN & MAEDER 534C NOW DETERMINE WHETHER OR NOT MU GRADIENTS ARE STEEP ENOUGH TO 535C INHIBIT TRANSPORT. MULTIPLY THE RESULTING VELOCITY ESTIMATES 536C BY THE USER DEFINED PARAMETERS FES AND FGSF. 537C IMU=3 KIPPENHAHN AND MOLLENHOF(1974)METHOD;IMU=2 LOCAL DAMPING 538C FACTOR METHOD. 539C AGAIN, OMIT OBSOLETE MU GRADIENT TREATMENTS. 540 IF(IMU.EQ.3)THEN 541 IF(IT.GT.1)THEN 542 DO I = 1,M 543 VMU(I) = 0.5D0*(VMU(I)+VMUP(I)) 544 END DO 545 ENDIF 546 DO I = IMIN,IMAX 547 FM = 1.0D0+VMU2(I) 548 RMID = 0.5D0*(HRU(I)+HRU(I-1)) 549 FCC = SQRT(FC*FESC) 550 VEST=MAX(0.0D0,FES*VES(I)-SQRT(FES*VES(I)*RMID) 551 * * *FCC*VMU(I)) 552 IF(VEST.LE.0.0D0)THEN 553 VES(I) = ABS(FES*VES(I)/FM) 554 ELSE 555 VES(I) = VEST 556 ENDIF 557 FCC = SQRT(FC*FGSFC) 558 IF(VGSF(I).GT.0.0D0)THEN 559 VGSFT=MAX(0.0D0,FGSF*VGSF(I)-SQRT(FGSF*VGSF(I)*RMID) 560 * * *FCC*VMU(I)) 561 ELSE 562 VGSFT = 0.0D0 563 ENDIF 564 IF(VGSFT.LE.0.0D0)THEN 565 VGSF(I) = ABS(FES*VGSF(I)/FM) 566 ELSE 567 VGSF(I) = VGSFT 568 ENDIF 569 VSS(I)=MAX(0.0D0,FSS*VSS(I)) 570 END DO 571C ELSE 572C IF(IMU.NE.2)THEN 573C FM = 1.0D0+VMU(I) 574C ELSE 575C FM = 1.0D0 576C ENDIF 577C VES(I)=ABS(FES*VES(I)/FM) 578C MHP 05/02 ONLY DO THIS IF MU GRADIENTS NOT 579C ALREADY ACCOUNTED FOR 580C VGSF(I)=ABS(FGSF*VGSF(I)/FM) 581C IF(IGSF.NE.0)THEN 582C VGSF(I)=ABS(FGSF*VGSF(I)/FM) 583C ELSE 584C VGSF(I) = FGSF*VGSF(I) 585C ENDIF 586C ALREADY INCLUDED - ONLY USE SCALE FACTOR 587C VSS(I)=ABS(FSS*VSS(I)/FM) 588C VSS(I)= FSS*VSS(I) 589C END DO 590C MHP 8/03 MULTIPLY VELOCITY ESTIMATES BY USER PARAMETER 591C SCALE FACTORS 592C MHP 4/26 SCALE FACTORS ALREADY INCLUDED FOR IMU=3 CASE, SO ONLY DO THIS IF NOT DONE BEFORE 593 ELSE 594 DO I = IMIN,IMAX 595 VES(I)=ABS(FES*VES(I)) 596 VGSF(I) = FGSF*VGSF(I) 597 VSS(I)= FSS*VSS(I) 598 END DO 599 ENDIF 600C MHP 11/94 601C REPEAT FOR DIF+AD 602 IF(LDIFAD)THEN 603 DO I = IMIN,IMAX 604 VESD(I)=ABS(FES*VESD(I)) 605 VESA(I)=FES*VESA(I) 606 END DO 607 ENDIF 608C IF(IMU.NE.2)THEN 609C FM = 1.0D0+VMU(I) 610C ELSE 611C FM = 1.0D0 612C ENDIF 613C FM = 1.0D0+VMU(I) 614C VESD(I)=ABS(FES*VESD(I)/FM) 615C VESA(I)=FES*VESA(I)/FM 616C END DO 617C ENDIF 618C ENDIF 619 IF(LDIFAD)THEN 620C MHP 05/02 621C CHANGED TO REFLECT THE DIFFERENT TREATMENT OF THE GSF INSTABILITY. 622C THE ENDAL AND SOFIA VARIANT HAD V=THDIF*DLNWDR*W**2/R WHICH WAS 623C THEN ADDED TO A DIFFUSION EQUATION D/DR(CON*V*R*DW/DR). 624C IN THE ZAHN FORMULISM, WE ARE SOLVING AN EQUATION OF THE FORM 625C D/DR(CON*V*W). BECAUSE THE VELOCITY ESTIMATE ITSELF INCLUDES 626C DW/DR THIS IS IN EFFECT A DIFFUSION TERM. SO...WE NEED TO 627C MULTIPLY THE ORIGINAL VELOCITY BY W AND DIVIDE BY DW/DR TO CAST 628C IT AS AN ENTRY IN THE DIFFUSION EQUATION. 629 DO I = IMIN,IMAX 630 VESD(I) = VESD(I)+RM(I)*(ABS(VGSF(I))+ABS(VSS(I))) 631 DSS(I) = RM(I)*ABS(VSS(I)) 632 DGSF(I) = RM(I)*ABS(VGSF(I)) 633C IF(VGSF(I).GT.0.0D0)THEN 634C D LN W/DR 635C DR = HRU(I) - HRU(I-1) 636C QLNWR = ABS(LOG(OMEGA(I))-LOG(OMEGA(I-1))/DR) 637C CEILING SET BY DYNAMICAL SHEAR 638C QLNWRMAX = ABS(QWRMAX(I)/(WM(I)*RM(I))) 639C TAKE THE SMALLER OF THE TWO 640C QLNWR = MIN(QLNWR,QLNWRMAX) 641C IF(QLNWR.GT.1.0D-32)VESD(I) = VESD(I)+ABS(VGSF(I)/QLNWR) 642C ENDIF 643C IF(VSS(I).GT.0.0D0)THEN 644C VESD(I) = VESD(I)+ABS(VSS(I)*RM(I)) 645C ENDIF 646 END DO 647 ENDIF 648C AVERAGE PREVIOUS AND NEW VELOCITY ESTIMATES AFTER THE FIRST ITERATION. 649 IF(IT.GT.1) THEN 650 DO 70 I = IMIN,IMAX 651 VGSF(I) = 0.5D0*(VGSF(I) + VGSF0(I)) 652 VES(I) = 0.5D0*(VES(I) + VES0(I)) 653 VSS(I) = 0.5D0*(VSS(I) + VSS0(I)) 654 70 70 CONTINUE 655C MHP 11/94 656 IF(LDIFAD)THEN 657 DO I = IMIN,IMAX 658 VESA(I) = 0.5D0*(VESA(I) + VESA0(I)) 659 VESD(I) = 0.5D0*(VESD(I) + VESD0(I)) 660 END DO 661 ENDIF 662 ENDIF 663 DO 80 I =IMIN,IMAX 664 HV(I) = VGSF(I) + VES(I) + VSS(I) 665 IF(HV(I).LT.1.0D-20)HV(I)=0.0D0 666 IF(HV(I).GT.0.0D0)LDO=.TRUE. 667 80 80 CONTINUE 668C 9/93 MIXING WITHOUT TRANSPORT ADDED. 669C ZERO OUT COEFFICIENTS IN CORE TO AVOID NUMERICAL PROBLEMS IN 670C THE H-BURNING SHELL. 671 IF(LNOJ) THEN 672 DO I = IMIN,IMAX 673 DELL = ALM(I)/CLTOT 674 IF(DELL.LT.9.9D-1)THEN 675 HV(I) = 0.0D0 676 VGSF(I) = 0.0D0 677 VES(I) = 0.0D0 678 VSS(I) = 0.0D0 679C MHP 11/94 680 VESA(I) = 0.0D0 681 VESD(I) = 0.0D0 682 ELSE 683 GOTO 81 684 ENDIF 685 END DO 686 81 81 CONTINUE 687 ENDIF 688 689 RETURN 690 END

1"""Python version of vcirc.f (YREC 7be4133): the velocities of the 2rotational mixing processes (meridional circulation, the GSF instability, 3secular and dynamical shear) at each interface of a radiative region.""" 4 5import math 6 7import numpy as np 8 9from pyrec import constants_mod 10from pyrec import difad2_mod 11from pyrec import difad3_mod 12from pyrec import dimensions_mod 13from pyrec import egridder_mod 14from pyrec import intvar_mod 15from pyrec import mdphy_mod 16from pyrec import nmlpar_mod 17from pyrec import oldrot2_mod 18from pyrec import quadd_mod 19from pyrec import quadru_mod 20from pyrec import taukh_mod 21from pyrec import temp2_mod 22from pyrec import vfact_mod 23 24JSON = dimensions_mod.json 25 26# Local arrays kept between calls (SAVE): the midpoint values (formerly 27# COMMON/INTVR2/), VMU of the previous iteration (formerly COMMON/PREVMU/) 28# and VMU2. 29viscmi = np.zeros(JSON + 1) 30thdifmi = np.zeros(JSON + 1) 31amumi = np.zeros(JSON + 1) 32wm = np.zeros(JSON + 1) 33vmup = np.zeros(JSON + 1) 34vmu2 = np.zeros(JSON + 1) 35 36# The speed of light. Suspected bug (see the page): its DATA statement is 37# commented out, so the Fortran's CC is undefined; 0.0 stands for it here 38# (TKHS is then zero, and the division by it raises ZeroDivisionError 39# where the Fortran would give an infinity). 40CC = 0.0 41 42 43def _interp4(f: np.ndarray, faci: np.ndarray, i: int, j0: int) -> float: 44 """4-point Lagrangian interpolation in the logarithm of f, from shells 45 j0..j0+3 with the weights FACI(1..4, I).""" 46 47 return math.exp(math.log(f[j0]) * faci[1, i] 48 + math.log(f[j0 + 1]) * faci[2, i] 49 + math.log(f[j0 + 2]) * faci[3, i] 50 + math.log(f[j0 + 3]) * faci[4, i]) 51 52 53def vcirc(hr: np.ndarray, hru: np.ndarray, imin: int, imax: int, it: int, lcz: np.ndarray, m: int, omega: np.ndarray, 54 ldo: bool, qwlnr: np.ndarray, qwrmax: np.ndarray, hv: np.ndarray, cltot: float, dt: float, hp: np.ndarray, 55 faci: np.ndarray) -> bool: 56 """Estimate the velocities of meridional circulation, the GSF 57 instability and secular and dynamical shear at each radiative 58 interface imin..imax (between shells i-1 and i), inhibited by mu 59 gradients and scaled by the user's factors; their sum hv sets the 60 diffusion coefficient for angular momentum and mixing. 61 62 Translates: vcirc.f#VCIRC (7be4133) 63 64 Arguments: log r hr and r hru, the range imin..imax, the iteration it, 65 the convective flags lcz, the number of shells m, omega, ldo (set), 66 the shear qwlnr and its dynamical limit qwrmax, hv (filled in place), 67 the total luminosity cltot, the time step dt, log P hp, and the 68 interpolation weights faci[1..4, i]. The velocities also go to 69 temp2_mod, difad2_mod, difad3_mod and quadd_mod.rat. 70 71 Returns ldo: True if any hv is positive. 72 73 """ 74 75 n = nmlpar_mod 76 c = constants_mod 77 cln = c.cln 78 t2 = temp2_mod 79 d2 = difad2_mod 80 d3 = difad3_mod 81 o2 = oldrot2_mod 82 tk = taukh_mod 83 vf = vfact_mod 84 iv = intvar_mod 85 rat = quadd_mod.rat 86 quad = quadru_mod.quad 87 gg = quadru_mod.gg 88 amum = mdphy_mod.amum 89 viscm = mdphy_mod.viscm 90 thdifm = mdphy_mod.thdifm 91 ves, vgsf, vss, vmu = t2.ves, t2.vgsf, t2.vss, t2.vmu 92 rm = iv.rm 93 94 ldo = False 95 96 # Store the old velocity estimates for later use (Fortran lines 136–149). 97 if it > 1: 98 for i in range(imin, imax + 1): # label 5 99 vmup[i] = vmu[i] 100 t2.vgsf0[i] = vgsf[i] 101 t2.ves0[i] = ves[i] 102 t2.vss0[i] = vss[i] 103 if n.ldifad: 104 for i in range(imin, imax + 1): 105 d2.vesa0[i] = d2.vesa[i] 106 d2.vesd0[i] = d2.vesd[i] 107 108 # Clear the velocities (Fortran lines 150–156). 109 for i in range(1, m + 1): # label 10 110 hv[i] = 0.0 111 ves[i] = 0.0 112 vgsf[i] = 0.0 113 vss[i] = 0.0 114 vmu[i] = 0.0 115 116 # Viscosity and thermal diffusivity at the interfaces, by 4-point 117 # Lagrangian interpolation in their logarithms (Fortran lines 163–191). 118 # First interface. 119 if imin < 3: 120 viscmi[2] = _interp4(viscm, faci, 2, 1) 121 thdifmi[2] = _interp4(thdifm, faci, 2, 1) 122 i0 = 3 123 else: 124 i0 = imin 125 # Last interface. 126 if imax == m: 127 viscmi[m] = _interp4(viscm, faci, m, m - 3) 128 thdifmi[m] = _interp4(thdifm, faci, m, m - 3) 129 i1 = m - 1 130 else: 131 i1 = imax 132 # General case. 133 for i in range(i0, i1 + 1): # label 20 134 viscmi[i] = _interp4(viscm, faci, i, i - 2) 135 thdifmi[i] = _interp4(thdifm, faci, i, i - 2) 136 137 # Linear interpolation for omega and mu (Fortran lines 193–196). 138 for i in range(2, m + 1): 139 amumi[i] = 0.5 * (amum[i] + amum[i - 1]) 140 wm[i] = 0.5 * (omega[i] + omega[i - 1]) 141 142 # Meridional circulation (Fortran lines 202–275). 143 if n.ies == 2: 144 # Kippenhahn & Moellenhoff (1974), restored (MHP 4/26). 145 # Suspected bug (see the page): with LDIFAD, VESA, VESD, VESD2, 146 # VESD3, THN and EQQCOD keep the previous call's values. 147 for i in range(imin, imax + 1): # label 30 148 if lcz[i] and lcz[i - 1]: 149 continue 150 dr = hru[i] - hru[i - 1] 151 dv = ((omega[i] * hru[i])**2 152 - (omega[i - 1] * hru[i - 1])**2) / dr 153 gmid = 0.5 * (gg[i] + gg[i - 1]) 154 phistd = 2.0 * c.cc23 * wm[i]**2 * rm[i] / gmid 155 phis2 = c.cc13 * wm[i]**2 * (hru[i]**2 / gg[i] 156 - hru[i - 1]**2 / gg[i - 1]) / dr 157 phi2 = (quad[i - 1] / gg[i - 1] - quad[i] / gg[i]) / dr 158 rat[i] = (phis2 + phi2) / phistd 159 # Question (see the page): overwritten at once, dropping the 160 # quadrupole term. 161 rat[i] = phis2 / phistd 162 ves[i] = abs(vf.fact1[i] * (vf.fact2[i] * rat[i] * wm[i]**2 163 + tk.fact6[i] * dv)) 164 else: 165 # Zahn (1991), with the quadrupole term of Zahn (1992). 166 for i in range(imin, imax + 1): # label 31 167 if lcz[i] and lcz[i - 1]: 168 continue 169 ves[i] = wm[i]**2 * (tk.fes1[i] + wm[i]**2 * tk.fes2[i]) 170 dr = hru[i] - hru[i - 1] 171 gmid = 0.5 * (gg[i] + gg[i - 1]) 172 phistd = 2.0 * c.cc23 * wm[i]**2 * rm[i] / gmid 173 phis2 = c.cc13 * wm[i]**2 * (hru[i]**2 / gg[i] 174 - hru[i - 1]**2 / gg[i - 1]) / dr 175 phi2 = (quad[i - 1] / gg[i - 1] - quad[i] / gg[i]) / dr 176 qw = wm[i] * (omega[i] - omega[i - 1]) / dr 177 rat[i] = (phis2 + phi2) / phistd 178 if rat[i] < 0.0: 179 print(f"{i:5d} RAT {rat[i]:10.3E} POT SPH,CYL,QUA" 180 f"{phistd:10.3E}{phis2:10.3E}{phi2:10.3E}\n" 181 f" QUA,G,R {quad[i]:12.3E}{quad[i - 1]:12.3E}" 182 f"{gg[i]:12.3E}{gg[i - 1]:12.3E}{hru[i]:12.3E}" 183 f"{hru[i - 1]:12.3E}") # FORMAT 303 184 if not n.ldifad: 185 ves[i] = abs(rat[i] * ves[i] - tk.fes3[i] * qw) 186 else: 187 # The d theta/dt term, and the factor 1/5 (MHP 06/02, 05/02). 188 o2.thn[i] = o2.festime[i] * (o2.thm[i] * qw - o2.thp[i]) / dt 189 d2.vesa[i] = 0.2 * (rat[i] * ves[i] + o2.thn[i]) 190 q1 = vf.fact7[i] * wm[i]**2 191 q2 = vf.fact8[i] * qw 192 q0 = tk.fes3[i] * wm[i]**2 193 d2.vesd[i] = 0.2 * (q0 + q1 + q2) 194 # Second- and third-order terms. 195 d3.vesd2[i] = 0.2 * d3.facd2[i] * wm[i] 196 d3.vesd3[i] = 0.2 * d3.facd3[i] * wm[i] 197 ves[i] = (rat[i] * ves[i] + o2.thn[i] 198 + (tk.fes3[i] + vf.fact7[i]) * qw 199 + vf.fact8[i] * ((omega[i] - omega[i - 1]) / dr)**2) 200 egridder_mod.eqqcod[i] = egridder_mod.qqcod[i] * wm[i] 201 202 # Inhibition by mu gradients (Fortran lines 280–337). 203 if n.imu == 3: 204 # Kippenhahn & Moellenhoff (1974): fictitious mu currents that oppose 205 # the circulation, restored (MHP 4/26). 206 for i in range(imin, imax + 1): 207 if lcz[i] and lcz[i - 1]: 208 continue # GOTO 32 209 # Square root of tKH* (Kippenhahn, IAU 66, p. 23); the factor 210 # of g is in FACT4. 211 tkhs = math.sqrt(viscmi[i] * CC / thdifmi[i]) 212 qmu = (amum[i] - amum[i - 1]) / (hru[i] - hru[i - 1]) 213 vmu[i] = abs(n.fmu * vf.fact1[i] * vf.fact4[i] * amumi[i] * qmu 214 / tkhs) 215 vmu2[i] = n.fmu * iv.hgm[i] * abs(qmu) / amumi[i] / wm[i]**2 216 else: 217 # Zahn & Maeder (1998): (del - del_ad) replaced by 218 # (del + del_mu - del_ad). 219 # Question (see the page): convective interfaces are not skipped. 220 for i in range(imin, imax + 1): 221 qmu = math.log(amum[i]) - math.log(amum[i - 1]) 222 qp = cln * (hp[i] - hp[i - 1]) 223 ddel = o2.ddelm[i] + qmu / qp 224 ddtest = max(o2.ddelm[i], 1.0e-3) 225 ddtest2 = max(ddel, 1.0e-3) 226 qqq = ddtest / ddtest2 227 ves[i] = ves[i] * qqq 228 if n.ldifad: 229 d2.vesa[i] = d2.vesa[i] * qqq 230 d2.vesd[i] = d2.vesd[i] * qqq 231 d3.vesd2[i] = d3.vesd2[i] * qqq 232 d3.vesd3[i] = d3.vesd3[i] * qqq 233 q1 = vf.fact7[i] * wm[i]**2 * qqq 234 q0 = tk.fes3[i] * wm[i]**2 * qqq 235 dr = hru[i] - hru[i - 1] 236 v2 = 0.2 * (q0 + q1) * (omega[i] - omega[i - 1]) / dr / wm[i] # noqa: F841 (unused in the Fortran too) 237 # Suspected bug (see the page): FESTIME's factor compounds over 238 # SECULR's iterations; FV0 is read nowhere. 239 o2.festime[i] = o2.festime[i] * qqq 240 d3.fv0[i] = d3.fv0[i] * qqq 241 242 # GSF instability: the Kippenhahn (1980) estimate. IGSF = 0: inhibited 243 # by viscosity and mu gradients; 1: not inhibited; 2: by viscosity only 244 # (Fortran lines 351–407). 245 for i in range(imin, imax + 1): # label 40 246 if lcz[i] and lcz[i - 1]: 247 continue 248 if n.igsf == 2 or n.igsf == 0: 249 qwrmx = 2.0 * math.sqrt(viscmi[i] / thdifmi[i]) * qwrmax[i] 250 if abs(qwlnr[i]) < qwrmx: 251 vgsf[i] = 0.0 252 continue 253 fxx = math.sqrt((abs(qwlnr[i]) - qwrmx) / qwrmx) 254 else: 255 fxx = 1.0 256 dr = hru[i] - hru[i - 1] 257 if n.igsf == 0: 258 qwrmx = 2.0 * math.sqrt(iv.hgm[i] * abs(amum[i] - amum[i - 1]) 259 / dr / amumi[i]) 260 if abs(qwlnr[i]) < qwrmx: 261 vgsf[i] = 0.0 262 continue 263 # Tests to avoid dividing by zero (MHP 05/02). 264 qq = 1.0e-2 * qwlnr[i] 265 if abs(qwrmx) < abs(qq): 266 fxx2 = 1.0 267 elif abs(qwrmx) < 1.0e-20: 268 fxx2 = 1.0 269 else: 270 fxx2 = math.sqrt((abs(qwlnr[i]) - qwrmx) / qwrmx) 271 fxx = min(fxx, fxx2) 272 rmid = 0.5 * (hru[i] + hru[i - 1]) 273 dlnwdr = abs(math.log(omega[i]) - math.log(omega[i - 1])) / dr 274 # If d ln(j)/dr < 0, use the larger of the two criteria; otherwise 275 # only d omega/dz nonzero. The factor 1/r: all the velocities are 276 # later multiplied by r. 277 if omega[i] * hru[i]**2 < omega[i - 1] * hru[i - 1]**2: 278 dlnjmdr = abs(2.0 / rmid + (math.log(omega[i]) 279 - math.log(omega[i - 1])) / dr) 280 fx = max(2.0 * dlnjmdr, 0.25 * dlnwdr) 281 vgsf[i] = tk.fgsfk[i] * fx * thdifmi[i] * wm[i]**2 / rmid 282 else: 283 vgsf[i] = 0.25 * tk.fgsfk[i] * thdifmi[i] * dlnwdr * wm[i]**2 \ 284 / rmid 285 vgsf[i] = abs(fxx * vgsf[i]) 286 287 # Secular and dynamical shear (Endal & Sofia, paper II) (Fortran lines 288 # 441–532). 289 for i in range(imin, imax + 1): # label 60 290 if abs(qwlnr[i]) > qwrmax[i]: 291 # Dynamical shear: a (large) velocity, and the GSF velocity 292 # rescaled to the gradient at the limit. 293 qwr = abs(qwlnr[i]) 294 vss[i] = 8.0 / 4.5e1 * thdifmi[i] * (qwr / qwrmax[i])**2 / rm[i] 295 dr = hru[i] - hru[i - 1] 296 rmid = 0.5 * (hru[i] + hru[i - 1]) 297 if omega[i] < omega[i - 1]: 298 wmid = omega[i] + qwrmax[i] * cln * (hr[i] - math.log10(rmid)) # noqa: F841 (unused in the Fortran too) 299 else: 300 wmid = omega[i] - qwrmax[i] * cln * (hr[i] - math.log10(rmid)) # noqa: F841 (unused in the Fortran too) 301 # WMIN: the OMEGA(I-1) at which the centred gradient QWLNR 302 # equals QWRMAX, between OMEGA(I) and OMEGA(I-1) (since 7be4133). 303 wmin = omega[i] + (omega[i - 1] - omega[i]) * qwrmax[i] / qwr 304 # DLNWDR0 > 0: the branch runs only if the two omegas differ. 305 dlnwdr0 = abs(math.log(omega[i]) - math.log(omega[i - 1])) / dr 306 dlnjmdr0 = abs(2.0 / rmid + (math.log(omega[i]) 307 - math.log(omega[i - 1])) / dr) 308 dlnjmdr = abs(2.0 / rmid + (math.log(omega[i]) 309 - math.log(wmin)) / dr) 310 dlnwdr = abs(math.log(omega[i]) - math.log(wmin)) / dr 311 # The branch from the actual profile (since 7be4133). 312 if omega[i] * hru[i]**2 < omega[i - 1] * hru[i - 1]**2: 313 fx = max(2.0 * dlnjmdr, 0.25 * dlnwdr) 314 fx2 = max(2.0 * dlnjmdr0, 0.25 * dlnwdr0) 315 vgsf[i] = vgsf[i] * fx / fx2 316 else: 317 vgsf[i] = vgsf[i] * dlnwdr / dlnwdr0 318 # To unit 6 at every such interface, at every call. 319 print(f" DYNAMICAL SHEAR-SHELL{i:5d} WTOP{omega[i]:11.3E}" 320 f" WBOT{omega[i - 1]:11.3E} LIMIT{wmin:11.3E}") # 9911 321 continue 322 # Secular shear: the maximum gradient it allows, with the square 323 # root of the Prandtl number; velocity from Zahn 1991 (MHP 3/92). 324 qwrmx = math.sqrt(viscmi[i] / thdifmi[i] * 1.25e-1 * n.rcrit) \ 325 * qwrmax[i] 326 if abs(qwlnr[i]) > qwrmx: 327 if abs((amum[i] - amum[i - 1]) / amumi[i]) < 1.0e-10: 328 qwrmx2 = 0.0 329 qwr = abs(qwlnr[i]) - qwrmx 330 vss[i] = 8.0 / 4.5e1 * thdifmi[i] * (qwr / qwrmax[i])**2 \ 331 / rm[i] 332 else: 333 # The limit set by the mu gradient; the larger one is used. 334 qwrmx2 = 2.0 * math.sqrt(max(1.0e-20, vf.fact5[i] * abs( 335 (amum[i] - amum[i - 1]) / amumi[i]))) 336 if abs(qwlnr[i]) > qwrmx2: 337 qwrmx = max(qwrmx2, qwrmx) 338 qwr = abs(qwlnr[i]) - qwrmx 339 vss[i] = 8.0 / 4.5e1 * thdifmi[i] \ 340 * (qwr / qwrmax[i])**2 / rm[i] 341 342 # The user's scale factors (and, with IMU = 3, the mu currents) 343 # (Fortran lines 540–599). 344 if n.imu == 3: 345 if it > 1: 346 for i in range(1, m + 1): 347 vmu[i] = 0.5 * (vmu[i] + vmup[i]) 348 for i in range(imin, imax + 1): 349 # VMU2 is stale at convective interfaces, where the velocities 350 # it divides are zero. 351 fm = 1.0 + vmu2[i] 352 rmid = 0.5 * (hru[i] + hru[i - 1]) 353 fcc = math.sqrt(n.fc * n.fesc) 354 vest = max(0.0, n.fes * ves[i] 355 - math.sqrt(n.fes * ves[i] * rmid) * fcc * vmu[i]) 356 if vest <= 0.0: 357 ves[i] = abs(n.fes * ves[i] / fm) 358 else: 359 ves[i] = vest 360 fcc = math.sqrt(n.fc * n.fgsfc) 361 if vgsf[i] > 0.0: 362 vgsft = max(0.0, n.fgsf * vgsf[i] 363 - math.sqrt(n.fgsf * vgsf[i] * rmid) * fcc 364 * vmu[i]) 365 else: 366 vgsft = 0.0 367 if vgsft <= 0.0: 368 # Suspected bug (see the page): FES, not FGSF. 369 vgsf[i] = abs(n.fes * vgsf[i] / fm) 370 else: 371 vgsf[i] = vgsft 372 vss[i] = max(0.0, n.fss * vss[i]) 373 else: 374 for i in range(imin, imax + 1): 375 ves[i] = abs(n.fes * ves[i]) 376 vgsf[i] = n.fgsf * vgsf[i] 377 vss[i] = n.fss * vss[i] 378 379 # Repeat for the diffusive and advective parts (MHP 11/94) (Fortran 380 # lines 602–607). 381 # Question (see the page): VES is made positive, VESA keeps its sign. 382 if n.ldifad: 383 for i in range(imin, imax + 1): 384 d2.vesd[i] = abs(n.fes * d2.vesd[i]) 385 d2.vesa[i] = n.fes * d2.vesa[i] 386 387 # The GSF and shear velocities, times r, join the diffusive part (MHP 388 # 05/02) (Fortran lines 619–647). 389 if n.ldifad: 390 for i in range(imin, imax + 1): 391 d2.vesd[i] = d2.vesd[i] + rm[i] * (abs(vgsf[i]) + abs(vss[i])) 392 d3.dss[i] = rm[i] * abs(vss[i]) 393 d3.dgsf[i] = rm[i] * abs(vgsf[i]) 394 395 # Average the previous and new estimates after the first iteration 396 # (Fortran lines 649–662). 397 if it > 1: 398 for i in range(imin, imax + 1): # label 70 399 vgsf[i] = 0.5 * (vgsf[i] + t2.vgsf0[i]) 400 ves[i] = 0.5 * (ves[i] + t2.ves0[i]) 401 vss[i] = 0.5 * (vss[i] + t2.vss0[i]) 402 if n.ldifad: 403 for i in range(imin, imax + 1): 404 d2.vesa[i] = 0.5 * (d2.vesa[i] + d2.vesa0[i]) 405 d2.vesd[i] = 0.5 * (d2.vesd[i] + d2.vesd0[i]) 406 407 # The total velocity (Fortran lines 663–667). 408 for i in range(imin, imax + 1): # label 80 409 hv[i] = vgsf[i] + ves[i] + vss[i] 410 if hv[i] < 1.0e-20: 411 hv[i] = 0.0 412 if hv[i] > 0.0: 413 ldo = True 414 415 # No transport in the core (L(r)/L < 0.99), to avoid numerical problems in 416 # the H-burning shell (9/93) (Fortran lines 671–687). 417 # Question (see the page): the threshold is hard-coded. 418 if n.lnoj: 419 for i in range(imin, imax + 1): 420 dell = iv.alm[i] / cltot 421 if dell < 9.9e-1: 422 hv[i] = 0.0 423 vgsf[i] = 0.0 424 ves[i] = 0.0 425 vss[i] = 0.0 426 d2.vesa[i] = 0.0 427 d2.vesd[i] = 0.0 428 else: 429 break # GOTO 81 430 431 return ldo

vcirc.f

Purpose

How it fits

The method

Design

From YREC's documentation

Suspected bugs, checked by running the code

1. (cannot be reached) The $\mu$ factor compounds on FESTIME.

2. (changes results, small) Checked by KC (Nov 10, 2026) CC is never set.

3. (never occurred) FES for the GSF velocity.

4. (cannot be reached) IES = 2 with LDIFAD reuses the previous call's parts.

Definitions

subroutineVCIRC

Change log

1. (cannot be reached) The $\mu$ factor compounds on `FESTIME`.

2. (changes results, small) Checked by KC (Nov 10, 2026) `CC` is never set.

3. (never occurred) `FES` for the GSF velocity.

4. (cannot be reached) `IES = 2` with `LDIFAD` reuses the previous call's parts.