calcad.f

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

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

Variables on this page (40)

atmbeg, atmmax, atmmin, atmstp, c4pil, cgl, cln, clsunl, crsunl, csigl, deladj, envbeg, envd, envmax, envmin, envp, envr, envstp, envt, eos, iclcd, index, ishort, json, lclcd, lopale06, mv, mx, nr, nt, numenv, spotf, spotx, stotal, tatmos, tauczn, tauhe, tcz, tnorm, whe

Purpose

This file computes acoustic depths: the sound travel time from the surface to the center, $t_\mathrm{norm}$, and from the surface to the base of the convection zone, $t_\mathrm{CZ}$. It joins the interior model and a fresh envelope integration into one profile, computes the sound speed at each point from the equation of state, and integrates $dr/c_s$.

How it fits

Inputs: the arguments, the interior model HR, HD, HP, HT ($\log r$, $\log\rho$, $\log P$, $\log T$), HCOMP (composition), M (last interior point), JENV (the base of the convection zone), RCZL ($\log r_\mathrm{CZ}/R_\odot$), BL, TEFFL, FP, FT, HSTOT, DAGE (the age), CMIXL (the mixing length, passed on to ENVINT).
Who calls it: WRTOUT, when the options LACOUT and LADON are on and the model has a convective envelope (JENV > 1) (wrtout.f, lines 696–700 ↗). LACOUT is in no namelist and set nowhere, so it keeps its initial value, false (acdpth_mod.f, line 55 ↗), and CALCAD never runs today (settled from the code, Nov 9, 2026; see ACDPTH_MOD).
What it calls: ENVINT (the envelope, which fills /ENVSTRUCT/); for the sound speed, ESAC06 (OPAL 2006 EOS) or EQSTAT2; SPLINJ and SPLINT to interpolate at $r_\mathrm{CZ}$; BOOLE to integrate.
Outputs: in /ACDPTH/, TNORM, TCZ, TAUCZN $= t_\mathrm{CZ}/t_\mathrm{norm}$, TATMOS (always 0, see below) and DELADJ ($\nabla_\mathrm{ad}$ at each point). WRTOUT prints them in its output columns (its comment: columns 68–70). A line to ISHORT at every call; with LCLCD, the full profile to unit ICLCD.
Modules: ACDPTH_MOD (with ENVINT, GETOPAL95, the main program, PARMIN, WRTOUT; it also holds the output switches LCLCD, LADON, LACOUT); ENVSTRUCT_MOD (the envelope profile, filled by ENVINT); NMLPAR_MOD (the integration step sizes ATMBEG … ENVMAX and ATMSTP, ENVSTP, the spot parameters SPOTF, SPOTX, and the switch LOPALE06); EOS06MISC_MOD (ESAC06's output EOS); COMP_MOD, CONSTANTS_MOD, IOUNITS_MOD (the units ISHORT, ICLCD) and DIMENSIONS_MOD.

The method

The envelope. Set the atmosphere and envelope step sizes to the fixed values ATMSTP, ENVSTP, then integrate the envelope with ENVINT (with the model's mixing length CMIXL), from the surface to the interior's outer point, with $$\log R = \tfrac12(\log L - \log 4\pi - \log\sigma - 4\log T_\mathrm{eff}), \qquad \log g = \log G + \log M - 2\log R.$$ For a spotted star (SPOTF $\ne 0$, SPOTX $\ne 1$) and only when JENV = M, the envelope uses the ambient temperature $$\log T_\mathrm{amb} = \log T_\mathrm{eff} - \tfrac14\log_{10}\!\left(f_s x_s^4 + 1 - f_s\right),$$ with $f_s$ = SPOTF (the spot filling factor) and $x_s$ = SPOTX (probably the spot-to-ambient temperature ratio).

JENV = M holds when no surface convection zone is found (findsh.f, lines 113–120 ↗), so the spot correction here applies only to stars without one, the opposite of the routines that test LC(M) (stitch.f, suspected bug 6).
One profile. The interior points $1..M$, then envelope points $2..N_\mathrm{env}$, in linear units, $T$ in $10^6$ K.
Sound speed. With the OPAL 2006 EOS, $\Gamma_1$ = EOS(8) and $\nabla_\mathrm{ad} = 1/$EOS(9) (EOS(9) is $\Gamma_2/(\Gamma_2-1)$). Otherwise, from EQSTAT2's derivatives (QDPJ, QDTJ, probably $(\partial\ln\rho/\partial\ln P)_T$ and $(\partial\ln\rho/\partial\ln T)_P$) and $c_P$, by Cox and Giuli's chapter 9: $$\chi_\rho = \frac{1}{\mathtt{QDP}}, \quad \chi_T = -\chi_\rho\,\mathtt{QDT}, \quad c_V = c_P - \frac{P}{\rho T}\frac{\chi_T^2}{\chi_\rho}, \quad \Gamma_1 = \frac{\chi_\rho\,c_P}{c_V}.$$ The code stores the inverse sound speed $1/c_s$, with $c_s = \sqrt{\Gamma_1 P/\rho}$.
Integrate with Boole's rule: $$t_\mathrm{norm} = \int_0^{R} \frac{dr}{c_s}, \qquad t_\mathrm{CZ} = \int_{r_\mathrm{CZ}}^{R} \frac{dr}{c_s},$$ in seconds. For $t_\mathrm{CZ}$, the point just below $r_\mathrm{CZ}$ is replaced by $r_\mathrm{CZ}$ itself, with $1/c_s$ from a cubic spline through 7 neighboring points.

The surface point of the profile is left out of both integrals (suspected bug 8).
$\tau_\mathrm{CZ,n} = t_\mathrm{CZ}/t_\mathrm{norm}$.

The atmosphere's acoustic depth (TATMOS) was computed by a section now commented out (lines 387–437 ↗), so it is always 0.

Design

The header lists IACAT, IJLAST and other output switches, some marked "obsolete?" (IACAT and the other unused unit numbers are no longer in the module). Many local variables and arrays remain from the commented-out sections.

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. (cannot be reached) Step sizes not restored.

(lines 196–207 ↗; lines 233–235 ↗; wrtmod.f, lines 237–242 ↗) The routine saves ATMBEG, ATMMIN, ATMMAX, ENVBEG, ENVMIN, ENVMAX in ABEG0 … EMAX0, sets them to ATMSTP/ENVSTP for its own envelope integration, and never restores them. Later envelope integrations in the same run (the surface boundary condition at every step) then use the fixed steps. WRTMOD, whose code the comment says CALCAD copied (lines 194–195 ↗), does restore them. Not reachable today (LACOUT is set nowhere, above).

A possible fix:

calcad.f

@@ lines 230–238 @@
       ELSE
             ATEFFL = TEFFL
       ENDIF
+C 11/09/26 KC via Claude
+C code changed to restore the integration step sizes after the
+C envelope integration, as WRTMOD does.
       CALL ENVINT(BJ,FPLJ,FTLJ,GLJ,HSTOT,IXXJ,LPRTJ,LSBC0J,
      *         PLIMJ,RLJ,ATEFFL,XJ,ZJ,DUM1,IDUMJ,KATMJ,KENVJ,KSAHAJ,
      *         DUM2,DUM3,DUM4,LPULPTJ,CMIXL)
+      ATMBEG = ABEG0
+      ATMMIN = AMIN0
+      ATMMAX = AMAX0
+      ENVBEG = EBEG0
+      ENVMIN = EMIN0
+      ENVMAX = EMAX0
 C      G Somers END

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): Not reached in any of 34 runs: 26 of YREC's cases, from a 0.04 $M_\odot$ brown dwarf and the solar models to 9 $M_\odot$, rotation, and helium burning, and 8 made-up ones; it cannot be reached at 7be4133, since CALCAD runs only when LACOUT is true, and LACOUT is false in the source and in no namelist (acdpth_mod.f, suspected bug 1). The fix was not run. Since nothing goes wrong until that switch is opened, it is fragile code, to be fixed before LACOUT is made settable. Intent: none found; WRTMOD, which the comment says this code copies, restores the step sizes. Related: STITCH sets the same six step sizes and does not restore them either (stitch.f, suspected bug 1); it runs at each store with LSTORE, LSTCH, and LSTENV on, and its check in the same round found a real but small effect ($\log L$ by up to 5.8e-5 dex, against 9.4e-9 dex for the null change).

2. (cannot be reached) The interior's $Z$.

(lines 246–247 ↗; line 283 ↗) STARZ(I) = HCOMP(3,M) for every interior point: the envelope's $Z$, not the local one, while STARX(I) takes the local $X$. It is used only on the non-OPAL path (EQSTAT2, when LOPALE06 is off, its default). With diffusion or burning that changes $Z$, the interior $Z$ differs. Not reachable today (LACOUT is set nowhere, above).

A possible fix:

calcad.f

@@ lines 244–250 @@
             STARD(I)=DEXP(CLN*HD(I))
             STARP(I)=DEXP(CLN*HP(I))
             STARX(I)=HCOMP(1,I)
-            STARZ(I)=HCOMP(3,M)
+C 11/09/26 KC via Claude
+C code changed to take the local Z at each interior point, as for X.
+            STARZ(I)=HCOMP(3,I)
 34      CONTINUE
       DO 35, I=1,NUMENV-1
             STARR(M+I)=DEXP(CLN*ENVR(I+1))

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): Not reached in the same 34 runs; it cannot be reached at 7be4133, since CALCAD runs only when LACOUT is true, and LACOUT is false in the source and in no namelist (acdpth_mod.f, suspected bug 1). Once it runs, this needs LOPALE06 off (as in testsuite/Test_solarGS98_yaleeos) and a $Z$ that varies inside, as with diffusion. The fix was not run. Since nothing goes wrong until that switch is opened, it is fragile code, to be fixed before LACOUT is made settable. Intent: none found; the envelope loop takes the outer point's $X$ and $Z$ throughout (lines 249–256 ↗), while this loop takes the local $X$.

