Discretization and Operators

ClimaAtmos discretizes the governing equations with a hybrid scheme: a Galerkin spectral element method in the horizontal and finite differences on a staggered grid in the vertical [1]. The horizontal method is continuous Galerkin (CG) today; ClimaCore also provides discontinuous Galerkin (DG) on the same element structure, and the two differ only in how element-local results are completed across element boundaries, so this page is written to cover both and marks the one place they part ways.

The page explains the choice of scheme, states the two rules that decide which operator each term uses, defines the operator symbols used throughout these docs, and shows how the pieces assemble into a discretized equation. Readers who need only the symbols can go straight to the operator reference; those adding a term to the model want Writing a new tendency.

The operators come from ClimaCore, whose Mathematical framework, Spectral elements: continuous and discontinuous Galerkin, Staggered vertical discretization, and DSS and numerical fluxes pages define the nodal basis, the strong and weak forms, the staggered vertical operators, and the completion step (direct stiffness summation for CG, interface fluxes for DG). This page states the choices ClimaAtmos makes among those operators — which form each term uses, and which reconstruction — and gives the symbol used for each one in the equations. For the short names the tendencies are written in, see Discrete operators in the API.

Why a hybrid scheme

The horizontal and vertical directions of an atmospheric grid have different requirements, and the scheme treats them differently.

Horizontally, spectral elements give high-order accuracy with one neighbor exchange per step; vertically, finite differences on a staggered grid couple only neighboring levels, so the fast vertical waves that fine vertical spacing admits can be solved implicitly column by column, and the staggering suppresses the computational modes of an unstaggered vertical grid. ClimaCore's Staggered vertical discretization gives the argument in full, including why Lorenz staggering is used. See Implicit Solver for the implicit solve.

The two parts are independent. The vertical staggering, the reconstruction rules, and the implicit solve are independent of how horizontal derivatives are computed, and the assembled equations at the end of this page are written in operator symbols rather than in any one method's terms. An alternative horizontal discretization — a discontinuous Galerkin method, for example — would supply its own realization of the horizontal operators and of the projection $\mathcal{P}$, and leave the rest of this page, and of the model, unchanged.

Grid layout and staggering

The domain is divided into $N_h$ horizontal elements, each extruded into $N_v$ vertical layers. Horizontal elements carry $(N_p + 1)^2$ Gauss–Lobatto–Legendre nodal points, so fields are polynomials of order $N_p$ within an element. On the sphere, the horizontal mesh is an equiangular cubed sphere [22, 23]; in Cartesian geometry the same machinery discretizes a box. See Grids for the constructors and Topography Representation for the terrain-following vertical coordinate.

The vertical arrangement is a Lorenz staggering: the covariant vertical velocity component $u_3$ is defined on element faces, and every other variable — including the horizontal velocity components $u_1$ and $u_2$ — is defined on element centers [1].

Throughout the docs and the code, a $ᶜ$ prefix marks a center field and a $ᶠ$ prefix marks a face field; see Notation and Symbols for the full convention.

The horizontal discretization: spectral elements

Within each element, fields are expanded in a nodal polynomial basis on the Gauss–Lobatto–Legendre points, and horizontal derivatives are computed by differentiating that expansion. The element-local operators, and the rule below for choosing between their strong and weak forms, are the same for continuous and discontinuous Galerkin; the two methods differ only in the completion step at the end of this section. The vertical discretization does not depend on either.

Strong and weak forms: which to use

Every horizontal derivative comes in a strong form, which differentiates the basis functions directly, and a weak form, constructed to be the negative adjoint of the strong form under the discrete inner product — a discrete integration by parts, from which discrete divergence and Stokes theorems follow. ClimaCore's Spectral elements page derives both; what matters here is which one each term uses.

The rule

Adjointness is a property of a pair. A term conserves what it should when the two operators acting in it are adjoint to one another, so that summing over the domain telescopes into boundary terms alone. That gives a short rule:

Mixing the two forms is therefore deliberate: using the weak divergence for the mass flux and the strong gradient for kinetic energy is what makes the two terms cancel exactly in the energy budget. Using the same form for both would break that cancellation and leave a spurious energy source. See Conservation Properties.

