cowind.f

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

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

Variables on this page (13)

c4pil, cln, clsunl, constfactor, crsun, csigl, exm, exmd, exr, exw, ishort, taucz, wmax

Purpose

This file gives the change of the surface convection zone's angular velocity over a time step due to magnetized-wind angular-momentum loss, with the older "Kawaler-type" wind law and a saturation threshold. It is used when the newer Matt et al. law (LMWIND) is off.

How it fits

Inputs: the arguments BL ($\log L/L_\odot$), TEFFL, SMASS ($M/M_\odot$), DT (the time step), HICZ (the moment of inertia of the surface convection zone), OMEGAS (the surface $\Omega$ after the previous iteration), WOLD (at the start of the step), IT (the iteration); WMAX, EXMD, EXW, EXR, EXM, CONSTFACTOR in NMLPAR_MOD (set in PARMIN from the input index "ALFA," per the comment); TAUCZ (the convective overturn time) in OVRTRN_MOD; the solar constants in CONSTANTS_MOD.
Who calls it: MCOWIND, when LMWIND is false ("run old windlaw").
Outputs: WIND1, WIND2, the change of $\Omega$ over the step computed with $\Omega$ at its start and at its end. The comment says they enter the diffusion equations and are averaged: SECULR, which calls MCOWIND, passes them to DCOEFT or DADCOEFT, where the surface boundary takes $-\tfrac12(\mathtt{WIND1}+\mathtt{WIND2})$ times the moment of inertia.

The method

Saturation. If WMAX $> 1$, it is read as $\Omega_\mathrm{sat,\odot}\,\tau_\odot$ and $\Omega_\mathrm{sat} = \mathtt{WMAX}/\tau_\mathrm{cz}$ (MHP 3/09; the run stops if $\tau_\mathrm{cz}$ is not defined, with exit code 1 and the message also in the .short file); otherwise $\Omega_\mathrm{sat} = \mathtt{WMAX}$. Unlike WIND, the reinterpretation does not depend on LINSTB.

Question 1 (lines 80–96 ↗). The test WMAX > 1. A real $\Omega_\mathrm{sat}$ is far below 1 rad s$^{-1}$, so the test works as a switch; but the Physics namelist page describes WMAX only as an angular velocity in rad/s (and WMAX_SUN, which PARMIN converts for the new law only), not this second meaning. Should the documentation describe it, or a separate switch replace it?
The radius from $L$ and $T_\mathrm{eff}$: $\log R = \tfrac12(\log(L/L_\odot) + \log L_\odot - \log 4\pi - \log\sigma - 4\log T_\mathrm{eff})$, in cm.
The comment's law, $dJ/dt \propto K\,(\dot M/10^{-14})^{a}\,\Omega^{n}\,(M/M_\odot)^{c}(R/R_\odot)^{b}$, gives $$C = \frac{\Delta t}{I_\mathrm{CZ}}\,K\left(\frac{\dot M}{10^{-14}\,M_\odot\,\mathrm{yr}^{-1}}\right)^{a}\left(\frac{R}{R_\odot}\right)^{b}\left(\frac{M}{M_\odot}\right)^{c},$$ $$\Delta\Omega = C\,\min(\Omega, \Omega_\mathrm{sat})^{\,n-1}\,\Omega,$$ with $K$ = CONSTFACTOR, $a$ = EXMD, $b$ = EXR, $c$ = EXM, $n$ = EXW. So $\Delta\Omega \propto \Omega^n$ below saturation and $\propto \Omega\,\Omega_\mathrm{sat}^{n-1}$ above.
WIND1 uses WOLD. WIND2 uses OMEGAS; from the second iteration on, it is the mean of its previous value and the new one, which damps the iteration.

Question 2 (lines 111–115 ↗). The damping of WIND2. Averaging the previous WIND2 with the new value at each iteration means the converged WIND2 mixes values from several iterations, weighting the latest most. Is this only for stability, and does it converge to the value at the end of the step?

Design