3. (cannot be reached) `ICZ` may be stale.

(lines 333–351 ↗; lines 370–373 ↗) ICZ, the point just below $r_\mathrm{CZ}$, is set only when a pair of points brackets $r_\mathrm{CZ}$, and is kept between calls (SAVE). Inside the loop, the replacement at I.EQ.ICZ is tested before the bracket is found, so a value from an earlier call that is below this call's bracket matches first: The spline replaces a point below $r_\mathrm{CZ}$, FLG is set, and the real point just below $r_\mathrm{CZ}$ is not replaced, so $t_\mathrm{CZ}$ starts from that point instead of $r_\mathrm{CZ}$. If no pair brackets $r_\mathrm{CZ}$, the old ICZ is used (undefined at the first call; 0 makes loop 27 read STARR(0)). The 7-point stencil ICZ-3 … ICZ+3 is also not checked against the end of the profile. Not reachable today (LACOUT is set nowhere, above).

A possible fix:

calcad.f

@@ lines 330–339 @@


 C      get the location of the convection zone, and count nonzero entries:
+C 11/09/26 KC via Claude
+C code changed to reset ICZ at each call and to find the bracket
+C before the replacement is tested.
                   IENDJ=0
                   FLG=0
+                  ICZ=0
                   DO 22, I=1,K
-                        IF (STARR(I).LE.RCZ(1) .AND. STARR(I+1).GE.RCZ(1)) THEN
+                        IF (ICZ.EQ.0 .AND. STARR(I).LE.RCZ(1) .AND.
+     *                      STARR(I+1).GE.RCZ(1)) THEN
                               ICZ = I
                         ENDIF

(A call with no bracket would still need a decision, such as TCZ = 0.)

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): Not reached in the same 34 runs; it cannot be reached at 7be4133, since CALCAD runs only when LACOUT is true, and LACOUT is false in the source and in no namelist (acdpth_mod.f, suspected bug 1). Once it runs, a stale ICZ needs a second call with a surface convection zone. The fix was not run. Since nothing goes wrong until that switch is opened, it is fragile code, to be fixed before LACOUT is made settable. Intent: none found.

4. (cannot be reached) `FXIONJ` is a scalar.

(line 82 ↗; line 292 ↗; eqstat2.f, line 54 ↗; eqstat2.f, lines 171–173 ↗) CALCAD passes the scalar FXIONJ where EQSTAT2 declares FXION(3) and writes three elements, so EQSTAT2 writes past it, into whatever CALCAD stores next. On the non-OPAL path (LOPALE06 off, its default). Not reachable today (LACOUT is set nowhere, above).

A possible fix:

calcad.f

@@ lines 79–87 @@
 C explicit declarations.
       IMPLICIT NONE
       REAL*8 AMUJ,ATEFFL,BET14J,BETAIJ,BETAJ,BL,CHRH,CHT,CMIXL,CVJ,DAGE,
-     *       DELAJ,EMUJ,ETAJ,FP,FPLJ,FT,FTLJ,FXIONJ,HCOMP,HD,HP,HR,
+     *       DELAJ,EMUJ,ETAJ,FP,FPLJ,FT,FTLJ,HCOMP,HD,HP,HR,
      *       HSTOT,HT,QAPJ,QATJ,QCPJ,QCPPJ,QCPTJ,QDPJ,QDTJ,QDTPJ,QDTTJ,
      *       RCZ,RCZL,RMUJ,STD,STDL,STP,STPL,STT,STTL,STX,STZ,TEFFL
+C 11/09/26 KC via Claude
+C code changed to declare FXIONJ with the three elements EQSTAT2
+C fills.
+      REAL*8 FXIONJ(3)
       INTEGER IORDER,IRAD,JENV,KHI,KLO,M,MV,MX,NR,NT
       LOGICAL*4 LATMOJ,LCALCJ,LDERIJ
 C 10/9/26 MHP via Claude

(with FXIONJ dropped from the REAL*8 list at line 82 ↗).

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): Not reached in the same 34 runs; it cannot be reached at 7be4133, since CALCAD runs only when LACOUT is true, and LACOUT is false in the source and in no namelist (acdpth_mod.f, suspected bug 1). Once it runs, this needs LOPALE06 off. The write past FXIONJ cannot be counted without reading the memory it overwrites, so only the call was to be counted. The fix was not run. Since nothing goes wrong until that switch is opened, it is fragile code, to be fixed before LACOUT is made settable. Intent: none found.

5. (cannot be reached) EOS errors ignored.

(lines 269–273 ↗) ESAC06's alternate return *999 jumps to the next statement, so a failed EOS call uses the EOS values of the previous point; ESAC06 writes a message to ISHORT, but CALCAD goes on. Only with LOPALE06, for points outside the OPAL table. Not reachable today (LACOUT is set nowhere, above).

A possible fix: a choice for the developers: fall back to EQSTAT2 for that point, as the EOS does elsewhere outside the table, or stop.

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): Not reached in the same 34 runs; it cannot be reached at 7be4133, since CALCAD runs only when LACOUT is true, and LACOUT is false in the source and in no namelist (acdpth_mod.f, suspected bug 1). Once it runs, this needs LOPALE06 on and points outside the OPAL table (esac06.f, lines 488–494 ↗). The fix (measurement only: recomputing such a point with EQSTAT2) was not run. Since nothing goes wrong until that switch is opened, it is fragile code, to be fixed before LACOUT is made settable. Intent: none found.

6. (cannot be reached) Output format.

(lines 445–446 ↗) Format 67 has no space before "Acoustic depth of atmopshere" (also misspelled), so the line runs the two together. Not reachable today (LACOUT is set nowhere, above).

A possible fix:

calcad.f

@@ lines 442–449 @@

 C Output acoustic depth info to ISHORT
       WRITE(ISHORT,67)TCZ,TNORM,TAUCZN,TATMOS
