mwind.f

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

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

Variables on this page (22)

c4pil, c_2, cgl, cln, clsunl, cmsun, constfactor, crsun, csigl, excen, exw, fracstep, json, lbscale, lmwind, lrossby, pmmsoltau, pmmsolw, taucz, taucz0, wmax, wnew

Purpose

Removes angular momentum from the surface convection zone by a magnetized stellar wind over one time step. If LMWIND is false it hands over to the older wind law in WIND; otherwise it uses the wind law of Matt et al. (2012), as G. Somers' comments (3/17) call it, with a centrifugal correction from Matt et al. (2012, ApJ 754, L26; MHP 8/17) and optional Rossby scaling.

How it fits

Inputs: $\log L/L_\odot$ (BL), the step DELTS (s), the wind switch LJDOT (no angular momentum is removed when it is false), the masses at the top and bottom of the convection zone (HSTOP, HSBOT), its shells JSTART..JEND, the surface angular velocity OMEGAS, the mass in solar units SMASS (settled from the code, Nov 9, 2026: $f$ divides by SMASS and CMSUN, line 140 ↗, and mdot.f, line 262 ↗ updates it in solar masses), $\log T_\mathrm{eff}$ (TEFFL), the zone's moment of inertia HICZ, and the specific angular momentum HJM. Wind parameters (LMWIND, LROSSBY, LBSCALE, WMAX, EXW, CONSTFACTOR, EXCEN, PMMSOLW, PMMSOLTAU) are namelist inputs in NMLPAR_MOD, read by PARMIN; the constant C_2 is in CONSTANTS_MOD (set in PARMIN); the convective turnover times at the start and end of the step and the fraction of the step (TAUCZ0, TAUCZ, FRACSTEP) come from OVRTRN_MOD (AMCALC also reads the photospheric pressures there).
Who calls it: GETW and SECULR.
What it calls: WIND (old law) or AMCALC, which returns the structure factor STRUCTFACTOR as an output argument.
Output: HJM of the convection zone's shells, reduced.
Modules: DIMENSIONS_MOD (JSON), OVRTRN_MOD, NMLPAR_MOD and CONSTANTS_MOD. NMLPAR_MOD holds WNEW, the name the routine uses for the new angular velocity, so MWIND overwrites the input WNEW (bug 1).

The method

The switch. With LJDOT false, nothing is removed.
Rossby scaling (LROSSBY): the turnover time is interpolated in the step, $\tau = \tau_0 + f\,(\tau_1 - \tau_0)$ (FRACSTEP). With LBSCALE, the rotation rate used in the saturation term is scaled, $\Omega' = \Omega\,\tau/\tau_\odot$; otherwise the saturation threshold is scaled, $\Omega_\mathrm{sat}' = \Omega_\mathrm{sat}\,\tau_\odot/\tau$. LBSCALE acts only with LROSSBY.

