Turbulence

Modules implementing the two-equation turbulence closures and algebraic stability functions. Each module corresponds to one GOTM Fortran subroutine.

Turbulence closure modules.

Dynamic transport equation for turbulent kinetic energy \(k\).

Implements GOTM Section 4.7.23 (tkeeq.F90) — solves the \(k\)-equation for a Boussinesq fluid in horizontally homogeneous flow (Eq. 150):

\[\dot{k} = \mathcal{D}_k + P + G + P_x + P_s - \varepsilon \comma\]

where \(P\) is shear production, \(G\) is buoyancy production (negative in stable stratification), \(P_s\) is Stokes shear production, \(P_x\) is extra production, and \(\varepsilon\) is the dissipation rate. The diffusive transport is (Eq. 151):

\[\mathcal{D}_k = \frac{\partial}{\partial z}\left( \frac{\nu_t}{\sigma_k} \frac{\partial k}{\partial z}\right) \comma\]

with Schmidt number \(\sigma_k\). The dissipation rate \(\varepsilon\) is supplied either from the dynamic dissipationeq or omegaeq, or from a prescribed length-scale closure (Eq. 153):

\[\varepsilon = (c_\mu^0)^3 \frac{k^{3/2}}{l} \point\]

Boundary conditions

At solid boundaries (surface and bottom), boundary conditions follow logarithmic-layer theory (Dirichlet):

\[k = \frac{u_\tau^2}{(c_\mu^0)^{1/2}} \comma\]

or wave-breaking injection following Craig and Banner (1994):

\[F_k = \eta u_{\tau s}^3 \comma\]

where \(\eta \approx 100\) is the Craig–Banner coefficient. The injection boundary condition results in a non-trivial k-profile near the surface (Dirichlet or Neumann depending on k_ubc/k_lbc).

Authors (original Fortran): Lars Umlauf (rewrite), Hans Burchard, Karsten Bolding.

class pygotm.turbulence.tkeeq.TKEEquationWorkspace(nlev, *, n_cols=None)[source]

Bases: ColumnWorkspace

Workspace arrays for the translated TKE equation.

Parameters:

  • nlev (int)

  • n_cols (int | None)

tke: ndarray
tkeo: ndarray
h: ndarray
P: ndarray
B: ndarray
Px: ndarray
PSTK: ndarray
num: ndarray
eps: ndarray
avh: ndarray
l_sour: ndarray
q_sour: ndarray
u_taus: ndarray
u_taub: ndarray
z0s: ndarray
z0b: ndarray
au: ndarray
bu: ndarray
cu: ndarray
du: ndarray
ru: ndarray
qu: ndarray
pygotm.turbulence.tkeeq.step_tkeeq(batch_size, nlev, dt, sig_k, k_min, k_ubc, k_lbc, ubc_type, lbc_type, cm0, cmsf, cw, gen_alpha, gen_l, tke, tkeo, h, P, B, Px, PSTK, num, eps, avh, l_sour, q_sour, u_taus, u_taub, z0s, z0b, au, bu, cu, du, ru, qu)[source]

Advance the dynamic k-equation for one or more columns.

Return type:

None

Parameters:
pygotm.turbulence.tkeeq.step_tkeeq_single(nlev, dt, sig_k, k_min, k_ubc, k_lbc, ubc_type, lbc_type, cm0, cmsf, cw, gen_alpha, gen_l, tke, tkeo, h, P, B, Px, PSTK, num, eps, avh, l_sour, q_sour, u_taus, u_taub, z0s, z0b, au, bu, cu, du, ru, qu)

Advance the dynamic k-equation for a single column.

Return type:

None

Parameters:

Dynamic transport equation for the dissipation rate \(\varepsilon\).

Implements GOTM Section 4.7.27 (dissipationeq.F90) — solves the \(k\)–\(\varepsilon\) dissipation equation (Eq. 163):

