surfp.f

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

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

Variables on this page (14)

ap, at, atmgl, atmpl, atmtl, imin, iminmax, iowr, ishort, istor, jp, kp, ng, nt

Purpose

SURFP gives the pressure at $T = T_\mathrm{eff}$ in a model atmosphere, by interpolating a table of $\log P$ in $\log T_\mathrm{eff}$ and $\log g$. It is the outer boundary of the envelope integration with Kurucz atmospheres (KTTAU = 3). With KTTAU = 4 SETUPS reads the same table, but SURFP is not called: a run with Allard atmospheres that leaves their table falls back to the gray atmosphere, not to this one (atmos2_mod.f, question 1).

How it fits

Inputs: $\log T_\mathrm{eff}$ (TEFFL) and $\log g$ (GL); the table in ATMOS2_MOD (ATMPL(57,11), the grids ATMTL, ATMGL), read by SETUPS; and in FAC_MOD, the lowest and highest valid gravity index for each temperature (IMIN, IMINMAX), also set by SETUPS.
Who calls it: ENVINT, for KTTAU = 3. Its siblings are KCSURFP (Kurucz–Castelli tables, KTTAU = 5) and ALSURFP (Allard tables, KTTAU = 4).
What it calls: INTER3 (3-point Lagrange weights), KSPLINE and KSPLINT (a natural cubic spline through four points).
Outputs: in ATMPRT_MOD, AP = $\log P$ and AT = $\log T_\mathrm{eff}$, read by ENVINT (and QATM, MASSLOSS); JP and KP in FAC_MOD, the table cell used. With LPRT, a line on ISHORT and ISTOR (units from IOUNITS_MOD).

The method

Stop the run if $\log T_\mathrm{eff} < 3.5$ or $\log g < -0.5$, with exit status 1.

Question 1 (line 44 ↗; lines 70–71 ↗). Hard-coded limits. 3.5, $-0.5$, 4.5 and 4.55 are written in the code and tie it to one table; a different table would need them changed (the Kurucz 1990 tables start at $\log T_\mathrm{eff}$ = 3.544, so between 3.5 and 3.544 the spline extrapolates). Should they come from the table's grid (ATMTL, ATMGL)?
Find four table temperatures around $T_\mathrm{eff}$.
At each of them, interpolate $\log P$ in $\log g$: for $\log T > 4.55$ linearly between the two highest gravities; for $4.5 < \log T \le 4.55$ with the 3-point Lagrange formula on the three highest gravities; otherwise with a natural cubic spline through the four valid gravities nearest $\log g$ (the top four, a bracketing four, or the bottom four).

Question 2 (lines 70–87 ↗). Hot temperatures. For $\log T > 4.5$ only the two or three highest gravities are used, whatever $\log g$ is: probably because the table has no low-gravity values there (the Kurucz 1990 tables have $-999$ below $\log g$ = 4.0 at $\log T \ge 4.50$). Is linear extrapolation to low gravity acceptable? SURFBC switches a Kurucz run to the gray atmosphere above $\log T_\mathrm{eff} = 3.95$ (surfbc.f, lines 83–95 ↗), so the four table temperatures stay below about 4.0 there; these branches would be reached only by a call of ENVINT with KTTAU = 3 for a hotter model before that switch (not traced).
Interpolate the four values in $\log T_\mathrm{eff}$ with a natural cubic spline through four points.

The four-point spline, for points $x_1 < \dots < x_4$, is the cubic-spline formula of splint.f with natural ends.

The header (lines 38–41 ↗) describes an older method: It says the routine interpolates "USING A 4-POINT LAGRANGIAN INTERPOLATOR" and "WILL QUIT IF THE DESIRED DATA POINT HAS TEFF OR LOG G MORE THAN ONE TABLE POINT FROM THE DATA," but the code uses 4-point splines (the Lagrange calls INTRP2 are commented out), except for table rows with $\log T_\mathrm{eff} > 4.5$, where gravity is interpolated linearly (above 4.55, lines 71–75 ↗) or with the 3-point Lagrange INTER3 (lines 77–82 ↗), and it checks only the lower limits of $T_\mathrm{eff}$ and $g$ (line 44 ↗).

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. (never occurred) Missing gravities used at the top of a row.

(lines 88–102 ↗; lines 104–123 ↗) SETUPS sets IMINMAX(J), the highest gravity with a pressure, below NG when a temperature's highest-gravity entries are $-999$ (setups.f, lines 270–299 ↗). G. Somers limited the general case's search to IMINMAX(J) "to prevent using -999," but the top-four branch (GL.GE.ATMGL(NG-1)) still takes the gravities NG-3 … NG, and KK = MIN(NG-3,KK) still allows them. The general search covers $\log g$ only up to ATMGL(IMINMAX(J)-1). In such a row, a $\log g$ from there up to ATMGL(NG-1) = 4.5 matches no interval and falls through to the lowest four valid gravities, a long extrapolation; a $\log g \ge 4.5$ takes the top four, $-999$ entries included, into the spline. The tables in input/atmos/kurucz/ with $[\mathrm{Fe/H}] \le -2.3$ (atmk1990m23.tab to atmk1990m50.tab) have such rows at $\log T_\mathrm{eff}$ = 3.53–3.68 (for example atmk1990m40.tab: no pressure above $\log g$ = 3.0 at 3.544), where cool dwarfs have $\log g \approx 4.5$–5; the other Kurucz 1990 tables have none below $\log T$ = 4.5, apart from an empty row at 4.041 in atmk1990p01.tab, above the switch to gray at 3.95. What follows: a wrong surface pressure AP for a cool metal-poor dwarf. Only with KTTAU = 3 (the default is 0, gray) and one of these tables.

A possible fix:

surfp.f

@@ lines 85–95 @@
             ENDIF
             GOTO 20
          ENDIF
-         IF(GL.GE.ATMGL(NG-1))THEN
+C 11/09/26 KC via Claude
+C code changed to use IMINMAX(J), the highest gravity with a pressure,
+C in place of NG, so that missing (-999) values at the top of a row
+C are not used, and the general search reaches the top of the row.
+         IF(GL.GE.ATMGL(IMINMAX(J)-1))THEN
 C DESIRED LOG G ABOVE SECOND TO TOP TABLE LOG G - USE TOP 4 LOG G VALUES.
             DO KK = 1,4
-               QG(KK)=ATMGL(NG-4+KK)
-               PTAB(KK) = ATMPL(J,NG-4+KK)
+               QG(KK)=ATMGL(IMINMAX(J)-4+KK)
+               PTAB(KK) = ATMPL(J,IMINMAX(J)-4+KK)
             END DO
 C            CALL INTRP2(QR,QG,GL)
             CALL KSPLINE(QG,PTAB,YG)
@@ lines 97–103 @@
 C            PP(N) = ATMPL(J,NG-3)*QG(1)+ATMPL(J,NG-2)*QG(2)+
 C     *              ATMPL(J,NG-1)*QG(3)+ATMPL(J,NG)*QG(4)
             PP(N) = Y0
-            KP(N) = NG-3
+            KP(N) = IMINMAX(J)-3
             GOTO 20
          ENDIF
 C GENERAL CASE - FIND 4 NEAREST POINTS IN GRAVITY THAT ARE IN THE TABLE.
@@ lines 106–112 @@
          DO K = IMINMAX(J)-3,IMIN(J),-1
             IF(GL.LT.ATMGL(K+2).AND.GL.GE.ATMGL(K+1))THEN
                KK = MAX(IMIN(J),K)
