microdiff_setup.f

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

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

Variables on this page (9)

cln, con_mass, con_rad, con_temp, con_time, ishort, json, xmin, ymin

Purpose

The first step of MICRODIFF: it decides which part of the star takes part in settling, and converts the model's variables to the units of Bahcall and Loeb (1989), as its comments call them.

How it fits

Inputs: the arguments of MICRODIFF (see there); XMIN and YMIN from NMLPAR_MOD (namelist variables read by PARMIN); the solar constants in CONSTANTS_MOD (not used: the routine has its own constants).
Who calls it: only MICRODIFF.
Outputs: IBEGIN, IEND, LALLCZ; the unlogged radius HRU and temperature HTU in the new units; the unlogged density ODEN and temperature OTEM in cgs units; and, in place, HS1, HQPR, DT and STOT, converted to the new units. MICRODIFF_ETM converts them back at the end.
Modules: it writes the conversion factors to CONFAC_MOD, which MICRODIFF_ETM reads (as do SETUP_GRSETT and EQUAL_TO_MODEL in the older settling path). It reads NMLPAR_MOD (XMIN, YMIN), CONSTANTS_MOD (CLN) and IOUNITS_MOD (ISHORT, for the messages), and takes JSON from DIMENSIONS_MOD. The unused COMMON blocks of rp_1 (/GRAVST/, /GRAVS3/, /SCRTCH/, /GSCOF/, /CMIXING/) are no longer declared.

The method

The diffusion region. If the centre is convective (LC(1)), IBEGIN is the last point of the convective core; if every point is convective, settling is skipped (LALLCZ = .TRUE.). IBEGIN then moves out to the first point with $X >$ XMIN. If the surface is convective, IEND is the innermost point of the surface convection zone, otherwise the last point M; it then moves inward to the outermost point with $Y >$ YMIN. A model with $X \le$ XMIN or $Y \le$ YMIN everywhere is skipped.

Question 1 (line 101 ↗; line 123 ↗). Overshoot. Two comments say "COMPUTE OVERSHOOT (TO BE ADDED)," so settling ignores overshoot at the convective boundaries: IBEGIN and IEND are the convective boundaries themselves. Is this still planned?

The second exit's comment is copied from the first. The exit taken when no helium is above YMIN (lines 133–137 ↗) is commented "HYDROGEN-FREE MODEL," as the hydrogen exit is (line 110 ↗); its message says "Y BELOW," so the comment should read "HELIUM-FREE MODEL."
Units. With $R_\odot = 6.9598\times10^{10}$ cm (CRSUN_BAH) and the year $3.1558\times10^{7}$ s (CSECYR_BAH): $$r' = \frac{r}{R_\odot},\quad T' = \frac{T}{10^7\ \mathrm{K}},\quad m' = \frac{10^{-2}\,m}{R_\odot^{3}},\quad \left(\frac{d\ln P}{dr}\right)' = R_\odot\,\frac{d\ln P}{dr},\quad \Delta t' = \frac{\Delta t}{2.7\times10^{13}\ \mathrm{yr}}.$$ That is, the mass unit is $100\ \mathrm{g\,cm^{-3}}\times R_\odot^3$. The comment says the time unit "includes factor of 2.2 from ln Lambda," the Coulomb logarithm that MICRODIFF_COD divides by.

Question 2 (lines 86–87 ↗; lines 144–148 ↗). Constants of its own. The routine uses its own $R_\odot$ and year (CRSUN_BAH, CSECYR_BAH) rather than CRSUN and CSECYR from CONSTANTS_MOD, which it uses (for CLN). Is this intended, to keep Bahcall and Loeb's values?

Design

HS1, HQPR, DT and STOT are converted in place, and MICRODIFF_ETM converts them back; since MICRODIFF passes copies of them (microdiff.f, lines 112–117 ↗), the caller's own values are not touched. When settling is skipped, the routine returns before converting, so the two stay consistent.

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. (cannot be reached) No check that `IEND > IBEGIN`.