\[\dot{\varepsilon} = \mathcal{D}_\varepsilon + \frac{\varepsilon}{k}\bigl( c_{\varepsilon 1}\,P + c_{\varepsilon 3}\,G + c_{\varepsilon x}\,P_x + c_{\varepsilon 4}\,P_s - c_{\varepsilon 2}\,\varepsilon \bigr) \comma\]

with diffusive transport (Eq. 164):

\[\mathcal{D}_\varepsilon = \frac{\partial}{\partial z} \left( \frac{\nu_t}{\sigma_\varepsilon} \frac{\partial \varepsilon}{\partial z} \right) \comma\]

and Schmidt number \(\sigma_\varepsilon\). After the transport step the turbulent length scale is recovered as:

\[l = c_{de} \frac{k^{3/2}}{\varepsilon} \point\]

Model constants (Rodi 1987)

The default constants for the \(k\)–\(\varepsilon\) model are:

Constant

Value

\(c_\mu^0\)

0.5477

\(\sigma_k\)

1.0

\(\sigma_\varepsilon\)

1.3

\(c_{\varepsilon 1}\)

1.44

\(c_{\varepsilon 2}\)

1.92

Length-scale limiter

When length_lim is active, Galperin et al. (1988) stability criterion:

\[\varepsilon \ge \frac{c_{de}}{\sqrt{2}\,\gamma} k \sqrt{\max(N^2, 0)}\]

is enforced after the transport step, preventing excessive length scales in stably stratified flows.

Boundary conditions

Boundary conditions follow logarithmic-layer theory (Dirichlet):

\[\varepsilon = \frac{c_{de}\, k^{3/2}}{\kappa\,(z + z_0)} \comma\]

or Craig–Banner wave-breaking injection.

Author (original Fortran): Lars Umlauf.

class pygotm.turbulence.dissipationeq.DissipationEquationWorkspace(nlev, *, n_cols=None)[source]

Bases: ColumnWorkspace

Workspace arrays for the translated epsilon equation.

Parameters:

  • nlev (int)

  • n_cols (int | None)

tke: ndarray
tkeo: ndarray
eps: ndarray
L: ndarray
h: ndarray
NN: ndarray
SS: ndarray
P: ndarray
B: ndarray
Px: ndarray
PSTK: ndarray
num: ndarray
avh: ndarray
sig_eff: ndarray
l_sour: ndarray
q_sour: ndarray
u_taus: ndarray
u_taub: ndarray
z0s: ndarray
z0b: ndarray
au: ndarray
bu: ndarray
cu: ndarray
du: ndarray
ru: ndarray
qu: ndarray
pygotm.turbulence.dissipationeq.step_dissipationeq(batch_size, nlev, dt, ce1, ce2, ce3plus, ce3minus, cex, ce4, cm0, cde, kappa, galp, sig_k, sig_e, sig_e0, sig_peps, length_lim, eps_min, psi_ubc, psi_lbc, ubc_type, lbc_type, cmsf, cw, gen_alpha, gen_l, tke, tkeo, eps, L, h, NN, SS, P, B, Px, PSTK, num, avh, sig_eff, l_sour, q_sour, u_taus, u_taub, z0s, z0b, au, bu, cu, du, ru, qu)[source]

Advance the dynamic epsilon-equation (batch).

Return type:

None

Parameters:
pygotm.turbulence.dissipationeq.step_dissipationeq_single(nlev, dt, ce1, ce2, ce3plus, ce3minus, cex, ce4, cm0, cde, kappa, galp, sig_k, sig_e, sig_e0, sig_peps, length_lim, eps_min, psi_ubc, psi_lbc, ubc_type, lbc_type, cmsf, cw, gen_alpha, gen_l, tke, tkeo, eps, L, h, NN, SS, P, B, Px, PSTK, num, avh, sig_eff, l_sour, q_sour, u_taus, u_taub, z0s, z0b, au, bu, cu, du, ru, qu)