Note the asymmetry between the two Laplacians. Both put one strong and one weak operator in the pair, but the scalar Laplacian takes the divergence weakly, while the vector Laplacian takes the outer operator weakly. Either arrangement satisfies the identity, because each is the adjoint of the other; what matters is that the pair contains one of each. Fourth-order hyperdiffusion applies the Laplacian twice with direct stiffness summation in between, which removes inter-element discontinuities without disturbing the inner product; see Hyperdiffusion.

In the code

The strong operators are divₕ, gradₕ, and curlₕ; the weak ones take a w prefix: wdivₕ, wgradₕ, wcurlₕ; the split divergence is split_divₕ. So the scalar and vector Laplacians read

ᶜ∇²s_d = @. wdivₕ(gradₕ(s_d - sd_r))                           # weak div ∘ strong grad
ᶜ∇²u = @. C123(wgradₕ(divₕ(ᶜu))) - C123(wcurlₕ(C123(curlₕ(ᶜu))))

The variable-resolution behavior of these operators is examined by [26].

Completing a tendency across elements

An element-local operator evaluates each element on its own, so at a node on an element boundary the weak-form result is incomplete: it holds only that element's share. A completion step, written $\mathcal{P}$ in these docs, supplies the rest, and it is the one place where continuous and discontinuous Galerkin differ; ClimaCore's DSS and numerical fluxes page develops both.

With continuous elements, which ClimaAtmos uses, $\mathcal{P}$ is direct stiffness summation (DSS): the volume-weighted average over the copies of each boundary node, which makes the field single-valued there. It preserves the discrete inner product, so it leaves the conservation properties of the weak operators intact, and it is the only operation that communicates horizontally between elements. ClimaAtmos applies it several times per step — at stage boundaries, around the implicit solve, and at the end of the step — which keeps round-off errors from accumulating as discontinuities across element boundaries.

With discontinuous elements, fields may jump at element boundaries, and $\mathcal{P}$ instead adds an interface numerical flux (a central or Rusanov flux, or a user-supplied one) to each conservative term and a lifting flux to each gradient or curl. The element-local operators and the choice of form are untouched; only the completion changes, and the discrete conservation statements above carry over because the numerical flux is single-valued and antisymmetric across each interface.

The vertical discretization: interpolating between centers and faces

Staggering forces interpolation. The covariant vertical velocity $u_3$ is defined on faces and everything else on centers, so any flux term needs one or the other moved. As with the strong and weak forms, the choice of average decides whether a discrete conservation law holds.

Three kinds appear.

Arithmetic mean ($I^c$, $I^f$; ᶜinterp, ᶠinterp). The plain average of the two neighbors. Use it for quantities that are not weighted by mass, such as the covariant velocity components, and wherever no conservation statement rides on the result.

Mass-weighted average ($WI^f$; ᶠwinterp). The average weighted by $\rho J$, that is, $WI^f(w, x) = I^f(w x) / I^f(w)$. Use it for a quantity that will be multiplied by a mass flux. The arithmetic and mass-weighted averages form an adjoint pair satisfying a density-weighted averaging-by-parts identity, so the vertical flux divergence telescopes; using two arithmetic means instead would leave a residual. In the code, it appears where the horizontal velocity is reconstructed onto faces, and in the hyperdiffusive momentum tendency. The face mass flux itself is not this operator: it is $I^f(\rho J) \tilde{\boldsymbol{u}} / J^f$, divided by the face Jacobian rather than by $I^f(J)$, so that it telescopes exactly.

Upwind or limited reconstruction ($U^f$; ᶠupwind1, ᶠupwind3, ᶠlin_vanleer). A biased or flux-corrected reconstruction [27–29]. Use it for advected scalars, where a centered average would produce dispersive oscillations and negative concentrations. The default for grid-mean energy and tracer transport is the van Leer limiter of [28], constrained by the local extrema of the neighboring center values.

Two consequences follow. Setting $\psi = 1$ in the scalar flux reconstruction recovers the mass flux divergence, so tracer transport stays consistent with mass transport, and a uniform tracer field stays uniform. At the domain boundaries, $I^f$ extrapolates by reusing the nearest interior value, while the vertical gradient and curl operators are set to zero there.

Operator reference

Each operator below has a short name in the ClimaAtmos source, documented under Discrete operators in the API.

Reconstruction between centers and faces

• The face-to-center interpolation $I^c$ is ClimaCore.Operators.InterpolateF2C, an arithmetic mean.