Question 1 (lines 90–98 ↗). $\tau_\mathrm{cz}$ = 0. With LROSSBY, WSAT = WMAX*PMMSOLTAU/TCZ, and with LROSSBY and LBSCALE, WNOW = OMEGAS*TCZ/PMMSOLTAU. GETTAU sets TAUCZ = 0 first (gettau.f, line 79 ↗), and TAUINTNEW (tauintnew.f, line 279 ↗) and MASSLOSS (massloss.f, line 120 ↗) set it to 0 when they find no surface convection zone. MWIND's callers call it only with a surface convection zone (LCZSUR), but GETTAU and the callers judge the zone differently. Can MWIND be reached with TAUCZ and TAUCZ0 both 0, for example as a surface convection zone appears or vanishes, giving an infinite saturation threshold (a division by zero) or $\Omega' = 0$ (no torque)? Should it guard against it?
The structure factor from AMCALC: $S = (M/M_\odot)^{a_M}\,(R/R_\odot)^{a_R}\,(L/L_\odot)^{a_L}\,(P_\mathrm{phot}/P_{\mathrm{phot},\odot})^{a_P}\,(\tau/\tau_\odot)^{a_\tau}$.
The centrifugal factor, relative to the Sun, with $f = \Omega^2 R^3/(2GM)$ and $R$ from $L = 4\pi R^2\sigma T_\mathrm{eff}^4$: $$F_\mathrm{cen} = \left(\frac{c_2^2 + f_\odot}{c_2^2 + f}\right)^{e}\quad(c_2 = \texttt{C\_2},\ e = \texttt{EXCEN}),$$ where $f_\odot = \Omega_\odot^2R_\odot^3/(2GM_\odot)$ with $\Omega_\odot$ = PMMSOLW. The constant $c_2$ = 0.0506 (set in PARMIN) and $f$, which is $\tfrac12 f_\mathrm{M}^2$ for Matt et al.'s $f_\mathrm{M} = \Omega R^{3/2}/(GM)^{1/2}$, follow the centrifugal correction of Matt et al. (2012), written in PARMIN's comment as $\{K_2^2/(K_2^2 + 0.5\,\Omega^2R^3/GM)\}$ raised to a power.
The step. The torque law is $dJ/dt = -K\,S\,\Omega\,\min(\Omega', \Omega_\mathrm{sat})^{a_\Omega-1}F_\mathrm{cen}$: $\propto\Omega^{a_\Omega}$ unsaturated ($a_\Omega$ = 3; the YREC paper, §2.6.4, says of the Kawaler law that "the torque for slow rotators therefore scales as $\Omega^3$"; the spin-down this gives is asked below) and $\propto\Omega$ saturated; with LBSCALE, $\propto \Omega\,(\Omega\tau/\tau_\odot)^{a_\Omega-1}$, a dependence on $\mathrm{Ro}^{-1}$. PARMIN sets EXW $= 1 + a - 2am + 4mb$ from the law's exponents (PMMA, PMMM, PMMB), which is 3 for both V13 and K97 (with $a = 2$, $b = 1$ the terms in $m$ cancel); only CUS can give another value. The loop's comment, "THE CONSTANT AND EXPONENTS ARE SET IN PARMIN BASED ON THE INPUT INDEX ALFA" (lines 158–159 ↗), describes the old law: PARMIN uses ALFA only when LMWIND is off (parmin.f, lines 1024–1040 ↗), and sets EXW for this law from PMMA, PMMM, PMMB (parmin.f, line 1092 ↗) (settled from the code, Nov 9, 2026). The convection zone spins down as $$\Omega_\mathrm{new} = \Omega_s - \frac{\Delta t}{I_\mathrm{cz}}\;K\,S\;\Omega\;\min(\Omega', \Omega_\mathrm{sat})^{\,a_\Omega - 1}\;F_\mathrm{cen},$$ with $K$ = CONSTFACTOR, $a_\Omega$ = EXW. If the first guess of the change exceeds $\Omega_s$, the step is split into $\lfloor \Delta\Omega/\Omega_s \rfloor + 1$ substeps, so $\Omega$ stays positive. In each, $\Omega$ in the loss term is iterated as the mean of the start and the previous end value (the mean of the explicit and implicit solutions), until the relative change is below $10^{-6}$ or 20 iterations. The new rate is stored in WNEW, which is also a namelist variable (bug 1).

Question 2 (lines 164–200 ↗). Silent non-convergence. After 21 evaluations without convergence (lines 194–197 ↗), the substep continues with the last value, without a warning. Should the routine warn, or cut the substep?

Question 3 (lines 142–143 ↗; lines 186–187 ↗; parmin.f, line 1092 ↗). The spin-down law. With $a_\Omega$ = 3, $dJ/dt \propto \Omega^3$ unsaturated; at fixed moment of inertia that integrates to $\Omega \propto t^{-1/2}$, the Skumanich law. Is matching the Skumanich spin-down the reason both presets give EXW = 3, and does a changing moment of inertia on the pre-main sequence or the subgiant branch spoil the match?
The loss. $\Delta j = (\Omega_s - \Omega_\mathrm{new})\,I_\mathrm{cz}/\Delta M_\mathrm{cz}$ is subtracted from every shell of the zone.

From YREC's documentation

Physics Options (Angular momentum loss & winds): LJDOT0 turns on angular-momentum loss and LMWIND the wind model of Matt et al. (2012); AWIND chooses the law: V13 (van Saders & Pinsonneault 2013), K97 (the legacy Krishnamurthi et al. 1997) or CUS, custom, which alone reads PMMA…PMMM. FK scales the loss, usually to reproduce the Sun's rotation at its age. LROSSBY (default true) scales the loss with the Rossby number, and should go with a saturation threshold: WMAX (rad/s), replaced by WMAX_SUN $\times\,\Omega_\odot$ when WMAX_SUN is below 1000. PMMSOLTAU is the Sun's convective overturn time, and PMMSOLP $\log P$ at $T_\mathrm{eff}$ in a solar model, for the field's calibration.

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. (never occurred) `WNEW` overwrites the namelist's solid-body rate.

