op953d.f

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

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

Variables on this page (15)

bd, bdd, bt, btt, bx, cappa, jd, jt, jx, lexh, lexth, lextl, rho_t63l, rho_t63ll, rlh

Purpose

Gets the OPAL95 opacity from the table pre-interpolated to the model's $Z$ (MHP 7/98: "at fixed Z"): a three-dimensional interpolation in $X$, $\log T$ and $\log R$ ($R = \rho/T_6^3$). The same as OP952D, with one more interpolation in $X$.

How it fits

Inputs: the table CAPPA(10,70,19) of $\log\kappa$ at the model's $Z$ in LLOT95_MOD, built by OP95ZTAB; the indices and weights in OP95INTERP_MOD (BX(1,·) for $X$; formerly /OP95INDX/, /OP95FACT/ and /OP95EXT/), from GETOPAL95.
Who calls it: GETOPAL95, when the point's $Z$ matches the table's ZTAB (to a relative $10^{-4}$) but its $X$ does not match the surface table's, and the acoustic depth is not being computed.
Outputs: $\kappa$, $\log\kappa$, QOD and QOT, as in OP952D.

The method

For each of the four nearest $X$ tables, the interpolation of OP952D in $\log R$ (with linear extrapolation beyond the table) and in $\log T$ gives $\log\kappa$ and its two derivatives; these three are then interpolated in $X$ with four-point weights BX(1,J), and $\partial\log\kappa/\partial\log T$ is converted to constant $\rho$ as in OP952D.

Definitions

subroutineOP953D

lines 6–89 · GitHub

Uses modules LLOT95_MOD, OP95INTERP_MOD

Called by GETOPAL95

Change log

The changes since rp_1

Changed since rp_1. Declarations only: COMMON blocks to modules (the three interpolation blocks merged into OP95INTERP_MOD, with RL, RLL renamed RHO_T63L, RHO_T63LL), the dimension parameters moved into LLOT95_MOD, IMPLICIT NONE.

1C 2C 3CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 4C OP953D 5CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 6 SUBROUTINE OP953D(O,OL,QOD,QOT) 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/LLOT95/ into LLOT95_MOD. 12 USE LLOT95_MOD 13C 10/20/26 MHP via Claude 14C code changed to replace implicit typing with IMPLICIT NONE and 15C explicit declarations. 16 IMPLICIT NONE 17 REAL*8 C,CC,CCD,CCT,DELTAR,O,OL,QC,QOD,QOT 18 INTEGER I,ID,IT,IX,J 19C MHP 7/98 GET OPAL95 OPACITY FROM TABLE AT FIXED Z 20C AS ABOVE FOR THE MODEL Z. 21C INDICES FOR INTERPOLATION IN Z,X,T, AND R; INTERPOLATION FACTORS FOR 22C Z,X,T, AND R (AND DERIVATIVE FACTORS FOR T AND RHO); AND DATA FOR 23C LINEAR EXTRAPOLATION WHEN OUTSIDE THE TABLE IN RHO. 24C 10/9/26 MHP via Claude 25C code changed to merge COMMON/OP95INDX/, COMMON/OP95FACT/, and 26C COMMON/OP95EXT/ into COMMON/OP95INTERP/ (see getopal95.f for the 27C full rationale). 28C 10/13/26 MHP via Claude 29C code changed to rename COMMON/OP95INTERP/'s RL and RLL to RHO_T63L and 30C RHO_T63LL. RL is log(rho/T6**3), set here as RL = DL - 3*(TL-6); the 31C old names clash with RL and RLL as used for radius throughout the rest 32C of the code, which is why main.f declared only a 10-member prefix of 33C this block to avoid importing them. With the clash gone main.f can 34C declare the full block, which a USE statement will require once this 35C becomes a module. 36 DIMENSION C(4),QC(4),CC(4),CCT(4),CCD(4) 37 SAVE 38C FIND OPACITY AT EACH OF THE 4 NEARBY VALUES OF X 39 DO J = 1,4 40 IX = JX(1,J) 41C FOR EACH X, GET CAPPA FOR 4 VALUES OF T 42 DO I = 1,4 43 IT = JT + I - 1 44C FOR EACH T, GET CAPPA FOR 4 VALUES OF R 45 ID = JD(I) 46C LOG CAPPA AT DESIRED RHO FOR EACH OF THE 4 DESIRED T. 47 C(I) = BD(I,1)*CAPPA(IX,IT,ID)+ BD(I,2)*CAPPA(IX,IT,ID+1) 48 * * + BD(I,3)*CAPPA(IX,IT,ID+2) + BD(I,4)*CAPPA(IX,IT,ID+3) 49C D LOG CAPPA/D LOG R FOR EACH OF THE 4 DESIRED T. 50 QC(I)=BDD(I,1)*CAPPA(IX,IT,ID)+BDD(I,2)*CAPPA(IX,IT,ID+1) 51 * * +BDD(I,3)*CAPPA(IX,IT,ID+2)+BDD(I,4)*CAPPA(IX,IT,ID+3) 52 END DO 53C CHECK ON WHETHER RHO IS OUTSIDE OF THE TABLE AND NEEDS EXTRAPOLATION 54 IF(LEXTL)THEN 55C FACTOR IN R 56 DELTAR = RHO_T63L-RHO_T63LL 57C CORRECT CAPPA BY USING THE DERIVATIVE AT THE BOUNDARY. 58 DO I = 1,4 59 C(I) = C(I)+DELTAR*QC(I) 60 END DO 61 ELSE IF(LEXTH)THEN 62 DO I = 1,4 63 IF(LEXH(I))THEN 64 DELTAR = RHO_T63L - RLH(I) 65 C(I) = C(I)+DELTAR*QC(I) 66 ENDIF 67 END DO 68 ENDIF 69C INTERPOLATE FOR LOG CAPPA IN T. 70 CC(J) = BT(1)*C(1) + BT(2)*C(2) + BT(3)*C(3) + BT(4)*C(4) 71C D LOG CAPPA/D LOG T 72 CCT(J) = BTT(1)*C(1) + BTT(2)*C(2) + BTT(3)*C(3) + BTT(4)*C(4) 73C D LOG CAPPA/D LOG R 74 CCD(J) = BT(1)*QC(1) + BT(2)*QC(2) + BT(3)*QC(3) + BT(4)*QC(4) 75 END DO 76C INTERPOLATE FOR LOG CAPPA IN X. 77 OL = BX(1,1)*CC(1) + BX(1,2)*CC(2) + BX(1,3)*CC(3) + 78 * * BX(1,4)*CC(4) 79C INTERPOLATE FOR QOT IN X. 80 QOT = BX(1,1)*CCT(1) + BX(1,2)*CCT(2) + BX(1,3)*CCT(3) + 81 * * BX(1,4)*CCT(4) 82C INTERPOLATE FOR QOD IN X. 83 QOD = BX(1,1)*CCD(1) + BX(1,2)*CCD(2) + BX(1,3)*CCD(3) + 84 * * BX(1,4)*CCD(4) 85C CORRECT FROM DERIVATE AT FIXED R TO DERIVATIVE AT FIXED RHO. 86 QOT = QOT - 3.0D0*QOD 87 O = 1.0D1**OL 88 RETURN 89 END

