kcsurfp.f

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

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

Variables on this page (14)

ap, at, atmglc, atmplc, atmtlc, imin2, iminmax2, iowr, ishort, istor, jp, kp, ngc, ntc

Purpose

KCSURFP finds the gas pressure at the point where $T = T_\mathrm{eff}$ by interpolating in a table of model-atmosphere pressures as a function of $\log T_\mathrm{eff}$ and $\log g$. It is the version of SURFP for the newer Kurucz/Castelli tables (JNT 06/2014), whose size is $76\times11$ (NTC, NGC).

How it fits

Who calls it: only ENVINT, the envelope integration, which needs the pressure at the base of the atmosphere.
Inputs: TEFFL = $\log T_\mathrm{eff}$, GL = $\log g$, LPRT (print or not); the table in ATMOS2C_MOD (ATMPLC = $\log P$, grids ATMTLC, ATMGLC), filled by SETUPS, the module's only other user; in FAC_MOD, the per-temperature range of valid gravities IMIN2, IMINMAX2 (table entries that are missing in the grid).
Outputs: in ATMPRT_MOD, AP = $\log P$ and AT = $\log T_\mathrm{eff}$; ATMPRT_MOD is read by ENVINT, QATM, MASSLOSS and the other surface-pressure routines. It also sets JP and KP(4) in FAC_MOD (the table corner used). When LPRT is true it writes the result to the short output (ISHORT) and the model-storage unit (ISTOR), both from IOUNITS_MOD.
Calls: KSPLINE and KSPLINT (4-point natural splines), INTER3 (3-point Lagrange).
Modules: ATMPRT_MOD, ATMOS2C_MOD (with SETUPS), FAC_MOD (with SETUPS, SURFP), IOUNITS_MOD.

The method

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

Question 1 (line 53 ↗; lines 79–80 ↗). Hard-coded table limits. The bounds 3.5 and −0.5 and the thresholds 4.5 and 4.55 in $\log T$ (above which only the highest gravities are used) are written into the code; they belong to one table. Are they right for every Kurucz/Castelli table the code may read? No upper bounds are checked, and SURFBC has no switch to gray for KTTAU = 5 (surfbc.f, question 1), so above the tables' last temperature (4.699 in input/atmos/CastelliKurucz/) the spline extrapolates silently.
Take four table temperatures around $\log T_\mathrm{eff}$.
At each of them, interpolate $\log P$ in $\log g$:
for $\log T > 4.55$, linearly between the two highest gravities;

Question 2 (lines 80–85 ↗). Linear extrapolation in gravity for hot stars. For $\log T > 4.55$ the linear formula extrapolates freely in $\log g$ beyond the two highest gravities (the tables in input/atmos/CastelliKurucz/ have no pressure below $\log g$ = 4.5 from $\log T$ = 4.602 up). Is that acceptable for the stars that reach these temperatures?
for $4.5 < \log T \le 4.55$, by 3-point Lagrange (INTER3) over the three highest gravities;
otherwise with a 4-point natural spline over the four valid gravities around $\log g$ (the top four if $\log g$ is above the second-highest one, the lowest four valid ones if it is below them).
Interpolate the four results in $\log T_\mathrm{eff}$ with a 4-point natural spline: AP.

The comment (lines 46–49 ↗) is from an older routine: It calls the routine "SURFPL," says it interpolates in temperature and gravity with a 4-point Lagrange interpolator (the routine uses natural splines), and that "IT WILL QUIT IF THE DESIRED DATA POINT HAS TEFF OR LOG G MORE THAN ONE TABLE POINT FROM THE DATA" (only the lower limits are checked, line 53 ↗).

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 pressures at the highest gravities.

(lines 97–111 ↗; lines 112–132 ↗) SETUPS sets IMINMAX2(J) below NGC when a temperature's highest-gravity entries are missing (pressure $\le 0$, e.g. −999) (setups.f, lines 327–358 ↗). The general search (loop over K from IMINMAX2(J)−3 down to IMIN2(J)) then covers $\log g$ only up to ATMGLC(IMINMAX2(J)-1). A $\log g$ between that and ATMGLC(NGC-1) matches no interval and falls through to the lowest four valid gravities, a long extrapolation; and a $\log g$ at or above ATMGLC(NGC-1) takes the top four gravities, missing entries included. Correct with the tables in input/atmos/CastelliKurucz/, which have all their highest gravities below $\log T$ = 4.5 (checked: every row of the 16 tables), but not with a table like the metal-poor Kurucz 1990 ones, where the same code in SURFP gives wrong pressures (surfp.f, suspected bug 1).

A possible fix:

kcsurfp.f

@@ lines 94–104 @@
             ENDIF
             GOTO 20
          ENDIF
