liburn2.f

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

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

Variables on this page (20)

cln, delam, delrm, hro, jenvos, json, lbliold, linstb, lovste, lrot, pscal0, pscal1, rbase0, rbe9, rbe90, rli6, rli60, rli7, rli70, sdel

Purpose

LIBURN2 is a cheaper version of LIBURN: it burns ⁶Li, ⁷Li and ⁹Be over a step in one go, with the rate taken at the geometric mean of its start and end values, without substeps, extrapolation or accretion. Its banner comment still says "LIBURN" (line 4 ↗).

How it fits

Called by GETW (rotating models, as an alternative to LIBURN) and BURSMIX. Same arguments (DELTS in seconds, as for LIBURN, whose comments' "in years" are stale: LIBURN), and the same modules as LIBURN without DEUTER_MOD and IOUNITS_MOD (it neither accretes nor writes): NEWRAT_MOD, OLDRAT_MOD (rates from LIRATE88), SCRTCH_MOD, MDPHY_MOD, OLDMOD_MOD, LIOV_MOD, PREVCZ_MOD, NMLPAR_MOD (LOVSTE, LROT, LINSTB), CONSTANTS_MOD, DIMENSIONS_MOD. No calls.

The method

The base of the convection zone is found as in LIBURN. Then, in each radiative point below both bases (skipping points where the ⁹Be rate, at the start or the end, is at most $10^{-32}$; destroying everything above $\log T = 7$), and in the convection zone as a whole,

$$x \leftarrow x\,\exp\!\left(-\Delta t\,\sqrt{\lambda_0\lambda_1}\right),$$

with $\lambda_0$, $\lambda_1$ the rates (mass-weighted means in the zone) at the start and end of the step; set to 0 when the exponent is 30 or more, or the result below $10^{-24}$. The zone's start rate is always computed from the start-of-step rates RLI60…, never taken from LIBURN's stored end rates, and LIBURN2 stores no end rates: It clears PREVCZ_MOD's flag LBLIOLD, so that LIBURN does not reuse a rate after it (the code's comment: BURSMIX restarts each extrapolation sequence from the start-of-step model and does not advance RLI60). Points left by a retreating zone get the weighted rate of LIBURN's step 4. A point whose ⁹Be rate is at most $10^{-32}$ is skipped in the radiative loop (line 158 ↗) but ends the last loop and the routine (line 334 ↗), as in both of LIBURN's loops; which is meant is LIBURN's question 1.

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) The sign of `FRAD`.

(lines 319–339 ↗) The same lines as LIBURN's (liburn.f, suspected bug 1): DMASS = HSBEG - HSEND is negative when the zone retreats, so FRAD lies between −1 and 0 and the blended rate extrapolates beyond the zone's rate; and the weights are the other way round from the comment's "a point 1/3 from the old to the new CZ base in mass spent 1/3 of the time in the CZ." Only with LEXCOM (off by default), for rotating models through GETW and BURSMIX, on steps where the zone retreats.

A possible fix: as for LIBURN,

liburn2.f

@@ lines 327–342 @@
          HSBEG = 0.0D0
       ENDIF
       HSEND = 0.5D0*(HS1(JENV)+HS1(JENV-1))
-      DMASS = HSBEG - HSEND
+C 11/09/26 KC via Claude
+C code changed to make FRAD the fraction of the step spent in the CZ,
+C between 0 at the old base and 1 at the new one, as the comment above
+C says, and to give it the weight of the convective rate.
+      DMASS = HSEND - HSBEG
       DO I = JENV0,JENV-1
 C MHP 9/91 CHANGE TO AVOID DIVISION BY ZERO.
 C SKIP IF SHELL TEMPERATURE DROPS BELOW BURNING THRESHOLD.
          IF(RBE9(I).LE.1.0D-32)GOTO 200
          FRAD = (HS1(I)-HSBEG)/DMASS
-C USE FRAD*RADIATIVE RATE AND (1-FRAD)*CONVECTIVE RATE.
-         DLI6 = DELTS*EXP(FRAD*LOG(RLI6(I))+(1.0D0-FRAD)*BLI60)
-         DLI7 = DELTS*EXP(FRAD*LOG(RLI7(I))+(1.0D0-FRAD)*BLI70)
-         DBE9 = DELTS*EXP(FRAD*LOG(RBE9(I))+(1.0D0-FRAD)*BBE90)
+C USE FRAD*CONVECTIVE RATE AND (1-FRAD)*RADIATIVE RATE.
+         DLI6 = DELTS*EXP((1.0D0-FRAD)*LOG(RLI6(I))+FRAD*BLI60)
+         DLI7 = DELTS*EXP((1.0D0-FRAD)*LOG(RLI7(I))+FRAD*BLI70)
+         DBE9 = DELTS*EXP((1.0D0-FRAD)*LOG(RBE9(I))+FRAD*BBE90)
 C***REMEMBER TO ADD FAILSAFES FOR LARGE DEPLETION***
          HCOMP(13,I) = HCOMP(13,I)/EXP(DLI6)
          IF(HCOMP(13,I).LT.1.0D-24)HCOMP(13,I)=0.0D0

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): Not reached in any of the 34 runs, and cannot be reached at 7be4133: LIBURN2 is called only from GETW and from BURSMIX, both only when LBURS is true (getw.f, lines 403–418 ↗; getw.f, lines 456–460 ↗), and GETW sets LBURS = .FALSE. on every call (getw.f, line 140 ↗), which no namelist setting changes. It is code that cannot run today; with LBURS switched on, the sign would matter as in LIBURN (liburn.f, suspected bug 1).

2. (cannot be reached) Points newly in a deepened zone are not burned.

(lines 155–157 ↗; lines 274–278 ↗; lines 308–312 ↗) As in LIBURN (liburn.f, suspected bug 2): When the zone deepens, its starting abundance includes the points JENV to JENV0−1, but the result is written only from MAXJ = JENV0 outward, and the radiative loop stops below JENV, so those points are not burned over the step. Only with LEXCOM, on steps where the zone deepens.

A possible fix:

liburn2.f

@@ lines 305–311 @@
             FBE9 = 0.0D0
          ENDIF
          IF(FBE9.LT.1.0D-24)FBE9 = 0.0D0
-      DO I = MAXJ,M
+C 11/09/26 KC via Claude
+C code changed to give the burned abundances to the whole zone at the
+C end of the step, including points that joined it as it deepened.
+      DO I = JENV,M
          HCOMP(13,I) = FLI6
          HCOMP(14,I) = FLI7
          HCOMP(15,I) = FBE9

3. (cannot be reached) The overshoot zone's width mixes two models.

(lines 135–142 ↗) With envelope overshoot, the offset FX of the new base is measured in a zone width DR = EXP(CLN*HRO(J+1)) - R, the start-of-step radius of the outer point minus the end-of-step radius of the inner one. LIBURN takes both from the end-of-step model since a8a5950 (liburn.f, lines 172–175 ↗). A wrong FX changes the mass and rates given to the bottom point of the zone. Only with LOVSTE (off by default) and LEXCOM.

A possible fix:

liburn2.f

@@ lines 136–142 @@
                R = EXP(CLN*HR(J))
                IF(R.LT.ROS)THEN
                   JENV = J + 1
-                  DR = EXP(CLN*HRO(J+1))-R
+C 11/09/26 KC via Claude
+C code changed to take the zone width from the end-of-step model alone,
+C as in LIBURN.
+                  DR = EXP(CLN*HR(J+1))-R
                   FX = 0.5D0-((ROS-R)/DR)
                   FX = MAX(-0.5D0,FX)
                   FX = MIN(0.5D0,FX)

4. (cannot be reached) `HRO(0)` read for a fully convective start.

(lines 113–115 ↗) With envelope overshoot, the starting base radius reads HRO(JNENV0-1), which is HRO(0), outside the array, when the start-of-step model is fully convective (JNENV0 = 1). The development build's -fcheck=all (src/Makefile line 61) stops the run there; otherwise the radius is computed from whatever lies before HRO in memory. LIBURN guards this case (liburn.f, lines 139–148 ↗). Only with LOVSTE (off by default) and LEXCOM, for a fully convective start.

A possible fix:

liburn2.f

@@ lines 111–118 @@
          ELSE
 C STARTING CZ DEPTH
             IF(RBASE0.EQ.0.0D0)THEN
+C 11/09/26 KC via Claude
+C code changed to take the radius below shell 1 as zero (the centre),
+C as in LIBURN.
+            IF(JNENV0.GT.1)THEN
                RBASE0 = 0.5D0*(EXP(CLN*HRO(JNENV0))
      *                  +EXP(CLN*HRO(JNENV0-1)))
+            ELSE
+               RBASE0 = 0.5D0*EXP(CLN*HRO(1))
+            ENDIF
                ROS = RBASE0 - PSCAL0
                DO J = JNENV0-1,1,-1
                   R = EXP(CLN*HRO(J))

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): Not reached in any of the 34 runs, and cannot be reached at 7be4133: LIBURN2 is called only from GETW and from BURSMIX, both only when LBURS is true (getw.f, lines 403–418 ↗; getw.f, lines 456–460 ↗), and GETW sets LBURS = .FALSE. on every call (getw.f, line 140 ↗), which no namelist setting changes. It would also need LOVSTE and a fully convective start.

5. (cannot be reached) No overflow guard in the last loop.

(lines 340–346 ↗) The last loop divides by EXP(DLI6), EXP(DLI7), EXP(DBE9) without the cut at 30 of the other loops or LIBURN's SAFEDIVEXP. For an exponent above 709.78 the default build gives 0, the intended result, and the development build (-ffpe-trap=overflow) stops the run.

A possible fix:

liburn2.f

@@ lines 338–348 @@
          DLI7 = DELTS*EXP(FRAD*LOG(RLI7(I))+(1.0D0-FRAD)*BLI70)
          DBE9 = DELTS*EXP(FRAD*LOG(RBE9(I))+(1.0D0-FRAD)*BBE90)
 C***REMEMBER TO ADD FAILSAFES FOR LARGE DEPLETION***
-         HCOMP(13,I) = HCOMP(13,I)/EXP(DLI6)
+C 11/09/26 KC via Claude
+C code changed to guard the division against overflow, as in LIBURN.
+         CALL SAFEDIVEXP(HCOMP(13,I),DLI6)
          IF(HCOMP(13,I).LT.1.0D-24)HCOMP(13,I)=0.0D0
-         HCOMP(14,I) = HCOMP(14,I)/EXP(DLI7)
+         CALL SAFEDIVEXP(HCOMP(14,I),DLI7)
          IF(HCOMP(14,I).LT.1.0D-24)HCOMP(14,I)=0.0D0
-         HCOMP(15,I) = HCOMP(15,I)/EXP(DBE9)
+         CALL SAFEDIVEXP(HCOMP(15,I),DBE9)
          IF(HCOMP(15,I).LT.1.0D-24)HCOMP(15,I)=0.0D0
       END DO
  200  CONTINUE

6. (cannot be reached) The bottom point's rates are left adjusted.

(lines 233–242 ↗; lines 270–271 ↗) As in LIBURN (liburn.f, suspected bug 4): Only HS2(JENV) is restored, and the adjusted RLI6, RLI7, RBE9 at JENV stay in NEWRAT_MOD. BURSMIX copies them into the start rates OLDRAT_MOD (bursmix.f, lines 138–142 ↗), where this routine's mean start rate (loop 68) weights the adjusted rate with the unadjusted mass.

A possible fix: as for LIBURN, save RLI6(JENV), RLI7(JENV), RBE9(JENV) with HS2SAV and restore them with it.

Definitions

subroutineLIBURN2

lines 8–350 · GitHub

One-step burning with the geometric-mean rate (above).

Uses modules CONSTANTS_MOD, DIMENSIONS_MOD, LIOV_MOD, MDPHY_MOD, NEWRAT_MOD, NMLPAR_MOD, OLDMOD_MOD, OLDRAT_MOD, PREVCZ_MOD, SCRTCH_MOD

Called by BURSMIX, GETW

Compiler warnings (1)

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

line 113: Equality comparison for REAL(8) [-Wcompare-reals]

Change log

The changes since rp_1

Changed since rp_1. Declarations only: COMMON blocks to modules (the namelist members of /DPMIX/ and /ROT/ to NMLPAR_MOD, /CONST1/ to CONSTANTS_MOD; the unused /LUOUT/ dropped), JSON from DIMENSIONS_MOD, IMPLICIT NONE.

Changed since a8a5950. The start rate is always computed, and the end rates are no longer stored in PREVCZ_MOD; the routine clears the new flag LBLIOLD instead. (Before, the test BLI6OLD.LE.0 was always true, so the start rate was always computed too.)

1C 2C 3C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 4C LIBURN 5C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 6C 11/91 HR ADDED TO CALL. 7C SUBROUTINE LIBURN2(DELTS,HCOMP,HD,HR,HS1,HS2,HT,JNENV,JNENV0,M) ! KC 2025-05-31 8 SUBROUTINE LIBURN2(DELTS,HCOMP,HR,HS1,HS2,HT,JNENV,JNENV0,M) 9C 10/15/26 MHP via Claude 10C code changed to take JSON from DIMENSIONS_MOD instead of declaring 11C PARAMETER(JSON=5000) here. 12 USE DIMENSIONS_MOD 13C 10/14/26 MHP via Claude 14C code changed to convert COMMON/SCRTCH/ into SCRTCH_MOD. 15 USE SCRTCH_MOD 16C 10/14/26 MHP via Claude 17C code changed to convert COMMON/PREVCZ/ into PREVCZ_MOD. 18 USE PREVCZ_MOD 19C 10/14/26 MHP via Claude 20C code changed to convert COMMON/OLDRAT/ into OLDRAT_MOD. 21 USE OLDRAT_MOD 22C 10/14/26 MHP via Claude 23C code changed to convert COMMON/OLDMOD/ into OLDMOD_MOD. 24 USE OLDMOD_MOD 25C 10/14/26 MHP via Claude 26C code changed to convert COMMON/NEWRAT/ into NEWRAT_MOD. 27 USE NEWRAT_MOD 28C 10/14/26 MHP via Claude 29C code changed to convert COMMON/MDPHY/ into MDPHY_MOD. 30 USE MDPHY_MOD 31C 10/14/26 MHP via Claude 32C code changed to convert COMMON/LIOV/ into LIOV_MOD. 33 USE LIOV_MOD 34C 10/14/26 MHP via Claude 35C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD. 36 USE NMLPAR_MOD 37C 10/14/26 MHP via Claude 38C code changed to convert COMMON/CONSTANTS/ into CONSTANTS_MOD. 39 USE CONSTANTS_MOD 40C 10/20/26 MHP via Claude 41C code changed to replace implicit typing with IMPLICIT NONE and 42C explicit declarations. 43 IMPLICIT NONE 44 REAL*8 BBE9,BBE90,BLI6,BLI60,BLI7,BLI70,CZBE9,CZLI6,CZLI7,DBE9, 45 * * DDJ,DDJM1,DELTS,DLI6,DLI7,DMASS,DR,FBE9,FBE90,FLI6,FLI60, 46 * * FLI7,FLI70,FM,FM0,FRAD,FX,HCOMP,HR,HS1,HS2,HS2SAV,HSBEG, 47 * * HSEND,HT,R,RATBE9,RATLI6,RATLI7,RBASE,ROS,TBE9,TLI6,TLI7 48 INTEGER I,J,JENV,JENV0,JNENV,JNENV0,M,MAXJ,MINJ 49 50C 10/13/26 MHP via Claude: removed the unused COMMON/IOUNITS/ declaration 51C (no member referenced here; residue of the obsolete-option removals). 52C 9/25/26 MHP via Claude 53C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/, 54C and COMMON/CONST3/ into a single COMMON/CONSTANTS/ 55C MHP 10/91 COMMON BLOCK ADDED FOR OVERSHOOT. 56C 10/12/26 MHP via Claude 57C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 58C namelist variables that the physics routines read from 27 COMMON 59C blocks which also held derived or working members. Those blocks keep 60C only their non-namelist members; five of them -- NEWMX, VMULT, 61C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 62C Derivation-only namelist variables from the same blocks went to 63C COMMON/NMLBYPASS/ instead, not here. 64C 65C Declared identically in all 121 files that hold any of these members. 66C Generated from one member list, so every copy is byte-identical; the 67C linker's 'size of symbol changed' warning is the check that they are. 68C 10/13/26 MHP via Claude: removed the unused COMMON/DPMIX/ declaration 69C (no member referenced here; residue of the obsolete-option removals). 70C MHP 11/91 COMMON BLOCK ADDED FOR OVERSHOOT. 71C 11/91 COMMON BLOCK ADDED FOR THE EXACT LOCATION OF THE BASE OF THE 72C SURFACE CZ AND THE INTEGRATED BURNING RATES THROUGHOUT IT. USED TO 73C REMEMBER THE PREVIOUS "END OF TIMESTEP VALUES" FOR THE NEW 74C "BEGINNING OF TIMESTEP VALUES". 75C DIMENSION HCOMP(15,JSON),HD(JSON),HS2(JSON),HT(JSON),HS1(JSON) ! KC 2025-05-31 76 DIMENSION HCOMP(15,JSON),HS2(JSON),HT(JSON),HS1(JSON) 77 DIMENSION RATLI6(JSON),RATLI7(JSON),RATBE9(JSON) 78 DIMENSION HR(JSON) 79C DATA TLIM/6.0D0/ 80 SAVE 81C DETERMINE LITHIUM AND BERYLLIUM BURNING. 82C THE BURNING RATES DEPEND ON THE LOCAL T AND RHO, THE ABUNDANCE, AND 83C IN THE CASE OF CONVECTION ZONES THE MIXING. 84C 85C IN RADIATIVE REGIONS BURNING IS DONE IMPLICITLY. 86C IN CONVECTIVE REGIONS THE PROGRAM ITERATES BETWEEN BURNING AND MIXING 87C TO CORRECT THE RATES FOR THE CHANGE IN ABUNDANCE AS A FUNCTION OF TIME. 88C THE DEGREE OF LITHIUM BURNING IN A SURFACE CZ DEPENDS SENSITIVELY 89C ON THE TEMPERATURE AT ITS BASE - SO ACCURATELY LOCATING IS IMPORTANT. 90C DETERMINE THE TRUE LOCATION (FX) OF THE BASE OF THE CZ AT THE END 91C OF THE TIMESTEP, AND THE LOCATION OF THE EDGE OF OVERSHOOT REGIONS 92C IF APPLICABLE. 93 IF(JNENV.GT.1.AND.JNENV.LT.M)THEN 94 IF(LROT.AND.LINSTB)THEN 95 DDJ = DELAM(JNENV) - DELRM(JNENV) 96 DDJM1 = DELAM(JNENV-1) - DELRM(JNENV-1) 97 ELSE 98C EVALUATE DEL(AD) - DEL(RAD) AT THE LAST CONVECTIVE POINT AND THE ONE 99C BELOW IT. 100 DDJ = SDEL(3,JNENV)-SDEL(1,JNENV) 101 DDJM1 = SDEL(3,JNENV-1)-SDEL(1,JNENV-1) 102 ENDIF 103C USE LINEAR INTERPOLATION TO FIND THE DISTANCE OF THE TRUE LOCATION 104C OF THE BASE FROM THE ZONE MIDPOINT. IF FX IS NEGATIVE,THEN THE TRUE 105C BASE IS HIGHER; IF IT IS POSITIVE, THE TRUE BASE IS LOWER. 106 FX = MAX(-0.5D0,0.5D0-DDJM1/(DDJM1-DDJ)) 107 FX = MIN(0.5D0,FX) 108 IF(.NOT.LOVSTE)THEN 109 JENV = JNENV 110 JENV0 = JNENV0 111 ELSE 112C STARTING CZ DEPTH 113 IF(RBASE0.EQ.0.0D0)THEN 114 RBASE0 = 0.5D0*(EXP(CLN*HRO(JNENV0)) 115 * * +EXP(CLN*HRO(JNENV0-1))) 116 ROS = RBASE0 - PSCAL0 117 DO J = JNENV0-1,1,-1 118 R = EXP(CLN*HRO(J)) 119 IF(R.LT.ROS)THEN 120 JENV0 = J + 1 121 GOTO 11 122 ENDIF 123 END DO 124 JENV0 = 1 125 11 11 CONTINUE 126 ELSE 127 JENV0 = JENVOS 128 ENDIF 129C ENDING CZ DEPTH : DETERMINE OVERSHOOT FROM TRUE CZ BASE. 130 DR = EXP(CLN*HR(JNENV))-EXP(CLN*HR(JNENV-1)) 131 RBASE = 0.5D0*(EXP(CLN*HR(JNENV))+EXP(CLN*HR(JNENV-1))) 132 * * -FX*DR 133 RBASE0 = RBASE 134 ROS = RBASE - PSCAL1 135 DO J = JNENV-1,1,-1 136 R = EXP(CLN*HR(J)) 137 IF(R.LT.ROS)THEN 138 JENV = J + 1 139 DR = EXP(CLN*HRO(J+1))-R 140 FX = 0.5D0-((ROS-R)/DR) 141 FX = MAX(-0.5D0,FX) 142 FX = MIN(0.5D0,FX) 143 GOTO 12 144 ENDIF 145 END DO 146 JENV = 1 147 12 12 CONTINUE 148 ENDIF 149 ELSE 150 JENV = JNENV 151 JENV0 = JNENV0 152 ENDIF 153 JENVOS = JENV 154C RADIATIVE INTERIOR. 155 MINJ = MIN(JENV,JENV0) 156 MAXJ = MAX(JENV,JENV0) 157 DO I = 1,MINJ-1 158 IF(RBE9(I).LE.1.0D-32 .OR. RBE90(I).LE.1.0D-32)GOTO 50 159 IF(HCOMP(13,I).LT.1.0D-24.AND.HCOMP(14,I).LT.1.0D-24 160 * * .AND.HCOMP(15,I).LT.1.0D-24)GOTO 50 161 IF(HT(I).GT.7.0D0)THEN 162 HCOMP(13,I) = 0.0D0 163 HCOMP(14,I) = 0.0D0 164 HCOMP(15,I) = 0.0D0 165 GOTO 50 166 ENDIF 167 TLI6 = 0.5D0*(LOG(RLI6(I)) + LOG(RLI60(I))) 168 TLI7 = 0.5D0*(LOG(RLI7(I)) + LOG(RLI70(I))) 169 TBE9 = 0.5D0*(LOG(RBE9(I)) + LOG(RBE90(I))) 170 RATLI6(I) = DELTS*EXP(TLI6) 171 RATLI7(I) = DELTS*EXP(TLI7) 172 RATBE9(I) = DELTS*EXP(TBE9) 173C THE REACTION RATES ARE OF THE FORM 174C DLI/DT = F(RHO,T,X)*LI =>DLN LI = DT*F(RHO,T,X) 175C SOLVE FOR D LN (SPECIES)/DT AND ZERO OUT IF THE DEPLETION IS TOO HIGH. 176 IF(RATLI6(I).LT.3.0D1)THEN 177 HCOMP(13,I) = HCOMP(13,I)/EXP(RATLI6(I)) 178 ELSE 179 HCOMP(13,I) = 0.0D0 180 ENDIF 181 IF(RATLI7(I).LT.3.0D1)THEN 182 HCOMP(14,I) = HCOMP(14,I)/EXP(RATLI7(I)) 183 ELSE 184 HCOMP(14,I) = 0.0D0 185 ENDIF 186 IF(RATBE9(I).LT.3.0D1)THEN 187 HCOMP(15,I) = HCOMP(15,I)/EXP(RATBE9(I)) 188 ELSE 189 HCOMP(15,I) = 0.0D0 190 ENDIF 191C WRITE NEW ABUNDANCES AND EXIT. 192 IF(HCOMP(13,I).LT.1.0D-24)HCOMP(13,I)=0.0D0 193 IF(HCOMP(14,I).LT.1.0D-24)HCOMP(14,I)=0.0D0 194 IF(HCOMP(15,I).LT.1.0D-24)HCOMP(15,I)=0.0D0 195 50 50 CONTINUE 196 END DO 197C CONVECTION ZONE. 198C 199C SKIP IF WHOLE CZ IS BELOW THE BURNING THRESHOLD. 200 IF(RBE90(JENV0).LE.1.0D-32.OR.RBE9(JENV).LE.1.0D-32)GOTO 200 201C FIND RATES AT THE BEGINNING OF THE TIMESTEP (USING THE DEPTH AT THE START). 202 FLI60 = 0.0D0 203 FLI70 = 0.0D0 204 FBE90 = 0.0D0 205 FM0 = 0.0D0 206 DO 65 I = JENV0,M 207 FLI60 = FLI60+HCOMP(13,I)*HS2(I) 208 FLI70 = FLI70+HCOMP(14,I)*HS2(I) 209 FBE90 = FBE90+HCOMP(15,I)*HS2(I) 210 FM0 = FM0 + HS2(I) 211 65 65 CONTINUE 212C 67 CONTINUE 213 FLI60 = FLI60/FM0 214 FLI70 = FLI70/FM0 215 FBE90 = FBE90/FM0 216C 11/02/26 MHP via Claude 217C code changed to always compute the start-of-step rates here. BURSMIX 218C restarts each extrapolation sequence from the start-of-step model and 219C does not advance RLI60, so an end-of-step rate from an earlier call 220C must not be reused (the old BLI6OLD.LE.0 test never allowed it). 221C MASS-WEIGHTED AVERAGE RATES AT THE START OF THE STEP. 222 BLI60 = 0.0D0 223 BLI70 = 0.0D0 224 BBE90 = 0.0D0 225 DO 68 I = JENV0,M 226 BLI60 = BLI60 + RLI60(I)*HS2(I) 227 BLI70 = BLI70 + RLI70(I)*HS2(I) 228 BBE90 = BBE90 + RBE90(I)*HS2(I) 229 68 68 CONTINUE 230 BLI60 = LOG(BLI60/FM0) 231 BLI70 = LOG(BLI70/FM0) 232 BBE90 = LOG(BBE90/FM0) 233C USE THE LOCATION OF THE TRUE EDGE OF THE CONVECTION ZONE (FX, FOUND AT 234C BEGINNING OF SR) TO ADJUST THE BURNING RATE AND MASS OF THE BOTTOM POINT 235C SUCH THAT IT INCLUDES THE ENTIRE C.Z. 236 HS2SAV = HS2(JENV) 237 IF(JENV.GT.1.AND.JENV.LT.M)THEN 238 HS2(JENV) = HS2(JENV)+FX*(HS1(JENV)-HS1(JENV-1)) 239 RLI6(JENV) = RLI6(JENV)+0.5D0*FX*(RLI6(JENV-1)-RLI6(JENV)) 240 RLI7(JENV) = RLI7(JENV)+0.5D0*FX*(RLI7(JENV-1)-RLI7(JENV)) 241 RBE9(JENV) = RBE9(JENV)+0.5D0*FX*(RBE9(JENV-1)-RBE9(JENV)) 242 ENDIF 243C FIND RATES AT THE END OF THE TIMESTEP (USING THE DEPTH AT THE END). 244C ALSO STORE INITIAL ABUNDANCES(FLI60,FLI70,FBE90). 245C FM IS THE TOTAL MASS IN THE CZ. 246 FLI6 = 0.0D0 247 FLI7 = 0.0D0 248 FBE9 = 0.0D0 249 FM = 0.0D0 250 BLI6 = 0.0D0 251 BLI7 = 0.0D0 252 BBE9 = 0.0D0 253 DO 70 I = JENV,M 254 BLI6 = BLI6 + RLI6(I)*HS2(I) 255 BLI7 = BLI7 + RLI7(I)*HS2(I) 256 BBE9 = BBE9 + RBE9(I)*HS2(I) 257 FLI6 = FLI6+HCOMP(13,I)*HS2(I) 258 FLI7 = FLI7+HCOMP(14,I)*HS2(I) 259 FBE9 = FBE9+HCOMP(15,I)*HS2(I) 260 FM = FM + HS2(I) 261 70 70 CONTINUE 262C 75 CONTINUE 263 FLI6 = FLI6/FM 264 FLI7 = FLI7/FM 265 FBE9 = FBE9/FM 266C MASS-WEIGHTED AVERAGE RATES AT THE END OF THE STEP. 267 BLI6 = LOG(BLI6/FM) 268 BLI7 = LOG(BLI7/FM) 269 BBE9 = LOG(BBE9/FM) 270C RESTORE ORIGINAL MASS TO BASE OF C.Z. 271 HS2(JENV) = HS2SAV 272C FOR THE STARTING ABUNDANCE, USE THE MIXED ABUNDANCE FOR THE 273C DEEPER CZ. 274 IF(JENV.LT.JENV0)THEN 275 FLI60 = FLI6 276 FLI70 = FLI7 277 FBE90 = FBE9 278 ENDIF 279C INITIALIZE ABUNDANCES. 280 FLI6 = FLI60 281 FLI7 = FLI70 282 FBE9 = FBE90 283C TIME-AVERAGED RATES USING LINEAR INTERPOLATION IN THE LOG. 284 CZLI6 = 0.5D0*(BLI6+BLI60) 285 CZLI7 = 0.5D0*(BLI7+BLI70) 286 CZBE9 = 0.5D0*(BBE9+BBE90) 287 DLI6 = DELTS*EXP(CZLI6) 288 DLI7 = DELTS*EXP(CZLI7) 289 DBE9 = DELTS*EXP(CZBE9) 290 IF(DLI6.LT.3.0D1 .AND. FLI6.GT.1.0D-24)THEN 291 FLI6 = FLI6/EXP(DLI6) 292 ELSE 293 FLI6 = 0.0D0 294 ENDIF 295 IF(FLI6.LT.1.0D-24)FLI6 = 0.0D0 296 IF(DLI7.LT.3.0D1 .AND. FLI7.GT.1.0D-24)THEN 297 FLI7 = FLI7/EXP(DLI7) 298 ELSE 299 FLI7 = 0.0D0 300 ENDIF 301 IF(FLI7.LT.1.0D-24)FLI7 = 0.0D0 302 IF(DBE9.LT.3.0D1 .AND. FBE9.GT.1.0D-24)THEN 303 FBE9 = FBE9/EXP(DBE9) 304 ELSE 305 FBE9 = 0.0D0 306 ENDIF 307 IF(FBE9.LT.1.0D-24)FBE9 = 0.0D0 308 DO I = MAXJ,M 309 HCOMP(13,I) = FLI6 310 HCOMP(14,I) = FLI7 311 HCOMP(15,I) = FBE9 312 END DO 313C 11/02/26 MHP via Claude 314C code changed to clear LBLIOLD instead of storing the ending rates, so 315C LIBURN does not reuse a rate from this routine (see above). 316 LBLIOLD = .FALSE. 317C NOW SOLVE FOR ABUNDANCES IN THE REGION WHICH BEGAN CONVECTIVE AND 318C ENDED RADIATIVE. 319 IF(JENV.LE.JENV0)GOTO 200 320C FIND STARTING AND ENDING LOCATION IN MASS OF THE CZ BASE, AND 321C ASSUME THAT THE FRACTION OF TIME SPENT IN THE CZ IS PROPORTIONAL TO 322C THE LOCATION IN MASS (I.E. A POINT 1/3 FROM THE OLD TO THE NEW CZ 323C BASE IN MASS SPENT 1/3 OF THE TIME IN THE CZ AND 2/3 OUT OF IT). 324 IF(JENV0.GT.1)THEN 325 HSBEG = 0.5D0*(HS1(JENV0)+HS1(JENV0-1)) 326 ELSE 327 HSBEG = 0.0D0 328 ENDIF 329 HSEND = 0.5D0*(HS1(JENV)+HS1(JENV-1)) 330 DMASS = HSBEG - HSEND 331 DO I = JENV0,JENV-1 332C MHP 9/91 CHANGE TO AVOID DIVISION BY ZERO. 333C SKIP IF SHELL TEMPERATURE DROPS BELOW BURNING THRESHOLD. 334 IF(RBE9(I).LE.1.0D-32)GOTO 200 335 FRAD = (HS1(I)-HSBEG)/DMASS 336C USE FRAD*RADIATIVE RATE AND (1-FRAD)*CONVECTIVE RATE. 337 DLI6 = DELTS*EXP(FRAD*LOG(RLI6(I))+(1.0D0-FRAD)*BLI60) 338 DLI7 = DELTS*EXP(FRAD*LOG(RLI7(I))+(1.0D0-FRAD)*BLI70) 339 DBE9 = DELTS*EXP(FRAD*LOG(RBE9(I))+(1.0D0-FRAD)*BBE90) 340C***REMEMBER TO ADD FAILSAFES FOR LARGE DEPLETION*** 341 HCOMP(13,I) = HCOMP(13,I)/EXP(DLI6) 342 IF(HCOMP(13,I).LT.1.0D-24)HCOMP(13,I)=0.0D0 343 HCOMP(14,I) = HCOMP(14,I)/EXP(DLI7) 344 IF(HCOMP(14,I).LT.1.0D-24)HCOMP(14,I)=0.0D0 345 HCOMP(15,I) = HCOMP(15,I)/EXP(DBE9) 346 IF(HCOMP(15,I).LT.1.0D-24)HCOMP(15,I)=0.0D0 347 END DO 348 200 200 CONTINUE 349 RETURN 350 END