(lines 107–116 ↗; lines 130–139 ↗) After the two searches, nothing checks that the region from IBEGIN to IEND is not empty. The search for $X >$ XMIN runs outward to M, past IEND, so if the hydrogen-exhausted centre reached the base of the convective envelope, IBEGIN would end inside the envelope, beyond IEND (but settling stops before the centre is exhausted: below). What follows in MICRODIFF_MTE: the total width DRTOT is zero or negative and is divided by the minimum spacing, which is the same number (microdiff_mte.f, lines 54–63 ↗); with IEND = IBEGIN this is 0/0 and DR = 0, and MICRODIFF_RUN divides by DR (microdiff_run.f, line 145 ↗); with IEND < IBEGIN the grid runs inward, and MICRODIFF_ETM applies both the central and the envelope change to the shells between them. Only with LNEWDIF and LDIFY, LDIFZ or LDIFLI (all off by default; the test suite's settling decks set them).

A possible fix:

microdiff_setup.f

@@ lines 135–142 @@
    17 FORMAT(1X,'Y BELOW ',F9.6,' IN WHOLE MODEL-NO SETTLING')
       LALLCZ = .TRUE.
       GOTO 9999
+C 11/09/26 KC via Claude
+C code changed to skip settling when the region between IBEGIN and
+C IEND is empty, as for a fully convective model.
    47 CONTINUE
       IEND = I
+      IF(IEND.LE.IBEGIN)THEN
+         WRITE(ISHORT,18)
+   18    FORMAT(1X,' NO RADIATIVE REGION LEFT - NO SETTLING')
+         LALLCZ = .TRUE.
+         GOTO 9999
+      ENDIF
 C     CON_MASS=CONVERSION FACTOR FOR MASS.
 C     CON_RADIUS=CONVERSION FACTOR FOR RADIUS.
 C     CON_TEMP=CONVERSION FACOTR FOR TEMPERATURE.

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): The check point was passed at every settling call in all seven runs with settling (testsuite/Test_solarGS98_dif, testsuite/Test_solarAAG21_dif, testsuite/Test_solarM22M_difkur, testsuite/Test_solarGS98_difrotmixper18p4d1p0, testsuite/Test_m0p8feh+0p0GS98_rotp8p0d0p1mag1e4std_TAMS, testsuite/Test_m1p4feh+0p0GS98_dif_TAMS, and testsuite/Test_solarAAG21_dif with a custom mixture), 690 to 3,703 times per run, with IBEGIN from 1 to 536 and IEND from 3 to 972, and IEND was above IBEGIN every time. It cannot be reached as described at 7be4133: MIX and ROTMIX switch settling off for the rest of the run once central $X$ falls below XMIN (mix.f, lines 412–416 ↗; rotmix.f, lines 151–155 ↗), so IBEGIN stays at 1 or at the edge of a mixed core, and the searches give IEND ≥ IBEGIN + 2 (lines 91–127 ↗); only the helium check (lines 130–139 ↗) could move IEND lower, with $Y$ below YMIN outside the core. It is fragile code. The fix was not run. Intent: none found.

Definitions

subroutineMICRODIFF_SETUP

lines 4–169 · GitHub

Uses modules CONFAC_MOD, CONSTANTS_MOD, DIMENSIONS_MOD, IOUNITS_MOD, NMLPAR_MOD

Called by MICRODIFF

Change log

The changes since rp_1

Changed since rp_1. Declarations only: COMMON blocks to modules (the unused ones dropped), IMPLICIT NONE.

1CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 2C MICRODIFF_SETUP 3CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 4 SUBROUTINE MICRODIFF_SETUP(DT,HQPR,HR,HD,HS1,HT,LC,M,STOT,HCOMP,HRU,HTU, 5 * * IBEGIN,IEND,LALLCZ,ODEN,OTEM) 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/CONFAC/ into CONFAC_MOD. 12 USE CONFAC_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/14/26 MHP via Claude 20C code changed to convert COMMON/IOUNITS/ into IOUNITS_MOD. 21 USE IOUNITS_MOD 22C 10/20/26 MHP via Claude 23C code changed to replace implicit typing with IMPLICIT NONE and 24C explicit declarations. 25 IMPLICIT NONE 26 REAL*8 CRSUN_BAH,CSECYR_BAH,DT,HCOMP,HD,HQPR,HR,HRU,HS1,HT,HTU, 27 * * ODEN,OTEM,STOT 28 INTEGER I,IBEGIN,IEND,M 29 LOGICAL*4 LALLCZ,LC 30 31 32C 9/25/26 MHP via Claude 33C code changed to consolidate all I/O logical unit numbers, previously 34C scattered across COMMON/LUOUT/, COMMON/LUNUM/, COMMON/ALEX06/, 35C COMMON/ACDPTH/, COMMON/NEWOPAC/, COMMON/LOPAL95/, COMMON/ATMOS2/, 36C COMMON/ALATM03/, COMMON/OPALEOS/, and COMMON/SCV2/, into a single 37C COMMON/IOUNITS/ 38C 9/25/26 MHP via Claude 39C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/, 40C and COMMON/CONST3/ into a single COMMON/CONSTANTS/ 41C 9/24/26 MHP via Claude 42C code changed to remove dead CMIXL from COMMON/CONST3/; it is now 43C passed as an explicit argument where actually used 44C 10/12/26 MHP via Claude 45C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 46C namelist variables that the physics routines read from 27 COMMON 47C blocks which also held derived or working members. Those blocks keep 48C only their non-namelist members; five of them -- NEWMX, VMULT, 49C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 50C Derivation-only namelist variables from the same blocks went to 51C COMMON/NMLBYPASS/ instead, not here. 52C 53C Declared identically in all 121 files that hold any of these members. 54C Generated from one member list, so every copy is byte-identical; the 55C linker's 'size of symbol changed' warning is the check that they are. 56C 10/13/26 MHP via Claude: removed the unused COMMON/SCRTCH/ declaration 57C (no member referenced here; residue of the obsolete-option removals). 58C MHP 8/94 ADDED I/O FOR DIFFUSION 59C 10/13/26 MHP via Claude 60C code changed to remove COMMON/GSCOF/: none of its four members is 61C referenced in executable code anywhere. 62C COMMON/GSCOF2/TAPP(JSON),TATP(JSON),TCLP(JSON),TAPZP(JSON), 63C * TATZP(JSON) 64CFD 10/09 ADDED COMMON BLOCK FOR EXTRA MIXING. IT MIMIC SOME MIXING BY AFFECTING 65C THE SETTLING COEFFICIENT DIRECTLY (in setup_grsett.f) 66C 67 DIMENSION HQPR(JSON),HR(JSON),HT(JSON),HS1(JSON),HD(JSON), 68 * * LC(JSON),HCOMP(15,JSON),HRU(JSON),HTU(JSON), 69 * * ODEN(JSON),OTEM(JSON) 70 SAVE 71C 72C OUTPUT VARIABLES : 73C 74C HRU - VECTOR OF UNLOGGED RADII IN BAHCALL AND LOEB UNITS. 75C HTU - VECTOR OF UNLOGGED TEMPERATURES IN BAHCALL AND LOEB UNITS. 76C IBEGIN - FIRST ZONE FOR DIFFUSION PURPOSES (EITHER THE FIRST MODEL POINT 77C OR THE OUTERMOST POINT OF A CENTRAL CONVECTION ZONE). 78C IEND - LAST ZONE FOR DIFFUSION PURPOSES (EITHER THE LAST MODEL POINT OR 79C THE INNERMOST POINT OF A SURFACE CONVECTION ZONE). 80C THE VECTORS HQPR AND HS1, AND THE SCALARS DT AND STOT, ARE ALSO CONVERTED 81C TO BAHCALL AND LOEB UNITS. 82C CONSTANTS DEFINED : 83C CLN = CONVERSION FACTOR FROM LN TO LOG10 84C CRSUN = SOLAR RADIUS (CM) 85C CSECYR = NUMBER OF SECONDS IN A YEAR. 86 CRSUN_BAH=6.9598D10 87 CSECYR_BAH=3.1558D7 88C LALLCZ=T FOR FULLY CONVECTIVE MODEL(AND IF TRUE, DIFFUSION IS SKIPPED). 89 LALLCZ=.FALSE. 90C CHECK FOR CONVECTIVE CORE. 91 IF(LC(1))THEN 92 DO 10 I=2,M 93 IF(.NOT.LC(I))GOTO 20 94 10 10 CONTINUE 95C DIFFUSION NOT COMPUTED FOR FULLY CONVECTIVE MODELS. 96 LALLCZ=.TRUE. 97 WRITE(ISHORT,15) 98 15 15 FORMAT(1X,' FULLY CONVECTIVE MODEL - NO SETTLING') 99 GOTO 9999 100 20 20 CONTINUE 101C COMPUTE OVERSHOOT (TO BE ADDED). 102 IBEGIN = I-1 103 ELSE 104 IBEGIN = 1 105 ENDIF 106C MHP 6/90 CHECK FOR HYDROGEN-EXHAUSTED CORE. 107 DO 23 I = IBEGIN,M 108 IF(HCOMP(1,I).GT.XMIN)GOTO 25 109 23 23 CONTINUE 110C HYDROGEN-FREE MODEL - EXIT. 111 WRITE(ISHORT,16)XMIN 112 16 16 FORMAT(1X,'X BELOW ',F9.6,' IN WHOLE MODEL-NO SETTLING') 113 LALLCZ = .TRUE. 114 GOTO 9999 115 25 25 CONTINUE 116 IBEGIN = I 117C CHECK FOR CONVECTIVE ENVELOPE. 118 IF(LC(M))THEN 119 DO 30 I=M-1,2,-1 120 IF(.NOT.LC(I))GOTO 40 121 30 30 CONTINUE 122 40 40 CONTINUE 123C COMPUTE OVERSHOOT (TO BE ADDED). 124 IEND = I+1 125 ELSE 126 IEND = M 127 ENDIF 128C CHECK FOR HELIUM-EXHAUSTED SURFACE. 129C OUTER POINT IS SET WHEREVER Y>YMIN. 130 DO 45 I=IEND,1,-1 131 IF(HCOMP(2,I).GT.YMIN) GOTO 47 132 45 45 CONTINUE 133C HYDROGEN-FREE MODEL - EXIT. 134 WRITE(ISHORT,17)YMIN 135 17 17 FORMAT(1X,'Y BELOW ',F9.6,' IN WHOLE MODEL-NO SETTLING') 136 LALLCZ = .TRUE. 137 GOTO 9999 138 47 47 CONTINUE 139 IEND = I 140C CON_MASS=CONVERSION FACTOR FOR MASS. 141C CON_RADIUS=CONVERSION FACTOR FOR RADIUS. 142C CON_TEMP=CONVERSION FACOTR FOR TEMPERATURE. 143C CON_TIME=CONVERSION FACTOR FOR TIME. 144 CON_RAD=1.0D0/CRSUN_BAH 145 CON_MASS=1.0D-2*CON_RAD**3 146 CON_TEMP=1.0D-7 147C INCLUDES FACTOR OF 2.2 FROM LN LAMBDA 148 CON_TIME=2.7D13*CSECYR_BAH 149C CONVERT LOG(RADIUS) AND LOG(TEMPERATURE) TO NATURAL UNITS. 150C ALSO CONVERT NATURAL UNITS TO BAHCALL AND LOEB UNITS. 151 DO 50 I=1,M 152 HRU(I)=EXP(CLN*HR(I))*CON_RAD 153 HTU(I)=EXP(CLN*HT(I))*CON_TEMP 154 HS1(I)=HS1(I)*CON_MASS 155 HQPR(I)=HQPR(I)/CON_RAD 156C SDEL(2,I)=0.4D0 !COMMENT OUT IN REAL CODE 157 50 50 CONTINUE 158 DT=DT/CON_TIME 159 STOT=STOT*CON_MASS 160C 161C COLLECT THE NECESSARY QUANTITIES (NAMELY RHO AND T) FOR LATER 162C TRANSFORMATION TO THE EQUALLY SPACED GRID. 163 DO I = 1,M 164 ODEN(I) = EXP(CLN*HD(I)) 165 OTEM(I) = EXP(CLN*HT(I)) 166 ENDDO 167 9999 9999 CONTINUE 168 RETURN 169 END

