momi.f

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

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

Variables on this page (6)

cc13, cc23, cgl, cln, json, walpcz

Purpose

Computes the moment of inertia HI of each shell of a rotating, distorted model, and its derivative with respect to $\ln A$, $A$ being the distortion parameter below, QIW $= dI/d\ln A$, which at fixed shape ($\eta_2$, $r_0$) is $\tfrac12\,\Omega\,\partial I/\partial\Omega$ (since $A \propto \Omega^2$; the comments name it a derivative in $\Omega$, line 50 ↗, and SOLID uses it as one, a suspected bug: solid.f, suspected bug 2). The comment: "Moment of inertia is calculated by the method used in Wai-Yuen Law's thesis (Yale, 1980) p.61," and "MOMI assumes that SHAPE has already been called."

How it fits

Inputs: the shape parameters ETA2 (Radau's $\eta_2$) and R0 (the characteristic radius $r_0$ of the equipotential, $r = r_0(1 - AP_2)$) per shell, from SHAPE; $\log r$ (HR), $\log m$ (HS), the shell masses HS2, the angular velocities OMEGA, and the shell range JSTART..JEND.
Who calls it: HPOINT, SOLID and STARIN.
What it calls: INTMOM, which evaluates the integral at each shell.
Outputs: HI and QIW per shell.
Modules: NMLPAR_MOD (WALPCZ), CONSTANTS_MOD (CLN, CC13, CC23, and $\log G$, CGL), DIMENSIONS_MOD (JSON). The unused /CCOUT/ and /LUOUT/ of rp_1 are no longer declared.

The method

Spherical shortcut. If WALPCZ $\ne 0$ (WALPCZ is the exponent of the rotation law $\Omega\, r^{\mathrm{WALPCZ}} =$ constant in the surface convection zone, as GETNEWENV's comment says; PARMIN keeps it between −2 and 0, so 0 is solid-body rotation), each shell is a thin spherical shell, $I = \tfrac23\,\Delta m\, r^2$, with $dI/d\Omega = 0$; the central shell is a sphere, $I = 0.4\,\Delta m\, r^2$, at the mean of the first two $\log r$ (lines 62–72 ↗).

Question 1 (lines 62–72 ↗). WALPCZ as a switch. The spherical shortcut is taken for every shell in JSTART..JEND, not only those in the convection zone, whenever the convection zone's rotation law is not solid-body (WALPCZ $\ne 0$). Why should that law switch off the distortion of all shells? - Distorted shells. Otherwise, for each shell, with $r_0$ = R0, the parameter $$A = \frac{5}{3}\,\frac{\Omega^2\, r_0^3}{G\,m\,(2 + \eta_2)}$$ (CON $=5/(3G)$), the ratio R0R $= r_0^5/r^3$, and QR0R $= d\ln r_0/d\ln r$ set to 1, go to INTMOM, which returns $I$ and QIWM per unit mass; these are multiplied by HS2. INTMOM sums a series to $A^5$ (Law 1980's coefficients $c_j$): $I = H_1\bigl(1 + 0.15\sum_j c_j(j\eta_2+5)A^j\bigr)$ and QIWM $= H_1\cdot0.15\sum_j j\,c_j(j\eta_2+5)A^j = A\,dI/dA$ at fixed $\eta_2$ and $r_0$ (lines 73–89 ↗). R0P and RPHIP (the previous shell's $\ln r_0$ and $\ln r$, lines 54–60 ↗, lines 87–88 ↗) and H0 (the thin-shell moment of inertia, line 77 ↗) are computed and never read (unused?), and QR0R is fixed at 1 (line 76 ↗).

Question 2 (lines 76–86 ↗; intmom.f, lines 26–31 ↗). A series to $A^5$ on a first-order shape. INTMOM expands $I$ to $A^5$, while SHAPE keeps the equipotentials to first order in $\Omega^2$ (only the $P_2$ term). Is the expansion the exact moment of a body $r = r_0(1-AP_2)$, kept for accuracy at large $A$, or should it be truncated consistently with the shape?

Definitions

subroutineMOMI

lines 3–92 · GitHub

Uses modules CONSTANTS_MOD, DIMENSIONS_MOD, NMLPAR_MOD

Called by HPOINT ×2, SOLID, STARIN

Flowchart of its calls, in order (1)

for I = JSTART,JEND (DO 20), lines 73–90

INTMOM 84

Compiler warnings (1)

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

line 62: Inequality comparison for REAL(8) [-Wcompare-reals]

Change log

The changes since rp_1

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

1C $$$$$$ 2C SUBROUTINE MOMI(ETA2,HD,HR,HS,HS2,JSTART,JEND,OMEGA,R0,HI,QIW,M) ! KC 2025-05-31 3 SUBROUTINE MOMI(ETA2,HR,HS,HS2,JSTART,JEND,OMEGA,R0,HI,QIW) 4C 10/15/26 MHP via Claude 5C code changed to take JSON from DIMENSIONS_MOD instead of declaring 6C PARAMETER(JSON=5000) here. 7 USE DIMENSIONS_MOD 8C 10/14/26 MHP via Claude 9C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD. 10 USE NMLPAR_MOD 11C 10/14/26 MHP via Claude 12C code changed to convert COMMON/CONSTANTS/ into CONSTANTS_MOD. 13 USE CONSTANTS_MOD 14C 10/20/26 MHP via Claude 15C code changed to replace implicit typing with IMPLICIT NONE and 16C explicit declarations. 17 IMPLICIT NONE 18 REAL*8 A,CON,ETA2,ETA2X,H0,HI,HIM,HR,HS,HS2,OMEGA,QIW,QIWM,QR0R, 19 * * R0,R03,R0P,R0R,RAD,RPHI3,RPHIP 20 INTEGER I,JEND,JSTART 21 22C 10/11/26 MHP via Claude 23C code changed to remove LSCRIB from COMMON/CCOUT/: it was declared 24C positionally in all 17 files sharing this block but never read or 25C written in executable code anywhere, and no namelist file sets it. 26C 10/12/26 MHP via Claude 27C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 28C namelist variables that the physics routines read from 27 COMMON 29C blocks which also held derived or working members. Those blocks keep 30C only their non-namelist members; five of them -- NEWMX, VMULT, 31C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 32C Derivation-only namelist variables from the same blocks went to 33C COMMON/NMLBYPASS/ instead, not here. 34C 35C Declared identically in all 121 files that hold any of these members. 36C Generated from one member list, so every copy is byte-identical; the 37C linker's 'size of symbol changed' warning is the check that they are. 38C 10/13/26 MHP via Claude: removed the unused COMMON/CCOUT/ declaration 39C (no member referenced here; residue of the obsolete-option removals). 40C 9/25/26 MHP via Claude 41C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/, 42C and COMMON/CONST3/ into a single COMMON/CONSTANTS/ 43C 10/13/26 MHP via Claude: removed the unused COMMON/IOUNITS/ declaration 44C (no member referenced here; residue of the obsolete-option removals). 45C DIMENSION HD(JSON),HR(JSON),HS2(JSON),R0(JSON),ETA2(JSON),HI(JSON) ! KC 2025-05-31 46 DIMENSION HR(JSON),HS2(JSON),R0(JSON),ETA2(JSON),HI(JSON) 47 DIMENSION OMEGA(JSON),HS(JSON),QIW(JSON) 48 SAVE 49 50C FIND THE MOMENT OF INERTIA (HI) AND DI/D(OMEGA) (QIW). 51C MOMI ASSUMES THAT SHAPE HAS ALREADY BEEN CALLED. 52C MOMENT OF INERTIA IS CALCULATED BY THE METHOD USED IN WAI-YUEN LAW'S 53C THESIS(YALE,1980) P.61. 54 IF(JSTART.EQ.1) THEN 55 R0P = 0.0D0 56 RPHIP = 0.0D0 57 ELSE 58 R0P = DLOG(R0(JSTART - 1)) 59 RPHIP = CLN*HR(JSTART - 1) 60 ENDIF 61 CON = CC13*5.0D0/DEXP(CLN*CGL) 62 IF(WALPCZ.NE.0.0D0)THEN 63 DO I = JSTART,JEND 64 HI(I) = CC23*HS2(I)*EXP(CLN*2.0D0*HR(I)) 65 QIW(I) = 0.0D0 66 END DO 67 IF(JSTART.EQ.1)THEN 68 RAD = 0.5D0*(HR(1)+HR(2)) 69 HI(1) = 0.4D0*HS2(1)*EXP(CLN*2.0D0*RAD) 70 ENDIF 71 RETURN 72 ENDIF 73 DO 20 I = JSTART,JEND 74C QR0R = D LN R0/ D LN R 75C 76 QR0R = 1.0D0 77 H0 = CC23*HS2(I)*DEXP(CLN*2.0D0*HR(I)) 78 R03 = R0(I)**3 79 RPHI3 = DEXP(CLN*3.0D0*HR(I)) 80 R0R = (R03/RPHI3)*R0(I)**2 81 A = CON*OMEGA(I)**2*R03/(DEXP(CLN*HS(I))*(2.0D0+ETA2(I))) 82 ETA2X = ETA2(I) 83C EVALUATE INTEGRAL AT THE ZONE CENTER 84 CALL INTMOM(A,ETA2X,QR0R,R0R,HIM,QIWM) 85 HI(I) = HIM*HS2(I) 86 QIW(I) = QIWM*HS2(I) 87 R0P = DLOG(R0(I)) 88 RPHIP = CLN*HR(I) 89 20 20 CONTINUE 90 91 RETURN 92 END

1"""Python version of momi.f (YREC 7be4133): the moment of inertia of each 2shell of a rotating, distorted model, and its derivative in A (called 3dI/dOmega in the comments).""" 4 5import math 6import numpy as np 7 8from pyrec import constants_mod 9from pyrec import nmlpar_mod 10from pyrec.intmom import intmom 11 12 13def momi(eta2: np.ndarray, hr: np.ndarray, hs: np.ndarray, hs2: np.ndarray, jstart: int, jend: int, omega: np.ndarray, 14 r0: np.ndarray, hi: np.ndarray, qiw: np.ndarray) -> None: 15 """Compute the moment of inertia hi and A dI/dA qiw of the shells 16 jstart..jend (Law's thesis, Yale 1980, p. 61), from the shapes found by 17 SHAPE; thin spherical shells when WALPCZ is not 0. 18 19 Translates: momi.f#MOMI (7be4133) 20 21 Arguments: the shape parameters eta2 and r0; log r (hr), log m (hs), 22 the shell masses hs2 and the angular velocities omega; the range 23 jstart..jend. Arrays have index 0 unused. 24 25 Sets in place: hi, qiw. Returns nothing. 26 27 """ 28 29 c = constants_mod 30 31 # The previous shell's ln r0 and ln r, and the constant of A 32 # (Fortran lines 54–61). r0p and rphip are unused (in the Fortran too). 33 if jstart == 1: 34 r0p = 0.0 # noqa: F841 (unused in the Fortran too) 35 rphip = 0.0 # noqa: F841 (unused in the Fortran too) 36 else: 37 r0p = math.log(r0[jstart - 1]) # noqa: F841 (unused in the Fortran too) 38 rphip = c.cln * hr[jstart - 1] # noqa: F841 (unused in the Fortran too) 39 con = c.cc13 * 5.0 / math.exp(c.cln * c.cgl) 40 41 # Thin spherical shells when WALPCZ is not 0 (Fortran lines 62–72). 42 # Question (see the page): WALPCZ as a switch for every shell. 43 if nmlpar_mod.walpcz != 0.0: 44 for i in range(jstart, jend + 1): 45 hi[i] = c.cc23 * hs2[i] * math.exp(c.cln * 2.0 * hr[i]) 46 qiw[i] = 0.0 47 if jstart == 1: 48 rad = 0.5 * (hr[1] + hr[2]) 49 hi[1] = 0.4 * hs2[1] * math.exp(c.cln * 2.0 * rad) 50 return 51 52 # Distorted shells, through INTMOM (Fortran lines 73–89). 53 for i in range(jstart, jend + 1): # loop 20 54 # QR0R = d ln r0 / d ln r. 55 qr0r = 1.0 56 h0 = c.cc23 * hs2[i] * math.exp(c.cln * 2.0 * hr[i]) # noqa: F841 (unused in the Fortran too) 57 r03 = r0[i]**3 58 rphi3 = math.exp(c.cln * 3.0 * hr[i]) 59 r0r = (r03 / rphi3) * r0[i]**2 60 a = con * omega[i]**2 * r03 / (math.exp(c.cln * hs[i]) 61 * (2.0 + eta2[i])) 62 eta2x = eta2[i] 63 # The integral at the zone centre. 64 him, qiwm = intmom(a, eta2x, qr0r, r0r, 0.0, 0.0) 65 hi[i] = him * hs2[i] 66 qiw[i] = qiwm * hs2[i] 67 r0p = math.log(r0[i]) # noqa: F841 (unused in the Fortran too) 68 rphip = c.cln * hr[i] # noqa: F841 (unused in the Fortran too)