findsh.f

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

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

Variables on this page (2)

atime, json

Purpose

FINDSH locates, in a model, the outer edge of a central convection zone, the inner edge of a surface convection zone, and, once the core's hydrogen falls below a threshold, the beginning, middle and end of the hydrogen-burning shell. Its comment lists a planned addition ("TO BE ADDED"): finding the exact edges of the convection zones and adding points there.

How it fits

Who calls it: the main program (main.f) and RSCALE.
Inputs: HCOMP (abundances; row 1 is $X$), HL (luminosity in solar units), LC (convective flags), M (number of points).
Outputs: JCORE (outermost point of the convective core; 1 if none), JENV (innermost point of the surface convection zone; M if none), LSHELL (a hydrogen shell exists), and JXBEG, JXMID, JXEND (the shell's start, middle and end).
Modules: NMLPAR_MOD, for ATIME(1) (formerly in /CTLIM/), the central $X$ below which the star has a shell rather than a core (set in PARMIN).

The method

Shell (if $X_\mathrm{c} \le$ ATIME(1)): going outward, JXBEG counts the points with $X \le$ ATIME(1), JXMID also those with $X \le X_\mathrm{surf}/2$; the shell ends (JXEND) at the first point where $L_i - L_{i-1} < 10^{-5}L_\mathrm{surf}$ or $X_\mathrm{surf} - X_i < 10^{-5}$. The difference HL(I) - HL(I-1) would read HL(0) at the centre, which the compiler flags, but at I = 1 the first test always holds, since $X_\mathrm{c} \le$ ATIME(1) there, so HL(0) is never read (settled from the code, Nov 9, 2026: line 80 ↗; lines 85–92 ↗).

Question 1 (lines 85–91 ↗). The shell's middle and end. JXMID counts points rather than recording an index; it equals the first point with $X > X_\mathrm{surf}/2$ only if the composition increases monotonically outward (and likewise JXBEG, the first point with $X >$ ATIME(1), as the header says). Is that guaranteed?

ENDSHL is a stale name. The comments speak of ENDSHL (lines 57–59 ↗; lines 72–73 ↗), but the code uses DLMIN, the fraction, and XLEND = DLMIN*HL(M), the threshold in $L$ (line 83 ↗), so the comments could name those. 2. Core: JCORE is the last point of the run of convective points from the centre. 3. Envelope: JENV is the first point of the run of convective points from the surface. 4. For a fully convective star (JENV = 1), JCORE is set to 1, "because JCORE is used in computing the nuclear timestep, and putting it at the surface does strange things."

Definitions

subroutineFINDSH

lines 6–126 · GitHub

Finds the convective core's edge, the surface convection zone's base, and the hydrogen shell's start, middle and end (method above).

Uses modules DIMENSIONS_MOD, NMLPAR_MOD

Called by (file level) ×2, RSCALE

Compiler warnings (1)

gfortran -Wall -Wextra, compiling only; nothing was run.

line 92: Array reference out of bounds (0 < 1) in loop beginning [-Wdo-subscript]

Change log

The changes since rp_1

Changed since rp_1. Declarations only: COMMON blocks to modules, IMPLICIT NONE.

1C 2C 3C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 4C FINDSH 5C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 6 SUBROUTINE FINDSH(HCOMP,HL,LC,M,JCORE,JENV,JXBEG,JXEND,JXMID, 7 * * LSHELL) 8C 10/15/26 MHP via Claude 9C code changed to take JSON from DIMENSIONS_MOD instead of declaring 10C PARAMETER(JSON=5000) here. 11 USE DIMENSIONS_MOD 12C 10/14/26 MHP via Claude 13C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD. 14 USE NMLPAR_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 DLMIN,DXMIN,HCOMP,HL,XLEND,XMID 20 INTEGER I,JCORE,JENV,JXBEG,JXEND,JXMID,M 21 LOGICAL*4 LC,LSHELL 22C 10/12/26 MHP via Claude 23C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 24C namelist variables that the physics routines read from 27 COMMON 25C blocks which also held derived or working members. Those blocks keep 26C only their non-namelist members; five of them -- NEWMX, VMULT, 27C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 28C Derivation-only namelist variables from the same blocks went to 29C COMMON/NMLBYPASS/ instead, not here. 30C 31C Declared identically in all 121 files that hold any of these members. 32C Generated from one member list, so every copy is byte-identical; the 33C linker's 'size of symbol changed' warning is the check that they are. 34C 10/13/26 MHP via Claude 35C code changed to drop COMMON/CTLIM/, declared here but never used. 36 DIMENSION HCOMP(15,JSON),HL(JSON),LC(JSON) 37 DATA DLMIN,DXMIN/1.0D-5,1.0D-5/ 38 SAVE 39C SR FINDSH LOCATES THE OUTER EDGE OF A CENTRAL CONVECTION ZONE AND THE INNER 40C EDGE OF A SURFACE C.Z. IT ALSO LOCATES THE BEGINNING,CENTER,AND END OF 41C A HYDROGEN-BURNING SHELL IF APPLICABLE. 42C ** TO BE ADDED ** 43C LOCATING THE EXACT EDGES OF CENTRAL AND SURFACE C.Z.'S AND ADDING 44C POINTS AT THE EDGES. 45C 46C INPUT VARIABLES: 47C DECLARED: 48C HCOMP : RUN OF MASS FRACTION OF THE VARIOUS CHEMICAL SPECIES. 49C HL : RUN OF LUMINOSITY AS A FUNCTION OF MASS (SOLAR UNITS). 50C LC : FLAG SET T IF A SHELL IS CONVECTIVE. 51C M : NUMBER OF POINTS IN THE MODEL. 52C IN COMMON BLOCKS : 53C ATIME(1) : USER PARAMETER; IF X < ATIME(1) THEN THE PROGRAM CONSIDERS 54C THE STAR TO HAVE A HYDROGEN-BURNING SHELL RATHER THAN A HYDROGEN- 55C BURNING CORE. 56C LOCAL VARIABLES : 57C DLMIN : THE OUTER EDGE OF THE H-BURNING SHELL IS DEFINED AS THE POINT 58C WHERE THE CHANGE IN LUMINOSITY FROM ONE SHELL TO THE NEXT 59C IS LESS THAN ENDSHL. 60C DXMIN : THE OUTER EDGE OF THE H-BURNING SHELL IS REACHED WHEN X DIFFERS 61C FROM THE SURFACE VALUE BY LESS THAN DXMIN REGARDLESS OF THE 62C LUMINOSITY TEST. 63C 64C OUTPUT VARIABLES : 65C JCORE : OUTERMOST POINT IN THE CONVECTIVE CORE (1=NO CONVECTIVE CORE) 66C JENV : INNERMOST POINT IN THE SURFACE C.Z. (M=NO SURFACE C.Z.) 67C LSHELL : FLAG SET T IF THE PROGRAM CONSIDERS THE MODEL TO HAVE A 68C HYDROGEN-BURNING SHELL. 69C IF LSHELL=T THEN THE FOLLOWING ARE COMPUTED : 70C JXBEG : FIRST SHELL WITH X > ATIME(1). 71C JXMID : FIRST SHELL WHERE X EXCEEDS 1/2 THE SURFACE VALUE. 72C JXEND : LAST SHELL WHERE L(I)-L(I-1) > ENDSHL OR X IS WITHIN 1.0E-5 73C OF THE SURFACE VALUE. 74CCC H-SHELL VALUES 75 JXBEG = 1 76 JXMID = 1 77 JXEND = 1 78 LSHELL = .FALSE. 79C IF CENTRAL X BELOW THRESHOLD THEN CALCULATE H SHELL VALUES 80 IF(HCOMP(1,1).LE.ATIME(1)) THEN 81 LSHELL = .TRUE. 82 XMID = 0.50D0*HCOMP(1,M) 83 XLEND = DLMIN*HL(M) 84C FIND BEGINNING(JXBEG), MIDDLE(JXMID) AND END(JXEND) OF H SHELL 85 DO 10 I = 1,M 86c IF(HCOMP(1,I).LE.1.0D-10) THEN ! Changed after discussion with Marc 87 IF(HCOMP(1,I).LE.ATIME(1)) THEN ! to force consistency with above LLP 9/24/08 88 JXBEG = JXBEG+1 89 JXMID = JXMID+1 90 ELSE IF(HCOMP(1,I).LE.XMID) THEN 91 JXMID = JXMID+1 92 ELSE IF(HL(I) - HL(I-1).LT.XLEND) THEN 93c write(*,*)'luminosity criteria' 94 GOTO 20 95 ELSE IF(HCOMP(1,M) - HCOMP(1,I).LT.DXMIN) THEN 96c write(*,*)'composition criteria' 97 GOTO 20 98 ENDIF 99 10 10 CONTINUE 100 I = M 101 20 20 JXEND = I 102 ENDIF 103CCC FIND BOUNDARY OF CENTRAL CONVECTION ZONE. 104 DO 30 I = 1,M 105 IF(.NOT.LC(I)) GOTO 40 106 30 30 CONTINUE 107 40 40 IF(I.GT.1) THEN 108 JCORE = I-1 109 ELSE 110 JCORE = 1 111 ENDIF 112CCC FIND BOUNDARY OF SURFACE C.Z. 113 DO 50 I = M,1,-1 114 IF(.NOT.LC(I)) GOTO 60 115 50 50 CONTINUE 116 60 60 IF(I.LT.M) THEN 117 JENV = I+1 118 ELSE 119 JENV = M 120 ENDIF 121C FOR A FULLY CONVECTIVE STAR (JENV=1),TURN THE CONVECTIVE CORE OFF(JCORE=1). 122C THIS IS DONE BECAUSE JCORE IS USED IN COMPUTING THE NUCLEAR TIMESTEP, 123C AND PUTTING IT AT THE SURFACE DOES STRANGE THINGS. 124 IF(JENV.EQ.1) JCORE = 1 125 RETURN 126 END

1"""Python version of findsh.f (YREC 7be4133): the edges of the 2convective core, the surface convection zone and the hydrogen shell.""" 3 4from pyrec import nmlpar_mod 5import numpy as np 6 7# DATA: the shell ends where L grows by less than dlmin * L_surf, or X is 8# within dxmin of its surface value. 9dlmin = 1.0e-5 10dxmin = 1.0e-5 11 12 13def findsh(hcomp: np.ndarray, hl: np.ndarray, lc: np.ndarray, m: int, jcore: int, jenv: int, jxbeg: int, jxend: int, 14 jxmid: int, lshell: bool) -> tuple[int, int, int, int, int, bool]: 15 """Find the outermost point of the convective core, the innermost 16 point of the surface convection zone and, once the central X is at 17 most nmlpar_mod.atime[1], the start, middle and end of the hydrogen 18 shell. 19 20 Translates: findsh.f#FINDSH (7be4133) 21 22 Arguments: hcomp (abundances, row 1 is X), hl (luminosities), lc 23 (convective flags), m (number of points); jcore, jenv, jxbeg, jxend, 24 jxmid, lshell are only set. 25 26 Returns (jcore, jenv, jxbeg, jxend, jxmid, lshell). 27 28 """ 29 30 # Fortran lines 74–78. 31 jxbeg = 1 32 jxmid = 1 33 jxend = 1 34 lshell = False 35 36 # The hydrogen shell, if the central X is below the threshold (Fortran 37 # lines 79–102). 38 if hcomp[1, 1] <= nmlpar_mod.atime[1]: 39 lshell = True 40 xmid = 0.50 * hcomp[1, m] 41 xlend = dlmin * hl[m] 42 # Question (see the page): jxmid counts points rather than 43 # recording an index. 44 i = m # if no end is found 45 for ii in range(1, m + 1): # label 10 46 if hcomp[1, ii] <= nmlpar_mod.atime[1]: 47 jxbeg = jxbeg + 1 48 jxmid = jxmid + 1 49 elif hcomp[1, ii] <= xmid: 50 jxmid = jxmid + 1 51 elif hl[ii] - hl[ii - 1] < xlend: 52 i = ii 53 break # GOTO 20 54 elif hcomp[1, m] - hcomp[1, ii] < dxmin: 55 i = ii 56 break # GOTO 20 57 jxend = i # label 20 58 59 # The core: the last point of the convective run from the centre 60 # (Fortran lines 103–111). 61 i = m + 1 # Fortran's index after a loop that runs out 62 for ii in range(1, m + 1): # label 30 63 if not lc[ii]: 64 i = ii 65 break # GOTO 40 66 if i > 1: # label 40 67 jcore = i - 1 68 else: 69 jcore = 1 70 71 # The envelope: the first point of the convective run from the surface 72 # (Fortran lines 112–120). 73 i = 0 74 for ii in range(m, 0, -1): # label 50 75 if not lc[ii]: 76 i = ii 77 break # GOTO 60 78 if i < m: # label 60 79 jenv = i + 1 80 else: 81 jenv = m 82 83 # Fully convective: JCORE is used in the nuclear time step, and 84 # putting it at the surface does strange things (Fortran lines 85 # 121–124). 86 if jenv == 1: 87 jcore = 1 88 89 return jcore, jenv, jxbeg, jxend, jxmid, lshell