Advance the dynamic epsilon-equation (single column).

Return type:

None

Parameters:

Dynamic transport equation for the specific dissipation rate \(\omega\).

Implements GOTM Section 4.7.28 (genericeq.F90, k–ω branch) — solves the \(k\)–\(\omega\) model of Umlauf et al. (2003). The inverse turbulence time scale is

\[\omega = \frac{\sqrt{k}}{c_\mu^0 l} \comma\]

and its transport equation (Eq. 168) is:

\[\dot{\omega} = \mathcal{D}_\omega + \frac{\omega}{k}\bigl( c_{\omega 1}\,P + c_{\omega 3}\,G + c_{\omega x}\,P_x + c_{\omega 4}\,P_s - c_{\omega 2}\,\varepsilon \bigr) \comma\]

with diffusive transport:

\[\mathcal{D}_\omega = \frac{\partial}{\partial z}\left( \frac{\nu_t}{\sigma_\omega} \frac{\partial \omega}{\partial z} \right) \comma\]

and Schmidt number \(\sigma_\omega\). After solving the \(\omega\) transport, the dissipation rate and turbulent length scale are recovered as:

\[\varepsilon = (c_\mu^0)^4 k\,\omega, \quad l = c_{de} \frac{k^{3/2}}{\varepsilon} \point\]

An optional Galperin et al. (1988) length-scale limiter can be activated via length_lim, restricting \(\varepsilon\) from below in stably stratified flows.

The k–ω model constants (Umlauf et al. 2003) correspond to GLS exponents \(p = 0\), \(m = 1/2\), \(n = -1\) in the GLS framework (Eq. 165: \(\psi = (c_\mu^0)^p k^m l^n\)).

Author (original Fortran): Lars Umlauf.

class pygotm.turbulence.omegaeq.OmegaEquationWorkspace(nlev, *, n_cols=None)[source]

Bases: ColumnWorkspace

Workspace arrays for the translated omega equation.

Parameters:

  • nlev (int)

  • n_cols (int | None)

tke: ndarray
tkeo: ndarray
eps: ndarray
L: ndarray
h: ndarray
NN: ndarray
SS: ndarray
P: ndarray
B: ndarray
Px: ndarray
PSTK: ndarray
num: ndarray
u_taus: ndarray
u_taub: ndarray
z0s: ndarray
z0b: ndarray
omega: ndarray
avh: ndarray
l_sour: ndarray
q_sour: ndarray
au: ndarray
bu: ndarray
cu: ndarray
du: ndarray
ru: ndarray
qu: ndarray
pygotm.turbulence.omegaeq.step_omegaeq(batch_size, nlev, dt, cw1, cw2, cw3plus, cw3minus, cwx, cw4, sig_w, cm0, kappa, cde, galp, length_lim, eps_min, psi_ubc, psi_lbc, ubc_type, lbc_type, sig_k, cmsf, cw, gen_alpha, gen_l, tke, tkeo, eps, L, h, NN, SS, P, B, Px, PSTK, num, u_taus, u_taub, z0s, z0b, omega, avh, l_sour, q_sour, au, bu, cu, du, ru, qu)[source]

Advance the dynamic omega-equation for one or more columns.

Return type:

None

Parameters:
pygotm.turbulence.omegaeq.step_omegaeq_single(nlev, dt, cw1, cw2, cw3plus, cw3minus, cwx, cw4, sig_w, cm0, kappa, cde, galp, length_lim, eps_min, psi_ubc, psi_lbc, ubc_type, lbc_type, sig_k, cmsf, cw, gen_alpha, gen_l, tke, tkeo, eps, L, h, NN, SS, P, B, Px, PSTK, num, u_taus, u_taub, z0s, z0b, omega, avh, l_sour, q_sour, au, bu, cu, du, ru, qu)

Advance the dynamic omega-equation for a single column.

Return type:

None

Parameters:

Turbulence shear and buoyancy production.

Implements GOTM Sections 4.3 and 4.7.20 — computes the mechanical shear production \(P\), buoyancy production \(G\), buoyancy-variance production \(P_b\), and extra production \(P_x\).

Shear production (Eq. 146):

\[P = \nu_t (M^2 + \alpha_w N^2) + P_{\mathrm{STK}} \comma\]

where \(\alpha_w\) is the angle between the wave and current vectors (internal wave model coupling, active only when iw_model == 1), and \(P_{\mathrm{STK}} = \nu_{\mathrm{cl}} \cdot \mathrm{SSCSTK}\) is the Stokes-drift shear production cross-term.

Buoyancy production (Eq. 147):

\[G = -\nu_t^B \left(N^2 - \tilde{\Gamma}_B\right) \comma\]

where \(\nu_t^B = \kappa_t\) is the turbulent scalar diffusivity and \(\tilde{\Gamma}_B\) is the counter-gradient buoyancy flux. This implementation sets \(\tilde{\Gamma}_B = 0\) (standard GOTM behaviour), so in practice \(G = -\kappa_t N^2\). \(G\) is negative (destruction) in stable stratification.

Buoyancy-variance production (Eq. 148):

\[P_b = -G N^2 = \nu_t^B N^4 \point\]

The total Stokes production (used in the turbulence transport equations) is:

\[P_{\mathrm{STK}} = \nu_t \cdot \mathrm{SSCSTK} + \nu_{\mathrm{cl}} \cdot \mathrm{SSSTK} \comma\]

where SSCSTK and SSSTK are the Stokes-drift cross-shear and Stokes-shear-squared terms from pygotm.meanflow.shear.

Author (original Fortran): Lars Umlauf.

class pygotm.turbulence.production.ProductionWorkspace(nlev, *, n_cols=None)[source]

Bases: ColumnWorkspace

Workspace arrays for production kernels.

Parameters:

  • nlev (int)

  • n_cols (int | None)

pygotm.turbulence.production.step_production(batch_size, nlev, iw_model, alpha, has_xP, has_sscstk, has_ssstk, NN, SS, xP, SSCSTK, SSSTK, num, nuh, nucl, P, B, Pb, Px, PSTK)[source]

Update turbulence production terms (batch).

Return type:

None

Parameters:
pygotm.turbulence.production.step_production_single(nlev, iw_model, alpha, has_xP, has_sscstk, has_ssstk, NN, SS, xP, SSCSTK, SSSTK, num, nuh, nucl, P, B, Pb, Px, PSTK)

Update turbulence production terms (single column).

Return type:

None

Parameters:

Dimensionless stability parameters \(\alpha_M\), \(\alpha_N\), \(\alpha_b\).

Implements GOTM Section 4.7.21 — updates the non-dimensional shear and buoyancy parameters used by the algebraic stress models (Eq. 149):

\[\alpha_M = \frac{k^2}{\varepsilon^2} M^2, \quad \alpha_N = \frac{k^2}{\varepsilon^2} N^2, \quad \alpha_b = \frac{k}{\varepsilon^2} \langle b'^2 \rangle \comma\]

where \(k/\varepsilon\) is the turbulence time scale \(\tau\). In the code, at stores \(\alpha_b = (k/\varepsilon)(k_b/\varepsilon)\).

When Stokes drift is active, two additional parameters are computed:

\[\alpha_v = \frac{k^2}{\varepsilon^2} \cdot \mathrm{SSCSTK}, \quad \alpha_w = \frac{k^2}{\varepsilon^2} \cdot \mathrm{SSSTK} \comma\]

stored in av and aw respectively.

All parameters are floored at \(10^{-10}\) to prevent numerical singularities in the stability-function evaluations.

Author (original Fortran): Lars Umlauf.

class pygotm.turbulence.alpha_mnb.AlphaMNBWorkspace(nlev, *, n_cols=None)[source]

Bases: ColumnWorkspace

Workspace arrays for alpha_M, alpha_N, alpha_B, and Stokes terms.

Parameters:

  • nlev (int)

  • n_cols (int | None)

pygotm.turbulence.alpha_mnb.step_alpha_mnb(batch_size, nlev, has_sscstk, has_ssstk, tke, eps, kb, NN, SS, SSCSTK, SSSTK, as_, an, at, av, aw)[source]

Update alpha_M, alpha_N, alpha_b, and optional Stokes terms (batch).

Return type:

None

Parameters:
pygotm.turbulence.alpha_mnb.step_alpha_mnb_single(nlev, has_sscstk, has_ssstk, tke, eps, kb, NN, SS, SSCSTK, SSSTK, as_, an, at, av, aw)

Update alpha_M, alpha_N, alpha_b, and optional Stokes terms (single column).

Return type:

None

Parameters:

Local weak-equilibrium stability functions \(c_\mu\) and \(c_\mu'\).

Implements GOTM Section 4.7.38 (cmue_c.F90) — computes the stability functions under the local weak-equilibrium assumption (Canuto et al. 2001; Cheng et al. 2002). Two simplifications are applied relative to the non-local variant:

  1. The buoyancy variance \(k_b = \varepsilon_b\) (local production–dissipation balance) eliminates the \(\bar{T}\) dependence from Eq. 71.

  2. The \(\Gamma\)-term in Eq. 74 drops out, so the solution is characterised by \(c_\mu\) and \(c_\mu'\) only.

The result gives the eddy viscosity and diffusivity:

\[\nu_t = c_\mu \frac{k^2}{\varepsilon}, \quad \kappa_t = c_\mu' \frac{k^2}{\varepsilon} \point\]

The denominator and numerators are polynomials (Eqs. 191–194):

\[\begin{split}D &= d_0 + d_1 \bar{N}^2 + d_2 \bar{S}^2 + d_3 \bar{N}^2 \bar{S}^2 + d_4 \bar{N}^4 + d_5 \bar{S}^4 \comma \\ N_n &= n_0 + n_1 \bar{N}^2 + n_2 \bar{S}^2 \comma \\ N_b &= n_{b0} + n_{b1} \bar{N}^2 + n_{b2} \bar{S}^2 \point\end{split}\]

Here \(\bar{S}^2 = \alpha_M\) and \(\bar{N}^2 = \alpha_N\). The coefficients \(d_i\), \(n_i\), \(n_{bi}\) are derived from the model constants \(a_1, a_2, a_3, a_5, a_{t1}, a_{t2}, a_{t3}, a_{t5}\) and \(\mathcal{N} = c_1/2\), \(\mathcal{N}_b = c_{t1}\) (see code for explicit expressions).

A clipping on \(\alpha_N\) is applied (anLimitFact = 0.5), and a corresponding limiter on \(\alpha_M\) prevents negative stability functions in strongly sheared flows.

Author (original Fortran): Lars Umlauf.

class pygotm.turbulence.cmue_c.CmueCWorkspace(nlev, *, n_cols=None)[source]

Bases: ColumnWorkspace

Workspace arrays for local weak-equilibrium stability functions.

Parameters:

  • nlev (int)

  • n_cols (int | None)

as_: ndarray
an: ndarray
cmue1: ndarray
cmue2: ndarray
pygotm.turbulence.cmue_c.step_cmue_c(batch_size, nlev, cm0, cc1, ct1, a1, a2, a3, a5, at1, at2, at3, at5, as_, an, cmue1, cmue2)[source]

Update the local weak-equilibrium stability functions (batch).

Return type:

None

Parameters:
pygotm.turbulence.cmue_c.step_cmue_c_single(nlev, cm0, cc1, ct1, a1, a2, a3, a5, at1, at2, at3, at5, as_, an, cmue1, cmue2)

Update the local weak-equilibrium stability functions (single column).

Return type:

None

Parameters:

Quasi-equilibrium stability functions \(c_\mu\) and \(c_\mu'\).

Implements GOTM Section 4.7.39 (cmue_d.F90) — computes the stability functions under the quasi-equilibrium (QE) assumption (Umlauf and Burchard 2005). Compared to pygotm.turbulence.cmue_c, the quasi-equilibrium closure retains an additional constraint: the TKE equilibrium condition (Eq. 195):

\[\frac{P + G}{\varepsilon} = \hat{c}_\mu(\alpha_M, \alpha_N)\,\alpha_M - \hat{c}_\mu'(\alpha_M, \alpha_N)\,\alpha_N = 1 \comma\]

where (149) has been used. This is an implicit relation for \(\alpha_M\) as a function of \(\alpha_N\). The solution \(\alpha_M(\alpha_N)\) is a quadratic polynomial that is solved analytically. The resulting \(\alpha_M\) is substituted back into the stability functions from Section 4.7.38 (pygotm.turbulence.cmue_c).

For negative \(\alpha_N\) (convection) the QE \(\alpha_M\) may become negative; the limiter anLimitFact = 0.5 prevents this from happening (see Umlauf and Burchard 2005).

The quasi-equilibrium variant is more accurate than the local variant in strongly stratified flows, but slightly more expensive due to the additional algebraic solve.

Author (original Fortran): Lars Umlauf.

class pygotm.turbulence.cmue_d.CmueDWorkspace(nlev, *, n_cols=None)[source]

Bases: ColumnWorkspace

Workspace arrays for quasi-equilibrium stability functions.

Parameters:

  • nlev (int)

  • n_cols (int | None)

as_: ndarray
an: ndarray
cmue1: ndarray
cmue2: ndarray
pygotm.turbulence.cmue_d.step_cmue_d(batch_size, nlev, cm0, cc1, ct1, a1, a2, a3, a5, at1, at2, at3, at5, as_, an, cmue1, cmue2)[source]

Update the quasi-equilibrium stability functions (batch).

Return type:

None

Parameters:
pygotm.turbulence.cmue_d.step_cmue_d_single(nlev, cm0, cc1, ct1, a1, a2, a3, a5, at1, at2, at3, at5, as_, an, cmue1, cmue2)

Update the quasi-equilibrium stability functions (single column).

Return type:

None

Parameters:

Algebraic closure for buoyancy variance \(k_b\).

Implements GOTM Section 4.7.30 (kbalgebraic.F90) — computes the buoyancy variance \(k_b = \langle b'^2 \rangle / 2\) under the algebraic equilibrium assumption.

The equilibrium condition \(P_b = \varepsilon_b\) (Eq. 171) is assumed, where \(P_b\) is the buoyancy-variance production. Using the definition of the time-scale ratio \(r = c_b\) (Eq. 66), this gives (Eq. 172):

\[k_b = \frac{k_{b\varepsilon}}{k\varepsilon} P_b = r \frac{k}{\varepsilon} P_b = c_b \frac{k}{\varepsilon} P_b \point\]

In the code, the coefficient \(c_b\) is passed as ctt (the scalar turbulence time-scale ratio).

The result is clipped to kb_min to prevent negative values.

Author (original Fortran): Lars Umlauf.

class pygotm.turbulence.kbalgebraic.KBAlgebraicWorkspace(nlev, *, n_cols=None)[source]

Bases: ColumnWorkspace

Workspace arrays for the algebraic buoyancy-variance closure.

Parameters:

  • nlev (int)

  • n_cols (int | None)

tke: ndarray
eps: ndarray
kb: ndarray
Pb: ndarray
pygotm.turbulence.kbalgebraic.step_kbalgebraic(batch_size, nlev, ctt, kb_min, tke, eps, kb, Pb)[source]

Advance the algebraic buoyancy-variance closure (batch).

Return type:

None

Parameters:
pygotm.turbulence.kbalgebraic.step_kbalgebraic_single(nlev, ctt, kb_min, tke, eps, kb, Pb)

Advance the algebraic buoyancy-variance closure (single column).

Return type:

None

Parameters:

Algebraic closure for buoyancy-variance dissipation \(\varepsilon_b\).

Implements GOTM Section 4.7.32 (epsbalgebraic.F90) — computes the molecular rate of destruction of buoyancy variance \(\varepsilon_b\) under the algebraic equilibrium assumption.

Assuming a constant time-scale ratio \(r = c_b\) (Eq. 66), the production–dissipation balance for \(k_b\) gives (Eq. 179):

\[\varepsilon_b = \frac{1}{c_b} \frac{\varepsilon}{k} k_b \point\]

In the code, \(c_b\) is passed as ctt (scalar turbulence time-scale ratio), so:

\[\varepsilon_b = \frac{1}{c_{tt}} \frac{\varepsilon}{k} k_b \point\]

The result is clipped to epsb_min.

Author (original Fortran): Lars Umlauf.

class pygotm.turbulence.epsbalgebraic.EpsBAlgebraicWorkspace(nlev, *, n_cols=None)[source]

Bases: ColumnWorkspace

Workspace arrays for the algebraic buoyancy-destruction closure.

Parameters:

  • nlev (int)

  • n_cols (int | None)

tke: ndarray
eps: ndarray
kb: ndarray
epsb: ndarray
pygotm.turbulence.epsbalgebraic.step_epsbalgebraic(batch_size, nlev, ctt, epsb_min, tke, eps, kb, epsb)[source]

Advance the algebraic buoyancy-destruction closure (batch).

Return type:

None

Parameters:
pygotm.turbulence.epsbalgebraic.step_epsbalgebraic_single(nlev, ctt, epsb_min, tke, eps, kb, epsb)

Advance the algebraic buoyancy-destruction closure (single column).

Return type:

None

Parameters:

Algebraic velocity-variance components \(\langle u'^2\rangle\), \(\langle v'^2\rangle\), \(\langle w'^2\rangle\).

Implements GOTM Section 4.7.33 (variances.F90) — derives the three diagonal Reynolds-stress components from algebraic expressions (Eq. 180) using the solution of the second-moment closure Eqs. 70 and 74:

\[\begin{split}\frac{\langle u'^2 \rangle}{k} &= \frac{2}{3} + \frac{1}{\mathcal{N}\varepsilon} \left[ \left(\frac{a_2}{3} + a_3\right) \nu_t \left(\frac{\partial U}{\partial z}\right)^2 - \frac{2 a_2}{3} \nu_t \left(\frac{\partial V}{\partial z}\right)^2 - \frac{4}{3} a_5 G \right] \comma \\ \frac{\langle v'^2 \rangle}{k} &= \frac{2}{3} + \frac{1}{\mathcal{N}\varepsilon} \left[ \left(\frac{a_2}{3} + a_3\right) \nu_t \left(\frac{\partial V}{\partial z}\right)^2 - \frac{2 a_2}{3} \nu_t \left(\frac{\partial U}{\partial z}\right)^2 - \frac{4}{3} a_5 G \right] \comma \\ \frac{\langle w'^2 \rangle}{k} &= \frac{2}{3} + \frac{1}{\mathcal{N}\varepsilon} \left[ \left(\frac{a_2}{3} - a_3\right)(P + P_x) + \frac{8}{3} a_5 G \right] \comma\end{split}\]

where \(\mathcal{N} = c_1/2\) is the rapid-return-to-isotropy parameter. In the code: n_value = 0.5 * cc1. The directional shear components SSU and SSV enter the \(u'\) and \(v'\) variance expressions.

The variances fall back to the isotropic value \(2k/3\) when \(\varepsilon \le 0\) or \(\mathcal{N} \le 0\).

Author (original Fortran): Lars Umlauf.

class pygotm.turbulence.variances.VariancesWorkspace(nlev, *, n_cols=None)[source]

Bases: ColumnWorkspace

Workspace arrays for algebraic velocity variances.

Parameters:

  • nlev (int)

  • n_cols (int | None)

tke: ndarray
eps: ndarray
P: ndarray
B: ndarray
Px: ndarray
num: ndarray
SSU: ndarray
SSV: ndarray
uu: ndarray
vv: ndarray
ww: ndarray
pygotm.turbulence.variances.step_variances(batch_size, nlev, cc1, ct1, a2, a3, a5, tke, eps, P, B, Px, num, SSU, SSV, uu, vv, ww)[source]

Update the algebraic velocity variances (batch).

Return type:

None

Parameters:
pygotm.turbulence.variances.step_variances_single(nlev, cc1, ct1, a2, a3, a5, tke, eps, P, B, Px, num, SSU, SSV, uu, vv, ww)

Update the algebraic velocity variances (single column).

Return type:

None

Parameters:

Internal-wave background mixing.

Implements GOTM Section 4.7.45 (internal_wave.F90) — imposes an eddy viscosity and diffusivity characteristic of internal wave activity and shear instability when the turbulence kinetic energy is at or below a threshold klimiw. Following Kantha and Clayson (1994), when \(k < k_{\mathrm{limiw}}\) the mixing coefficients are set to empirical values (Eq. 204):

\[\nu_t = \nu_t^{IW} + \nu_t^{SI}, \quad \kappa_t = \kappa_t^{IW} + \kappa_t^{SI} \comma\]

where the internal-wave background values are (Eq. 205):

\[\nu_t^{IW} = 10^{-4}\,\mathrm{m^2\,s^{-1}}, \quad \kappa_t^{IW} = 5 \times 10^{-5}\,\mathrm{m^2\,s^{-1}} \comma\]

and the shear-instability parts depend on the gradient Richardson number \(R_i = N^2 / M^2\):

\[\nu_t^{SI} = \kappa_t^{SI} = 0, \qquad R_i > 0.7 \quad (\text{Eq. 206})\]
\[\nu_t^{SI} = \kappa_t^{SI} = 5 \times 10^{-3} \left(1 - \left(\frac{R_i}{0.7}\right)^2\right)^3, \quad 0 < R_i < 0.7 \quad (\text{Eq. 207})\]
\[\nu_t^{SI} = \kappa_t^{SI} = 5 \times 10^{-3}, \qquad R_i < 0 \quad (\text{Eq. 208})\]

All diffusivities are in \(\mathrm{m^2\,s^{-1}}\).

This model is activated by iw_model = 2. The contribution is added to \(\nu_t\) and \(\kappa_t\) after the turbulence-closure update.

Authors (original Fortran): Karsten Bolding, Hans Burchard, Manuel Ruiz Villarreal.

class pygotm.turbulence.internal_wave.InternalWaveWorkspace(nlev, *, n_cols=None)[source]

Bases: ColumnWorkspace

Workspace arrays for internal-wave mixing.

Parameters:

  • nlev (int)

  • n_cols (int | None)

tke: ndarray
num: ndarray
nuh: ndarray
NN: ndarray
SS: ndarray
pygotm.turbulence.internal_wave.step_internal_wave(batch_size, nlev, iw_model, klimiw, rich_cr, numiw, nuhiw, numshear, tke, num, nuh, NN, SS)[source]

Apply the Kantha-Clayson internal-wave mixing fallback (batch).

Return type:

None

Parameters:
pygotm.turbulence.internal_wave.step_internal_wave_single(nlev, iw_model, klimiw, rich_cr, numiw, nuhiw, numshear, tke, num, nuh, NN, SS)

Apply the Kantha-Clayson internal-wave mixing fallback (single column).

Return type:

None

Parameters: