getopal95.f

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

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

Variables on this page (31)

bd, bdd, bt, btt, bx, bz, cappa, cappa2, ishort, jd, jt, jx, jz, lacout, ladon, lexh, lexth, lextl, numd, numt, numx, numz, r0gr, rho_t63l, rho_t63ll, rlh, tgr, xop95, xxg, ztab, zzg

Purpose

GETOPAL95 is the driver for the OPAL 1995 opacity tables ("7/98 MHP"): given $\rho$, $T$, $X$, $Z$, it finds the interpolation stencil and weights in $T$, $R = \rho/T_6^3$, $X$ and $Z$, handles the table's empty corner at high density, and calls a 2D, 3D or 4D interpolator depending on how far the composition is from the tables already prepared.

How it fits

Who calls it: only GETOPAC, when LOPAL95 is set (with $Z$ = 0.1 in metal-rich helium-burning cores).
Inputs: DL, TL, X, Z.
What it calls: INTERP (four-point Lagrange weights with derivative weights, for $T$ and $R$; the file's header speaks of "CUBIC SPLINE INTERPLOATION," which is stale: Settled from the code, Nov 9, 2026: lines 1–2 ↗, interp.f, lines 16–18 ↗), INTRP2 (weights only, for $X$ and $Z$), and one of OP952D (2D, in the table for the surface $X$ and $Z$), OP953D (3D, in the tables for the model's $Z$) or OP954D (4D, in the full set).
Outputs: O, OL, QOD, QOT, as from GETOPAC.
Modules (COMMON blocks until rp_1): the tables LLOT95A_MOD (the full set CAPPA2, indexed by a combined X–Z table number, and the grids TGR, R0GR, XXG, ZZG; read by LL95TBL), LLOT95_MOD (CAPPA at the model's $Z$, ZTAB, made by OP95ZTAB), LLOT95E_MOD (CAPPAE at the surface $X$, XOP95 (formerly XTAB), made by OP95XTAB); OP95INTERP_MOD, which merges the former /OP95INDX/ (indices), /OP95FACT/ (weights) and /OP95EXT/ (extrapolation flags), shared with the OP95*D routines, $\log R$ now named RHO_T63L (formerly RL, which clashed with the radius elsewhere) and its edge value RHO_T63LL (formerly RLL); ACDPTH_MOD for the flags LADON, LACOUT (acoustic-depth output; shared with CALCAD, ENVINT, WRTOUT); IOUNITS_MOD for ISHORT, the .short file that the stop reports to.

Why three interpolators. SETUPOPAC, through LL95TBL (which calls OP95ZTAB at ZOPAL951) and OP95XTAB (at XENV), interpolates the full table set ahead of time: CAPPA at one $Z$ for every $X$ node, and CAPPAE from it at the surface $X$. Most calls (the envelope, the atmosphere, and the outer convection zone, which share the surface composition) then need only the 2D interpolation in $T$ and $R$; the 3D path is needed where hydrogen burning has changed $X$, and the 4D path where $Z$ differs too, for example after settling or in a helium-burning core.

The method

$\log R = \log\rho - 3(\log T - 6)$ (RHO_T63L). The run stops if $\log T$ is outside the table. Below the smallest $\log R$, the weights are computed at the table edge and the value is extrapolated linearly later (LEXTL).
Stencils in $\log T$ (JT) and $\log R$ (JD), each starting from the last one found.
Empty corner. Values of 9.9 or more mark missing entries. If the stencil's high-density, high-temperature corner is missing, each of the four temperatures gets its own density stencil (LSHIFT), and where even that is outside the data, its weights are computed at the last valid density and the value is extrapolated linearly (LEXH, RLH). The test reads the $X$ table JX(1,1) of the previous call, since the current $X$ stencil is found only later; this does not matter, since every table of a file has the same empty region: the last one to five densities of the temperatures from TGR(58) up, set to 9.999 by LL95TBL, and in the OP files of the recent mixtures also every density above $\log T = 8.0$ (bug 2). For the same reason the search for valid densities stops at the 14th density or above, so the start of a shifted stencil, JD(I) = J - 3, is at least 11 and never zero or negative; in a row with no data it finds nothing, and JD(I) keeps its old value (settled from the code, Nov 9, 2026: lines 177–209 ↗, ll95tbl.f, lines 140–171 ↗). At low $X$ and the lowest temperatures (X < 0.2, JT ≤ 5) the comment postpones a check of the low-density, low-temperature empty region to the 2D, 3D or 4D branch, but no branch makes it, and the flag LCHK set there is unused? (lines 162–172 ↗; bug 1 below).
Which interpolation. If $Z$ is within a relative $10^{-4}$ of the model table's $Z$: 2D if $X$ is within $10^{-4}$ of the surface table's $X$, else 3D in $X$, $T$, $R$. Otherwise 4D. The 4D path is forced when acoustic-depth output is on (LADON and LACOUT, "JVS 04/11").
Missing high-$X$ tables. The $X$ nodes are 0, 0.1, 0.2, 0.35, 0.5, 0.7, 0.8, 0.9, 0.95, 1. For $X \ge 0.8$: if $0.04 \le Z < 0.1$, the $X = 0.95$ table is absent and the stencil is 0.7, 0.8, 0.9 and the last node, set to $1 - Z$; if $Z \ge 0.1$, the last two are absent and the stencil is 0.5, 0.7, 0.8, 0.9. (In 3D the test is on the model table's $Z$, in 4D on each $Z$ node. OP95ZTAB stops the run for $Z > 0.08$, so the 3D case $Z \ge 0.1$ is never reached; the 4D path can still meet the node $Z = 0.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. (not reached) `LSHIFT` and `LEXTH` left over from the last call.

(lines 167–176 ↗, line 214 ↗) When X < 0.2 and JT ≤ 5, GOTO 60 jumps past the lines that reset LSHIFT and LEXTH to false, so their values from the previous call (the routine has SAVE) decide the density stencil; on the very first call LSHIFT has no value at all. If the previous call was in the empty corner (LSHIFT true), the density weights are computed on the stencils JD(2..4) and with the flags LEXH of that call, at other temperatures, and the opacity is interpolated or extrapolated from the wrong densities. This needs a hydrogen-poor gas ($X$ < 0.2) at the table's first temperatures, which GETOPAC reaches only in the ramp between TMOLMIN and TMOLMAX (by default $\log T$ = 4.0 to 4.1), so a hydrogen-poor surface; hydrogen-rich runs do not take it.

A possible fix:

getopal95.f

@@ lines 164–169 @@
 C 4X4 MATRIX OF T,RHO AROUND IT OUTSIDE THE TABLE.
 C CHECK LOW RHO, LOW T REGIME; SINCE THE EMPTY REGION DEPENDS ON Z AT
 C LOW X, POSTPONE UNTIL YOU KNOW IF IT IS A 2D, 3D, OR 4D INTERPOLATION.
+C 11/09/26 KC via Claude
+C code changed to reset LSHIFT and LEXTH before the low-X, low-T
+C branch, which skips the corner check.
+      LSHIFT = .FALSE.
+      LEXTH = .FALSE.
       IF(X.LT.0.2D0 .AND. JT. LE. 5)THEN
          LCHK = .TRUE.
          GOTO 60
@@ lines 172–179 @@
       ENDIF
 C CHECK HIGH RHO, HIGH T REGIME(BECAUSE OF THE TABLE GEOMETRY, THE
 C POINT WILL FALL OUTSIDE THE TABLE FIRST AT THE HIGHEST T).
-      LSHIFT = .FALSE.
-      LEXTH = .FALSE.
       IF(CAPPA(JX(1,1),JT+3,JD(1)+3) .GE. 9.9D0) THEN
          LSHIFT = .TRUE.
 C NOW THERE ARE TWO POSSIBLE SOLUTIONS: WE CAN USE A DIFFERENT SET

(the two resets at lines 175–176 go).

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 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 them the AAG21 case of suspected bug 2, YREC's 1 $M_\odot$ helium-burning case run with FLIV95 changed to op05/AAG2021zx0.0187.tab). YREC's cases and documentation hold no hydrogen-poor surface, so no run reaches it yet. Intent: The comment postpones a check that no branch makes, and LCHK, set here, is read nowhere in src/, so the skipped reset looks like a side effect of the unfinished check rather than a design.

2. (changes results, large) Checked by KC (Nov 10, 2026) Empty temperatures in the OP tables above $\log T = 8.0$.

(lines 113–131 ↗, lines 177–209 ↗; op952d.f, lines 39–67 ↗) In the OP tables of the recent mixtures (AAG2021zx0.0187.tab, Magg22_met_zx0.0226.tab, Magg22_phot_zx0.0225.tab, and the other AAG21*, Magg22* and Lodders25* files in input/opacity/atomic/op05/), the last four temperatures of every table ($\log T$ = 8.1, 8.3, 8.5, 8.7) hold 9.999 at all densities; GS98.OP17 has data there. From $\log T$ = 7.9 up, the $T$ stencil includes such a row. The empty-corner search finds no valid density in it, so its JD(I) keeps an old value, LEXH(I) is set, and the interpolators take $\log\kappa \approx 9.999$ at that temperature (the extrapolation adds nothing, since the slope of a constant row is zero). Between $\log T$ = 7.9 and 8.0 the Lagrange weight of the empty row is down to about $-0.064$, so $\log\kappa$ is too low by up to about 0.7 (and its $T$ derivative is off); above 8.0, with two or more empty rows in the stencil, the value is meaningless. The run neither stops nor warns, since the range check uses TGR(NUMT) = 8.7. The README of the tables says the hydrogen-burning tables are also required for helium-burning models (input/opacity/atomic/README, line 1), so a run with these tables reaches this path in a helium-burning core, and at the tip of the red giant branch; how much conduction hides it in degenerate cores is not traced. The test suite uses these tables only in three solar models (Test_solarAAG21_dif, Test_solarM22M_dif, Test_solarM22M_difkur), which stay below this temperature, and GS98.OP17 in the others.

A possible fix, a choice for the developers: stop with a message when the stencil's last temperature has no data (CAPPA(JX(1,1),JT+3,1) ≥ 9.9), as below (not run), or switch to another table above $\log T = 8.0$. The run below measured a third version, kept within the rows that hold data and extrapolating in $T$ above 8.0, used only to measure.

getopal95.f

@@ lines 170–175 @@
       ELSE
          LCHK = .FALSE.
       ENDIF
+C 11/10/26 KC via Claude
+C code changed to stop, with a message, when the temperature stencil
+C reaches a row that holds no data even at the lowest density, as the
+C recent OP tables do above log T = 8.0, instead of interpolating 9.999.
+      IF(CAPPA(JX(1,1),JT+3,1) .GE. 9.9D0) THEN
+         WRITE(*,6)TL
+         WRITE(ISHORT,6)TL
+ 6       FORMAT(' LOG T OF',F11.6,' NEEDS A ROW OF THE OPAL 95 TABLE WITH NO DATA')
+         STOP 1
+      ENDIF
 C CHECK HIGH RHO, HIGH T REGIME(BECAUSE OF THE TABLE GEOMETRY, THE
 C POINT WILL FALL OUTSIDE THE TABLE FIRST AT THE HIGHEST T).
       LSHIFT = .FALSE.

How the code itself shows it. The stencil search assumes that every temperature has data at the lowest densities. Its comment says the point "WILL FALL OUTSIDE THE TABLE FIRST AT THE HIGHEST T" (lines 173–174 ↗), and the search for a valid density walks down from IDD+2 (lines 198–204 ↗) with no branch for a row that has none, so JD(I) keeps its previous value. The reader LL95TBL builds the empty region as that corner: In the last 13 temperatures it reads one to five fewer densities and sets only the highest ones to 9.999D0 (ll95tbl.f, lines 140–171 ↗). In AAG2021zx0.0187.tab, the rows at $\log T$ = 8.10, 8.30, 8.50, and 8.70 hold 9.999 in every column of all 126 tables, where GS98.OP17 has values (−0.782 at 8.10 in the first table). For a temperature outside the table the routine stops: IF(TL.LT.TGR(1) .OR. TL.GT.TGR(NUMT)) writes "OUT OF OPAL 95 TABLE RANGE" and ends with STOP 1 (lines 90–96 ↗). A row with no data is outside the table in the same way, but these files give TGR(NUMT) = 8.70, so the check passes.

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): YREC's 1 $M_\odot$ helium-burning case (examples/run_from_zahb_to_tahb/ZAHB_m1p0feh+0p0_TAHB, which uses GS98.OP17), run with only FLIV95 changed to op05/AAG2021zx0.0187.tab, the table that input/opacity/atomic/README recommends for the AAG21 mixture, takes the path from the start: An empty row was in the $T$ stencil about 113 million times, at $\log T$ = 7.90 to 8.34, with Lagrange weights from $-0.064$ to 1.03. The run neither stopped nor warned. Rerun with the stencil kept within the rows that hold data (one possible fix, extrapolating in $T$ above 8.0; used here only to measure), the same case changes well beyond the test suite's 2% tolerance. Without the fix, the convective core (the track's Mco_core) is 0.059 $M_\odot$ when central $Y$ falls below 0.9 and 0.267 $M_\odot$ below 0.01; with it, 0.123 and 0.151 $M_\odot$. Central $Y$ reaches 0.5 after 0.038 Gyr without the fix and 0.106 Gyr with it; the phase lasts 0.232 Gyr (2,256 models) without it and 0.222 Gyr (1,824 models) with it; the core ends with C/O = 0.21/0.78 without it and 0.34/0.65 with it. With the fix, the run follows the unchanged GS98.OP17 case closely (core 0.119 and 0.151 $M_\odot$, central $Y$ = 0.5 at 0.105 Gyr; these three from the pilot runs of Nov 9). In YREC's own runs, which use GS98.OP17, the stencil reaches $\log T \ge 7.9$ in six cases (up to 31 million times, in examples/run_from_zahb_to_tahb/ZAHB_m0p8feh+0p0_TAHB; also testsuite/Test_m1p0feh+0p0GS98_base_HeIgnite, up the giant branch to helium ignition) but never meets an empty row, since that table has data there. On the AAG21 case, the fix moves $\log L$ by up to 0.27 dex, against $1.4 \times 10^{-6}$ dex for a null change (nearly every opacity, all with $3 \le \log T \le 9$ and $-6 \le \log\rho \le 9.75$, times $1 + 10^{-10}$), the noise floor.

3. (never occurred) The wrong array in the empty-corner search.

(lines 193–199 ↗) The search for valid densities at each temperature tests CAPPA2(IXX,ITT,J), with IXX = JX(1,1), an index into the $X$ tables of CAPPA, while the lines above it test CAPPA. CAPPA2's first index is the combined X–Z table number (LLOT95A_MOD; op95ztab.f, lines 56–64 ↗), so IXX picks the $Z = 0$ table with that $X$ index. The result is right today, since every table has the same empty corner (ll95tbl.f, lines 140–171 ↗); a table set whose empty region depends on the composition would get the stencil of the wrong table.

A possible fix:

getopal95.f

@@ lines 196–202 @@
             ITT = JT+I-1
             IF(CAPPA(IXX,ITT,IDD+3).GT.9.9D0)THEN
                DO J = IDD+2,1,-1
-                  IF(CAPPA2(IXX,ITT,J).LE.9.9D0)THEN
+C 11/09/26 KC via Claude
+C code changed to search CAPPA, the model-Z tables that IXX indexes,
+C instead of the full set CAPPA2.
+                  IF(CAPPA(IXX,ITT,J).LE.9.9D0)THEN
                      JD(I) = J - 3
                      IDD = J - 3
                      GOTO 50

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 search ran in five of 34 runs, all helium-burning or giant cores ($\log T$ = 7.3 to 8.34): 113 million times in the AAG21 case of suspected bug 2 (examples/run_from_zahb_to_tahb/ZAHB_m1p0feh+0p0_TAHB with FLIV95 changed to op05/AAG2021zx0.0187.tab), 25 million in testsuite/Test_m1p0feh+0p0GS98_base_HeIgnite, and fewer in three more of YREC's cases; CAPPA2 never gave a different answer from CAPPA, so results are unchanged today, as above. Intent: This is the routine's only read of CAPPA2, and the test two lines above reads CAPPA; elsewhere CAPPA2 is always indexed by the combined X–Z table number (op95ztab.f, lines 62–64 ↗, op954d.f, lines 55–58 ↗), never by an $X$ index such as IXX, so this looks like a slip.

Definitions

subroutineGETOPAL95

lines 17–430 · GitHub

Driver for the OPAL 1995 opacities: stencils and weights in $T$, $R$, $X$, $Z$, the empty-corner handling, and the choice of 2D, 3D or 4D interpolation (method above).

Uses modules ACDPTH_MOD, IOUNITS_MOD, LLOT95A_MOD, LLOT95E_MOD, LLOT95_MOD, OP95INTERP_MOD

Called by GETOPAC ×2

Flowchart of its calls, in order (11)

INTERP 137

if (.NOT.LSHIFT), lines 214–230

if (.NOT.LEXTL), lines 218–219

INTERP 219

else, lines 220–223

INTERP 222

else, lines 231–250

for J = 1,4, lines 232–249

if (.NOT.LEXH(J)), lines 238–239

INTERP 239

else, lines 240–244

INTERP 242

if (ABS(Z-ZTAB)/MAX(ZTAB,1.0D-6).LE.1.0D-4), lines 255–326

if (ABS(X-XOP95).LE.1.0D-4), lines 256–259

OP952D 258

else, lines 260–325

INTRP2 319→OP953D 323

INTRP2 358

for K = 1,4, lines 363–426

INTRP2 422

OP954D 427

Change log

The changes since rp_1

Changed since rp_1. Declarations and names only: COMMON blocks to modules (the three interpolation blocks merged into OP95INTERP_MOD), IMPLICIT NONE, RL, RLL and XTAB renamed RHO_T63L, RHO_T63LL and XOP95, and the unused parameters NUMXZ and JSON removed. The computation is unchanged; no question settled.

Changed since a8a5950. The stop for $\log T$ outside the table now writes its message to the .short file as well as the terminal (the routine now uses IOUNITS_MOD, for ISHORT), and exits with status 1, the code for run-time errors, instead of a bare STOP (KNOWN_ISSUES.md, issue 38, which this stop prompted).

1C THIS IS THE INTERPOLATION FACILITY FOR THE LIVERMORE OPACITY TABLES 2C USING CUBIC SPLINE INTERPLOATION SCHEME. (FOR OPAL95) 3C -------------------------------------------------------------- 4C SUBROUTINE GETOPAL95(DL,TL,X,Z,O,OL,QOD,QOT) DRIVER- GET CAPPA 5C SUBROUTINE OP952D(O,OL,QOD,QOT) 2D INTERPOLATION IN (T,RHO) 6C SUBROUTINE OP953D(O,OL,QOD,QOT) 3D INTERPOLATION IN (X,T,RHO) 7C SUBROUTINE OP954D(O,OL,QOD,QOT) 4D INTERPOLATION IN (Z,X,T,RHO) 8C SUBROUTINE LL95TBL READ IN ALL TABLES 9C SUBROUTINE OP95ZTAB(Z) GENERATE TABLE AT FIXED Z 10C SUBROUTINE OP95XTAB(X) GENERATE TABLE AT FIXED X,Z 11C 12C 13C 14CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 15C GETOPAL95 16CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 17 SUBROUTINE GETOPAL95(DL,TL,X,Z,O,OL,QOD,QOT) 18C 10/14/26 MHP via Claude 19C code changed to convert COMMON/OP95INTERP/ into OP95INTERP_MOD. 20 USE OP95INTERP_MOD 21C 10/14/26 MHP via Claude 22C code changed to convert COMMON/LLOT95A/ into LLOT95A_MOD. 23 USE LLOT95A_MOD 24C 10/14/26 MHP via Claude 25C code changed to convert COMMON/ACDPTH/ into ACDPTH_MOD. 26 USE ACDPTH_MOD 27C 10/14/26 MHP via Claude 28C code changed to convert COMMON/LLOT95E/ into LLOT95E_MOD. 29 USE LLOT95E_MOD 30C 10/14/26 MHP via Claude 31C code changed to convert COMMON/LLOT95/ into LLOT95_MOD. 32 USE LLOT95_MOD 33C 10/20/26 MHP via Claude 34C code changed to replace implicit typing with IMPLICIT NONE and 35C explicit declarations. 36C 10/30/26 MHP via Claude 37C code changed to USE IOUNITS_MOD for ISHORT, so the error stop below 38C also reports to the .short file (issue 38). 39 USE IOUNITS_MOD 40 IMPLICIT NONE 41 REAL*8 A,B,C,DL,O,OL,QOD,QOT,RR,TL,X,Z,ZZ 42 INTEGER I,IDD,ITT,IXX,J,JXX,K,NUMD,NUMT,NUMX,NUMZ 43 LOGICAL*4 LCHK,LSHIFT 44C 7/98 MHP DRIVER ROUTINE FOR OPACITY GIVEN RHO,T, X, AND Z. 45 PARAMETER (NUMT=70) 46 PARAMETER (NUMD=19) 47 PARAMETER (NUMX=10) 48 PARAMETER (NUMZ=13) 49C JVS Need this one too: 50C FULL SET OF TABLES: OPACITY AS A FUNCTION OF Z AND X, T, RHO/T6**3 51C TABLES ARE INCREMENTED IN SETS OF NZ*NX. SO THE TABLES FOR THE 52C THIRD METAL ABUNDANCE (3 X 10**-4)BEGIN AT TABLE 21 AND END AT TABLE 30. 53C FOR THE HIGH VALUES OF Z, THE NUMBER OF X TABLES IS NOT THE SAME (I.E. 54C X<0.9 FOR Z=0.1). 55C FOR EACH COMPOSITION A FULL GRID IN (T,RHO/T6**3) IS RETAINED. 56C AS ABOVE FOR THE MODEL Z. 57C AS ABOVE FOR DESIRED SURFACE VALUE OF X. 58C INDICES FOR INTERPOLATION IN Z,X,T, AND R; INTERPOLATION FACTORS FOR 59C Z,X,T, AND R (AND DERIVATIVE FACTORS FOR T AND RHO); AND DATA FOR 60C LINEAR EXTRAPOLATION WHEN OUTSIDE THE TABLE IN RHO. 61C 10/9/26 MHP via Claude 62C code changed to merge COMMON/OP95INDX/, COMMON/OP95FACT/, and 63C COMMON/OP95EXT/ into COMMON/OP95INTERP/: all three are small 64C interpolation-support arrays (not table data, unlike LLOT95/ 65C LLOT95E/LLOT95A, which hold the actual 2D/3D/4D opacity tables and 66C are intentionally left separate), with no member name colliding 67C with any dummy argument across the 6 files that declare any of 68C them. 69C 10/13/26 MHP via Claude 70C code changed to rename COMMON/OP95INTERP/'s RL and RLL to RHO_T63L and 71C RHO_T63LL. RL is log(rho/T6**3), set here as RL = DL - 3*(TL-6); the 72C old names clash with RL and RLL as used for radius throughout the rest 73C of the code, which is why main.f declared only a 10-member prefix of 74C this block to avoid importing them. With the clash gone main.f can 75C declare the full block, which a USE statement will require once this 76C becomes a module. 77C 10/13/26 MHP via Claude: removed the unused COMMON/IOUNITS/ declaration 78C (no member referenced here; residue of the obsolete-option removals). 79C 9/25/26 MHP via Claude 80C code changed to remove dead IACAT, IJVS, IJENT, IJDEL unit 81C reservations from COMMON/ACDPTH/ (declared but never opened anywhere) 82C *LCLCD, AGEOUT(5), IACAT, IJLAST, LJLAST, LJWRT, LADON,LAOLY, IJVS, 83C END JVS 84 DIMENSION A(4),B(4),C(4) 85 SAVE 86C ENSURE THAT WE ARE WITHIN THE OPAL 95 TABLES. 87C COMPUTE LOG R = RHO/T6**3 88 RHO_T63L = DL - 3.0D0*(TL-6.0D0) 89C CHECK T 90 IF(TL.LT.TGR(1) .OR. TL.GT.TGR(NUMT))THEN 91 WRITE(*,5)TL 92 WRITE(ISHORT,5)TL 93 5 5 FORMAT(' LOG T OF',F11.6,'OUT OF OPAL 95 TABLE RANGE') 94C 10/30/26 MHP via Claude: bare STOP -> STOP 1, the run-time error code (issue 38). 95 STOP 1 96 ENDIF 97C CHECK TO SEE IF EXTRAPOLATION BELOW THE FIRST TABLE ELEMENT 98C IN DENSITY IS NEEDED. 99 IF(RHO_T63L.LT.R0GR(1))THEN 100 LEXTL = .TRUE. 101 ELSE 102 LEXTL = .FALSE. 103 ENDIF 104C THIS SECTION IS REPLACED BY LINEAR EXTRAPOLATION IN R FROM THE 105C LAST TWO OPACITY VALUES. 106C IF(RHO_T63L.LT.(R0GR(1)-4.0D0) .OR. RHO_T63L.GT.(R0GR(NUMD)+1.0D0))THEN 107C WRITE(*,7)RHO_T63L,TL 108C 7 FORMAT(' LOG R OF',F11.6,' OUT OF OPAL 95 TABLE RANGE;TL' 109C * ,F11.6) 110C STOP 111C ENDIF 112C GET INTERPOLATION FACTORS IN T 113 IF(TL.GE.TGR(JT+2))THEN 114 DO I = JT+3,NUMT-1 115 IF(TL.LT.TGR(I))THEN 116 JT = I - 2 117 GOTO 10 118 ENDIF 119 END DO 120 JT = NUMT - 3 121 10 10 CONTINUE 122 ELSE 123 DO I = JT+1,2,-1 124 IF(TL.GT.TGR(I))THEN 125 JT = I - 1 126 GOTO 20 127 ENDIF 128 END DO 129 JT = 1 130 20 20 CONTINUE 131 ENDIF 132C T INTERPOLATION FACTORS 133 DO I = 1,4 134 A(I) = TGR(JT+I-1) 135 END DO 136C GET INTERPOLATION FACTORS IN T 137 CALL INTERP(A,B,C,TL) 138 DO I = 1,4 139 BT(I) = B(I) 140 BTT(I) = C(I) 141 END DO 142C GET INDICES IN RHO FOR EACH OF THE 4 T VALUES. 143 IF(RHO_T63L.GE.R0GR(JD(1)+2))THEN 144 DO I = JD(1)+3,NUMD-1 145 IF(RHO_T63L.LT.R0GR(I))THEN 146 JD(1) = I - 2 147 GOTO 30 148 ENDIF 149 END DO 150 JD(1) = NUMD - 3 151 30 30 CONTINUE 152 ELSE 153 DO I = JD(1)+1,2,-1 154 IF(RHO_T63L.GT.R0GR(I))THEN 155 JD(1) = I - 1 156 GOTO 40 157 ENDIF 158 END DO 159 JD(1) = 1 160 40 40 CONTINUE 161 ENDIF 162C NOW CHECK IF WITHIN PORTION OF TABLE WITH DATA. 163C NOTE THAT A POINT CAN BE WITHIN THE TABLE AND STILL HAVE THE 164C 4X4 MATRIX OF T,RHO AROUND IT OUTSIDE THE TABLE. 165C CHECK LOW RHO, LOW T REGIME; SINCE THE EMPTY REGION DEPENDS ON Z AT 166C LOW X, POSTPONE UNTIL YOU KNOW IF IT IS A 2D, 3D, OR 4D INTERPOLATION. 167 IF(X.LT.0.2D0 .AND. JT. LE. 5)THEN 168 LCHK = .TRUE. 169 GOTO 60 170 ELSE 171 LCHK = .FALSE. 172 ENDIF 173C CHECK HIGH RHO, HIGH T REGIME(BECAUSE OF THE TABLE GEOMETRY, THE 174C POINT WILL FALL OUTSIDE THE TABLE FIRST AT THE HIGHEST T). 175 LSHIFT = .FALSE. 176 LEXTH = .FALSE. 177 IF(CAPPA(JX(1,1),JT+3,JD(1)+3) .GE. 9.9D0) THEN 178 LSHIFT = .TRUE. 179C NOW THERE ARE TWO POSSIBLE SOLUTIONS: WE CAN USE A DIFFERENT SET 180C OF DENSITIES FOR EACH TEMPERATURE AND STAY WITHIN THE TABLE, OR 181C THE DESIRED DENSITY COULD BE OUTSIDE THE TABLE AT ONE OR MORE 182C TEMPERATURES. IF THE DESIRED DENSITY IS OUTSIDE THE TABLE, 183C LINEAR EXTRAPOLATION IN R AT FIXED T IS USED. 184 DO K = 1,4 185 IF(CAPPA(JX(1,1),JT+K-1,JD(1)+1).GE.9.9D0)THEN 186 LEXTH = .TRUE. 187 LEXH(K) = .TRUE. 188 ELSE 189 LEXH(K) = .FALSE. 190 ENDIF 191 END DO 192C FIND 4 NEAREST ELEMENTS WITHIN THE TABLE AT EACH OF THE 4 TEMPERATURES. 193 IXX = JX(1,1) 194 IDD = JD(1) 195 DO I = 1,4 196 ITT = JT+I-1 197 IF(CAPPA(IXX,ITT,IDD+3).GT.9.9D0)THEN 198 DO J = IDD+2,1,-1 199 IF(CAPPA2(IXX,ITT,J).LE.9.9D0)THEN 200 JD(I) = J - 3 201 IDD = J - 3 202 GOTO 50 203 ENDIF 204 END DO 205 50 50 CONTINUE 206 ELSE 207 JD(I) = IDD 208 ENDIF 209 END DO 210 ENDIF 211 60 60 CONTINUE 212C GET INTERPOLATION FACTORS IN RHO FOR EACH OF THE 4 T VALUES. 213C ASSIGN ALL THE SAME VALUES IF GRID IN R IS SAME FOR ALL T. 214 IF(.NOT.LSHIFT)THEN 215 DO I = 1,4 216 A(I) = R0GR(JD(1)+I-1) 217 END DO 218 IF(.NOT.LEXTL)THEN 219 CALL INTERP(A,B,C,RHO_T63L) 220 ELSE 221 RHO_T63LL = R0GR(1) 222 CALL INTERP(A,B,C,RHO_T63LL) 223 ENDIF 224 DO J = 1,4 225 JD(J) = JD(1) 226 DO I = 1,4 227 BD(J,I) = B(I) 228 BDD(J,I) = C(I) 229 END DO 230 END DO 231 ELSE 232 DO J = 1,4 233 DO I = 1,4 234 A(I) = R0GR(JD(J)+I-1) 235 END DO 236C IF EXTRAPOLATION IS BEING PERFORMED, ASSIGN R=REDGE 237C AND EXTRAPOLATE LINEARLY LATER 238 IF(.NOT.LEXH(J))THEN 239 CALL INTERP(A,B,C,RHO_T63L) 240 ELSE 241 RR = R0GR(JD(J)+3) 242 CALL INTERP(A,B,C,RR) 243 RLH(J) = RR 244 ENDIF 245 DO I = 1,4 246 BD(J,I) = B(I) 247 BDD(J,I) = C(I) 248 END DO 249 END DO 250 ENDIF 251C DETERMINE WHETHER A 2D (RHO,T); 3D (X,RHO,T); OR 4D (Z,X,RHO,T) 252C INTERPOLATION IS NEEDED TO GET THE OPACITY. 253C JVS 04/11 force 4d interpolation for acoustic depth calculations 254 IF (LADON .AND. LACOUT) GOTO 153 ! If we're worrying about the acoustic depth, default to 4d interp 255 IF(ABS(Z-ZTAB)/MAX(ZTAB,1.0D-6).LE.1.0D-4)THEN 256 IF(ABS(X-XOP95).LE.1.0D-4)THEN 257C 2D INTERPOLATION IN SURFACE TABLE 258 CALL OP952D(O,OL,QOD,QOT) 259 GOTO 9999 260 ELSE 261C 3D INTERPOLATION IN FIXED Z TABLE (X,T,RHO) 262C GET INTERPOLATION FACTORS IN X. 263 IF(X.GT.XXG(JX(1,1)+2))THEN 264 DO I = JX(1,1)+3,NUMX-1 265 IF(X.LT.XXG(I))THEN 266 JX(1,1) = I - 2 267 GOTO 70 268 ENDIF 269 END DO 270 JX(1,1) = NUMX - 3 271 70 70 CONTINUE 272 ELSE 273 DO I = JX(1,1)+1,2,-1 274 IF(X.GE.XXG(I))THEN 275 JX(1,1) = I - 1 276 GOTO 80 277 ENDIF 278 END DO 279 JX(1,1) = 1 280 80 80 CONTINUE 281 ENDIF 282C CHECK FOR ABSENT TABLES AT HIGH X. 283 IF(ZTAB.GE.0.04D0)THEN 284 IF(X.GE.0.8D0)THEN 285 IF(ZTAB.LT.0.1D0)THEN 286C NO TABLE 9 (X=0.95) 287 JX(1,4) = NUMX 288 JX(1,3) = NUMX - 2 289 JX(1,2) = NUMX - 3 290 JX(1,1) = NUMX - 4 291 ELSE 292C NO TABLE 9 OR TABLE 10 293 JX(1,4) = NUMX - 2 294 JX(1,3) = NUMX - 3 295 JX(1,2) = NUMX - 4 296 JX(1,1) = NUMX - 5 297 ENDIF 298 ELSE 299 JXX = JX(1,1) - 1 300 DO I = 1,4 301 JX(1,I) = JXX + I 302 END DO 303 ENDIF 304 ELSE 305 JXX = JX(1,1) - 1 306 DO I = 1,4 307 JX(1,I) = JXX + I 308 END DO 309 ENDIF 310C X INTERPOLATION FACTORS 311 DO I = 1,4 312 IXX = JX(1,I) 313 A(I) = XXG(IXX) 314 END DO 315 IF(JX(1,4).EQ.NUMX)THEN 316 A(4) = 1.0D0 - ZTAB 317 ENDIF 318C GET INTERPOLATION FACTORS IN X 319 CALL INTRP2(A,B,X) 320 DO J = 1,4 321 BX(1,J) = B(J) 322 END DO 323 CALL OP953D(O,OL,QOD,QOT) 324 GOTO 9999 325 ENDIF 326 ENDIF 327153 153 CONTINUE 328C 4D INTERPOLATION IN Z,X,T,RHO 329C GET NEAREST TABLES IN Z. 330 IF(Z.GT.ZZG(JZ+2))THEN 331c mhp 7/12 corrected typo 332C DO I = JZ+3,NUMT-1 333 DO I = JZ+3,NUMZ-1 334 IF(Z.LT.ZZG(I))THEN 335 JZ = I - 2 336 GOTO 110 337 ENDIF 338 END DO 339 JZ = NUMZ - 3 340 110 110 CONTINUE 341 ELSE 342 DO I = JZ+1,2,-1 343 IF(Z.GE.ZZG(I))THEN 344 JZ = I - 1 345 GOTO 120 346 ENDIF 347 END DO 348 JZ = 1 349 120 120 CONTINUE 350 ENDIF 351C mhp 7/12 added failsafe 352 JZ = MIN(JZ,NUMZ-3) 353C Z INTERPOLATION FACTORS 354 DO I = 1,4 355 A(I) = ZZG(JZ+I-1) 356 END DO 357C INTERPOLATION FACTORS IN Z 358 CALL INTRP2(A,B,Z) 359 DO I = 1,4 360 BZ(I) = B(I) 361 END DO 362C GET NEAREST TABLES IN X AT EACH VALUE OF Z 363 DO K = 1,4 364 ZZ = ZZG(JZ+K-1) 365 IXX = JX(K,1) 366 IF(X.GT.XXG(IXX+2))THEN 367 DO I = JX(K,1)+3,NUMX-1 368 IF(X.LT.XXG(I))THEN 369 JX(K,1) = I - 2 370 GOTO 75 371 ENDIF 372 END DO 373 JX(K,1) = NUMX - 3 374 75 75 CONTINUE 375 ELSE 376 DO I = JX(K,1)+1,2,-1 377 IF(X.GE.XXG(I))THEN 378 JX(K,1) = I - 1 379 GOTO 85 380 ENDIF 381 END DO 382 JX(K,1) = 1 383 85 85 CONTINUE 384 ENDIF 385C CHECK FOR ABSENT TABLES AT HIGH X. 386 IF(ZZ.GE.0.04D0)THEN 387 IF(X.GE.0.8D0)THEN 388 IF(ZZ.LT.0.1D0)THEN 389C NO TABLE 9 (X=0.95) 390 JX(K,4) = NUMX 391 JX(K,3) = NUMX - 2 392 JX(K,2) = NUMX - 3 393 JX(K,1) = NUMX - 4 394 ELSE 395C NO TABLE 9 OR TABLE 10 396 JX(K,4) = NUMX - 2 397 JX(K,3) = NUMX - 3 398 JX(K,2) = NUMX - 4 399 JX(K,1) = NUMX - 5 400 ENDIF 401 ELSE 402 JXX = JX(K,1) - 1 403 DO I = 1,4 404 JX(K,I) = JXX + I 405 END DO 406 ENDIF 407 ELSE 408 JXX = JX(K,1) - 1 409 DO I = 1,4 410 JX(K,I) = JXX + I 411 END DO 412 ENDIF 413C X INTERPOLATION FACTORS 414 DO I = 1,4 415 IXX = JX(K,I) 416 A(I) = XXG(IXX) 417 END DO 418 IF(JX(K,4).EQ.NUMX)THEN 419 A(4) = 1.0D0 - ZZ 420 ENDIF 421C GET INTERPOLATION FACTORS IN X 422 CALL INTRP2(A,B,X) 423 DO J = 1,4 424 BX(K,J) = B(J) 425 END DO 426 END DO 427 CALL OP954D(O,OL,QOD,QOT) 428 9999 9999 CONTINUE 429 RETURN 430 END

1"""Python version of getopal95.f (YREC development, 7be4133): the driver for 2the OPAL 1995 opacity tables.""" 3 4import numpy as np 5 6from pyrec import acdpth_mod 7from pyrec import iounits_mod 8from pyrec import llot95_mod 9from pyrec import llot95a_mod 10from pyrec import llot95e_mod 11from pyrec import op95interp_mod 12from pyrec.interp import interp 13from pyrec.intrp2 import intrp2 14from pyrec.op952d import op952d 15from pyrec.op953d import op953d 16from pyrec.op954d import op954d 17 18NUMT = 70 19NUMD = 19 20NUMX = 10 21NUMZ = 13 22 23# GETOPAL95 has SAVE: its locals keep their values between calls. Only 24# LSHIFT is read before it is set (a suspected bug, see the page); the work arrays A, B, C are 25# kept here as well, as in the Fortran. 26a = np.zeros(5) 27b = np.zeros(5) 28c = np.zeros(5) 29lshift = False 30lchk = False 31 32 33def _x_stencil(k, x, zz): 34 """The X stencil JX(K,1..4) and its weights BX(K,1..4) at the Z value zz. 35 36 The same block appears twice in the Fortran: for the 3D path (K = 1, 37 with ZZ = ZTAB, labels 70 and 80) and for each Z node of the 4D path 38 (labels 75 and 85). It is one helper here; the statements are those of 39 the Fortran. 40 41 """ 42 43 jx = op95interp_mod.jx 44 xxg = llot95a_mod.xxg 45 46 # Stencil in X (Fortran lines 263–281, 366–384). 47 if x > xxg[jx[k, 1] + 2]: 48 for i in range(jx[k, 1] + 3, NUMX): # I = JX(K,1)+3, NUMX-1 49 if x < xxg[i]: 50 jx[k, 1] = i - 2 51 break 52 else: 53 jx[k, 1] = NUMX - 3 54 else: 55 for i in range(jx[k, 1] + 1, 1, -1): # I = JX(K,1)+1, 2, -1 56 if x >= xxg[i]: 57 jx[k, 1] = i - 1 58 break 59 else: 60 jx[k, 1] = 1 61 62 # Check for absent tables at high X (Fortran lines 283–309, 386–412). 63 if zz >= 0.04 and x >= 0.8: 64 if zz < 0.1: 65 # No table 9 (X = 0.95). 66 jx[k, 4] = NUMX 67 jx[k, 3] = NUMX - 2 68 jx[k, 2] = NUMX - 3 69 jx[k, 1] = NUMX - 4 70 else: 71 # No table 9 or table 10. 72 jx[k, 4] = NUMX - 2 73 jx[k, 3] = NUMX - 3 74 jx[k, 2] = NUMX - 4 75 jx[k, 1] = NUMX - 5 76 else: 77 jxx = jx[k, 1] - 1 78 for i in range(1, 5): 79 jx[k, i] = jxx + i 80 81 # X interpolation factors (Fortran lines 311–322, 414–425). 82 for i in range(1, 5): 83 a[i] = xxg[jx[k, i]] 84 if jx[k, 4] == NUMX: 85 a[4] = 1.0 - zz 86 intrp2(a, b, x) 87 for j in range(1, 5): 88 op95interp_mod.bx[k, j] = b[j] 89 90 91def getopal95(dl: float, tl: float, x: float, z: float, o: float, ol: float, qod: float, qot: float) -> tuple[float, 92 float, float, float]: 93 """Opacity from the OPAL 1995 tables, given rho, T, X and Z. 94 95 Finds the interpolation stencils and weights in T, R = rho/T6**3, X and 96 Z, handles the tables' empty high-density corner, and calls the 2D, 3D 97 or 4D interpolator. 98 99 Translates: getopal95.f#GETOPAL95 (7be4133) 100 101 Parameters 102 ---------- 103 dl, tl : float 104 log10 density and log10 temperature. 105 x, z : float 106 Hydrogen and metal mass fractions. 107 o, ol, qod, qot : float 108 Set here (passed in with any value). 109 110 Returns 111 ------- 112 o, ol, qod, qot 113 The opacity, its log10, and dlog(kappa)/dlog(rho) and 114 dlog(kappa)/dlog(T), as from the interpolators. 115 116 """ 117 118 global lshift, lchk 119 im = op95interp_mod 120 tgr = llot95a_mod.tgr 121 r0gr = llot95a_mod.r0gr 122 zzg = llot95a_mod.zzg 123 cappa = llot95_mod.cappa 124 cappa2 = llot95a_mod.cappa2 125 jd = im.jd 126 jx = im.jx 127 128 # log R = log rho - 3 (log T - 6), and the check of T (Fortran lines 88–96). 129 im.rho_t63l = dl - 3.0 * (tl - 6.0) 130 if tl < tgr[1] or tl > tgr[NUMT]: 131 text = f" LOG T OF{tl:11.6f}OUT OF OPAL 95 TABLE RANGE" # format 5 132 print(text) 133 iounits_mod.write(iounits_mod.ishort, text) 134 raise SystemExit(1) # STOP 1, the run-time error code 135 136 # Extrapolation below the first density of the table? (Fortran lines 99–103.) 137 im.lextl = im.rho_t63l < r0gr[1] 138 139 # Stencil in T, starting from the last one (Fortran lines 113–131). 140 # Suspected bug (see the page): with the OP tables of the recent 141 # mixtures, the stencil reaches empty rows from log T = 7.9 up. 142 if tl >= tgr[im.jt + 2]: 143 for i in range(im.jt + 3, NUMT): 144 if tl < tgr[i]: 145 im.jt = i - 2 146 break # label 10 147 else: 148 im.jt = NUMT - 3 149 else: 150 for i in range(im.jt + 1, 1, -1): 151 if tl > tgr[i]: 152 im.jt = i - 1 153 break # label 20 154 else: 155 im.jt = 1 156 jt = im.jt 157 158 # T interpolation factors (Fortran lines 133–141). 159 for i in range(1, 5): 160 a[i] = tgr[jt + i - 1] 161 interp(a, b, c, tl) 162 for i in range(1, 5): 163 im.bt[i] = b[i] 164 im.btt[i] = c[i] 165 166 # Stencil in R (for the first temperature) (Fortran lines 143–161). 167 if im.rho_t63l >= r0gr[jd[1] + 2]: 168 for i in range(jd[1] + 3, NUMD): 169 if im.rho_t63l < r0gr[i]: 170 jd[1] = i - 2 171 break # label 30 172 else: 173 jd[1] = NUMD - 3 174 else: 175 for i in range(jd[1] + 1, 1, -1): 176 if im.rho_t63l > r0gr[i]: 177 jd[1] = i - 1 178 break # label 40 179 else: 180 jd[1] = 1 181 182 # The empty corner (Fortran lines 167–211). At low X and low T the check is 183 # postponed (and never made): GOTO 60 skips the reset of LSHIFT and LEXTH. 184 if x < 0.2 and jt <= 5: 185 # Suspected bug (see the page): LSHIFT and LEXTH keep their values 186 # from the last call on this path; LCHK is never used. 187 lchk = True 188 else: 189 lchk = False 190 lshift = False 191 im.lexth = False 192 if cappa[jx[1, 1], jt + 3, jd[1] + 3] >= 9.9: 193 lshift = True 194 for k in range(1, 5): 195 if cappa[jx[1, 1], jt + k - 1, jd[1] + 1] >= 9.9: 196 im.lexth = True 197 im.lexh[k] = True 198 else: 199 im.lexh[k] = False 200 # The four nearest densities with data at each temperature. 201 ixx = jx[1, 1] 202 idd = jd[1] 203 for i in range(1, 5): 204 itt = jt + i - 1 205 if cappa[ixx, itt, idd + 3] > 9.9: 206 for j in range(idd + 2, 0, -1): 207 # Suspected bug (see the page): CAPPA2 here, CAPPA above. 208 if cappa2[ixx, itt, j] <= 9.9: 209 jd[i] = j - 3 210 idd = j - 3 211 break # label 50 212 else: 213 jd[i] = idd 214 # label 60 215 216 # Weights in R for each of the four temperatures (Fortran lines 214–250). 217 if not lshift: 218 for i in range(1, 5): 219 a[i] = r0gr[jd[1] + i - 1] 220 if not im.lextl: 221 interp(a, b, c, im.rho_t63l) 222 else: 223 im.rho_t63ll = r0gr[1] 224 interp(a, b, c, im.rho_t63ll) 225 for j in range(1, 5): 226 jd[j] = jd[1] 227 for i in range(1, 5): 228 im.bd[j, i] = b[i] 229 im.bdd[j, i] = c[i] 230 else: 231 for j in range(1, 5): 232 for i in range(1, 5): 233 a[i] = r0gr[jd[j] + i - 1] 234 # Outside the data: weights at the edge, linear extrapolation 235 # later. 236 if not im.lexh[j]: 237 interp(a, b, c, im.rho_t63l) 238 else: 239 rr = r0gr[jd[j] + 3] 240 interp(a, b, c, rr) 241 im.rlh[j] = rr 242 for i in range(1, 5): 243 im.bd[j, i] = b[i] 244 im.bdd[j, i] = c[i] 245 246 # 2D, 3D or 4D? The 4D path is forced for acoustic-depth output (JVS 04/11) 247 # (Fortran lines 254–326). 248 if not (acdpth_mod.ladon and acdpth_mod.lacout): 249 ztab = llot95_mod.ztab 250 if abs(z - ztab) / max(ztab, 1.0e-6) <= 1.0e-4: 251 if abs(x - llot95e_mod.xop95) <= 1.0e-4: 252 # 2D in the surface table. 253 o, ol, qod, qot = op952d(o, ol, qod, qot) 254 return o, ol, qod, qot 255 # 3D in the table at fixed Z (X, T, rho). 256 _x_stencil(1, x, ztab) 257 o, ol, qod, qot = op953d(o, ol, qod, qot) 258 return o, ol, qod, qot 259 260 # label 153: 4D in Z, X, T, rho. Nearest tables in Z (Fortran lines 330–361). 261 if z > zzg[im.jz + 2]: 262 for i in range(im.jz + 3, NUMZ): # up to NUMZ, a typo fixed (mhp 7/12) 263 if z < zzg[i]: 264 im.jz = i - 2 265 break # label 110 266 else: 267 im.jz = NUMZ - 3 268 else: 269 for i in range(im.jz + 1, 1, -1): 270 if z >= zzg[i]: 271 im.jz = i - 1 272 break # label 120 273 else: 274 im.jz = 1 275 im.jz = min(im.jz, NUMZ - 3) # failsafe: four Z points (mhp 7/12) 276 for i in range(1, 5): 277 a[i] = zzg[im.jz + i - 1] 278 intrp2(a, b, z) 279 for i in range(1, 5): 280 im.bz[i] = b[i] 281 282 # Nearest tables in X at each value of Z, then the 4D interpolation 283 # (Fortran lines 363–427). 284 for k in range(1, 5): 285 _x_stencil(k, x, zzg[im.jz + k - 1]) 286 o, ol, qod, qot = op954d(o, ol, qod, qot) 287 return o, ol, qod, qot

getopal95.f

Purpose

How it fits

The method

Suspected bugs, checked by running the code

1. (not reached) LSHIFT and LEXTH left over from the last call.

2. (changes results, large) Checked by KC (Nov 10, 2026) Empty temperatures in the OP tables above $\log T = 8.0$.

3. (never occurred) The wrong array in the empty-corner search.

Definitions

subroutineGETOPAL95

Change log

1. (not reached) `LSHIFT` and `LEXTH` left over from the last call.