taucal.f

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

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

Variables on this page (6)

cgl, cln, cmsun, json, lrot, taucz

Purpose

TAUCAL computes the convective overturn time scale $\tau_\mathrm{cz}$ (TAUCZ) near the base of the surface convection zone, for a convection zone that lies within the envelope (JVS 02/12). Its comments say the code was taken from MIDMOD, and that the routine is called only from ENVINT.

How it fits

Inputs: the envelope's profile, passed by ENVINT as arguments: the shell masses HS2, unlogged masses HS1, convective flags LCZ, $\log r$ (HR), $\log P$ (HP), $\log\rho$ (HD), gravity HG, the convective velocity SVEL, and two temperature gradients DEL1, DEL2 (EDEL1 and EDEL2 in ENVINT, copied from EDELS(1,·) and EDELS(2,·), which STITCH's column list names the radiative and the actual gradient; their difference changes sign at the zone's edge).
Who calls it: only ENVINT.
Output: TAUCZ in OVRTRN_MOD, shared with GETTAU, ENVINT, MASSLOSS, GETW and others (the wind and mass-loss routines, by their names).
Modules: OVRTRN_MOD, NMLPAR_MOD (LROT) and CONSTANTS_MOD.

The method

Going down from the surface, find the first convective shell, then the base of that zone, IMAX; note whether a second, deeper convective zone exists (LTWOCZ, its top ICZTOP); with two zones, IMAX ends at the base of the deeper one. For a single zone ICZTOP is not set (bug 1). The routine also sets HSTOP, HSBOT and LCZSUR, which nothing reads since the test on LCZSUR was commented out (JVS 10/11/13; unused?, lines 106–128 ↗); the test (HSTOP-HSBOT)/CMSUN .GT. 0 would be true for any zone.
Locate the base $r_b$ where $\nabla_1 - \nabla_2$ changes sign, by linear interpolation between shells IMAX-1 and IMAX (factor $f$), and the pressure scale height there, $$H_p = \frac{P}{\rho\,g},\qquad g = \frac{G m}{r^2}\ \text{(computed here only without rotation)},$$ interpolated with the same $f$.
If the zone is thin, that is, with $r_\mathrm{top}$ the radius of shell ICZTOP, $$\left|1 - \frac{r_b + H_p - r_\mathrm{IMAX}}{r_\mathrm{top} - r_\mathrm{IMAX}}\right| < 0.5,$$ so that the point $r_b + H_p$ lies between the zone's middle and half a width above its top, or there are two zones, the mass-weighted mean velocity over shells IMAX to ICZTOP is used: $$\bar v = \frac{\sum_j \Delta m_j\,\tfrac12 (v_j + v_{j+1})}{\sum_j \Delta m_j},\qquad \tau_\mathrm{cz} = \frac{|r_\mathrm{top} - r_\mathrm{IMAX}|}{\bar v}.$$

Question 1 (lines 148–153 ↗). .AND. became .OR. The commented-out test required a thin zone and two zones; the active one takes either. Which is intended?

Question 2 (line 159 ↗; lines 139–140 ↗). Two widths. In the averaging branch, the width is $r_\mathrm{top} - r_\mathrm{IMAX}$, from the shells, rather than from the interpolated base $r_b$ (ENVR). Is the shell's radius meant, or should the width start at $r_b$?
Otherwise (the original method), the velocity at one pressure scale height above the base, $v(r_b + H_p)$, by linear interpolation in $\log r$, gives $$\tau_\mathrm{cz} = \frac{H_p}{v(r_b + H_p)} .$$
ENVINT calls TAUCAL only with the old method (LNEWTCZ false; the default is true) and only when the envelope's deepest point is radiative and a point above it convective (JENV > 1, envint.f, lines 937–971 ↗). So LCZ(1) is false: when a convective point lies below M, the loop leaves the zone and sets IMAX ≥ 2; when only the surface point M is convective (the loop starts at M-1), LCZEX stays false and IMAX = M. IMAX is thus always set in the call, LALLCZ is never set, and the branch TAUCZ = 0 is not reached (settled from the code, Nov 9, 2026: lines 87–101 ↗, line 122 ↗, lines 223–224 ↗).

