tauint.f

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

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

Variables on this page (9)

cgl, cln, cmsun, envr, evels, json, lrot, numenv, taucz

Purpose

TAUINT computes the convective overturn time scale $\tau_\mathrm{cz}$ (TAUCZ) at the base of the surface convection zone, on the whole star's profile, interior and envelope joined (JVS 02/12, "snagged from MIDMOD"). It is the older method; tauintnew.f is the newer one (G. Somers 3/17).

How it fits

Inputs: the profile built by GETTAU: the interior shells only, since GETTAU stitches the envelope onto the profile only when LCALCENV and LNEWTCZ are both set, and calls TAUINT only when LNEWTCZ is false. HS1 (mass), LCZ, HR, HP, HD, HG, SVEL, and DEL1, DEL2, which GETTAU takes from SDEL(1,·) and SDEL(3,·), the radiative and adiabatic gradients; and, through ENVSTRUCT_MOD, the envelope from the last ENVINT call.
Who calls it: only GETTAU, when LNEWTCZ is false and TAUCZ is still 0 after its ENVINT call (which sets it through TAUCAL when the top of the interior model is radiative).
Output: TAUCZ in OVRTRN_MOD.
Modules: OVRTRN_MOD, ENVSTRUCT_MOD, NMLPAR_MOD (LROT) and CONSTANTS_MOD.

The method

If the top shell or the one below it is convective, the zone's base IMAX is the lowest shell of the convective run from the surface. The overturn time is computed only if LCZSUR is set, which the test (HSTOP-HSBOT)/CMSUN .GT. 0 makes true for any such zone.

Partly convective star: as in taucal.f, step 2 and 4: the base $r_b$ where $\nabla_1 - \nabla_2$ changes sign, $H_p = P/(\rho g)$ there, and $$\tau_\mathrm{cz} = \frac{H_p}{v(r_b + H_p)} ,$$ with $v$ interpolated linearly in $\log r$. If $r_b + H_p$ lies beyond the profile, the envelope in ENVSTRUCT_MOD is searched for the first convective point above it ("stitch on the envelope for more room"; the profile holds the interior only, so TAUINT searches the envelope itself), and $v$ is interpolated between the last shell and that point.

Question 2 (lines 126–132 ↗). The envelope interpolation. Beyond the profile, $v$ is interpolated between the last shell M and envelope point K, with FX = (RTESTL-HR(M))/(ENVR(K)-HR(M)), skipping envelope points 2..K-1, so a velocity profile with structure between the fitting point and point K is replaced by a straight line. Should the interpolation use the two envelope points that bracket $r_b + H_p$?
Fully convective star: the point where $H_p = r$, by linear interpolation of $H_p - r$; then $\tau_\mathrm{cz} = H_p/v$ there. If already $H_p \le r$ at the first shell, it takes $H_p \approx \sqrt{H_{p,1}\,r_1}$ and $v_1$ (the comment: "assume v constant inside and $H_p = k/r$" near the centre).
No surface zone: $\tau_\mathrm{cz} = 0$.

Definitions

subroutineTAUINT

lines 3–196 · GitHub

The older overturn time: $H_p/v$ one pressure scale height above the base of the surface convection zone, or at $H_p = r$ for a fully convective star.

Uses modules CONSTANTS_MOD, DIMENSIONS_MOD, ENVSTRUCT_MOD, NMLPAR_MOD, OVRTRN_MOD

Called by GETTAU

Change log

The changes since rp_1

Changed since rp_1. Declarations only: COMMON blocks to modules (/OVRTRN/ to OVRTRN_MOD, /ENVSTRUCT/ to ENVSTRUCT_MOD, /ROT/ to NMLPAR_MOD, the constants to CONSTANTS_MOD; the unused /COMP/, /DEUTER/ and /JTEST/ dropped, which settles the former question 3), the array size JSON from DIMENSIONS_MOD, IMPLICIT NONE.

