Vertical coordinates

SpeedyWeather.jl supports two families of terrain-following vertical coordinates. Both are implemented as subtypes of AbstractVerticalCoordinates and are interchangeable in the model constructor. For more mathematical background on coordinate transformations see Vertical coordinates in the Primitive equation model documentation.

Sigma coordinates

Sigma coordinates use the fraction of surface pressure as the vertical coordinate

with   at the model top and   at the surface. Sigma levels are terrain-following because   is always at surface pressure, bending around every mountain. One specifies the half levels and the full levels are obtained as midpoints

Creating sigma coordinates

using SpeedyWeather
spectral_grid = SpectralGrid(nlayers = 8)
SigmaCoordinates(spectral_grid)

8-layer SigmaCoordinates{Int64, Vector{Float32}}
 ├─ 0.0000  k = 0.5
 │× 0.0625  k = 1
 ├─ 0.1250  k = 1.5
 │× 0.1875  k = 2
 ├─ 0.2500  k = 2.5
 │× 0.3125  k = 3
 ├─ 0.3750  k = 3.5
 │× 0.4375  k = 4
 ├─ 0.5000  k = 4.5
 │× 0.5625  k = 5
 ├─ 0.6250  k = 5.5
 │× 0.6875  k = 6
 ├─ 0.7500  k = 6.5
 │× 0.8125  k = 7
 ├─ 0.8750  k = 7.5
 │× 0.9375  k = 8
 └─ 1.0000  k = 8.5

From a custom vector of half levels (nlayers is inferred from the length):

SigmaCoordinates([0.0, 0.1, 0.25, 0.45, 0.62, 0.78, 0.9, 0.97, 1.0])

8-layer SigmaCoordinates{Int64, Vector{Float32}}
 ├─ 0.0000  k = 0.5
 │× 0.0500  k = 1
 ├─ 0.1000  k = 1.5
 │× 0.1750  k = 2
 ├─ 0.2500  k = 2.5
 │× 0.3500  k = 3
 ├─ 0.4500  k = 3.5
 │× 0.5350  k = 4
 ├─ 0.6200  k = 4.5
 │× 0.7000  k = 5
 ├─ 0.7800  k = 5.5
 │× 0.8400  k = 6
 ├─ 0.9000  k = 6.5
 │× 0.9350  k = 7
 ├─ 0.9700  k = 7.5
 │× 0.9850  k = 8
 └─ 1.0000  k = 8.5

From a range:

SigmaCoordinates(0:0.2:1)

5-layer SigmaCoordinates{Int64, Vector{Float32}}
 ├─ 0.0000  k = 0.5
 │× 0.1000  k = 1
 ├─ 0.2000  k = 1.5
 │× 0.3000  k = 2
 ├─ 0.4000  k = 2.5
 │× 0.5000  k = 3
 ├─ 0.6000  k = 3.5
 │× 0.7000  k = 4
 ├─ 0.8000  k = 4.5
 │× 0.9000  k = 5
 └─ 1.0000  k = 5.5

Frierson (2006) spacing with finer resolution near the surface and stratosphere:

FriersonSigmaCoordinates(spectral_grid)

8-layer SigmaCoordinates{Int64, Vector{Float32}}
 ├─ 0.0000  k = 0.5
 │× 0.0167  k = 1
 ├─ 0.0333  k = 1.5
 │× 0.0726  k = 2
 ├─ 0.1118  k = 2.5
 │× 0.1900  k = 3
 ├─ 0.2682  k = 3.5
 │× 0.3778  k = 4
 ├─ 0.4874  k = 4.5
 │× 0.5981  k = 5
 ├─ 0.7088  k = 5.5
 │× 0.7905  k = 6
 ├─ 0.8722  k = 6.5
 │× 0.9162  k = 7
 ├─ 0.9603  k = 7.5
 │× 0.9801  k = 8
 └─ 1.0000  k = 8.5

Using sigma coordinates in a model

Pass the coordinate to Geometry and then to the model constructor:

σ = FriersonSigmaCoordinates(spectral_grid)
geometry = Geometry(spectral_grid; vertical_coordinates = σ)
model = PrimitiveDryModel(spectral_grid; geometry)
simulation = initialize!(model)

Hybrid sigma-pressure coordinates

Work in progress

Hybrid sigma-pressure coordinates are implemented for the coordinate geometry, the pressure, pressure_thickness, and sigma functions, and all parameterizations, but the dynamical core still uses the pure sigma formulation throughout (vertical advection, surface pressure tendency, geopotential, and adiabatic conversion). See the TODO comments in those source files for details. Using SigmaPressureCoordinates with a transition close to pure sigma (e.g. the default  ) therefore remains consistent with the dynamics.

Pure sigma coordinates tilt sharply with orography even at high altitudes where terrain-following levels are not needed. Hybrid sigma-pressure coordinates solve this by blending constant-pressure surfaces near the model top with terrain-following sigma surfaces near the surface.

The pressure at layer is

where is a constant reference pressure and , are layer coefficients. Given a nominal sigma level and a transition function   the coefficients are

so that    always holds. When   the level is a pure constant-pressure surface; when   it is a pure sigma surface. The layer thickness in pressure is

with    and equivalently for .

Creating hybrid sigma-pressure coordinates

The default transition is  , giving a linear blend:

SigmaPressureCoordinates(spectral_grid)

8-layer SigmaPressureCoordinates{Float32, Vector{Float32}}
    p_ref = 100000.0 Pa   ░ pressure │ sigma █
 ├─ 0.0000  k = 0.5
 │× 0.0625  k = 1   █░░░░░░░░░  A=0.0547  B=0.0078
 ├─ 0.1250  k = 1.5
 │× 0.1875  k = 2   ██░░░░░░░░  A=0.1484  B=0.0391
 ├─ 0.2500  k = 2.5
 │× 0.3125  k = 3   ███░░░░░░░  A=0.2109  B=0.1016
 ├─ 0.3750  k = 3.5
 │× 0.4375  k = 4   ████░░░░░░  A=0.2422  B=0.1953
 ├─ 0.5000  k = 4.5
 │× 0.5625  k = 5   ██████░░░░  A=0.2422  B=0.3203
 ├─ 0.6250  k = 5.5
 │× 0.6875  k = 6   ███████░░░  A=0.2109  B=0.4766
 ├─ 0.7500  k = 6.5
 │× 0.8125  k = 7   ████████░░  A=0.1484  B=0.6641
 ├─ 0.8750  k = 7.5
 │× 0.9375  k = 8   █████████░  A=0.0547  B=0.8828
 └─ 1.0000  k = 8.5

Pure pressure levels everywhere ( ):

SigmaPressureCoordinates(spectral_grid; transition = _ -> 0.0)

8-layer SigmaPressureCoordinates{Float32, Vector{Float64}}
    p_ref = 100000.0 Pa   ░ pressure │ sigma █
 ├─ 0.0000  k = 0.5
 │× 0.0625  k = 1   ░░░░░░░░░░  A=0.0625  B=0.0000
 ├─ 0.1250  k = 1.5
 │× 0.1875  k = 2   ░░░░░░░░░░  A=0.1875  B=0.0000
 ├─ 0.2500  k = 2.5
 │× 0.3125  k = 3   ░░░░░░░░░░  A=0.3125  B=0.0000
 ├─ 0.3750  k = 3.5
 │× 0.4375  k = 4   ░░░░░░░░░░  A=0.4375  B=0.0000
 ├─ 0.5000  k = 4.5
 │× 0.5625  k = 5   ░░░░░░░░░░  A=0.5625  B=0.0000
 ├─ 0.6250  k = 5.5
 │× 0.6875  k = 6   ░░░░░░░░░░  A=0.6875  B=0.0000
 ├─ 0.7500  k = 6.5
 │× 0.8125  k = 7   ░░░░░░░░░░  A=0.8125  B=0.0000
 ├─ 0.8750  k = 7.5
 │× 0.9375  k = 8   ░░░░░░░░░░  A=0.9375  B=0.0000
 └─ 1.0000  k = 8.5

Pure sigma levels everywhere ( , equivalent to SigmaCoordinates):

SigmaPressureCoordinates(spectral_grid; transition = _ -> 1.0)

8-layer SigmaPressureCoordinates{Float32, Vector{Float64}}
    p_ref = 100000.0 Pa   ░ pressure │ sigma █
 ├─ 0.0000  k = 0.5
 │× 0.0625  k = 1   ██████████  A=0.0000  B=0.0625
 ├─ 0.1250  k = 1.5
 │× 0.1875  k = 2   ██████████  A=0.0000  B=0.1875
 ├─ 0.2500  k = 2.5
 │× 0.3125  k = 3   ██████████  A=0.0000  B=0.3125
 ├─ 0.3750  k = 3.5
 │× 0.4375  k = 4   ██████████  A=0.0000  B=0.4375
 ├─ 0.5000  k = 4.5
 │× 0.5625  k = 5   ██████████  A=0.0000  B=0.5625
 ├─ 0.6250  k = 5.5
 │× 0.6875  k = 6   ██████████  A=0.0000  B=0.6875
 ├─ 0.7500  k = 6.5
 │× 0.8125  k = 7   ██████████  A=0.0000  B=0.8125
 ├─ 0.8750  k = 7.5
 │× 0.9375  k = 8   ██████████  A=0.0000  B=0.9375
 └─ 1.0000  k = 8.5

A ready-made CubicSigmaPressureCoordinates uses a cubic smoothstep transition: pure pressure above (in the vertical sense) a given sigma value called pressure_only_above, typically in the upper atmosphere, pure sigma below sigma_only_below typically between the boundary layer and the mid-troposphere, and a smooth cubic interpolation in between. The bar on each full level shows the A/B split (█ = sigma fraction, ░ = pressure fraction):

CubicSigmaPressureCoordinates(spectral_grid)