+C 11/09/26 KC via Claude
+C code changed to separate and correct the last label of format 67.
 67       FORMAT(1X,'Acoustic depth to CZ:',F14.8,2X,'Acoustic depth to center',
-     *F13.7,2X,'Normalized taucz:',F11.9, 'Acoustic depth of atmopshere:',F16.8,2X)
+     *F13.7,2X,'Normalized taucz:',F11.9,2X,'Acoustic depth of atmosphere:',
+     *F16.8,2X)

 C 555            CONTINUE
 C--------------------------------------------------------------

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): Not reached in the same 34 runs; it cannot be reached at 7be4133, since CALCAD runs only when LACOUT is true, and LACOUT is false in the source and in no namelist (acdpth_mod.f, suspected bug 1). The fix was not run. Since nothing goes wrong until that switch is opened, it is a defect only in output that is never written, to be fixed before LACOUT is made settable. Intent: none found.

7. (cannot be reached) A short profile.

(lines 317–322 ↗; line 441 ↗; boole.f, lines 56–63 ↗) If K ≤ 3, NN = 2 and BOOLE returns 0, so TAUCZN divides by zero. A profile that short does not occur in a real model. Not reachable today (LACOUT is set nowhere, above).

A possible fix:

calcad.f

@@ lines 437–444 @@
 c      CALL BOOLE(ATMOSR,ATMOSC,V,NN,TATMOS)

 C Output normalized acoustic depth
+C 11/09/26 KC via Claude
+C code changed to leave TAUCZN at zero when TNORM is zero.
       TAUCZN=0.0D0
-      TAUCZN=TCZ/TNORM
+      IF (TNORM.GT.0.0D0) TAUCZN=TCZ/TNORM

 C Output acoustic depth info to ISHORT
       WRITE(ISHORT,67)TCZ,TNORM,TAUCZN,TATMOS

8. (cannot be reached) The surface point is left out of the integrals.

(lines 310–322 ↗; lines 241–256 ↗; boole.f, lines 41–52 ↗) Loops 34 and 35 build a profile of $M + N_\mathrm{env} - 1$ points, and the sound speed is computed at all of them (line 264 ↗; line 278 ↗), but K = NUMENV-1+M-1 (line 310 ↗), one point fewer, is passed to BOOLE as the number of points, and BOOLE integrates over its first N points only (boole.f, lines 41–52 ↗). The profile for $t_\mathrm{CZ}$ ends at point K too (lines 370–373 ↗). What follows: Both $t_\mathrm{norm}$ and $t_\mathrm{CZ}$ leave out the outermost interval, between the envelope's last two points, although the sound speed at the surface point is computed as at the others. Not reachable today (LACOUT is set nowhere, above).

A possible fix:

calcad.f

@@ lines 307–314 @@


 C Integrate for the full acoustic depth for normalization
-      K=NUMENV-1+M-1
-      U=MOD((NUMENV-1+M-1),4)
+C 11/10/26 KC via Claude
+C code changed to integrate over all M+NUMENV-1 points of the profile,
+C the surface point included.
+      K=NUMENV-1+M
+      U=MOD(K,4)
       IF (U.NE.0) THEN
             NN=K+4-U+1
       ELSE
@@ lines 332–338 @@
 C      get the location of the convection zone, and count nonzero entries:
                   IENDJ=0
                   FLG=0
-                  DO 22, I=1,K
+C 11/10/26 KC via Claude
+C code changed to stop the search one point short of K, so that
+C STARR(I+1) stays inside the profile.
+                  DO 22, I=1,K-1
                         IF (STARR(I).LE.RCZ(1) .AND. STARR(I+1).GE.RCZ(1)) THEN
                               ICZ = I
                         ENDIF

Checked by running YREC (Nov 10, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): Not reached in five of YREC's cases, testsuite/Test_m1p0feh+0p0GS98_base_TAMS, testsuite/Test_m0p3feh+0p0GS98_base_ZAMS, testsuite/Test_m1p4feh+0p0GS98_dif_TAMS (with diffusion), testsuite/Test_m0p8feh+0p0GS98_rotp8p0d0p1mag1e4std_TAMS (with rotation), and examples/run_standard_solar_model/Test_solarcal_noGS_norot (the solar calibration): A line at line 310 ↗ never counted a call. It cannot be reached at 7be4133, since CALCAD runs only when LACOUT is true, and LACOUT is false in the source and in no namelist (acdpth_mod.f, suspected bug 1). The fix was not run; since nothing goes wrong until that switch is opened, it is fragile code, to be fixed before LACOUT is made settable. Intent: none found.

Definitions

subroutineCALCAD

lines 48–509 · GitHub

Integrates the envelope, joins it to the interior, computes $1/c_s$ from the EOS and integrates it with BOOLE to get the acoustic depth of the star and of its convection zone (TNORM, TCZ, TAUCZN in /ACDPTH/).

Uses modules ACDPTH_MOD, COMP_MOD, CONSTANTS_MOD, DIMENSIONS_MOD, ENVSTRUCT_MOD, EOS06MISC_MOD, IOUNITS_MOD, NMLPAR_MOD

Called by WRTOUT

Flowchart of its calls, in order (9)

ENVINT 233

if (LOPALE06), lines 261–275

for I=1,M+NUMENV-1 (DO 36), lines 264–274

ESAC06 269

if ( .NOT. LOPALE06), lines 277–304

for I=1,M+NUMENV-1 (DO 41), lines 278–303

EQSTAT2 292

BOOLE 322

for I=1,K (DO 22), lines 335–369

if (I.EQ.ICZ .AND. FLG.EQ.0 .AND. ICZ.GT.3), lines 341–351

SPLINJ 347→SPLINT 348

if (I.EQ.ICZ .AND. FLG.EQ.0 .AND. ICZ.LE.3), lines 352–367

else, lines 357–366

SPLINJ 362→SPLINT 363

BOOLE 382

Compiler warnings (6)

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

line 85: Unused parameter ‘mv’ declared [-Wunused-parameter]
line 85: Unused parameter ‘mx’ declared [-Wunused-parameter]
line 85: Unused parameter ‘nr’ declared [-Wunused-parameter]
line 85: Unused parameter ‘nt’ declared [-Wunused-parameter]
line 228: Inequality comparison for REAL(8) [-Wcompare-reals]
line 228: Inequality comparison for REAL(8) [-Wcompare-reals]

Change log

The changes since rp_1

Changed since rp_1. The new argument CMIXL, passed on to ENVINT; otherwise declarations only: COMMON blocks to modules, IMPLICIT NONE, unused blocks removed (among them a wrong /LLOT95/ declaration that aliased ZTAB to an opacity value, harmless since nothing here reads it).