1"""Python version of liburn2.f (YREC development, 7be4133).""" 2 3import math 4 5import numpy as np 6 7from pyrec import constants_mod 8from pyrec import dimensions_mod 9from pyrec import liov_mod 10from pyrec import mdphy_mod 11from pyrec import newrat_mod 12from pyrec import nmlpar_mod 13from pyrec import oldmod_mod 14from pyrec import oldrat_mod 15from pyrec import prevcz_mod 16from pyrec import scrtch_mod 17 18 19def liburn2(delts: float, hcomp: np.ndarray, hr: np.ndarray, hs1: np.ndarray, hs2: np.ndarray, ht: np.ndarray, 20 jnenv: int, jnenv0: int, m: int) -> None: 21 """Burn Li6, Li7 and Be9 over one step in one go (cheap LIBURN). 22 23 The same as LIBURN without substeps, extrapolation or accretion: each 24 radiative point, and the convection zone as a whole, decays with the 25 geometric mean of the start and end rates. Points left behind by a 26 retreating zone get LIBURN's blended rate. 27 28 Translates: liburn2.f#LIBURN2 (7be4133) 29 30 Parameters 31 ---------- 32 delts : float 33 The time step. 34 hcomp : ndarray, shape (16, JSON + 1) 35 Composition; rows 13, 14, 15 (Li6, Li7, Be9) updated in place. 36 hr : ndarray 37 log r at each point. 38 hs1, hs2 : ndarray 39 Mass inside each point, and each point's mass (g). ``hs2[jenv]`` 40 is changed and restored. 41 ht : ndarray 42 log T at each point. 43 jnenv, jnenv0 : int 44 The base of the surface convection zone at the end and the start 45 of the step. 46 m : int 47 Number of points. 48 49 Also changes ``newrat_mod`` (the bottom point's rates) and 50 ``prevcz_mod`` (``rbase0``, ``lbliold``, ``jenvos``). Nothing is 51 returned. 52 53 """ 54 55 # Module data (Fortran lines 12–39). 56 json = dimensions_mod.json 57 cln = constants_mod.cln 58 rli6 = newrat_mod.rli6 59 rli7 = newrat_mod.rli7 60 rbe9 = newrat_mod.rbe9 61 rli60 = oldrat_mod.rli60 62 rli70 = oldrat_mod.rli70 63 rbe90 = oldrat_mod.rbe90 64 hro = oldmod_mod.hro 65 sdel = scrtch_mod.sdel 66 67 # Local arrays (Fortran lines 76–80). SAVE in the Fortran: no local is 68 # read before it is set in a call. 69 ratli6 = np.zeros(json + 1) 70 ratli7 = np.zeros(json + 1) 71 ratbe9 = np.zeros(json + 1) 72 fx = 0.0 73 74 # The base of the convection zone, as in LIBURN (Fortran lines 88–153). 75 if 1 < jnenv < m: 76 if nmlpar_mod.lrot and nmlpar_mod.linstb: 77 ddj = mdphy_mod.delam[jnenv] - mdphy_mod.delrm[jnenv] 78 ddjm1 = mdphy_mod.delam[jnenv - 1] - mdphy_mod.delrm[jnenv - 1] 79 else: 80 ddj = sdel[3, jnenv] - sdel[1, jnenv] 81 ddjm1 = sdel[3, jnenv - 1] - sdel[1, jnenv - 1] 82 fx = max(-0.5, 0.5 - ddjm1 / (ddjm1 - ddj)) 83 fx = min(0.5, fx) 84 if not nmlpar_mod.lovste: 85 jenv = jnenv 86 jenv0 = jnenv0 87 else: 88 # Envelope overshoot: the start of the step (Fortran lines 89 # 112–128). 90 if prevcz_mod.rbase0 == 0.0: 91 # Suspected bug (see the page): reads hro[0] when jnenv0 = 1 92 # (0 here, outside the array in the Fortran); LIBURN guards 93 # this. 94 prevcz_mod.rbase0 = 0.5 * (math.exp(cln * hro[jnenv0]) 95 + math.exp(cln * hro[jnenv0 - 1])) 96 ros = prevcz_mod.rbase0 - liov_mod.pscal0 97 jenv0 = 1 98 for j in range(jnenv0 - 1, 0, -1): 99 r = math.exp(cln * hro[j]) 100 if r < ros: 101 jenv0 = j + 1 102 break # to label 11 103 else: 104 jenv0 = prevcz_mod.jenvos 105 106 # The end of the step, from the true base (Fortran lines 107 # 129–147). 108 dr = math.exp(cln * hr[jnenv]) - math.exp(cln * hr[jnenv - 1]) 109 rbase = (0.5 * (math.exp(cln * hr[jnenv]) 110 + math.exp(cln * hr[jnenv - 1])) - fx * dr) 111 prevcz_mod.rbase0 = rbase 112 ros = rbase - liov_mod.pscal1 113 jenv = 1 114 for j in range(jnenv - 1, 0, -1): 115 r = math.exp(cln * hr[j]) 116 if r < ros: 117 jenv = j + 1 118 # Suspected bug (see the page): HRO, the old radius, here. 119 dr = math.exp(cln * hro[j + 1]) - r 120 fx = 0.5 - ((ros - r) / dr) 121 fx = max(-0.5, fx) 122 fx = min(0.5, fx) 123 break # to label 12 124 else: 125 jenv = jnenv 126 jenv0 = jnenv0 127 prevcz_mod.jenvos = jenv 128 129 # The radiative interior, below both bases: skip, not leave, points 130 # with a negligible Be9 rate (label 50 is inside the loop) (Fortran 131 # lines 154–196). 132 minj = min(jenv, jenv0) 133 maxj = max(jenv, jenv0) 134 for i in range(1, minj): 135 if rbe9[i] <= 1.0e-32 or rbe90[i] <= 1.0e-32: 136 continue 137 if (hcomp[13, i] < 1.0e-24 and hcomp[14, i] < 1.0e-24 138 and hcomp[15, i] < 1.0e-24): 139 continue 140 if ht[i] > 7.0: 141 hcomp[13, i] = 0.0 142 hcomp[14, i] = 0.0 143 hcomp[15, i] = 0.0 144 continue 145 146 # The geometric mean of the start and end rates, and the burning 147 # (Fortran lines 167–194). 148 tli6 = 0.5 * (math.log(rli6[i]) + math.log(rli60[i])) 149 tli7 = 0.5 * (math.log(rli7[i]) + math.log(rli70[i])) 150 tbe9 = 0.5 * (math.log(rbe9[i]) + math.log(rbe90[i])) 151 ratli6[i] = delts * math.exp(tli6) 152 ratli7[i] = delts * math.exp(tli7) 153 ratbe9[i] = delts * math.exp(tbe9) 154 if ratli6[i] < 3.0e1: 155 hcomp[13, i] = hcomp[13, i] / math.exp(ratli6[i]) 156 else: 157 hcomp[13, i] = 0.0 158 if ratli7[i] < 3.0e1: 159 hcomp[14, i] = hcomp[14, i] / math.exp(ratli7[i]) 160 else: 161 hcomp[14, i] = 0.0 162 if ratbe9[i] < 3.0e1: 163 hcomp[15, i] = hcomp[15, i] / math.exp(ratbe9[i]) 164 else: 165 hcomp[15, i] = 0.0 166 if hcomp[13, i] < 1.0e-24: 167 hcomp[13, i] = 0.0 168 if hcomp[14, i] < 1.0e-24: 169 hcomp[14, i] = 0.0 170 if hcomp[15, i] < 1.0e-24: 171 hcomp[15, i] = 0.0 172 173 # The convection zone; skipped if it does not burn (Fortran lines 174 # 197–200). 175 if rbe90[jenv0] <= 1.0e-32 or rbe9[jenv] <= 1.0e-32: 176 return # label 200 177 178 # Mean abundances over the zone at the start of the step (loop 65, 179 # Fortran lines 201–215). 180 fli60 = 0.0 181 fli70 = 0.0 182 fbe90 = 0.0 183 fm0 = 0.0 184 for i in range(jenv0, m + 1): 185 fli60 = fli60 + hcomp[13, i] * hs2[i] 186 fli70 = fli70 + hcomp[14, i] * hs2[i] 187 fbe90 = fbe90 + hcomp[15, i] * hs2[i] 188 fm0 = fm0 + hs2[i] 189 fli60 = fli60 / fm0 190 fli70 = fli70 / fm0 191 fbe90 = fbe90 / fm0 192 193 # Always the zone's mean start rates: BURSMIX restarts from the 194 # start-of-step model, so no end rate from an earlier call is reused 195 # (Fortran lines 216–232). 196 bli60 = 0.0 197 bli70 = 0.0 198 bbe90 = 0.0 199 for i in range(jenv0, m + 1): # loop 68 200 bli60 = bli60 + rli60[i] * hs2[i] 201 bli70 = bli70 + rli70[i] * hs2[i] 202 bbe90 = bbe90 + rbe90[i] * hs2[i] 203 bli60 = math.log(bli60 / fm0) 204 bli70 = math.log(bli70 / fm0) 205 bbe90 = math.log(bbe90 / fm0) 206 207 # The bottom point's mass and rates adjusted for the offset fx 208 # (Fortran lines 233–242). 209 # Suspected bug (see the page): the rates are not restored afterwards. 210 hs2sav = hs2[jenv] 211 if 1 < jenv < m: 212 hs2[jenv] = hs2[jenv] + fx * (hs1[jenv] - hs1[jenv - 1]) 213 rli6[jenv] = rli6[jenv] + 0.5 * fx * (rli6[jenv - 1] - rli6[jenv]) 214 rli7[jenv] = rli7[jenv] + 0.5 * fx * (rli7[jenv - 1] - rli7[jenv]) 215 rbe9[jenv] = rbe9[jenv] + 0.5 * fx * (rbe9[jenv - 1] - rbe9[jenv]) 216 217 # Mean abundances and rates over the zone at the end, then the mass 218 # restored (loop 70, Fortran lines 243–271). 219 fli6 = 0.0 220 fli7 = 0.0 221 fbe9 = 0.0 222 fm = 0.0 223 bli6 = 0.0 224 bli7 = 0.0 225 bbe9 = 0.0 226 for i in range(jenv, m + 1): 227 bli6 = bli6 + rli6[i] * hs2[i] 228 bli7 = bli7 + rli7[i] * hs2[i] 229 bbe9 = bbe9 + rbe9[i] * hs2[i] 230 fli6 = fli6 + hcomp[13, i] * hs2[i] 231 fli7 = fli7 + hcomp[14, i] * hs2[i] 232 fbe9 = fbe9 + hcomp[15, i] * hs2[i] 233 fm = fm + hs2[i] 234 fli6 = fli6 / fm 235 fli7 = fli7 / fm 236 fbe9 = fbe9 / fm 237 bli6 = math.log(bli6 / fm) 238 bli7 = math.log(bli7 / fm) 239 bbe9 = math.log(bbe9 / fm) 240 hs2[jenv] = hs2sav 241 242 # A deepened zone starts from the end-of-step mean (Fortran lines 243 # 272–278). 244 if jenv < jenv0: 245 fli60 = fli6 246 fli70 = fli7 247 fbe90 = fbe9 248 249 # Burn the zone once with the geometric mean of its start and end 250 # rates (Fortran lines 279–307). 251 fli6 = fli60 252 fli7 = fli70 253 fbe9 = fbe90 254 czli6 = 0.5 * (bli6 + bli60) 255 czli7 = 0.5 * (bli7 + bli70) 256 czbe9 = 0.5 * (bbe9 + bbe90) 257 dli6 = delts * math.exp(czli6) 258 dli7 = delts * math.exp(czli7) 259 dbe9 = delts * math.exp(czbe9) 260 if dli6 < 3.0e1 and fli6 > 1.0e-24: 261 fli6 = fli6 / math.exp(dli6) 262 else: 263 fli6 = 0.0 264 if fli6 < 1.0e-24: 265 fli6 = 0.0 266 if dli7 < 3.0e1 and fli7 > 1.0e-24: 267 fli7 = fli7 / math.exp(dli7) 268 else: 269 fli7 = 0.0 270 if fli7 < 1.0e-24: 271 fli7 = 0.0 272 if dbe9 < 3.0e1 and fbe9 > 1.0e-24: 273 fbe9 = fbe9 / math.exp(dbe9) 274 else: 275 fbe9 = 0.0 276 if fbe9 < 1.0e-24: 277 fbe9 = 0.0 278 279 # The zone's new abundances, and LBLIOLD cleared (Fortran lines 280 # 308–316). 281 # Suspected bug (see the page): from maxj, so points newly in a 282 # deepened zone keep their old abundances. 283 for i in range(maxj, m + 1): 284 hcomp[13, i] = fli6 285 hcomp[14, i] = fli7 286 hcomp[15, i] = fbe9 287 prevcz_mod.lbliold = False 288 289 # The region left by a retreating zone: burned once, blended rate 290 # (Fortran lines 317–347). 291 if jenv <= jenv0: 292 return # label 200 293 if jenv0 > 1: 294 hsbeg = 0.5 * (hs1[jenv0] + hs1[jenv0 - 1]) 295 else: 296 hsbeg = 0.0 297 hsend = 0.5 * (hs1[jenv] + hs1[jenv - 1]) 298 # Suspected bug (see the page): dmass < 0, so frad lies in (-1, 0), 299 # and the weights are the other way round from the comment. 300 dmass = hsbeg - hsend 301 for i in range(jenv0, jenv): 302 # Question (see the page): ends the routine, not the point. 303 if rbe9[i] <= 1.0e-32: 304 return # label 200 305 frad = (hs1[i] - hsbeg) / dmass 306 dli6 = delts * math.exp(frad * math.log(rli6[i]) + (1.0 - frad) * bli60) 307 dli7 = delts * math.exp(frad * math.log(rli7[i]) + (1.0 - frad) * bli70) 308 dbe9 = delts * math.exp(frad * math.log(rbe9[i]) + (1.0 - frad) * bbe90) 309 # Suspected bug (see the page): no overflow guard here (math.exp 310 # raises OverflowError where the Fortran's EXP overflows). 311 hcomp[13, i] = hcomp[13, i] / math.exp(dli6) 312 if hcomp[13, i] < 1.0e-24: 313 hcomp[13, i] = 0.0 314 hcomp[14, i] = hcomp[14, i] / math.exp(dli7) 315 if hcomp[14, i] < 1.0e-24: 316 hcomp[14, i] = 0.0 317 hcomp[15, i] = hcomp[15, i] / math.exp(dbe9) 318 if hcomp[15, i] < 1.0e-24: 319 hcomp[15, i] = 0.0

liburn2.f

Purpose

How it fits

The method

Suspected bugs, checked by running the code

1. (cannot be reached) The sign of FRAD.

2. (cannot be reached) Points newly in a deepened zone are not burned.

3. (cannot be reached) The overshoot zone's width mixes two models.

4. (cannot be reached) HRO(0) read for a fully convective start.

5. (cannot be reached) No overflow guard in the last loop.

6. (cannot be reached) The bottom point's rates are left adjusted.

Definitions

subroutineLIBURN2

Change log

1. (cannot be reached) The sign of `FRAD`.

4. (cannot be reached) `HRO(0)` read for a fully convective start.