1"""Python version of op953d.f (YREC 7be4133): the OPAL95 opacity from the 2table at the model's Z.""" 3 4import numpy as np 5 6from pyrec import llot95_mod 7from pyrec import op95interp_mod 8 9# The Fortran's local arrays (SAVE). 10c = np.zeros(5) 11qc = np.zeros(5) 12cc = np.zeros(5) 13cct = np.zeros(5) 14ccd = np.zeros(5) 15 16 17def op953d(o: float, ol: float, qod: float, qot: float) -> tuple[float, float, float, float]: 18 """Get the OPAL95 opacity from the table pre-interpolated to the 19 model's Z (MHP 7/98): as OP952D (log R, then log T) for each of the four 20 nearest X tables, then four-point interpolation in X, with the indices 21 and weights GETOPAL95 left in op95interp_mod. 22 23 Translates: op953d.f#OP953D (7be4133) 24 25 Arguments: o, ol, qod, qot, passed in with any value. 26 27 Returns (o, ol, qod, qot): kappa, log kappa, dlog kappa/dlog rho and 28 dlog kappa/dlog T (at constant rho). 29 30 """ 31 32 im = op95interp_mod 33 cappa = llot95_mod.cappa 34 35 # Find the opacity at each of the 4 nearby values of X (Fortran lines 39–75). 36 for j in range(1, 5): 37 ix = im.jx[1, j] 38 # For each X, log kappa and its slope in log R at four values of T 39 # (Fortran lines 42–52). 40 for i in range(1, 5): 41 it = im.jt + i - 1 42 id_ = im.jd[i] 43 c[i] = (im.bd[i, 1] * cappa[ix, it, id_] 44 + im.bd[i, 2] * cappa[ix, it, id_ + 1] 45 + im.bd[i, 3] * cappa[ix, it, id_ + 2] 46 + im.bd[i, 4] * cappa[ix, it, id_ + 3]) 47 qc[i] = (im.bdd[i, 1] * cappa[ix, it, id_] 48 + im.bdd[i, 2] * cappa[ix, it, id_ + 1] 49 + im.bdd[i, 3] * cappa[ix, it, id_ + 2] 50 + im.bdd[i, 4] * cappa[ix, it, id_ + 3]) 51 52 # Outside the table in rho: extrapolate with the slope at the edge 53 # (Fortran lines 54–68). 54 if im.lextl: 55 deltar = im.rho_t63l - im.rho_t63ll 56 for i in range(1, 5): 57 c[i] = c[i] + deltar * qc[i] 58 elif im.lexth: 59 for i in range(1, 5): 60 if im.lexh[i]: 61 deltar = im.rho_t63l - im.rlh[i] 62 c[i] = c[i] + deltar * qc[i] 63 64 # Interpolate in T (Fortran lines 70–74). 65 cc[j] = (im.bt[1] * c[1] + im.bt[2] * c[2] + im.bt[3] * c[3] 66 + im.bt[4] * c[4]) 67 cct[j] = (im.btt[1] * c[1] + im.btt[2] * c[2] + im.btt[3] * c[3] 68 + im.btt[4] * c[4]) 69 ccd[j] = (im.bt[1] * qc[1] + im.bt[2] * qc[2] + im.bt[3] * qc[3] 70 + im.bt[4] * qc[4]) 71 72 # Interpolate in X (Fortran lines 77–84). 73 bx = im.bx 74 ol = bx[1, 1] * cc[1] + bx[1, 2] * cc[2] + bx[1, 3] * cc[3] + bx[1, 4] * cc[4] 75 qot = (bx[1, 1] * cct[1] + bx[1, 2] * cct[2] + bx[1, 3] * cct[3] 76 + bx[1, 4] * cct[4]) 77 qod = (bx[1, 1] * ccd[1] + bx[1, 2] * ccd[2] + bx[1, 3] * ccd[3] 78 + bx[1, 4] * ccd[4]) 79 80 # From the derivative at fixed R to the derivative at fixed rho 81 # (Fortran lines 86–87). 82 qot = qot - 3.0 * qod 83 o = 10.0**ol 84 return o, ol, qod, qot