tridia.f

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

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

Variables on this page (7)

a, b, c, d, gama, json, u

Purpose

TRIDIA solves the tridiagonal system of the implicit angular-momentum diffusion step (the Thomas algorithm), and returns the change of angular momentum at each point.

How it fits

Inputs: TRIDI_MOD, filled by DCOEFT just before the call: the three diagonals A, B, C, the right-hand side D (by DCOEFT's comment, the angular velocities at the start of the time step, the outermost one already reduced by the wind), and GAMA(N), the wind term at the surface. GAMA is also the elimination's work array (loop 1 sets GAMA(2..N), line 29 ↗), so TRIDIA reads the wind term first (line 22 ↗). The argument EI holds the points' moments of inertia (by SECULR's use, probably).
Who calls it: only SECULR.
Outputs: U in TRIDI_MOD, the angular velocities at the end of the step; DJ, the change of angular momentum at each point; and SUMDJ, their sum. SECULR passes DJ to EQ2MOD to return to the model's grid.
Modules: TRIDI_MOD (formerly COMMON/TRIDI/), shared with DCOEFT, CCOEFT, CTRIDI, GETQUA and MIXCOM.

The method

Forward elimination and back substitution for $a_j u_{j-1} + b_j u_j + c_j u_{j+1} = d_j$; then $$\Delta J_j = (\Omega_j^\mathrm{new} - \Omega_j^\mathrm{old})\,I_j,$$ with, at the surface, the wind term added: $\Delta J_N = \gamma_N + (\Omega_N^\mathrm{new} - \Omega_N^\mathrm{old})\,I_N$. A zero pivot stops the run ("#TRIDIA:SINGULAR MATRIX").

Definitions

subroutineTRIDIA

lines 3–45 · GitHub

Thomas algorithm for the implicit angular-momentum diffusion, returning each point's change of angular momentum and their sum.

Uses modules DIMENSIONS_MOD, TRIDI_MOD

Called by SECULR

Compiler warnings (1)

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

line 31: Equality comparison for REAL(8) [-Wcompare-reals]

Change log

The changes since rp_1

Changed since rp_1. Declarations only: COMMON/TRIDI/ to TRIDI_MOD, JSON from DIMENSIONS_MOD, IMPLICIT NONE.

1C $$$$$$ 2C SUBROUTINE TRIDIA(N,EI,EJ,DJ,SUMDJ) ! KC 2025-05-31 3 SUBROUTINE TRIDIA(N,EI,DJ,SUMDJ) 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/TRIDI/ into TRIDI_MOD. 10 USE TRIDI_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 BET,DJ,EI,FJ,SUMDJ,U0 16 INTEGER I,J,N 17 18C DIMENSION EI(JSON),EJ(JSON),DJ(JSON),U0(JSON) ! KC 2025-05-31 19 DIMENSION EI(JSON),DJ(JSON),U0(JSON) 20 SAVE 21 22 DJ(N) = GAMA(N) 23 DO I = 1,N 24 U0(I) = D(I) 25 END DO 26 BET = B(1) 27 U(1) = D(1)/BET 28 DO 1 J=2,N 29 GAMA(J) = C(J-1)/BET 30 BET = B(J) - A(J)*GAMA(J) 31 IF( BET.EQ.0 ) STOP '#TRIDIA:SINGULAR MATRIX' 32 U(J) = (D(J) - A(J)*U(J-1))/BET 33 1 1 CONTINUE 34 DJ(N) = DJ(N)+(U(N)-U0(N))*EI(N) 35 SUMDJ = DJ(N) 36 FJ = 1.0D0+(U(N)-U0(N))/U0(N) 37 U(N) = U0(N)*FJ 38 DO 2 J=N-1,1,-1 39 U(J) = U(J) - GAMA(J+1)*U(J+1) 40 DJ(J) = (U(J)-U0(J))*EI(J) 41 SUMDJ = SUMDJ + DJ(J) 42 2 2 CONTINUE 43 44 RETURN 45 END

1"""Python version of tridia.f (YREC 7be4133): the tridiagonal solver of the 2implicit angular-momentum diffusion step, with each point's change of 3angular momentum.""" 4 5import numpy as np 6 7from pyrec import dimensions_mod 8from pyrec import tridi_mod 9 10# The Fortran's local array of starting angular velocities (SAVE). 11u0 = np.zeros(dimensions_mod.json + 1) 12 13 14def tridia(n: int, ei: np.ndarray, dj: np.ndarray, sumdj: float) -> float: 15 """Solve the system a, b, c, d of tridi_mod (set by dcoeft) for the new 16 angular velocities tridi_mod.u, and return each point's change of 17 angular momentum dj (the surface's including the wind term, which 18 dcoeft leaves in tridi_mod.gama[n]) and their sum. 19 20 Translates: tridia.f#TRIDIA (7be4133) 21 22 Arguments: n, the number of points; ei, the moments of inertia; dj, set 23 in place; sumdj is passed in with any value. 24 25 Returns sumdj. 26 27 """ 28 29 a = tridi_mod.a 30 b = tridi_mod.b 31 c = tridi_mod.c 32 d = tridi_mod.d 33 u = tridi_mod.u 34 gama = tridi_mod.gama 35 36 # The wind term, read from gama[n] before the elimination reuses gama 37 # as work space; the starting angular velocities (Fortran lines 22–25). 38 dj[n] = gama[n] 39 for i in range(1, n + 1): 40 u0[i] = d[i] 41 42 # Forward elimination (Fortran lines 26–33). 43 bet = b[1] 44 u[1] = d[1] / bet 45 for j in range(2, n + 1): # loop 1 46 gama[j] = c[j - 1] / bet 47 bet = b[j] - a[j] * gama[j] 48 if bet == 0: 49 raise SystemExit('#TRIDIA:SINGULAR MATRIX') 50 u[j] = (d[j] - a[j] * u[j - 1]) / bet 51 52 # The surface's change of angular momentum, wind included (Fortran 53 # lines 34–35). 54 dj[n] = dj[n] + (u[n] - u0[n]) * ei[n] 55 sumdj = dj[n] 56 57 # Question (see the page): a no-op, u[n] unchanged up to rounding 58 # (Fortran lines 36–37). 59 fj = 1.0 + (u[n] - u0[n]) / u0[n] 60 u[n] = u0[n] * fj 61 62 # Back substitution, and each point's change of angular momentum 63 # (Fortran lines 38–42). 64 for j in range(n - 1, 0, -1): # loop 2 65 u[j] = u[j] - gama[j + 1] * u[j + 1] 66 dj[j] = (u[j] - u0[j]) * ei[j] 67 sumdj = sumdj + dj[j] 68 69 return sumdj