• The center-to-face interpolation $I^f$ is ClimaCore.Operators.InterpolateC2F, an arithmetic mean with constant extrapolation to the domain boundaries.

• The center-to-face weighted interpolation $WI^f$ is ClimaCore.Operators.WeightedInterpolateC2F, with $WI^f(w, x) = I^f(w x) / I^f(w)$ for a weight $w$. It appears with the weight $\rho J$, to reconstruct the horizontal velocity onto faces and in the hyperdiffusive vertical-momentum tendency.

• The face mass flux $\mathcal{M}^f(\rho)$ is

  \[\mathcal{M}^f(\rho) = \frac{I^f(\rho J)}{J^f} \, \tilde{\boldsymbol{u}} ,\]

  with $J^f$ the Jacobian on faces. Dividing by $J^f$ rather than by $I^f(J)$ makes the flux divergence telescope: $D^c$ multiplies its argument by $J^f$ before differencing, so the $J^f$ cancels and the flux that leaves one cell is exactly the flux that enters the next. The two denominators agree only where $J^f = I^f(J)$, which holds on a flat or linearly warped grid but not under a SLEVE warp; see Topography Representation.

• The center-to-face upwind product $U^f$: first order ClimaCore.Operators.UpwindBiasedProductC2F, third order ClimaCore.Operators.Upwind3rdOrderBiasedProductC2F, or the van Leer limiter ClimaCore.Operators.LinVanLeerC2F. The van Leer limiter is the default for grid-mean energy and tracer vertical transport, set by the energy_q_tot_upwinding and tracer_upwinding configuration keys.

Horizontal differential operators

These are the element-local strong and weak forms of the previous section, and they are the same operators for continuous and discontinuous elements; only $\mathcal{P}$ below changes between the two.

• The strong horizontal spectral divergence $\mathcal{D}_h$ is ClimaCore.Operators.Divergence.

• The weak horizontal spectral divergence $\hat{\mathcal{D}}_h$ is Divergence{WeakForm} (see ClimaCore.Operators.Divergence).

• The split, skew-symmetric horizontal divergence $\mathcal{D}^{split}_h$ is ClimaCore.Operators.SplitDivergence,

  \[\mathcal{D}^{split}_h(\rho \boldsymbol{u}, \psi) = \tfrac{1}{2} \hat{\mathcal{D}}_h(\rho \boldsymbol{u} \psi) + \tfrac{1}{2} \left[ \psi \, \hat{\mathcal{D}}_h(\rho \boldsymbol{u}) + \rho \boldsymbol{u} \cdot \mathcal{G}_h \psi \right].\]

  The horizontal advective fluxes of energy, moisture, and tracers use this entropy-stable split form; for $\psi = 1$ it reduces to the weak divergence, which is why the semi-discrete mass equation below is written with $\hat{\mathcal{D}}_h$ (in the code, the mass flux is split_divₕ(ρu, 1)). The horizontal pressure-gradient term uses an analogous split form.

• The strong horizontal spectral gradient $\mathcal{G}_h$ is ClimaCore.Operators.Gradient.

• The weak horizontal spectral gradient $\hat{\mathcal{G}}_h$ is Gradient{WeakForm} (see ClimaCore.Operators.Gradient), the outer gradient of the vector Laplacian.

• The curl $\mathcal{C}_h$ acts on the components involving horizontal derivatives ClimaCore.Operators.Curl. Applied to $\boldsymbol{u}_h$ it returns a vector with only vertical contravariant components; applied to $\boldsymbol{u}_v$ it returns a vector with only horizontal contravariant components.

• The corresponding weak curl $\hat{\mathcal{C}}_h$ is Curl{WeakForm} (see ClimaCore.Operators.Curl).

• The completion step $\mathcal{P}$ across element boundaries is direct stiffness summation on the continuous elements ClimaAtmos uses; see Completing a tendency across elements.

Vertical differential operators

• The face-to-center vertical divergence $D^c$ is ClimaCore.Operators.DivergenceF2C. Separate variants impose the boundary conditions for advective, precipitation, and diffusive fluxes.

• The center-to-face vertical gradient $G^f$ is ClimaCore.Operators.GradientC2F, set to zero at the top and bottom boundaries. The boundary values are placeholders: the vertical velocity at the boundaries is fixed instead by the impenetrability condition described above.

• The face-to-center vertical gradient $G^c$ is ClimaCore.Operators.GradientF2C, used where a face-defined field such as the geopotential is differentiated onto centers.

• The center-to-face curl $C^f$ acts on the components involving vertical derivatives ClimaCore.Operators.CurlC2F, set to zero at the top and bottom boundaries. Applied to $\boldsymbol{u}_h$ it returns a vector with only a horizontal contravariant component.

Reconstructions and conservation

The reconstructions above are the ones that make the global conservation laws hold [1]; they resemble those of [30] in some respects.

• Density is reconstructed onto faces as $I^f(\rho J) / J^f$, the Jacobian-weighted average that the face mass flux needs.
• Velocity covariant components use unweighted averages; the contravariant vertical component on faces uses a mass-weighted average. There is no unique reconstruction of the contravariant vertical velocity onto centers.
• Vorticity contravariant components are computed with a weak horizontal curl.
• Momentum advection uses the vector-invariant form, with a strong horizontal gradient of kinetic energy and a weighted average of the vorticity term. In this form, momentum advection conserves kinetic energy and vorticity globally, and avoids the curvature terms that appear in advection terms in non-orthogonal coordinates.

Because total energy is separately conserved, any numerical conversion between kinetic and non-kinetic energy comes from the discretized pressure-gradient term and the physical sources, and not from the advection scheme.

Writing a new tendency

Putting the two rules together, a new term added to src/prognostic_equations/ generally follows this pattern.

1. Decide where the result lives. A tendency for a center variable must end on centers, one for u₃ on faces.
2. For a horizontal advective flux of a conserved scalar, use the split form split_divₕ; for other horizontal flux divergences, such as diffusive ones, use wdivₕ. For a gradient in the momentum equation, use gradₕ.
3. For the vertical part, use ᶜdivᵥ on a face flux, and build that face flux with ᶠwinterp if it multiplies a mass flux, or with an upwind operator if it transports a scalar.
4. If the term diffuses energy, diffuse the dry static energy and the effective total water separately rather than a lumped total enthalpy, so the decomposition stays energetically consistent. See Thermodynamics and the Working Fluid.
5. If the term is vertical and fast, add it to the implicit tendency and to the Jacobian; see Implicit Solver.

The semi-discrete equations at the end of this page show these patterns applied to each governing equation.

Timestepping

Tendencies are split into an explicit part and an implicit part and advanced with a horizontally explicit, vertically implicit (HEVI) additive Runge–Kutta method [24, 31]. The implicit part takes the vertical terms responsible for sound and gravity waves, falling and sedimenting condensate, the Rayleigh damping of $u_3$, and, when implicit_diffusion is enabled, the vertical diffusion; the viscous sponge stays explicit. Because the implicit part involves no horizontal derivatives, it can be solved independently in each column, with no horizontal communication.

This lifts the timestep restriction from fast vertical dynamics and leaves the horizontal propagation of sound waves as the limit, so the maximum timestep scales as $\delta t \sim (\delta x)_{\min} / c_s$ with the minimum horizontal distance between nodal points and the speed of sound.

Implicit Solver documents the Newton solve, the Jacobian approximation, and the available Jacobian algorithms. Integer Time (ITime) explains how time itself is represented.

Assembling a discretized equation

The operators above combine in the same way in every prognostic equation. Two reconstructed velocities and the kinetic energy built from them appear throughout; all follow from the staggering, since the covariant vertical component lives on faces and everything else on centers:

\[\tilde{\boldsymbol{u}} = WI^f(\rho J, \boldsymbol{u}_h) + \boldsymbol{u}_v , \qquad \bar{\boldsymbol{u}} = \boldsymbol{u}_h + I^c(\boldsymbol{u}_v) , \qquad K = \tfrac{1}{2} \left( \boldsymbol{u}_h \cdot \boldsymbol{u}_h + 2 \boldsymbol{u}_h \cdot I^c(\boldsymbol{u}_v) + I^c(\boldsymbol{u}_v \cdot \boldsymbol{u}_v) \right) .\]

$\tilde{\boldsymbol{u}}$ is the face velocity (see compute_ᶠuₕ³ in src/cache/precomputed_quantities.jl), $\bar{\boldsymbol{u}}$ the center velocity, and $K$ the specific kinetic energy at centers. The no-flux condition at the surface and the model top fixes the covariant vertical component from the horizontal one, $u_3 = -(g^{31} u_1 + g^{32} u_2)/g^{33}$, so that the contravariant $\tilde{u}^3$ vanishes there; see Topography.

