Utilities

Shared numerical utilities used across all pyGOTM physics kernels. Each submodule is a direct Python translation of the corresponding GOTM Fortran source file under src/.

Submodules at a glance

Module

Purpose

tridiagonal

Thomas algorithm (tridiagonal) solver

adv_center

Vertical advection — cell-centred tracers

diff_center

Vertical diffusion — cell-centred variables

diff_face

Vertical diffusion — face-centred variables

density

Equation of state (density, α, β)

gsw

Numba-compiled TEOS-10 GSW library

time

Time management (Julian day + seconds)

convert_fluxes

Flux conversion for the KPP turbulence model

gridinterpol

Grid interpolation from observation to model grid

lagrange

Lagrangian particle random walk

ode_solvers

ODE solvers for biogeochemical models (11 methods)

ode_solvers_template

Alternate RK2/RK4 implementations (template variant)

gotm_version

GOTM version constants

compilation

Compilation metadata stub

Selector Constants

util defines the integer selector constants shared by the advection, diffusion, and advection-boundary-condition routines.

Advection schemes (method argument to adv_center()):

Constant

Value

Scheme

CENTRAL

-1

Central differencing

UPSTREAM

1

First-order upwind

P1

2

First-order upwind (P1)

P2

3

Second-order unbounded

Superbee

4

Superbee limiter (Roe 1986)

MUSCL

5

MUSCL (van Leer 1979)

P2_PDM

6

P2 with Positive Definite Method (Pietrzak 1998)

SPLMAX13

13

SPLMAX13 (Pietrzak 1998)

Diffusion boundary conditions (bc_up / bc_down):

Constant

Value

Meaning

Dirichlet

0

Prescribe value at boundary

Neumann

1

Prescribe flux at boundary

Advection boundary conditions (bc_up / bc_down):

Constant

Value

Meaning

flux

1

Prescribed flux [tracer·m/s]

value

2

Prescribed tracer value

oneSided

3

One-sided upwind (outflow only)

zeroDivergence

4

Zero-divergence extrapolation

Tridiagonal Solver

Translation of mtridiagonal.F90. Provides a Numba-JIT Thomas algorithm (simplified Gaussian elimination) for tridiagonal linear systems of the form:

\[b_i\,y_i + c_i\,y_{i+1} + a_i\,y_{i-1} = d_i, \quad i = f_i,\ldots,l_t\]

The main diagonal is stored in bu, upper diagonal in au, lower diagonal in cu, and the right-hand side in du. Work arrays ru and qu hold intermediate values.

All diffusion, momentum, and turbulence-closure routines in pyGOTM call tridiagonal() as their final linear-system step.

Workspace classes pre-allocate all six arrays (shape (nlev+1,) for a single column, shape (batch_size, nlev+1) for batch runs):

Tridiagonal (Thomas algorithm) solver — translation of mtridiagonal.F90.

Provides a Numba-JIT Thomas algorithm solver for tridiagonal linear systems, used by the diffusion, turbulence-closure, and momentum routines throughout pyGOTM.

The main diagonal is stored in bu, the upper diagonal in au, the lower diagonal in cu, and the right-hand side in du. Work arrays ru and qu hold intermediate values during forward substitution and back substitution respectively.

Public interface: init_tridiagonal(), tridiagonal(), clean_tridiagonal(), TridiagonalWorkspace, TridiagonalBatchWorkspace.

Original FORTRAN authors: Hans Burchard, Karsten Bolding.

class pygotm.util.tridiagonal.TridiagonalBatchWorkspace(nlev, batch_size)[source]

Bases: ColumnWorkspace

Batch tridiagonal workspace — arrays shape (batch_size, nlev+1).

Parameters:
class pygotm.util.tridiagonal.TridiagonalWorkspace(nlev)[source]

Bases: ColumnWorkspace

Single-column tridiagonal workspace — arrays shape (nlev+1,).

Parameters:

nlev (int)

pygotm.util.tridiagonal.clean_tridiagonal(workspace)[source]

No-op — NumPy workspaces are garbage-collected normally.

Return type:

None

Parameters:

workspace (TridiagonalWorkspace)

pygotm.util.tridiagonal.init_tridiagonal(nlev)[source]

Allocate the tridiagonal coefficients and Thomas work arrays.

Return type:

TridiagonalWorkspace

Parameters:

nlev (int)

pygotm.util.tridiagonal.tridiagonal(au, bu, cu, du, ru, qu, value, fi, lt)[source]

Solve a tridiagonal system with the Thomas algorithm (simplified Gaussian elimination).

A linear equation with tridiagonal matrix structure is solved here. The main diagonal is stored in bu, the upper diagonal in au, the lower diagonal in cu, and the right-hand side in du. Indices run from fi to lt inclusive.

Return type:

None

Parameters:

Vertical Advection

Translation of adv_center.F90. Solves the one-dimensional advection equation for tracers defined at cell centres.

Two conservation modes are supported:

  • Conservative (CONSERVATIVE = 1): \(\partial_t Y = -\partial_z(wY)\). Used for settling or rising tracers such as sediment or phytoplankton.

  • Non-conservative (NON_CONSERVATIVE = 0): \(\partial_t Y = -w\,\partial_z Y\). Used when the water column has a prescribed net vertical velocity.

The face flux is reconstructed with one of the eight slope-limiter schemes listed in the selector-constants table above (UPSTREAM through SPLMAX13). Sub-stepping is applied automatically when \(\max(|w|\Delta t / \Delta z) > 1\), up to 100 sub-steps per time step.

Workspace classes:

Vertical advection for cell-centred tracers — translation of adv_center.F90.

Solves the one-dimensional advection equation for variables defined at cell centres. Two formulations are supported:

  • Conservative (CONSERVATIVE = 1): \(\partial_t Y = -\partial_z(wY)\). Used for settling or rising tracers (e.g. sediment, phytoplankton).

  • Non-conservative (NON_CONSERVATIVE = 0): \(\partial_t Y = -w\partial_z Y\). Used when the water column has a prescribed net vertical velocity.

The advective face flux is reconstructed with one of seven slope-limiter schemes selected by the method integer:

Constant

Value

Scheme

UPSTREAM

1

First-order upwind

P1

2

First-order upwind (P1 — same as UPSTREAM in this implementation)

P2

3

Second-order unbounded (may produce over/undershoots)

Superbee

4

Superbee limiter (Roe 1986)

MUSCL

5

MUSCL (van Leer 1979)

P2_PDM

6

P2 with Positive Definite Method limiter (Pietrzak 1998)

SPLMAX13

13

SPLMAX13 (Pietrzak 1998)

Sub-stepping is applied when \(\max(|w|\Delta t / \Delta z) > 1\), up to _ITMAX = 100 sub-steps per timestep.

Original FORTRAN author: Lars Umlauf.

class pygotm.util.adv_center.AdvectionBatchWorkspace(nlev, batch_size)[source]

Bases: ColumnWorkspace

Batch advection workspace — cu has shape (batch_size, nlev+1).

Parameters:
class pygotm.util.adv_center.AdvectionWorkspace(nlev)[source]

Bases: ColumnWorkspace

Single-column advection workspace — cu has shape (nlev+1,).

Parameters:

nlev (int)

pygotm.util.adv_center.adv_center(nlev, dt, h, ho, ww, bc_up, bc_down, y_up, y_down, method, mode, y, cu)[source]

Advance a single-column tracer on cell centers using GOTM advection.

Return type:

None

Parameters:
pygotm.util.adv_center.adv_center_batch(batch_size, nlev, dt, h, ho, ww, bc_up, bc_down, y_up, y_down, method, mode, y, cu)[source]

Batch variant: process batch_size columns in parallel with numba.prange.

Return type:

None

Parameters:
pygotm.util.adv_center.clean_adv_center(workspace)[source]

No-op — NumPy workspaces are garbage-collected normally.

Return type:

None

Parameters:

workspace (AdvectionWorkspace)

pygotm.util.adv_center.init_adv_center(nlev)[source]

Allocate the temporary face-flux array used by adv_center.

Return type:

AdvectionWorkspace

Parameters:

nlev (int)

Vertical Diffusion (Cell Centres)

Translation of diff_center.F90. Solves the one-dimensional diffusion equation with optional source terms and relaxation for variables defined at cell centres:

\[\frac{\partial Y}{\partial t} = \frac{\partial}{\partial z}\!\left(\nu_Y \frac{\partial Y}{\partial z}\right) - \frac{Y - Y_{\mathrm{obs}}}{\tau_R} + Y\,L_{\mathrm{sour}} + Q_{\mathrm{sour}}\]

The diffusivity \(\nu_Y\) is defined at cell faces. The diffusion term, linear source \(L_{\mathrm{sour}}\), and relaxation are treated implicitly with Crank–Nicolson implicitness cnpar; the constant source \(Q_{\mathrm{sour}}\) is explicit. Relaxation is only applied where tau_r[i] < 1e10.

For non-negative concentrations (posconc = 1), negative Neumann boundary fluxes are linearised following Patankar (1980) to preserve positivity.

Vertical diffusion for cell-centred variables — translation of diff_center.F90.

Solves the one-dimensional diffusion equation with optional source terms and relaxation toward observed values:

\[\frac{\partial Y}{\partial t} = \frac{\partial}{\partial z}\!\left(\nu_Y \frac{\partial Y}{\partial z}\right) - \frac{Y - Y_{\mathrm{obs}}}{\tau_R} + Y\,L_{\mathrm{sour}} + Q_{\mathrm{sour}}\]

The diffusivity \(\nu_Y\) is defined at cell faces. The diffusion term, linear source \(L_{\mathrm{sour}}\), and relaxation term are treated implicitly with Crank–Nicolson implicitness cnpar; the constant source \(Q_{\mathrm{sour}}\) is explicit. Relaxation is only applied where tau_r[i] < 1e10.

Boundary conditions (bc_up, bc_down) are Dirichlet (Dirichlet = 0, prescribes the value) or Neumann (Neumann = 1, prescribes the flux). Fluxes entering a boundary cell are positive by convention. For non-negative concentrations (posconc = 1), negative Neumann boundary fluxes are linearised following Patankar (1980) to preserve positivity.

The Thomas algorithm (tridiagonal()) solves the resulting banded system.

Original FORTRAN author: Lars Umlauf.

pygotm.util.diff_center.diff_center(nlev, dt, cnpar, posconc, h, bc_up, bc_down, y_up, y_down, nu_y, l_sour, q_sour, tau_r, y_obs, y, au, bu, cu, du, ru, qu)[source]

Solve diffusion for a cell-centred variable using Crank–Nicolson.

Advances y at cell centres (indices 1 to nlev) one time step. Solves the one-dimensional diffusion equation including linear and explicit source terms, optional relaxation toward observed values, and Dirichlet or Neumann boundary conditions.

Return type:

None

Parameters:
pygotm.util.diff_center.diff_center_batch(batch_size, nlev, dt, cnpar, posconc, h, bc_up, bc_down, y_up, y_down, nu_y, l_sour, q_sour, tau_r, y_obs, y, au, bu, cu, du, ru, qu)[source]

Batch variant: process batch_size columns in parallel with numba.prange.

Return type:

None

Parameters:

Vertical Diffusion (Cell Faces)

Translation of diff_face.F90. Solves the same diffusion equation as diff_center but for variables defined at grid faces (interfaces) rather than cell centres. The solved range spans interior face indices 1 to nlev − 1.

A bug fix attributed to Georg Umgiesser handles the nlev == 2 edge case: when only two layers are present, boundary diffusivities and values are replicated from the single interior face to stabilise the tridiagonal system.

Vertical diffusion for face-centred variables — translation of diff_face.F90.

Solves the one-dimensional diffusion equation for variables defined at grid faces (interfaces) rather than cell centres. Uses the same Crank–Nicolson approach as diff_center but with the stencil appropriate for face-centred quantities. Source terms l_sour (linear, implicit) and q_sour (constant, explicit) are supported; relaxation is not included.

Boundary conditions (bc_up, bc_down) are Dirichlet (Dirichlet = 0) or Neumann (Neumann = 1). The solved range spans indices 1 to nlev - 1 (interior faces only).

Includes a bug-fix for nlev == 2 attributed to Georg Umgiesser: when only two layers are present, boundary diffusivities and values are replicated from the single interior face to stabilise the system.

The Thomas algorithm (tridiagonal()) solves the resulting banded system.

Original FORTRAN author: Lars Umlauf.

pygotm.util.diff_face.diff_face(nlev, dt, cnpar, h, bc_up, bc_down, y_up, y_down, nu_y, l_sour, q_sour, y, au, bu, cu, du, ru, qu)[source]

Solve diffusion for a face-centred variable using Crank–Nicolson.

Advances y at interior faces (indices 1 to nlev−1) one time step. Implicitness is controlled by cnpar; optional linear (l_sour) and explicit constant (q_sour) source terms are included. Boundary conditions bc_up and bc_down are DIRICHLET (prescribe value) or NEUMANN (prescribe flux).

Return type:

None

Parameters:
pygotm.util.diff_face.diff_face_batch(batch_size, nlev, dt, cnpar, h, bc_up, bc_down, y_up, y_down, nu_y, l_sour, q_sour, y, au, bu, cu, du, ru, qu)[source]

Batch variant: process batch_size columns in parallel with numba.prange.

Return type:

None

Parameters:

Equation of State

Translation of density.F90. Computes in-situ density \(\bar{\rho}(S,\Theta,P)\), potential density, and buoyancy expansion coefficients from salinity, temperature, and pressure.

Three methods are selected via DensityState.density_method:

Code

Constant

Description

1

METHOD_TEOS10

Full TEOS-10 EOS via the gsw package: gsw_rho, gsw_sigma0, gsw_alpha, gsw_beta.

2

METHOD_LINEAR_TEOS10

Linearised EOS; \(\alpha_0\) and \(\beta_0\) computed from TEOS-10 at the user reference point \((S_0, T_0, p_0)\).

3

METHOD_LINEAR_USER

Linearised EOS with user-supplied \(\rho_0\), \(\alpha_0\), \(\beta_0\).

The TEOS-10 standard specific heat of seawater at \(S_A = 0\), \(C_T = 0\), \(p = 0\) is CP0 = 3991.86795711963 J/(kg·K), matching gsw_mod_teos10_constants::gsw_cp0 in the Fortran source.

The linear EOS (methods 2 and 3) computes density as:

\[\rho = \rho_b\,(1 - \alpha_0\,(T - T_0) + \beta_0\,(S - S_0))\]

No pressure dependence is applied in the linear methods (matching the Fortran comment from Lars Umlauf).

DensityState holds all module-level state. Call init_density() once to allocate arrays, then do_density() each time step.

Equation of state (density) — translation of density.F90.

Computes density \(\bar{\rho}(S, \Theta, P)\), potential density, and buoyancy expansion coefficients from salinity, temperature, and pressure. Three methods are supported, selected via DensityState.density_method:

  • METHOD_TEOS10 (1) — Full TEOS-10 equation of state via the gsw package: gsw_rho, gsw_sigma0, gsw_alpha, gsw_beta.

  • METHOD_LINEAR_TEOS10 (2) — Linearised EOS; thermal expansion \(\alpha_0\) and haline contraction \(\beta_0\) are computed from TEOS-10 at the user reference point \((S_0, T_0, p_0)\).

  • METHOD_LINEAR_USER (3) — Linearised EOS with user-supplied \(\rho_0\), \(\alpha_0\), \(\beta_0\).

The GOTM Fortran source (density.F90) uses the TEOS-10 Gibbs SeaWater toolbox (gsw) rather than the original UNESCO 1983 polynomial described in the historical module header. This Python translation wraps the same TEOS-10 library via the gsw Python package (McDougall & Barker 2011).

TEOS-10 reference: IOC, SCOR and IAPSO (2010). The International Thermodynamic Equation of Seawater 2010. Intergovernmental Oceanographic Commission, Manuals and Guides No. 56.

Public interface: init_density(), do_density(), get_rho(), get_alpha(), get_beta(), clean_density(), DensityState.

Original FORTRAN authors: Hans Burchard, Karsten Bolding.

class pygotm.util.density.DensityState[source]

Bases: object

State for the GOTM density / EOS module.

Mirrors the module-level public and private variables in density.F90. Set configuration attributes, call init_density(), then call do_density() each time step.

Attributes

density_methodint

EOS method: METHOD_TEOS10 (1), METHOD_LINEAR_TEOS10 (2), or METHOD_LINEAR_USER (3).

T0, S0, p0float

Reference temperature [°C], salinity [g/kg], pressure [dbar]. Used by methods 2 and 3 for the linearisation reference point.

rho0float

Reference density [kg/m³], default 1027. Used by method 3.

alpha0, beta0float

Thermal expansion coefficient [1/K] and haline contraction coefficient [kg/g]. For method 2, computed from TEOS-10 at (S0, T0, p0). For method 3, supplied by the caller.

cpfloat

Specific heat capacity of seawater [J/(kg·K)]; set by init_density.

alpha, betanp.ndarray or None, shape (nlev+1,)

Per-interface thermal expansion and haline contraction.

rho, rho_pnp.ndarray or None, shape (nlev+1,)

In-situ and potential density at cell centres. Index 0 is seabed (GOTM convention); index nlev is surface.

pygotm.util.density.clean_density(state)[source]

Release the density arrays held by state.

Sets alpha, beta, rho, and rho_p to None, mirroring Fortran deallocate. Call init_density() again before the next do_density() call.

Return type:

None

Parameters:

state (DensityState)

pygotm.util.density.do_density(state, nlev, S, T, p, pi)[source]

Compute density, potential density, and expansion coefficients.

Interface salinity and temperature are computed as arithmetic means of adjacent cell-centre values. Boundary faces (seabed index 0 and surface index nlev) use the single adjacent cell value.

Method-specific computation:

  • METHOD_TEOS10 (1): computes rho[1:nlev+1] via gsw_rho, rho_p[1:nlev+1] via gsw_sigma0 + 1000, and the full interface arrays alpha and beta via gsw_alpha / gsw_beta evaluated at the interface salinity, temperature, and pressure.

  • METHOD_LINEAR_TEOS10 / METHOD_LINEAR_USER (2, 3): applies the linear EOS \(\rho_p = \rho_b(1 - \alpha_0(T-T_0) + \beta_0(S-S_0))\) at all cell centres. In-situ density is set equal to potential density (no pressure dependence implemented, following the Fortran source).

Parameters

stateDensityState

Initialised density state.

nlevint

Number of vertical layers.

Snp.ndarray, shape (nlev+1,)

Absolute salinity at cell centres [g/kg].

Tnp.ndarray, shape (nlev+1,)

Conservative temperature at cell centres [°C].

pnp.ndarray, shape (nlev+1,)

In-situ pressure at cell centres [dbar].

pinp.ndarray, shape (nlev+1,)

In-situ pressure at cell interfaces [dbar].

rtype:

None

Parameters:
Return type:

None

pygotm.util.density.get_alpha(state, S, T, p)[source]

Compute thermal expansion coefficient for a single water parcel.

  • METHOD_TEOS10 (1): returns gsw_alpha(S, T, p) from TEOS-10.

  • METHOD_LINEAR_TEOS10 / METHOD_LINEAR_USER (2, 3): returns the constant reference value alpha0 set during initialisation.

Parameters

stateDensityState

Initialised density state.

Sfloat

Absolute salinity [g/kg].

Tfloat

Conservative temperature [°C].

pfloat

In-situ pressure [dbar].

rtype:

float

Parameters:
Return type:

float

pygotm.util.density.get_beta(state, S, T, p)[source]

Compute haline contraction coefficient for a single water parcel.

  • METHOD_TEOS10 (1): returns gsw_beta(S, T, p) from TEOS-10.

  • METHOD_LINEAR_TEOS10 / METHOD_LINEAR_USER (2, 3): returns the constant reference value beta0 set during initialisation.

Parameters

stateDensityState

Initialised density state.

Sfloat

Absolute salinity [g/kg].

Tfloat

Conservative temperature [°C].

pfloat

In-situ pressure [dbar].

rtype:

float

Parameters:
Return type:

float

pygotm.util.density.get_rho(state, S, T, p=None)[source]

Compute density for a single water parcel.

  • METHOD_TEOS10 (1): returns gsw_rho(S, T, p) when p is given, or gsw_sigma0(S, T) + 1000 (potential density) when p is omitted.

  • METHOD_LINEAR_TEOS10 / METHOD_LINEAR_USER (2, 3): returns \(\rho_b(1 - \alpha_0(T-T_0) + \beta_0(S-S_0))\). No pressure dependence is applied (following the Fortran source).

Parameters

stateDensityState

Initialised density state.

Sfloat

Absolute salinity [g/kg].

Tfloat

Conservative temperature [°C].

pfloat, optional

In-situ pressure [dbar]. If omitted under method 1, returns potential density (referenced to 0 dbar).

rtype:

float

Parameters:
Return type:

float

pygotm.util.density.init_density(state, nlev)[source]

Initialise the density module state for a column of nlev layers.

Method-specific initialisation:

  • METHOD_TEOS10 (1): sets cp = CP0; all other coefficients are computed per call in do_density().

  • METHOD_LINEAR_TEOS10 (2): computes _rhob = gsw_sigma0(S0,T0) + 1000, alpha0 = gsw_alpha(S0,T0,p0), beta0 = gsw_beta(S0,T0,p0), and cp = CP0 from TEOS-10 at the user reference point.

  • METHOD_LINEAR_USER (3): sets _rhob = rho0 from the caller-supplied reference density; alpha0 and beta0 must be set by the caller before calling this function.

After all methods, allocates alpha, beta (shape nlev+1, filled with alpha0 / beta0), rho, and rho_p (shape nlev+1, zero-initialised).

Parameters

stateDensityState

Pre-configured state (set density_method and reference values first).

nlevint

Number of vertical layers.

rtype:

None

Parameters:
Return type:

None

TEOS-10 GSW Library

A Numba-compiled subset of the TEOS-10 Gibbs SeaWater toolbox, mirroring the Fortran files in extern/gsw/. The equation-of-state functions used from the time loop are scalar @numba.njit functions; the salinity conversion helpers are vectorised Python wrappers around GOTM’s bundled SAAR data.

Exported functions

Function

Description

gsw_rho

In-situ density [kg/m³] from \((S_A, C_T, p)\)

gsw_sigma0

Potential density anomaly [kg/m³] at \(p = 0\)

gsw_alpha

Thermal expansion coefficient \(\alpha\) [1/K]

gsw_beta

Haline contraction coefficient \(\beta\) [kg/g]

gsw_specvol

Specific volume [m³/kg]

gsw_sa_from_sp

Absolute salinity from practical salinity

gsw_sp_from_sa

Practical salinity from absolute salinity

gsw_saar

Salinity anomaly ratio (SAAR) from the GOTM lookup table

Numba-compiled TEOS-10 GSW library.

Mirrors the Fortran folder layout:

extern/gsw/
modules/

gsw_mod_teos10_constants.f90 → modules/gsw_mod_teos10_constants.py gsw_mod_specvol_coefficients.f90 → modules/gsw_mod_specvol_coefficients.py

toolbox/

gsw_saar.f90 → toolbox/gsw_saar.py gsw_sa_from_sp.f90 → toolbox/gsw_sa_from_sp.py gsw_sp_from_sa.f90 → toolbox/gsw_sp_from_sa.py gsw_specvol.f90 → toolbox/gsw_specvol.py gsw_sigma0.f90 → toolbox/gsw_sigma0.py gsw_alpha.f90 → toolbox/gsw_alpha.py gsw_beta.f90 → toolbox/gsw_beta.py gsw_rho.f90 → toolbox/gsw_rho.py

The equation-of-state functions used from the compiled time loop are scalar @numba.njit functions. The salinity conversion helpers are vectorized Python wrappers around GOTM’s bundled GSW SAAR data and are used during setup/output.

pygotm.util.gsw.gsw_alpha(sa, ct, p)[source]

Thermal expansion coefficient [1/K] w.r.t. Conservative Temperature.

Translated from gsw_alpha.f90 (Roquet et al. 2015).

Return type:

float

Parameters:

Parameters

safloat

Absolute Salinity [g/kg]

ctfloat

Conservative Temperature [°C]

pfloat

Sea pressure [dbar]

Returns

float

alpha [1/K]

pygotm.util.gsw.gsw_beta(sa, ct, p)[source]

Haline contraction coefficient [kg/g] at constant Conservative Temperature.

Translated from gsw_beta.f90 (Roquet et al. 2015).

Return type:

float

Parameters:

Parameters

safloat

Absolute Salinity [g/kg]

ctfloat