-               KK = MIN(NG-3,KK)
+               KK = MIN(IMINMAX(J)-3,KK)
                DO KKK = 1,4
                   QG(KKK) = ATMGL(KK+KKK-1)
                   PTAB(KKK) = ATMPL(J,KK+KKK-1)

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 branch ran in the two runs of 34 runs (26 of YREC's cases, from a 0.04 $M_\odot$ brown dwarf and the solar models to 9 $M_\odot$, rotation, and helium burning, and 8 made-up ones; one of these, with FCORE = 1.0, stopped in its first second) with KTTAU = 3, testsuite/Test_solarGS98_kurucz (10,264 times) and testsuite/Test_solarM22M_difkur (14,236 times), both with atmk1990p00.tab, at $\log T_\mathrm{eff}$ = 3.63 to 3.76 and $\log g$ = 2.70 to 4.57; no row with missing highest gravities was met, as expected for that table. A cool dwarf with one of the metal-poor tables ($[\mathrm{Fe/H}] \le -2.3$) would reach it; none was run. No fix was run. Intent: The comment (lines 104–105 ↗) says the change to IMINMAX was made "to prevent using -999"; only the loop's start was changed.

Definitions

subroutineSURFP

lines 1–155 · GitHub

$\log P$ at $T = T_\mathrm{eff}$ from the Kurucz atmosphere table, by interpolation in $\log g$ then $\log T_\mathrm{eff}$.

Uses modules ATMOS2_MOD, ATMPRT_MOD, FAC_MOD, IOUNITS_MOD

Called by ENVINT

Flowchart of its calls, in order (9)

for J = JJ,JJ+3 (DO 20), lines 66–138

if (ATMTL(J).GT.4.5D0), lines 70–87

else, lines 76–85

INTER3 81

if (GL.GE.ATMGL(NG-1)), lines 88–102

KSPLINE 95→KSPLINT 96

for K = IMINMAX(J)-3,IMIN(J),-1, lines 106–123

if (GL.LT.ATMGL(K+2).AND.GL.GE.ATMGL(K+1)), lines 107–122

KSPLINE 114→KSPLINT 115

KSPLINE 129→KSPLINT 130

KSPLINE 139→KSPLINT 140

Change log

The changes since rp_1

Changed since rp_1. Declarations only: COMMON blocks to modules (/ATMOS2/, /FAC/, /ATMPRT/ and /LUOUT/ to ATMOS2_MOD, FAC_MOD, ATMPRT_MOD and IOUNITS_MOD; the unused parameters NTC, NGC dropped), IMPLICIT NONE.

Changed since a8a5950. The out-of-table stop ends the run with exit status 1 (STOP 1) instead of a bare STOP, which exited with 0 (KNOWN_ISSUES.md, issue 38).

1 SUBROUTINE SURFP(TEFFL,GL,LPRT) 2C 10/14/26 MHP via Claude 3C code changed to convert COMMON/FAC/ into FAC_MOD. 4 USE FAC_MOD 5C 10/14/26 MHP via Claude 6C code changed to convert COMMON/ATMPRT/ into ATMPRT_MOD. 7 USE ATMPRT_MOD 8C 10/14/26 MHP via Claude 9C code changed to convert COMMON/ATMOS2/ into ATMOS2_MOD. 10 USE ATMOS2_MOD 11C 10/14/26 MHP via Claude 12C code changed to convert COMMON/IOUNITS/ into IOUNITS_MOD. 13 USE IOUNITS_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 DUM,FX,GL,PP,PTAB,QG,QGG,QS,QT,TEFFL,Y0,YG,YT 19 INTEGER J,JJ,K,KK,KKK,N,NG,NT 20 LOGICAL*4 LPRT 21 22C PARAMETERS NT AND NG FOR TABULATED SURFACE PRESSURES. 23 PARAMETER(NT=57,NG=11) 24C CHARACTER*256 FATM 25 26C MHP 8/25 Removed file names from common block 27C COMMON/ATMOS2/ATMPL(NT,NG),ATMTL(NT), 28C * ATMGL(NG),ATMZ,IOATM,FATM 29C 9/25/26 MHP via Claude 30C code changed to consolidate all I/O logical unit numbers, previously 31C scattered across COMMON/LUOUT/, COMMON/LUNUM/, COMMON/ALEX06/, 32C COMMON/ACDPTH/, COMMON/NEWOPAC/, COMMON/LOPAL95/, COMMON/ATMOS2/, 33C COMMON/ALATM03/, COMMON/OPALEOS/, and COMMON/SCV2/, into a single 34C COMMON/IOUNITS/ 35 DIMENSION QT(4),QG(4),PP(4),DUM(3),QGG(3),QS(3),PTAB(4), 36 * * YT(4),YG(4) 37 SAVE 38C SURFPL INTERPOLATES IN TEMPERATURE USING A 4-POINT LAGRANGIAN 39C INTERPOLATOR, AND INTERPOLATES IN GRAVITY THE SAME WAY IF 4 OR 40C MORE POINTS ARE AVAILABLE. IT WILL QUIT IF THE DESIRED DATA POINT 41C HAS TEFF OR LOG G MORE THAN ONE TABLE POINT FROM THE DATA. 42C 43C CHECK TO ENSURE THAT DATA IS WITHIN TABLE. 44 IF(TEFFL.LT.3.5D0 .OR. GL.LT.-0.5D0)THEN 45 WRITE(IOWR,911)TEFFL,GL 46 WRITE(ISHORT,911)TEFFL,GL 47 911 911 FORMAT(1X,'DESIRED ATMOSPHERE OUTSIDE TABLE RANGE'/ 48 * * ' LOG TEFF',F10.6,' LOG G',F10.6/' RUN STOPPED') 49C 10/30/26 MHP via Claude: bare STOP -> STOP 1, the run-time error code (issue 38). 50 STOP 1 51 ENDIF 52C TEMPERATURE INTERPOLATION FACTORS. 53 DO J = 1,NT 54 IF(TEFFL.LE.ATMTL(J))GOTO 10 55 END DO 56 J = NT 57 10 10 CONTINUE 58 JJ = MAX(1,J-2) 59 JJ = MIN(NT-3,JJ) 60 DO K = 1,4 61 QT(K) = ATMTL(JJ+K-1) 62 END DO 63C CALL INTRP2(QR,QT,TEFFL) 64 JP = JJ 65C GRAVITY INTERPOLATION FACTORS. 66 DO 20 J = JJ,JJ+3 67 N = J-JJ+1 68C CHECK IF 4 LOG VALUES AVAILABLE - OTHERWISE, USE 3 POINT LAGRANGIAN 69C OR LINEAR INTERPOLATION. 70 IF(ATMTL(J).GT.4.5D0)THEN 71 IF(ATMTL(J).GT.4.55D0)THEN 72C LINEAR INTERPOLATION 73 FX = (GL-ATMGL(NG-1))/(ATMGL(NG)-ATMGL(NG-1)) 74 PP(N)=ATMPL(J,NG-1)+FX*(ATMPL(J,NG)-ATMPL(J,NG-1)) 75 KP(N) = NG-1 76 ELSE 77C 3-POINT LAGRANGIAN INTERPOLATION. 78 DO K = 1,3 79 QS(K) = ATMGL(NG-3+K) 80 END DO 81 CALL INTER3(QS,QGG,DUM,GL) 82 PP(N)=ATMPL(J,NG-2)*QGG(1)+ATMPL(J,NG-1)*QGG(2)+ 83 * * ATMPL(J,NG)*QGG(3) 84 KP(N) = NG-2 85 ENDIF 86 GOTO 20 87 ENDIF 88 IF(GL.GE.ATMGL(NG-1))THEN 89C DESIRED LOG G ABOVE SECOND TO TOP TABLE LOG G - USE TOP 4 LOG G VALUES. 90 DO KK = 1,4 91 QG(KK)=ATMGL(NG-4+KK) 92 PTAB(KK) = ATMPL(J,NG-4+KK) 93 END DO 94C CALL INTRP2(QR,QG,GL) 95 CALL KSPLINE(QG,PTAB,YG) 96 CALL KSPLINT(QG,PTAB,YG,GL,Y0) 97C PP(N) = ATMPL(J,NG-3)*QG(1)+ATMPL(J,NG-2)*QG(2)+ 98C * ATMPL(J,NG-1)*QG(3)+ATMPL(J,NG)*QG(4) 99 PP(N) = Y0 100 KP(N) = NG-3 101 GOTO 20 102 ENDIF 103C GENERAL CASE - FIND 4 NEAREST POINTS IN GRAVITY THAT ARE IN THE TABLE. 104C G Somers, I changed NG to IMINMAX in the next line. This prevents the 105C code from using -999 to interpolate in some instances. 106 DO K = IMINMAX(J)-3,IMIN(J),-1 107 IF(GL.LT.ATMGL(K+2).AND.GL.GE.ATMGL(K+1))THEN 108 KK = MAX(IMIN(J),K) 109 KK = MIN(NG-3,KK) 110 DO KKK = 1,4 111 QG(KKK) = ATMGL(KK+KKK-1) 112 PTAB(KKK) = ATMPL(J,KK+KKK-1) 113 END DO 114 CALL KSPLINE(QG,PTAB,YG) 115 CALL KSPLINT(QG,PTAB,YG,GL,Y0) 116 PP(N) = Y0 117C CALL INTRP2(QR,QG,GL) 118C PP(N) = ATMPL(J,KK)*QG(1)+ATMPL(J,KK+1)*QG(2)+ 119C * ATMPL(J,KK+2)*QG(3)+ATMPL(J,KK+3)*QG(4) 120 KP(N) = KK 121 GOTO 20 122 ENDIF 123 END DO 124C DESIRED LOG G BELOW 2ND TABLE ENTRY -USE FIRST 4 POINTS. 125 DO K = 1,4 126 QG(K) = ATMGL(K+IMIN(J)-1) 127 PTAB(K) = ATMPL(J,K+IMIN(J)-1) 128 END DO 129 CALL KSPLINE(QG,PTAB,YG) 130 CALL KSPLINT(QG,PTAB,YG,GL,Y0) 131 PP(N) = Y0 132C CALL INTRP2(QR,QG,GL) 133C II = IMIN(J) 134C PP(N) = ATMPL(J,II)*QG(1)+ATMPL(J,II+1)*QG(2)+ 135C * ATMPL(J,II+2)*QG(3)+ATMPL(J,II+3)*QG(4) 136 KP(N) = IMIN(J) 137 20 20 CONTINUE 138C INTERPOLATE IN TEMPERATURE TO FIND CORRECT LOG P. 139 CALL KSPLINE(QT,PP,YT) 140 CALL KSPLINT(QT,PP,YT,TEFFL,Y0) 141C AP = PP(1)*QT(1)+PP(2)*QT(2)+PP(3)*QT(3)+PP(4)*QT(4) 142 AP = Y0 143 AT = TEFFL 144C WRITE OUT INFORMATION TO THE MODEL FILE. 145 IF (LPRT) THEN 146 WRITE(ISHORT,70) 147 WRITE(ISTOR,70) 14870 70 FORMAT('********PRESSURE AT T=TEFF INTERPOLATED FROM TABULATED' 149 * * , ' VALUES********') 150 WRITE(ISHORT,71) TEFFL,AP 151 WRITE(ISTOR,71) TEFFL,AP 15271 71 FORMAT(' ',20X,'LOG (Teff) =',F10.5,' LOG P =',F10.5) 153 ENDIF 154 RETURN 155 END

1"""Python version of surfp.f (YREC 7be4133): the pressure at T = Teff from 2the Kurucz model-atmosphere table (KTTAU = 3), the outer boundary of the 3envelope integration.""" 4 5import numpy as np 6 7from pyrec import atmos2_mod 8from pyrec import atmprt_mod 9from pyrec import fac_mod 10from pyrec import iounits_mod 11from pyrec.inter3 import inter3 12from pyrec.kspline import kspline 13from pyrec.ksplint import ksplint 14 15# PARAMETER (NT=57, NG=11): the table's temperatures and gravities 16# (Fortran line 23). 17NT = 57 18NG = 11 19 20 21def _spline4(qg, ptab, x): 22 """KSPLINE then KSPLINT on four points: the spline's value at x.""" 23 24 yg = np.zeros(5) 25 kspline(qg, ptab, yg) 26 return ksplint(qg, ptab, yg, x, 0.0) 27 28 29def surfp(teffl: float, gl: float, lprt: bool) -> None: 30 """log P at T = Teff, interpolated in log g at four table temperatures 31 around log Teff (linearly, by 3-point Lagrange or by a 4-point natural 32 spline, by temperature), then in log Teff by a 4-point natural spline. 33 34 Translates: surfp.f#SURFP (7be4133) 35 36 Arguments: teffl (log Teff), gl (log g); lprt, to print the result on 37 ISHORT and ISTOR. Sets atmprt_mod.ap (log P) and atmprt_mod.at 38 (log Teff), and fac_mod.jp and fac_mod.kp[1..4] (the table cell used). 39 None is returned. 40 41 """ 42 43 # The local arrays (Fortran lines 35–36). 44 qt = np.zeros(5) 45 pp = np.zeros(5) 46 qs = np.zeros(4) 47 qgg = np.zeros(4) 48 dum = np.zeros(4) 49 qg = np.zeros(5) 50 ptab = np.zeros(5) 51 atmpl = atmos2_mod.atmpl 52 atmtl = atmos2_mod.atmtl 53 atmgl = atmos2_mod.atmgl 54 55 # Question (see the page): only the lower limits are checked, and they 56 # are hard-coded (Fortran lines 43–51). 57 if teffl < 3.5 or gl < -0.5: 58 for unit in (iounits_mod.iowr, iounits_mod.ishort): 59 iounits_mod.write(unit, " DESIRED ATMOSPHERE OUTSIDE TABLE RANGE") 60 iounits_mod.write(unit, f" LOG TEFF{teffl:10.6f} LOG G{gl:10.6f}") 61 iounits_mod.write(unit, " RUN STOPPED") # format 911 62 raise SystemExit(1) # STOP 1, the run-time error code 63 64 # The four table temperatures around teffl (Fortran lines 52–64). 65 for j in range(1, NT + 1): 66 if teffl <= atmtl[j]: 67 break # to label 10 68 else: 69 j = NT 70 jj = max(1, j - 2) # label 10 71 jj = min(NT - 3, jj) 72 for k in range(1, 5): 73 qt[k] = atmtl[jj + k - 1] 74 fac_mod.jp = jj 75 76 # log P in log g at each of the four temperatures (loop 20, Fortran 77 # lines 65–137). 78 for j in range(jj, jj + 4): 79 n = j - jj + 1 80 if atmtl[j] > 4.5: 81 # Question (see the page): hot temperatures use only the 82 # highest gravities, without checking IMINMAX(J) (Fortran lines 83 # 70–87). 84 if atmtl[j] > 4.55: 85 fx = (gl - atmgl[NG - 1]) / (atmgl[NG] - atmgl[NG - 1]) 86 pp[n] = atmpl[j, NG - 1] + fx * (atmpl[j, NG] 87 - atmpl[j, NG - 1]) 88 fac_mod.kp[n] = NG - 1 89 else: 90 for k in range(1, 4): 91 qs[k] = atmgl[NG - 3 + k] 92 inter3(qs, qgg, dum, gl) 93 pp[n] = (atmpl[j, NG - 2] * qgg[1] 94 + atmpl[j, NG - 1] * qgg[2] 95 + atmpl[j, NG] * qgg[3]) 96 fac_mod.kp[n] = NG - 2 97 continue 98 99 if gl >= atmgl[NG - 1]: 100 # The top four gravities (Fortran lines 88–102). Suspected bug 101 # (see the page): IMINMAX(J) is not checked here. 102 for kk in range(1, 5): 103 qg[kk] = atmgl[NG - 4 + kk] 104 ptab[kk] = atmpl[j, NG - 4 + kk] 105 pp[n] = _spline4(qg, ptab, gl) 106 fac_mod.kp[n] = NG - 3 107 continue 108 109 # General case: the four valid gravities bracketing gl (Fortran 110 # lines 103–123). Suspected bug (see the page): NG-3, not 111 # IMINMAX(J)-3, bounds kk, and the search stops below IMINMAX(J)-1. 112 found = False 113 for k in range(fac_mod.iminmax[j] - 3, fac_mod.imin[j] - 1, -1): 114 if gl < atmgl[k + 2] and gl >= atmgl[k + 1]: 115 kk = max(fac_mod.imin[j], k) 116 kk = min(NG - 3, kk) 117 for kkk in range(1, 5): 118 qg[kkk] = atmgl[kk + kkk - 1] 119 ptab[kkk] = atmpl[j, kk + kkk - 1] 120 pp[n] = _spline4(qg, ptab, gl) 121 fac_mod.kp[n] = kk 122 found = True 123 break 124 if found: 125 continue 126 127 # The lowest four valid gravities (Fortran lines 124–136). 128 for k in range(1, 5): 129 qg[k] = atmgl[k + fac_mod.imin[j] - 1] 130 ptab[k] = atmpl[j, k + fac_mod.imin[j] - 1] 131 pp[n] = _spline4(qg, ptab, gl) 132 fac_mod.kp[n] = fac_mod.imin[j] 133 134 # log P in log Teff (Fortran lines 138–143). 135 atmprt_mod.ap = _spline4(qt, pp, teffl) 136 atmprt_mod.at = teffl 137 138 139 # Print to ISHORT and ISTOR if requested (Fortran lines 144–153). 140 if lprt: 141 for unit in (iounits_mod.ishort, iounits_mod.istor): 142 iounits_mod.write(unit, "********PRESSURE AT T=TEFF INTERPOLATED " 143 "FROM TABULATED VALUES********") # format 70 144 for unit in (iounits_mod.ishort, iounits_mod.istor): 145 iounits_mod.write(unit, f" {'':20s}LOG (Teff) ={teffl:10.5f}" 146 f" LOG P ={atmprt_mod.ap:10.5f}") # format 71