(line 186 ↗; line 199 ↗; line 208 ↗) As in wind.f, suspected bug 1: The iterated new rate is assigned to WNEW, which USE NMLPAR_MOD makes the namelist variable WNEW (nmlpar_mod.f, line 129 ↗), the rate at which STARIN starts a model in solid-body rotation when LWNEW is set (starin.f, lines 365–371 ↗). Every call of MWIND replaces it with the latest surface rate. STARIN is called again for each run in the main program's DO 200 NK = 1, NUMRUN loop (main.f, line 397 ↗), and applies WNEW to a model it reads from file (LFIRST, starin.f, lines 255–263 ↗), so a later run that reads its model from file with LWNEW would start from the wind's last value, not from the input. Only with LROT, LWNEW (off by default), LJDOT, and more than one run (NUMRUN > 1).

A possible fix:

mwind.f

@@ lines 20–28 @@
 C code changed to replace implicit typing with IMPLICIT NONE and
 C explicit declarations.
       IMPLICIT NONE
+C 11/09/26 KC via Claude
+C code changed to iterate on the local WNEWL instead of WNEW, which
+C is the namelist's solid-body rate in NMLPAR_MOD.
       REAL*8 BL,DELTS,DM,DT,DW,DWTEST,FCEN,FCORR,FJDOM,FSUN,HICZ,HJM,
      *       HSBOT,HSTOP,OMEGAS,RL,SMASS,STRUCTFACTOR,TCZ,TEFFL,W,WNEW0,
-     *       WNOW,WS,WSAT
+     *       WNEWL,WNOW,WS,WSAT
       INTEGER II,JEND,JSTART,K,N,NSTEP
       LOGICAL*4 LJDOT

@@ lines 183–202 @@
 C NOTE THAT THIS IS IMPLEMENTED HERE RELATIVE TO THE SUN
          FCORR = 0.5*W**2*EXP(CLN*(3.0*RL-CGL))/SMASS/CMSUN
          FCEN = ((C_2**2+FSUN)/(C_2**2+FCORR))**EXCEN
-         WNEW = WS - (DT/HICZ)*CONSTFACTOR*STRUCTFACTOR*W
+         WNEWL = WS - (DT/HICZ)*CONSTFACTOR*STRUCTFACTOR*W
      *          *MIN(WNOW,WSAT)**(EXW-1.0D0)*FCEN
 C         WNEW = WS - (DT/HICZ)*CONSTFACTOR*STRUCTFACTOR*W
 C     *          *MIN(W,WSAT)**(EXW-1.0D0)
 C G Somers END
-         DW = 2.0D0*ABS((WNEW0-WNEW)/(WNEW0+WNEW))
+         DW = 2.0D0*ABS((WNEW0-WNEWL)/(WNEW0+WNEWL))
 C       WRITE(*,4)WS,W,WNEW,DW,HICZ
 C    4 FORMAT(1X,1P5E14.6)
          IF(DW.GT.1.0D-6)THEN
-            W = 0.5D0*(WS+WNEW)
-            WNEW0 = WNEW
+            W = 0.5D0*(WS+WNEWL)
+            WNEW0 = WNEWL
             IF(K.LE.20)GOTO 5
          ENDIF
-         WS = WNEW
+         WS = WNEWL
   100 CONTINUE
 C     CON = DELTS*FACTOR*(DMDOT/1.0D-14)**EXMD*OMEGAS**(EXW-1.0D0)
 C    *        *(RTOT/CRSUN)**EXR*SMASS**EXM
@@ lines 205–211 @@
       DM = HSTOP - HSBOT
 C FIND CHANGE IN ANGULAR MOMENTUM PER UNIT MASS AND SUBTRACT THIS
 C NUMBER FROM THE J/M OF EACH SHELL IN THE SURFACE CONVECTION ZONE.
-      FJDOM = (OMEGAS-WNEW)*HICZ/DM
+      FJDOM = (OMEGAS-WNEWL)*HICZ/DM
 C      WRITE(*,11)FJDOM,HJM(JSTART),HJM(JEND)
 C   11 FORMAT(5X,1P3E14.6)
 C     FJDOM=FJDOT/DM

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 overwrite happens at every step with the new law: 360 times in testsuite/Test_m0p8feh+0p0GS98_rotp8p0d0p1mag1e4std_TAMS and 80 times in testsuite/Test_solarGS98_difrotmixper18p4d1p0 (WNEW from $3.9\times10^{-6}$ to $4.9\times10^{-5}$). No later run used it: The line added where STARIN applies LWNEW never saw a changed WNEW, since every multi-run case of YREC with LWNEW sets LFIRST false for its later runs, which then start from the model in memory. The fix was not run. What would reach it: a later run (LFIRST(k) true, k > 1) with LROT, LWNEW, and LMWIND. No comment says the shared name is meant.

