Mean Flow

Modules implementing the 1D mean-flow equations. Each module corresponds to one GOTM Fortran subroutine and is compiled with @numba.njit.

Mean-flow physics modules.

Coriolis rotation of horizontal velocity.

Implements GOTM Section 3.2.4 — applies an exact rotation of the horizontal velocity vector \((U, V)\) through angle \(\omega = f \Delta t\) at each vertical level, where \(f\) is the Coriolis parameter and \(\Delta t\) is the time step. This exactly satisfies the homogeneous Coriolis system

\[\frac{\partial U}{\partial t} = fV, \quad \frac{\partial V}{\partial t} = -fU \comma\]

by applying the rotation matrix

\[\begin{split}\begin{pmatrix} U^{n+1} \\ V^{n+1} \end{pmatrix} = \begin{pmatrix} \cos\omega & \sin\omega \\ -\sin\omega & \cos\omega \end{pmatrix} \begin{pmatrix} U^n \\ V^n \end{pmatrix} \comma\end{split}\]

where \(\omega = f \Delta t\).

When Stokes drift profiles \((u_s, v_s)\) are supplied, they are added to the Eulerian velocity before rotation and subtracted after, so that only the Eulerian part of the horizontal velocity is rotated:

\[\tilde{U} = U + u_s, \quad \tilde{V} = V + v_s \comma\]

and the rotation is applied to \((\tilde{U}, \tilde{V})\).

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

pygotm.meanflow.coriolis.coriolis(state, nlev, dt, usprof=None, vsprof=None)[source]

Apply the Coriolis rotation to horizontal velocities for one column.

Return type:

None

Parameters:

Parameters

state:

MeanflowState with u, v (horizontal velocities), cori (Coriolis parameter f [rad/s]).

nlev:

Number of model layers.

dt:

Time step [s]. Rotation angle is omega = state.cori * dt.

usprof:

Stokes drift profile in x, shape (nlev+1,) [m/s]. If None, zero.

vsprof:

Stokes drift profile in y, shape (nlev+1,) [m/s]. If None, zero.

pygotm.meanflow.coriolis.step_coriolis_batch(batch_size, nlev, cosomega, sinomega, u, v, usprof, vsprof)[source]

Batch Coriolis rotation: process batch_size columns in parallel.

Return type:

None

Parameters:
pygotm.meanflow.coriolis.step_coriolis_single(nlev, cosomega, sinomega, u, v, usprof, vsprof)
Return type:

None

Parameters:

U-momentum (east–west) equation.

Implements GOTM Section 3.2.5 — advances the depth-varying eastward velocity \(U\) by one timestep using a Crank–Nicolson scheme.

Governing equation

The horizontally-averaged, incompressible U-momentum equation reads (Eq. 12):

