ll95tbl.f

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

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

Variables on this page (13)

cappa2, iliv95, ishort, numd, numx, numxz, numz, nz, r0gr, tgr, xxg, zopal951, zzg

Purpose

LL95TBL reads the full set of OPAL 1995 opacity tables from a file into memory: $\log\kappa$ as a function of $Z$, $X$, $\log T$ and $\log R$ with $R = \rho/T_6^3$. Then it builds the table at the run's metallicity with OP95ZTAB.

How it fits

Who calls it: only SETUPOPAC, at set-up, when the run uses the OPAL95 tables.
Input: the file name FLIV95 (argument) and its unit ILIV95 (IOUNITS_MOD, formerly COMMON LOPAL95).
Output: LLOT95A_MOD (the grids TGR, R0GR read here, and XXG, ZZG, NUMXM, NZ, initialized in the module, which formerly were DATA here; the tables CAPPA2(126,70,19)), read by GETOPAL95, OP95XTAB, OP95ZTAB and OP954D; and OP95INTERP_MOD (formerly COMMON OP95INDX), whose starting interpolation indices JZ, JX, JT, JD (all 1) are now initialized in the module.
Calls: OP95ZTAB with $Z$ = ZOPAL951 (from NMLPAR_MOD, formerly COMMON NEWOPAC). Its error messages go to ISHORT (IOUNITS_MOD).

The method

Skip the file's header until the first line starting with "TABLE"; read its $X$ and $Z$.
For each table: check that its $X$ and $Z$ are the expected ones (from the grids XXG, 10 values from 0 to 1, and ZZG, 13 values from 0 to 0.1; the last $X$ is $1 - Z$), else stop (exit status 3, the code for input errors); read the 19 values of $\log R$; read 70 rows in $\log T$. The last 13 rows have fewer densities (the corner of the table without data), and the missing entries are set to 9.999.

Question 1 (lines 136–171 ↗). Hard-coded table shape. The number of complete rows (57), the row ranges with 18, 17, 16, 15, 14 densities, and the filler 9.999 are written into the code.
The file is ordered by $X$ and, within each $X$, by $Z$; fewer $Z$ values exist for the high $X$ (10 for $X = 0.95$, 12 for $X = 1-Z$). Table $(Z_i, X_j)$ is stored at index $N = \texttt{NZ}(i) + j$ (line 117 ↗; for the tables at $X = 1 - Z$ this overwrites two others, bug 1). It stops after 126 tables. An end of file before then, in the header, in a table's $\log R$ line or in its rows, stops the run (exit status 3) with a message to the terminal and to ISHORT naming the file and the table being read (NN of NUMXZ; NN = 0 while the header is skipped). An end of file at the line that gives the next table's $X$ and $Z$ (format 900, which has no END=) is left to the compiler's run-time error.
Build the fixed-$Z$ table: OP95ZTAB(ZOPAL951).

From YREC's documentation

Control Options: FLIV95 is the path of the atomic opacity table; ZOPAL951 makes it build a table for the surface mixture at this $Z$, which speeds up models without diffusion.

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. (changes results, large) Checked by KC (Nov 10, 2026) The $X = 1 - Z$ tables at $Z$ = 0.06 and 0.08 overwrite the $X = 0$ tables at $Z$ = 0.08 and 0.10.

(line 117 ↗, lines 176–191 ↗; llot95a_mod.f, lines 46–47 ↗) NZ gives $Z$ = 0.06 and 0.08 nine slots each (NZ = 100, 109, 118 for the last three $Z$, as NUMXM = 9, 9, 8 tables), but the table at $X = 1 - Z$ is stored at NZ(IZ)+10. For $Z$ = 0.06 that is slot 110, the slot of the $X = 0$ table at $Z$ = 0.08; for $Z$ = 0.08 it is slot 119, the $X = 0$ table at $Z$ = 0.10. The file is ordered by $X$ (in GS98.OP17, tables 125 and 126 are $X$ = 0.94, $Z$ = 0.06 and $X$ = 0.92, $Z$ = 0.08), so these two are read last and overwrite the $X = 0$ tables, while slots 109 and 118 stay empty. The readers use the same NZ and the index JX, so the $X = 1 - Z$ tables are found where they were stored, but a reader asking for $X = 0$ at $Z$ = 0.08 or 0.10 gets the opacity of a gas with $X \approx 0.93$: OP954D (op954d.f, lines 45–48 ↗) whenever the $Z$ stencil reaches 0.08 ($Z \ge 0.04$) and the $X$ stencil starts at $X = 0$ ($X < 0.2$), as in a helium-burning core, and always in GETOPAC's helium-burning branch, which calls GETOPAL95 at $Z = 0.1$ (getopac.f, lines 117–120 ↗): There, for $X = 0$, the Lagrange weights in $Z$ and $X$ put all the weight on slot 119, so the interior opacity used for a metal-rich helium-burning core is that of $X = 0.92$, $Z = 0.08$; and OP95ZTAB (op95ztab.f, lines 54–67 ↗) when ZOPAL951 > 0.04. Taken with every table file in this format (they all hold the same 126 compositions), in every helium-burning core whose $Z$ reaches 0.04.

