amcalc.f

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

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

Variables on this page (18)

c4pil, cln, clsunl, crsun, csigl, exl, exm, expr, exr, extau, fracstep, lrossby, pmmsolp, pmmsoltau, pphot, pphot0, taucz, taucz0

Purpose

This file computes the structure-dependent factor of the magnetic-braking torque, "according to the Matt et al. (2012) formulation" (comment by G. Somers, 6/16). The factor scales the angular-momentum loss of the wind with the star's mass, radius, luminosity, photospheric pressure and, optionally, its convective overturn time.

How it fits

Inputs: the arguments SMASS (the star's mass in solar masses: MDOT updates it as SMASS + DELM/CMSUN, and COWIND's header says "stellar mass (solar units)"), BL ($\log L/L_\odot$) and TEFFL ($\log T_\mathrm{eff}$); the exponents EXM, EXR, EXL, EXPR, EXTAU, the solar normalizations PMMSOLP and PMMSOLTAU and the switch LROSSBY, all namelist parameters in NMLPAR_MOD (set in PARMIN); the photospheric pressure PPHOT, PPHOT0, the overturn time TAUCZ, TAUCZ0 and FRACSTEP in OVRTRN_MOD.
Who calls it: MWIND and MCOWIND.
Output: the argument STRUCTFACTOR, which MWIND and MCOWIND multiply with CONSTFACTOR and the time step to get the spin-down.
Modules: CONSTANTS_MOD (CLSUNL, C4PIL, CSIGL, CLN, CRSUN), NMLPAR_MOD and OVRTRN_MOD.

The method

The radius from the Stefan–Boltzmann law, in cgs units: $$\log R = \tfrac12\left(\log L - \log 4\pi - \log\sigma - 4\log T_\mathrm{eff}\right),$$ with $\log L$ = BL + CLSUNL; then $R/R_\odot$ (RRAD).
$L/L_\odot = 10^{\mathtt{BL}}$.
The photospheric pressure, interpolated within the step between its values at the start (PPHOT0) and end (PPHOT) of the step, both $\log P$, and normalized to the solar value: $$\frac{P_\mathrm{phot}}{P_{\mathrm{phot},\odot}} = \frac{10^{\,p_0 + f\,(p - p_0)}}{10^{\,p_\odot}},$$ with $f$ = FRACSTEP.
If LROSSBY, the overturn time likewise, $\tau/\tau_\odot = [\tau_0 + f(\tau - \tau_0)]/\tau_\odot$; otherwise 1.
The factor: $$\mathtt{STRUCTFACTOR} = M^{a_M}\,\left(\frac{R}{R_\odot}\right)^{a_R}\left(\frac{L}{L_\odot}\right)^{a_L}\left(\frac{P_\mathrm{phot}}{P_{\mathrm{phot},\odot}}\right)^{a_P}\left(\frac{\tau}{\tau_\odot}\right)^{a_\tau}.$$

Design

The radius comes from $L$ and $T_\mathrm{eff}$, not from the model's outer radius. The mass is used as given, with no solar normalization in this routine.

From YREC's documentation

Physics Options: PMMSOLP (default 4.9304) is $\log P$ at $T_\mathrm{eff}$ in the solar model, for the magnetic field's calibration, and PMMSOLTAU (9.4805D5 s) the Sun's convective overturn time; LROSSBY scales the loss with the Rossby number.

Definitions

subroutineAMCALC

lines 5–89 · GitHub

Sets STRUCTFACTOR, the product of $M$, $R/R_\odot$, $L/L_\odot$, $P_\mathrm{phot}/P_{\mathrm{phot},\odot}$ and optionally $\tau/\tau_\odot$, each to its exponent from NMLPAR_MOD (Matt et al. 2012), and returns it in STRUCTFACTOR.

Uses modules CONSTANTS_MOD, NMLPAR_MOD, OVRTRN_MOD

Called by MCOWIND, MWIND

Change log

The changes since rp_1

Changed since rp_1. STRUCTFACTOR is now an output argument (/CWIND/ retired); COMMON blocks to modules, IMPLICIT NONE; the unused constants and AWIND are gone, which settles the former question on them.