-         IF(GL.GE.ATMGLC(NGC-1))THEN
+C 11/09/26 KC via Claude
+C code changed to use IMINMAX2(J), the highest gravity with a pressure,
+C in place of NGC, 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.ATMGLC(IMINMAX2(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)=ATMGLC(NGC-4+KK)
-               PTAB(KK) = ATMPLC(J,NGC-4+KK)
+               QG(KK)=ATMGLC(IMINMAX2(J)-4+KK)
+               PTAB(KK) = ATMPLC(J,IMINMAX2(J)-4+KK)
             END DO
 C            CALL INTRP2(QR,QG,GL)
             CALL KSPLINE(QG,PTAB,YG)
@@ lines 106–112 @@
 C            PP(N) = ATMPLC(J,NGC-3)*QG(1)+ATMPLC(J,NGC-2)*QG(2)+
 C     *              ATMPLC(J,NGC-1)*QG(3)+ATMPLC(J,NGC)*QG(4)
             PP(N) = Y0
-            KP(N) = NGC-3
+            KP(N) = IMINMAX2(J)-3
             GOTO 20
          ENDIF
 C GENERAL CASE - FIND 4 NEAREST POINTS IN GRAVITY THAT ARE IN THE TABLE.
@@ lines 115–121 @@
          DO K = IMINMAX2(J)-3,IMIN2(J),-1
             IF(GL.LT.ATMGLC(K+2).AND.GL.GE.ATMGLC(K+1))THEN
                KK = MAX(IMIN2(J),K)
-               KK = MIN(NGC-3,KK)
+               KK = MIN(IMINMAX2(J)-3,KK)
                DO KKK = 1,4
                   QG(KKK) = ATMGLC(KK+KKK-1)
                   PTAB(KKK) = ATMPLC(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): No case of YREC's uses KTTAU = 5, so a made-up run was added: testsuite/Test_solarGS98_kurucz with KTTAU = 5 and atmk2004p00.tab, as input/atmos/CastelliKurucz/README asks (line 15: "You NEED TO CHANGE KTTAU to 5"). The branch ran 10,156 times, at $\log T_\mathrm{eff}$ = 3.63 to 3.76 and $\log g$ = 2.70 to 4.57, and never met a row with missing highest gravities, as the check of the tables above found. A table with such rows would reach it; none is shipped. No fix was run. Intent: The comment (lines 113–114 ↗) says the change to IMINMAX2 was made "to prevent using -999"; only the loop's start was changed.

Definitions

subroutineKCSURFP

lines 1–164 · GitHub

Pressure at $T = T_\mathrm{eff}$ from the Kurucz/Castelli grid: splines in $\log g$ at four temperatures, then in $\log T_\mathrm{eff}$ (above).

Uses modules ATMOS2C_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 75–147

if (ATMTLC(J).GT.4.5D0), lines 79–96

else, lines 85–94

INTER3 90

if (GL.GE.ATMGLC(NGC-1)), lines 97–111

KSPLINE 104→KSPLINT 105

for K = IMINMAX2(J)-3,IMIN2(J),-1, lines 115–132

if (GL.LT.ATMGLC(K+2).AND.GL.GE.ATMGLC(K+1)), lines 116–131

KSPLINE 123→KSPLINT 124

KSPLINE 138→KSPLINT 139

KSPLINE 148→KSPLINT 149

Change log

The changes since rp_1

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

1"""Python version of kcsurfp.f (YREC development, 7be4133).""" 2 3import numpy as np 4 5from pyrec import atmos2c_mod 6from pyrec import atmprt_mod 7from pyrec import fac_mod 8from pyrec import iounits_mod 9from pyrec.inter3 import inter3 10from pyrec.kspline import kspline 11from pyrec.ksplint import ksplint 12 13# Parameters NTC and NGC for the tabulated surface pressures (Fortran 14# line 27). 15NTC = 76 16NGC = 11 17 18 19def _spline4(qg, ptab, x): 20 """KSPLINE then KSPLINT on four points: the spline's value at x.""" 21 22 yg = np.zeros(5) 23 kspline(qg, ptab, yg) 24 return ksplint(qg, ptab, yg, x, 0.0) 25 26 27def kcsurfp(teffl: float, gl: float, lprt: bool) -> None: 28 """Gas pressure at T = Teff from the Kurucz/Castelli atmosphere table. 29 30 The version of SURFP for the newer Kurucz/Castelli tables (76 x 11): 31 log P is interpolated in log g at four table temperatures around 32 log Teff (linearly, by 3-point Lagrange or by a 4-point natural spline, 33 depending on the temperature), then in log Teff by a 4-point natural 34 spline. Called only by ENVINT. 35 36 Translates: kcsurfp.f#KCSURFP (7be4133) 37 38 Parameters 39 ---------- 40 teffl : float 41 log Teff. 42 gl : float 43 log g. 44 lprt : bool 45 Write the result to the short output and the model-storage unit. 46 47 Sets ``atmprt_mod.ap`` (log P) and ``atmprt_mod.at`` (log Teff), and 48 ``fac_mod.jp`` and ``fac_mod.kp[1..4]`` (the table corner used). 49 Nothing is returned. 50 51 """ 52 53 # The local arrays (Fortran lines 43–44). 54 qt = np.zeros(5) 55 pp = np.zeros(5) 56 qs = np.zeros(4) 57 qgg = np.zeros(4) 58 dum = np.zeros(4) 59 qg = np.zeros(5) 60 ptab = np.zeros(5) 61 atmplc = atmos2c_mod.atmplc 62 atmtlc = atmos2c_mod.atmtlc 63 atmglc = atmos2c_mod.atmglc 64 65 # Check that the point is within the table (Fortran lines 52–60). 66 # Question (see the page): hard-coded bounds; no upper bounds checked. 67 if teffl < 3.5 or gl < -0.5: 68 for unit in (iounits_mod.iowr, iounits_mod.ishort): 69 iounits_mod.write(unit, " DESIRED ATMOSPHERE OUTSIDE TABLE RANGE") 70 iounits_mod.write(unit, f" LOG TEFF{teffl:10.6f} LOG G{gl:10.6f}") 71 iounits_mod.write(unit, " RUN STOPPED") # format 911 72 raise SystemExit(1) # STOP 1, the run-time error code 73 74 # Temperature interpolation: the four table points around teffl 75 # (Fortran lines 61–73). 76 for j in range(1, NTC + 1): 77 if teffl <= atmtlc[j]: 78 break # to label 10 79 else: 80 j = NTC 81 jj = max(1, j - 2) # label 10 82 jj = min(NTC - 3, jj) 83 for k in range(1, 5): 84 qt[k] = atmtlc[jj + k - 1] 85 fac_mod.jp = jj 86 87 # Gravity interpolation at each of the four temperatures (loop 20, 88 # Fortran lines 74–146). 89 # Hot temperatures: the highest gravities only (Fortran lines 79–96). 90 for j in range(jj, jj + 4): 91 n = j - jj + 1 92 if atmtlc[j] > 4.5: 93 if atmtlc[j] > 4.55: 94 # Linear interpolation (extrapolation) in the two highest 95 # gravities. Question (see the page): unbounded in log g. 96 fx = (gl - atmglc[NGC - 1]) / (atmglc[NGC] - atmglc[NGC - 1]) 97 pp[n] = (atmplc[j, NGC - 1] 98 + fx * (atmplc[j, NGC] - atmplc[j, NGC - 1])) 99 fac_mod.kp[n] = NGC - 1 100 else: 101 # 3-point Lagrange interpolation in the three highest. 102 for k in range(1, 4): 103 qs[k] = atmglc[NGC - 3 + k] 104 inter3(qs, qgg, dum, gl) 105 pp[n] = (atmplc[j, NGC - 2] * qgg[1] 106 + atmplc[j, NGC - 1] * qgg[2] 107 + atmplc[j, NGC] * qgg[3]) 108 fac_mod.kp[n] = NGC - 2 109 continue 110 111 # Suspected bug (see the page): this branch and the search below use 112 # NGC, not IMINMAX2(J), so missing top entries can enter (Fortran 113 # lines 97–111). 114 if gl >= atmglc[NGC - 1]: 115 # Above the second-highest gravity: the top four. 116 for kk in range(1, 5): 117 qg[kk] = atmglc[NGC - 4 + kk] 118 ptab[kk] = atmplc[j, NGC - 4 + kk] 119 pp[n] = _spline4(qg, ptab, gl) 120 fac_mod.kp[n] = NGC - 3 121 continue 122 123 # General case: the four nearest valid gravities (Fortran lines 124 # 112–132). 125 found = False 126 for k in range(fac_mod.iminmax2[j] - 3, fac_mod.imin2[j] - 1, -1): 127 if atmglc[k + 2] > gl >= atmglc[k + 1]: 128 kk = max(fac_mod.imin2[j], k) 129 kk = min(NGC - 3, kk) 130 for kkk in range(1, 5): 131 qg[kkk] = atmglc[kk + kkk - 1] 132 ptab[kkk] = atmplc[j, kk + kkk - 1] 133 pp[n] = _spline4(qg, ptab, gl) 134 fac_mod.kp[n] = kk 135 found = True 136 break 137 if found: 138 continue 139 140 # Below the second valid gravity: the lowest four valid ones 141 # (Fortran lines 133–145). 142 for k in range(1, 5): 143 qg[k] = atmglc[k + fac_mod.imin2[j] - 1] 144 ptab[k] = atmplc[j, k + fac_mod.imin2[j] - 1] 145 pp[n] = _spline4(qg, ptab, gl) 146 fac_mod.kp[n] = fac_mod.imin2[j] 147 148 # Interpolate in temperature to find log P (Fortran lines 147–152). 149 atmprt_mod.ap = _spline4(qt, pp, teffl) 150 atmprt_mod.at = teffl 151 152 # Write the result to the model file (Fortran lines 153–162). 153 if lprt: 154 for unit in (iounits_mod.ishort, iounits_mod.istor): 155 iounits_mod.write(unit, "********PRESSURE AT T=TEFF INTERPOLATED " 156 "FROM TABULATED VALUES********") # format 70 157 for unit in (iounits_mod.ishort, iounits_mod.istor): 158 iounits_mod.write(unit, f" {'':20s}LOG (Teff) ={teffl:10.5f}" 159 f" LOG P ={atmprt_mod.ap:10.5f}") # format 71