oeqos06.f

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

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

Variables on this page (7)

cgas, eos, ishort, mv, mx, nr, nt

Purpose

The interface to the 2006 OPAL equation of state ("LLP, October 17, 2006"), with the same arguments and outputs as OEQOS and OEQOS01 (removed since rp_1); it is now the only OPAL equation of state. Its inline comments document the entries of EOS it uses: EOS(5) is $c_V = dE/dT_6$ at constant volume, EOS(6) is $\partial\log P/\partial\log\rho$ at constant $T_6$, EOS(7) is $\partial\log P/\partial\log T_6$ at constant $\rho$, EOS(8) is $\Gamma_1$, and EOS(9) is $\Gamma_2/(\Gamma_2-1)$.

How it fits

Inputs: T10 $=\log T$, P10 $=\log P$, XC, Z; the tables read by READCOEOS06.
Who calls it: EQSTAT2.
What it calls: RHOOFP06, ESAC06 (IORDER = 9) and MU.
Outputs: $T$, $P$, $\log\rho$ and $\rho$, BETA, BETAI, BETA14, RMU, AMU, EMU, QDT, QDP, QCP, DELA, and the alternate return.
Modules: EOS06MISC_MOD (the result EOS(10), formerly /EEOS06/; shared with ESAC06, RADSUB06, READCOEOS06, RHOOFP06, T6RINTEOS06 and CALCAD), CONSTANTS_MOD (CGAS), IOUNITS_MOD (ISHORT). The unused /COMP/, /CTLIM/, /CONST1/ and /CCOUT2/ of rp_1 are no longer declared. The parameters MX, MV, NR, NT (the table's dimensions, line 51 ↗) are declared but not used here (unused?), and ZTAB is set to Z (line 91 ↗) but no longer passed on (unused?).

The method

Density. From $T$ and $P$ (given as logs T10, P10), with $T_6 = T/10^6$ K and $P_{12} = P/10^{12}$ dyn cm⁻², the table routine finds $\rho(P, T, X)$ (RHOOFP…), with radiation included (IRAD = 1). A value $\le -998$ means failure.
Thermodynamics. The table routine ESAC… fills the array EOS at $(X, T_6, \rho)$. From it: $$\texttt{QDP} = \left(\frac{\partial\ln\rho}{\partial\ln P}\right)_T = \frac{1}{\chi_\rho},\qquad \texttt{QDT} = \left(\frac{\partial\ln\rho}{\partial\ln T}\right)_P = -\frac{\chi_T}{\chi_\rho},$$ $$c_P = 10^{6}\,\frac{c_V\,\Gamma_1}{\chi_\rho}\ \ (\texttt{QCP};\ c_V\ \text{per}\ 10^6\ \mathrm{K}\ \text{in the table}),\qquad \nabla_\mathrm{ad} = \frac{\Gamma_2 - 1}{\Gamma_2}\ \ (\texttt{DELA}),$$ where $\chi_\rho = (\partial\ln P/\partial\ln\rho)_T$, $\chi_T = (\partial\ln P/\partial\ln T)_\rho$.

Question 1 (lines 104–109 ↗). The pressure check. The check that the table reproduces the input pressure, EOS(1) against P12 to $5\times10^{-7}$, is commented out (as it was in OEQOS01, while OEQOS stopped on it; both removed since rp_1). Since RHOOFP06 now accepts a relative error of $10^{-5}$ (rhoofp06.f, line 120 ↗), the check as written would stop runs whose density has converged within that tolerance. Is it left out on purpose, or should it come back as a warning with the looser tolerance?
Radiation. $1-\beta = aT^4/(3P)$, coded as BETA14 $= 2.521971383\times10^{-3}\,T_6^4/P_{12}$; BETA $=\beta$, BETAI $= 1/\beta$.
Mean molecular weight. RMU $= \beta P/(\rho T)$; MU gives the ion and electron terms AMU and EMU; if EMU $\le 0$, it is set to zero and AMU $=$ RMU$/\mathcal{R}$.

The alternate return (RETURN 1) tells the caller to fall back on another equation of state when the point is outside the table or the table fails.

Range. $0.00187 \le T_6 \le 200$.

Definitions

subroutineOEQOS06

lines 30–132 · GitHub

Uses modules CONSTANTS_MOD, EOS06MISC_MOD, IOUNITS_MOD

Called by EQSTAT2

Flowchart of its calls, in order (3)

RHOOFP06 94→ESAC06 102→MU 122

Compiler warnings (4)

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

line 48: Unused parameter ‘mv’ declared [-Wunused-parameter]
line 48: Unused parameter ‘mx’ declared [-Wunused-parameter]
line 48: Unused parameter ‘nr’ declared [-Wunused-parameter]
line 48: Unused parameter ‘nt’ declared [-Wunused-parameter]

Change log

The changes since rp_1

Changed since rp_1. Declarations only: COMMON blocks to modules, IMPLICIT NONE; the unused parameters IVARX, CNVS and ZERO removed.

1CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 2C 3C 2006 OPAL Equation of State 4C 5C LLP October 17, 2006 6C 7CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 8C 9C Routines associated with 2001 OPAL EOS 10C 11C SUBROUTINE OEOS06 12C SUBROUTINE ESAC06 13C SUBROUTINE T6RINTEOS06 14C SUBROUTINE READCOEOS06 15C FUNCTION QUADEOS06 16C FUNCTION GMASS06 17C SUBROUTINE RADSUB06 18C FUNCTION RHOOFP06 19C SUBROUTINE EQBOUND06 20 21C SUBROUTINE MU is also used 22 23 24C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 25C 26C OEQOS06 27C 28C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 29 30 SUBROUTINE OEQOS06(T10,T,P10,P,D10,D,XC,Z,BETA,BETAI, 31 1 1 BETA14,RMU,AMU,EMU,QDT,QDP,QCP,DELA,*) 32C 10/14/26 MHP via Claude 33C code changed to convert COMMON/CONSTANTS/ into CONSTANTS_MOD. 34 USE CONSTANTS_MOD 35C 10/14/26 MHP via Claude 36C code changed to convert COMMON/IOUNITS/ into IOUNITS_MOD. 37 USE IOUNITS_MOD 38C 10/20/26 MHP via Claude 39C code changed to convert COMMON/EOS06MISC/ into 40C EOS06MISC_MOD. 41 USE EOS06MISC_MOD 42C 10/20/26 MHP via Claude 43C code changed to replace implicit typing with IMPLICIT NONE and 44C explicit declarations. 45 IMPLICIT NONE 46 REAL*8 AMU,BETA,BETA14,BETAI,D,D10,DELA,EMU,P,P10,P12,QCP,QDP,QDT, 47 * * R,RHOOFP06,RMU,RMU2,T,T10,T6,X,XC,Z,ZTAB 48 INTEGER IORDER,IRAD,MV,MX,NR,NT 49 50 51 parameter (mx=5,mv=10,nr=169,nt=197) 52 53C DIMENSION ATOMWT(4) 54C 10/13/26 MHP via Claude: removed the unused COMMON/COMP/ declaration 55C (no member referenced here; residue of the obsolete-option removals). 56C 10/13/26 MHP via Claude: removed the unused COMMON/NMLPAR/ declaration 57C (no member referenced here; residue of the obsolete-option removals). 58C 10/13/26 MHP via Claude 59C code changed to drop COMMON/CTLIM/, declared here but never used. 60C 9/25/26 MHP via Claude 61C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/, 62C and COMMON/CONST3/ into a single COMMON/CONSTANTS/ 63C 9/25/26 MHP via Claude 64C code changed to consolidate all I/O logical unit numbers, previously 65C scattered across COMMON/LUOUT/, COMMON/LUNUM/, COMMON/ALEX06/, 66C COMMON/ACDPTH/, COMMON/NEWOPAC/, COMMON/LOPAL95/, COMMON/ATMOS2/, 67C COMMON/ALATM03/, COMMON/OPALEOS/, and COMMON/SCV2/, into a single 68C COMMON/IOUNITS/ 69 70C 10/9/26 MHP via Claude 71C code changed to merge COMMON/EEOS06/, COMMON/AAEOS06/, COMMON/BBEOS06/, 72C and COMMON/BEOS06/ into COMMON/EOS06MISC/ (see esac06.f for the full 73C rationale). 74 75C DATA NZP1/12/ 76C DATA ATOMWT/0.9921D0,0.24975D0,0.08322D0,0.4995D0/ 77 78 SAVE 79C 80 iorder=9 ! gives all 1st and 2nd order data. See instructions 81c in esac01. 82c NOTE: irad=0 does not add radiation; irad=1 adds radiation 83 irad=1 ! does add radiation corrections 84 85 T=10.0D0**T10 86 P=10.0D0**P10 87 T6=T/1.0D6 88 P12=P/1.0D12 89 IF(T6.LT.0.001870D0.OR.T6.GT.200.0D0)RETURN 1 90 X=XC 91 ZTAB=Z 92 93C R=RHOOFP06(X,ZTAB,T6,P12,IRAD) ! KC 2025-05-31 94 R=RHOOFP06(X,T6,P12,IRAD) 95 IF(R.LE.-998.0D0)THEN 96 RETURN 1 97 ENDIF 98 D=R 99 D10=DLOG10(D) 100 101C CALL ESAC06(X,ZTAB,T6,R,IORDER,IRAD,*999) ! KC 2025-05-31 102 CALL ESAC06(X,T6,R,IORDER,IRAD,*999) 103 104c IF(ABS((P12-EOS(1))/P12).GT.0.5D-6)THEN 105c WRITE(ISHORT,*)'***** RUN TERMINATED --ERROR IN OEQOS06 PTOT' 106c WRITE(ISHORT,*) 'P12,EOS(1) Differ: P12,EOS(1),T6,R,X,ZTAB=', 107c * P12,EOS(1),T6,R,X,ZTAB 108c STOP ' ERROR IN OEQOS06 PTOT' 109c ENDIF 110 111 QDP=1.0D0/EOS(6) ! O2006 EOS(6) is dlogP/dlogRho at const T6 112 QDT= -EOS(7)/ EOS(6) ! O2006 EOS(7) is dlogp/dlogT6 at const Rho 113 QCP=1.0D6*EOS(5)*EOS(8)/EOS(6) ! O2006 EOS(5) is the specific heat. dE/dT6 114 ! at const Vol 115 ! O2006 EOS(8) is gamma1 116 DELA=1.0D0/EOS(9) ! O2006 EOS(9) is gamma2/(gamma2-1) 117 118 BETA14=(2.521971383D-3*T6*T6)*(T6*T6/P12) 119 BETA=1.0D0-BETA14 120 BETAI=1.0D0/BETA 121 RMU=P*BETA/(D*T) 122 CALL MU(T,P,D,XC,Z,RMU2,AMU,EMU,BETA) 123 IF(EMU.LE.0.0D0)THEN 124 EMU=0.0D0 125 AMU=RMU/CGAS-EMU 126 ENDIF 127 128 RETURN 129 999 999 WRITE(ISHORT, *)'WARNNING... OPAL TBL FAIL' 130 131 RETURN 1 132 END

1"""Python version of oeqos06.f (YREC 7be4133): the interface to the 2006 2OPAL equation of state (LLP, October 17, 2006).""" 3 4import math 5 6from pyrec import constants_mod 7from pyrec import eos06misc_mod 8from pyrec import iounits_mod 9from pyrec.esac06 import esac06 10from pyrec.mu import mu 11from pyrec.rhoofp06 import rhoofp06 12 13# The table's dimensions (unused here). 14MX, MV, NR, NT = 5, 10, 169, 197 15 16 17def oeqos06(t10: float, t: float, p10: float, p: float, d10: float, d: float, xc: float, z: float, beta: float, 18 betai: float, beta14: float, rmu: float, amu: float, emu: float, qdt: float, qdp: float, qcp: float, 19 dela: float) -> tuple[float, float, float, float, float, float, float, float, float, float, 20 float, float, float, float, bool]: 21 """From log T and log P, find the density in the 2006 OPAL tables 22 (RHOOFP06, radiation included), the thermodynamic derivatives (ESAC06), 23 the radiation pressure fraction and the mean-molecular-weight terms 24 (MU). Fails (alt1) outside 0.00187 <= T6 <= 200 or where the table 25 fails, so that the caller falls back on another equation of state. 26 27 Translates: oeqos06.f#OEQOS06 (7be4133) 28 29 Arguments: t10 = log T, p10 = log P, the composition xc (X) and z; the 30 others are passed in with any value. 31 32 Returns (t, p, d10, d, beta, betai, beta14, rmu, amu, emu, qdt, qdp, 33 qcp, dela, alt1): T and P unlogged, log rho and rho, beta = 1 - aT^4/3P, 34 1/beta and 1 - beta, the mean-molecular-weight terms, dln rho/dln T at 35 constant P, dln rho/dln P at constant T, c_P, nabla_ad, and alt1, True 36 when the Fortran takes its alternate return (RETURN 1). 37 38 """ 39 40 def out(alt1): 41 return (t, p, d10, d, beta, betai, beta14, rmu, amu, emu, qdt, qdp, 42 qcp, dela, alt1) 43 44 # The options, the units, and the range (Fortran lines 80–91). 45 iorder = 9 # all first and second order data 46 irad = 1 # adds radiation corrections 47 t = 10.0**t10 48 p = 10.0**p10 49 t6 = t / 1.0e6 50 p12 = p / 1.0e12 51 if t6 < 0.001870 or t6 > 200.0: 52 return out(True) # RETURN 1 53 x = xc 54 ztab = z # noqa: F841 (unused in the Fortran too) 55 56 # The density (Fortran lines 93–99). 57 r = rhoofp06(x, t6, p12, irad) 58 if r <= -998.0: 59 return out(True) # RETURN 1 60 d = r 61 d10 = math.log10(d) 62 63 # The thermodynamic quantities (Fortran lines 101–102, 129–131). 64 if esac06(x, t6, r, iorder, irad): 65 iounits_mod.write(iounits_mod.ishort, " WARNNING... OPAL TBL FAIL") 66 return out(True) # label 999 67 # Question (see the page): the pressure check is commented out 68 # (Fortran lines 104–109). 69 70 # The derivatives (Fortran lines 111–116). 71 eos = eos06misc_mod.eos 72 qdp = 1.0 / eos[6] # EOS(6): dlogP/dlogRho at const T6 73 qdt = -eos[7] / eos[6] # EOS(7): dlogP/dlogT6 at const Rho 74 qcp = 1.0e6 * eos[5] * eos[8] / eos[6] # EOS(5): dE/dT6, EOS(8): Gamma1 75 dela = 1.0 / eos[9] # EOS(9): Gamma2/(Gamma2-1) 76 77 # Radiation and the mean molecular weight (Fortran lines 118–126). 78 beta14 = (2.521971383e-3 * t6 * t6) * (t6 * t6 / p12) 79 beta = 1.0 - beta14 80 betai = 1.0 / beta 81 rmu = p * beta / (d * t) 82 rmu2, amu, emu = mu(t, p, d, xc, z, 0.0, amu, emu, beta) 83 if emu <= 0.0: 84 emu = 0.0 85 amu = rmu / constants_mod.cgas - emu 86 return out(False)