One equation in full: total energy

Total energy shows every pattern the other equations reuse. Its continuous form is

\[\frac{\partial}{\partial t} \rho e_{tot} = - \nabla \cdot \left[ (\rho e_{tot} + p) \boldsymbol{u} + \boldsymbol{F}_R \right] + \rho S_e ,\]

stabilized by fourth-order hyperdiffusion of the total enthalpy, decomposed into dry static energy and a water-enthalpy piece (see Hyperdiffusion). Discretized, it reads

\[\frac{\partial}{\partial t} \rho e_{tot} \approx - \mathcal{D}^{split}_h \left[ \rho \bar{\boldsymbol{u}}, \tfrac{\rho e_{tot} + p}{\rho} \right] - D^c \left[ \mathcal{M}^f(\rho) \, I^f \left( \tfrac{\rho e_{tot} + p}{\rho} \right) \right] - D^c \left[ \boldsymbol{F}_R \right] - \nu_h \left[ \hat{\mathcal{D}}_h \left( \rho \, \mathcal{G}_h(\psi_{s_d}) \right) + \hat{\mathcal{D}}_h \left( \rho \, h_{\mathrm{tot,cl}} \, \mathcal{G}_h(\psi_{q_t^{\mathrm{eff}}}) \right) \right], \qquad \psi_x = \mathcal{P} \left[ \hat{\mathcal{D}}_h \left( \mathcal{G}_h (x - x_r) \right) \right],\]

with $x_r$ the hydrostatic reference profile subtracted before the first Laplacian. Reading the terms in order:

• The horizontal advective flux uses the split divergence $\mathcal{D}^{split}_h$ of the center velocity, the entropy-stable form from the rule above.
• The vertical advective flux is the face mass flux $\mathcal{M}^f(\rho)$ times the centrally interpolated enthalpy, differenced with $D^c$. This term is implicit. When upwinding is enabled for energy (the van Leer limiter by default), the implicit solve still uses this central form, and the difference between the upwinded and central fluxes is applied after the Newton solve as a T_post_imp! correction (correct_implicit_advection_tendency!), so the upwind flux sees the Newton-solved velocity.
• The radiative flux divergence is explicit and applied separately.
• Hyperdiffusion is two Laplacians with the completion $\mathcal{P}$ in between, each a weak divergence of a strong gradient.

The other equations

The remaining equations follow the same template; the table gives the operator each term uses and where it lives. Everything in the "Implicit" column is solved column by column as described under Timestepping; everything else is explicit.

Two details of the momentum rows are worth spelling out. The pressure-gradient and geopotential terms in the momentum equation are written as a departure from a hydrostatic reference state $(\theta_{v,r}, \Phi_r)$, with $\theta_v' = \theta_v - \theta_{v,r}$; the horizontal part uses the split form ``\tfrac{c{pd}}{2}{\thetav' \mathcal{G}h[\Pi] + \mathcal{G}h[\theta_v' \Pi]

• \Pi \mathcal{G}h[\thetav']}$, and the reference profile itself is on the [Governing Equations](equations.md) page. And the momentum hyperviscosity is the vector Laplacian of the rule table applied twice,$-\nuu { \delta{div} \hat{\mathcal{G}}h(\mathcal{D}h \boldsymbol{\psi})
• \hat{\mathcal{C}}h(\mathcal{C}h \boldsymbol{\psi}) }$with$\boldsymbol{\psi} = \mathcal{P}[\hat{\mathcal{G}}h(\mathcal{D}h \bar{\boldsymbol{u}})
• \hat{\mathcal{C}}h(\mathcal{C}h \bar{\boldsymbol{u}})]$, projected onto the horizontal covariant directions for$\boldsymbol{u}h$and onto the vertical one, after a$\rho J$-weighted interpolation to faces, for$u3$; the divergence damping factor$\delta_{div}`` is described under Hyperdiffusion.

The fully expanded forms of each tendency are in the source files linked above and in Appendix B of [1].

Where this is implemented

The upwinding and limiter choices are exposed as configuration keys; see Configuration Options for energy_q_tot_upwinding, tracer_upwinding, and the corresponding PROPHET keys.