The mass-loss rate is fixed: DMDOT = 2.0D-14 ($M_\odot$ yr$^{-1}$), so the factor $(\dot M/10^{-14})^a$ is the constant $2^a$.

Definitions

subroutineCOWIND

lines 2–118 · GitHub

Older wind law: $\Delta\Omega = (\Delta t/I_\mathrm{CZ})\,K\,(\dot M/10^{-14})^{a}(R/R_\odot)^{b}(M/M_\odot)^{c}\,\min(\Omega,\Omega_\mathrm{sat})^{n-1}\,\Omega$, at the start (WIND1) and end (WIND2) of the step; saturation optionally scaled by the overturn time.

Uses modules CONSTANTS_MOD, IOUNITS_MOD, NMLPAR_MOD, OVRTRN_MOD

Called by MCOWIND

Change log

The changes since rp_1

Changed since rp_1. Declarations only: COMMON blocks to modules, IMPLICIT NONE; the unused /DEUTER/ is gone, and the wind parameters this file does not use (EXTAU, EXL, EXPR, EXCEN, C_2) are now namelist variables in NMLPAR_MOD, shared by all its users, which settles the former question 3 on unused inputs.

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$$$$$$ 2 SUBROUTINE COWIND(BL,DT,HICZ,IT,OMEGAS,SMASS, 3 * * TEFFL,WOLD,WIND1,WIND2) 4C 10/14/26 MHP via Claude 5C code changed to convert COMMON/OVRTRN/ into OVRTRN_MOD. 6 USE OVRTRN_MOD 7C 10/14/26 MHP via Claude 8C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD. 9 USE NMLPAR_MOD 10C 10/14/26 MHP via Claude 11C code changed to convert COMMON/CONSTANTS/ into CONSTANTS_MOD. 12 USE CONSTANTS_MOD 13C 10/20/26 MHP via Claude 14C code changed to replace implicit typing with IMPLICIT NONE and 15C explicit declarations. 16C 10/30/26 MHP via Claude 17C code changed to USE IOUNITS_MOD for ISHORT, so the error stop below 18C also reports to the .short file (issue 38). 19 USE IOUNITS_MOD 20 IMPLICIT NONE 21 REAL*8 BL,C,DMDOT,DT,HICZ,OMEGAS,RL,RTOT,SMASS,TEFFL,TEMP,WIND1, 22 * * WIND2,WN,WOLD,WP,WSAT 23 INTEGER IT 24 25C 9/25/26 MHP via Claude 26C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/, 27C and COMMON/CONST3/ into a single COMMON/CONSTANTS/ 28C MHP 8/17 ADDED EXCEN, C_2 TO COMMON BLOCK FOR MATT ET AL. 2012 CENT. TERM 29C 10/12/26 MHP via Claude 30C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 31C namelist variables that the physics routines read from 27 COMMON 32C blocks which also held derived or working members. Those blocks keep 33C only their non-namelist members; five of them -- NEWMX, VMULT, 34C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 35C Derivation-only namelist variables from the same blocks went to 36C COMMON/NMLBYPASS/ instead, not here. 37C 38C Declared identically in all 121 files that hold any of these members. 39C Generated from one member list, so every copy is byte-identical; the 40C linker's 'size of symbol changed' warning is the check that they are. 41C 10/13/26 MHP via Claude: removed the unused COMMON/CWIND/ declaration 42C (no member referenced here; residue of the obsolete-option removals). 43C COMMON/CWIND/WMAX,EXMD,EXW,EXTAU,EXR,EXM,CONSTFACTOR,STRUCTFACTOR,LJDOT0 44C MHP 3/09 ADDED OPTION TO SCALE THE SATURATION RATE BY THE OVERTURN TIMESCALE 45C 10/13/26 MHP via Claude: removed the unused COMMON/DEUTER/ declaration 46C (no member referenced here; residue of the obsolete-option removals). 47 SAVE 48C COWIND RETURNS THE CHANGE IN THE ANGULAR VELOCITY OF THE SURFACE 49C CONVECTION ZONE OF A STAR WHICH IS EXPERIENCING ANGULAR MOMENTUM LOSS. 50C 51C INPUT VARIABLES: 52C 53C DECLARED 54C BL : TOTAL LUMINOSITY 55C DT : DIFFUSION TIMESTEP. 56C HICZ : TOTAL MOMENT OF INERTIA OF THE SURFACE C.Z. 57C IT : ITERATION NUMBER (IN THE DIFFUSION ROUTINES). 58C OMEGAS : SURFACE ANGULAR VELOCITY AT THE END OF THE PREVIOUS ITERATION. 59C SMASS : STELLAR MASS (SOLAR UNITS) 60C TEFFL : STELLAR EFFECTIVE TEMPERATURE. 61C WOLD : SURFACE ANGULAR VELOCITY AT THE BEGINNING OF THE TIMESTEP. 62C IN COMMON BLOCKS 63C 64C EXMD,EXW,EXR,EXM,CONSTFACTOR : USER PARAMETERS WHICH VARY THE STRENGTH 65C OF ANGULAR MOMENTUM LOSS AND ITS DEPENDENCE ON SURFACE ROTATION RATE 66C AND OTHER STELLAR PARAMETERS. 67C 68C OUTPUT VARIABLES 69C 70C WIND1,WIND2:CHANGE IN ANGULAR VELOCITY OF THE SURFACE C.Z. FOR THE TIMESTEP 71C BASED ON THE ANGULAR VELOCITY AT THE BEGINNING AND END OF THE TIMESTEP 72C RESPECTIVELY. 73C 74C WIND1 AND WIND2 ARE INCORPORATED INTO THE DIFFUSION EQUATIONS; 75C THEY ARE AVERAGED TO FIND THE TOTAL ANGULAR MOMENTUM LOSS FOR THE 76C TIMESTEP. 77C FIND TOTAL RADIUS OF STAR. 78C mhp 10/02 cgrav not used; omit 79C CGRAV = EXP(CLN*CGL) 80C MHP 3/09 IF WMAX > 1 THEN ASSUME THAT THE PARAMETER WMAX IS DEFINED BY 81C WMAX = WMAX(SUN)*TAUCZ(SUN) AND THE SATURATION THRESHOLD WSAT = WMAX/TAUCZ(STAR) 82 IF(WMAX.GT.1.0D0)THEN 83 IF(TAUCZ.GT.1.0D0)THEN 84 WSAT = WMAX/TAUCZ 85c WRITE(*,912)WSAT,TAUCZ 86c 912 FORMAT('Omega sat, Tau',1p2e12.3) 87 ELSE 88 WRITE(*,911)WMAX,TAUCZ 89 WRITE(ISHORT,911)WMAX,TAUCZ 90 911 911 FORMAT('ERROR IN WIND - TAUCZ NOT DEFINED ',1P2E12.3,'STOPPED') 91C 10/30/26 MHP via Claude: bare STOP -> STOP 1, the run-time error code (issue 38). 92 STOP 1 93 ENDIF 94 ELSE 95 WSAT = WMAX 96 ENDIF 97 RL=0.5D0*(BL+CLSUNL-C4PIL-CSIGL-4.D0*TEFFL) 98 RTOT = EXP(CLN*RL) 99C DMDOT IS THE MASS LOSS RATE IN SOLAR MASSES PER YEAR. 100 DMDOT = 2.0D-14 101C DJ/DT = DT*CONSTFACTOR*(DMDOT/1.0D-14)**EXMD*OMEGA**EXW*(M/MSUN)**EXM 102C *(R/RSUN)**EXR 103C THE CONSTANT AND EXPONENTS ARE SET IN PARMIN BASED ON THE INPUT 104C INDEX ALFA;SEE PARMIN FOR DETAILS ON THE DEPENDENCE OF EACH ON ALFA. 105 C = DT/HICZ*CONSTFACTOR*(DMDOT/1.0D-14)**EXMD 106 * * *(RTOT/CRSUN)**EXR*SMASS**EXM 107 WP = MIN(WOLD,WSAT) 108 WN = MIN(OMEGAS,WSAT) 109 WIND1 = C*WP**(EXW-1.0D0)*WOLD 110 TEMP = C*WN**(EXW-1.0D0)*OMEGAS 111 IF(IT.EQ.1) THEN 112 WIND2 = TEMP 113 ELSE 114 WIND2 = 0.5D0*(WIND2+TEMP) 115 ENDIF 116 117 RETURN 118 END