Conservative Temperature [°C]

pfloat

Sea pressure [dbar]

Returns

float

beta [kg/g]

pygotm.util.gsw.gsw_rho(sa, ct, p)[source]

In-situ density from Absolute Salinity, Conservative Temperature, pressure.

Translated from gsw_rho.f90.

Return type:

float

Parameters:

Parameters

safloat

Absolute Salinity [g/kg]

ctfloat

Conservative Temperature [°C]

pfloat

Sea pressure [dbar]

Returns

float

In-situ density [kg/m³]

pygotm.util.gsw.gsw_sa_from_sp(sp, p, long, lat)[source]

Calculate Absolute Salinity from Practical Salinity.

Return type:

ndarray | float

Parameters:
pygotm.util.gsw.gsw_saar(p, long, lat)[source]

Calculate the Absolute Salinity Anomaly Ratio.

Return type:

ndarray | float

Parameters:
pygotm.util.gsw.gsw_sigma0(sa, ct)[source]

Potential density anomaly [kg/m³] with reference pressure 0 dbar.

Translated from gsw_sigma0.f90. Uses only the p=0 terms of the specific-volume polynomial (Roquet et al. 2015).

Return type:

float

Parameters:

Parameters

safloat

Absolute Salinity [g/kg]

ctfloat

Conservative Temperature [°C]

Returns

float

sigma0 = potential_density - 1000 [kg/m³]

pygotm.util.gsw.gsw_sp_from_sa(sa, p, long, lat)[source]

Calculate Practical Salinity from Absolute Salinity.

Return type:

ndarray | float

Parameters:
pygotm.util.gsw.gsw_specvol(sa, ct, p)[source]

Specific volume [m³/kg] from Absolute Salinity, Conservative Temperature, pressure.

Translated from gsw_specvol.f90 (Roquet et al. 2015, 75-term polynomial).

Return type:

float

Parameters:

Parameters

safloat

Absolute Salinity [g/kg]

ctfloat

Conservative Temperature [°C]

pfloat

Sea pressure [dbar] (absolute pressure minus 10.1325 dbar)

Returns

float

Specific volume [m³/kg]

Time Management

Translation of time.F90. Time is represented as a pair of integers — true Julian day and integer seconds since midnight — so all arithmetic is simple integer operations.

timefmt modes (configured on GotmTime):

Mode

Behaviour

1

MaxN given directly; fake start date 2000-01-01 00:00:00 is used.

2

start and stop strings given; MaxN computed from total duration divided by timestep.

3

start string and MaxN given; stop is computed.

Time strings use the format 'yyyy-mm-dd hh:mm:ss'.

GotmTime is a Python dataclass that mirrors the Fortran module-level variables. Call init_time() once, then update_time() each step.

Time management — translation of time.F90.

Time is represented as a pair of integers: true Julian day and integer seconds since midnight. This makes all time arithmetic simple integer operations.

Three timefmt modes control how start/stop/MaxN are resolved:

  • 1 — MaxN given directly; a fake start date of 2000-01-01 is used.

  • 2 — start and stop strings given; MaxN is computed from the total duration divided by timestep.

  • 3 — start string and MaxN given; stop is computed.

Original FORTRAN authors: Karsten Bolding, Hans Burchard.

Public interface: GotmTime, calendar_date(), julian_day(), time_diff(), sunrise_sunset(), in_time_interval(), read_time_string(), write_time_string().

class pygotm.util.time.GotmTime(timestep=3600.0, timefmt=2, start='2000-01-01 00:00:00', stop='2001-01-01 00:00:00')[source]

Bases: object

GOTM time-management state, mirroring the Fortran ‘time’ module.

Time is represented as true Julian day + integer seconds since midnight. All public fields mirror the Fortran module-level variables.

timefmt controls how start/stop/MaxN are resolved:

1 — MaxN given directly; fake start date 2000-01-01 00:00:00 2 — start and stop strings given; MaxN computed from duration 3 — start string + MaxN given; stop computed

Parameters:
timestep: float = 3600.0
timefmt: int = 2
start: str = '2000-01-01 00:00:00'
stop: str = '2001-01-01 00:00:00'
timestr: str = ''
julianday: int = 0
secondsofday: int = 0
fsecs: float = 0.0
fsecondsofday: float = 0.0
simtime: float = 0.0
yearday: int = 0
MinN: int = 1
MaxN: int = 0
jul1: int = 0
secs1: int = 0
jul2: int = 0
secs2: int = 0
init_time()[source]

Initialise the time system.

Resolves start/stop/MaxN according to timefmt, then calls update_time(0) to set all state for the initial instant. Mirrors init_time() from time.F90.

Return type:

None

update_time(n)[source]

Update all time-tracking variables for time step n.

Parameters:

n (int) – Time step index (0 = start of simulation).

Return type:

None

Mirrors update_time(n) from time.F90.

pygotm.util.time.calendar_date(julian)[source]

Convert true Julian day to calendar date — year, month and day. Based on a similar routine in Numerical Recipes.

Parameters:

julian (int) – True Julian day number.

Return type:

tuple[int, int, int]

Returns:

(yyyy, mm, dd) calendar date.

pygotm.util.time.julian_day(yyyy, mm, dd)[source]

Convert a calendar date to true Julian day. Based on a similar routine in Numerical Recipes.

Parameters:

  • yyyy (int) – Year (use negative integers for BCE; 0 is skipped — 1 BCE is -1).

  • mm (int) – Month (1–12).

  • dd (int) – Day (1–31).

Return type:

int

Returns:

True Julian day number.

pygotm.util.time.time_diff(jul1, secs1, jul2, secs2)[source]

Return the time difference between two dates in seconds.

The dates are given as Julian day and seconds of that day. Returns (jul1, secs1) − (jul2, secs2) in seconds.

Return type:

float

Parameters:
pygotm.util.time.sunrise_sunset(latitude, declination)[source]

Return the times of sunrise and sunset.

Parameters:

  • latitude (float) – Geographic latitude in degrees.

  • declination (float) – Solar declination in degrees.

Return type:

tuple[float, float]

Returns:

(sunrise, sunset) in decimal hours (UTC).

pygotm.util.time.in_time_interval(j1, s1, j, s, j2, s2)[source]

Return True if (j, s) lies within the closed interval [(j1, s1), (j2, s2)].

Times are expressed as (Julian day, seconds of that day).

Return type:

bool

Parameters:
pygotm.util.time.read_time_string(timestr)[source]

Convert a time string to Julian day and seconds of that day.

The format of the time string must be: ‘yyyy-mm-dd hh:mm:ss’.

Return type:

tuple[int, int]

Returns:

(julian_day, seconds_of_day)

Parameters:

timestr (str)

pygotm.util.time.write_time_string(jul, secs)[source]

Format Julian day and seconds of that day to a time string.

Output format: ‘yyyy-mm-dd hh:mm:ss’.

Return type:

str

Parameters:

Flux Conversion

Translation of convert_fluxes.F90. Converts surface heat, salinity, and shortwave radiative fluxes to the temperature and buoyancy flux forms required by the KPP turbulence closure.

convert_fluxes() returns a 6-tuple:

t_flux, s_flux, bt_flux, bs_flux, t_rad, b_rad = convert_fluxes(
    state, nlev, gravity, shf, ssf, rad, T_srf, S_srf
)
Return values

Name

Units

Description

t_flux

K·m/s

Temperature flux from surface heat flux

s_flux

psu·m/s

Salinity flux from P − E

bt_flux

m²/s³

Buoyancy flux from surface heat flux

bs_flux

m²/s³

Buoyancy flux from surface salinity flux

t_rad

K·m/s, shape (nlev+1,)

Temperature flux profile from shortwave radiation

b_rad

m²/s³, shape (nlev+1,)

Buoyancy flux profile from shortwave radiation

Only called when the KPP turbulence model is active.

Flux conversion for KPP turbulence model — translation of convert_fluxes.F90.

Converts surface heat, salinity, and shortwave radiative fluxes into the temperature and buoyancy flux forms required by the KPP turbulence closure. Three surface contributions are handled:

  1. Surface heat flux → temperature flux and thermal buoyancy flux.

  2. Surface salinity flux (P − E) → salinity flux and haline buoyancy flux.

  3. Shortwave radiative flux profile → radiative temperature and buoyancy flux profiles.

Only called when the KPP turbulence model is active. Callers outside GOTM must supply fluxes consistent with their own equation of state.

Original FORTRAN author: Lars Umlauf.

Public interface: convert_fluxes().

pygotm.util.convert_fluxes.convert_fluxes(state, nlev, gravity, shf, ssf, rad, T_srf, S_srf)[source]

Convert surface and radiative fluxes to temperature and buoyancy fluxes.

The algorithm proceeds in three steps:

  1. Evaluate the thermal expansion coefficient \(\alpha_0\) and haline contraction coefficient \(\beta_0\) at the surface using the configured equation of state.

  2. Convert the surface heat flux and salinity flux to temperature and buoyancy fluxes:

    \[t_{\mathrm{flux}} = -\frac{q_{\mathrm{shf}}}{\rho_0 c_p},\quad bt_{\mathrm{flux}} = g\,\alpha_0\,t_{\mathrm{flux}}\]
    \[s_{\mathrm{flux}} = -q_{\mathrm{ssf}},\quad bs_{\mathrm{flux}} = -g\,\beta_0\,s_{\mathrm{flux}}\]

  3. Convert the shortwave radiative profile to temperature and buoyancy flux profiles:

    \[t_{\mathrm{rad}} = \frac{\mathrm{rad}}{\rho_0 c_p},\quad b_{\mathrm{rad}} = g\,\alpha\,t_{\mathrm{rad}}\]

Parameters

stateDensityState

Initialised density state providing rho0, cp, alpha, beta.

nlevint

Number of vertical layers.

gravityfloat

Gravitational acceleration [m/s²].

shffloat

Surface heat flux [W/m²]. Positive = into ocean.

ssffloat

Surface salinity flux [psu·m/s or kg/(m²·s)]. Sign: positive = freshwater input → salinity decrease.

radnp.ndarray, shape (nlev+1,)

Shortwave radiative heat flux profile [W/m²].

T_srffloat

Surface temperature [°C] — used to evaluate α, β at the surface.

S_srffloat

Surface salinity [g/kg] — used to evaluate α, β at the surface.

Returns

t_fluxfloat

Temperature flux from surface heat flux [K·m/s].

s_fluxfloat

Salinity flux from P-E [psu·m/s].

bt_fluxfloat

Buoyancy flux from surface heat flux [m²/s³].

bs_fluxfloat

Buoyancy flux from surface salinity flux [m²/s³].

t_radnp.ndarray, shape (nlev+1,)

Temperature flux profile from shortwave radiation [K·m/s].

b_radnp.ndarray, shape (nlev+1,)

Buoyancy flux profile from shortwave radiation [m²/s³].

rtype:

tuple[float, float, float, float, ndarray, ndarray]

Parameters:
Return type:

tuple[float, float, float, float, ndarray, ndarray]

Grid Interpolation

Translation of gridinterpol.F90. gridinterpol() linearly interpolates (and extrapolates) observational data defined on an arbitrary structured depth grid to the (potentially moving) GOTM model grid:

  • Model levels above the topmost observation receive the topmost observed value.

  • Model levels below the deepest observation receive the deepest observed value.

  • Interior levels are interpolated linearly between the two nearest observations.

Original FORTRAN authors: Karsten Bolding, Hans Burchard.

Grid interpolation from observation space to model grid — translation of gridinterpol.F90.

Linearly interpolates (and extrapolates) observational data defined on an arbitrary structured depth grid to the (potentially moving) GOTM model grid. Extrapolation clamps values: model levels above the topmost observation receive the topmost observed value; model levels below the deepest observation receive the deepest observed value.

Original FORTRAN authors: Karsten Bolding, Hans Burchard.

Public interface: gridinterpol().

pygotm.util.gridinterpol.gridinterpol(obs_z, obs_prof, model_z, nlev)[source]

Linearly interpolate/extrapolate observational data to the model grid.

Observational data on an arbitrary structured grid (obs_z, obs_prof) are mapped to the (moving) model grid (model_z). Linear extrapolation is used outside the observational range: surface levels above obs_z[N] receive the topmost observed value, and bottom levels below obs_z[1] receive the lowest observed value.

Return type:

ndarray

Parameters:

Parameters

obs_znp.ndarray, shape (N+1,)

Observation depth levels [m]. Indices 1..N are the actual levels; index 0 is a padding element (never used by this routine, following GOTM Fortran convention with DIMENSION(0:N)).

obs_profnp.ndarray, shape (N+1, cols)

Observed profile values at each obs level and column.

model_znp.ndarray, shape (nlev+1,)

Model depth levels [m], indexed 0..nlev.

nlevint

Number of model layers.

Returns

model_profnp.ndarray, shape (nlev+1, cols)

Interpolated profile on the model grid. Index 0 is not set by this routine (consistent with the Fortran loop bounds 1..nlev).

Lagrangian Particle Tracking

Translation of lagrange.F90. Implements the Visser (1997) random-walk scheme for Lagrangian particles in spatially inhomogeneous turbulence. Each particle is advanced by:

\[z^{n+1} = z^n + \partial_z \nu_t(z^n)\,\Delta t + R\left[2\,r^{-1}\,\nu_t\!\left(z^n + \tfrac{1}{2}\partial_z\nu_t(z^n)\,\Delta t\right) \Delta t\right]^{1/2}\]

where \(R\) is a zero-mean random variable with variance \(\langle R^2 \rangle = r\).

Fixed parameters:

  • VISC_BACK = 0.0e-6 m²/s — background viscosity floor.

  • RND_VAR = 0.333333333 — variance \(r\) of the random walk.

Reflective boundary conditions are applied at the surface (\(z = 0\)) and the bottom (\(z = -\mathrm{depth}\)). Pass a seeded np.random.Generator for reproducible results.

Original FORTRAN authors: Hans Burchard, Karsten Bolding.

Lagrangian particle random walk — translation of lagrange.F90.

Implements the Visser (1997) random-walk scheme for spatially inhomogeneous turbulence. Each active particle is advanced from \(z^n\) to \(z^{n+1}\) by:

\[z^{n+1} = z^n + \partial_z \nu_t(z^n)\,\Delta t + R\left[2\,r^{-1}\,\nu_t\!\left(z^n + \tfrac{1}{2}\partial_z\nu_t(z^n)\,\Delta t\right) \Delta t\right]^{1/2}\]

where \(R\) is a zero-mean random variable with variance \(\langle R^2 \rangle = r\).

Fixed parameters: background viscosity VISC_BACK = 0.0e-6 m²/s, random-walk variance RND_VAR = 1/3. Semi-implicit viscosity correction (visc_corr) is disabled, matching the Fortran default. Reflective boundary conditions are applied at the surface and bottom.

Original FORTRAN authors: Hans Burchard, Karsten Bolding.

Public interface: lagrange(), VISC_BACK, RND_VAR.

pygotm.util.lagrange.lagrange(nlev, dt, zlev, nuh, w, npar, active, zi, zp, rng=None)[source]

Lagrangian particle random walk for spatially inhomogeneous turbulence.

Implements the Visser (1997) random-walk scheme. Each particle position zp[n] and its enclosing layer index zi[n] are updated in-place.

Return type:

None

Parameters:
The step formula is:
z^{n+1} = z^n + (dzn[i] + w) * dt

  • sqrt(2 * rnd_var_inv * dt_inv * visc) * rnd * dt

Reflective boundary conditions are applied at the surface (z=0) and bottom (z=-depth). Particle indices zi are 1-based (matching Fortran DIMENSION(0:nlev) convention).

Parameters

nlevint

Number of model layers.

dtfloat

Time step [s].

zlevnp.ndarray, shape (nlev+1,)

Layer interface depths [m], indexed 0..nlev. zlev[0] is the bottom.

nuhnp.ndarray, shape (nlev+1,)

Eddy diffusivity [m²/s] at layer interfaces, indexed 0..nlev.

wfloat

Vertical velocity [m/s], positive upward.

nparint

Number of particles.

activenp.ndarray of bool, shape (npar,)

Active flag per particle (passed for interface compatibility; not used inside the walk loop, matching GOTM Fortran behaviour).

zinp.ndarray of int, shape (npar,)

Layer index (1-based) enclosing each particle. Modified in-place.

zpnp.ndarray of float, shape (npar,)

Particle vertical position [m]. Modified in-place.

rngnp.random.Generator, optional