1C 2C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 3C CALCAD: links together and integrates sound speed profiles, for acoustic 4C depth to convection zone, acoustic depth of atmosphere 5C INPUTS: 6C HR: array of radii for the interiors solution 7C RCZ: physical radius of the convection zone (logged), calculated in wrtout 8C M: last interior point 9C HD: interior's density 10C HP: interior pressure 11C HT: interior temperature 12C BL: logged luminosity 13C FP: Rotational distortion, needed for envelope integration 14C FT: Rotational distortion, needed for envelope integration 15C HSTOT: needed for ENVINT 16C LPRT: needed for ENVINT 17C TEFFL: logged effective temperature 18C HCOMP: vector of compositions for interior 19C NKK: obsolete? 20C DAGE: Age 21C DDAGE: obsolete? 22C 23C 24C 25C 26C OUTPUTS: 27C Adjusts the values in the common block for acoustic depth: 28C TAUCZN = acoustic depth to CZ / sound travel time from surface-center 29C TCZ = acoustic depth to CZ (seconds) 30C TNORM = sound travel time from surface to center (seconds) 31C DELADJ = adiabatic gradient from OPAL 2006 EOS or SCV EOS 32C 33C 34C 35C 36C 37C COMMON BLOCKED THINGS: 38C 39C LCLCD = Logical flag for output of CALCAD files at AGEOUT ages 40C AGEOUT(5) = Ages for which you want output. Hard coded in parmin 41C IACAT = 42C IJLAST = Again, for output 43C LJLAST = Again, for output 44C LJWRT = ouput on/off toggle 45C LADON = toggle on calcad imposed 4d interpolation in opacity 46C 47C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 48 SUBROUTINE CALCAD(HR, RCZL, M, HD, HP, HT, BL, FP, FT, HSTOT, 49C * LPRT, TEFFL, HCOMP, NKK, DAGE, DDAGE, JENV) ! KC 2025-05-31 50 * * TEFFL, HCOMP, DAGE, JENV, CMIXL) 51C 10/15/26 MHP via Claude 52C code changed to take JSON from DIMENSIONS_MOD instead of declaring 53C PARAMETER(JSON=5000) here. 54 USE DIMENSIONS_MOD 55C 10/14/26 MHP via Claude 56C code changed to convert COMMON/ACDPTH/ into ACDPTH_MOD. 57 USE ACDPTH_MOD 58C 10/14/26 MHP via Claude 59C code changed to convert COMMON/ENVSTRUCT/ into ENVSTRUCT_MOD. 60 USE ENVSTRUCT_MOD 61C 10/14/26 MHP via Claude 62C code changed to convert COMMON/COMP/ into COMP_MOD. 63 USE COMP_MOD 64C 10/14/26 MHP via Claude 65C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD. 66 USE NMLPAR_MOD 67C 10/14/26 MHP via Claude 68C code changed to convert COMMON/CONSTANTS/ into CONSTANTS_MOD. 69 USE CONSTANTS_MOD 70C 10/14/26 MHP via Claude 71C code changed to convert COMMON/IOUNITS/ into IOUNITS_MOD. 72 USE IOUNITS_MOD 73C 10/20/26 MHP via Claude 74C code changed to convert COMMON/EOS06MISC/ into 75C EOS06MISC_MOD. 76 USE EOS06MISC_MOD 77C 10/20/26 MHP via Claude 78C code changed to replace implicit typing with IMPLICIT NONE and 79C explicit declarations. 80 IMPLICIT NONE 81 REAL*8 AMUJ,ATEFFL,BET14J,BETAIJ,BETAJ,BL,CHRH,CHT,CMIXL,CVJ,DAGE, 82 * * DELAJ,EMUJ,ETAJ,FP,FPLJ,FT,FTLJ,FXIONJ,HCOMP,HD,HP,HR, 83 * * HSTOT,HT,QAPJ,QATJ,QCPJ,QCPPJ,QCPTJ,QDPJ,QDTJ,QDTPJ,QDTTJ, 84 * * RCZ,RCZL,RMUJ,STD,STDL,STP,STPL,STT,STTL,STX,STZ,TEFFL 85 INTEGER IORDER,IRAD,JENV,KHI,KLO,M,MV,MX,NR,NT 86 LOGICAL*4 LATMOJ,LCALCJ,LDERIJ 87C 10/9/26 MHP via Claude 88C code changed to add mx,mv,nr,nt (matching every other file that 89C declares an OPAL 2006 EOS common block), needed below for 90C COMMON/EOS06MISC/'s array dimensions; this file previously only 91C needed COMMON/EEOS06/, which it hardcoded as eos(10) instead. 92 PARAMETER (mx=5,mv=10,nr=169,nt=197) 93 REAL*8 HOOD1(7), HOOD2(7), D(7), ANS(1) 94 REAL*8 DUM1(4),DUM2(3),DUM3(3),DUM4(3) 95 REAL*8 XJ,ZJ,GLJ,RLJ,PLIMJ, BJ 96 REAL*8 ABEG0,AMIN0,AMAX0,EBEG0,EMIN0,EMAX0 97C 10/15/26 MHP via Claude 98C code changed to drop TNORM, TCZ and TAUCZN from this REAL*8 statement. 99C All three are COMMON/ACDPTH/ members, declared at line 99 below, and 100C explicitly typing a COMMON member is legal while doing the same to a 101C USE-associated name is rejected, so they had to go before ACDPTH 102C becomes a module. Removing them changes nothing: IMPLICIT 103C REAL*8(A-H,O-Z) already gives all three that type, T being in O-Z. 104 REAL*8 STARR(JSON),START(JSON),STARD(JSON),STARX(JSON),STARP(JSON), 105 * * STARZ(JSON),STARC(JSON),STARCZ(JSON),STARRZ(JSON) 106C REAL*8 ATMOSX, ATMOST(JSON),ATMOSD(JSON),ATMOSP(JSON),ADELAD(JSON),ATMOSC(JSON), 107C * ATGAM1(JSON), ATMOSR(JSON),GM1(JSON) 108 REAL*8 ATMOSX, GM1(JSON) 109C KC 2025-05-30 addressed warning messages from Makefile.legacy 110C INTEGER I, K, N, U, ICZ, NN, IENDJ, FLG, V, IATCNT, KTSAV 111 INTEGER I, K, P, U, ICZ, NN, IENDJ, FLG, V ! IATCNT, KTSAV 112 INTEGER IDUMJ, KATMJ,KENVJ,KSAHAJ,IXXJ 113 DIMENSION RCZL(1), RCZ(1),HR(JSON), HD(JSON), HP(JSON), BJ(1),FPLJ(1),FTLJ(1), 114 * * FT(JSON),FP(JSON), HCOMP(15,JSON), HT(JSON) ! HL(JSON) 115 LOGICAL LPULPTJ,LPRTJ,LSBC0J ! LTEST 116 117C KC 2025-05-30 reordered common block elements 118C COMMON/ACDPTH/TAUCZN,DELADJ(JSON),TAUHE, TNORM, TCZ, WHE, ICLCD, 119C 9/25/26 MHP via Claude 120C code changed to consolidate all I/O logical unit numbers, previously 121C scattered across COMMON/LUOUT/, COMMON/LUNUM/, COMMON/ALEX06/, 122C COMMON/ACDPTH/, COMMON/NEWOPAC/, COMMON/LOPAL95/, COMMON/ATMOS2/, 123C COMMON/ALATM03/, COMMON/OPALEOS/, and COMMON/SCV2/, into a single 124C COMMON/IOUNITS/ 125C 9/25/26 MHP via Claude 126C code changed to remove dead IACAT, IJVS, IJENT, IJDEL unit 127C reservations from COMMON/ACDPTH/ (declared but never opened anywhere) 128C *LCLCD, AGEOUT(5), IACAT, IJLAST, LJLAST, LJWRT, LADON, LAOLY, IJVS, 129C 10/12/26 MHP via Claude 130C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 131C namelist variables that the physics routines read from 27 COMMON 132C blocks which also held derived or working members. Those blocks keep 133C only their non-namelist members; five of them -- NEWMX, VMULT, 134C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 135C Derivation-only namelist variables from the same blocks went to 136C COMMON/NMLBYPASS/ instead, not here. 137C 138C Declared identically in all 121 files that hold any of these members. 139C Generated from one member list, so every copy is byte-identical; the 140C linker's 'size of symbol changed' warning is the check that they are. 141C 9/25/26 MHP via Claude 142C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/, 143C and COMMON/CONST3/ into a single COMMON/CONSTANTS/ 144C 10/13/26 MHP via Claude: removed the unused COMMON/SCRTCH/ declaration 145C (no member referenced here; residue of the obsolete-option removals). 146C 10/13/26 MHP via Claude: removed the unused COMMON/SOUND/ declaration 147C (no member referenced here; residue of the obsolete-option removals). 148C 10/9/26 MHP via Claude 149C code changed to merge COMMON/EEOS06/, COMMON/AAEOS06/, COMMON/BBEOS06/, 150C and COMMON/BEOS06/ into COMMON/EOS06MISC/ (see esac06.f for the full 151C rationale). 152C 10/11/26 MHP via Claude 153C code changed to remove an incorrect COMMON /LLOT95/ declaration from 154C this file. Every other file declares it as 155C COMMON /LLOT95/CAPPA(NUMX,NUMT,NUMD),ZTAB 156C so declaring only ZTAB here aliased it to CAPPA(1,1,1) -- an opacity 157C value, not the table's Z. Harmless to date only because this file's 158C two references to ZTAB are both in commented-out CALL ESAC06 lines 159C (see below); nothing reads it. NOTE: if those calls are ever 160C restored they need the real metallicity, which must come from the 161C full block declaration or (better) be passed in as an argument. 162C MHP 10/24 ADDED STOP CRITERIA FOR CENTRAL H,D,AND HE4 163C 10/13/26 MHP via Claude: removed the unused COMMON/SETT/ declaration 164C (no member referenced here; residue of the obsolete-option removals). 165C COMMON/SETT/ENDAGE(50),SETDT(50),LENDAG(50),LSETDT(50) 166C 10/13/26 MHP via Claude: removed the unused COMMON/ATMOS/ declaration 167C (no member referenced here; residue of the obsolete-option removals). 168C 10/13/26 MHP via Claude: removed the unused COMMON/SCVEOS/ declaration 169C (no member referenced here; residue of the obsolete-option removals). 170C MHP 8/25 Remove file names from common blocks 171C 9/24/26 MHP via Claude 172C code changed to remove the OPAL 1995/2001 equations of state 173C Added for experimantal mass loss section 174C G Somers 10/14, Add spot common block 175C G Somers END 176 SAVE 177 178 179C Initialize values: 180 TAUCZN=0.0D0 181 TCZ= 0.0D0 182 TNORM=0.0D0 183 WHE= 0.0D0 184 TAUHE = 0.0D0 185 ATMOSX=0.0D0 186 TATMOS = 0.0D0 187C KTSAV = 3 188 189C Need to calculate the sound speed in the envelope: first, stitch 190C envelope and interior together (so that T, P, R all agree at fitting 191C point. Then calculate sound speed in the envelope and interior 192C and integrate. 193 194C First, do the envelope integration to make sure the values in ENVSTRUCT 195C are up-to-date (this section of code grabbed from wrtmod). First, integrate the atm and store values with grey atm 196 ABEG0 = ATMBEG 197 AMIN0 = ATMMIN 198 AMAX0 = ATMMAX 199 EBEG0 = ENVBEG 200 EMIN0 = ENVMIN 201 EMAX0 = ENVMAX 202 ATMBEG = ATMSTP 203 ATMMIN = ATMSTP 204 ATMMAX = ATMSTP 205 ENVBEG = ENVSTP 206 ENVMIN = ENVSTP 207 ENVMAX = ENVSTP 208 IDUMJ = 0 209 BJ = DEXP(CLN*BL) 210 FPLJ(1) = FP(M) 211 FTLJ(1) = FT(M) 212 KATMJ = 0 213 KENVJ = 0 214 KSAHAJ = 0 215 LPULPTJ=.TRUE. 216 IXXJ=0 217 LPRTJ=.TRUE. 218 LSBC0J = .FALSE. 219 XJ = HCOMP(1,M) 220 ZJ = HCOMP(3,M) 221 RLJ = 0.5D0*(BL + CLSUNL - 4.0D0*TEFFL - C4PIL - CSIGL) 222 GLJ = CGL + STOTAL - RLJ - RLJ 223 PLIMJ = HP(M) 224c IF (KTTAU .EQ. 0) LAOLY = .TRUE. 225c IF (KTTAU .EQ. 3) LAOLY = .FALSE. ! for grey atm intergration: stores values in common block 226C G Somers 10/14, FOR SPOTTED RUNS, FIND THE PRESSURE AT 227C THE AMBIENT TEMPERATURE ATEFFL 228 IF(JENV.EQ.M.AND.SPOTF.NE.0.0.AND.SPOTX.NE.1.0)THEN 229 ATEFFL = TEFFL - 0.25*LOG10(SPOTF * SPOTX**4.0 + 1.0 - SPOTF) 230 ELSE 231 ATEFFL = TEFFL 232 ENDIF 233 CALL ENVINT(BJ,FPLJ,FTLJ,GLJ,HSTOT,IXXJ,LPRTJ,LSBC0J, 234 * * PLIMJ,RLJ,ATEFFL,XJ,ZJ,DUM1,IDUMJ,KATMJ,KENVJ,KSAHAJ, 235 * * DUM2,DUM3,DUM4,LPULPTJ,CMIXL) 236C G Somers END 237 238 239 240C Stitch interior and envelope together, convert into ESAC06 units 241 DO 34, I=1,M 242 STARR(I)=DEXP(CLN*HR(I)) ! get unlogged units 243 START(I)=DEXP(CLN*HT(I))/1.0D6 ! ESAC06 takes T in units of 10^6K 244 STARD(I)=DEXP(CLN*HD(I)) 245 STARP(I)=DEXP(CLN*HP(I)) 246 STARX(I)=HCOMP(1,I) 247 STARZ(I)=HCOMP(3,M) 24834 34 CONTINUE 249 DO 35, I=1,NUMENV-1 250 STARR(M+I)=DEXP(CLN*ENVR(I+1)) 251 START(M+I)=DEXP(CLN*ENVT(I+1))/1.0D6 252 STARD(M+I)=DEXP(CLN*ENVD(I+1)) 253 STARP(M+I)=DEXP(CLN*ENVP(I+1)) 254 STARX(M+I)=HCOMP(1,M) 255 STARZ(M+I)=HCOMP(3,M) 25635 35 CONTINUE 257 258 259C Call EOS interpolator: Takes OPAL2006 EOS if possible. If OPAL06 is off, uses SCV 260C For OPAL 2006 EOS: 261 IF (LOPALE06) THEN 262 IORDER=9 263 IRAD=1 264 DO 36, I=1,M+NUMENV-1 265 STX=STARX(I) 266 STT=START(I) 267 STD=STARD(I) 268C CALL ESAC06(STX,ZTAB,STT,STD,IORDER,IRAD,*999) ! KC 2025-05-31 269 CALL ESAC06(STX,STT,STD,IORDER,IRAD,*999) 270999 999 CONTINUE 271 STARC(I)=1.0D0/SQRT(eos(8)*STARP(I)/STARD(I)) 272 DELADJ(I)=1.0D0/eos(9) 273 GM1(I) = eos(8) 27436 36 CONTINUE 275 ENDIF 276C For SCV EOS: (when OPAL06 is turned off, SCV on for backup, Yale if SCV is off) 277 IF ( .NOT. LOPALE06) THEN 278 DO 41, I=1,M+NUMENV-1 279 STX=STARX(I) 280 STT=START(I)*1.0D6 281 STD=STARD(I) 282 STDL=DLOG10(STARD(I)) 283 STZ=STARZ(I) 284 STTL=DLOG10(STT) 285 STP=STARP(I) 286 STPL=DLOG10(STP) 287 LCALCJ = .TRUE. 288C CALL EQSCVE(STTL,STT,STPL,STP,STDL,STD,STX,STZ,BETAJ,BETAIJ,BET14J,FXIONJ,RMUJ, 289C * AMUJ,EMUJ,ETAJ,QDTJ,QDPJ,QCPJ,DELAJ,LCALCJ) 290 LDERIJ = .FALSE. 291 LATMOJ = .TRUE. 292 CALL EQSTAT2(STTL,STT,STPL,STP,STDL,STD,STX,STZ,BETAJ,BETAIJ,BET14J,FXIONJ, 293 * * RMUJ,AMUJ,EMUJ,ETAJ,QDTJ,QDPJ,QCPJ,DELAJ,QDTTJ,QDTPJ,QATJ,QAPJ,QCPTJ,QCPPJ, 294 * * LDERIJ,LATMOJ,KSAHAJ) 295 296C calculate gamma1 (taken from envint example, eqns in chap. 9 of Cox& Giuli): 297 CHRH = 1.0D0/QDPJ 298 CHT = -CHRH*QDTJ 299 CVJ = QCPJ - EXP(CLN*(STPL-STDL-STTL))*CHT**2/CHRH 300 GM1(I) = CHRH*QCPJ/CVJ 301 STARC(I)=1.0D0/SQRT(GM1(I)*STARP(I)/STARD(I)) 302 DELADJ(I) = DELAJ 30341 41 CONTINUE 304 ENDIF 305 306 307 308 309C Integrate for the full acoustic depth for normalization 310 K=NUMENV-1+M-1 311 U=MOD((NUMENV-1+M-1),4) 312 IF (U.NE.0) THEN 313 NN=K+4-U+1 314 ELSE 315 NN=K+1 316 ENDIF 317 IF (K.LE.3) THEN 318 K=4 319 NN=1+1 320 ENDIF 321 322 CALL BOOLE(STARR,STARC,K,NN,TNORM) 323 324 325 326 327C Find the location of the CZ, but only if one exists 328 IF (JENV.LT.M) THEN 329 RCZ=DEXP(CLN*(RCZL+CRSUNL)) 330 331 332C get the location of the convection zone, and count nonzero entries: 333 IENDJ=0 334 FLG=0 335 DO 22, I=1,K 336 IF (STARR(I).LE.RCZ(1) .AND. STARR(I+1).GE.RCZ(1)) THEN 337 ICZ = I 338 ENDIF 339 340C Interpolate to find Cs at the exact CZ radius 341 IF (I.EQ.ICZ .AND. FLG.EQ.0 .AND. ICZ.GT.3) THEN 342 FLG=1 343 DO 26 P=1,7 344 HOOD1(P)=STARR(ICZ-4+P) 345 HOOD2(P)=STARC(ICZ-4+P) 34626 26 CONTINUE 347 CALL SPLINJ(HOOD1,HOOD2,D,7) 348 CALL SPLINT(HOOD1,HOOD2,7,D,RCZ,ANS,KLO,KHI) 349 STARR(ICZ)=RCZ(1) 350 STARC(ICZ)=ANS(1) 351 ENDIF 352 IF (I.EQ.ICZ .AND. FLG.EQ.0 .AND. ICZ.LE.3) THEN 353 FLG=1 354 IF (ICZ .EQ. 1) THEN 355 RCZ(1) = STARR(ICZ) 356 ANS(1) = STARC(ICZ) 357 ELSE 358 DO 66 P=1,7 359 HOOD1(P)=STARR(ICZ-2+P) 360 HOOD2(P)=STARC(ICZ-2+P) 36166 66 CONTINUE 362 CALL SPLINJ(HOOD1,HOOD2,D,7) 363 CALL SPLINT(HOOD1,HOOD2,7,D,RCZ,ANS,KLO,KHI) 364 STARR(ICZ)=RCZ(1) 365 STARC(ICZ)=ANS(1) 366 ENDIF 367 ENDIF 368 36922 22 CONTINUE 370 DO 27, I=1,K-ICZ+1 371 STARRZ(I)=STARR(ICZ-1+I) 372 STARCZ(I)=STARC(ICZ-1+I) 37327 27 CONTINUE 374C Then call Boole from surface to cz: 375 U=MOD((K-ICZ+1),4) 376 V=K-ICZ+1 377 IF (U.NE.0) THEN 378 NN=K-ICZ+1+4-U+1 379 ELSE 380 NN=K-ICZ+1+1 381 ENDIF 382 CALL BOOLE(STARRZ,STARCZ,V,NN,TCZ) 383 ELSE 384 TCZ = 0.0d0 385 ENDIF 386 387C Deal with the atmosphere seperately 388C Invert vectors and shift units 389C Adopt the atmosphere metallicity to be the same as the last envelope pt. 390C ATMOSX=STARX(M+NUMENV-1) 391C Determine length of Atmosphere vector: 392C IATCNT=0 393C DO 56, I=1,JSON-1 394C IF(ACATMR(I) .NE. 0.0 .AND. ACATMR(I+1) .LE. 1.0D-7) THEN 395C IATCNT=I 396C GOTO 56 397C END IF 398C56 CONTINUE 399 400 401C Put the atmosphere in actual units, not differential 402C DO 39, I=1,IATCNT-1 403C ACATMR(IATCNT-I)=ACATMR(IATCNT-I+1)+ACATMR(IATCNT-I) 404C39 CONTINUE 405C DO 40, I=1,IATCNT 406C ACATMR(I)=ACATMR(I)+STARR(M+NUMENV-1) 407C40 CONTINUE 408 409C DO 37,I=1,IATCNT 410C ATMOSD(I)=DEXP(CLN*ACATMD(IATCNT+1-I)) 411C ATMOST(I)=DEXP(CLN*ACATMT(IATCNT+1-I))/1.0D6 412C ATMOSP(I)=DEXP(CLN*ACATMP(IATCNT+1-I)) 413C ATMOSR(I)=ACATMR(IATCNT+1-I) 414C37 CONTINUE 415 416C Pass atm vectors to ESAC06 417C IORDER=9 418C IRAD=1 419C DO 38, I=1,IATCNT 420C CALL ESAC06(ATMOSX,ZTAB,ATMOST(I),ATMOSD(I),IORDER,IRAD,*998) 421C998 CONTINUE 422C ATMOSC(I)=1.0D0/SQRT(eos(8)*ATMOSP(I)/ATMOSD(I)) 423C ADELAD(I)=1.0D0/eos(9) 424C ATGAM1(I)= eos(8) 425C38 CONTINUE 426 427 428C Then call Boole to do the integration: 429C U=MOD((IATCNT),4) 430c V=IATCNT 431c IF (U.NE.0) THEN 432c NN=V+4-U+1 433c ELSE 434c NN=V+1 435c ENDIF 436c 437c CALL BOOLE(ATMOSR,ATMOSC,V,NN,TATMOS) 438 439C Output normalized acoustic depth 440 TAUCZN=0.0D0 441 TAUCZN=TCZ/TNORM 442 443C Output acoustic depth info to ISHORT 444 WRITE(ISHORT,67)TCZ,TNORM,TAUCZN,TATMOS 44567 67 FORMAT(1X,'Acoustic depth to CZ:',F14.8,2X,'Acoustic depth to center', 446 * *F13.7,2X,'Normalized taucz:',F11.9, 'Acoustic depth of atmopshere:',F16.8,2X) 447 448C 555 CONTINUE 449C-------------------------------------------------------------- 450C Save all vectors of interest when the end of a kind card is reached. 451 IF(LCLCD)THEN 452 453 WRITE(UNIT=ICLCD,FMT=1506) NUMENV, M, ICZ 4541506 1506 FORMAT(1X, 'Number of points in envelope:',I5,2X, 455 * * 'Number of points in interior:',I5,2X,'Index near Rcz:' 456 * * ,I5,2X) 457 458 DO 1505 I=1,M+NUMENV-1 459 IF (I .LE. K-ICZ+1) THEN 460C KC 2025-05-30 addressed warning messages from Makefile.legacy 461C WRITE(UNIT=ICLCD,FMT=1504),DAGE, STARR(I), STARC(I), 462 WRITE(UNIT=ICLCD,FMT=1504) DAGE, STARR(I), STARC(I), 463 * * STARRZ(I), STARCZ(I), DELADJ(I), GM1(I), STARP(I), 464 * * START(I),STARD(I),STARX(I) 465 466 ELSE 467C KC 2025-05-30 addressed warning messages from Makefile.legacy 468C WRITE(UNIT=ICLCD,FMT=1504) DAGE, STARR(I), STARC(I), 469 WRITE(UNIT=ICLCD,FMT=1504) DAGE, STARR(I), STARC(I), 470 * * 0.0d0, 0.0d0, DELADJ(I), GM1(I), STARP(I), 471 * * START(I),STARD(I),STARX(I) 472 473 ENDIF 474 475 4761505 1505 CONTINUE 4771504 1504 FORMAT(1X,11E16.8) 478c DO 1520 I=1,IATCNT 479c WRITE(UNIT=IACAT,FMT=1521),DAGE, ATMOSR(I), ATMOSC(I), 480c * DELADJ(I), ATGAM1(I), ATMOSP(I), ATMOST(I), ATMOSX 481c 482c1520 CONTINUE 483c1521 FORMAT(1X,8E16.8) 484 485 486 ENDIF 487 488 489 490c DO 444, I=1,JSON 491c ACATMR(I) =0.0d0 492c ACATMP(I) =0.0d0 493c ACATMT(I) =0.0d0 494c ACATMD(I) =0.0d0 495c444 CONTINUE 496 497 498 499 500 501 502 503 504 505 506C 50 CONTINUE 507 RETURN 508 509 END 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542

1"""Python version of calcad.f (YREC 7be4133): acoustic depths of the star 2and of its convection zone.""" 3 4import math 5 6import numpy as np 7 8from pyrec import acdpth_mod 9from pyrec import comp_mod 10from pyrec import constants_mod 11from pyrec import dimensions_mod 12from pyrec import envstruct_mod 13from pyrec import eos06misc_mod 14from pyrec import iounits_mod 15from pyrec import nmlpar_mod 16from pyrec.boole import boole 17from pyrec.envint import envint 18from pyrec.eqstat2 import eqstat2 19from pyrec.esac06 import esac06 20from pyrec.splinj import splinj 21from pyrec.splint import splint 22 23JSON = dimensions_mod.json 24 25# CALCAD has SAVE. Kept between calls because later calls can read them: 26# ICZ (when no pair of points brackets the base of the convection zone, or 27# when JENV >= M) and the convection-zone profile written with LCLCD. 28icz = 0 29starrz = np.zeros(JSON + 1) 30starcz = np.zeros(JSON + 1) 31# The step sizes saved on entry and never restored. Suspected bug (see the 32# page). 33abeg0 = amin0 = amax0 = ebeg0 = emin0 = emax0 = 0.0 34 35 36def calcad(hr: np.ndarray, rczl: float, m: int, hd: np.ndarray, hp: np.ndarray, ht: np.ndarray, bl: float, 37 fp: np.ndarray, ft: np.ndarray, hstot: float, teffl: float, hcomp: np.ndarray, dage: float, jenv: int, 38 cmixl: float) -> None: 39 """Join the interior model and a fresh envelope integration into one 40 profile, compute 1/c_s from the equation of state, and integrate it to 41 get the acoustic depth of the star (TNORM) and of its convection zone 42 (TCZ), and their ratio TAUCZN. 43 44 Translates: calcad.f#CALCAD (7be4133) 45 46 Arguments: the interior model hr, hd, hp, ht (log r, log rho, log P, 47 log T; arrays, index 0 unused); rczl, log r_CZ / R_sun (RCZL(1) in the 48 Fortran, a scalar here); m, the last interior point; bl, log L; fp, ft, 49 the rotational distortion; hstot; teffl, log Teff; hcomp[species, 50 shell]; dage, the age; jenv, the base of the convection zone; cmixl, 51 the mixing length (passed on to ENVINT). 52 53 Returns nothing: the results go to acdpth_mod (tnorm, tcz, tauczn, 54 tatmos, deladj). 55 56 """ 57 58 global icz, abeg0, amin0, amax0, ebeg0, emin0, emax0 59 c = constants_mod 60 nml = nmlpar_mod 61 acd = acdpth_mod 62 env = envstruct_mod 63 eos = eos06misc_mod.eos 64 65 starr = np.zeros(JSON + 1) 66 start = np.zeros(JSON + 1) 67 stard = np.zeros(JSON + 1) 68 starx = np.zeros(JSON + 1) 69 starp = np.zeros(JSON + 1) 70 starz = np.zeros(JSON + 1) 71 starc = np.zeros(JSON + 1) 72 gm1 = np.zeros(JSON + 1) 73 hood1 = np.zeros(8) 74 hood2 = np.zeros(8) 75 d = np.zeros(8) 76 dum1 = np.zeros(5) 77 dum2 = np.zeros(4) 78 dum3 = np.zeros(4) 79 dum4 = np.zeros(4) 80 fxionj = np.zeros(4) 81 82 # Initialize the results (Fortran lines 180–186). 83 acd.tauczn = 0.0 84 acd.tcz = 0.0 85 acd.tnorm = 0.0 86 acd.whe = 0.0 87 acd.tauhe = 0.0 88 acd.tatmos = 0.0 89 90 # The envelope, integrated with fixed step sizes (Fortran lines 91 # 196–207). Suspected bug (see the page): the step sizes are never 92 # restored. 93 abeg0, amin0, amax0 = nml.atmbeg, nml.atmmin, nml.atmmax 94 ebeg0, emin0, emax0 = nml.envbeg, nml.envmin, nml.envmax 95 nml.atmbeg = nml.atmmin = nml.atmmax = nml.atmstp 96 nml.envbeg = nml.envmin = nml.envmax = nml.envstp 97 98 # The arguments of ENVINT (Fortran lines 208–223). 99 idumj = 0 100 bj = math.exp(c.cln * bl) 101 fplj = fp[m] 102 ftlj = ft[m] 103 katmj = 0 104 kenvj = 0 105 ksahaj = 0 106 lpulptj = True 107 ixxj = 0 108 lprtj = True 109 lsbc0j = False 110 xj = hcomp[1, m] 111 zj = hcomp[3, m] 112 rlj = 0.5 * (bl + c.clsunl - 4.0 * teffl - c.c4pil - c.csigl) 113 glj = c.cgl + comp_mod.stotal - rlj - rlj 114 plimj = hp[m] 115 116 # For a spotted star, the ambient temperature (Fortran lines 228–232). 117 # Question (see the page): spots only when JENV = M. 118 if jenv == m and nml.spotf != 0.0 and nml.spotx != 1.0: 119 ateffl = teffl - 0.25 * math.log10( 120 nml.spotf * nml.spotx**4.0 + 1.0 - nml.spotf) 121 else: 122 ateffl = teffl 123 124 # The envelope integration, which fills envstruct_mod (Fortran lines 125 # 233–235). 126 rlj, idumj, katmj, kenvj, ksahaj = envint( 127 bj, fplj, ftlj, glj, hstot, ixxj, lprtj, lsbc0j, plimj, rlj, ateffl, 128 xj, zj, dum1, idumj, katmj, kenvj, ksahaj, dum2, dum3, dum4, 129 lpulptj, cmixl) 130 numenv = env.numenv 131 132 # One profile: interior points 1..M, then envelope points 2..NUMENV, 133 # in linear units, T in 10^6 K (labels 34, 35; Fortran lines 241–256). 134 for i in range(1, m + 1): 135 starr[i] = math.exp(c.cln * hr[i]) 136 start[i] = math.exp(c.cln * ht[i]) / 1.0e6 137 stard[i] = math.exp(c.cln * hd[i]) 138 starp[i] = math.exp(c.cln * hp[i]) 139 starx[i] = hcomp[1, i] 140 # Suspected bug (see the page): the envelope's Z at every interior 141 # point. 142 starz[i] = hcomp[3, m] 143 for i in range(1, numenv): 144 starr[m + i] = math.exp(c.cln * env.envr[i + 1]) 145 start[m + i] = math.exp(c.cln * env.envt[i + 1]) / 1.0e6 146 stard[m + i] = math.exp(c.cln * env.envd[i + 1]) 147 starp[m + i] = math.exp(c.cln * env.envp[i + 1]) 148 starx[m + i] = hcomp[1, m] 149 starz[m + i] = hcomp[3, m] 150 151 # Inverse sound speed and the adiabatic gradient, from the OPAL 2006 152 # EOS (label 36; Fortran lines 261–275). 153 if nml.lopale06: 154 iorder = 9 155 irad = 1 156 for i in range(1, m + numenv): 157 # The Fortran's alternate return (*999) lands on the next 158 # statement, so a failed call goes on with the previous EOS. 159 # Suspected bug (see the page): EOS errors ignored. 160 esac06(starx[i], start[i], stard[i], iorder, irad) 161 starc[i] = 1.0 / math.sqrt(eos[8] * starp[i] / stard[i]) 162 acd.deladj[i] = 1.0 / eos[9] 163 gm1[i] = eos[8] 164 165 # Otherwise from EQSTAT2, with Gamma_1 from Cox and Giuli, chapter 9 166 # (label 41; Fortran lines 277–304). 167 if not nml.lopale06: 168 for i in range(1, m + numenv): 169 stx = starx[i] 170 stt = start[i] * 1.0e6 171 std = stard[i] 172 stdl = math.log10(stard[i]) 173 stz = starz[i] 174 sttl = math.log10(stt) 175 stp = starp[i] 176 stpl = math.log10(stp) 177 lcalcj = True # noqa: F841 (unused in the Fortran too) 178 lderij = False 179 latmoj = True 180 # Suspected bug (see the page): FXIONJ is a scalar in CALCAD but 181 # FXION(3) in EQSTAT2, which writes past it; here an array, so 182 # that EQSTAT2 can fill it. 183 (stt, stp, stdl, std, betaj, betaij, bet14j, rmuj, amuj, emuj, 184 etaj, qdtj, qdpj, qcpj, delaj, qdttj, qdtpj, qatj, qapj, qcptj, 185 qcppj, ksahaj) = eqstat2( 186 sttl, stt, stpl, stp, stdl, std, stx, stz, 0.0, 0.0, 0.0, 187 fxionj, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 188 0.0, 0.0, 0.0, 0.0, lderij, latmoj, ksahaj) 189 chrh = 1.0 / qdpj 190 cht = -chrh * qdtj 191 cvj = qcpj - math.exp(c.cln * (stpl - stdl - sttl)) * cht**2 / chrh 192 gm1[i] = chrh * qcpj / cvj 193 starc[i] = 1.0 / math.sqrt(gm1[i] * starp[i] / stard[i]) 194 acd.deladj[i] = delaj 195 196 # t_norm: Boole's rule on 1 + 4p points (Fortran lines 310–322). 197 # Question (see the page): K leaves out the outermost point. 198 k = numenv - 1 + m - 1 199 u = (numenv - 1 + m - 1) % 4 200 if u != 0: 201 nn = k + 4 - u + 1 202 else: 203 nn = k + 1 204 # Suspected bug (see the page): a short profile gives TNORM = 0. 205 if k <= 3: 206 k = 4 207 nn = 1 + 1 208 acd.tnorm = boole(starr, starc, k, nn, acd.tnorm) 209 210 # t_CZ: from r_CZ to the surface, if there is a convection zone 211 # (Fortran lines 328–385). 212 if jenv < m: 213 rcz = math.exp(c.cln * (rczl + c.crsunl)) 214 ans = 0.0 215 klo = khi = 0 216 flg = 0 217 218 # Find the point below r_CZ and put r_CZ there, with 1/c_s from a 219 # spline through 7 points (loop 22; Fortran lines 333–369). 220 for i in range(1, k + 1): 221 if starr[i] <= rcz and starr[i + 1] >= rcz: 222 icz = i 223 # Suspected bug (see the page): ICZ may be stale, from an 224 # earlier call, and then can match I before the bracket is 225 # found. 226 if i == icz and flg == 0 and icz > 3: 227 flg = 1 228 for p in range(1, 8): 229 hood1[p] = starr[icz - 4 + p] 230 hood2[p] = starc[icz - 4 + p] 231 splinj(hood1, hood2, d, 7) 232 ans, klo, khi = splint(hood1, hood2, 7, d, rcz, ans, klo, khi) 233 starr[icz] = rcz 234 starc[icz] = ans 235 if i == icz and flg == 0 and icz <= 3: 236 flg = 1 237 if icz == 1: 238 rcz = starr[icz] 239 ans = starc[icz] 240 else: 241 for p in range(1, 8): 242 hood1[p] = starr[icz - 2 + p] 243 hood2[p] = starc[icz - 2 + p] 244 splinj(hood1, hood2, d, 7) 245 ans, klo, khi = splint(hood1, hood2, 7, d, rcz, ans, 246 klo, khi) 247 starr[icz] = rcz 248 starc[icz] = ans 249 250 # The profile from r_CZ outward (loop 27; Fortran lines 370–373). 251 for i in range(1, k - icz + 2): # label 27 252 starrz[i] = starr[icz - 1 + i] 253 starcz[i] = starc[icz - 1 + i] 254 255 # Boole's rule from r_CZ to the surface (Fortran lines 375–382). 256 u = (k - icz + 1) % 4 257 v = k - icz + 1 258 if u != 0: 259 nn = k - icz + 1 + 4 - u + 1 260 else: 261 nn = k - icz + 1 + 1 262 acd.tcz = boole(starrz, starcz, v, nn, acd.tcz) 263 else: 264 acd.tcz = 0.0 265 266 # The normalized acoustic depth (Fortran lines 440–441). Suspected bug 267 # (see the page): TNORM = 0 for a short profile; Python raises 268 # ZeroDivisionError where the Fortran gives Inf or NaN. 269 acd.tauczn = 0.0 270 acd.tauczn = acd.tcz / acd.tnorm 271 272 # The line to ISHORT (Fortran lines 444–446). Suspected bug (see the 273 # page): format 67. Question (see the page): TATMOS is always 0. 274 iounits_mod.write( 275 iounits_mod.ishort, 276 f" Acoustic depth to CZ:{acd.tcz:14.8f} Acoustic depth to center" 277 f"{acd.tnorm:13.7f} Normalized taucz:{acd.tauczn:11.9f}" 278 f"Acoustic depth of atmopshere:{acd.tatmos:16.8f} ") 279 280 # With LCLCD, the full profile to ICLCD (Fortran lines 451–486). 281 if acd.lclcd: 282 iounits_mod.write( 283 iounits_mod.iclcd, 284 f" Number of points in envelope:{numenv:5d} " 285 f"Number of points in interior:{m:5d} Index near Rcz:{icz:5d} ") 286 for i in range(1, m + numenv): # label 1505 287 if i <= k - icz + 1: 288 row = (dage, starr[i], starc[i], starrz[i], starcz[i], 289 acd.deladj[i], gm1[i], starp[i], start[i], stard[i], 290 starx[i]) 291 else: 292 row = (dage, starr[i], starc[i], 0.0, 0.0, acd.deladj[i], 293 gm1[i], starp[i], start[i], stard[i], starx[i]) 294 iounits_mod.write(iounits_mod.iclcd, 295 " " + "".join(f"{x:16.8E}" for x in row))

calcad.f

Purpose

How it fits

The method

Design

Suspected bugs, checked by running the code

1. (cannot be reached) Step sizes not restored.

2. (cannot be reached) The interior's $Z$.

3. (cannot be reached) ICZ may be stale.

4. (cannot be reached) FXIONJ is a scalar.

5. (cannot be reached) EOS errors ignored.

6. (cannot be reached) Output format.

7. (cannot be reached) A short profile.

8. (cannot be reached) The surface point is left out of the integrals.

Definitions

subroutineCALCAD

Change log

3. (cannot be reached) `ICZ` may be stale.

4. (cannot be reached) `FXIONJ` is a scalar.