readcoeos06.f

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

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

Variables on this page (19)

dfs, dfsr, dfsx, eos, index, iopale, ishort, midx, mv, mx, nr, nt, nta, rho, t6a, t6list, xx, xz, zz

Purpose

Reads the OPAL 2006 equation-of-state tables (five hydrogen abundances) into memory. It follows READCOEOS01 (removed since rp_1, with the OPAL 2001 equation of state) closely; see also READCO (removed since rp_1, with the OPAL 1995 equation of state).

How it fits

Who calls it: SETEOS06, once at set-up (from SETUPS).
Input: the file on unit IOPALE (from IOUNITS_MOD), opened by PARMIN; closed here. The close (2) before the reading (line 95 ↗), as in READCOEOS01, closes a unit that no part of YREC opens (no unit number in IOUNITS_MOD or PARMIN is 2), so it does nothing (unused?).
Outputs, through modules: AEOS06_MOD (the tables XZ(5, 10, 197, 169), grids, inverse spacings; shared with ESAC06, RHOOFP06, T6RINTEOS06, EQBOUND06; its MIDX, formerly M, is the loop index over the compositions here) and EOS06MISC_MOD (ZZ read here; INDEX, NTA used; EOS(1..10) set to 1 on the first call). The per-composition headers and the extra columns AMU_M, ALOGNE, formerly /EEEEOS06/, are saved local arrays; AMU_M and ALOGNE are read nowhere (unused?), nor is the array ALOGR (line 64 ↗, unused?).
Other modules: EEEOS06_MOD (XX, for the inverse spacings DFSX); ISHORT from IOUNITS_MOD.

The method

As in READCOEOS01, with 197 temperatures, and rows that hold $T_6$, a mean molecular weight (AMU_M), $\log n_e$ (ALOGNE, by its name) and nine quantities. Unlike the 2001 version, it does not set up the ramp's edge RMX; EQBOUND06 uses NRA in /BEOS06/ (now in EOS06MISC_MOD) instead.

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) The 2006 tables called "OEOS96" in a message.

(lines 114–118 ↗) When a density row of the file has fewer temperatures than NTA requires, the run stops with "Problem with OEOS96 data files," naming the 1996 tables.

A possible fix:

readcoeos06.f

@@ lines 112–118 @@
       read(IOPALE,'(a)') blank
       read(IOPALE,'(a)') blank
       if (icycuse(midx,jcs) .lt. nta(jcs)) then
-         write (ISHORT,'("Problem with OEOS96 data files: X=",f6.4,
+C 11/09/26 KC via Claude
+C code changed to name the 2006 tables in the message.
+         write (ISHORT,'("Problem with OEOS06 data files: X=",f6.4,
      X         " density=",e14.4)') xin(midx), rhogr(midx,jcs)
          stop
       endif

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 message was never reached in the 31 runs with the OPAL 2006 tables, since every density row of the shipped file is long enough; only a damaged or different file reaches it, and the run then stops. No fix was run, since it changes only a message. Intent: none found.

2. (harmless in the runs) `ALOGNE` read into one element.

(lines 124–126 ↗; line 81 ↗, line 89 ↗) The data read is alogNe(jcs,j), but the loop variable is i: j is left from the initialization loops (11 after them), so every row writes to alogNe(jcs,11), and ALOGNE holds only the last value read for each density. ALOGNE is read nowhere, so the results are not changed today.

A possible fix:

readcoeos06.f

@@ lines 121–128 @@
          read (IOPALE,'(a)') blank
          go to 4
       endif
+C 11/09/26 KC via Claude
+C code changed to read ALOGNE into the row of the loop.
       read (IOPALE,'(f11.6,1x,f6.4,e11.4,2e13.6,2e11.3,5f10.6)')
-     x t6list(jcs,i),amu_M(jcs,i),alogNe(jcs,j),
+     x t6list(jcs,i),amu_M(jcs,i),alogNe(jcs,i),
      x (xz(midx,index(iv),i,jcs),iv=1,9)
     4 continue
       enddo

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged). Each of the 31 runs with the OPAL 2006 tables read 117,440 data rows, and 116,595 of them went to alogNe(jcs,11) instead of their own row (j = 11 throughout); the others are the rows with i = 11. No fix was run: ALOGNE is read nowhere, so no output can change. Intent: none found.

Definitions

subroutineREADCOEOS06

lines 7–163 · GitHub

Uses modules AEOS06_MOD, EEEOS06_MOD, EOS06MISC_MOD, IOUNITS_MOD

Called by SETEOS06

Compiler warnings (5)

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

line 21: Unused variable ‘alogr’ declared [-Wunused-variable]
line 28: Unused variable ‘lopale06’ declared [-Wunused-variable]
line 28: Unused variable ‘lnumderiv’ declared [-Wunused-variable]
line 75: Legacy Extension: Blanket SAVE statement follows previous SAVE statement
line 137: Equality comparison for REAL(8) [-Wcompare-reals]

Change log

The changes since rp_1

Changed since rp_1. Declarations only: COMMON blocks to modules (/BEOS06/, /EEOS06/ and /AAEOS06/ merged into EOS06MISC_MOD; /EEEEOS06/ made local arrays, with SAVE; /OPALEOS/ dropped with the OPAL 1995 and 2001 equations of state), IMPLICIT NONE; AEOS06's M renamed MIDX.

1C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 2C 3C READCOEOS06 4C 5C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 6 7 subroutine readcoeos06 8C 10/14/26 MHP via Claude 9C code changed to convert COMMON/IOUNITS/ into IOUNITS_MOD. 10 USE IOUNITS_MOD 11C 10/20/26 MHP via Claude 12C code changed to convert COMMON/AEOS06/, COMMON/EOS06MISC/, COMMON/EEEOS06/ into 13C AEOS06_MOD, EOS06MISC_MOD, EEEOS06_MOD. 14 USE AEOS06_MOD 15 USE EEEOS06_MOD 16 USE EOS06MISC_MOD 17C 10/20/26 MHP via Claude 18C code changed to replace implicit typing with IMPLICIT NONE and 19C explicit declarations. 20 IMPLICIT NONE 21 REAL*8 ALOGNE,ALOGR,AMU_M,DUM,FRAC,RHOGR,TMASS,XIN 22 INTEGER I,ICYCUSE,ITIMECO,IV,J,JCS,K,L,MV,MX,NR,NT,NUMTOT 23 24c..... The purpose of this subroutine is to read the OPAL 2006 EOS data tables 25 26 parameter (mx=5,mv=10,nr=169,nt=197) 27 real*8 moles 28 LOGICAL*4 lopale06,LNumDeriv 29C MHP 8/25 Remove unused variables 30C 9/24/26 MHP via Claude 31C code changed to remove the OPAL 1995/2001 equations of state 32 character*1 blank 33C 10/9/26 MHP via Claude 34C code changed to rename AEOS06's M member to MIDX throughout: M is 35C such a universally-used name elsewhere in this codebase (e.g. shell 36C count in calcad.f) that keeping it here was just asking for a 37C future collision. Used in this file as the composition/mixture-axis 38C loop index below. 39C KC 2025-05-30 reordered common block elements 40C common/beos06/ iri(10),index(10),nta(nr),zz(mx),nra(nt) 41C 10/9/26 MHP via Claude 42C code changed to merge COMMON/EEOS06/, COMMON/AAEOS06/, COMMON/BBEOS06/, 43C and COMMON/BEOS06/ into COMMON/EOS06MISC/ (see esac06.f for the full 44C rationale). AEOS06 and EEEOS06 are left untouched. 45C Amu_M and alogNe were incorrectly dimensoned Amu_M(mx,nr) and alogNe(mx,nr). 46C the required dimensions are Amu_M(nr,nt) alogNe(nr,nt), which is much larger. 47C Something was being overwritten. llp 8/19/08 48C x amu_M(mx,nr),alogNe(mx,nr),rhogr(mx,nr),frac(mx,6), 49C KC 2025-05-30 reordered common block elements 50C common/eeeeos06/moles(mx),xin(mx),tmass(mx),icycuse(mx,nr), 51C 10/9/26 MHP via Claude 52C code changed to convert COMMON/EEEEOS06/ to plain local arrays: 53C declared only in this file (confirmed via a codebase-wide scan), so 54C there's nothing to share via COMMON -- and nothing here reads these 55C values back after they're populated either, so this subroutine's 56C value-across-calls behavior is unaffected either way (it's only 57C ever called once, from esac06.f). Added SAVE below anyway, which 58C costs nothing (a SAVE'd local array and a COMMON-block array are 59C both ordinary static storage, same footprint) and protects against 60C any future refactor that calls this more than once. 61 dimension moles(mx),xin(mx),tmass(mx), 62 x x amu_M(nr,nt),alogNe(nr,nt),rhogr(mx,nr),frac(mx,6), 63C x alogr(nr,nt) 64 x x alogr(nr,nt),icycuse(mx,nr) 65 save 66C 9/25/26 MHP via Claude 67C code changed to consolidate all I/O logical unit numbers, previously 68C scattered across COMMON/LUOUT/, COMMON/LUNUM/, COMMON/ALEX06/, 69C COMMON/ACDPTH/, COMMON/NEWOPAC/, COMMON/LOPAL95/, COMMON/ATMOS2/, 70C COMMON/ALATM03/, COMMON/OPALEOS/, and COMMON/SCV2/, into a single 71C COMMON/IOUNITS/ 72C 10/13/26 MHP via Claude: removed the unused COMMON/NMLPAR/ declaration 73C (no member referenced here; residue of the obsolete-option removals). 74 75 save 76 blank=' ' 77 78 79 if (itimeco .ne. 12345678) then 80 do i=1,mx 81 do j=1,mv 82 do k=1,nt 83 do l=1,nr 84 xz(i,j,k,l)=1.D+35 85 enddo 86 enddo 87 enddo 88 enddo 89 do j = 1,mv 90 eos(j) = 1D0 91 enddo 92 itimeco=12345678 93 endif 94 95 close (2) 96c.....read tables 97C MHP 8/25 Moved opening of file to parmin 98 do 3 midx=1,mx 99 read (IOPALE,'(3x,f6.4,3x,f12.9,11x,f10.7,17x,f10.7)') 100 x x xin(midx),zz(midx),moles(midx),tmass(midx) 101 read (IOPALE,'(21x,e14.7,4x,e14.7,3x,e11.4,3x,e11.4,3x,e11.4, 102 x x 4x,e11.4)') (frac(midx,i),i=1,6) 103 read (IOPALE,'(a)') blank 104 do 2 jcs=1,nr 105 read (IOPALE,'(2i5,2f12.7,17x,e15.7)') numtot,icycuse(midx,jcs) 106 x x ,dum,dum,rhogr(midx,jcs) 107 if(numtot .ne. jcs) then 108 write (ISHORT,'(" OEOS06 Data file incorrect: numtot,jcs= " 109 x x , 2i5)') numtot,jcs 110 stop 111 endif 112 read(IOPALE,'(a)') blank 113 read(IOPALE,'(a)') blank 114 if (icycuse(midx,jcs) .lt. nta(jcs)) then 115 write (ISHORT,'("Problem with OEOS96 data files: X=",f6.4, 116 X X " density=",e14.4)') xin(midx), rhogr(midx,jcs) 117 stop 118 endif 119 do i=1,icycuse(midx,jcs) 120 if (i .gt. nta(jcs)) then 121 read (IOPALE,'(a)') blank 122 go to 4 123 endif 124 read (IOPALE,'(f11.6,1x,f6.4,e11.4,2e13.6,2e11.3,5f10.6)') 125 x x t6list(jcs,i),amu_M(jcs,i),alogNe(jcs,j), 126 x x (xz(midx,index(iv),i,jcs),iv=1,9) 127 4 4 continue 128 enddo 129 read(IOPALE,'(a)') blank 130 read(IOPALE,'(a)') blank 131 read(IOPALE,'(a)') blank 132 2 2 continue 133 read(IOPALE,'(a)') blank 134 3 3 continue 135 136 do i=1,nt 137 if(t6list(1,i) .eq. 0D0) then 138 write(ISHORT,'("READCOEOS06: Error:",i4, 139 $ $ "-th T6 value is zero")') i 140 stop 141 endif 142 t6a(i)=t6list(1,i) 143 enddo 144 do 12 i=2,nt 145C KC 2025-05-30 fixed "DO termination statement which is not END DO or CONTINUE" 146C 12 dfs(i)=1D0/(t6a(i)-t6a(i-1)) 147 dfs(i)=1D0/(t6a(i)-t6a(i-1)) 148 12 12 continue 149 rho(1)=rhogr(1,1) 150 do 13 i=2,nr 151C KC 2025-05-30 fixed "DO termination statement which is not END DO or CONTINUE" 152C rho(i)=rhogr(1,i) 153C 13 dfsr(i)=1D0/(rho(i)-rho(i-1)) 154 rho(i)=rhogr(1,i) 155 dfsr(i)=1D0/(rho(i)-rho(i-1)) 156 13 13 continue 157 do i=2,mx 158 dfsx(i)=1D0/(xx(i)-xx(i-1)) 159 enddo 160 CLOSE (IOPALE) 161 162 return 163 end

1"""Python version of readcoeos06.f (YREC 7be4133): reads the OPAL 2006 2equation-of-state tables.""" 3 4import numpy as np 5 6from pyrec import aeos06_mod 7from pyrec import eeeos06_mod 8from pyrec import eos06misc_mod 9from pyrec import iounits_mod 10 11MX = 5 12MV = 10 13NR = 169 14NT = 197 15 16# Locals kept between calls (SAVE), formerly in part COMMON/EEEEOS06/. 17moles = np.zeros(MX + 1) 18xin = np.zeros(MX + 1) 19tmass = np.zeros(MX + 1) 20amu_m = np.zeros((NR + 1, NT + 1)) 21alogne = np.zeros((NR + 1, NT + 1)) 22rhogr = np.zeros((MX + 1, NR + 1)) 23frac = np.zeros((MX + 1, 7)) 24alogr = np.zeros((NR + 1, NT + 1)) 25icycuse = np.zeros((MX + 1, NR + 1), dtype=int) 26itimeco = 0 27j = 0 # saved: the leftover index used in alogne[jcs, j] below 28 29 30def _fields(line: str, widths: list, start: int = 0, ints: tuple = ()) -> list: 31 """Split a fixed-format record into numbers (a blank field reads as 0, 32 as in Fortran); the fields numbered in ints are integers. 33 34 """ 35 36 out, pos = [], start 37 for n, w in enumerate(widths): 38 text = line[pos:pos + w].strip() 39 if n in ints: 40 out.append(int(text) if text else 0) 41 else: 42 out.append(float(text.replace('D', 'E')) if text else 0.0) 43 pos += w 44 return out 45 46 47def readcoeos06() -> None: 48 """Read the OPAL 2006 EOS tables (five hydrogen abundances) from unit 49 IOPALE into aeos06_mod (xz, t6list, and the grids t6a, rho with their 50 inverse spacings dfs, dfsr, dfsx) and eos06misc_mod (zz); on the 51 first call, fill xz with 1e35 and set eos06misc_mod.eos to 1. 52 53 Translates: readcoeos06.f#READCOEOS06 (7be4133) 54 55 No arguments; nothing is returned. 56 57 """ 58 59 global itimeco, j 60 a = aeos06_mod 61 e = eos06misc_mod 62 iopale = iounits_mod.iopale 63 ishort = iounits_mod.ishort 64 65 # On the first call, fill the tables with 1e35 and EOS with 1 66 # (Fortran lines 79–93). 67 if itimeco != 12345678: 68 for i in range(1, MX + 1): 69 for j in range(1, MV + 1): 70 for k in range(1, NT + 1): 71 for l in range(1, NR + 1): 72 a.xz[i, j, k, l] = 1.0e35 73 for j in range(1, MV + 1): 74 e.eos[j] = 1.0 75 j = MV + 1 # the value a Fortran DO leaves 76 itimeco = 12345678 77 78 # Close unit 2, which nothing opens (Fortran line 95). 79 iounits_mod.close(2) 80 81 # Read the tables (the file is opened in parmin) (Fortran lines 96–134). 82 for midx in range(1, MX + 1): 83 84 # The composition's header (Fortran lines 99–103). 85 a.midx = midx 86 xin[midx], e.zz[midx], moles[midx], tmass[midx] = _fields( 87 iounits_mod.read(iopale), [6, 3, 12, 11, 10, 17, 10], 88 start=3)[0::2] 89 frac[midx, 1:7] = _fields( 90 iounits_mod.read(iopale), 91 [14, 4, 14, 3, 11, 3, 11, 3, 11, 4, 11], start=21)[0::2] 92 iounits_mod.read(iopale) 93 94 for jcs in range(1, NR + 1): 95 96 # The density's header (Fortran lines 105–118). 97 numtot, icycuse[midx, jcs], dum, dum, _, rhogr[midx, jcs] = \ 98 _fields(iounits_mod.read(iopale), [5, 5, 12, 12, 17, 15], 99 ints=(0, 1)) 100 if numtot != jcs: 101 iounits_mod.write( 102 ishort, f" OEOS06 Data file incorrect: numtot,jcs= " 103 f"{numtot:5d}{jcs:5d}") 104 raise SystemExit 105 iounits_mod.read(iopale) 106 iounits_mod.read(iopale) 107 if icycuse[midx, jcs] < e.nta[jcs]: 108 # Suspected bug (see the page): "OEOS96" for the 2006 109 # tables. 110 iounits_mod.write( 111 ishort, f"Problem with OEOS96 data files: X=" 112 f"{xin[midx]:6.4f} density={rhogr[midx, jcs]:14.4E}") 113 raise SystemExit 114 115 # The rows, the temperatures beyond nta skipped (Fortran lines 116 # 119–131). 117 for i in range(1, icycuse[midx, jcs] + 1): 118 if i > e.nta[jcs]: 119 iounits_mod.read(iopale) 120 continue # go to 4 121 vals = _fields(iounits_mod.read(iopale), 122 [11, 1, 6, 11, 13, 13, 11, 11, 123 10, 10, 10, 10, 10]) 124 a.t6list[jcs, i] = vals[0] 125 amu_m[jcs, i] = vals[2] 126 # Suspected bug (see the page): j, not i (j is left over). 127 alogne[jcs, j] = vals[3] 128 for iv in range(1, 10): 129 a.xz[midx, e.index[iv], i, jcs] = vals[3 + iv] 130 iounits_mod.read(iopale) 131 iounits_mod.read(iopale) 132 iounits_mod.read(iopale) 133 iounits_mod.read(iopale) 134 135 # The grids and their inverse spacings (Fortran lines 136–160). 136 for i in range(1, NT + 1): 137 if a.t6list[1, i] == 0.0: 138 iounits_mod.write( 139 ishort, f"READCOEOS06: Error:{i:4d}-th T6 value is zero") 140 raise SystemExit 141 a.t6a[i] = a.t6list[1, i] 142 for i in range(2, NT + 1): 143 a.dfs[i] = 1.0 / (a.t6a[i] - a.t6a[i - 1]) 144 a.rho[1] = rhogr[1, 1] 145 for i in range(2, NR + 1): 146 a.rho[i] = rhogr[1, i] 147 a.dfsr[i] = 1.0 / (a.rho[i] - a.rho[i - 1]) 148 for i in range(2, MX + 1): 149 a.dfsx[i] = 1.0 / (eeeos06_mod.xx[i] - eeeos06_mod.xx[i - 1]) 150 iounits_mod.close(iopale)

readcoeos06.f

Purpose

How it fits

The method

Suspected bugs, checked by running the code

1. (never occurred) The 2006 tables called "OEOS96" in a message.

2. (harmless in the runs) ALOGNE read into one element.

Definitions

subroutineREADCOEOS06

Change log

2. (harmless in the runs) `ALOGNE` read into one element.