1"""Python version of cowind.f (YREC 7be4133): the change of the surface 2convection zone's angular velocity by the older (Kawaler-type) wind law.""" 3 4import math 5 6from pyrec import constants_mod 7from pyrec import iounits_mod 8from pyrec import nmlpar_mod 9from pyrec import ovrtrn_mod 10 11 12def cowind(bl: float, dt: float, hicz: float, it: int, omegas: float, smass: float, teffl: float, wold: float, 13 wind1: float, wind2: float) -> tuple[float, float]: 14 """The change of the surface convection zone's Omega over the step by 15 magnetized-wind angular-momentum loss, with a saturation threshold: 16 delta Omega = (dt/I_cz) K (Mdot/1e-14)**a (R/R_sun)**b (M/M_sun)**c 17 min(Omega, Omega_sat)**(n-1) Omega, for Omega at the start of the 18 step (wind1) and at its end (wind2). 19 20 Translates: cowind.f#COWIND (7be4133) 21 22 Arguments: bl, log L/L_sun; dt, the time step; hicz, the moment of 23 inertia of the surface convection zone; it, the iteration; omegas, 24 the surface Omega after the previous iteration; smass, M/M_sun; 25 teffl, log Teff; wold, Omega at the start of the step; wind1, passed 26 in (with any value); wind2, its value from the previous iteration. 27 28 Returns (wind1, wind2). 29 30 """ 31 32 c = constants_mod 33 n = nmlpar_mod 34 taucz = ovrtrn_mod.taucz 35 36 # WMAX > 1 is read as WMAX(sun)*TAUCZ(sun), and the saturation threshold is 37 # WMAX/TAUCZ(star) (MHP 3/09) (Fortran lines 82–96). 38 # Question (see the page): the test WMAX > 1 works as a switch. 39 if n.wmax > 1.0: 40 if taucz > 1.0: 41 wsat = n.wmax / taucz 42 else: 43 msg = (f"ERROR IN WIND - TAUCZ NOT DEFINED " 44 f"{n.wmax:12.3E}{taucz:12.3E}STOPPED") 45 print(msg) 46 iounits_mod.write(iounits_mod.ishort, msg) 47 raise SystemExit(1) # STOP 1, the run-time error code 48 else: 49 wsat = n.wmax 50 51 # The total radius of the star (cm) (Fortran lines 97–98). 52 rl = 0.5 * (bl + c.clsunl - c.c4pil - c.csigl - 4.0 * teffl) 53 rtot = math.exp(c.cln * rl) 54 55 # The mass-loss rate, in solar masses per year (Fortran line 100). 56 # Question (see the page): fixed, not the model's. 57 dmdot = 2.0e-14 58 59 # The changes at the start and end of the step; cc is the Fortran's C 60 # (c here names constants_mod) (Fortran lines 105–110). 61 cc = (dt / hicz * n.constfactor * (dmdot / 1.0e-14) ** n.exmd 62 * (rtot / c.crsun) ** n.exr * smass ** n.exm) 63 wp = min(wold, wsat) 64 wn = min(omegas, wsat) 65 wind1 = cc * wp ** (n.exw - 1.0) * wold 66 temp = cc * wn ** (n.exw - 1.0) * omegas 67 68 # WIND2 averaged with its previous value (Fortran lines 111–115). 69 # Question (see the page): the mean with the previous WIND2. 70 if it == 1: 71 wind2 = temp 72 else: 73 wind2 = 0.5 * (wind2 + temp) 74 return wind1, wind2