A possible fix, in LLOT95A_MOD, giving every $Z$ ten slots (LL95TBL keeps its own NUMXZ = 126, the number of tables in the file):

llot95a_mod.f (another file)

@@ lines 34–40 @@
       INTEGER*4, PARAMETER, PRIVATE :: NUMD = 19
       INTEGER*4, PARAMETER, PRIVATE :: NUMT = 70
       INTEGER*4, PARAMETER, PRIVATE :: NUMX = 10
-      INTEGER*4, PARAMETER, PRIVATE :: NUMXZ = 126
+C 11/09/26 KC via Claude
+C code changed to give every Z ten table slots, so that the X=1-Z
+C tables at Z = 0.06 and 0.08 (slot NZ+10) no longer overwrite the
+C X=0 tables at Z = 0.08 and 0.10.
+      INTEGER*4, PARAMETER, PRIVATE :: NUMXZ = 130
       INTEGER*4, PARAMETER, PRIVATE :: NUMZ = 13

       REAL*8    :: TGR(NUMT)
@@ lines 44–49 @@
      *  0.08D0, 0.10D0/)
       REAL*8    :: CAPPA2(NUMXZ,NUMT,NUMD)
       INTEGER*4 :: NUMXM(NUMZ) = (/10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 9, 9, 8/)
-      INTEGER*4 :: NZ(NUMZ) = (/0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 109, 118/)
+      INTEGER*4 :: NZ(NUMZ) = (/0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110,
+     *  120/)

       END MODULE LLOT95A_MOD

How the code itself shows it. The reader stores each table by its position in the $X$ grid, not by its count: N = NZ(IZ)+IX (line 117 ↗), where IX = 10 is the table at $X = 1 - Z$ (lines 118–122 ↗). The loop at lines 177–190 ↗ skips IX = 9 ($X$ = 0.95) for $Z$ = 0.06, 0.08, and 0.10 (NMAX = NUMZ - 3) but still reads IX = 10 for $Z$ = 0.06 and 0.08 (NMAX = NUMZ - 1), so each of these two $Z$ holds nine tables spread over ten positions. OP95ZTAB reads the slots by $X$ position for every $Z$ as well: NZ(IDZ)+I for the tables 1–8 (op95ztab.f, lines 54–59 ↗), NZ(IDZ2)+9 for $X$ = 0.95 (op95ztab.f, lines 85–88 ↗), and NZ(IDZ2)+10 for $X = 1 - Z$ (op95ztab.f, lines 114–117 ↗); so does OP954D, with IZ = NZ(JZ+K-1) and IXZ = JX(K,J)+IZ (op954d.f, lines 45–48 ↗). NZ's steps of 9 after $Z$ = 0.06 and 0.08 (100, 109, and 118; llot95a_mod.f, line 47 ↗) follow NUMXM, the number of tables at each $Z$ (llot95a_mod.f, line 46 ↗), which nothing else reads, while every index into CAPPA2 counts positions. The fix makes the steps 10 throughout, as for the first ten $Z$.

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): 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 (each one of YREC's cases with a few settings changed), of which one stopped in its first second. In all 33 others, the tables at $X$ = 0.94, $Z$ = 0.06 and $X$ = 0.92, $Z$ = 0.08 were stored over the $X = 0$ slots at the start. The 4-D interpolation then took one of them as the $X = 0$ table in six runs, all with helium-burning or giant cores, with the full weight in $X$ on it: 85 million times in YREC's examples/run_from_zahb_to_tahb/ZAHB_m1p0feh+0p0_TAHB run with the AAG21 table of GETOPAL95's bug 2, 64 million in examples/run_from_zahb_to_tahb/ZAHB_m0p8feh+0p0_TAHB, 31 million in testsuite/Test_m9p0feh+0p0GS98_base_TAHB, 29 million in testsuite/Test_m0p95feh+0p0_GN93_TAHB, and fewer in two more. With the fix, testsuite/Test_m0p95feh+0p0_GN93_TAHB changes by up to 0.36 dex in $\log L$, 0.083 dex in $\log T_c$, and 0.43 in central $Y$, against $1.7\times10^{-7}$ dex, $3.8\times10^{-8}$ dex, and $8\times10^{-9}$ for the null change used as a noise floor (nearly every opacity, all with $3 \le \log T \le 9$ and $-6 \le \log\rho \le 9.75$, times $1 + 10^{-10}$): real and large; the track ends at 0.0943 Gyr instead of 0.1293 (854 models instead of 966), both runs ending normally. On the 3 $M_\odot$ giant testsuite/Test_3Msolar_GS_rot_gray_OPALSCV_SF2_GS98_OP_CF10_RGp2diffrotI3cz2, the change is within the noise. Intent: none found.

Definitions

subroutineLL95TBL

lines 6–219 · GitHub

Reads the OPAL95 opacity file into LLOT95A_MOD and builds the fixed-$Z$ table (above).

Uses modules IOUNITS_MOD, LLOT95A_MOD, NMLPAR_MOD, OP95INTERP_MOD

Called by SETUPOPAC

Flowchart of its calls, in order (1)

OP95ZTAB 206

Compiler warnings (2)

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

line 123: Inequality comparison for REAL(8) [-Wcompare-reals]
line 123: Inequality comparison for REAL(8) [-Wcompare-reals]

Change log

The changes since rp_1

Changed since rp_1. Declarations only: COMMON blocks to modules, IMPLICIT NONE; the grids' and indices' DATA moved into LLOT95A_MOD and OP95INTERP_MOD; NEWOPAC dropped, ZOPAL951 now from NMLPAR_MOD, which settles former question 1 (the authors' change notes confirm the layout mismatch).

Changed since a8a5950. An end of file no longer ends the read silently: every END= now goes to a message (the file, and the table reached of 126) and a stop, where before it went to a bare RETURN, leaving CAPPA2 partly filled, the file open, and the fixed-$Z$ table unbuilt; an end of file at a table's $\log R$ line, which went to the normal exit, does the same (this settles former question 2). That stop and the stop on a table of the wrong $X$ or $Z$ exit with status 3, the code for input errors (KNOWN_ISSUES.md, issue 38).

1C 2C 3CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 4C LL95TBL 5CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 6 SUBROUTINE LL95TBL(FLIV95) 7C 10/14/26 MHP via Claude 8C code changed to convert COMMON/OP95INTERP/ into OP95INTERP_MOD. 9 USE OP95INTERP_MOD 10C 10/14/26 MHP via Claude 11C code changed to convert COMMON/LLOT95A/ into LLOT95A_MOD. 12 USE LLOT95A_MOD 13C 10/14/26 MHP via Claude 14C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD. 15 USE NMLPAR_MOD 16C 10/14/26 MHP via Claude 17C code changed to convert COMMON/IOUNITS/ into IOUNITS_MOD. 18 USE IOUNITS_MOD 19C 10/20/26 MHP via Claude 20C code changed to replace implicit typing with IMPLICIT NONE and 21C explicit declarations. 22 IMPLICIT NONE 23 REAL*8 XX,XXT,Z,ZZ 24 INTEGER I,IX,IZ,J,N,NFMT,NMAX,NN,NUMD,NUMX,NUMXZ,NUMZ 25c MHP 7/98 MODIFIED TO READ IN ALL METAL ABUNDACES FOR OPACITY TABLES. 26 PARAMETER (NUMD=19) 27 PARAMETER (NUMX=10) 28 PARAMETER (NUMZ=13) 29 PARAMETER (NUMXZ=126) 30 CHARACTER*256 FLIV95 31C MHP 8/25 Removed character file names from common block 32C 9/25/26 MHP via Claude 33C code changed to consolidate all I/O logical unit numbers, previously 34C scattered across COMMON/LUOUT/, COMMON/LUNUM/, COMMON/ALEX06/, 35C COMMON/ACDPTH/, COMMON/NEWOPAC/, COMMON/LOPAL95/, COMMON/ATMOS2/, 36C COMMON/ALATM03/, COMMON/OPALEOS/, and COMMON/SCV2/, into a single 37C COMMON/IOUNITS/ 38 SAVE 39C FULL SET OF TABLES: OPACITY AS A FUNCTION OF Z AND X, T, RHO/T6**3 40C TABLES ARE INCREMENTED IN SETS OF NZ*NX. SO THE TABLES FOR THE 41C THIRD METAL ABUNDANCE (3 X 10**-4)BEGIN AT TABLE 21 AND END AT TABLE 30. 42C FOR THE HIGH VALUES OF Z, THE NUMBER OF X TABLES IS NOT THE SAME (I.E. 43C X<0.9 FOR Z=0.1). 44C FOR EACH COMPOSITION A FULL GRID IN (T,RHO/T6**3) IS RETAINED. 45C INDICES FOR INTERPOLATION IN Z,X,T, AND R; INTERPOLATION FACTORS FOR 46C Z,X,T, AND R (AND DERIVATIVE FACTORS FOR T AND RHO); AND DATA FOR 47C LINEAR EXTRAPOLATION WHEN OUTSIDE THE TABLE IN RHO. 48C 10/9/26 MHP via Claude 49C code changed to merge COMMON/OP95INDX/ into COMMON/OP95INTERP/ (see 50C getopal95.f for the full rationale). This file declares the FULL 51C block, not just the JZ..JD prefix it originally needed: it has its 52C own DATA statement initializing JZ,JX,JT,JD below, and mixing a 53C DATA-initialized short declaration of a COMMON block with a larger 54C declaration of the same name elsewhere is a genuine linker/ABI 55C hazard (confirmed: this caused NaN corruption when first tried as a 56C prefix -- see the linker's "size of symbol changed" warning). 57C Declaring the full block here avoids any size mismatch. 58C 10/13/26 MHP via Claude 59C code changed to rename COMMON/OP95INTERP/'s RL and RLL to RHO_T63L and 60C RHO_T63LL. RL is log(rho/T6**3), set here as RL = DL - 3*(TL-6); the 61C old names clash with RL and RLL as used for radius throughout the rest 62C of the code, which is why main.f declared only a 10-member prefix of 63C this block to avoid importing them. With the clash gone main.f can 64C declare the full block, which a USE statement will require once this 65C becomes a module. 66C 9/23/26 MHP via Claude 67C code changed to remove the OPAL92 opacity tables 68C 9/23/26 MHP via Claude 69C code changed to remove the Kurucz 1990 opacity tables 70C 9/23/26 MHP via Claude 71C code changed to remove the LAOL89 opacity tables 72C 9/24/26 MHP via Claude 73C code changed to remove the Alexander 1995 opacity tables 74C 9/24/26 MHP via Claude 75C code changed to remove dead opacity-table Z-value inputs 76C (this also fixes a pre-existing COMMON /NEWOPAC/ layout mismatch: 77C this file previously omitted TMOLMIN,TMOLMAX,LALEX06, which every 78C other file sharing this block declares) 79C 9/24/26 MHP via Claude 80C code changed to merge COMMON /MISCOPAC/ and COMMON/NWLAOL/ into 81C COMMON /NEWOPAC/ 82C 10/12/26 MHP via Claude 83C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 84C namelist variables that the physics routines read from 27 COMMON 85C blocks which also held derived or working members. Those blocks keep 86C only their non-namelist members; five of them -- NEWMX, VMULT, 87C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 88C Derivation-only namelist variables from the same blocks went to 89C COMMON/NMLBYPASS/ instead, not here. 90C 91C Declared identically in all 121 files that hold any of these members. 92C Generated from one member list, so every copy is byte-identical; the 93C linker's 'size of symbol changed' warning is the check that they are. 94C 10/13/26 MHP via Claude: removed the unused COMMON/NEWOPAC/ declaration 95C (no member referenced here; residue of the obsolete-option removals). 96C NUMBER OF COMPOSITION TABLES AT LOWER Z FOR EACH ABUNDANCE 97C NUMBER OF X TABLES AT EACH Z 98C TABULTED SET OF X 99C TABULATED SET OF Z 100c INDICES FOR INTERPOLATION 101C LOCAL VECTOR, USED TO SKIP LONG HEADER. 102 CHARACTER*132 DUMMY 103C READ IN OPACITY TABLES, SKIPPING HEADER. HEY, IT WORKS. 104 OPEN(ILIV95,FILE=FLIV95,STATUS='OLD',ACCESS='SEQUENTIAL') 105 NFMT=1 106 NN = 0 107 108 10 10 READ(ILIV95,1,END=9999)DUMMY 109 IF( DUMMY(NFMT:NFMT+4).NE.'TABLE' ) GO TO 10 110 1 1 FORMAT(A) 111 READ(DUMMY,'(36X,F7.4,11X,F7.4)')XX,ZZ 112 IZ = 1 113 IX = 1 114 NN = 1 115C ENTRY POINT FOR GETTING NEW TABLES. 116 15 15 CONTINUE 117 N = NZ(IZ)+IX 118 IF(IX.LT.NUMX)THEN 119 XXT = XXG(IX) 120 ELSE 121 XXT = 1.0D0 - ZZG(IZ) 122 ENDIF 123 IF(ZZG(IZ).NE.ZZ .OR. XXT.NE.XX)THEN 124 WRITE(ISHORT,*)' OPAL95: Z ERROR INCOMPATIBLE TABLE' 125C 10/30/26 MHP via Claude: bare STOP -> STOP 3, the input-error code (issue 38). 126 STOP 3 127 ENDIF 128 129C READ IN HEADER INFO: GRID IN RHO/T6**3 130C 10/30/26 MHP via Claude: END=1000 -> END=9999. An end of file at a 131C table header is an error; 1000 is the exit after all NUMXZ tables. 132 READ(ILIV95,20,END=9999) (R0GR(I),I=1,NUMD) 133 20 20 FORMAT(///,4X,19F7.1,/) 134C READ IN FULL TABLE: LOG CAPPA AS A FUNCTION OF LOG T AND 135C LOG R = RHO/T6**3 136 DO I = 1,57 137 READ(ILIV95,30,END=9999)TGR(I),(CAPPA2(N,I,J),J=1,NUMD) 138 END DO 139 30 30 FORMAT(F4.2,19F7.3) 140C MHP 12/97 NOW TREAT CORNER WITHOUT DATA. 141 READ(ILIV95,31,END=9999)TGR(58),(CAPPA2(N,58,J),J=1,NUMD-1) 142 31 31 FORMAT(F4.2,18F7.3) 143 144 CAPPA2(N,58,19) = 9.999D0 145 DO I = 59,60 146 READ(ILIV95,32,END=9999)TGR(I),(CAPPA2(N,I,J),J=1,NUMD-2) 147 32 32 FORMAT(F4.2,17F7.3) 148 DO J = 18,19 149 CAPPA2(N,I,J) = 9.999D0 150 END DO 151 END DO 152 DO I = 61,64 153 READ(ILIV95,33,END=9999)TGR(I),(CAPPA2(N,I,J),J=1,NUMD-3) 154 33 33 FORMAT(F4.2,16F7.3) 155 DO J = 17,19 156 CAPPA2(N,I,J) = 9.999D0 157 END DO 158 END DO 159 DO I = 65,69 160 READ(ILIV95,34,END=9999)TGR(I),(CAPPA2(N,I,J),J=1,NUMD-4) 161 34 34 FORMAT(F4.2,15F7.3) 162 DO J = 16,19 163 CAPPA2(N,I,J) = 9.999D0 164 END DO 165 END DO 166 I = 70 167 READ(ILIV95,35,END=9999)TGR(I),(CAPPA2(N,I,J),J=1,NUMD-5) 168 35 35 FORMAT(F4.2,14F7.3) 169 DO J = 15,19 170 CAPPA2(N,I,J) = 9.999D0 171 END DO 172 173C EXIT IF CORRECT NUMBER OF TABLES READ IN. 174 IF(NN.GE.NUMXZ)GOTO 1000 175 NN = NN + 1 176C NEED TO ACCOUNT FOR FEWER X VALUES AT HIGHER Z. 177 IF(IX.LE.8)THEN 178 NMAX = NUMZ 179 ELSE IF(IX.EQ.9)THEN 180 NMAX = NUMZ - 3 181 ELSE 182 NMAX = NUMZ - 1 183 ENDIF 184C READ IN NEXT METAL ABUNDANCE (AT FIXED X) OR READ IN FIRST METAL 185C ABUNDANCE AT NEXT X. 186 IF(IZ.LT.NMAX)THEN 187 IZ = IZ + 1 188 ELSE 189 IZ = 1 190 IX = IX + 1 191 ENDIF 192 193C RETURN TO READ IN NEXT TABLE. 194 READ(ILIV95,900)XX,ZZ 195 900 900 FORMAT(/36X,F7.4,11X,F7.4) 196 GOTO 15 197 1000 1000 CONTINUE 198 199 CLOSE(ILIV95) 200C NOW GENERATE A TABLE AT A FIXED VALUE OF Z. 201C NOTE THAT FOR METAL DIFFUSION A 4-D INTERPOLATION (IN X,Z,T,RHO) 202C IS PERFORMED RATHER THAN A LINEAR INTERPOLATION BETWEEN TWO FIXED 203C Z TABLES. 204 Z = ZOPAL951 205 206 CALL OP95ZTAB(Z) 207 RETURN 208C 10/30/26 MHP via Claude 209C code changed to stop on a truncated table file. Every END= exit came 210C here, to a bare RETURN, so a short file left CAPPA2 partly filled, the 211C unit open and the fixed-Z table never built, with no message. 212 9999 9999 CONTINUE 213 WRITE(*,9000)TRIM(FLIV95),NN,NUMXZ 214 WRITE(ISHORT,9000)TRIM(FLIV95),NN,NUMXZ 215 9000 9000 FORMAT(' OPAL95: UNEXPECTED END OF FILE IN ',A/ 216 * * ' WHILE READING TABLE',I4,' OF',I4,'. RUN STOPPED') 217C 10/30/26 MHP via Claude: STOP 3, the input-error code (issue 38). 218 STOP 3 219 END

1"""Python version of ll95tbl.f (YREC 7be4133): read the OPAL 1995 opacity 2tables for all metallicities and build the table at the run's Z.""" 3 4from pyrec import iounits_mod # unit numbers; Python opens files by name 5from pyrec import llot95a_mod 6from pyrec import nmlpar_mod 7from pyrec.op95ztab import op95ztab 8 9NUMD = 19 # densities, log R = log(rho/T6**3) 10NUMX = 10 # hydrogen abundances 11NUMZ = 13 # metallicities 12NUMXZ = 126 # tables in the file 13 14 15class _EndOfFile(Exception): 16 """An END= branch of a Fortran READ.""" 17 18 19def _record(f): 20 """One record, as a READ consumes it; _EndOfFile at the end.""" 21 22 line = f.readline() 23 if line == '': 24 raise _EndOfFile 25 return line.rstrip('\n') 26 27 28def _fields(line, widths, start=0): 29 """Split a fixed-format record into floats (a blank field reads as 0, 30 as in Fortran).""" 31 32 out, pos = [], start 33 for w in widths: 34 text = line[pos:pos + w].strip() 35 out.append(float(text.replace('D', 'E')) if text else 0.0) 36 pos += w 37 return out 38 39 40def _xz(line): 41 """X and Z from a table's title line, format (36X,F7.4,11X,F7.4): the 42 11 columns between them ("Y=... Z=") are skipped, not read.""" 43 44 return _fields(line, [7], start=36)[0], _fields(line, [7], start=54)[0] 45 46 47def _read_row(f, n, i, ncol): 48 """A row in log T: F4.2 for TGR(I), then NCOL values F7.3 into 49 CAPPA2(N,I,1..NCOL) (formats 30 to 35).""" 50 51 a = llot95a_mod 52 vals = _fields(_record(f), [4] + [7] * ncol) 53 a.tgr[i] = vals[0] 54 a.cappa2[n, i, 1:ncol + 1] = vals[1:] 55 56 57def _stop_eof(fliv95, nn): 58 """Label 9999: an unexpected end of file; format 9000 to the terminal 59 and the .short file, then STOP 3, the input-error code.""" 60 61 text = (f" OPAL95: UNEXPECTED END OF FILE IN {fliv95.strip()}\n" 62 f" WHILE READING TABLE{nn:4d} OF{NUMXZ:4d}. RUN STOPPED") 63 print(text) 64 iounits_mod.write(iounits_mod.ishort, text) 65 raise SystemExit(3) 66 67 68def ll95tbl(fliv95: str) -> None: 69 """Read the full set of OPAL95 opacity tables (log kappa against Z, X, 70 log T and log R) into llot95a_mod, checking each table's X and Z 71 against the grids, filling the corner without data with 9.999; then 72 build the fixed-Z table with OP95ZTAB at Z = ZOPAL951. 73 74 Translates: ll95tbl.f#LL95TBL (7be4133) 75 76 Arguments: fliv95, the file name. Nothing is returned. 77 78 """ 79 80 a = llot95a_mod 81 82 # Read in the opacity tables, skipping the header. "Hey, it works." 83 # (Fortran: OPEN(ILIV95,...,STATUS='OLD').) Every END= goes to label 84 # 9999, which stops the run (Fortran lines 104–111). 85 with open(fliv95) as f: 86 nfmt = 1 87 nn = 0 88 try: 89 while True: # label 10 90 dummy = _record(f) 91 if dummy[nfmt - 1:nfmt + 4] == 'TABLE': 92 break 93 except _EndOfFile: 94 _stop_eof(fliv95, nn) # label 9999 95 xx, zz = _xz(dummy) 96 iz = 1 97 ix = 1 98 nn = 1 99 try: 100 while True: 101 # Entry point for getting new tables (label 15): the slot, 102 # and the check of X and Z (Fortran lines 116–127). 103 # Suspected bug (see the page): the X = 1-Z tables at Z = 0.06 104 # and 0.08 overwrite the X = 0 tables at Z = 0.08 and 0.10. 105 n = a.nz[iz] + ix 106 if ix < NUMX: 107 xxt = a.xxg[ix] 108 else: 109 xxt = 1.0 - a.zzg[iz] 110 if a.zzg[iz] != zz or xxt != xx: 111 iounits_mod.write(iounits_mod.ishort, 112 ' OPAL95: Z ERROR INCOMPATIBLE TABLE') 113 raise SystemExit(3) # STOP 3, the input-error code 114 115 # Header: the grid in rho/T6**3 (format 20: ///,4X,19F7.1,/) 116 # (Fortran lines 132–133). 117 # An end of file here goes to label 9999 (the outer except). 118 for _ in range(3): 119 _record(f) 120 a.r0gr[1:NUMD + 1] = _fields(_record(f), [7] * NUMD, 121 start=4) 122 _record(f) 123 124 # Full table: log kappa against log T and log R (lines 125 # 136–171). 126 # Question (see the page): hard-coded table shape. 127 for i in range(1, 58): 128 _read_row(f, n, i, NUMD) # format 30 129 # The corner of the table without data (MHP 12/97). 130 _read_row(f, n, 58, NUMD - 1) # format 31 131 a.cappa2[n, 58, 19] = 9.999 132 for i in range(59, 61): 133 _read_row(f, n, i, NUMD - 2) # format 32 134 a.cappa2[n, i, 18:20] = 9.999 135 for i in range(61, 65): 136 _read_row(f, n, i, NUMD - 3) # format 33 137 a.cappa2[n, i, 17:20] = 9.999 138 for i in range(65, 70): 139 _read_row(f, n, i, NUMD - 4) # format 34 140 a.cappa2[n, i, 16:20] = 9.999 141 i = 70 142 _read_row(f, n, i, NUMD - 5) # format 35 143 a.cappa2[n, i, 15:20] = 9.999 144 145 # Exit if the correct number of tables was read in (lines 146 # 174–175). 147 if nn >= NUMXZ: 148 break # go to 1000 149 nn = nn + 1 150 151 # Fewer Z values at the highest X; the next metallicity at 152 # fixed X, or the first at the next X (Fortran lines 177–191). 153 if ix <= 8: 154 nmax = NUMZ 155 elif ix == 9: 156 nmax = NUMZ - 3 157 else: 158 nmax = NUMZ - 1 159 if iz < nmax: 160 iz = iz + 1 161 else: 162 iz = 1 163 ix = ix + 1 164 165 # The next table's X and Z (format 900: /36X,F7.4,11X,F7.4; 166 # no END=, so an end of file here is a Fortran run-time 167 # error, not label 9999) (Fortran lines 194–196). 168 try: 169 _record(f) 170 line = _record(f) 171 except _EndOfFile: 172 raise EOFError('LL95TBL: end of file in READ 900 ' 173 '(Fortran run-time error)') from None 174 xx, zz = _xz(line) 175 except _EndOfFile: 176 _stop_eof(fliv95, nn) # label 9999 177 178 # Label 1000 (Fortran lines 197–206): CLOSE(ILIV95), then a table at a fixed value of Z. (For 179 # metal diffusion a 4-D interpolation in X, Z, T, rho is done instead.) 180 z = nmlpar_mod.zopal951 181 op95ztab(z)