\[\frac{\partial U}{\partial t} - f V = -g \frac{\partial \zeta}{\partial x} + \frac{\partial}{\partial z} \left[ (\nu_t + \nu) \frac{\partial U}{\partial z} - \tilde{\Gamma}_U \right] + \int_z^\eta \frac{\partial B}{\partial x}\,dz' - \frac{1}{\tau_R}(U - U_{\mathrm{obs}}) - C_f U \sqrt{U^2 + V^2} \comma\]

where \(f\) is the Coriolis parameter (handled separately by pygotm.meanflow.coriolis), \(\nu_t\) is the turbulent eddy viscosity, \(\nu\) is the molecular viscosity, \(\tilde{\Gamma}_U\) is the counter-gradient momentum flux (Langmuir/Stokes correction), \(B = -g(\rho - \rho_0)/\rho_0\) is the buoyancy, and \(C_f\) is a quadratic bottom drag coefficient.

The external (barotropic) pressure gradient from sea-surface slope is \(-g \partial\zeta/\partial x\) (active when ext_method = 0). Internal (baroclinic) pressure gradients from horizontal density gradients enter through idpdx.

Numerics

The Crank–Nicolson scheme (implicitness \(\sigma\), set to 1 for fully implicit) discretises the vertical diffusion term, producing an unconditionally stable tridiagonal system solved by the Thomas algorithm at each time step (see pygotm.util.diff_center).

Bottom drag is treated as a linearised source term \(l_1 = -C_f \sqrt{U_1^2 + V_1^2} / h_1\) applied at the lowest layer.

When w_adv_active is True, vertical advection is applied using the scheme given by w_adv_discr (upwind, Superbee, etc.) from pygotm.util.adv_center.

Author (original Fortran): Lars Umlauf.

pygotm.meanflow.uequation.uequation(state, nlev, dt, cnpar, tx, num, nucl, gamu, ext_method=0, dpdx=0.0, idpdx=None, dusdz=None, w_adv_active=False, w_adv_discr=4, vel_relax_tau=None, vel_relax_ramp=1000000000000000.0, uprof=None, plume_active=False, seagrass_active=False)[source]

Advance the U-momentum equation for one column.

Return type:

None

Parameters:

Parameters

state:

MeanflowState with h, u, uo, v, w, drag, avh, gravity, avmolu.

nlev:

Number of model layers.

dt:

Time step [s].

cnpar:

Crank–Nicolson implicitness (0=explicit, 1=fully implicit).

tx:

Surface wind stress in x [m² s⁻²] (Neumann BC at surface).

num:

Turbulent eddy viscosity [m² s⁻¹], shape (nlev+1,).

nucl:

Non-local eddy viscosity (Stokes/Langmuir correction) [m² s⁻¹], shape (nlev+1,).

gamu:

Counter-gradient momentum flux (not currently used; reserved).

ext_method:

External pressure treatment (0 = barotropic slope from dpdx).

dpdx:

Barotropic pressure gradient \(g\partial\zeta/\partial x\) [m s⁻²].

idpdx:

Baroclinic pressure gradient profile [m s⁻²], shape (nlev+1,).

dusdz:

Stokes drift shear in x [s⁻¹], shape (nlev+1,).

pygotm.meanflow.uequation.step_uequation(batch_size, nlev, dt, cnpar, avmolu, gravity, ext_method, w_adv_active, w_adv_discr, seagrass_active, plume_active, tx, dzetadx, u, uo, v, h, w, drag, num, nucl, dusdz, idpdx, uprof, tau_r, avh, q_sour, l_sour, au, bu, cu, du, ru, qu, adv_cu)[source]

Batch variant: process batch_size columns in parallel.

Return type:

None

Parameters:
pygotm.meanflow.uequation.step_uequation_single(nlev, dt, cnpar, avmolu, gravity, ext_method, w_adv_active, w_adv_discr, seagrass_active, plume_active, tx_val, dzetadx_val, u, uo, v, h, w, drag, num, nucl, dusdz, idpdx, uprof, tau_r, avh, q_sour, l_sour, au, bu, cu, du, ru, qu, adv_cu)
Return type:

None

Parameters:

V-momentum (north–south) equation.

Implements GOTM Section 3.2.6 — advances the depth-varying northward velocity \(V\) by one timestep using a Crank–Nicolson scheme identical to that of pygotm.meanflow.uequation.

Governing equation

The northward momentum equation reads (Eq. 14):

\[\frac{\partial V}{\partial t} + f U = \mathcal{D}_V -g \frac{\partial \zeta}{\partial y} + \int_z^\eta \frac{\partial B}{\partial y}\,dz' - \frac{1}{\tau_R}(V - V_{\mathrm{obs}}) - C_f V \sqrt{U^2 + V^2} \comma\]

where the diffusion operator is (Eq. 15):

\[\mathcal{D}_V = \frac{\partial}{\partial z} \left[ (\nu_t + \nu) \frac{\partial V}{\partial z} - \tilde{\Gamma}_V \right] \point\]

The Coriolis term \(+fU\) is handled separately by pygotm.meanflow.coriolis before this equation is solved. All other terms are identical in structure to the U-equation; see pygotm.meanflow.uequation for details.

Author (original Fortran): Lars Umlauf.

pygotm.meanflow.vequation.vequation(state, nlev, dt, cnpar, ty, num, nucl, gamv, ext_method=0, dpdy=0.0, idpdy=None, dvsdz=None, w_adv_active=False, w_adv_discr=4, vel_relax_tau=None, vel_relax_ramp=1000000000000000.0, vprof=None, plume_active=False)[source]

Advance the V-momentum equation for one column.

Return type:

None

Parameters:
pygotm.meanflow.vequation.step_vequation(batch_size, nlev, dt, cnpar, avmolu, gravity, ext_method, w_adv_active, w_adv_discr, plume_active, ty, dzetady, v, vo, u, h, w, drag, num, nucl, dvsdz, idpdy, vprof, tau_r, avh, q_sour, l_sour, av, bv, cv, dv, rv, qv, adv_cv)[source]

Batch variant: process batch_size columns in parallel.

Return type:

None

Parameters:
pygotm.meanflow.vequation.step_vequation_single(nlev, dt, cnpar, avmolu, gravity, ext_method, w_adv_active, w_adv_discr, plume_active, ty_val, dzetady_val, v, vo, u, h, w, drag, num, nucl, dvsdz, idpdy, vprof, tau_r, avh, q_sour, l_sour, av, bv, cv, dv, rv, qv, adv_cv)
Return type:

None

Parameters:

Temperature (potential) transport equation.

Implements GOTM Section 3.2.10 — advances potential temperature \(\Theta\) by one timestep using a Crank–Nicolson advection–diffusion solver.

Transport equation

The temperature satisfies (Eq. 27):

\[\dot{\Theta} = \mathcal{D}_\Theta - \frac{1}{\tau_R^\Theta}(\Theta - \Theta_{\mathrm{obs}}) + \frac{1}{C_p \rho_0} \frac{\partial I}{\partial z} \comma\]

where \(\mathcal{D}_\Theta\) is the turbulent diffusion operator (Eq. 28):

\[\mathcal{D}_\Theta = \frac{\partial}{\partial z} \left( (\nu_t^\Theta + \nu^\Theta) \frac{\partial \Theta}{\partial z} - \tilde{\Gamma}_\Theta \right) \comma\]

with turbulent scalar diffusivity \(\nu_t^\Theta\), molecular diffusivity \(\nu^\Theta = \nu_{\mathrm{mol}}^T\) (avmolT), and counter-gradient flux \(\tilde{\Gamma}_\Theta\) (gamh).

Short-wave radiation absorption

The radiative heating term uses the Paulson and Simpson (1977) double-exponential profile (Eq. 29):

\[I(z) = I_0 \left[ A\,e^{-z/\eta_1} + (1 - A)\,e^{-z/\eta_2}\,\mathrm{bioshade} \right] \comma\]

where \(I_0\) is the surface short-wave irradiance [W m⁻²], \(A = 0.58\) is the fraction in the red/near-infrared band, \(\eta_1 = 0.35\,\mathrm{m}\) and \(\eta_2 = 23.0\,\mathrm{m}\) are the extinction depths (defaults from Paulson & Simpson 1977), and bioshade is an optional biological shading factor.

The divergence \(\partial I/\partial z\) contributes to q_sour at each layer. The surface layer receives the residual flux absorbed in the topmost cell.

Surface boundary condition

At the sea surface a Neumann (prescribed flux) condition is applied:

\[(\nu_t^\Theta + \nu^\Theta) \frac{\partial \Theta}{\partial z}\Bigg|_{z=\zeta} = \frac{Q_{\mathrm{net}}}{C_p \rho_0} \comma\]

where \(Q_{\mathrm{net}}\) is the net non-penetrative heat flux (positive into the ocean). In the code, hflux carries the opposite sign — it follows the atmospheric convention (positive = heat loss from ocean, negative = heat gain) — so the Neumann value passed to the diffusion solver is diff_t_up = -hflux / (rho0 * cp), which is positive (warming) when hflux < 0. This matches the calling convention in both pygotm.gotm.gotm and Fortran GOTM v7, where shf = -heat_input%value before calling temperature().

A sea-ice correction suppresses the warming flux when the sea-surface temperature is at or below the saline freezing point (\(T \le -0.0575\,S\)).

Author (original Fortran): Lars Umlauf.

pygotm.meanflow.temperature.temperature(state, nlev, dt, cnpar, I_0, wflux, hflux, nuh, gamh, *, rho0, cp, A=0.58, g1=0.35, g2=23.0, Tobs=None, tau_r=None, dtdx=None, dtdy=None, w_adv_active=False, w_adv_discr=4, t_adv=False, apply_simple_ice_correction=False)[source]

Advance the temperature equation for one column.

Return type:

None

Parameters:

Parameters

state:

MeanflowState with T, S, h, w, u, v, bioshade, rad, Tobs, avh, avmolT.

nlev:

Number of model layers.

dt:

Time step [s].

cnpar:

Crank–Nicolson implicitness (1 = fully implicit).

I_0:

Surface short-wave irradiance [W m⁻²].

wflux:

Freshwater flux [m s⁻¹] (reserved; not currently used in heat equation).

hflux:

Net non-penetrative surface heat flux [W m⁻²], atmospheric convention (positive = heat loss from ocean; negative = heat gain into ocean). Identical convention to Fortran GOTM v7 shf = -heat_input%value. Converted to Neumann BC as diff_t_up = -hflux / (rho0 * cp); negative hflux produces a positive (warming) Neumann value.

nuh:

Turbulent heat diffusivity profile [m² s⁻¹], shape (nlev+1,).

gamh:

Counter-gradient heat flux \(\tilde{\Gamma}_\Theta\), shape (nlev+1,).

rho0:

Reference density [kg m⁻³].

cp:

Specific heat capacity [J kg⁻¹ K⁻¹].

A:

Fraction of short-wave in the red/near-IR band (Paulson & Simpson 1977). Default 0.58.

g1:

Extinction depth for band 1 [m]. Default 0.35 m.

g2:

Extinction depth for band 2 [m]. Default 23.0 m.

pygotm.meanflow.temperature.step_temperature(batch_size, nlev, dt, cnpar, avmolT, rho0, cp, A, g1, g2, w_adv_active, w_adv_discr, t_adv, T, S, h, w, u, v, nuh, gamh, bioshade, rad, Tobs, tau_r, i_0, diff_t_up, apply_freezing_correction, dtdx, dtdy, avh, q_sour, l_sour, au, bu, cu, du, ru, qu, adv_cu)[source]

Batch variant: process batch_size columns in parallel.

Return type:

None

Parameters:
pygotm.meanflow.temperature.step_temperature_single(nlev, dt, cnpar, avmolT, rho0, cp, A, g1, g2, w_adv_active, w_adv_discr, t_adv, T, S, h, w, u, v, nuh, gamh, bioshade, rad, Tobs, tau_r, i_0, diff_t_up, apply_freezing_correction, dtdx, dtdy, avh, q_sour, l_sour, au, bu, cu, du, ru, qu, adv_cu)
Return type:

None

Parameters:

Salinity transport equation.

Implements GOTM Section 3.2.11 — advances salinity \(S\) by one timestep using an advection–diffusion solver with Neumann surface and bottom boundary conditions.

The salinity satisfies the general scalar transport equation (Eq. 30):

\[\dot{S} = \mathcal{D}_S - \frac{1}{\tau_R^S}(S - S_{\mathrm{obs}}) \comma\]

where \(\mathcal{D}_S\) is the turbulent diffusion operator (Eq. 31):

\[\mathcal{D}_S = \frac{\partial}{\partial z} \left( (\nu_t^S + \nu^S) \frac{\partial S}{\partial z} - \tilde{\Gamma}_S \right) \comma\]

with turbulent scalar diffusivity \(\nu_t^S = \nu_t^{\Theta}\) (the same as for temperature in the current GOTM implementation), molecular diffusivity \(\nu^S\), and counter-gradient flux \(\tilde{\Gamma}_S\).

The optional \(-\tau_R^{-1}(S - S_{\mathrm{obs}})\) term relaxes the modelled salinity towards an observed profile \(S_{\mathrm{obs}}\) with time scale \(\tau_R^S\).

Surface boundary condition

At the sea surface \(z = \zeta\) a Neumann (prescribed flux) condition is applied (Eq. 32):

\[(\nu_t^S + \nu^S) \frac{\partial S}{\partial z} = S (P - E) \comma\]

where \(P - E\) is precipitation minus evaporation [m s⁻¹]. In the code this is passed as sflux = S_sfc * (P - E) [psu m s⁻¹] — positive when precipitation exceeds evaporation (P > E, dilution regime), negative when evaporation dominates (E > P, concentration regime).

Note

Sign convention vs. the GOTM manual (Eq. 32). The manual writes the surface BC as \((\nu_t^S + \nu^S)\,\partial S/\partial z = S(P-E)\) at \(z = \zeta\), where the right-hand side is the diffusive flux in the upward (+z) direction — positive means salt leaves the water column (P > E, dilution). The diff_center kernel, however, uses the opposite sign convention: the upper Neumann argument Yup is the flux entering the top cell (positive = salt gains the surface layer). Because these two conventions have opposite signs, the code negates sflux before passing it to the kernel:

diff_s_up = -sflux = -S*(P-E) = S*(E-P)

For evaporation (E > P): diff_s_up > 0 → flux entering the top cell → surface salinity increases. ✓ For precipitation (P > E): diff_s_up < 0 → flux leaving the top cell → surface salinity decreases. ✓

This matches the Fortran verbatim (DiffSup = -sflux in salinity.F90).

Both surface and bottom boundary conditions are of Neumann type; the bottom flux is zero.

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

pygotm.meanflow.salinity.salinity(state, nlev, dt, cnpar, wflux, sflux, nus, gams, *, Sobs=None, tau_r=None, dsdx=None, dsdy=None, w_adv_active=False, w_adv_discr=4, s_adv=False)[source]

Advance the salinity equation for one column.

Return type:

None

Parameters:

Parameters

state:

MeanflowState with S, h, w, u, v, Sobs, avh, avmolS.

nlev:

Number of model layers.

dt:

Time step [s].

cnpar:

Crank–Nicolson implicitness (0 = explicit, 1 = fully implicit).

wflux:

Freshwater flux at the surface [m s⁻¹] (positive = precipitation). Not used directly; kept for API symmetry with temperature().

sflux:

Virtual salinity flux S_surface · (P − E) [psu m s⁻¹]. Positive when P > E (precipitation-dominated, dilution); negative when E > P (evaporation-dominated, concentration). The kernel receives -sflux as the upper Neumann boundary condition; see the module docstring for the sign-convention note.

nus:

Turbulent salinity diffusivity profile [m² s⁻¹], shape (nlev+1,).

gams:

Counter-gradient salinity flux \(\tilde{\Gamma}_S\), shape (nlev+1,).

pygotm.meanflow.salinity.step_salinity(batch_size, nlev, dt, cnpar, avmolS, w_adv_active, w_adv_discr, s_adv, S, h, w, u, v, nus, gams, Sobs, tau_r, diff_s_up, dsdx, dsdy, avh, q_sour, l_sour, au, bu, cu, du, ru, qu, adv_cu)[source]

Batch variant: process batch_size columns in parallel.

Return type:

None

Parameters:
pygotm.meanflow.salinity.step_salinity_single(nlev, dt, cnpar, avmolS, w_adv_active, w_adv_discr, s_adv, S, h, w, u, v, nus, gams, Sobs, tau_r, diff_s_up, dsdx, dsdy, avh, q_sour, l_sour, au, bu, cu, du, ru, qu, adv_cu)
Return type:

None

Parameters:

Vertical friction and boundary-layer drag.

Implements GOTM Section 3.2.9 — computes bottom friction velocity \(u_{\tau b}\) from a logarithmic drag law and surface stress from wind forcing, then applies them as Neumann boundary conditions for the momentum equations.

Bottom roughness and friction velocity

The bottom roughness length combines hydraulically smooth, hydraulically rough, and moveable-bed (sediment) contributions (Eq. 23):

\[z_{0b} = \frac{0.1\nu}{u_{\tau b}} + 0.03 h_{0b} + z_a \comma\]

where \(\nu\) is the molecular viscosity, \(h_{0b}\) is the physical roughness height, and \(z_a\) is a moveable bed roughness. When \(\nu \le 0\) the smooth-wall term is omitted.

The bottom friction velocity is (Eq. 24):

\[u_{\tau b} = r \sqrt{U_1^2 + V_1^2} \comma\]

where \(U_1, V_1\) are the velocities in the lowest model layer, and the drag coefficient \(r\) is (Eq. 25):

\[r = \frac{\kappa}{\ln\bigl((z_{0b} + h_1/2) / z_{0b}\bigr)} \comma\]

with \(h_1\) the thickness of the lowest model layer and \(\kappa = 0.4\) the von Kármán constant. Because \(z_{0b}\) itself depends on \(u_{\tau b}\) through the smooth-wall term, the equations are iterated up to MaxItz0b times.

Surface roughness

When Charnock’s formula is activated (Eq. 26):

\[z_{0s} = \frac{\alpha_c\, u_{\tau s}^2}{g} \comma\]

where \(\alpha_c\) is the Charnock coefficient (charnock_val) and \(g\) is gravitational acceleration. Otherwise \(z_{0s}\) is set to the fixed minimum value z0s_min.

When plume_type = 1 the surface friction velocity is computed analogously to the bottom, from the top-layer velocity. Otherwise it is derived from the prescribed wind-stress components \((\tau_x, \tau_y)\):

\[u_{\tau s} = (\tau_x^2 + \tau_y^2)^{1/4} \point\]

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

class pygotm.meanflow.friction.FrictionWorkspace(nlev, *, batch_size=1)[source]

Bases: ColumnWorkspace

Batch friction workspace — profile arrays shape (batch_size, nlev+1), scalar arrays shape (batch_size,).

Parameters:
pygotm.meanflow.friction.friction(state, nlev, *, kappa=0.4, avmolu=None, tx=0.0, ty=0.0, plume_type=0, rho0=1027.0, _first=None)[source]

Update bottom roughness and compute friction velocities and drag.

Return type:

None

Parameters:
pygotm.meanflow.friction.step_friction_batch(batch_size, nlev, kappa, avmolu, rho0, gravity, h0b, z0s_min, charnock, charnock_val, calc_bottom_stress, MaxItz0b, plume_type, first, h, u, v, drag, z0b, z0s, za, u_taub, u_taubo, u_taus, taub, tx, ty)[source]

Batch variant: process batch_size columns in parallel.

Return type:

None

Parameters:
pygotm.meanflow.friction.step_friction_single(nlev, kappa, avmolu, rho0, gravity, h0b, z0s_min, charnock, charnock_val, calc_bottom_stress, MaxItz0b, plume_type, first, h, u, v, drag, z0b, z0s, za, u_taub, u_taubo, u_taus, taub, tx, ty)

Single-column friction kernel. Scalar outputs use length-1 arrays.

Return type:

None

Parameters:

Vertical shear-frequency squared.

Implements GOTM Section 3.2.13. The squared shear frequency is defined as (Eq. 36):

\[M^2 = \left( \frac{\partial U}{\partial z} \right)^2 + \left( \frac{\partial V}{\partial z} \right)^2 \point\]

Energy-conserving discretisation

The \(U\)- and \(V\)-contributions are discretised with the energy-conserving scheme of Burchard (2002) (Eq. 37). Defining the averaged velocities \(\tilde{U}_j = \tfrac{1}{2}(\hat{U}_j + U_j^n)\) (where \(\hat{U}\) is the updated and \(U^n\) is the old velocity), the discrete shear is

\[M^2_{U,j} = \frac{1}{2} \left[ \frac{(\hat{U}_{j+1}-\hat{U}_j)(\hat{U}_{j+1}-U_j^n)}{h_{j+1/2}\,h_j} + \frac{(\hat{U}_{j+1}-\hat{U}_j)(U^n_{j+1}-\hat{U}_j)}{h_{j+1/2}\,h_{j+1}} \right] \comma\]

and analogously for \(M^2_V\). The scheme guarantees that no spurious mean kinetic energy is generated in the conversion to turbulent kinetic energy.

Stokes drift contributions

When Stokes drift shear profiles \((\partial u_s/\partial z, \partial v_s/\partial z)\) are supplied, two additional arrays are computed:

  • SSCSTK: cross-production term \((\partial u_s/\partial z)(\partial U/\partial z) + (\partial v_s/\partial z)(\partial V/\partial z)\), used in the Stokes-drift shear production in the turbulence equations.

  • SSSTK: squared Stokes shear \((\partial u_s/\partial z)^2 + (\partial v_s/\partial z)^2\), entering the extra Stokes production term.

The resulting \(M^2\) drives shear production in all two-equation turbulence closures.

Author (original Fortran): Lars Umlauf.

pygotm.meanflow.shear.shear(state, nlev, cnpar, dusdz=None, dvsdz=None)[source]

Compute the shear-frequency squared (M²) at layer interfaces.

Updates state.SS, state.SSU, state.SSV, state.SSCSTK, and state.SSSTK in-place using the energy-conserving Burchard (2002) discretisation.

Return type:

None

Parameters:

Parameters

state:

MeanflowState with h, u, v, uo, vo, SS, SSU, SSV, SSCSTK, SSSTK. All arrays have shape (nlev+1,); index 0 = seabed, nlev = surface.

nlev:

Number of model layers.

cnpar:

Numerical implicitness parameter (0 = explicit, 1 = fully implicit; 0.5 for Crank-Nicolson). Controls weighting between old and new time-level velocity differences in the shear production formula.

dusdz:

Stokes drift shear in x, shape (nlev+1,) [s⁻¹]. Defaults to zeros.

dvsdz:

Stokes drift shear in y, shape (nlev+1,) [s⁻¹]. Defaults to zeros.

pygotm.meanflow.shear.step_shear_single(nlev, cnpar, h, u, v, uo, vo, dusdz, dvsdz, SS, SSU, SSV, SSCSTK, SSSTK)[source]

Compute shear-frequency squared for one 0:nlev column.

Return type:

None

Parameters: