wind.f

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

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

Variables on this page (17)

c4pil, cgl, cln, clsunl, constfactor, crsun, csigl, exm, exmd, exr, exw, ishort, json, linstb, taucz, wmax, wnew

Purpose

WIND removes angular momentum from the surface convection zone by a magnetized wind, "a Weber–Davis magnetic wind model" (its comment), over one time step. It is the older loss law; MWIND calls it when LMWIND is off, and otherwise uses its own law.

How it fits

Inputs: $\log L$ (BL), the time step DELTS (s), the convection zone's mass range (HSBOT, HSTOP) and shells (JSTART, JEND), the switch LJDOT, the surface rate OMEGAS, the mass SMASS (in solar masses, probably), $\log T_\mathrm{eff}$, the zone's moment of inertia HICZ. Through modules: the law's constant and exponents (NMLPAR_MOD, formerly COMMON/CWIND/: CONSTFACTOR, EXMD, EXW, EXR, EXM, the saturation WMAX, set by PARMIN); the overturn time TAUCZ (OVRTRN_MOD, from GETTAU); LINSTB (NMLPAR_MOD, formerly COMMON/ROT/).
Who calls it: only MWIND.
Output: HJM, each shell's specific angular momentum in the zone, reduced.
Modules: NMLPAR_MOD, OVRTRN_MOD, CONSTANTS_MOD, DIMENSIONS_MOD, IOUNITS_MOD (for ISHORT).

The method

Saturation. $\Omega_\mathrm{sat}$ = WMAX; but with angular momentum transport (LINSTB) and WMAX > 1, WMAX is read as $\Omega_\mathrm{sat}\tau_\mathrm{cz}$ for the Sun, and $\Omega_\mathrm{sat} = $ WMAX$/\tau_\mathrm{cz}$ (MHP 3/09); the run stops if $\tau_\mathrm{cz} \le 1$, with exit code 1 and the message also in the .short file.

Question 1 (lines 64–80 ↗). The saturation's units. WMAX > 1 switches the meaning of WMAX (rad/s, or rad/s × s). A rate above 1 rad/s is unphysical, so the test works as a switch, but the Physics namelist page describes WMAX only in rad/s. Would a separate switch, or a sentence in the documentation, be clearer?
The loss law. With $R = \sqrt{L/(4\pi\sigma T_\mathrm{eff}^4)}$ and the fixed mass-loss rate $\dot M = 2\times10^{-14}\,M_\odot$/yr, $$I\,\frac{d\Omega}{dt} = -K\left(\frac{\dot M}{10^{-14}}\right)^{a}\left(\frac{R}{R_\odot}\right)^{c} M^{d}\;\Omega\,\min(\Omega, \Omega_\mathrm{sat})^{\,b-1},$$ with $K$ = CONSTFACTOR, $a$ = EXMD, $b$ = EXW, $c$ = EXR, $d$ = EXM.

Question 2 (line 87 ↗; lines 97–99 ↗). A fixed mass-loss rate. DMDOT = 2.0D-14 is hard-coded, so the factor $(\dot M/10^{-14})^{a}$ is the constant $2^{a}$, whatever the model's own rate (DMDT0). Is this meant, as a part of the calibration through CONSTFACTOR?
Integration. If one explicit step would remove more than $\Omega$, the step is split into NSTEP substeps. In each, $\Omega_\mathrm{new}$ is found by iteration: the loss is evaluated at the mean of the start and the previous estimate, until the relative change is below $10^{-6}$ (at most 21 evaluations). The new rate is stored in WNEW, which is also a namelist variable (bug 1).

Question 3 (lines 117–134 ↗). No convergence check. After 21 evaluations without convergence, the last estimate is used without a message. Should the routine warn, or cut the substep?
Removal. $\Delta j = (\Omega_\mathrm{s} - \Omega_\mathrm{new})\,I_\mathrm{cz}/\Delta M_\mathrm{cz}$, with $\Delta M_\mathrm{cz}$ = HSTOP - HSBOT, is subtracted from every shell of the zone.

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

(line 122 ↗; line 133 ↗; line 142 ↗) 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 ↗). Each call of WIND therefore replaces the input with the last wind rate. The main program calls STARIN again for each run of its DO 200 NK = 1, NUMRUN loop (main.f, line 397 ↗), and STARIN applies WNEW to a model it reads from file (LFIRST, starin.f, lines 255–263 ↗), so a later run of the same job that reads its model from file with LWNEW would start from the previous run's last surface rate, not from the input. Only with LROT, LWNEW (off by default), the old wind law (LMWIND off, the default), and more than one run (NUMRUN > 1); mwind.f, suspected bug 1 is the same for the new law.

A possible fix: a local name for the iterate, here and in MWIND:

wind.f

@@ lines 24–31 @@
 C also reports to the .short file (issue 38).
       USE IOUNITS_MOD
       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,CGRAV,DELTS,DM,DMDOT,DT,DW,DWTEST,FJDOM,HICZ,HJM,HSBOT,
-     *       HSTOP,OMEGAS,RL,RTOT,SMASS,TEFFL,W,WNEW0,WS,WSAT
+     *       HSTOP,OMEGAS,RL,RTOT,SMASS,TEFFL,W,WNEW0,WNEWL,WS,WSAT
       INTEGER II,JEND,JSTART,K,N,NSTEP
       LOGICAL*4 LJDOT

@@ lines 119–136 @@
             WNEW0 = WS
     5       CONTINUE
             K = K + 1
-            WNEW = WS - (DT/HICZ)*CONSTFACTOR*(DMDOT/1.0D-14)**EXMD
+            WNEWL = WS - (DT/HICZ)*CONSTFACTOR*(DMDOT/1.0D-14)**EXMD
      *             *W*(RTOT/CRSUN)**EXR*SMASS**EXM
      *             *MIN(W,WSAT)**(EXW-1.0D0)
-            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*CONSTFACTOR*(DMDOT/1.0D-14)**EXMD*OMEGAS**(EXW-1.0D0)
 C    *           *(RTOT/CRSUN)**EXR*SMASS**EXM
@@ lines 139–145 @@
          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): Not reached: Every rotating case of YREC uses the new law (LMWIND on), so WIND never ran in the 33 runs; and, as for MWIND, YREC's multi-run cases with LWNEW start their later runs from the model in memory, so no run would use the overwritten value. The fix was not run. What would reach it: LROT with LMWIND off, and a later run that reads its model from file (LFIRST(k) true, k > 1) with LWNEW. No comment says the shared name is meant.

Definitions

subroutineWIND

lines 3–152 · GitHub

The older wind braking law: $I\,d\Omega/dt \propto \Omega\,\min(\Omega,\Omega_\mathrm{sat})^{b-1}$ with fixed exponents, integrated over the step and taken from the surface convection zone's specific angular momentum.

Uses modules CONSTANTS_MOD, DIMENSIONS_MOD, IOUNITS_MOD, NMLPAR_MOD, OVRTRN_MOD

Called by MWIND

Change log

The changes since rp_1

Changed since rp_1. Declarations only: COMMON blocks to modules (/CWIND/ and /ROT/ to NMLPAR_MOD, /OVRTRN/ to OVRTRN_MOD, the constants to CONSTANTS_MOD; the unused /DEUTER/ dropped), IMPLICIT NONE; former question 1 (now bug 1) stands, WNEW now being NMLPAR_MOD's.

Changed since a8a5950. The error stop when $\tau_\mathrm{cz}$ is not defined now also writes its message to the .short file (ISHORT, through IOUNITS_MOD), and ends with STOP 1 instead of a bare STOP, so that a failed run no longer exits with status 0 (KNOWN_ISSUES.md, issue 38).

1C$$$$$$ 2C MHP 10/02 UNUSED LFIRST REMOVED FROM CALL 3 SUBROUTINE WIND(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. 22C 10/30/26 MHP via Claude 23C code changed to USE IOUNITS_MOD for ISHORT, so the error stop below 24C also reports to the .short file (issue 38). 25 USE IOUNITS_MOD 26 IMPLICIT NONE 27 REAL*8 BL,CGRAV,DELTS,DM,DMDOT,DT,DW,DWTEST,FJDOM,HICZ,HJM,HSBOT, 28 * * HSTOP,OMEGAS,RL,RTOT,SMASS,TEFFL,W,WNEW0,WS,WSAT 29 INTEGER II,JEND,JSTART,K,N,NSTEP 30 LOGICAL*4 LJDOT 31 32 DIMENSION HJM(JSON) 33C 9/25/26 MHP via Claude 34C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/, 35C and COMMON/CONST3/ into a single COMMON/CONSTANTS/ 36C 9/24/26 MHP via Claude 37C code changed to remove dead, unused CMIXZ (never referenced) from 38C COMMON/CONST3/ 39C MHP 8/17 ADDED EXCEN, C_2 TO COMMON BLOCK FOR MATT ET AL. 2012 CENT. TERM 40C 10/12/26 MHP via Claude 41C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 42C namelist variables that the physics routines read from 27 COMMON 43C blocks which also held derived or working members. Those blocks keep 44C only their non-namelist members; five of them -- NEWMX, VMULT, 45C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 46C Derivation-only namelist variables from the same blocks went to 47C COMMON/NMLBYPASS/ instead, not here. 48C 49C Declared identically in all 121 files that hold any of these members. 50C Generated from one member list, so every copy is byte-identical; the 51C linker's 'size of symbol changed' warning is the check that they are. 52C 10/13/26 MHP via Claude: removed the unused COMMON/CWIND/ declaration 53C (no member referenced here; residue of the obsolete-option removals). 54C COMMON/CWIND/WMAX,EXMD,EXW,EXTAU,EXR,EXM,CONSTFACTOR,STRUCTFACTOR,LJDOT0 55C MHP 3/09 ADDED OPTION TO SCALE THE SATURATION RATE BY THE OVERTURN TIMESCALE 56C 10/13/26 MHP via Claude: removed the unused COMMON/DEUTER/ declaration 57C (no member referenced here; residue of the obsolete-option removals). 58 SAVE 59C IF DESIRED, REMOVE ANGULAR MOMENTUM FROM OUTER CONVECTION ZONE 60C USING A WEBER-DAVIS MAGNETIC WIND MODEL 61C MHP 3/09 IF WMAX > 1 THEN ASSUME THAT THE PARAMETER WMAX IS DEFINED BY 62C WMAX = WMAX(SUN)*TAUCZ(SUN) AND THE SATURATION THRESHOLD WSAT = WMAX/TAUCZ(STAR) 63C ONLY APPLY THIS IF ANGULAR MOMENTUM TRANSPORT ENABLED 64 IF(.NOT.LINSTB)THEN 65 WSAT = WMAX 66 ELSE IF(WMAX.GT.1.0D0)THEN 67 IF(TAUCZ.GT.1.0D0)THEN 68 WSAT = WMAX/TAUCZ 69c WRITE(*,912)WSAT,TAUCZ 70c 912 FORMAT('Omega sat, Tau',1p2e12.3) 71 ELSE 72 WRITE(*,911)WMAX,TAUCZ 73 WRITE(ISHORT,911)WMAX,TAUCZ 74 911 911 FORMAT('ERROR IN WIND - TAUCZ NOT DEFINED ',1P2E12.3,'STOPPED') 75C 10/30/26 MHP via Claude: bare STOP -> STOP 1, the run-time error code (issue 38). 76 STOP 1 77 ENDIF 78 ELSE 79 WSAT = WMAX 80 ENDIF 81 IF (LJDOT) THEN 82C FIND TOTAL RADIUS OF STAR. 83 CGRAV = DEXP(CLN*CGL) 84 RL=0.5D0*(BL+CLSUNL-C4PIL-CSIGL-4.D0*TEFFL) 85 RTOT = DEXP(CLN*RL) 86C DMDOT IS THE MASS LOSS RATE IN SOLAR MASSES PER YEAR. 87 DMDOT = 2.0D-14 88C DJ/DT = DT*CONSTFACTOR*(DMDOT/1.0D-14)**EXMD*OMEGA**EXW*(M/MSUN)**EXM 89C *(R/RSUN)**EXR 90C TEST : THE LOSS RATE DEPENDS ON OMEGA, AND FOR TIMESTESP THAT ARE 91C TOO LARGE, ROTATION RATES THAT ARE TOO HIGH, AND THIN SURFACE C.Z.S 92C NEGATIVE SURFACE ANGULAR VELOCITIES CAN BE PRODUCED. 93C TO AVOID THIS, CHECK THAT THE TIMESTEP IS SMALL ENOUGH TO ALLOW 94C A POSITIVE SOLUTION FOR OMEGA IN THE FIRST GUESS AT THE LOSS RATE. 95C IF NOT, USE A SERIES OF SMALL STEPS. 96C MHP 12/91 CAP LOSS RATE AT WSAT. 97 DWTEST = (DELTS/HICZ)*CONSTFACTOR*(DMDOT/1.0D-14)**EXMD 98 * * *OMEGAS*(RTOT/CRSUN)**EXR*SMASS**EXM 99 * * *MIN(OMEGAS,WSAT)**(EXW-1.0D0) 100 IF(DWTEST.GT.OMEGAS)THEN 101 NSTEP = INT(DWTEST/OMEGAS)+1 102 DT = DELTS/DFLOAT(NSTEP) 103 ELSE 104 NSTEP = 1 105 DT = DELTS 106 ENDIF 107C WRITE(*,3)NSTEP 108C 3 FORMAT(5X,I5) 109 WS = OMEGAS 110 DO 100 N = 1,NSTEP 111C THE CONSTANT AND EXPONENTS ARE SET IN PARMIN BASED ON THE INPUT 112C INDEX ALFA;SEE PARMIN FOR DETAILS ON THE DEPENDENCE OF EACH ON ALFA. 113C ITERATIVE SOLUTION : FOR FIRST GUESS, USE OMEGA=INITIAL OMEGA IN 114C COMPUTING LOSS RATE. TO COMPUTE SUBSEQUENT RATES, USE THE AVERAGE OF THE 115C STARTING OMEGA AND THE *PREVIOUS* ENDING OMEGA (I.E. AVERAGE EXPLICIT 116C AND IMPLICIT SOLUTIONS). 117 K = 0 118 W = WS 119 WNEW0 = WS 120 5 5 CONTINUE 121 K = K + 1 122 WNEW = WS - (DT/HICZ)*CONSTFACTOR*(DMDOT/1.0D-14)**EXMD 123 * * *W*(RTOT/CRSUN)**EXR*SMASS**EXM 124 * * *MIN(W,WSAT)**(EXW-1.0D0) 125 DW = 2.0D0*ABS((WNEW0-WNEW)/(WNEW0+WNEW)) 126C WRITE(*,4)WS,W,WNEW,DW,HICZ 127C 4 FORMAT(1X,1P5E14.6) 128 IF(DW.GT.1.0D-6)THEN 129 W = 0.5D0*(WS+WNEW) 130 WNEW0 = WNEW 131 IF(K.LE.20)GOTO 5 132 ENDIF 133 WS = WNEW 134 100 100 CONTINUE 135C CON = DELTS*CONSTFACTOR*(DMDOT/1.0D-14)**EXMD*OMEGAS**(EXW-1.0D0) 136C * *(RTOT/CRSUN)**EXR*SMASS**EXM 137C FJDOT = CON*OMEGAS/(1.0D0+(EXW*CON/HICZ)) 138C DM IS THE TOTAL MASS IN THE CONVECTION ZONE. 139 DM = HSTOP - HSBOT 140C FIND CHANGE IN ANGULAR MOMENTUM PER UNIT MASS AND SUBTRACT THIS 141C NUMBER FROM THE J/M OF EACH SHELL IN THE SURFACE CONVECTION ZONE. 142 FJDOM = (OMEGAS-WNEW)*HICZ/DM 143C WRITE(*,11)FJDOM,HJM(JSTART),HJM(JEND) 144C 11 FORMAT(5X,1P3E14.6) 145C FJDOM=FJDOT/DM 146C TAUJ=SJTOT/(FJDOT/DELTS)/CSECYR 147 DO 10 II = JSTART,JEND 148 HJM(II) = HJM(II) - FJDOM 149 10 10 CONTINUE 150 ENDIF 151 RETURN 152 END

1"""Python version of wind.f (YREC 7be4133): the older magnetized-wind loss 2of angular momentum from the surface convection zone.""" 3 4import math 5import numpy as np 6 7from pyrec import constants_mod 8from pyrec import iounits_mod 9from pyrec import nmlpar_mod 10from pyrec import ovrtrn_mod 11 12 13def wind(bl: float, delts: float, hsbot: float, hstop: float, jstart: int, jend: int, ljdot: bool, omegas: float, 14 smass: float, teffl: float, hicz: float, hjm: np.ndarray) -> None: 15 """Remove angular momentum from the surface convection zone over one 16 step by the Weber-Davis wind law 17 I dOmega/dt = -K (Mdot/1e-14)^a (R/Rsun)^c M^d Omega min(Omega, Omega_sat)^(b-1), 18 integrated in substeps by iteration, and subtract the loss of specific 19 angular momentum from every shell of the zone. 20 21 Translates: wind.f#WIND (7be4133) 22 23 Arguments: log L bl; the step delts (s); the zone's mass range hsbot, 24 hstop and shells jstart..jend; ljdot, the switch; the surface rate 25 omegas; smass; log Teff teffl; the zone's moment of inertia hicz; hjm, 26 the specific angular momentum (index 0 unused), changed in place. 27 None is returned (but nmlpar_mod.wnew is overwritten). 28 29 """ 30 31 c = constants_mod 32 p = nmlpar_mod 33 34 # The saturation rate, scaled by the overturn time if WMAX > 1 (Fortran 35 # lines 64–80). 36 if not p.linstb: 37 wsat = p.wmax 38 elif p.wmax > 1.0: 39 # Question (see the page): wmax > 1 switches its meaning. 40 if ovrtrn_mod.taucz > 1.0: 41 wsat = p.wmax / ovrtrn_mod.taucz 42 else: 43 msg = (f'ERROR IN WIND - TAUCZ NOT DEFINED ' 44 f'{p.wmax:12.3E}{ovrtrn_mod.taucz:12.3E}STOPPED') 45 print(msg) # WRITE(*,911) 46 iounits_mod.write(iounits_mod.ishort, msg) 47 raise SystemExit(1) # STOP 1, the run-time error code 48 else: 49 wsat = p.wmax 50 51 if ljdot: 52 # The radius, and the fixed mass-loss rate (Fortran lines 83–87). 53 cgrav = math.exp(c.cln * c.cgl) # noqa: F841 (unused, as in Fortran) 54 rl = 0.5 * (bl + c.clsunl - c.c4pil - c.csigl - 4.0 * teffl) 55 rtot = math.exp(c.cln * rl) 56 # Question (see the page): a fixed mass-loss rate. 57 dmdot = 2.0e-14 58 59 # Substeps, if one explicit step would remove more than Omega 60 # (Fortran lines 97–106). 61 dwtest = ((delts / hicz) * p.constfactor * (dmdot / 1.0e-14)**p.exmd 62 * omegas * (rtot / c.crsun)**p.exr * smass**p.exm 63 * min(omegas, wsat)**(p.exw - 1.0)) 64 if dwtest > omegas: 65 nstep = int(dwtest / omegas) + 1 66 dt = delts / float(nstep) 67 else: 68 nstep = 1 69 dt = delts 70 71 # The new rate in each substep, by iteration (Fortran lines 72 # 109–134). 73 ws = omegas 74 for n in range(1, nstep + 1): # loop 100 75 k = 0 76 w = ws 77 wnew0 = ws 78 while True: # label 5 79 k = k + 1 80 # Suspected bug (see the page): wnew is nmlpar_mod's 81 # solid-body rate, overwritten here. 82 p.wnew = (ws - (dt / hicz) * p.constfactor 83 * (dmdot / 1.0e-14)**p.exmd 84 * w * (rtot / c.crsun)**p.exr * smass**p.exm 85 * min(w, wsat)**(p.exw - 1.0)) 86 dw = 2.0 * abs((wnew0 - p.wnew) / (wnew0 + p.wnew)) 87 if dw > 1.0e-6: 88 w = 0.5 * (ws + p.wnew) 89 wnew0 = p.wnew 90 if k <= 20: 91 continue 92 # Question (see the page): no message when it fails. 93 break 94 ws = p.wnew 95 96 # Subtract the loss of specific angular momentum from every shell 97 # of the zone (Fortran lines 139–149). 98 dm = hstop - hsbot 99 fjdom = (omegas - p.wnew) * hicz / dm 100 for ii in range(jstart, jend + 1): # loop 10 101 hjm[ii] = hjm[ii] - fjdom

wind.f

Purpose

How it fits

The method

Suspected bugs, checked by running the code

1. (not reached) WNEW overwrites the namelist's solid-body rate.

Definitions

subroutineWIND

Change log

1. (not reached) `WNEW` overwrites the namelist's solid-body rate.