1C$$$$$$ 2C SUBROUTINE TAUINT(HCOMP,HS2,HS1,LCZ,HR,HP,HD,HG,M,SVEL,DEL1,DEL2) ! KC 2025-05-31 3 SUBROUTINE TAUINT(HS1,LCZ,HR,HP,HD,HG,M,SVEL,DEL1,DEL2) 4C 10/15/26 MHP via Claude 5C code changed to take JSON from DIMENSIONS_MOD instead of declaring 6C PARAMETER(JSON=5000) here. 7 USE DIMENSIONS_MOD 8C 10/14/26 MHP via Claude 9C code changed to convert COMMON/OVRTRN/ into OVRTRN_MOD. 10 USE OVRTRN_MOD 11C 10/14/26 MHP via Claude 12C code changed to convert COMMON/ENVSTRUCT/ into ENVSTRUCT_MOD. 13 USE ENVSTRUCT_MOD 14C 10/14/26 MHP via Claude 15C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD. 16 USE NMLPAR_MOD 17C 10/14/26 MHP via Claude 18C code changed to convert COMMON/CONSTANTS/ into CONSTANTS_MOD. 19 USE CONSTANTS_MOD 20C 10/20/26 MHP via Claude 21C code changed to replace implicit typing with IMPLICIT NONE and 22C explicit declarations. 23 IMPLICIT NONE 24 REAL*8 CVEL,DD1,DD2,DEL1,DEL2,ENVRC,ENVRCL,FX,HD,HG,HP,HR,HS1, 25 * * HSBOT,HSTOP,PS1,PS2,PSCA,PSCA1,PSCA2,RTEST1,RTEST2,RTESTL, 26 * * SVEL 27 INTEGER I,IMAX,K,M 28 LOGICAL*4 LALLCZ,LCZ,LCZSUR 29C 10/13/26 MHP via Claude: removed the unused COMMON/COMP/ declaration 30C (no member referenced here; residue of the obsolete-option removals). 31C 10/13/26 MHP via Claude: removed the unused COMMON/DEUTER/ declaration 32C (no member referenced here; residue of the obsolete-option removals). 33C 10/12/26 MHP via Claude 34C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 35C namelist variables that the physics routines read from 27 COMMON 36C blocks which also held derived or working members. Those blocks keep 37C only their non-namelist members; five of them -- NEWMX, VMULT, 38C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 39C Derivation-only namelist variables from the same blocks went to 40C COMMON/NMLBYPASS/ instead, not here. 41C 42C Declared identically in all 121 files that hold any of these members. 43C Generated from one member list, so every copy is byte-identical; the 44C linker's 'size of symbol changed' warning is the check that they are. 45C 9/25/26 MHP via Claude 46C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/, 47C and COMMON/CONST3/ into a single COMMON/CONSTANTS/ 48C 10/13/26 MHP via Claude 49C code changed to remove COMMON/JTEST/: IMAX1, IMAX2 and LJVS are 50C referenced nowhere in executable code. The only code that used LJVS is 51C a commented-out block in wrtout.f, removed with it. 52C Needs access to this common block: grants knowledge of envellope 53C KC 2025-05-30 reordered common block elements 54C JvS 08/25 Updated with new elements 55C G Somers 3/17, ADDING NEW TAUCZ COMMON BLOCK 56C DIMENSION HCOMP(15,JSON),HS2(JSON),LCZ(JSON), 57 DIMENSION LCZ(JSON), 58 * * HR(JSON), HP(JSON), HD(JSON), HG(JSON), HS1(JSON), ! XSUM(15) 59 * * DEL1(JSON), DEL2(JSON), SVEL(JSON) 60 SAVE 61 62 63C JVS 02/12 CALCULATE THE LOCAL CONVECTIVE OVERTURN TIMESCALE AT THE BASE 64C OF THE CZ. IN OLDER VERSIONS THIS WAS ONLY DONE FOR ROTATING MODELS; 65C THIS MAKES IT SO TAUCZ IS CALCULATED FOR ALL MODELS. 66C THIS CODE SNAGGED FROM MIDMOD 67C DETERMINE EXTENT OF SURFACE CONVECTION ZONE. 68 LALLCZ = .FALSE. 69C JVS Allows the last point to be stable. 70 IF(LCZ(M) .OR. LCZ(M-1))THEN 71C SURFACE C.Z. EXISTS. FIND LOWEST SHELL (IMAX), WHICH IS ALSO THE 72C UPPERMOST ZONE CONSIDERED FOR STABILITY AGAINST ROTATIONALLY INDUCED MIXING. 73 DO 71 I = M-1,1,-1 74 IF(.NOT.LCZ(I)) GOTO 81 75 71 71 CONTINUE 76 LALLCZ = .TRUE. 77 I = 0 78 81 81 IMAX = I + 1 79C HSTOP IS THE MASS AT THE TOP OF THE C.Z. 80C HSBOT IS THE MASS AT THE BOTTOM OF THE C.Z. 81 HSTOP = HS1(M) 82 IF(IMAX.GT.1)THEN 83 HSBOT = 0.5D0*(HS1(IMAX)+HS1(IMAX-1)) 84 ELSE 85 HSBOT = 0.0D0 86 ENDIF 87C LCZSUR=T IF A SURFACE C.Z.DEEP ENOUGH FOR ANGULAR MOMENTUM LOSS EXISTS 88 IF((HSTOP-HSBOT)/CMSUN.GT.0.0D0)THEN 89 LCZSUR = .TRUE. 90 ELSE 91 LCZSUR= .FALSE. 92 ENDIF 93 ELSE 94C NO SURFACE C.Z. 95 IMAX = M 96 LCZSUR = .FALSE. 97 ENDIF 98 99 IF(LCZSUR)THEN 100 IF(.NOT.LALLCZ)THEN 101 IF(.NOT.LROT)THEN 102 HG(IMAX)=HS1(IMAX)*EXP(CLN*(CGL-2.0D0*HR(IMAX))) 103 HG(IMAX-1)=HS1(IMAX-1)*EXP(CLN*(CGL-2.0D0*HR(IMAX-1))) 104 ENDIF 105C PINPOINT RCZ 106 DD2 = DEL1(IMAX-1)-DEL2(IMAX-1) 107 DD1 = DEL1(IMAX)-DEL2(IMAX) 108 FX = DD2/(DD2-DD1) 109C INFER HP 110 ENVRCL = HR(IMAX-1)+FX*(HR(IMAX)-HR(IMAX-1)) 111 ENVRC = EXP(CLN*ENVRCL) 112 PS2 = EXP(CLN*(HP(IMAX)-HD(IMAX)))/HG(IMAX) 113 PS1 = EXP(CLN*(HP(IMAX-1)-HD(IMAX-1)))/HG(IMAX-1) 114 PSCA = PS1 + FX*(PS2-PS1) 115 RTESTL = DLOG10(ENVRC+PSCA) 116C FIND V 117 DO K = IMAX+1,M 118 IF(HR(K).GT.RTESTL)THEN 119 FX = (RTESTL-HR(K-1))/(HR(K)-HR(K-1)) 120 CVEL = SVEL(K-1)+FX*(SVEL(K)-SVEL(K-1)) 121 GOTO 85 122 ENDIF 123 END DO 124 ! One pressure scale height overshoots edge of interior 125 ! calculation. Stitch on the envelope for more room 126 DO K = 2,NUMENV 127 IF(ENVR(K).GT.RTESTL .AND. EVELS(K) .GT. 0.0)THEN 128 FX = (RTESTL-HR(M))/(ENVR(K)-HR(M)) 129 CVEL = SVEL(M)+FX*(EVELS(K)-SVEL(M)) 130 GOTO 85 131 ENDIF 132 END DO 133 CVEL = SVEL(M) 134 85 85 CONTINUE 135C DEFINE TAUCZ 136 TAUCZ = PSCA/CVEL 137 138 ELSE 139C INFER HP 140 IF(.NOT.LROT)THEN 141 HG(1)=HS1(1)*EXP(CLN*(CGL-2.0D0*HR(1))) 142 ENDIF 143 PSCA2 = EXP(CLN*(HP(1)-HD(1)))/HG(1) 144 RTEST2 = EXP(CLN*HR(1)) 145 IF(PSCA2.LE.RTEST2)THEN 146C HP < R AT THE FIRST POINT. ASSUME V CONSTANT INSIDE AND HP = K/R FOR 147C SLOWLY VARYING DENSITY AND PRESSURE NEAR THE CENTER. 148 CVEL = SVEL(1) 149 PSCA = (PSCA2*RTEST2)**0.5D0 150 TAUCZ = PSCA/CVEL 151 ELSE 152 DO K = 2,M 153 PSCA1 = PSCA2 154 RTEST1 = RTEST2 155 IF(.NOT.LROT)THEN 156 HG(K)=HS1(K)*EXP(CLN*(CGL-2.0D0*HR(K))) 157 ENDIF 158 PSCA2 = EXP(CLN*(HP(K)-HD(K)))/HG(K) 159 RTEST2 = EXP(CLN*HR(K)) 160C FIND LOCATION WHERE HP = R 161 IF(PSCA2.LE.RTEST2)THEN 162 FX = (RTEST1-PSCA1)/((PSCA2-RTEST2)-(PSCA1-RTEST1)) 163C FIND V 164 CVEL = SVEL(K-1)+FX*(SVEL(K)-SVEL(K-1)) 165 PSCA = PSCA1+FX*(PSCA2-PSCA1) 166C DEFINE TAUCZ 167 TAUCZ = PSCA/CVEL 168 GOTO 95 169 ENDIF 170 END DO 171 K = M 172 CVEL = SVEL(M) 173 PSCA = PSCA2 174 TAUCZ = PSCA/CVEL 175 95 95 CONTINUE 176 ENDIF 177 ENDIF 178 ELSE 179 TAUCZ = 0.0D0 180 ENDIF 181 182C-------------------------------------------------------------- 183C OPEN(UNIT=100,FILE='diagnostic.out',STATUS='OLD') 184C DO 1505 I=1,M 185C IF (LCZ(I)) HOLD = 1.0 186C IF (.NOT.LCZ(I)) HOLD = 0.0 187C WRITE(100,1504), HS2(I), HS1(I), HR(I),HOLD, 188C * DEL1(I), DEL2(I), SVEL(I), 189C1505 CONTINUE 190C1504 FORMAT(1X,7E16.8) 191C CLOSE(100) 192C---------------------------------------------------------------- 193 194C END JVS 195 RETURN 196 END

1"""Python version of tauint.f (YREC 7be4133): the convective overturn time 2at the base of the surface convection zone, on the interior profile (the 3older method; TAUINTNEW is the newer).""" 4 5import math 6import numpy as np 7 8from pyrec import constants_mod 9from pyrec import envstruct_mod 10from pyrec import nmlpar_mod 11from pyrec import ovrtrn_mod 12 13 14def tauint(hs1: np.ndarray, lcz: np.ndarray, hr: np.ndarray, hp: np.ndarray, hd: np.ndarray, hg: np.ndarray, m: int, 15 svel: np.ndarray, del1: np.ndarray, del2: np.ndarray) -> None: 16 """Set ovrtrn_mod.taucz: H_P / v one pressure scale height above the 17 base of the surface convection zone (searching the envelope of the 18 last ENVINT call when the profile ends first); for a fully convective 19 star, H_P / v where H_P = r; 0 with no surface zone. 20 21 Translates: tauint.f#TAUINT (7be4133) 22 23 Arguments: hs1 (masses), lcz (convective flags), hr, hp, hd (log r, 24 log P, log rho), hg (gravity; changed in place where computed, when 25 there is no rotation), m, svel (convective velocities), del1, del2 26 (the radiative and adiabatic gradients). None is returned. 27 28 """ 29 30 c = constants_mod 31 env = envstruct_mod 32 lrot = nmlpar_mod.lrot 33 lallcz = False 34 35 # IMAX, the zone's mass range and LCZSUR (Fortran lines 68–97). 36 if lcz[m] or lcz[m - 1]: 37 # IMAX: the lowest shell of the convective run from the surface. 38 for i in range(m - 1, 0, -1): 39 if not lcz[i]: 40 break # to label 81 41 else: 42 lallcz = True 43 i = 0 44 imax = i + 1 # label 81 45 hstop = hs1[m] 46 if imax > 1: 47 hsbot = 0.5 * (hs1[imax] + hs1[imax - 1]) 48 else: 49 hsbot = 0.0 50 # Question (see the page): true for any zone. 51 if (hstop - hsbot) / c.cmsun > 0.0: 52 lczsur = True 53 else: 54 lczsur = False 55 else: 56 imax = m 57 lczsur = False 58 59 # No surface zone (Fortran lines 99, 178–179). 60 if not lczsur: 61 ovrtrn_mod.taucz = 0.0 62 return 63 64 if not lallcz: 65 # A partly convective star: the base r_b and H_P there (Fortran 66 # lines 101–115). 67 if not lrot: 68 hg[imax] = hs1[imax] * math.exp(c.cln * (c.cgl - 2.0 * hr[imax])) 69 hg[imax - 1] = hs1[imax - 1] * math.exp( 70 c.cln * (c.cgl - 2.0 * hr[imax - 1])) 71 dd2 = del1[imax - 1] - del2[imax - 1] 72 dd1 = del1[imax] - del2[imax] 73 fx = dd2 / (dd2 - dd1) 74 envrcl = hr[imax - 1] + fx * (hr[imax] - hr[imax - 1]) 75 envrc = math.exp(c.cln * envrcl) 76 ps2 = math.exp(c.cln * (hp[imax] - hd[imax])) / hg[imax] 77 ps1 = math.exp(c.cln * (hp[imax - 1] - hd[imax - 1])) / hg[imax - 1] 78 psca = ps1 + fx * (ps2 - ps1) 79 rtestl = math.log10(envrc + psca) 80 81 # The velocity one scale height above the base (Fortran lines 82 # 117–136). 83 found = False 84 for k in range(imax + 1, m + 1): 85 if hr[k] > rtestl: 86 fx = (rtestl - hr[k - 1]) / (hr[k] - hr[k - 1]) 87 cvel = svel[k - 1] + fx * (svel[k] - svel[k - 1]) 88 found = True 89 break # to label 85 90 if not found: 91 # One pressure scale height overshoots the interior: stitch on 92 # the envelope for more room. 93 # Question (see the page): between shell M and envelope point 94 # K, skipping the points before K. 95 for k in range(2, env.numenv + 1): 96 if env.envr[k] > rtestl and env.evels[k] > 0.0: 97 fx = (rtestl - hr[m]) / (env.envr[k] - hr[m]) 98 cvel = svel[m] + fx * (env.evels[k] - svel[m]) 99 break # to label 85 100 else: 101 cvel = svel[m] 102 ovrtrn_mod.taucz = psca / cvel # label 85 103 return 104 105 # A fully convective star: where H_P = r (Fortran lines 140–175). 106 if not lrot: 107 hg[1] = hs1[1] * math.exp(c.cln * (c.cgl - 2.0 * hr[1])) 108 psca2 = math.exp(c.cln * (hp[1] - hd[1])) / hg[1] 109 rtest2 = math.exp(c.cln * hr[1]) 110 if psca2 <= rtest2: 111 cvel = svel[1] 112 psca = (psca2 * rtest2) ** 0.5 113 ovrtrn_mod.taucz = psca / cvel 114 return 115 for k in range(2, m + 1): 116 psca1 = psca2 117 rtest1 = rtest2 118 if not lrot: 119 hg[k] = hs1[k] * math.exp(c.cln * (c.cgl - 2.0 * hr[k])) 120 psca2 = math.exp(c.cln * (hp[k] - hd[k])) / hg[k] 121 rtest2 = math.exp(c.cln * hr[k]) 122 if psca2 <= rtest2: 123 fx = (rtest1 - psca1) / ((psca2 - rtest2) - (psca1 - rtest1)) 124 cvel = svel[k - 1] + fx * (svel[k] - svel[k - 1]) 125 psca = psca1 + fx * (psca2 - psca1) 126 ovrtrn_mod.taucz = psca / cvel 127 return # to label 95 128 cvel = svel[m] 129 psca = psca2 130 ovrtrn_mod.taucz = psca / cvel