1"""Python version of microdiff_setup.f (YREC 7be4133): the settling 2region and the conversion to Bahcall and Loeb's units, for MICRODIFF.""" 3 4import math 5 6import numpy as np 7 8from pyrec import confac_mod 9from pyrec import constants_mod 10from pyrec import iounits_mod 11from pyrec import nmlpar_mod 12 13 14def microdiff_setup(dt: float, hqpr: np.ndarray, hr: np.ndarray, hd: np.ndarray, hs1: np.ndarray, ht: np.ndarray, 15 lc: np.ndarray, m: int, stot: float, hcomp: np.ndarray, hru: np.ndarray, htu: np.ndarray, 16 ibegin: int, iend: int, lallcz: bool, oden: np.ndarray, otem: np.ndarray) -> tuple[float, float, 17 int, int, bool]: 18 """Find the region that takes part in settling (IBEGIN to IEND) and 19 convert the model to Bahcall and Loeb's units; skip settling 20 (lallcz = True) for a fully convective model, or one with X <= XMIN or 21 Y <= YMIN everywhere. 22 23 Translates: microdiff_setup.f#MICRODIFF_SETUP (7be4133) 24 25 Arguments: the step dt (s); the model's hqpr (dlnP/dr), hr (log r), 26 hd (log rho), hs1 (mass, g), ht (log T), lc (convective flags) and 27 hcomp[k, j], with m shells; the total mass stot (g). ibegin, iend and 28 lallcz are passed in with any value. Arrays have index 0 unused. 29 30 Sets in place: hru, htu (unlogged r and T in the new units), oden, otem 31 (unlogged rho and T, cgs), and converts hs1 and hqpr; writes the 32 conversion factors to confac_mod. 33 34 Returns (dt, stot, ibegin, iend, lallcz): dt and stot in the new units 35 (unchanged when settling is skipped), the region, and the skip flag. 36 37 """ 38 39 cln = constants_mod.cln 40 ishort = iounits_mod.ishort 41 42 # Bahcall and Loeb's solar radius and year (Fortran lines 86–89). 43 crsun_bah = 6.9598e10 44 csecyr_bah = 3.1558e7 45 lallcz = False 46 47 # Check for a convective core (Fortran lines 90–105). 48 if lc[1]: 49 i = 2 50 while i <= m: # loop 10 51 if not lc[i]: 52 break # GOTO 20 53 i += 1 54 else: 55 lallcz = True 56 iounits_mod.write(ishort, " FULLY CONVECTIVE MODEL - NO SETTLING") 57 return dt, stot, ibegin, iend, lallcz # GOTO 9999 58 # Compute overshoot (to be added). 59 ibegin = i - 1 60 else: 61 ibegin = 1 62 63 # Check for a hydrogen-exhausted core, MHP 6/90 (Fortran lines 106–116). 64 for i in range(ibegin, m + 1): # loop 23 65 if hcomp[1, i] > nmlpar_mod.xmin: 66 break # GOTO 25 67 else: 68 # Question (see the page): "HYDROGEN-FREE MODEL" comment, both exits. 69 iounits_mod.write(ishort, f" X BELOW {nmlpar_mod.xmin:9.6f}" 70 " IN WHOLE MODEL-NO SETTLING") 71 lallcz = True 72 return dt, stot, ibegin, iend, lallcz 73 ibegin = i # label 25 74 75 # Check for a convective envelope (Fortran lines 117–127). 76 if lc[m]: 77 i = m - 1 78 while i >= 2: # loop 30 79 if not lc[i]: 80 break # GOTO 40 81 i -= 1 82 # On completion the Fortran's I is 1. 83 iend = i + 1 # label 40 84 else: 85 iend = m 86 87 # Check for a helium-exhausted surface: the outer point is the last 88 # one with Y > YMIN (Fortran lines 128–139). 89 for i in range(iend, 0, -1): # loop 45 90 if hcomp[2, i] > nmlpar_mod.ymin: 91 break # GOTO 47 92 else: 93 iounits_mod.write(ishort, f" Y BELOW {nmlpar_mod.ymin:9.6f}" 94 " IN WHOLE MODEL-NO SETTLING") 95 lallcz = True 96 return dt, stot, ibegin, iend, lallcz 97 iend = i # label 47 98 # Suspected bug (see the page): no check that iend > ibegin. 99 100 # Conversion factors; the time unit includes a factor of 2.2 from 101 # ln Lambda (Fortran lines 140–148). 102 confac_mod.con_rad = 1.0 / crsun_bah 103 confac_mod.con_mass = 1.0e-2 * confac_mod.con_rad**3 104 confac_mod.con_temp = 1.0e-7 105 confac_mod.con_time = 2.7e13 * csecyr_bah 106 107 # Convert to Bahcall and Loeb's units (Fortran lines 149–159). 108 for i in range(1, m + 1): # loop 50 109 hru[i] = math.exp(cln * hr[i]) * confac_mod.con_rad 110 htu[i] = math.exp(cln * ht[i]) * confac_mod.con_temp 111 hs1[i] = hs1[i] * confac_mod.con_mass 112 hqpr[i] = hqpr[i] / confac_mod.con_rad 113 dt = dt / confac_mod.con_time 114 stot = stot * confac_mod.con_mass 115 116 # Rho and T for the later transformation to the equal grid (Fortran 117 # lines 160–166). 118 for i in range(1, m + 1): 119 oden[i] = math.exp(cln * hd[i]) 120 otem[i] = math.exp(cln * ht[i]) 121 122 return dt, stot, ibegin, iend, lallcz

microdiff_setup.f

Purpose

How it fits

The method

Design

Suspected bugs, checked by running the code

1. (cannot be reached) No check that IEND > IBEGIN.

Definitions

subroutineMICRODIFF_SETUP

Change log

1. (cannot be reached) No check that `IEND > IBEGIN`.