Vertical interpolation onto pressure layers

SpeedyWeather.jl integrates on terrain-following vertical coordinates, either Sigma coordinates or Hybrid sigma-pressure coordinates, see Vertical coordinates. Model layers therefore sit at a different pressure in every column, bending around the orography. For analysis and for comparison with observations or reanalyses one usually wants the data on fixed pressure layers instead, e.g. temperature at 850 hPa everywhere. This page describes how to interpolate from the model's vertical layers onto pressure layers.

Interpolating a field onto pressure layers

SpeedyWeather.interpolate_pressure_layers! interpolates a field on model layers onto pressure layers, writing into an output field that you allocate yourself

SpeedyWeather.interpolate_pressure_layers!(
    out_field,          # OUTPUT: (horizontal, npressure)
    in_field,           # INPUT: (horizontal, nlayers) on model layers
    surface_pressure,   # INPUT: (horizontal,) surface pressure [Pa]
    p,                  # pressure layers [Pa]
    coordinates,        # model.geometry.vertical_coordinates
    interpolation,      # optional, LinearInLogPressure() by default
    extrapolation,      # optional, ConstantExtrapolation() by default
)

The pressure of every model layer in every column follows from the surface pressure and the vertical coordinates, so the same function covers sigma and hybrid sigma-pressure coordinates. The interpolation is column-local and therefore launched as a kernel over horizontal grid points and pressure layers, running on CPU and GPU. Consequently the pressure layers p have to be on the same architecture (see GPU and Architectures) as the fields, and ideally of the same number format. They do not have to be sorted.

The number of pressure layers is independent of the number of model layers, it is the second dimension of out_field that decides. Note that in_field must not have a time step dimension, so for prognostic variables select the step first, see Step dimension.

Interpolation methods

Interpolating a variable between two model layers at pressures    uses the weight of the lower layer

with the two available choices for being linear in pressure

or linear in the logarithm of pressure

SpeedyWeather.LinearInLogPressure is the default as most variables vary more linearly with than with .

Extrapolation beyond the model layers

A pressure layer can also lie outside the range spanned by the full model layers: above the top-most layer, or below the bottom-most layer. The latter splits into two cases, because the lowest full model layer is not the surface: a layer below the lowest full layer can still be above ground (p < pₛ), or genuinely below ground (p > pₛ). The following are available

• ConstantExtrapolation (default) holds the outer-most model layer constant.

• DryAdiabaticExtrapolation descends dry-adiabatically below the lowest model layer,  , which is what you want for a temperature at, say, 1000 hPa below a lowest model layer that sits above it. It is the same adiabat that the mean sea-level pressure output uses.

• SubsurfaceMask masks everything below the surface with a missing_value (NaN by default) and uses another extrapolation between the lowest model layer and the surface.

Above the model top every method currently holds the top-most layer constant.

Example

using SpeedyWeather
spectral_grid = SpectralGrid(truncation = 31, nlayers = 8)
model = PrimitiveDryModel(spectral_grid)
simulation = initialize!(model)
run!(simulation, period = Hour(6))

Take the temperature on model layers (selecting the time step to interpolate) and the surface pressure

(; variables) = simulation
temp = SpeedyWeather.get_prognostic_step(variables.grid.temperature, model.time_stepping, model.output)
pₛ = variables.parameterizations.surface_pressure
summary(temp)

"3168×8, 48-ring (XYZ) OctahedralGaussianField{Float32, 2} as Array on CPU"

allocate an output field for the pressure layers of interest, and interpolate

p = spectral_grid.NF[850e2, 500e2, 200e2]       # pressure layers in Pa
temp_p = zeros(spectral_grid.NF, spectral_grid.grid, length(p))

SpeedyWeather.interpolate_pressure_layers!(temp_p, temp, pₛ, p, model.geometry.vertical_coordinates)
[sum(temp_p[:, k]) / length(pₛ) for k in eachindex(p)]   # mean temperature [K] per layer

3-element Vector{Float32}:
 271.41107
 238.04094
 201.59813

A 1000 hPa layer is below ground over much of the globe. Extrapolating dry-adiabatically down to it while masking the points that are genuinely below the surface

p = spectral_grid.NF[1000e2]
temp_1000 = zeros(spectral_grid.NF, spectral_grid.grid, length(p))

SpeedyWeather.interpolate_pressure_layers!(
    temp_1000, temp, pₛ, p, model.geometry.vertical_coordinates,
    SpeedyWeather.LinearInLogPressure(),
    SpeedyWeather.SubsurfaceMask(
        above_surface = SpeedyWeather.DryAdiabaticExtrapolation(model.atmosphere.κ),
    ),
)
count(isnan, temp_1000), length(pₛ)     # masked points of total

(1869, 3168)

Output on pressure layers

You normally don't need to call any of this yourself: every output writer takes a layers keyword argument to write its 3D variables on pressure layers instead of model layers

output = NetCDFOutput(spectral_grid, PrimitiveWet, layers = PressureLayers(spectral_grid, [850, 500, 200] .* 100))

See Output layers for the details, including how a variable chooses its extrapolation below the lowest model layer.