eqbound06.f

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

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

Variables on this page (10)

k3, l2, mv, mx, nr, nra, nt, nta, rho, t6a

Purpose

EQBOUND06 decides whether $(T,\rho)$ lies inside the OPAL 2006 equation-of-state table and, near its edges, returns the weight for blending OPAL with the Yale/SCV equation of state ("LLP 10/16/06"). Unlike the 1995 and 2001 versions, it follows the table's irregular edge in both temperature and density, and combines the ramps.

How it fits

Who calls it: only EQSTAT2, after OEQOS06, when LOPALE06 is set and $1870 \le T \le 2\times10^8$ K and $\log\rho \le 7$ (EQSTAT2's own density).
Inputs and outputs: T (K) and DL ($\log_{10}\rho$) in; LTAB, LRAMP and FAC out.
Modules: AEOS06_MOD (formerly /AEOS06/) holds the table (T6A, RHO), shared with ESAC06, READCOEOS06, RHOOFP06 and T6RINTEOS06. EOS06MISC_MOD (formerly /BEOS06/ and /BBEOS06/) holds NTA (for each density, the index of the lowest allowed temperature) and NRA (for each temperature, the index of the highest allowed density), set by DATA statements in the module (formerly in ESAC06), and the interpolation indices K3 and L2 that ESAC06 found on its last call; this routine starts its search from them.

The method

With $T_6 = T/10^6$, $\rho = 10^{\mathrm{DL}}$, T6A decreasing and RHO increasing:

Out of table if $\rho < \rho_1$, $\rho \ge \rho_{N}$, $T_6 > T_{6,1}$ or $T_6 \le T_{6,N}$.
Ramps at the "easy borders": at high temperature, $f_T = (T_6 - T_{6,1})/(T_{6,2} - T_{6,1})$ when $T_6 > T_{6,2}$; at low density, $f_\rho = (\rho - \rho_1)/(\rho_2 - \rho_1)$ when $\rho < \rho_2$.
Find $k$ and $l$ with $T_{6,k-1} \ge T_6 > T_{6,k}$ and $\rho_{l-1} < \rho \le \rho_l$, by a linear search from K3 and L2. EQSTAT2 calls this routine only after OEQOS06 has returned normally, so ESAC06 has just set K3 and L2 for the same point (settled from the code, Nov 9, 2026: eqstat2.f, lines 313–317 ↗, oeqos06.f, lines 94–102 ↗).

Question 1 (lines 69–88 ↗). Hidden dependence on ESAC06. The search starts from K3 and L2 in EOS06MISC_MOD, left by the last call of ESAC06. That holds at the one call today, but a call without a preceding ESAC06 would start from unset indices (T6A(0) or RHO(0) could be read). Should the indices be passed explicitly, or the search start from a fixed index?
The irregular edge: out of table if $l > $ NRA(k) or $k > $ NTA(l). On the last allowed row or column, a ramp: $f_\rho = (\rho_l - \rho)/(\rho_l - \rho_{l-1})$ if $l = $ NRA(k); $f_T = (T_6 - T_{6,k})/(T_{6,k-1} - T_{6,k})$ if $k = $ NTA(l).
The weight is the product, $f = f_\rho f_T$.

The ramps are linear in $\rho$ and $T_6$, not in their logarithms.

A high-temperature ramp of step 2 ($k = 2$) and an edge ramp in temperature of step 4 ($k = $ NTA(l)) never meet, nor do the two density ramps ($l = 2$ and $l = $ NRA(k)): The smallest NTA is 26 and the smallest NRA is 87 (settled from the code, Nov 9, 2026: eos06misc_mod.f, lines 61–70 ↗), so step 4's assignment of $f_T$ or $f_\rho$ never replaces one from step 2.

Definitions

subroutineEQBOUND06

lines 5–130 · GitHub

In-table test and blending weight $f = f_\rho f_T$ for the OPAL 2006 table, following its irregular edge (NTA, NRA).

Uses modules AEOS06_MOD, EOS06MISC_MOD

Called by EQSTAT2

Compiler warnings (2)

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

line 16: Unused parameter ‘mv’ declared [-Wunused-parameter]
line 16: Unused parameter ‘mx’ declared [-Wunused-parameter]

Change log

The changes since rp_1

Changed since rp_1. Declarations only: COMMON blocks to modules, IMPLICIT NONE.

1C LLP 10/16/06 SUBROUTINE TO CHECK IF POINT IS WITHIN 2C OPAL 2006 EOS TABLE. 3 4C SUBROUTINE EQBOUND06(T,DL,DL0,FAC,LTAB,LRAMP) ! KC 2025-05-31 5 SUBROUTINE EQBOUND06(T,DL,FAC,LTAB,LRAMP) 6C 10/20/26 MHP via Claude 7C code changed to convert COMMON/AEOS06/, COMMON/EOS06MISC/ into 8C AEOS06_MOD, EOS06MISC_MOD. 9 USE AEOS06_MOD 10 USE EOS06MISC_MOD 11C 10/20/26 MHP via Claude 12C code changed to replace implicit typing with IMPLICIT NONE and 13C explicit declarations. 14 IMPLICIT NONE 15 REAL*8 D,DL,DRAMPF,FAC,T,T6,TRAMPF 16 INTEGER MV,MX,NR,NT 17 LOGICAL*4 LRAMP,LTAB 18 Integer*4 k,l 19 parameter (mx=5,mv=10,nr=169,nt=197) 20C 10/9/26 MHP via Claude 21C code changed to rename AEOS06's M member to MIDX: M is such a 22C universally-used name elsewhere in this codebase (e.g. shell count 23C in calcad.f) that keeping it here was just asking for a future 24C collision. M is genuinely used (as a composition/mixture-axis 25C index) in esac06.f and readcoeos06.f; declared but never read or 26C written by name in this file. 27C KC 2025-05-30 reordered common block elements 28C common/beos06/ iri(10),index(10),nta(nr),zz(mx),nra(nt) 29C 10/9/26 MHP via Claude 30C code changed to merge COMMON/EEOS06/, COMMON/AAEOS06/, COMMON/BBEOS06/, 31C and COMMON/BEOS06/ into COMMON/EOS06MISC/ (see esac06.f for the full 32C rationale). AEOS06 (above) is left untouched. 33 34 SAVE 35 T6 = T*1.0D-6 36 D = 10D0**DL 37 38C EXIT IF OUTSIDE TABLE IN RHO 39 IF((D.LT.rho(1)) .OR. (D.GE.rho(nr)))THEN 40 GOTO 9999 ! Out of Table in density. Go to out of table exit. 41 ENDIF 42 43C Exit if outside table in T6 44 IF((T6 .GT. T6a(1)) .OR. (T6 .LE. T6a(nt))) THEN 45 GOTO 9999 ! Out of Table in temperature. Go to out of table exit. 46 ENDIF 47 48C Initialize 49 DRampF = 1D0 ! Initialize density ramp factor to 1 50 TRampF = 1D0 ! Initialize temperature ramp factor to 1 51 LTAB = .TRUE. ! Presume that T6,D point is in OPAL table 52 LRAMP = .FALSE. ! presume ramping is not needed, 53 ! i.e. we are not in the border of the table. 54 55C Check for easy borders 56 IF(T6 .GT. T6a(2)) THEN 57 TRampF = (T6-T6a(1))/(t6a(2)-T6a(1)) ! Ramp in temperature 58 LRAMP = .TRUE. 59 ENDIF 60 IF(D .LT. rho(2)) THEN 61 DRampF = (D - rho(1))/(rho(2)-rho(1)) ! Ramp in density 62 LRAMP = .TRUE. 63 ENDIF 64 65C Find a k such that T6a(k-1) >= T6 > T6a(k) 66C and an l such that rho(l-1) < D <= rho(l) 67C The ESAC06 variables (k=)k3 and (l-)l2 are close to what we need. 68C A linear search will work fine. 69 k = k3 70 10 10 IF(T6 .LE. T6a(k)) then 71 k = k+1 72 goto 10 73 ENDIF 74 20 20 If( T6 .GT. T6a(k-1)) THEN 75 k = k-1 76 goto 20 77 ENDIF 78C We now have: T6a(k-1) >= T6 > T6a(k) 79 80 l = l2 81 30 30 IF(D .GT. rho(l)) THEN 82 l = l+1 83 goto 30 84 ENDIF 85 40 40 IF(D .LE. rho(l-1)) THEN 86 l = l-1 87 goto 40 88 ENDIF 89C We now have: rho(l-1) < D <= rho(l) 90 91C For a given temperature in array T6a with index k, e.e., T6a(k), element 92C nra(k) of array nra contains the index to the max allowed density in 93C array rho. So, for a given T6a(k), the associated D must be less than 94C rho(nra(k)). A valid index in array rho must be less than or equal to 95C nra(k) 96 97 IF(l .GT. nra(k)) then 98 goto 9999 ! Out of table exit in density. 99 ENDIF 100 IF(l .EQ. nra(k)) THEN 101 LRAMP = .TRUE. 102 DRampF = (rho(l)-D)/(rho(l)-rho(l-1)) 103 ENDIF 104 105C For a given density in array rho with index l, i.e., rho(l), element 106C nta(l) of array nta contains the index to the min allowed temperature 107C in array T6a. So for a given rho(l), the associated T6 must be 108C greater than T6a(nta(l)). (Note that T6a is a decreasing array.) A 109C valid index in T6a must be less than or equal to nta(l). 110 111 IF(k .GT. nta(l)) then 112 goto 9999 ! Out of table exit in temperature 113 ENDIF 114 IF(k .EQ. nta(l)) THEN 115 LRAMP = .TRUE. 116 TRampF = (T6-T6a(k))/(T6a(k-1)-T6A(k)) 117 ENDIF 118 119 FAC = DRampF * TRampF 120 return 121 122C OUT OF TABLE EXIT 123 9999 9999 CONTINUE 124 LTAB = .FALSE. ! Not in table 125 LRAMP = .TRUE. ! Turn on ramping 126 FAC = 0D0 ! Set ramping factor to zero 127C This way, out of table results are ramped to zero. 128 RETURN 129 130 END 131c

1"""Python version of eqbound06.f (YREC 7be4133): whether (T, rho) lies in 2the OPAL 2006 equation-of-state table, and the blending weight near its 3edges.""" 4 5from pyrec import aeos06_mod 6from pyrec import eos06misc_mod 7 8# PARAMETER (mx=5, mv=10, nr=169, nt=197): the table's dimensions. 9MX = 5 10MV = 10 11NR = 169 12NT = 197 13 14 15def eqbound06(t: float, dl: float, fac: float, ltab: bool, lramp: bool) -> tuple[float, bool, bool]: 16 """Test whether (T, log rho) lies inside the OPAL 2006 table, following 17 its irregular edge (nta, nra), and give the weight fac = f_rho * f_T for 18 blending OPAL with the Yale/SCV equation of state near the edges. 19 20 Translates: eqbound06.f#EQBOUND06 (7be4133) 21 22 Arguments: t (K), dl (log10 rho); fac, ltab, lramp (set only). 23 24 Returns (fac, ltab, lramp): the weight (0 out of table), whether the 25 point is in the table, and whether ramping is on (True out of table). 26 27 """ 28 29 rho = aeos06_mod.rho 30 t6a = aeos06_mod.t6a 31 nta = eos06misc_mod.nta 32 nra = eos06misc_mod.nra 33 t6 = t * 1.0e-6 34 d = 10.0**dl 35 36 while True: # one pass; break goes to the out-of-table exit, label 9999 37 38 # Out of the table's rectangle (Fortran lines 38–46). 39 if d < rho[1] or d >= rho[NR]: 40 break # out of table in density 41 if t6 > t6a[1] or t6 <= t6a[NT]: 42 break # out of table in temperature 43 44 # Initialize (Fortran lines 48–53). 45 drampf = 1.0 46 trampf = 1.0 47 ltab = True 48 lramp = False 49 50 # Ramps at the easy borders, linear in T6 and rho (Fortran lines 51 # 55–63). 52 # Question (see the page): linear, where the older OPAL tables 53 # ramped in log rho. 54 if t6 > t6a[2]: 55 trampf = (t6 - t6a[1]) / (t6a[2] - t6a[1]) 56 lramp = True 57 if d < rho[2]: 58 drampf = (d - rho[1]) / (rho[2] - rho[1]) 59 lramp = True 60 61 # Find k, l by a linear search from ESAC06's last indices (Fortran 62 # lines 65–89). 63 # Question (see the page): k3, l2 are left by ESAC06's last call. 64 k = eos06misc_mod.k3 65 while t6 <= t6a[k]: # label 10 66 k = k + 1 67 while t6 > t6a[k - 1]: # label 20 68 k = k - 1 69 l = eos06misc_mod.l2 70 while d > rho[l]: # label 30 71 l = l + 1 72 while d <= rho[l - 1]: # label 40 73 l = l - 1 74 75 # The irregular edge (Fortran lines 91–117). 76 if l > nra[k]: 77 break # out of table in density 78 if l == nra[k]: 79 lramp = True 80 drampf = (rho[l] - d) / (rho[l] - rho[l - 1]) 81 if k > nta[l]: 82 break # out of table in temperature 83 if k == nta[l]: 84 lramp = True 85 trampf = (t6 - t6a[k]) / (t6a[k - 1] - t6a[k]) 86 87 # The weight (Fortran lines 119–120). 88 fac = drampf * trampf 89 return fac, ltab, lramp 90 91 # label 9999: out of table (Fortran lines 122–128). 92 ltab = False 93 lramp = True 94 fac = 0.0 95 return fac, ltab, lramp