Definitions

subroutineMWIND

lines 3–218 · GitHub

Uses modules CONSTANTS_MOD, DIMENSIONS_MOD, NMLPAR_MOD, OVRTRN_MOD

Called by GETW ×2, SECULR ×2

Flowchart of its calls, in order (2)

if (.NOT.LMWIND), lines 76–81

WIND 77

AMCALC 120

Change log

The changes since rp_1

Changed since rp_1. Declarations only, apart from one call: COMMON blocks to modules (/CWIND/ retired, its namelist members and those of /PMMWIND/ to NMLPAR_MOD, C_2 to CONSTANTS_MOD; /CONST/…/CONST3/ merged into CONSTANTS_MOD; /DEUTER/ and the character AWIND dropped), IMPLICIT NONE; AMCALC now returns STRUCTFACTOR as a fourth argument instead of setting it in /CWIND/. /ROT/ is no longer declared, but WNEW now comes from NMLPAR_MOD, so question 5 (now bug 1) stands.

Changed since a8a5950. With the Matt et al. law, the routine now returns without removing angular momentum when LJDOT is false, as WIND does (SECULR calls it whatever LJDOT0 is, so before, a run without the wind still lost angular momentum here); the extra field scaling by TCZ in the iteration now acts only with LROSSBY, as in the scaling above it (before, LBSCALE without LROSSBY used TCZ, which is set only under LROSSBY, and gave no torque), which settles former question 1; former questions 2–7 are now 1–6.