8-layer SigmaPressureCoordinates{Float32, Vector{Float32}}
    p_ref = 100000.0 Pa   ░ pressure │ sigma █
 ├─ 0.0000  k = 0.5
 │× 0.0625  k = 1   ░░░░░░░░░░  A=0.0625  B=0.0000
 ├─ 0.1250  k = 1.5
 │× 0.1875  k = 2   ░░░░░░░░░░  A=0.1850  B=0.0025
 ├─ 0.2500  k = 2.5
 │× 0.3125  k = 3   █░░░░░░░░░  A=0.2715  B=0.0410
 ├─ 0.3750  k = 3.5
 │× 0.4375  k = 4   ████░░░░░░  A=0.2740  B=0.1635
 ├─ 0.5000  k = 4.5
 │× 0.5625  k = 5   ███████░░░  A=0.1892  B=0.3733
 ├─ 0.6250  k = 5.5
 │× 0.6875  k = 6   █████████░  A=0.0716  B=0.6159
 ├─ 0.7500  k = 6.5
 │× 0.8125  k = 7   ██████████  A=0.0074  B=0.8051
 ├─ 0.8750  k = 7.5
 │× 0.9375  k = 8   ██████████  A=0.0000  B=0.9375
 └─ 1.0000  k = 8.5

Custom thresholds can be set:

CubicSigmaPressureCoordinates(spectral_grid; pressure_only_above = 0.1, σ_only_below = 0.9)

8-layer SigmaPressureCoordinates{Float32, Vector{Float32}}
    p_ref = 100000.0 Pa   ░ pressure │ sigma █
 ├─ 0.0000  k = 0.5
 │× 0.0625  k = 1   ░░░░░░░░░░  A=0.0623  B=0.0002
 ├─ 0.1250  k = 1.5
 │× 0.1875  k = 2   █░░░░░░░░░  A=0.1758  B=0.0117
 ├─ 0.2500  k = 2.5
 │× 0.3125  k = 3   ██░░░░░░░░  A=0.2497  B=0.0628
 ├─ 0.3750  k = 3.5
 │× 0.4375  k = 4   ████░░░░░░  A=0.2613  B=0.1762
 ├─ 0.5000  k = 4.5
 │× 0.5625  k = 5   ██████░░░░  A=0.2104  B=0.3521
 ├─ 0.6250  k = 5.5
 │× 0.6875  k = 6   ████████░░  A=0.1200  B=0.5675
 ├─ 0.7500  k = 6.5
 │× 0.8125  k = 7   ██████████  A=0.0359  B=0.7766
 ├─ 0.8750  k = 7.5
 │× 0.9375  k = 8   ██████████  A=0.0013  B=0.9362
 └─ 1.0000  k = 8.5

Or custom sigma spacing too, to have, regardless sigma or sigma-pressure coordinates more layers near the surface and the top of the atmosphere:

(; σ_half) = FriersonSigmaCoordinates(spectral_grid)
CubicSigmaPressureCoordinates(spectral_grid;  σ_half)

8-layer SigmaPressureCoordinates{Float32, Vector{Float32}}
    p_ref = 100000.0 Pa   ░ pressure │ sigma █
 ├─ 0.0000  k = 0.5
 │× 0.0167  k = 1   ░░░░░░░░░░  A=0.0167  B=0.0000
 ├─ 0.0333  k = 1.5
 │× 0.0726  k = 2   ░░░░░░░░░░  A=0.0726  B=0.0000
 ├─ 0.1118  k = 2.5
 │× 0.1900  k = 3   ░░░░░░░░░░  A=0.1852  B=0.0048
 ├─ 0.2682  k = 3.5
 │× 0.3778  k = 4   ███░░░░░░░  A=0.2588  B=0.1190
 ├─ 0.4874  k = 4.5
 │× 0.5981  k = 5   ███████░░░  A=0.1516  B=0.4464
 ├─ 0.7088  k = 5.5
 │× 0.7905  k = 6   ██████████  A=0.0221  B=0.7684
 ├─ 0.8722  k = 6.5
 │× 0.9162  k = 7   ██████████  A=0.0000  B=0.9162
 ├─ 0.9603  k = 7.5
 │× 0.9801  k = 8   ██████████  A=0.0000  B=0.9801
 └─ 1.0000  k = 8.5

Using hybrid coordinates in a model

S = CubicSigmaPressureCoordinates(spectral_grid)
geometry = Geometry(spectral_grid; vertical_coordinates = S)
model = PrimitiveDryModel(spectral_grid; geometry)
simulation = initialize!(model)

Reference pressure and atmosphere

If the reference_pressure of the SigmaPressureCoordinates differs from the reference_pressure of the atmosphere (default 1e5 Pa), a warning is issued during initialize!. Make sure they match for physical consistency.

Coordinate functions

Three functions provide a coordinate-agnostic interface for evaluating the vertical coordinate at a given full level , working for both SigmaCoordinates and SigmaPressureCoordinates:

Function	Returns
SpeedyWeather.pressure(k, pₛ, coord)	Pressure [Pa] at full level
SpeedyWeather.pressure_thickness(k, pₛ, coord)	Pressure thickness [Pa] of layer
SpeedyWeather.sigma(k, coord)	Nominal sigma level (surface-pressure independent)

For SigmaCoordinates these reduce to , , and respectively. For SigmaPressureCoordinates they use the full   formula. Note that sigma always equals    regardless of the transition.