1C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 2C AMCALC 3C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 4 5 SUBROUTINE AMCALC(SMASS,BL,TEFFL,STRUCTFACTOR) 6C 10/14/26 MHP via Claude 7C code changed to convert COMMON/OVRTRN/ into OVRTRN_MOD. 8 USE OVRTRN_MOD 9C 10/14/26 MHP via Claude 10C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD. 11 USE NMLPAR_MOD 12C 10/14/26 MHP via Claude 13C code changed to convert COMMON/CONSTANTS/ into CONSTANTS_MOD. 14 USE CONSTANTS_MOD 15C 10/20/26 MHP via Claude 16C code changed to replace implicit typing with IMPLICIT NONE and 17C explicit declarations. 18 IMPLICIT NONE 19 REAL*8 BL,RL,RLUM,RMASS,RPHOT,RRAD,RTAU,RTOT,SMASS,STRUCTFACTOR, 20 * * TEFFL 21 22 23C 9/25/26 MHP via Claude 24C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/, 25C and COMMON/CONST3/ into a single COMMON/CONSTANTS/ 26C 9/24/26 MHP via Claude 27C code changed to remove dead CMIXL from COMMON/CONST3/; it is now 28C passed as an explicit argument where actually used 29C MHP 8/17 ADDED EXCEN, C_2 TO COMMON BLOCK FOR MATT ET AL. 2012 CENT. TERM 30C 10/12/26 MHP via Claude 31C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 32C namelist variables that the physics routines read from 27 COMMON 33C blocks which also held derived or working members. Those blocks keep 34C only their non-namelist members; five of them -- NEWMX, VMULT, 35C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 36C Derivation-only namelist variables from the same blocks went to 37C COMMON/NMLBYPASS/ instead, not here. 38C 39C Declared identically in all 121 files that hold any of these members. 40C Generated from one member list, so every copy is byte-identical; the 41C linker's 'size of symbol changed' warning is the check that they are. 42C 10/13/26 MHP via Claude 43C code changed to retire COMMON/CWIND/. Its namelist members moved to the 44C consolidated block earlier, leaving two. C_2 is a constant, set once in 45C parmin.f, so it is now the last member of COMMON/CONSTANTS/ -- appended 46C there rather than inserted, and added to every declaration so the block 47C keeps one size; parmin.f DATA-initializes CLSUN and CRSUN, so a short 48C declaration anywhere would be the OP95INTERP hazard again. It may 49C eventually become a namelist variable. STRUCTFACTOR, the loss-law 50C prefactor, is computed by AMCALC and read by its own two callers, 51C mwind.f and mcowind.f, so it is now an AMCALC output argument. 52C COMMON/CWIND/WMAX,EXMD,EXW,EXTAU,EXR,EXM,CONSTFACTOR,STRUCTFACTOR,LJDOT0 53C KC 2025-05-30 reordered common block elements 54C COMMON/PMMWIND/AWIND,LMWIND,LROSSBY,LBSCALE,PMMA,PMMB,PMMC,PMMD,PMMM, 55C * PMMJD,PMMMD,PMMSOLP,PMMSOLW,PMMSOLTAU 56C 10/12/26 MHP via Claude 57C code changed to remove the trailing CHARACTER member AWIND from 58C COMMON/PMMWIND/. Character data does not belong in a COMMON block, and 59C this file never referenced AWIND. It is now a local of parmin.f, the 60C only routine that reads it. Being last, its removal moves no other 61C member, and the block shrinks identically in all five files. 62 SAVE 63 64C G Somers, 6/16 65C THIS SHORT SUBROUTINE WILL USE THE DELIVERED STRUCTURAL VARIABLES TO DETERMINE 66C THE STRENGTH OF THE TORQUE, ACCORDING TO THE MATT ET AL. (2012) FORMULATION. 67C 68C COLLECT THE RELEVANT STRUCTURE VARIABLES, AND 69C MASS 70 RMASS = SMASS 71C RADIUS 72 RL = 0.5D0*(BL+CLSUNL-C4PIL-CSIGL-4.D0*TEFFL) 73 RTOT = DEXP(CLN*RL) 74 RRAD = RTOT/CRSUN 75C LUMINOSITY 76 RLUM = 10.**BL 77C PHOTOSPHERIC PRESSURE 78 RPHOT = 10.**(PPHOT0+FRACSTEP*(PPHOT-PPHOT0))/(10.**PMMSOLP) 79C CONVECTIVE OVERTURN TIMESCALE 80 IF(LROSSBY)THEN 81 RTAU = (TAUCZ0+FRACSTEP*(TAUCZ-TAUCZ0))/PMMSOLTAU 82 ELSE 83 RTAU = 1. 84 ENDIF 85C COMBINE THEM ALL 86 STRUCTFACTOR = RMASS**EXM * RRAD**EXR * RLUM**EXL * RPHOT**EXPR 87 * * * RTAU**EXTAU 88 RETURN 89 END

1"""Python version of amcalc.f (YREC 7be4133): the structure factor of the 2magnetic-braking torque (Matt et al. 2012).""" 3 4import math 5 6from pyrec import constants_mod 7from pyrec import nmlpar_mod 8from pyrec import ovrtrn_mod 9 10 11def amcalc(smass: float, bl: float, teffl: float, structfactor: float) -> float: 12 """Compute the structure factor of the wind torque: the mass, R/Rsun, 13 L/Lsun, the photospheric pressure over the solar one and, if LROSSBY, 14 the convective overturn time over the solar one, each raised to its 15 exponent from nmlpar_mod. 16 17 Translates: amcalc.f#AMCALC (7be4133) 18 19 Arguments: smass (mass, solar units), bl (log L/Lsun), teffl 20 (log Teff); structfactor is output only (the value passed is unused). 21 22 Returns structfactor. 23 24 """ 25 26 # Mass (Fortran line 70). 27 rmass = smass 28 29 # Radius, from L and Teff (Stefan-Boltzmann), in solar units (Fortran 30 # lines 72–74). 31 rl = 0.5 * (bl + constants_mod.clsunl - constants_mod.c4pil 32 - constants_mod.csigl - 4.0 * teffl) 33 rtot = math.exp(constants_mod.cln * rl) 34 rrad = rtot / constants_mod.crsun 35 36 # Luminosity (Fortran line 76). 37 rlum = 10.0 ** bl 38 39 # Photospheric pressure, interpolated within the step (Fortran line 78). 40 rphot = (10.0 ** (ovrtrn_mod.pphot0 + ovrtrn_mod.fracstep 41 * (ovrtrn_mod.pphot - ovrtrn_mod.pphot0)) 42 / 10.0 ** nmlpar_mod.pmmsolp) 43 44 # Convective overturn timescale (Fortran lines 80–84). 45 if nmlpar_mod.lrossby: 46 rtau = (ovrtrn_mod.taucz0 + ovrtrn_mod.fracstep 47 * (ovrtrn_mod.taucz - ovrtrn_mod.taucz0)) / nmlpar_mod.pmmsoltau 48 else: 49 rtau = 1.0 50 51 # Combine them all (Fortran lines 86–87). 52 structfactor = (rmass ** nmlpar_mod.exm * rrad ** nmlpar_mod.exr 53 * rlum ** nmlpar_mod.exl * rphot ** nmlpar_mod.expr 54 * rtau ** nmlpar_mod.extau) 55 return structfactor