t6rinteos06.f

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

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

Variables on this page (23)

dfs, dfsr, esact, esk, h, ip, iq, ishort, k1, k2, k3, k4, l1, l2, l3, l4, mv, mx, nr, nt, q, rho, t6a

Purpose

T6RINTEOS06 interpolates one equation-of-state quantity in temperature T6 and density for the OPAL 2006 EOS (its comment's words). At rp_1 it had two twins, line for line apart from table sizes and COMMON names: t6rinteos01.f for the 2001 tables and t6rinterp.f for the original OPAL tables (both removed since rp_1); this version has NT = 197 temperatures.

How it fits

Who calls it: only ESAC06.
What it calls: QUADEOS06.
Inputs: SLR and SLT, the density and temperature; from AEOS06_MOD the grids RHO, T6A, the inverse spacings DFS, DFSR and ESK, the values already interpolated in $X$; from EOS06MISC_MOD the grid indices K1–K4, L1–L4 and the switches IP, IQ, set by ESAC06, and the work arrays H, Q.
Output: ESACT in EOS06MISC_MOD, read by ESAC06; the module is also used by OEQOS06, RADSUB06, READCOEOS06, RHOOFP06, EQBOUND06 and CALCAD.

The method

QUADEOS06 gives the quadratic through three points, evaluated at $x$. For each of the temperatures $T_{k_1}, \dots$, a quadratic in density through the rows $l_1$–$l_3$ gives $h_k$, and, if IQ = 3, another through $l_2$–$l_4$ gives $q_k$. Quadratics in temperature through $h_{k_1..k_3}$ and $h_{k_2..k_4}$ (if IP = 3) are then blended linearly: $$f = w\,f_{123} + (1 - w)\,f_{234},\qquad w = (T_{k_3} - T)\,/\,(T_{k_3} - T_{k_2}),$$ the weight from DFS(K3), the inverse grid spacing. The same is done for the $q$'s, and, when IP = 3, the two results are blended in density with $w_\rho = (\rho_{l_3} - \rho)\,/\,(\rho_{l_3} - \rho_{l_2})$ (DFSR(L3)). A result above $10^{15}$ stops the run ("interpolation indices out of range").

Definitions

subroutineT6RINTEOS06

lines 7–104 · GitHub

Interpolation in temperature and density for the OPAL 2006 EOS tables; result in ESACT (EOS06MISC_MOD).

Uses modules AEOS06_MOD, EEEOS06_MOD, EOS06MISC_MOD, IOUNITS_MOD

Called by ESAC06

Flowchart of its calls, in order (6)

for kx=k1,k1+ip, lines 50–61

QUADEOS06 53

if (iq. eq. 3), lines 55–59

QUADEOS06 57

QUADEOS06 66

if (iq. eq. 3), lines 68–72

QUADEOS06 70

if (ip .eq. 3), lines 73–88

QUADEOS06 75

if (iq .eq. 3), lines 81–87

QUADEOS06 84

Compiler warnings (4)

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

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

Change log

The changes since rp_1

Changed since rp_1. Declarations only: COMMON blocks to modules (/AEOS06/ to AEOS06_MOD, /EEEOS06/ to EEEOS06_MOD, /AAEOS06/, /BBEOS06/ and /EEOS06/ to EOS06MISC_MOD, /LUOUT/ to IOUNITS_MOD), IMPLICIT NONE. The question on the two releases is settled: the 2001 version, t6rinteos01.f, was removed.

1C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 2C 3C T6RINTEOS01 4C 5C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 6 7 recursive subroutine t6rinteos06(slr,slt) 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/EEEOS06/, COMMON/AEOS06/, COMMON/EOS06MISC/ into 13C EEEOS06_MOD, AEOS06_MOD, EOS06MISC_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 DIX,DIX2,ESACT2,ESACTQ,ESACTQ2,QUADEOS06,SLR,SLT 22 INTEGER IS,IU,IW,KX,MV,MX,NR,NT 23 24c The purpose of this OPAL 2006 EOS subroutine is to interpolate in 25c T6 and rho 26 27 parameter (mx=5,mv=10,nr=169,nt=197) 28C 10/9/26 MHP via Claude 29C code changed to rename AEOS06's M member to MIDX: M is such a 30C universally-used name elsewhere in this codebase (e.g. shell count 31C in calcad.f) that keeping it here was just asking for a future 32C collision. M is genuinely used (as a composition/mixture-axis 33C index) in esac06.f and readcoeos06.f; declared but never read or 34C written by name in this file. 35C 10/9/26 MHP via Claude 36C code changed to merge COMMON/EEOS06/, COMMON/AAEOS06/, COMMON/BBEOS06/, 37C and COMMON/BEOS06/ into COMMON/EOS06MISC/ (see esac06.f for the full 38C rationale). AEOS06 and EEEOS06 are left untouched. 39C 9/25/26 MHP via Claude 40C code changed to consolidate all I/O logical unit numbers, previously 41C scattered across COMMON/LUOUT/, COMMON/LUNUM/, COMMON/ALEX06/, 42C COMMON/ACDPTH/, COMMON/NEWOPAC/, COMMON/LOPAL95/, COMMON/ATMOS2/, 43C COMMON/ALATM03/, COMMON/OPALEOS/, and COMMON/SCV2/, into a single 44C COMMON/IOUNITS/ 45c 46 save 47 iu=0 48 is=0 49 50 do kx=k1,k1+ip 51 iw=1 52 iu=iu+1 53 h(iu)=quadeos06(is,iw,slr,esk(kx,l1),esk(kx,l2),esk(kx,l3), 54 x x rho(l1),rho(l2),rho(l3)) 55 if(iq. eq. 3) then 56 iw=2 57 q(iu)=quadeos06(is,iw,slr,esk(kx,l2),esk(kx,l3), 58 x x esk(kx,l4),rho(l2),rho(l3),rho(l4)) 59 endif 60 is=1 61 enddo 62c 63 is=0 64 iw=1 65c..... eos(i) in lower-right 3x3(i=i1,i1+2 j=j1,j1+2) 66 esact=quadeos06(is,iw,slt,h(1),h(2),h(3),t6a(k1),t6a(k2), 67 x x t6a(k3)) 68 if (iq. eq. 3) then 69c..... eos(i) upper-right 3x3(i=i1+1,i1+3 j=j1,j1+2) 70 esactq=quadeos06(is,iw,slt,q(1),q(2),q(3),t6a(k1), 71 x x t6a(k2),t6a(k3)) 72 endif 73 if(ip .eq. 3) then 74c..... eos(i) in lower-left 3x3. 75 esact2=quadeos06(is,iw,slt,h(2),h(3),h(4),t6a(k2),t6a(k3), 76 x x t6a(k4)) 77c..... eos(i) smoothed in left 3x4 78 dix=(t6a(k3)-slt)*dfs(k3) 79 esact=esact*dix+esact2*(1D0-dix) 80c endif ! moved to loc a 81 if(iq .eq. 3) then 82 83c..... eos(i) in upper-right 3x3. 84 esactq2=quadeos06(is,iw,slt,q(2),q(3),q(4),t6a(k2),t6a(k3), 85 x x t6a(k4)) 86 esactq=esactq*dix+esactq2*(1D0-dix) 87 endif 88 endif ! loc a 89c 90 if(iq .eq. 3) then 91 dix2=(rho(l3)-slr)*dfsr(l3) 92 if(ip .eq. 3) then 93c..... eos(i) smoothed in both log(T6) and log(R) 94 esact=esact*dix2+esactq*(1-dix2) 95 endif 96 endif 97 if (esact .gt. 1.D+15) then 98 write(ISHORT,'("T6RINTEOS06: Interpolation indices out", 99 x x " of range;please report conditions.")') 100 stop 101 endif 102 103 return 104 end 105

1"""Python version of t6rinteos06.f (YREC 7be4133): interpolation of one 2OPAL 2006 equation-of-state quantity in temperature and density.""" 3 4from pyrec import aeos06_mod 5from pyrec import eos06misc_mod 6from pyrec import iounits_mod 7from pyrec.quadeos06 import quadeos06 8 9 10def t6rinteos06(slr: float, slt: float) -> None: 11 """Interpolate one EOS quantity at density slr and temperature slt on 12 the 4 x 4 grid ESAC06 has chosen, by blending overlapping quadratics; 13 the result goes to eos06misc_mod.esact. 14 15 Translates: t6rinteos06.f#T6RINTEOS06 (7be4133) 16 17 Arguments: slr, slt. Reads the grid indices K1-K4, L1-L4 and the 18 switches IP, IQ from eos06misc_mod, the tables from aeos06_mod. Uses 19 eos06misc_mod.h and .q as work arrays. None is returned. 20 21 """ 22 23 # Fortran lines 46–48. 24 a = aeos06_mod 25 e = eos06misc_mod 26 l1, l2, l3, l4 = e.l1, e.l2, e.l3, e.l4 27 k1, k2, k3, k4 = e.k1, e.k2, e.k3, e.k4 28 ip, iq = e.ip, e.iq 29 iu = 0 30 is_ = 0 31 32 # Quadratics in density at each temperature: h through rows l1-l3, and, 33 # if IQ = 3, q through l2-l4 (Fortran lines 50–61). 34 for kx in range(k1, k1 + ip + 1): 35 iw = 1 36 iu = iu + 1 37 e.h[iu] = quadeos06(is_, iw, slr, a.esk[kx, l1], a.esk[kx, l2], 38 a.esk[kx, l3], a.rho[l1], a.rho[l2], a.rho[l3]) 39 if iq == 3: 40 iw = 2 41 e.q[iu] = quadeos06(is_, iw, slr, a.esk[kx, l2], a.esk[kx, l3], 42 a.esk[kx, l4], a.rho[l2], a.rho[l3], 43 a.rho[l4]) 44 is_ = 1 45 46 # Quadratics in temperature, blended if IP = 3 (Fortran lines 63–88). 47 is_ = 0 48 iw = 1 49 e.esact = quadeos06(is_, iw, slt, e.h[1], e.h[2], e.h[3], a.t6a[k1], 50 a.t6a[k2], a.t6a[k3]) 51 if iq == 3: 52 esactq = quadeos06(is_, iw, slt, e.q[1], e.q[2], e.q[3], a.t6a[k1], 53 a.t6a[k2], a.t6a[k3]) 54 if ip == 3: 55 esact2 = quadeos06(is_, iw, slt, e.h[2], e.h[3], e.h[4], a.t6a[k2], 56 a.t6a[k3], a.t6a[k4]) 57 dix = (a.t6a[k3] - slt) * a.dfs[k3] 58 e.esact = e.esact * dix + esact2 * (1.0 - dix) 59 if iq == 3: 60 esactq2 = quadeos06(is_, iw, slt, e.q[2], e.q[3], e.q[4], 61 a.t6a[k2], a.t6a[k3], a.t6a[k4]) 62 esactq = esactq * dix + esactq2 * (1.0 - dix) 63 64 # The blend in density (Fortran lines 90–96). 65 if iq == 3: 66 dix2 = (a.rho[l3] - slr) * a.dfsr[l3] 67 if ip == 3: 68 e.esact = e.esact * dix2 + esactq * (1 - dix2) 69 70 # A result out of range stops the run (Fortran lines 97–101). 71 if e.esact > 1.0e15: 72 iounits_mod.write(iounits_mod.ishort, 73 "T6RINTEOS06: Interpolation indices out of " 74 "range;please report conditions.") 75 raise SystemExit()