Random number generator. If None, a default generator is used. Pass a seeded generator for reproducible results.

ODE Solvers

Translation of ode_solvers.F90. Provides 11 numerical solvers for the biogeochemical reaction-term ODE:

\[\partial_t c_i = P_i(\mathbf{c}) - D_i(\mathbf{c}), \quad i = 1, \ldots, I\]

where \(c_i\) are species concentrations and \(P_i\), \(D_i\) are the production and destruction terms (Burchard et al. 2003).

Available solvers

Code

Method

Conservative

Positive

1

First-order explicit (Euler)

No

No

2

Second-order explicit Runge-Kutta

No

No

3

Fourth-order explicit Runge-Kutta

No

No

4

First-order Patankar

No

Yes

5

Second-order Patankar-Runge-Kutta

No

Yes

6

Fourth-order Patankar-Runge-Kutta (non-functional)

—

—

7

First-order Modified Patankar

Yes

Yes

8

Second-order Modified Patankar-Runge-Kutta

Yes

Yes

9

Fourth-order Modified Patankar-Runge-Kutta (non-functional)

—

—

10

First-order Extended Modified Patankar (EMP)

Stoichiometric

Yes

11

Second-order EMP-Runge-Kutta

Stoichiometric

Yes

Solvers 6 and 9 are not yet developed (matching the Fortran status). EMP schemes 10–11 were developed by Bruggeman et al. (2005) to extend Modified Patankar to full stoichiometric conservation with multiple limiting nutrients.

ODE solvers for biogeochemical models — translation of ode_solvers.F90.

Provides 11 numerical solvers for the reaction-term ODEs that arise in biogeochemical models:

\[\partial_t c_i = P_i(\mathbf{c}) - D_i(\mathbf{c}), \quad i = 1, \ldots, I\]

where \(c_i\) are species concentrations and \(P_i\), \(D_i\) are the production and destruction terms (Burchard et al. 2003).

Available solvers

Code

Method

Conservative

Positive

1

First-order explicit (Euler)

No

No

2

Second-order explicit Runge-Kutta

No

No

3

Fourth-order explicit Runge-Kutta

No

No

4

First-order Patankar

No

Yes

5

Second-order Patankar-Runge-Kutta

No

Yes

6

Fourth-order Patankar-Runge-Kutta (non-functional)

—

—

7

First-order Modified Patankar

Yes

Yes

8

Second-order Modified Patankar-Runge-Kutta

Yes

Yes

9

Fourth-order Modified Patankar-Runge-Kutta (non-functional)

—

—

10

First-order Extended Modified Patankar (EMP)

Stoichiometric

Yes

11

Second-order Extended Modified Patankar-Runge-Kutta (EMP)

Stoichiometric

Yes

Schemes 6 and 9 are not yet developed in Fortran GOTM and are non-functional here. EMP schemes (10–11) — developed by Bruggeman et al. (2005) — extend Modified Patankar to full stoichiometric conservation with multiple limiting nutrients.

Original FORTRAN authors: Hans Burchard, Karsten Bolding.

Public interface: ode_solver(), RhsCallback, PpddCallback.

pygotm.util.ode_solvers.ode_solver(solver, numc, nlev, dt, cc, get_rhs=None, get_ppdd=None)[source]

Dispatch to one of 11 ODE solvers for biogeochemical reaction equations.

Solvers 1–3 and 10–11 require get_rhs; solvers 4–9 require get_ppdd. See the module docstring for the full solver table.

Return type:

None

Parameters:

Parameters

solverint

Solver identifier (1–11).

numcint

Number of biogeochemical state variables.

nlevint

Number of vertical levels (cc has shape (numc, nlev+1)).

dtfloat

Time step [s].

ccnp.ndarray, shape (numc, nlev+1)

State variable array, modified in-place. Level 0 is not updated.

get_rhscallable, optional

RHS callback get_rhs(first, numc, nlev, cc) → rhs. Required for solvers 1, 2, 3, 10, 11.

get_ppddcallable, optional

Production-destruction callback get_ppdd(first, numc, nlev, cc) → (pp, dd). Required for solvers 4–9.

pygotm.util.ode_solvers.euler_forward(dt, numc, nlev, cc, get_rhs)[source]

First-order Euler-forward (E1) scheme.

One evaluation of the right-hand side per time step:

\[c_i^{n+1} = c_i^n + \Delta t\,\bigl(P_i(c^n) - D_i(c^n)\bigr)\]

Updates cc in-place. Level 0 (index 0) is never modified.

Return type:

None

Parameters:
pygotm.util.ode_solvers.runge_kutta_2(dt, numc, nlev, cc, get_rhs)[source]

Second-order Runge-Kutta — Heun’s method (explicit trapezoidal).

Two evaluations of the right-hand side per time step:

\[\begin{split}c_i^{(1)} &= c_i^n + \Delta t\,f_i(c^n) \\ c_i^{n+1} &= c_i^n + \tfrac{\Delta t}{2}\bigl(f_i(c^n) + f_i(c^{(1)})\bigr)\end{split}\]

where \(f_i = P_i - D_i\). Updates cc in-place. Level 0 is never modified.

Return type:

None

Parameters:
pygotm.util.ode_solvers.runge_kutta_4(dt, numc, nlev, cc, get_rhs)[source]

Fourth-order Runge-Kutta scheme (ode_solvers.F90 variant — full-dt steps).

Four evaluations of the right-hand side per time step using full-\(\Delta t\) intermediates (differs from the standard half-step RK4 in ode_solvers_template):

\[\begin{split}c^{(1)} &= c^n + \Delta t\,f(c^n) \\ c^{(2)} &= c^n + \Delta t\,f(c^{(1)}) \\ c^{(3)} &= c^n + \Delta t\,f(c^{(2)}) \\ c^{n+1} &= c^n + \tfrac{\Delta t}{3}\bigl( \tfrac{1}{2}f(c^n) + f(c^{(1)}) + f(c^{(2)}) + \tfrac{1}{2}f(c^{(3)})\bigr)\end{split}\]

Updates cc in-place. Level 0 is never modified.

Return type:

None

Parameters:
pygotm.util.ode_solvers.patankar(dt, numc, nlev, cc, get_ppdd)[source]

First-order Patankar-Euler (PE1) scheme.

One evaluation of the production/destruction terms per time step. Not conservative; unconditionally positive:

\[c_i^{n+1} = \frac{c_i^n + \Delta t\,P_i(c^n)}{1 + \Delta t\,D_i(c^n)/c_i^n}\]

Updates cc in-place. Level 0 is never modified.

Return type:

None

Parameters:
pygotm.util.ode_solvers.patankar_runge_kutta_2(dt, numc, nlev, cc, get_ppdd)[source]

Second-order Patankar-Runge-Kutta (PRK2) scheme.

Two evaluations of the production/destruction terms per time step. Not conservative; unconditionally positive (Burchard et al. 2003).

Stage 1 — Patankar-Euler predictor:

\[c_i^{(1)} = \frac{c_i^n + \Delta t\,P_i(c^n)}{1 + \Delta t\,D_i(c^n)/c_i^n}\]

Stage 2 — second-order corrector:

\[c_i^{n+1} = \frac{c_i^n + \tfrac{\Delta t}{2}(P_i(c^n)+P_i(c^{(1)}))} {1 + \tfrac{\Delta t}{2}(D_i(c^n)+D_i(c^{(1)}))/c_i^{(1)}}\]

Updates cc in-place. Level 0 is never modified.

Return type:

None

Parameters:
pygotm.util.ode_solvers.patankar_runge_kutta_4(dt, numc, nlev, cc, get_ppdd)[source]

Fourth-order Patankar-Runge-Kutta (PRK4) scheme — does not work.

This scheme has not yet been developed in GOTM Fortran or here; the implementation is a placeholder. Updates cc in-place. Level 0 is never modified.

Return type:

None

Parameters:
pygotm.util.ode_solvers.modified_patankar(dt, numc, nlev, cc, get_ppdd)[source]

First-order Modified Patankar-Euler (MPE1) scheme.

One evaluation of the production/destruction terms per time step. Conservative and unconditionally positive (Burchard et al. 2003):

\[c_i^{n+1} = c_i^n + \Delta t\left[ \sum_{j \neq i} p_{ij}(c^n)\,\frac{c_j^{n+1}}{c_j^n} + p_{ii}(c^n) - \sum_j d_{ij}(c^n)\,\frac{c_i^{n+1}}{c_i^n} \right]\]

