microdiff_etm.f

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

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

Variables on this page (7)

con_mass, con_rad, con_temp, con_time, json, ldifli, ldifz

Purpose

"Equal to model": the last step of MICRODIFF. It interpolates the changes in abundance from the equally spaced grid back to the model points, applies them to HCOMP, recomputes helium, and converts the variables changed by MICRODIFF_SETUP back to cgs units.

How it fits

Inputs: the equal grid radii ER, the changes EDELX, EDELZ, EDELLI from MICRODIFF_RUN, IBEGIN, IEND, NPT, and the conversion factors in CONFAC_MOD (written by MICRODIFF_SETUP).
Who calls it: only MICRODIFF.
What it calls: INTRP2, for four-point Lagrange weights.
Outputs: HCOMP updated in place; HRU, HTU, HS1, HQPR, DT and STOT converted back (lines 165–172 ↗), though MICRODIFF passes its local HRU and HTU and copies of the other four, and returns right after the call (microdiff.f, lines 112–117 ↗; microdiff.f, lines 242–244 ↗): unused?
Modules: CONFAC_MOD, NMLPAR_MOD (LDIFZ, LDIFLI), DIMENSIONS_MOD (JSON).

The method

Central region (points 1..IBEGIN, the convective core or the hydrogen-exhausted centre): the change at the first grid point is added to each point: $X \to \max(X + \Delta X_1, 0)$, the light elements likewise floored at zero, and $Z \to \min(Z + \Delta Z_1,\ 1 - X - X_3)$.
Radiative interior (points IBEGIN+1..IEND-1): the change is interpolated with four-point Lagrange weights $w_k$ over the grid points around the model radius, $\Delta X(r) = \sum_k w_k\,\Delta X_{k_0+k-1}$, and $X \to \min(X + \Delta X, 1 - Z - X_3)$.
Envelope (points IEND..M): the change at the last grid point is added, with $X$ capped as above and $Z$ floored at zero.

Question 1 (lines 47–48 ↗; lines 104–111 ↗; lines 135–142 ↗). The guard comment. The comment says the failsafes guard "against negative X and Z>1 respectively," and so they do in the centre ($X \ge 0$, $Z \le 1 - X - X_3$); but in the interior $X$ and $Z$ are both capped from above, and in the envelope $X$ is capped from above and $Z$ floored at zero. Which guards are intended where?
Metals are scaled as a whole: with $Z' = Z + \Delta Z$ (capped at $1 - X - X_3$ in the centre and the interior), the individual metal rows 5–11 of HCOMP (¹²C, ¹³C, ¹⁴N, ¹⁵N, ¹⁶O, ¹⁷O, ¹⁸O, as STARIN names them) are multiplied by $Z'/Z$.
Helium follows from closure: $Y = 1 - X - Z - X_3$, where $X_3$ is ³He, HCOMP(4,·) (STARIN fills it from XHE3_INI).

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. (harmless in the runs) Light elements in the interior and the envelope are not floored at zero.

(lines 66–74 ↗; lines 122–131 ↗; lines 156–163 ↗) The light elements (HCOMP(13..15,·)) are floored at zero only in the centre (MAX(…,XMINIM)); in the interior (HCOMP(II,I) + DXMOD, with DXMOD from four-point Lagrange interpolation, which can overshoot) and in the envelope the change is added with no floor. What follows: where an element is nearly depleted, its abundance can become negative until MIX and ROTMIX floor every species at zero after the settling substeps (mix.f, line 470 ↗; rotmix.f, line 244 ↗), so only the later substeps of the same call see it. Only with LNEWDIF and LDIFLI (both off by default; the test suite's settling decks set them).

A possible fix:

microdiff_etm.f

@@ lines 126–132 @@
      *                 FACINTERP(2)*EDELLI(KK,K0+1)+
      *                 FACINTERP(3)*EDELLI(KK,K0+2)+
      *                 FACINTERP(4)*EDELLI(KK,K0+3)
-               HCOMP(II,I) = HCOMP(II,I) + DXMOD
+C 11/09/26 KC via Claude
+C code changed to floor the light elements at zero in the interior
+C and the envelope too, as in the centre.
+               HCOMP(II,I) = MAX(HCOMP(II,I) + DXMOD,0.0D0)
             ENDDO
          ENDIF
    20 CONTINUE
@@ lines 157–163 @@
          DO KK = 1,NLIGHT
             DO 65 I = IEND,M
                II = CSLIGHT(KK)
-               HCOMP(II,I)=HCOMP(II,I) + EDELLI(KK,NPT)
+               HCOMP(II,I)=MAX(HCOMP(II,I) + EDELLI(KK,NPT),0.0D0)
    65       CONTINUE
          ENDDO
       ENDIF

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): taken 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), up to 13.2 million times per run (testsuite/Test_solarM22M_difkur). A light element became negative 7,781 to 92,130 times per run, by at most $3.8\times10^{-13}$. With the fix, testsuite/Test_m1p4feh+0p0GS98_dif_TAMS and testsuite/Test_solarM22M_difkur give identical tracks and stored models, with the same numbers of models and ages, since MIX floors the values at zero anyway (above). The negative values never reach the model, and the code relies on that later floor. Intent: The comment guards X and Z only ("GUARDING AGAINST NEGATIVE X AND Z>1," lines 46–48 ↗), and the light elements have a floor only in the centre (lines 68–71 ↗).

2. (never occurred) Division by zero where $Z = 0$.

(line 53 ↗; line 112 ↗; line 143 ↗) ZZ2 = ZZ/HCOMP(3,I) divides by zero at a point with $Z = 0$; with IEEE arithmetic it gives an infinity or NaN rather than a stop, and the metal rows 5–11 (zero there) become NaN. Only with LNEWDIF and LDIFZ and a model with $Z = 0$ somewhere (a metal-free model).

A possible fix (the same at the three places; where $Z = 0$ the rows 5–11 stay zero):

microdiff_etm.f

@@ lines 50–56 @@
          DO I = IBEGIN,1,-1
             ZMAX = 1.0D0 - HCOMP(1,I) - HCOMP(4,I)
             ZZ=MIN(HCOMP(3,I)+EDELZ(1),ZMAX)
-            ZZ2 = ZZ/HCOMP(3,I)
+C 11/09/26 KC via Claude
+C code changed to avoid a division by zero where Z = 0.
+            IF(HCOMP(3,I).GT.0.0D0)THEN
+               ZZ2 = ZZ/HCOMP(3,I)
+            ELSE
+               ZZ2 = 1.0D0
+            ENDIF
             HCOMP(3,I) = ZZ
             DO J = 5,11
                HCOMP(J,I) = ZZ2*HCOMP(J,I)
@@ lines 109–115 @@
             DZMOD = FACINTERP(1)*EDELZ(K0)+FACINTERP(2)*EDELZ(K0+1)+
      *             FACINTERP(3)*EDELZ(K0+2)+FACINTERP(4)*EDELZ(K0+3)
             ZZ = MIN(HCOMP(3,I)+DZMOD,ZMAX)
-            ZZ2 = ZZ/HCOMP(3,I)
+C 11/09/26 KC via Claude
+C code changed to avoid a division by zero where Z = 0.
+            IF(HCOMP(3,I).GT.0.0D0)THEN
+               ZZ2 = ZZ/HCOMP(3,I)
+            ELSE
+               ZZ2 = 1.0D0
+            ENDIF
             HCOMP(3,I)=ZZ
             DO J = 5,11
                HCOMP(J,I) = ZZ2*HCOMP(J,I)
@@ lines 140–146 @@
          ZMIN = 0.0D0
          DO I = IEND,M
             ZZ = MAX(HCOMP(3,I)+EDELZ(NPT),ZMIN)
-            ZZ2 = ZZ/HCOMP(3,I)
+C 11/09/26 KC via Claude
+C code changed to avoid a division by zero where Z = 0.
+            IF(HCOMP(3,I).GT.0.0D0)THEN
+               ZZ2 = ZZ/HCOMP(3,I)
+            ELSE
+               ZZ2 = 1.0D0
+            ENDIF
             HCOMP(3,I) = ZZ
             DO J = 5,11
                HCOMP(J,I) = ZZ2*HCOMP(J,I)

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): taken in all seven runs with settling (six of the test suite's and a made-up one with a custom mixture, listed in suspected bug 1), up to 4.7 million times per run; $Z$ was never zero (at least 0.0020). Not reached: It needs LDIFZ and a model with $Z = 0$ somewhere, and no case has a metal-free model. The fix was not run. Intent: none found.

Definitions

subroutineMICRODIFF_ETM

lines 2–174 · GitHub

Uses modules CONFAC_MOD, DIMENSIONS_MOD, NMLPAR_MOD

Called by MICRODIFF

Flowchart of its calls, in order (1)

for I=IBEGIN+1,IEND-1 (DO 20), lines 77–133

INTRP2 100

Change log

The changes since rp_1

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

Changed since a8a5950. The file is unchanged, but MICRODIFF_MTE now gives the grid at least 4 points, so K0 stays within the grid; this settles former question 5.

1C $$$$$$ 2 SUBROUTINE MICRODIFF_ETM(DT,ER,EDELX,EDELZ,EDELLI,IBEGIN,IEND,NPT, 3 * * HCOMP,HQPR,HRU,HS1,HTU,M,STOT,NLIGHT,CSLIGHT) 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/CONFAC/ into CONFAC_MOD. 10 USE CONFAC_MOD 11C 10/14/26 MHP via Claude 12C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD. 13 USE NMLPAR_MOD 14C 10/20/26 MHP via Claude 15C code changed to replace implicit typing with IMPLICIT NONE and 16C explicit declarations. 17 IMPLICIT NONE 18 REAL*8 DT,DXMOD,DZMOD,EDELLI,EDELX,EDELZ,ER,FACINTERP,HCOMP,HQPR, 19 * * HRU,HS1,HTU,RADMOD,STOT,TABLER,XMAX,XMINIM,ZMAX,ZMIN,ZZ,ZZ2 20 INTEGER I,IBEGIN,IEND,II,J,JMIN,K,K0,KK,M,NLIGHT,NPT 21C 10/12/26 MHP via Claude 22C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 23C namelist variables that the physics routines read from 27 COMMON 24C blocks which also held derived or working members. Those blocks keep 25C only their non-namelist members; five of them -- NEWMX, VMULT, 26C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 27C Derivation-only namelist variables from the same blocks went to 28C COMMON/NMLBYPASS/ instead, not here. 29C 30C Declared identically in all 121 files that hold any of these members. 31C Generated from one member list, so every copy is byte-identical; the 32C linker's 'size of symbol changed' warning is the check that they are. 33 DIMENSION ER(JSON),EDELX(JSON),EDELZ(JSON),EDELLI(NLIGHT,JSON), 34 * * HCOMP(15,JSON),HQPR(JSON),HRU(JSON),HS1(JSON), 35C KC 2025-05-30 fixed -Wconversion 36C * HTU(JSON),TABLER(4),FACINTERP(4),CSLIGHT(NLIGHT) 37 * * HTU(JSON),TABLER(4),FACINTERP(4) 38 INTEGER CSLIGHT(NLIGHT) 39 SAVE 40C TRANSFORM BACK TO ORIGINAL GRID OF MODEL POINTS FROM EQUALLY 41C SPACED GRID. 42 XMINIM = 0.0D0 43 DO 10 I = IBEGIN,1,-1 44 HCOMP(1,I)=MAX(HCOMP(1,I) + EDELX(1),XMINIM) 45 10 10 CONTINUE 46C MHP 3/94 ADDED METAL DIFFUSION 47C NOTE THAT BECAUSE METALS SINK, AND HYDROGEN RISES, THE FAILSAFES 48C ARE OPPOSITE (GUARDING AGAINST NEGATIVE X AND Z>1 RESPECTIVELY). 49 IF(LDIFZ)THEN 50 DO I = IBEGIN,1,-1 51 ZMAX = 1.0D0 - HCOMP(1,I) - HCOMP(4,I) 52 ZZ=MIN(HCOMP(3,I)+EDELZ(1),ZMAX) 53 ZZ2 = ZZ/HCOMP(3,I) 54 HCOMP(3,I) = ZZ 55 DO J = 5,11 56 HCOMP(J,I) = ZZ2*HCOMP(J,I) 57 END DO 58 HCOMP(2,I)=1.0D0-HCOMP(1,I)-HCOMP(3,I)-HCOMP(4,I) 59 END DO 60 ELSE 61 DO I = IBEGIN,1,-1 62 HCOMP(2,I)=1.0D0-HCOMP(1,I)-HCOMP(3,I)-HCOMP(4,I) 63 END DO 64 ENDIF 65C G SOMERS 5/15; ADD LIGHT ELEMENT DIFFUSION 66 IF(LDIFLI)THEN 67 DO KK = 1,NLIGHT 68 XMINIM = 0.0D0 69 DO 15 I = IBEGIN,1,-1 70 II = CSLIGHT(KK) 71 HCOMP(II,I)=MAX(HCOMP(II,I) + EDELLI(KK,1),XMINIM) 72 15 15 CONTINUE 73 END DO 74 ENDIF 75C 76 JMIN=2 77 DO 20 I=IBEGIN+1,IEND-1 78 DO 30 J=JMIN,NPT 79C FIND EQUALLY SPACED GRID POINTS CLOSEST TO THE MODEL POINT. 80 IF(ER(J).GE.HRU(I))THEN 81C ENSURE THAT FIRST INTERP. POINT NO LESS THAN FIRST EQUALLY SPACED POINT. 82 K0 = MAX(J-2,1) 83C ENSURE THAT LAST INTERP. POINT NO GREATER THAN LAST EQUALLY SPACED POINT. 84 K0 = MIN(K0,NPT-3) 85C JVS fix for NPT = 3? 86 IF (K0 .EQ. 0) K0=1 87 JMIN=J 88 GOTO 40 89 ENDIF 90 30 30 CONTINUE 91 K0 = NPT-3 92 JMIN=NPT 93 40 40 CONTINUE 94 95 DO 50 K=1,4 96 TABLER(K)=ER(K0+K-1) 97 50 50 CONTINUE 98 RADMOD=HRU(I) 99C FIND 4 POINT LAGRANGIAN INTERPOLATION FACTORS. 100 CALL INTRP2(TABLER,FACINTERP,RADMOD) 101C PERFORM 4 POINT LAGRANGIAN INTERPOLATION FOR CHANGE IN X. 102 DXMOD = FACINTERP(1)*EDELX(K0)+FACINTERP(2)*EDELX(K0+1)+ 103 * * FACINTERP(3)*EDELX(K0+2)+FACINTERP(4)*EDELX(K0+3) 104 XMAX = 1.0D0 - HCOMP(3,I) - HCOMP(4,I) 105 HCOMP(1,I)=MIN(HCOMP(1,I) + DXMOD,XMAX) 106C MHP 3/94 ADDED METAL DIFFUSION 107 IF(LDIFZ)THEN 108 ZMAX = 1.0D0 - HCOMP(1,I) - HCOMP(4,I) 109 DZMOD = FACINTERP(1)*EDELZ(K0)+FACINTERP(2)*EDELZ(K0+1)+ 110 * * FACINTERP(3)*EDELZ(K0+2)+FACINTERP(4)*EDELZ(K0+3) 111 ZZ = MIN(HCOMP(3,I)+DZMOD,ZMAX) 112 ZZ2 = ZZ/HCOMP(3,I) 113 HCOMP(3,I)=ZZ 114 DO J = 5,11 115 HCOMP(J,I) = ZZ2*HCOMP(J,I) 116 ENDDO 117 HCOMP(2,I)=1.0D0-HCOMP(1,I)-HCOMP(3,I)-HCOMP(4,I) 118 ELSE 119 HCOMP(2,I)=1.0D0-HCOMP(1,I)-HCOMP(3,I)-HCOMP(4,I) 120 ENDIF 121C GES 5/15 ADDED LIGHT ELEMENT DIFFUSION 122 IF(LDIFLI)THEN 123 DO KK = 1,NLIGHT 124 II = CSLIGHT(KK) 125 DXMOD = FACINTERP(1)*EDELLI(KK,K0)+ 126 * * FACINTERP(2)*EDELLI(KK,K0+1)+ 127 * * FACINTERP(3)*EDELLI(KK,K0+2)+ 128 * * FACINTERP(4)*EDELLI(KK,K0+3) 129 HCOMP(II,I) = HCOMP(II,I) + DXMOD 130 ENDDO 131 ENDIF 132 20 20 CONTINUE 133C 134 DO 60 I = IEND,M 135 XMAX = 1.0D0 - HCOMP(3,I) - HCOMP(4,I) 136 HCOMP(1,I)=MIN(HCOMP(1,I) + EDELX(NPT),XMAX) 137 60 60 CONTINUE 138C MHP 3/94 ADDED METAL DIFFUSION 139 IF(LDIFZ)THEN 140 ZMIN = 0.0D0 141 DO I = IEND,M 142 ZZ = MAX(HCOMP(3,I)+EDELZ(NPT),ZMIN) 143 ZZ2 = ZZ/HCOMP(3,I) 144 HCOMP(3,I) = ZZ 145 DO J = 5,11 146 HCOMP(J,I) = ZZ2*HCOMP(J,I) 147 END DO 148 HCOMP(2,I)=1.0D0-HCOMP(1,I)-HCOMP(3,I)-HCOMP(4,I) 149 END DO 150 ELSE 151 DO I = IEND,M 152 HCOMP(2,I)=1.0D0-HCOMP(1,I)-HCOMP(3,I)-HCOMP(4,I) 153 END DO 154 ENDIF 155C GES 5/15 LIGHT ELEMENT DIFFUSION 156 IF(LDIFLI)THEN 157 DO KK = 1,NLIGHT 158 DO 65 I = IEND,M 159 II = CSLIGHT(KK) 160 HCOMP(II,I)=HCOMP(II,I) + EDELLI(KK,NPT) 161 65 65 CONTINUE 162 ENDDO 163 ENDIF 164C 165 DO 70 I=1,M 166 HRU(I)=HRU(I)/CON_RAD 167 HTU(I)=HTU(I)/CON_TEMP 168 HS1(I)=HS1(I)/CON_MASS 169 HQPR(I)=HQPR(I)*CON_RAD 170 70 70 CONTINUE 171 DT=DT*CON_TIME 172 STOT=STOT/CON_MASS 173 RETURN 174 END

1"""Python version of microdiff_etm.f (YREC 7be4133): "equal to model", 2the changes in abundance mapped back from the equal grid to the model.""" 3 4import numpy as np 5 6from pyrec import confac_mod 7from pyrec import nmlpar_mod 8from pyrec.intrp2 import intrp2 9 10# The Fortran's local arrays (SAVE). 11tabler = np.zeros(5) 12facinterp = np.zeros(5) 13 14 15def _rescale_metals(hcomp, i, zz): 16 # HCOMP(3,I) set to ZZ, the metal rows 5..11 scaled by ZZ/Z, and Y by 17 # closure: the same block at the centre, the interior and the envelope 18 # (Fortran lines 53–58, 112–117, 143–148). 19 # Suspected bug (see the page): divides by zero if hcomp[3, i] = 0. 20 zz2 = zz / hcomp[3, i] 21 hcomp[3, i] = zz 22 for j in range(5, 12): 23 hcomp[j, i] = zz2 * hcomp[j, i] 24 hcomp[2, i] = 1.0 - hcomp[1, i] - hcomp[3, i] - hcomp[4, i] 25 26 27def microdiff_etm(dt: float, er: np.ndarray, edelx: np.ndarray, edelz: np.ndarray, edelli: np.ndarray, ibegin: int, 28 iend: int, npt: int, hcomp: np.ndarray, hqpr: np.ndarray, hru: np.ndarray, hs1: np.ndarray, 29 htu: np.ndarray, m: int, stot: float, nlight: int, cslight: np.ndarray) -> tuple[float, float]: 30 """Interpolate the changes in X, Z and the light elements from the 31 equal grid back to the model points, apply them to hcomp (helium by 32 closure), and convert the variables MICRODIFF_SETUP changed back to cgs. 33 34 Translates: microdiff_etm.f#MICRODIFF_ETM (7be4133) 35 36 Arguments: dt (Bahcall and Loeb units); the grid radii er and the 37 changes edelx, edelz, edelli[kk, i] on npt points; the region ibegin, 38 iend; the model hcomp[k, j], hqpr, hru, hs1, htu with m shells; stot; 39 nlight light elements at the hcomp rows cslight[1..nlight]. Arrays 40 have index 0 unused. 41 42 Sets in place: hcomp; hru, htu, hs1, hqpr converted back. 43 44 Returns (dt, stot), converted back to seconds and grams. 45 46 """ 47 48 ldifz = nmlpar_mod.ldifz 49 ldifli = nmlpar_mod.ldifli 50 51 # The centre (points ibegin..1): the change at the first grid point 52 # (Fortran lines 40–45). 53 xminim = 0.0 54 for i in range(ibegin, 0, -1): # loop 10 55 hcomp[1, i] = max(hcomp[1, i] + edelx[1], xminim) 56 57 # Metals and helium at the centre (Fortran lines 46–64). 58 # Question (see the page): which guards are intended where. 59 if ldifz: 60 for i in range(ibegin, 0, -1): 61 zmax = 1.0 - hcomp[1, i] - hcomp[4, i] 62 zz = min(hcomp[3, i] + edelz[1], zmax) 63 _rescale_metals(hcomp, i, zz) 64 else: 65 for i in range(ibegin, 0, -1): 66 hcomp[2, i] = 1.0 - hcomp[1, i] - hcomp[3, i] - hcomp[4, i] 67 68 # Light elements at the centre, floored at zero (Fortran lines 65–74). 69 if ldifli: 70 for kk in range(1, nlight + 1): 71 xminim = 0.0 72 for i in range(ibegin, 0, -1): # loop 15 73 ii = cslight[kk] 74 hcomp[ii, i] = max(hcomp[ii, i] + edelli[kk, 1], xminim) 75 76 # The radiative interior: four-point Lagrange interpolation (Fortran 77 # lines 75–132). 78 jmin = 2 79 for i in range(ibegin + 1, iend): # loop 20 80 # The grid points around the model point (Fortran lines 76–93). 81 for j in range(jmin, npt + 1): # loop 30 82 if er[j] >= hru[i]: 83 k0 = max(j - 2, 1) 84 k0 = min(k0, npt - 3) 85 # Handles the case NPT = 3 (JVS). 86 if k0 == 0: 87 k0 = 1 88 jmin = j 89 break # GOTO 40 90 else: 91 k0 = npt - 3 92 jmin = npt 93 # label 40 94 95 # The interpolation factors, and X (Fortran lines 95–105). 96 for k in range(1, 5): # loop 50 97 tabler[k] = er[k0 + k - 1] 98 radmod = hru[i] 99 intrp2(tabler, facinterp, radmod) 100 dxmod = (facinterp[1] * edelx[k0] + facinterp[2] * edelx[k0 + 1] 101 + facinterp[3] * edelx[k0 + 2] + facinterp[4] * edelx[k0 + 3]) 102 xmax = 1.0 - hcomp[3, i] - hcomp[4, i] 103 hcomp[1, i] = min(hcomp[1, i] + dxmod, xmax) 104 105 # Metals and helium (Fortran lines 106–120). 106 if ldifz: 107 zmax = 1.0 - hcomp[1, i] - hcomp[4, i] 108 dzmod = (facinterp[1] * edelz[k0] + facinterp[2] * edelz[k0 + 1] 109 + facinterp[3] * edelz[k0 + 2] 110 + facinterp[4] * edelz[k0 + 3]) 111 zz = min(hcomp[3, i] + dzmod, zmax) 112 _rescale_metals(hcomp, i, zz) 113 else: 114 hcomp[2, i] = 1.0 - hcomp[1, i] - hcomp[3, i] - hcomp[4, i] 115 116 # Light elements (Fortran lines 121–131). 117 if ldifli: 118 for kk in range(1, nlight + 1): 119 ii = cslight[kk] 120 dxmod = (facinterp[1] * edelli[kk, k0] 121 + facinterp[2] * edelli[kk, k0 + 1] 122 + facinterp[3] * edelli[kk, k0 + 2] 123 + facinterp[4] * edelli[kk, k0 + 3]) 124 # Suspected bug (see the page): no floor at zero here. 125 hcomp[ii, i] = hcomp[ii, i] + dxmod 126 127 # The envelope (points iend..m): the change at the last grid point 128 # (Fortran lines 134–137). 129 for i in range(iend, m + 1): # loop 60 130 xmax = 1.0 - hcomp[3, i] - hcomp[4, i] 131 hcomp[1, i] = min(hcomp[1, i] + edelx[npt], xmax) 132 133 # Metals and helium in the envelope (Fortran lines 138–154). 134 if ldifz: 135 zmin = 0.0 136 for i in range(iend, m + 1): 137 zz = max(hcomp[3, i] + edelz[npt], zmin) 138 _rescale_metals(hcomp, i, zz) 139 else: 140 for i in range(iend, m + 1): 141 hcomp[2, i] = 1.0 - hcomp[1, i] - hcomp[3, i] - hcomp[4, i] 142 143 # Light elements in the envelope (Fortran lines 155–163). 144 # Suspected bug (see the page): no floor at zero here. 145 if ldifli: 146 for kk in range(1, nlight + 1): 147 for i in range(iend, m + 1): # loop 65 148 ii = cslight[kk] 149 hcomp[ii, i] = hcomp[ii, i] + edelli[kk, npt] 150 151 # Back to cgs units; MICRODIFF does not use these (Fortran lines 152 # 165–172). 153 for i in range(1, m + 1): # loop 70 154 hru[i] = hru[i] / confac_mod.con_rad 155 htu[i] = htu[i] / confac_mod.con_temp 156 hs1[i] = hs1[i] / confac_mod.con_mass 157 hqpr[i] = hqpr[i] * confac_mod.con_rad 158 dt = dt * confac_mod.con_time 159 stot = stot / confac_mod.con_mass 160 161 return dt, stot