1C $$$$$$
2 SUBROUTINE NDIFCOM(DT,COD2,ECOD2,EM,HD,HL,HP,HR,HS,HS1,HS2,
3C * HSTOT,HV,IBEG,IEND,IMAX,IMIN,LCZ,LOK,M,NTOT, ! KC 2025-05-31
4 * * HSTOT,HV,IBEG,IEND,IMAX,IMIN,LCZ,LOK,M,
5 * * HCOMP,JBEG,JEND)
6C 10/15/26 MHP via Claude
7C code changed to take JSON from DIMENSIONS_MOD instead of declaring
8C PARAMETER(JSON=5000) here.
9 USE DIMENSIONS_MOD
10C 10/20/26 MHP via Claude
11C code changed to replace implicit typing with IMPLICIT NONE and
12C explicit declarations.
13 IMPLICIT NONE
14 REAL*8 COD2,DT,ECOD2,EM,HCOMP,HD,HL,HP,HR,HS,HS1,HS2,HSTOT,HV
15 INTEGER IBEG,IEND,IMAX,IMIN,ISTART,J,JBEG,JEND,M
16 LOGICAL*4 LCZ,LDUM,LDUM2,LOK,LTEST
17 DIMENSION COD2(JSON),ECOD2(JSON),EM(JSON),HD(JSON),HL(JSON),
18 * * HP(JSON),HR(JSON),
19 * * HCOMP(15,JSON),HS(JSON),HS1(JSON),HS2(JSON),HV(JSON),
20 * * LCZ(JSON)
21 SAVE
22C DIFCOM CALCULATES THE DIFFUSION OF COMPOSITION DUE TO ANGULAR MOMENTUM
23C TRANSPORT. THIS IS DONE BY TRANSFORMING TO AN EQUALLY SPACED GRID IN
24C RADIUS, SOLVING A DIFFUSION EQUATION, AND TRANSFORMING BACK.
25C
26C INPUT VARIABLES:
27C DT : DIFFUSION TIMESTEP(SEC).
28C COD2 : DIFFUSION COEFFICIENTS FOR COMPOSITION TRANSPORT AT THE ORIGINAL
29C MODEL POINTS.
30C ECOD2 : DIFFUSION COEFFICIENTS FOR COMPOSITION TRANSPORT AT THE EQUALLY
31C SPACED GRID POINTS.
32C EM : MASSES OF THE EQUALLY SPACED GRID POINTS(GM).
33C NOTE: FOR CONVECTIVE BOUNDARIES THE MASS OF THE LAST GRID POINT IS
34C THE MASS OF THE ENTIRE CONVECTION ZONE.
35C HCOMP : ARRAY OF MASS FRACTION OF ALL OF THE SPECIES AT THE ORIGINAL
36C MODEL POINTS.
37C HS1 : LOCATION IN MASS(UNLOGGED) OF THE ORIGINAL MODEL POINTS.
38C HS2 : MASSES OF THE ORIGINAL MODEL POINTS(UNLOGGED).
39C HV : RUN OF DIFFUSION VELOCITIES.
40C IBEG,IEND :THE FIRST/LAST UNSTABLE POINTS IN THE REGION.
41C NOTE: FOR CONVECTIVE BOUNDARIES THESE ARE ONLY THE FIRST CONVECTIVE
42C POINTS ADJACENT TO AN UNSTABLE RADIATIVE REGION.
43C IMIN : THE INNERMOST RADIATIVE ZONE OUTSIDE OF ZONE 1.
44C LCZ : FLAG WHICH TELLS WHICH OF THE ORIGINAL MODEL POINTS ARE
45C CONVECTIVE FOR ANGULAR MOMENTUM TRANSPORT PURPOSES(I.E. INCLUDES
46C OVERSHOOT REGIONS). LCZ=T IF CONVECTIVE.
47C LOK : FLAG SET T WHEN FINAL ITERATION IS BEING PERFORMED.
48C M : NUMBER OF MODEL POINTS.
49
50C OUTPUT VARIABLES:
51C HCOMP IS UPDATED IN DIFCOM TO GIVE THE NEW RUN OF COMPOSITION AFTER
52C ANGULAR MOMENTUM TRANSPORT.
53C
54C BEFORE THE LAST ITERATION(LOK=F),ONLY DIFFUSION OF H,HE4,HE3 CALCULATED
55C TO CALCULATE CHANGE IN MU GRADIENTS CAUSED BY DIFFUSION.
56C
57C WHEN LOK IS T, ONLY DIFFUSION OF SPECIES HEAVIER THAN HE4 IS PERFORMED
58C (AS DIFFUSION OF LIGHTER SPECIES HAS ALREADY BEEN DONE).
59C THE EQUALLY SPACED GRID POINTS FROM THE LAST UNSTABLE REGION SOLVED
60C ARE USED INITIALLY, AND THEN THE PROGRAM CHECKS FOR ADDITIONAL UNSTABLE
61C REGIONS INTERIOR TO THAT.
62C CALL MIXCOM FOR PREVIOUSLY IDENTIFIED UNSTABLE REGION IF THIS IS NOT
63C THE CONVERGED CALL
64 IF(.NOT.LOK) THEN
65 CALL MIXCOM(DT,ECOD2,EM,HS2,IBEG,IEND,LCZ,LOK,M,HCOMP,
66 * * JBEG,JEND)
67 ELSE
68C FIND UNSTABLE REGIONS IN ORDER, AND CALL MIXGRID TO SET UP THE
69C EQUALLY SPACED GRID AND MIXCOM TO MIX THEM IN ORDER
70C EACH UNSTABLE REGION IS SOLVED SEPARATELY STARTING HERE.
71C LTEST IS SET T IF A NON-ZERO VELOCITY IS ENCOUNTERED.
72C IBEG IS THE ZONE BELOW THE FIRST NON-ZERO V;IEND IS THE ZONE ABOVE
73C THE LAST NON-ZERO V.
74 LDUM = .FALSE.
75 ISTART = IMIN
76 60 60 CONTINUE
77 LTEST = .FALSE.
78 DO J = ISTART,IMAX
79 IF(HV(J).GT.0.0D0) THEN
80 LTEST = .TRUE.
81 IF(.NOT.LDUM) THEN
82C START OF UNSTABLE REGION
83 LDUM = .TRUE.
84 IBEG = J - 1
85 ENDIF
86 ELSE IF(LDUM) THEN
87C END OF UNSTABLE REGION
88 IEND = J - 1
89 LDUM = .FALSE.
90 ISTART = J + 1
91 GOTO 80
92 ENDIF
93 END DO
94C IF THE LAST INTERFACE IS UNSTABLE (NON-ZERO V) ENSURE THAT IEND IS SET
95C PROPERLY.
96 IF(LDUM) IEND = IMAX
97 ISTART = IMAX + 1
98 80 80 CONTINUE
99C IF NO NON-ZERO V'S ENCOUNTERED, EXIT.
100 IF(.NOT.LTEST) GOTO 90
101C TRANSFORM TO EQUAL GRID SPACING IN R FOR THE REGION.
102 CALL MIXGRID(COD2,HD,HL,HP,HR,HS,HS1,HS2,
103 * * HSTOT,IBEG,IEND,LCZ,M,ECOD2,EM,LDUM2)
104C LDUM2=T IF TWO ZONES IN UNSTABLE REGION;
105C SKIP IF THIS OCCURS.
106 IF(LDUM2) THEN
107 IF(ISTART.LE.IMAX) THEN
108 GOTO 60
109 ELSE
110 GOTO 90
111 ENDIF
112 ENDIF
113C PERFORM COMPOSITION DIFFUSION.
114C UNTIL THE FINAL ITERATION, ONLY COMPOSITION DIFFUSION OF SPECIES WHICH
115C AFFECT GRADIENTS IN MEAN MOLECULAR WEIGHT IS COMPUTED (H,HE3,HE4).
116C ON THE FINAL ITERATION, DIFFUSION OF ALL SPECIES IS PERFORMED.
117 CALL MIXCOM(DT,ECOD2,EM,HS2,IBEG,IEND,LCZ,LOK,M,HCOMP,
118 * * JBEG,JEND)
119C RETURN FOR NEXT REGION IF APPLICABLE
120 IF(ISTART.LE.IMAX) GOTO 60
121 90 90 CONTINUE
122 ENDIF
123 RETURN
124 END
1"""Python version of ndifcom.f (YREC 7be4133): composition mixed across
2each region made unstable by rotation, through MIXGRID and MIXCOM."""
3
4from pyrec.mixcom import mixcom
5import numpy as np
6from pyrec.mixgrid import mixgrid
7
8# The Fortran's locals (SAVE).
9ldum = False
10ldum2 = False
11
12
13def ndifcom(dt: float, cod2: np.ndarray, ecod2: np.ndarray, em: np.ndarray, hd: np.ndarray, hl: np.ndarray,
14 hp: np.ndarray, hr: np.ndarray, hs: np.ndarray, hs1: np.ndarray, hs2: np.ndarray, hstot: float,
15 hv: np.ndarray, ibeg: int, iend: int, imax: int, imin: int, lcz: np.ndarray, lok: bool, m: int,
16 hcomp: np.ndarray, jbeg: int, jend: int) -> tuple[int, int]:
17 """Before the final iteration (lok False), mix rows jbeg..jend over the
18 region found earlier (ibeg..iend, with its grid ecod2, em). On the
19 final iteration, find each run of interfaces with hv > 0 between imin
20 and imax, set up its grid (MIXGRID) and mix it (MIXCOM); regions of a
21 single unstable interface are skipped.
22
23 Translates: ndifcom.f#NDIFCOM (7be4133)
24
25 Arguments: the step dt (s); the coefficients cod2 at the model points;
26 ecod2 and em on the equal grid (kept between calls by the caller); the
27 model hd, hl, hp, hr, hs, hs1, hs2, hstot; the velocities hv; ibeg,
28 iend (the last region); the search range imin..imax; the flags lcz; lok;
29 m shells; hcomp[k, j]; the rows jbeg..jend. Arrays have index 0 unused.
30
31 Sets in place: hcomp, and ecod2, em for the last region.
32
33 Returns (ibeg, iend), the last region.
34
35 """
36
37 global ldum, ldum2
38
39 # Before the final call: the region the caller found (Fortran lines 64–66).
40 if not lok:
41 mixcom(dt, ecod2, em, hs2, ibeg, iend, lcz, lok, m, hcomp, jbeg, jend)
42 return ibeg, iend
43
44 # Find the unstable regions in order: ibeg is the zone below the first
45 # non-zero v, iend the zone above the last (Fortran lines 74–100).
46 ldum = False
47 istart = imin
48 while True: # label 60
49 ltest = False
50 ended = False
51 for j in range(istart, imax + 1):
52 if hv[j] > 0.0:
53 ltest = True
54 if not ldum:
55 # Start of an unstable region.
56 ldum = True
57 ibeg = j - 1
58 elif ldum:
59 # End of an unstable region.
60 iend = j - 1
61 ldum = False
62 istart = j + 1
63 ended = True
64 break # GOTO 80
65 if not ended:
66 # The last interface is unstable: iend is IMAX.
67 if ldum:
68 iend = imax
69 istart = imax + 1
70 # label 80: no non-zero v, exit.
71 if not ltest:
72 break # GOTO 90
73
74 # The region's grid; a single unstable interface is skipped (Fortran lines 102–112).
75 # Question (see the page): a single unstable interface is skipped.
76 ldum2 = mixgrid(cod2, hd, hl, hp, hr, hs, hs1, hs2, hstot, ibeg,
77 iend, lcz, m, ecod2, em, ldum2)
78 if ldum2:
79 if istart <= imax:
80 continue # GOTO 60
81 break # GOTO 90
82
83 # Mix the region, then back for the next one, if any (Fortran lines 117–120).
84 mixcom(dt, ecod2, em, hs2, ibeg, iend, lcz, lok, m, hcomp, jbeg, jend)
85 if istart > imax:
86 break
87
88 # label 90
89 return ibeg, iend
PyREC scaffolding (from YREC development, 7be4133)