Solved as a linear system per vertical level using matrix_solve(). Updates cc in-place. Level 0 is never modified.

Return type:

None

Parameters:
pygotm.util.ode_solvers.modified_patankar_2(dt, numc, nlev, cc, get_ppdd)[source]

Second-order Modified Patankar-Runge-Kutta (MPRK2) scheme.

Two evaluations of the production/destruction terms per time step. Conservative and unconditionally positive (Burchard et al. 2003).

Stage 1: Modified Patankar-Euler step to obtain the predictor \(c^{(1)}\) (same system as modified_patankar()).

Stage 2: second-order corrector using averaged production/destruction terms \(\frac{1}{2}(pp + pp^{(1)})\) and \(\frac{1}{2}(dd + dd^{(1)})\), with \(c^{(1)}\) as the denominator state.

Updates cc in-place. Level 0 is never modified.

Return type:

None

Parameters:
pygotm.util.ode_solvers.modified_patankar_4(dt, numc, nlev, cc, get_ppdd)[source]

Fourth-order Modified Patankar-Runge-Kutta (MPRK4) — does not work.

This scheme has not yet been developed in GOTM Fortran or here; the implementation is a placeholder. Updates cc in-place. Level 0 is never modified.

Return type:

None

Parameters:
pygotm.util.ode_solvers.emp_1(dt, numc, nlev, cc, get_rhs)[source]

First-order Extended Modified Patankar (EMP-1) scheme.

One evaluation of the right-hand side per time step. Stoichiometrically conservative and positive (Bruggeman et al. 2005):

\[c^{n+1} = c^n + \Delta t\,f(t^n, c^n)\prod_{j \in J^n} \frac{c_j^{n+1}}{c_j^n}, \quad J^n = \{i : f_i(t^n, c^n) < 0\}\]

The product term \(p\) is found via findp_bisection(). Updates cc in-place. Level 0 is never modified.

Return type:

None

Parameters:
pygotm.util.ode_solvers.emp_2(dt, numc, nlev, cc, get_rhs)[source]

Second-order Extended Modified Patankar (EMP-2) scheme.

Two evaluations of the right-hand side per time step. Stoichiometrically conservative and positive (Bruggeman et al. 2005).

Step 1: identical to EMP-1 — advance cc to a midpoint cc_med using the product term from findp_bisection().

Step 2: average \(f(c^n)\) and \(f(c^{\mathrm{med}})\), correct for state variables in the negative-flux set \(J\), then apply a second findp_bisection() step to update cc.

Updates cc in-place. Level 0 is never modified.

Return type:

None

Parameters:
pygotm.util.ode_solvers.findp_bisection(numc, cc, derivative, dt, accuracy)[source]

Find the EMP product term p via bisection.

Solves the non-linear problem:

c^{n+1} = c^n + dt * f(t^n, c^n) * prod_{j in J} (c_j^{n+1} / c_j^n)

with J = {i : f_i < 0} by reducing it to a polynomial root problem and applying 20 bisection iterations. It has been proved that there exists exactly one p for which the above is true (Bruggeman et al. 2005).

Original FORTRAN author: Jorn Bruggeman.

Return type:

float

Parameters:

Parameters

numcint

Number of state variables.

ccnp.ndarray, shape (numc,)

Current concentrations at one vertical level.

derivativenp.ndarray, shape (numc,)

Right-hand side f evaluated at current state.

dtfloat

Time step [s].

accuracyfloat

Convergence criterion on the relative pi interval width.

Returns

pifloat

EMP product term p in (0, 1].

pygotm.util.ode_solvers.matrix_solve(n, a, r)[source]

Gaussian forward-elimination + back-substitution solver.

Solves the n×n system A·c = r. Modifies working copies of a and r internally; does not alter the caller’s arrays.

Return type:

ndarray

Parameters:

Parameters

nint

System dimension.

anp.ndarray, shape (n, n)

Coefficient matrix (modified internally; caller’s copy is unchanged).

rnp.ndarray, shape (n,)

Right-hand side vector (modified internally).

Returns

cnp.ndarray, shape (n,)

Solution vector.

ODE Solver Template

Translation of ode_solvers_template.F90. Provides alternate implementations of the RK2 and RK4 solvers that differ from ode_solvers:

  • RK2 (code 2) — explicit midpoint method:

    \[c^{\mathrm{mid}} = c^n + \tfrac{\Delta t}{2}\,f(c^n),\quad c^{n+1} = c^n + \Delta t\,f(c^{\mathrm{mid}})\]

  • RK4 (code 3) — standard four-stage RK with half-step intermediates:

    \[k_1 = f(c^n),\quad k_2 = f\!\left(c^n + \tfrac{\Delta t}{2}k_1\right),\quad k_3 = f\!\left(c^n + \tfrac{\Delta t}{2}k_2\right),\quad k_4 = f(c^n + \Delta t\,k_3)\]
    \[c^{n+1} = c^n + \tfrac{\Delta t}{3}\!\left(\tfrac{k_1}{2} + k_2 + k_3 + \tfrac{k_4}{2}\right)\]

All other solvers are re-exported unchanged from ode_solvers.

ODE solver template — translation of ode_solvers_template.F90.

Alternate implementations of the Runge-Kutta solvers (codes 2 and 3) that differ from ode_solvers:

  • RK2 (code 2): explicit midpoint method.

    \[c^{\mathrm{mid}} = c^n + \tfrac{\Delta t}{2}\,f(c^n),\quad c^{n+1} = c^n + \Delta t\,f(c^{\mathrm{mid}})\]

  • RK4 (code 3): standard four-stage RK with half-step intermediates.

    \[k_1 = f(c^n),\quad k_2 = f(c^n + \tfrac{\Delta t}{2}k_1),\quad k_3 = f(c^n + \tfrac{\Delta t}{2}k_2),\quad k_4 = f(c^n + \Delta t\,k_3)\]
    \[c^{n+1} = c^n + \tfrac{\Delta t}{3}\left(\tfrac{k_1}{2} + k_2 + k_3 + \tfrac{k_4}{2}\right)\]

All other solvers (Euler, Patankar, Modified Patankar, EMP — codes 1, 4–11) are re-exported unchanged from ode_solvers.

Original FORTRAN authors: Hans Burchard, Karsten Bolding.

pygotm.util.ode_solvers_template.ode_solver(solver, numc, nlev, dt, cc, get_rhs=None, get_ppdd=None)[source]

Dispatch to one of 11 ODE solvers (template version).

Identical to pygotm.util.ode_solvers.ode_solver() except that solvers 2 and 3 use the template-version algorithms: midpoint RK2 and standard half-step RK4. Solvers 4–11 delegate to the shared implementations in ode_solvers.

Return type:

None

Parameters:

Parameters

solverint

Solver identifier (1–11).

numcint

Number of biogeochemical state variables.

nlevint

Number of vertical levels.

dtfloat

Time step [s].

ccnp.ndarray, shape (numc, nlev+1)

State variable array, modified in-place. Level 0 is not updated.

get_rhscallable, optional

Required for solvers 1, 2, 3, 10, 11.

get_ppddcallable, optional

Required for solvers 4–9.

pygotm.util.ode_solvers_template.euler_forward(dt, numc, nlev, cc, get_rhs)[source]

First-order Euler-forward (E1) scheme.

One evaluation of the right-hand side per time step:

\[c_i^{n+1} = c_i^n + \Delta t\,\bigl(P_i(c^n) - D_i(c^n)\bigr)\]

Updates cc in-place. Level 0 (index 0) is never modified.

Return type:

None

Parameters:
pygotm.util.ode_solvers_template.runge_kutta_2(dt, numc, nlev, cc, get_rhs)[source]

Second-order Runge-Kutta — explicit midpoint method (template version).

Two evaluations of the right-hand side per time step:

\[\begin{split}c^{\mathrm{mid}} &= c^n + \tfrac{\Delta t}{2}\,f(c^n) \\ c^{n+1} &= c^n + \Delta t\,f(c^{\mathrm{mid}})\end{split}\]

This is the ode_solvers_template.F90 variant. The non-template ode_solvers.F90 uses Heun’s trapezoidal method instead. Updates cc in-place. Level 0 is never modified.

Return type:

None

Parameters:
pygotm.util.ode_solvers_template.runge_kutta_4(dt, numc, nlev, cc, get_rhs)[source]

Fourth-order Runge-Kutta — standard half-step variant (template version).

Four evaluations of the right-hand side per time step with standard half-step intermediates (achieves true 4th-order accuracy):

\[\begin{split}k_1 &= f(c^n),\quad k_2 = f\!\left(c^n + \tfrac{\Delta t}{2}k_1\right),\quad k_3 = f\!\left(c^n + \tfrac{\Delta t}{2}k_2\right),\quad k_4 = f(c^n + \Delta t\,k_3) \\ c^{n+1} &= c^n + \tfrac{\Delta t}{3}\bigl(\tfrac{k_1}{2} + k_2 + k_3 + \tfrac{k_4}{2}\bigr)\end{split}\]

This is the ode_solvers_template.F90 variant. The non-template ode_solvers.F90 uses full-\(\Delta t\) intermediate steps instead. Updates cc in-place. Level 0 is never modified.

Return type:

None

Parameters:
pygotm.util.ode_solvers_template.patankar(dt, numc, nlev, cc, get_ppdd)[source]

First-order Patankar-Euler (PE1) scheme.

One evaluation of the production/destruction terms per time step. Not conservative; unconditionally positive:

\[c_i^{n+1} = \frac{c_i^n + \Delta t\,P_i(c^n)}{1 + \Delta t\,D_i(c^n)/c_i^n}\]

Updates cc in-place. Level 0 is never modified.

Return type:

None

Parameters:
pygotm.util.ode_solvers_template.patankar_runge_kutta_2(dt, numc, nlev, cc, get_ppdd)[source]

Second-order Patankar-Runge-Kutta (PRK2) scheme.

Two evaluations of the production/destruction terms per time step. Not conservative; unconditionally positive (Burchard et al. 2003).

Stage 1 — Patankar-Euler predictor:

\[c_i^{(1)} = \frac{c_i^n + \Delta t\,P_i(c^n)}{1 + \Delta t\,D_i(c^n)/c_i^n}\]

Stage 2 — second-order corrector:

\[c_i^{n+1} = \frac{c_i^n + \tfrac{\Delta t}{2}(P_i(c^n)+P_i(c^{(1)}))} {1 + \tfrac{\Delta t}{2}(D_i(c^n)+D_i(c^{(1)}))/c_i^{(1)}}\]

Updates cc in-place. Level 0 is never modified.

Return type:

None

Parameters:
pygotm.util.ode_solvers_template.patankar_runge_kutta_4(dt, numc, nlev, cc, get_ppdd)[source]

Fourth-order Patankar-Runge-Kutta (PRK4) scheme — does not work.

This scheme has not yet been developed in GOTM Fortran or here; the implementation is a placeholder. Updates cc in-place. Level 0 is never modified.

Return type:

None

Parameters:
pygotm.util.ode_solvers_template.modified_patankar(dt, numc, nlev, cc, get_ppdd)[source]

First-order Modified Patankar-Euler (MPE1) scheme.

One evaluation of the production/destruction terms per time step. Conservative and unconditionally positive (Burchard et al. 2003):

\[c_i^{n+1} = c_i^n + \Delta t\left[ \sum_{j \neq i} p_{ij}(c^n)\,\frac{c_j^{n+1}}{c_j^n} + p_{ii}(c^n) - \sum_j d_{ij}(c^n)\,\frac{c_i^{n+1}}{c_i^n} \right]\]

Solved as a linear system per vertical level using matrix_solve(). Updates cc in-place. Level 0 is never modified.

Return type:

None

Parameters:
pygotm.util.ode_solvers_template.modified_patankar_2(dt, numc, nlev, cc, get_ppdd)[source]

Second-order Modified Patankar-Runge-Kutta (MPRK2) scheme.

Two evaluations of the production/destruction terms per time step. Conservative and unconditionally positive (Burchard et al. 2003).

Stage 1: Modified Patankar-Euler step to obtain the predictor \(c^{(1)}\) (same system as modified_patankar()).

Stage 2: second-order corrector using averaged production/destruction terms \(\frac{1}{2}(pp + pp^{(1)})\) and \(\frac{1}{2}(dd + dd^{(1)})\), with \(c^{(1)}\) as the denominator state.

Updates cc in-place. Level 0 is never modified.

Return type:

None

Parameters:
pygotm.util.ode_solvers_template.modified_patankar_4(dt, numc, nlev, cc, get_ppdd)[source]

Fourth-order Modified Patankar-Runge-Kutta (MPRK4) — does not work.

This scheme has not yet been developed in GOTM Fortran or here; the implementation is a placeholder. Updates cc in-place. Level 0 is never modified.

Return type:

None

Parameters:
pygotm.util.ode_solvers_template.emp_1(dt, numc, nlev, cc, get_rhs)[source]

First-order Extended Modified Patankar (EMP-1) scheme.

One evaluation of the right-hand side per time step. Stoichiometrically conservative and positive (Bruggeman et al. 2005):

\[c^{n+1} = c^n + \Delta t\,f(t^n, c^n)\prod_{j \in J^n} \frac{c_j^{n+1}}{c_j^n}, \quad J^n = \{i : f_i(t^n, c^n) < 0\}\]

The product term \(p\) is found via findp_bisection(). Updates cc in-place. Level 0 is never modified.

Return type:

None

Parameters:
pygotm.util.ode_solvers_template.emp_2(dt, numc, nlev, cc, get_rhs)[source]

Second-order Extended Modified Patankar (EMP-2) scheme.

Two evaluations of the right-hand side per time step. Stoichiometrically conservative and positive (Bruggeman et al. 2005).

Step 1: identical to EMP-1 — advance cc to a midpoint cc_med using the product term from findp_bisection().

Step 2: average \(f(c^n)\) and \(f(c^{\mathrm{med}})\), correct for state variables in the negative-flux set \(J\), then apply a second findp_bisection() step to update cc.

Updates cc in-place. Level 0 is never modified.

Return type:

None

Parameters:
pygotm.util.ode_solvers_template.findp_bisection(numc, cc, derivative, dt, accuracy)[source]

Find the EMP product term p via bisection.

Solves the non-linear problem:

c^{n+1} = c^n + dt * f(t^n, c^n) * prod_{j in J} (c_j^{n+1} / c_j^n)

with J = {i : f_i < 0} by reducing it to a polynomial root problem and applying 20 bisection iterations. It has been proved that there exists exactly one p for which the above is true (Bruggeman et al. 2005).

Original FORTRAN author: Jorn Bruggeman.

Return type:

float

Parameters:

Parameters

numcint

Number of state variables.

ccnp.ndarray, shape (numc,)

Current concentrations at one vertical level.

derivativenp.ndarray, shape (numc,)

Right-hand side f evaluated at current state.

dtfloat

Time step [s].

accuracyfloat

Convergence criterion on the relative pi interval width.

Returns

pifloat

EMP product term p in (0, 1].

pygotm.util.ode_solvers_template.matrix_solve(n, a, r)[source]

Gaussian forward-elimination + back-substitution solver.

Solves the n×n system A·c = r. Modifies working copies of a and r internally; does not alter the caller’s arrays.

Return type:

ndarray

Parameters:

Parameters

nint

System dimension.

anp.ndarray, shape (n, n)

Coefficient matrix (modified internally; caller’s copy is unchanged).

rnp.ndarray, shape (n,)

Right-hand side vector (modified internally).

Returns

cnp.ndarray, shape (n,)

Solution vector.

Version and Compilation

Two stub modules expose build-time metadata from the original Fortran GOTM:

  • gotm_version — git_commit_id = "4.1.0" and git_branch_name = "master".

  • compilation — compiler, compiler_id, compiler_version (all empty strings; injected by CMake in the Fortran build but not applicable in Python).

GOTM version constants — translation of gotm_version.F90.

Exposes the GOTM source version identifiers that were embedded at compile time in the original Fortran build. Values are fixed strings matching the GOTM release this Python translation targets.

Public interface: git_commit_id, git_branch_name.

Compilation metadata stub — translation of gotm_compilation.F90.

In the original Fortran build, compiler, compiler_id, and compiler_version were injected by the CMake build system. This Python stub exposes the same names as empty strings; they carry no functional meaning at runtime.

Public interface: compiler, compiler_id, compiler_version.