Question 3 (lines 120–124 ↗; lines 129–144 ↗). Only the surface point convective. When only point M is convective, IMAX = M, and $\tau_\mathrm{cz}$ is computed from shells M-1 and M as if a zone ended there: The search for $r_b + H_p$ then finds no shell (the loop from IMAX+1 to M is empty), and CVEL = SVEL(M). Should TAUCZ be 0 in that case?

Design

The routine prints "TauCal" to the screen at every call. A large commented-out block handled fully convective cases.

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. (not reached) `ICZTOP` is not set for a single zone.

(lines 87–98 ↗; lines 152–153 ↗; lines 159–163 ↗) ICZTOP, the top of the deeper zone, is set only when the loop finds two zones (line 94 ↗). But the test that chooses the averaging method evaluates HR(ICZTOP) for every zone (lines 152–153 ↗), and the averaging itself, which the first condition can select for a single zone, sums over shells IMAX to ICZTOP and takes the width from HR(ICZTOP). For a single zone, ICZTOP keeps the value SAVEd from an earlier call with two zones, and is undefined if there was none (an index of 0 would read outside HR). What follows: The choice of method and, if averaging is chosen, TAUCZ depend on an earlier model's zone. Only with the old method (LNEWTCZ false; the default is true), when the envelope holds the zone's base.

A possible fix: take the top of a single zone as ICZTOP:

taucal.f

@@ lines 94–99 @@
                   ICZTOP = I
             ELSE IF( LCZ(I) .AND. .NOT. LCZEX) THEN !First convective point from surface
                   LCZEX = .TRUE.
+C 11/09/26 KC via Claude
+C code changed to set ICZTOP for a single zone too: the thin-zone test
+C and the average below use it.  A second zone overwrites it.
+                  ICZTOP = I
             ENDIF
    71    CONTINUE
          IF(IMAX .LT. 1)THEN

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: TAUCAL runs only with LNEWTCZ off (envint.f, line 937 ↗), and no case of YREC turns it off. The fix was not run. What would reach it: LNEWTCZ = .FALSE., with the base of a single surface convection zone inside the envelope. No comment says a single zone should use an earlier ICZTOP.

Definitions

subroutineTAUCAL

lines 3–231 · GitHub

The convective overturn time near the base of a surface convection zone within the envelope, $H_p/v$ one scale height above the base, or a zone's width over its mean velocity for thin or double zones.

Uses modules CONSTANTS_MOD, DIMENSIONS_MOD, NMLPAR_MOD, OVRTRN_MOD

Called by ENVINT

Change log

The changes since rp_1

Changed since rp_1. Declarations only: COMMON blocks to modules (/OVRTRN/ to OVRTRN_MOD, /ROT/ to NMLPAR_MOD, the constants to CONSTANTS_MOD; the unused /COMP/ and /DEUTER/ dropped), the array size JSON from DIMENSIONS_MOD, IMPLICIT NONE.

1C$$$$$$ 2C SUBROUTINE TAUCAL(HCOMP,HS2,HS1,LCZ,HR,HP,HD,HG,M,SVEL,DEL1,DEL2,MINT) ! KC 2025-05-31 3 SUBROUTINE TAUCAL(HS2,HS1,LCZ,HR,HP,HD,HG,M,SVEL,DEL1,DEL2) 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/OVRTRN/ into OVRTRN_MOD. 10 USE OVRTRN_MOD 11C 10/14/26 MHP via Claude 12C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD. 13 USE NMLPAR_MOD 14C 10/14/26 MHP via Claude 15C code changed to convert COMMON/CONSTANTS/ into CONSTANTS_MOD. 16 USE CONSTANTS_MOD 17C 10/20/26 MHP via Claude 18C code changed to replace implicit typing with IMPLICIT NONE and 19C explicit declarations. 20 IMPLICIT NONE 21 REAL*8 CVEAVG,CVEL,CZWIDE,DD1,DD2,DEL1,DEL2,DENOM,ENVR,ENVRL,FX, 22 * * HD,HG,HP,HR,HS1,HS2,HSBOT,HSTOP,PS1,PS2,PSCA,RTESTL,SVEL, 23 * * TOL,UMER 24 INTEGER I,ICZTOP,J,K,M 25 LOGICAL*4 LALLCZ,LCZ,LCZEX,LCZSUR,LFLAG,LTWOCZ 26 INTEGER*4 IMAX 27C 10/13/26 MHP via Claude: removed the unused COMMON/COMP/ declaration 28C (no member referenced here; residue of the obsolete-option removals). 29C 10/13/26 MHP via Claude: removed the unused COMMON/DEUTER/ declaration 30C (no member referenced here; residue of the obsolete-option removals). 31C 10/12/26 MHP via Claude 32C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 33C namelist variables that the physics routines read from 27 COMMON 34C blocks which also held derived or working members. Those blocks keep 35C only their non-namelist members; five of them -- NEWMX, VMULT, 36C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 37C Derivation-only namelist variables from the same blocks went to 38C COMMON/NMLBYPASS/ instead, not here. 39C 40C Declared identically in all 121 files that hold any of these members. 41C Generated from one member list, so every copy is byte-identical; the 42C linker's 'size of symbol changed' warning is the check that they are. 43C 9/25/26 MHP via Claude 44C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/, 45C and COMMON/CONST3/ into a single COMMON/CONSTANTS/ 46C G Somers 3/17, ADDING NEW TAUCZ COMMON BLOCK 47C DIMENSION HCOMP(15,JSON),HS2(JSON),LCZ(JSON), ! KC 2025-05-31 48 DIMENSION HS2(JSON),LCZ(JSON), 49 * * HR(JSON), HP(JSON), HD(JSON), HG(JSON), HS1(JSON), ! XSUM(15) 50 * * DEL1(JSON), DEL2(JSON), SVEL(JSON) 51 SAVE 52 53 54 55C JVS 02/12 CALCULATE THE LOCAL CONVECTIVE OVERTURN TIMESCALE AT THE BASE 56C OF THE CZ. IN OLDER VERSIONS THIS WAS ONLY DONE FOR ROTATING MODELS; 57C THIS MAKES IT SO TAUCZ IS CALCULATED FOR ALL MODELS. THIS PARTICULAR SUBROUTINE 58C DOES CALCULATES THE OVERTURN TIMESCALES IN MODELS WHERE THE CZ IS WITHIN 59C THE ENVELOPE 60 61C THIS CODE SNAGGED FROM MIDMOD 62 63C DETERMINE EXTENT OF SURFACE CONVECTION ZONE. 64 LALLCZ = .FALSE. 65C JVS Allows the last point to be stable. 66C JVS 05/14 removed requirement for last point to be convective, since this routine 67C is only called from ENVINT, the near surface points in the envelope are often 68C radiative, despite the fact that a convective envelope exists. Code now checks 69C to make sure there is a CZ somewhere in the envelope. 70C IF(LCZ(M) .OR. LCZ(M-1) .OR. LCZ(M-2))THEN 71 72 73 74C SURFACE C.Z. EXISTS. FIND LOWEST SHELL (IMAX), WHICH IS ALSO THE 75C UPPERMOST ZONE CONSIDERED FOR STABILITY AGAINST ROTATIONALLY INDUCED MIXING. 76C DO 71 I = M-1,1,-1 77C IF(.NOT.LCZ(I)) GOTO 81 78C 71 CONTINUE 79C LALLCZ = .TRUE. 80C I = 0 81C 81 IMAX = I + 1 82 83C DEAL THE POSSIBILITY OF TWO OR MORE SURFACE CONVECTION ZONES IN THE ENVELOPE 84 LFLAG = .FALSE. 85 LTWOCZ = .FALSE. 86 LCZEX = .FALSE. ! Flag for the existence of a CZ 87 DO 71 I = M-1,1,-1 88 IF(.NOT.LCZ(I) .AND. .NOT. LFLAG .AND. LCZEX) THEN ! EXITING CZ 89 IMAX = I+1 90 LFLAG = .TRUE. 91 ELSE IF (LFLAG .AND. LCZ(I)) THEN ! TWO CZS WITH RADIATIVE ZONE BETWEEN 92 LFLAG = .FALSE. 93 LTWOCZ = .TRUE. 94 ICZTOP = I 95 ELSE IF( LCZ(I) .AND. .NOT. LCZEX) THEN !First convective point from surface 96 LCZEX = .TRUE. 97 ENDIF 98 71 71 CONTINUE 99 IF(IMAX .LT. 1)THEN 100 LALLCZ = .TRUE. 101 ENDIF 102 103 IF(LCZEX) THEN ! REPLACES "IF(LCZ(M)" THEN JVS 05/14 104 105 106C HSTOP IS THE MASS AT THE TOP OF THE C.Z. 107C HSBOT IS THE MASS AT THE BOTTOM OF THE C.Z. 108 HSTOP = HS1(M) 109 IF(IMAX.GT.1)THEN 110 HSBOT = 0.5D0*(HS1(IMAX)+HS1(IMAX-1)) 111 ELSE 112 HSBOT = HS1(1) 113 ENDIF 114C LCZSUR=T IF A SURFACE C.Z.DEEP ENOUGH FOR ANGULAR MOMENTUM LOSS EXISTS 115 IF((HSTOP-HSBOT)/CMSUN.GT.0.0D0)THEN 116 LCZSUR = .TRUE. 117 ELSE 118 LCZSUR= .FALSE. 119 ENDIF 120 ELSE 121C NO SURFACE C.Z. 122 IMAX = M 123 LCZSUR = .FALSE. 124 ENDIF 125 126C JVS 10/11/13 Concerned that this is inappropriate for LCZSUR = FALSE. 127C Should never be the case, but removed it. 128C IF(LCZSUR)THEN 129 IF(.NOT.LALLCZ)THEN 130 IF(.NOT.LROT)THEN 131 HG(IMAX)=HS1(IMAX)*EXP(CLN*(CGL-2.0D0*HR(IMAX))) 132 HG(IMAX-1)=HS1(IMAX-1)*EXP(CLN*(CGL-2.0D0*HR(IMAX-1))) 133 ENDIF 134C PINPOINT RCZ 135 DD2 = DEL1(IMAX-1)-DEL2(IMAX-1) 136 DD1 = DEL1(IMAX)-DEL2(IMAX) 137 FX = DD2/(DD2-DD1) 138C INFER HP 139 ENVRL = HR(IMAX-1)+FX*(HR(IMAX)-HR(IMAX-1)) 140 ENVR = EXP(CLN*ENVRL) 141 PS2 = EXP(CLN*(HP(IMAX)-HD(IMAX)))/HG(IMAX) 142 PS1 = EXP(CLN*(HP(IMAX-1)-HD(IMAX-1)))/HG(IMAX-1) 143 PSCA = PS1 + FX*(PS2-PS1) 144 RTESTL = DLOG10(ENVR+PSCA) 145C JVS 03/14 WHEN THERE ARE TWO CZS, THE DEEPER ONE CAN BE ~HP THICK. IN 146C THIS CASE, WE"LL WANT AND AVERAGE CONVECTIVE VELOCITY ACROSS THE REGION 147C-------------------------------------------------------------- 148 TOL = 0.5 149C IF(ABS(1.0 - ((DEXP(CLN*(RTESTL-HR(IMAX))) - 1.0)/ 150C * (DEXP(CLN*(HR(ICZTOP)-HR(IMAX))) - 1.0))) .LT. TOL 151C * .AND. LTWOCZ) THEN 152 IF((ABS(1.0 - ((DEXP(CLN*(RTESTL-HR(IMAX))) - 1.0)/ 153 * * (DEXP(CLN*(HR(ICZTOP)-HR(IMAX))) - 1.0))) .LT. TOL) .OR. LTWOCZ) THEN 154 ! TAKE AVERAGE CONV VELOCITY 155 CVEAVG = 0.0 156C DENOM = ABS(HS1(ICZTOP)-HS1(IMAX))/DEXP(CLN*STOTAL) ! MASS IN CZ 157 DENOM = 0.0 158 UMER = 0.0 159 CZWIDE = ABS(DEXP(CLN*HR(ICZTOP))-DEXP(CLN*HR(IMAX))) 160 DO J=IMAX,ICZTOP,1 161 DENOM = DENOM + HS2(J) 162 UMER = UMER+ 0.5d0*(SVEL(J)+SVEL(J+1))*HS2(J) 163 ENDDO 164 CVEAVG = UMER/DENOM 165 TAUCZ = CZWIDE/CVEAVG 166 ELSE 167C (ORIGINAL ROUTINE) 168C FIND V 169 DO K = IMAX+1,M 170 IF(HR(K).GT.RTESTL)THEN 171 FX = (RTESTL-HR(K-1))/(HR(K)-HR(K-1)) 172 CVEL = SVEL(K-1)+FX*(SVEL(K)-SVEL(K-1)) 173 GOTO 85 174 ENDIF 175 END DO 176 CVEL = SVEL(M) 177 85 85 CONTINUE 178C DEFINE TAUCZ 179 TAUCZ = PSCA/CVEL 180 ENDIF 181 182C JVS 10/11/13 183C ELSE 184CC INFER HP 185C IF(.NOT.LROT)THEN 186C HG(1)=HS1(1)*EXP(CLN*(CGL-2.0D0*HR(1))) 187C ENDIF 188C PSCA2 = EXP(CLN*(HP(1)-HD(1)))/HG(1) 189C RTEST2 = EXP(CLN*HR(1)) 190C IF(PSCA2.LE.RTEST2)THEN 191CC HP < R AT THE FIRST POINT. ASSUME V CONSTANT INSIDE AND HP = K/R FOR 192CC SLOWLY VARYING DENSITY AND PRESSURE NEAR THE CENTER. 193C CVEL = SVEL(1) 194C PSCA = (PSCA2*RTEST2)**0.5D0 195C TAUCZ = PSCA/CVEL 196C ELSE 197C DO K = 2,M 198C PSCA1 = PSCA2 199C RTEST1 = RTEST2 200C IF(.NOT.LROT)THEN 201C HG(K)=HS1(K)*EXP(CLN*(CGL-2.0D0*HR(K))) 202C ENDIF 203C PSCA2 = EXP(CLN*(HP(K)-HD(K)))/HG(K) 204C RTEST2 = EXP(CLN*HR(K)) 205CC FIND LOCATION WHERE HP = R 206C IF(PSCA2.LE.RTEST2)THEN 207C FX = (RTEST1-PSCA1)/((PSCA2-RTEST2)-(PSCA1-RTEST1)) 208CC FIND V 209C CVEL = SVEL(K-1)+FX*(SVEL(K)-SVEL(K-1)) 210C PSCA = PSCA1+FX*(PSCA2-PSCA1) 211CC DEFINE TAUCZ 212C TAUCZ = PSCA/CVEL 213C GOTO 95 214C ENDIF 215C END DO 216C K = M 217C CVEL = SVEL(M) 218C PSCA = PSCA2 219C TAUCZ = PSCA/CVEL 220C 95 CONTINUE 221C ENDIF 222C ENDIF 223 ELSE 224 TAUCZ = 0.0D0 225 ENDIF 226 print*, 'TauCal' 227 228C END JVS 229 230 RETURN 231 END

1"""Python version of taucal.f (YREC 7be4133): the convective overturn time 2near the base of a surface convection zone lying within the envelope.""" 3 4import math 5import numpy as np 6 7from pyrec import constants_mod 8from pyrec import nmlpar_mod 9from pyrec import ovrtrn_mod 10 11# Locals kept between calls (SAVE). Suspected bug (see the page): ICZTOP 12# is used for a single zone without being set in this call; 0 stands for 13# the compiler's usual zero initial value. 14imax = 0 15icztop = 0 16 17 18def taucal(hs2: np.ndarray, hs1: np.ndarray, lcz: np.ndarray, hr: np.ndarray, hp: np.ndarray, hd: np.ndarray, 19 hg: np.ndarray, m: int, svel: np.ndarray, del1: np.ndarray, del2: np.ndarray) -> None: 20 """Set ovrtrn_mod.taucz, the convective overturn time near the base of 21 the surface convection zone of the envelope: H_P / v one pressure scale 22 height above the base or, for a thin zone or two zones, the zone's 23 width over its mass-weighted mean velocity. 24 25 Translates: taucal.f#TAUCAL (7be4133) 26 27 Arguments: hs2 (shell masses), hs1 (masses), lcz (convective flags), 28 hr, hp, hd (log r, log P, log rho), hg (gravity; changed in place at 29 the base when there is no rotation), m, svel (convective velocities), 30 del1, del2 (two temperature gradients). None is returned. 31 32 """ 33 34 global imax, icztop 35 c = constants_mod 36 lallcz = False 37 lflag = False 38 ltwocz = False 39 lczex = False # a convection zone exists 40 41 # Down from the surface: the first convective shell, the zone's base 42 # IMAX, and a second, deeper zone (top ICZTOP) (Fortran lines 64–101). 43 for i in range(m - 1, 0, -1): 44 if not lcz[i] and not lflag and lczex: # leaving the zone 45 imax = i + 1 46 lflag = True 47 elif lflag and lcz[i]: # two zones with a radiative one between 48 lflag = False 49 ltwocz = True 50 icztop = i 51 elif lcz[i] and not lczex: # first convective shell from the top 52 lczex = True 53 if imax < 1: # never true where ENVINT calls TAUCAL (see the page) 54 lallcz = True 55 56 # The zone's mass range and LCZSUR, which nothing reads (Fortran lines 57 # 103–124). 58 if lczex: 59 hstop = hs1[m] 60 if imax > 1: 61 hsbot = 0.5 * (hs1[imax] + hs1[imax - 1]) 62 else: 63 hsbot = hs1[1] 64 if (hstop - hsbot) / c.cmsun > 0.0: 65 lczsur = True # noqa: F841 (unused in the Fortran too) 66 else: 67 lczsur = False # noqa: F841 (unused in the Fortran too) 68 else: 69 # Question (see the page): only the surface point convective, but 70 # a turnover time below. 71 imax = m 72 lczsur = False # noqa: F841 (unused in the Fortran too) 73 74 if not lallcz: 75 # Gravity at the base, without rotation (Fortran lines 130–133). 76 if not nmlpar_mod.lrot: 77 hg[imax] = hs1[imax] * math.exp(c.cln * (c.cgl - 2.0 * hr[imax])) 78 hg[imax - 1] = hs1[imax - 1] * math.exp( 79 c.cln * (c.cgl - 2.0 * hr[imax - 1])) 80 81 # The base r_b, where DEL1 - DEL2 changes sign, and H_P there 82 # (Fortran lines 135–144). 83 dd2 = del1[imax - 1] - del2[imax - 1] 84 dd1 = del1[imax] - del2[imax] 85 fx = dd2 / (dd2 - dd1) 86 envrl = hr[imax - 1] + fx * (hr[imax] - hr[imax - 1]) 87 envr = math.exp(c.cln * envrl) 88 ps2 = math.exp(c.cln * (hp[imax] - hd[imax])) / hg[imax] 89 ps1 = math.exp(c.cln * (hp[imax - 1] - hd[imax - 1])) / hg[imax - 1] 90 psca = ps1 + fx * (ps2 - ps1) 91 rtestl = math.log10(envr + psca) 92 93 # A thin zone or two zones: the mean velocity (Fortran lines 94 # 148–165). 95 tol = 0.5 96 # Question (see the page): .OR., not .AND. 97 # Suspected bug (see the page): ICZTOP is unset for a single zone. 98 if (abs(1.0 - ((math.exp(c.cln * (rtestl - hr[imax])) - 1.0) 99 / (math.exp(c.cln * (hr[icztop] - hr[imax])) - 1.0))) 100 < tol or ltwocz): 101 # The mass-weighted mean convective velocity. 102 # Question (see the page): the width from the shells. 103 denom = 0.0 104 umer = 0.0 105 czwide = abs(math.exp(c.cln * hr[icztop]) 106 - math.exp(c.cln * hr[imax])) 107 for j in range(imax, icztop + 1): 108 denom = denom + hs2[j] 109 umer = umer + 0.5 * (svel[j] + svel[j + 1]) * hs2[j] 110 cveavg = umer / denom 111 ovrtrn_mod.taucz = czwide / cveavg 112 else: 113 # The velocity one scale height above the base (Fortran lines 114 # 169–179). 115 for k in range(imax + 1, m + 1): 116 if hr[k] > rtestl: 117 fx = (rtestl - hr[k - 1]) / (hr[k] - hr[k - 1]) 118 cvel = svel[k - 1] + fx * (svel[k] - svel[k - 1]) 119 break # to label 85 120 else: 121 cvel = svel[m] 122 ovrtrn_mod.taucz = psca / cvel # label 85 123 else: 124 ovrtrn_mod.taucz = 0.0 # Fortran line 224. 125 126 # Question (see the page): a debugging print at every call (Fortran 127 # line 226). 128 print(" TauCal")

taucal.f

Purpose

How it fits

The method

Design

Suspected bugs, checked by running the code

1. (not reached) ICZTOP is not set for a single zone.

Definitions

subroutineTAUCAL

Change log

1. (not reached) `ICZTOP` is not set for a single zone.