dcoeft.f

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

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

Variables on this page (9)

a, b, c, d, ecod3, ecod4, gama, json, ldifad

Purpose

This file sets up the tridiagonal system for one implicit step of angular-momentum transport by diffusion (and, with LDIFAD, a simple advective term), on the even grid in the coordinate $\chi$ (DR is DCHI, set by ROTGRID). It is the solver setup used when the band solver of DADCOEFT is not.

How it fits

Inputs: the arguments ECOD (the diffusion coefficients at the grid points), DR, DT, EI (moments of inertia), EJ (angular momenta), EW ($\Omega$), NTOT, WIND1, WIND2 (the wind's change of $\Omega$, explicit and implicit), LFIX (fixed surface $\Omega$); with LDIFAD (a setting in NMLPAR_MOD, formerly /VARFC/), the advective and diffusive coefficients ECOD3, ECOD4 from DIFAD_MOD. ROTGRID has already multiplied ECOD by $\rho r^2\,(d\chi/dr)$ and by the shell's $I/\Delta m$ ($4\pi$ is in CODIFF), and with LDIFAD added the composition coefficient ECOD2 to ECOD4.
Who calls it: SECULR.
Outputs: A, B, C, D in TRIDI_MOD, solved next by TRIDIA; also GAMA(NTOT), set to the wind's loss of angular momentum, $\Delta J = -\tfrac12(W_1 + W_2)\,I_N$, which TRIDIA reads before overwriting it (line 82 ↗; tridia.f, line 22 ↗). GAMA is the solver's work array, so this is a message passed through work space, which an explicit argument would make clearer.

The method

The header gives $$\frac{\partial\Omega}{\partial t} = \frac{1}{4\pi\rho r^2}\,\frac{1}{I/M}\,\frac{\partial}{\partial r}\left(D\frac{\partial\Omega}{\partial r}\right),$$ differenced implicitly with $f_i = \Delta t/(\Delta r\,I_i)$: $$a_i = -f_i D_i, \qquad b_i = 1 + f_i(D_i + D_{i+1}), \qquad c_i = -f_i D_{i+1}, \qquad d_i = \Omega_i^n,$$ with no flux at the centre and at the surface (or a fixed surface $\Omega$ if LFIX). The wind enters the surface: $d_N = \Omega_N(1 + \Delta J/J_N)$. With LDIFAD, $D$ is ECOD4 and centred advective terms $\pm\tfrac12\Delta t\,\mathtt{ECOD3}/I_i$ are added, in conservative form: each interface's flux, $\tfrac12\mathtt{ECOD3}_i(\Omega_{i-1}+\Omega_i)$, enters both neighbouring rows, and none crosses the centre or the surface. This LDIFAD branch (lines 108–139 ↗) is unused?: SECULR calls DCOEFT only when LDIFAD is off and calls DADCOEFT otherwise (seculr.f, lines 390–404 ↗; settled from the code, Nov 9, 2026). Its advective factor, FACTA = 0.5D0*DT/EI(I) (line 116 ↗, line 122 ↗, line 130 ↗), has no $1/\Delta r$, unlike DADCOEFT's 0.5D0*DT/EI(1)/DR (dadcoeft.f, line 305 ↗).

Definitions

subroutineDCOEFT

lines 2–142 · GitHub

Builds the implicit tridiagonal system for $\Omega$ on the even grid (diffusion, optionally a centred advective term), with zero-flux or fixed-$\Omega$ boundaries and the wind's loss at the surface.

Uses modules DIFAD_MOD, DIMENSIONS_MOD, NMLPAR_MOD, TRIDI_MOD

Called by SECULR

Change log

The changes since rp_1

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

1C $$$$$$ 2 SUBROUTINE DCOEFT(ECOD,DR,DT,EI,EJ,EW,NTOT,WIND1,WIND2,LFIX) 3C 10/15/26 MHP via Claude 4C code changed to take JSON from DIMENSIONS_MOD instead of declaring 5C PARAMETER(JSON=5000) here. 6 USE DIMENSIONS_MOD 7C 10/14/26 MHP via Claude 8C code changed to convert COMMON/TRIDI/ into TRIDI_MOD. 9 USE TRIDI_MOD 10C 10/14/26 MHP via Claude 11C code changed to convert COMMON/DIFAD/ into DIFAD_MOD. 12 USE DIFAD_MOD 13C 10/14/26 MHP via Claude 14C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD. 15 USE NMLPAR_MOD 16C 10/20/26 MHP via Claude 17C code changed to replace implicit typing with IMPLICIT NONE and 18C explicit declarations. 19 IMPLICIT NONE 20 REAL*8 DR,DT,ECOD,EI,EJ,EW,FACT,FACT0,FACTA,WIND1,WIND2 21 INTEGER I,NTOT 22 LOGICAL*4 LFIX 23 24C MHP 7/93 25C 10/12/26 MHP via Claude 26C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 27C namelist variables that the physics routines read from 27 COMMON 28C blocks which also held derived or working members. Those blocks keep 29C only their non-namelist members; five of them -- NEWMX, VMULT, 30C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 31C Derivation-only namelist variables from the same blocks went to 32C COMMON/NMLBYPASS/ instead, not here. 33C 34C Declared identically in all 121 files that hold any of these members. 35C Generated from one member list, so every copy is byte-identical; the 36C linker's 'size of symbol changed' warning is the check that they are. 37C 10/13/26 MHP via Claude 38C code changed to drop COMMON/VARFC/, declared here but never used. 39 DIMENSION ECOD(JSON),EI(JSON),EW(JSON),EJ(JSON) 40 SAVE 41C DCOEFT SETS UP THE COEFFICIENTS FOR THE DIFFUSION DIFFERENCE EQUATION. 42C 43C INPUT VARIABLES: 44C ECOD : DIFFUSION COEFFICIENTS AT THE EQUALLY SPACED GRID POINTS. 45C DR : GRID SPACING. 46C DT : TIMESTEP. 47C EI : RUN OF MOMENTS OF INERTIA OF EQUALLY SPACED GRID POINTS. 48C EW : RUN OF ANGULAR VELOCITY OF EQUALLY SPACED GRID POINTS. 49C NTOT : NUMBER OF EQUALLY SPACED GRID POINTS. 50C WIND1 : THE ANGULAR MOMENTUM LOSS FROM A SURFACE C.Z. COMPUTED 51C EXPLICITILY. 52C WIND2 : AS WIND1, BUT COMPUTED IMPLICITLY. 53C *NOTE: IF NOT APPLICABLE, WIND1 AND WIND2 ARE ZEROED OUT. 54C 55C OUTPUT VARIABLES : 56C THE ANGULAR VELOCITY OF SHELL I AT TIME N+1 (W(I,N+1)) IS A FUNCTION OF 57C W(I-1),W(I),AND W(I+1) AS DISCUSSED BELOW. THIS CAN BE EXPRESSED 58C AS A TRIDIAGONAL MATRIX. 59C A(I) : CONTAINS ALL TERMS INVOLVING OMEGA(I-1). 60C B(I) : CONTAINS ALL TERMS INVOLVING OMEGA(I). 61C C(I) : CONTAINS ALL TERMS INVOLVING OMEGA(I+1) 62C D(I) : THE ANGULAR VELOCITY TERMS AT THE BEGINNING OF THE TIMESTEP. 63C *NOTE: ANGULAR MOMENTUM LOSS FROM A SURFACE C.Z. IS SUBTRACTED FROM 64C THE LAST ELEMENT OF ARRAY D IF APPLICABLE. 65C 66C THE DIFFUSION EQUATION WE ARE SOLVING IS 67C dW/dT = 1/(4pi*RHO*R**2) 1/(I/M) d/dR(D*dW/dR) 68C WHERE D IS OUR DIFFUSION COEFFICIENT,W IS THE ANGULAR VELOCITY, 69C I IS THE MOMENT OF INERTIA,M IS MASS,T IS TIME,R IS RADIUS, 70C AND RHO IS DENSITY. 71C USING M = 4pi*RHO*R**2*DR THIS IS DIFFERENCED FOR SHELL I AS 72C W(I,N+1)-W(I,N)=(DT/DR)*(1/I)*(D(I+1/2)*(W(I+1,N+1)-W(I,N+1)) - 73C D(I-1/2)*(W(I,N+1)-W(I-1,N+1))) 74C WHERE INDEX N+1 REFERS TO VALUES AT THE END OF THE TIMESTEP,AND 75C INDEX N REFERS TO VALUES AT THE BEGINNING OF THE TIMESTEP. 76C THIS SYSTEM IS SUPPLEMENTED WITH APPROPRIATE BOUNDARY CONDITIONS. 77C 78 FACT0 = DT/DR 79 DO 10 I = 1,NTOT 80 D(I) = EW(I) 81 10 10 CONTINUE 82 GAMA(NTOT) = -0.5D0*(WIND1+WIND2)*EI(NTOT) 83 D(NTOT) = D(NTOT)*(1.0D0+GAMA(NTOT)/EJ(NTOT)) 84 IF(.NOT.LDIFAD)THEN 85C FIRST SHELL BOUNDARY CONDITIONS: NO FLOW BELOW THE BOUNDARY 86C I.E. THE ANGULAR MOMENTUM TRANSPORT AT THE FIRST SHELL DOES NOT 87C DEPEND ON THE SHELLS BELOW IT. 88 FACT = FACT0/EI(1) 89 A(1) = 0.0D0 90 B(1) = 1.0D0 + FACT*ECOD(2) 91 C(1) = -FACT*ECOD(2) 92 DO 20 I = 2,NTOT-1 93 FACT = FACT0/EI(I) 94 A(I) = -FACT*ECOD(I) 95 B(I) = 1.0D0 + FACT*(ECOD(I)+ECOD(I+1)) 96 C(I) = -FACT*ECOD(I+1) 97 20 20 CONTINUE 98C LAST SHELL B.C. : SAME AS FIRST SHELL B.C. 99 FACT = FACT0/EI(NTOT) 100 IF(.NOT.LFIX)THEN 101 A(NTOT) = -FACT*ECOD(NTOT) 102 B(NTOT) = 1.0D0 + FACT*ECOD(NTOT) 103 ELSE 104 A(NTOT) = 0.0D0 105 B(NTOT) = 1.0D0 106 ENDIF 107 C(NTOT) = 0.0D0 108 ELSE 109C IF LDIFAD=T, WE ARE SOLVING A COMBINED DIFFUSION/ADVECTION EQUATION. 110C THIS ADDS A TERM D/DR(RHO*R**4*V*W) TO THE ORIGINAL D/DR(RHO*R**4* 111C V*R*DW/DR) EQUATION. 112C FIRST SHELL BOUNDARY CONDITIONS: NO FLOW BELOW THE BOUNDARY 113C I.E. THE ANGULAR MOMENTUM TRANSPORT AT THE FIRST SHELL DOES NOT 114C DEPEND ON THE SHELLS BELOW IT. 115 FACT = FACT0/EI(1) 116 FACTA = 0.5D0*DT/EI(1) 117 A(1) = 0.0D0 118 B(1) = 1.0D0 + FACT*ECOD4(2) - FACTA*ECOD3(2) 119 C(1) = -FACT*ECOD4(2) - FACTA*ECOD3(2) 120 DO I = 2,NTOT-1 121 FACT = FACT0/EI(I) 122 FACTA = 0.5D0*DT/EI(I) 123 A(I) = -FACT*ECOD4(I) + FACTA*ECOD3(I) 124 B(I) = 1.0D0 + FACT*(ECOD4(I)+ECOD4(I+1)) + 125 * * FACTA*(ECOD3(I)-ECOD3(I+1)) 126 C(I) = -FACT*ECOD4(I+1) - FACTA*ECOD3(I+1) 127 END DO 128C LAST SHELL B.C. : SAME AS FIRST SHELL B.C. 129 FACT = FACT0/EI(NTOT) 130 FACTA = 0.5D0*DT/EI(NTOT) 131 IF(.NOT.LFIX)THEN 132 A(NTOT) = -FACT*ECOD4(NTOT) + FACTA*ECOD3(NTOT) 133 B(NTOT) = 1.0D0 + FACT*ECOD4(NTOT) + FACTA*ECOD3(NTOT) 134 ELSE 135 A(NTOT) = 0.0D0 136 B(NTOT) = 1.0D0 137 ENDIF 138 C(NTOT) = 0.0D0 139 ENDIF 140 141 RETURN 142 END

1"""Python version of dcoeft.f (YREC 7be4133): the tridiagonal system for 2one implicit step of angular-momentum transport on the even grid.""" 3 4import numpy as np 5 6from pyrec import difad_mod 7from pyrec import nmlpar_mod 8from pyrec import tridi_mod 9 10 11def dcoeft(ecod: np.ndarray, dr: float, dt: float, ei: np.ndarray, ej: np.ndarray, ew: np.ndarray, ntot: int, 12 wind1: float, wind2: float, lfix: bool) -> None: 13 """Set up the implicit tridiagonal system for Omega (diffusion, and with 14 ldifad a centred advective term), with no flux at the centre, no flux 15 or a fixed Omega at the surface, and the wind's loss of angular 16 momentum in the surface row. 17 18 Translates: dcoeft.f#DCOEFT (7be4133) 19 20 Arguments: ecod, the diffusion coefficients; dr, the grid step in chi; 21 dt, the step; ei, ej, ew, the moments of inertia, angular momenta and 22 Omega; ntot, the number of points; wind1, wind2, the wind's change of 23 Omega (explicit, implicit); lfix, a fixed surface Omega. 24 25 Returns nothing: sets tridi_mod.a, b, c, d and gama[ntot] (the wind's 26 Delta J, read by TRIDIA). 27 28 """ 29 30 a = tridi_mod.a 31 b = tridi_mod.b 32 c = tridi_mod.c 33 d = tridi_mod.d 34 35 # The starting Omega, and the wind's Delta J in the surface row 36 # (passed to TRIDIA in gama[ntot]) (Fortran lines 78–83). 37 fact0 = dt / dr 38 for i in range(1, ntot + 1): # label 10 39 d[i] = ew[i] 40 tridi_mod.gama[ntot] = -0.5 * (wind1 + wind2) * ei[ntot] 41 d[ntot] = d[ntot] * (1.0 + tridi_mod.gama[ntot] / ej[ntot]) 42 43 if not nmlpar_mod.ldifad: 44 # Diffusion only; no flux below the first shell (Fortran lines 45 # 84–91). 46 fact = fact0 / ei[1] 47 a[1] = 0.0 48 b[1] = 1.0 + fact * ecod[2] 49 c[1] = -fact * ecod[2] 50 51 # The interior rows (Fortran lines 92–97). 52 for i in range(2, ntot): # label 20 53 fact = fact0 / ei[i] 54 a[i] = -fact * ecod[i] 55 b[i] = 1.0 + fact * (ecod[i] + ecod[i + 1]) 56 c[i] = -fact * ecod[i + 1] 57 58 # The surface row: no flux, or a fixed Omega (Fortran lines 59 # 98–107). 60 fact = fact0 / ei[ntot] 61 if not lfix: 62 a[ntot] = -fact * ecod[ntot] 63 b[ntot] = 1.0 + fact * ecod[ntot] 64 else: 65 a[ntot] = 0.0 66 b[ntot] = 1.0 67 c[ntot] = 0.0 68 69 else: 70 # Diffusion (ecod4) and advection (ecod3); unused? (SECULR calls 71 # DCOEFT only without ldifad). The first row (Fortran lines 72 # 108–119). 73 ecod3 = difad_mod.ecod3 74 ecod4 = difad_mod.ecod4 75 fact = fact0 / ei[1] 76 facta = 0.5 * dt / ei[1] 77 a[1] = 0.0 78 b[1] = 1.0 + fact * ecod4[2] - facta * ecod3[2] 79 c[1] = -fact * ecod4[2] - facta * ecod3[2] 80 81 # The interior rows (Fortran lines 120–127). 82 for i in range(2, ntot): 83 fact = fact0 / ei[i] 84 facta = 0.5 * dt / ei[i] 85 a[i] = -fact * ecod4[i] + facta * ecod3[i] 86 b[i] = (1.0 + fact * (ecod4[i] + ecod4[i + 1]) 87 + facta * (ecod3[i] - ecod3[i + 1])) 88 c[i] = -fact * ecod4[i + 1] - facta * ecod3[i + 1] 89 90 # The surface row (Fortran lines 128–138). 91 fact = fact0 / ei[ntot] 92 facta = 0.5 * dt / ei[ntot] 93 if not lfix: 94 a[ntot] = -fact * ecod4[ntot] + facta * ecod3[ntot] 95 b[ntot] = 1.0 + fact * ecod4[ntot] + facta * ecod3[ntot] 96 else: 97 a[ntot] = 0.0 98 b[ntot] = 1.0 99 c[ntot] = 0.0