1C$$$$$$ 2C MHP 10/02 UNUSED LFIRST REMOVED FROM CALL 3 SUBROUTINE MWIND(BL,DELTS,HSBOT,HSTOP,JSTART,JEND,LJDOT,OMEGAS, 4C * SJTOT,SMASS,TEFFL,HICZ,HJM,LFIRST) ! KC 2025-05-31 5 * * SMASS,TEFFL,HICZ,HJM) 6C 10/15/26 MHP via Claude 7C code changed to take JSON from DIMENSIONS_MOD instead of declaring 8C PARAMETER(JSON=5000) here. 9 USE DIMENSIONS_MOD 10C 10/14/26 MHP via Claude 11C code changed to convert COMMON/OVRTRN/ into OVRTRN_MOD. 12 USE OVRTRN_MOD 13C 10/14/26 MHP via Claude 14C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD. 15 USE NMLPAR_MOD 16C 10/14/26 MHP via Claude 17C code changed to convert COMMON/CONSTANTS/ into CONSTANTS_MOD. 18 USE CONSTANTS_MOD 19C 10/20/26 MHP via Claude 20C code changed to replace implicit typing with IMPLICIT NONE and 21C explicit declarations. 22 IMPLICIT NONE 23 REAL*8 BL,DELTS,DM,DT,DW,DWTEST,FCEN,FCORR,FJDOM,FSUN,HICZ,HJM, 24 * * HSBOT,HSTOP,OMEGAS,RL,SMASS,STRUCTFACTOR,TCZ,TEFFL,W,WNEW0, 25 * * WNOW,WS,WSAT 26 INTEGER II,JEND,JSTART,K,N,NSTEP 27 LOGICAL*4 LJDOT 28 29 DIMENSION HJM(JSON) 30C 9/25/26 MHP via Claude 31C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/, 32C and COMMON/CONST3/ into a single COMMON/CONSTANTS/ 33C 9/24/26 MHP via Claude 34C code changed to remove dead, unused CMIXZ (never referenced) from 35C COMMON/CONST3/ 36C MHP 8/17 ADDED EXCEN, C_2 TO COMMON BLOCK FOR MATT ET AL. 2012 CENT. TERM 37C 10/12/26 MHP via Claude 38C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 39C namelist variables that the physics routines read from 27 COMMON 40C blocks which also held derived or working members. Those blocks keep 41C only their non-namelist members; five of them -- NEWMX, VMULT, 42C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 43C Derivation-only namelist variables from the same blocks went to 44C COMMON/NMLBYPASS/ instead, not here. 45C 46C Declared identically in all 121 files that hold any of these members. 47C Generated from one member list, so every copy is byte-identical; the 48C linker's 'size of symbol changed' warning is the check that they are. 49C 10/13/26 MHP via Claude 50C code changed to retire COMMON/CWIND/. Its namelist members moved to the 51C consolidated block earlier, leaving two. C_2 is a constant, set once in 52C parmin.f, so it is now the last member of COMMON/CONSTANTS/ -- appended 53C there rather than inserted, and added to every declaration so the block 54C keeps one size; parmin.f DATA-initializes CLSUN and CRSUN, so a short 55C declaration anywhere would be the OP95INTERP hazard again. It may 56C eventually become a namelist variable. STRUCTFACTOR, the loss-law 57C prefactor, is computed by AMCALC and read by its own two callers, 58C mwind.f and mcowind.f, so it is now an AMCALC output argument. 59C COMMON/CWIND/WMAX,EXMD,EXW,EXTAU,EXR,EXM,CONSTFACTOR,STRUCTFACTOR,LJDOT0 60C MHP 3/09 ADDED OPTION TO SCALE THE SATURATION RATE BY THE OVERTURN TIMESCALE 61C 10/13/26 MHP via Claude: removed the unused COMMON/DEUTER/ declaration 62C (no member referenced here; residue of the obsolete-option removals). 63C G Somers 3/17, ADDING NEW TAUCZ COMMON BLOCK 64C KC 2025-05-30 reordered common block elements 65C COMMON/PMMWIND/AWIND,LMWIND,LROSSBY,LBSCALE,PMMA,PMMB,PMMC,PMMD,PMMM, 66C * PMMJD,PMMMD,PMMSOLP,PMMSOLW,PMMSOLTAU 67C 10/12/26 MHP via Claude 68C code changed to remove the trailing CHARACTER member AWIND from 69C COMMON/PMMWIND/. Character data does not belong in a COMMON block, and 70C this file never referenced AWIND. It is now a local of parmin.f, the 71C only routine that reads it. Being last, its removal moves no other 72C member, and the block shrinks identically in all five files. 73 SAVE 74C IF DESIRED, REMOVE ANGULAR MOMENTUM FROM OUTER CONVECTION ZONE 75C USING A WEBER-DAVIS MAGNETIC WIND MODEL 76 IF(.NOT.LMWIND)THEN 77 CALL WIND(BL,DELTS,HSBOT,HSTOP,JSTART,JEND,LJDOT,OMEGAS, 78C * SJTOT,SMASS,TEFFL,HICZ,HJM) ! KC 2025-05-31 79 * * SMASS,TEFFL,HICZ,HJM) 80 GOTO 9999 81 ENDIF 82C 10/30/26 MHP via Claude 83C code changed to honour the wind switch, as WIND does. SECULR calls 84C MWIND whatever LJDOT0 is, so a run with LJDOT0 = F still lost J here. 85 IF(.NOT.LJDOT) GOTO 9999 86C 87C ADD ROSSBY SCALING IF DESIRED. 88 IF(LROSSBY)THEN 89C MHP 8/17 CORRECTED TAUCZ CALCULATION TO INTERPOLATE PROPERLY IN TIMESTEP 90 TCZ = TAUCZ0+FRACSTEP*(TAUCZ-TAUCZ0) 91 IF(LBSCALE)THEN 92C G Somers 8/17 CREATE ROTATION DUMMY VARIABLES. 93 WNOW = OMEGAS*TCZ/PMMSOLTAU 94 WSAT = WMAX 95 ELSE 96 WNOW = OMEGAS 97 WSAT = WMAX*PMMSOLTAU/TCZ 98 ENDIF 99C If not scaling, just set dummies to the original values. 100 ELSE 101 WNOW = OMEGAS 102 WSAT = WMAX 103 ENDIF 104C 105C G Somers 3/17, CHANGING WINDLAW TO NEW MATT+2012 METHOD. 106C 107C FIND TOTAL RADIUS OF STAR. 108C G Somers - CGRAV, RTOT not used, so blacking out. 109C CGRAV = DEXP(CLN*CGL) 110C RL=0.5D0*(BL+CLSUNL-C4PIL-CSIGL-4.D0*TEFFL) 111C RTOT = DEXP(CLN*RL) 112C DMDOT IS THE MASS LOSS RATE IN SOLAR MASSES PER YEAR. 113C G Somers - Don't need DMDOT anymore. 114C DMDOT = 2.0D-14 115C DJ/DT = DT*FACTOR*(DMDOT/1.0D-14)**EXMD*OMEGA**EXW*(M/MSUN)**EXM 116C *(R/RSUN)**EXR 117C 118C G Somers - New calculation. Call AMCALC to set "STRUCTFACTOR". 119C 120 CALL AMCALC(SMASS,BL,TEFFL,STRUCTFACTOR) 121C 122C TEST : THE LOSS RATE DEPENDS ON OMEGA, AND FOR TIMESTESP THAT ARE 123C TOO LARGE, ROTATION RATES THAT ARE TOO HIGH, AND THIN SURFACE C.Z.S 124C NEGATIVE SURFACE ANGULAR VELOCITIES CAN BE PRODUCED. 125C TO AVOID THIS, CHECK THAT THE TIMESTEP IS SMALL ENOUGH TO ALLOW 126C A POSITIVE SOLUTION FOR OMEGA IN THE FIRST GUESS AT THE LOSS RATE. 127C IF NOT, USE A SERIES OF SMALL STEPS. 128C MHP 12/91 CAP LOSS RATE AT WSAT. 129C 130C G Somers, Commenting out the old DWTEST, adding the PMM version. 131C DWTEST = (DELTS/HICZ)*FACTOR*(DMDOT/1.0D-14)**EXMD 132C * *OMEGAS*(RTOT/CRSUN)**EXR*SMASS**EXM 133C * *MIN(OMEGAS,WSAT)**(EXW-1.0D0) 134C MHP 8/17 ADDED CENTRIFUGAL REDUCTION TERM FROM MATT+2012 ApJ 754, L26 135C NOTE THAT THIS IS IMPLEMENTED HERE RELATIVE TO THE SUN 136C C_2 = 0.0506 137 FSUN = 0.5*PMMSOLW**2*CRSUN**3/EXP(CLN*CGL)/CMSUN 138C RADIUS 139 RL = 0.5D0*(BL+CLSUNL-C4PIL-CSIGL-4.D0*TEFFL) 140 FCORR = 0.5*OMEGAS**2*EXP(CLN*(3.0*RL-CGL))/SMASS/CMSUN 141 FCEN = ((C_2**2+FSUN)/(C_2**2+FCORR))**EXCEN 142 DWTEST = (DELTS/HICZ)*CONSTFACTOR*STRUCTFACTOR*OMEGAS 143 * * *MIN(WNOW,WSAT)**(EXW-1.0D0)*FCEN 144C DWTEST = (DELTS/HICZ)*CONSTFACTOR*STRUCTFACTOR*OMEGAS 145C * *MIN(OMEGAS,WSAT)**(EXW-1.0D0) 146C G Somers END 147 IF(DWTEST.GT.OMEGAS)THEN 148 NSTEP = INT(DWTEST/OMEGAS)+1 149 DT = DELTS/DFLOAT(NSTEP) 150 ELSE 151 NSTEP = 1 152 DT = DELTS 153 ENDIF 154C WRITE(*,3)NSTEP 155C 3 FORMAT(5X,I5) 156 WS = OMEGAS 157 DO 100 N = 1,NSTEP 158C THE CONSTANT AND EXPONENTS ARE SET IN PARMIN BASED ON THE INPUT 159C INDEX ALFA;SEE PARMIN FOR DETAILS ON THE DEPENDENCE OF EACH ON ALFA. 160C ITERATIVE SOLUTION : FOR FIRST GUESS, USE OMEGA=INITIAL OMEGA IN 161C COMPUTING LOSS RATE. TO COMPUTE SUBSEQUENT RATES, USE THE AVERAGE OF THE 162C STARTING OMEGA AND THE *PREVIOUS* ENDING OMEGA (I.E. AVERAGE EXPLICIT 163C AND IMPLICIT SOLUTIONS). 164 K = 0 165 W = WS 166 WNEW0 = WS 167 5 5 CONTINUE 168C G Somers 08/17 IF ADDING ADDITIONAL B SCALING, ADD ADDITIONAL TAUCZ TERM. 169C 10/30/26 MHP via Claude 170C code changed to scale by TCZ only under LROSSBY, as above. TCZ is set 171C only there, so LBSCALE without LROSSBY used TCZ = 0 and gave no torque. 172 IF(LROSSBY.AND.LBSCALE)THEN 173 WNOW = W*TCZ/PMMSOLTAU 174 ELSE 175 WNOW = W 176 ENDIF 177 K = K + 1 178C G Somers, Commenting out the old WNEW, adding the PMM version. 179C WNEW = WS - (DT/HICZ)*FACTOR*(DMDOT/1.0D-14)**EXMD 180C * *W*(RTOT/CRSUN)**EXR*SMASS**EXM 181C * *MIN(W,WSAT)**(EXW-1.0D0) 182C MHP 8/17 ADDED CENTRIFUGAL REDUCTION TERM FROM MATT+2012 ApJ 754, L26 183C NOTE THAT THIS IS IMPLEMENTED HERE RELATIVE TO THE SUN 184 FCORR = 0.5*W**2*EXP(CLN*(3.0*RL-CGL))/SMASS/CMSUN 185 FCEN = ((C_2**2+FSUN)/(C_2**2+FCORR))**EXCEN 186 WNEW = WS - (DT/HICZ)*CONSTFACTOR*STRUCTFACTOR*W 187 * * *MIN(WNOW,WSAT)**(EXW-1.0D0)*FCEN 188C WNEW = WS - (DT/HICZ)*CONSTFACTOR*STRUCTFACTOR*W 189C * *MIN(W,WSAT)**(EXW-1.0D0) 190C G Somers END 191 DW = 2.0D0*ABS((WNEW0-WNEW)/(WNEW0+WNEW)) 192C WRITE(*,4)WS,W,WNEW,DW,HICZ 193C 4 FORMAT(1X,1P5E14.6) 194 IF(DW.GT.1.0D-6)THEN 195 W = 0.5D0*(WS+WNEW) 196 WNEW0 = WNEW 197 IF(K.LE.20)GOTO 5 198 ENDIF 199 WS = WNEW 200 100 100 CONTINUE 201C CON = DELTS*FACTOR*(DMDOT/1.0D-14)**EXMD*OMEGAS**(EXW-1.0D0) 202C * *(RTOT/CRSUN)**EXR*SMASS**EXM 203C FJDOT = CON*OMEGAS/(1.0D0+(EXW*CON/HICZ)) 204C DM IS THE TOTAL MASS IN THE CONVECTION ZONE. 205 DM = HSTOP - HSBOT 206C FIND CHANGE IN ANGULAR MOMENTUM PER UNIT MASS AND SUBTRACT THIS 207C NUMBER FROM THE J/M OF EACH SHELL IN THE SURFACE CONVECTION ZONE. 208 FJDOM = (OMEGAS-WNEW)*HICZ/DM 209C WRITE(*,11)FJDOM,HJM(JSTART),HJM(JEND) 210C 11 FORMAT(5X,1P3E14.6) 211C FJDOM=FJDOT/DM 212C TAUJ=SJTOT/(FJDOT/DELTS)/CSECYR 213 DO 10 II = JSTART,JEND 214 HJM(II) = HJM(II) - FJDOM 215 10 10 CONTINUE 216 9999 9999 CONTINUE 217 RETURN 218 END

1"""Python version of mwind.f (YREC development, 7be4133).""" 2 3import math 4import numpy as np 5 6from pyrec import constants_mod 7from pyrec import nmlpar_mod 8from pyrec import ovrtrn_mod 9from pyrec.amcalc import amcalc 10from pyrec.wind import wind 11 12# Fortran SAVEs its locals. TCZ is set only when LROSSBY is true, and is 13# read only then. 14_tcz = 0.0 15 16 17def mwind(bl: float, delts: float, hsbot: float, hstop: float, jstart: int, jend: int, ljdot: bool, omegas: float, 18 smass: float, teffl: float, hicz: float, hjm: np.ndarray) -> None: 19 """Remove angular momentum from the surface convection zone by a wind. 20 21 Uses the magnetized-wind law of Matt et al. (2012) with their 22 centrifugal correction and optional Rossby scaling, over one time step 23 (split into substeps if needed); hands over to WIND (the older law) 24 if LMWIND is false. The loss of specific angular momentum is 25 subtracted from every shell of the zone. 26 27 Translates: mwind.f#MWIND (7be4133) 28 29 Parameters 30 ---------- 31 bl : float 32 log10 L/Lsun. 33 delts : float 34 Time step (s). 35 hsbot, hstop : float 36 Masses at the bottom and top of the convection zone. 37 jstart, jend : int 38 First and last shells of the zone. 39 ljdot : bool 40 The wind switch: no angular momentum is removed when it is false. 41 omegas : float 42 Surface angular velocity at the start of the step. 43 smass : float 44 Mass in solar units. 45 teffl : float 46 log10 Teff. 47 hicz : float 48 Moment of inertia of the zone. 49 hjm : ndarray, shape (JSON + 1,) 50 Specific angular momentum, changed in place for shells 51 ``jstart..jend``. 52 53 Notes 54 ----- 55 Sets ``nmlpar_mod.wnew`` (the new surface angular velocity), as the 56 Fortran does through the module. 57 58 """ 59 60 global _tcz 61 c = constants_mod 62 nml = nmlpar_mod 63 64 # If desired, remove angular momentum using a Weber-Davis wind model; 65 # otherwise use the old law (Fortran lines 76–81). 66 if not nml.lmwind: 67 wind(bl, delts, hsbot, hstop, jstart, jend, ljdot, omegas, 68 smass, teffl, hicz, hjm) 69 return # label 9999 70 71 # Honour the wind switch, as WIND does (Fortran line 85). 72 if not ljdot: 73 return # label 9999 74 75 # Rossby scaling, if desired; tau_cz interpolated within the step 76 # (Fortran lines 88–103). 77 # Question (see the page): tcz may be 0. 78 if nml.lrossby: 79 _tcz = (ovrtrn_mod.taucz0 80 + ovrtrn_mod.fracstep * (ovrtrn_mod.taucz - ovrtrn_mod.taucz0)) 81 if nml.lbscale: 82 wnow = omegas * _tcz / nml.pmmsoltau 83 wsat = nml.wmax 84 else: 85 wnow = omegas 86 wsat = nml.wmax * nml.pmmsoltau / _tcz 87 else: 88 wnow = omegas 89 wsat = nml.wmax 90 91 # The structure factor of the Matt et al. (2012) law (Fortran line 92 # 120). 93 structfactor = amcalc(smass, bl, teffl, 0.0) 94 95 # Centrifugal reduction term (Matt et al. 2012, ApJ 754, L26), 96 # relative to the Sun. C_2 = 0.0506 (Fortran lines 137–141). 97 fsun = (0.5 * nml.pmmsolw ** 2 * c.crsun ** 3 98 / math.exp(c.cln * c.cgl) / c.cmsun) 99 # Radius: 100 rl = 0.5 * (bl + c.clsunl - c.c4pil - c.csigl - 4.0 * teffl) 101 fcorr = (0.5 * omegas ** 2 * math.exp(c.cln * (3.0 * rl - c.cgl)) 102 / smass / c.cmsun) 103 fcen = ((c.c_2 ** 2 + fsun) / (c.c_2 ** 2 + fcorr)) ** nml.excen 104 105 # First guess of the change in Omega; if it exceeds Omega, use 106 # several smaller steps so that Omega stays positive (Fortran lines 107 # 142–153). 108 dwtest = ((delts / hicz) * nml.constfactor * structfactor * omegas 109 * min(wnow, wsat) ** (nml.exw - 1.0) * fcen) 110 if dwtest > omegas: 111 nstep = int(dwtest / omegas) + 1 112 dt = delts / float(nstep) 113 else: 114 nstep = 1 115 dt = delts 116 117 # Iterative solution in each substep: first guess uses the starting 118 # Omega; then the average of the starting and the previous ending 119 # Omega (Fortran lines 156–200). 120 ws = omegas 121 for n in range(1, nstep + 1): # label 100 122 k = 0 123 w = ws 124 wnew0 = ws 125 while True: # label 5 126 if nml.lrossby and nml.lbscale: 127 wnow = w * _tcz / nml.pmmsoltau 128 else: 129 wnow = w 130 k = k + 1 131 fcorr = (0.5 * w ** 2 * math.exp(c.cln * (3.0 * rl - c.cgl)) 132 / smass / c.cmsun) 133 fcen = ((c.c_2 ** 2 + fsun) / (c.c_2 ** 2 + fcorr)) ** nml.excen 134 # Suspected bug (see the page): WNEW is the namelist input in 135 # nmlpar_mod, overwritten here. 136 nml.wnew = (ws - (dt / hicz) * nml.constfactor * structfactor 137 * w * min(wnow, wsat) ** (nml.exw - 1.0) * fcen) 138 dw = 2.0 * abs((wnew0 - nml.wnew) / (wnew0 + nml.wnew)) 139 if dw > 1.0e-6: 140 w = 0.5 * (ws + nml.wnew) 141 wnew0 = nml.wnew 142 # Question (see the page): after 21 passes the loop ends 143 # without a warning. 144 if k <= 20: 145 continue 146 break 147 ws = nml.wnew 148 149 # DM is the total mass in the convection zone; subtract the change in 150 # angular momentum per unit mass from every shell of the zone (Fortran 151 # lines 205–215). 152 dm = hstop - hsbot 153 fjdom = (omegas - nml.wnew) * hicz / dm 154 for ii in range(jstart, jend + 1): # label 10 155 hjm[ii] = hjm[ii] - fjdom

mwind.f

Purpose

How it fits

The method

From YREC's documentation

Suspected bugs, checked by running the code

1. (never occurred) WNEW overwrites the namelist's solid-body rate.

Definitions

subroutineMWIND

Change log

1. (never occurred) `WNEW` overwrites the namelist's solid-body rate.