Computing values on hybrid levels

import earthkit.meteo.vertical.array as vertical

This notebook gives an overview about computing various quantities on hybrid levels.

Hybrid levels

Hybrid levels divide the atmosphere into layers. These layers are defined by the pressures at the interfaces between them, which are the half-levels. The half-levels are defined by a set A and B coefficients in such a way that at the top of the atmosphere the first half level pressure is a constant, while at the bottom the last one is the surface pressure. The full-level pressure associated with each model level is defined as the middle of the layer. As a consequence, the bottom-most (full) model level is always above the surface. Please note that by convention the model level numbering starts at 1 at the top of the atmosphere and increases towards the surface.

For more details about the hybrid levels see IFS Documentation CY47R3 - Part IV Physical processes, Chapter 2, Section 2.2.1.

Getting the data

First, we get some sample data containing 137 IFS (full) model levels for two points. This data is in ascending order with regards to the model level number. So the first level is model level 1 (i.e. the top), while the last one is model level 137 (i.e. the bottom).

from earthkit.meteo.utils.sample import get_sample

DATA = get_sample("vertical_hybrid_data")
sp = DATA.p_surf  # surface pressure [Pa]
zs = DATA.z_surf  # surface geopotential [m2/s2]
t = DATA.t  # temperature [K]
q = DATA.q  # specific humidity [kg/kg]

sp.shape, zs.shape, t.shape, q.shape

((2,), (2,), (137, 2), (137, 2))

Next, get the hybrid level definitions using hybrid_level_parameters().

A, B = vertical.hybrid_level_parameters(137, model="ifs")
A.shape, B.shape

((138,), (138,))

The A and B arrays contain coefficients for each hybrid half-level in ascending model level number order.

Computing pressure

We can compute the pressure on all the full-levels from the surface pressure with pressure_on_hybrid_levels().

p_full = vertical.pressure_on_hybrid_levels(sp, A=B, B=B, output="full")
p_full.shape

(137, 2)

The output is ascending order with regards to the model level number, so levels go from the top of the atmosphere towards the surface.

# the lowest 3 levels above the surface
p_full[-3:]

array([[ 94560.03063869, 101992.82288719],
       [ 94828.72712554, 102282.63998019],
       [ 95066.55610599, 102539.16325875]])

By using the output kwarg other pressure related quantities can be computed.

p_full, p_half, alpha, delta = vertical.pressure_on_hybrid_levels(
    sp, A=A, B=B, output=("full", "half", "alpha", "delta")
)

p_full.shape, p_half.shape, alpha.shape, delta.shape

((137, 2), (138, 2), (137, 2), (137, 2))

It is possible to use only part of the levels via the levels option. The level indices start at 1 and go up to 137 for the current data. The example below only computes the pressure on the lowest 3 levels.

p_full = vertical.pressure_on_hybrid_levels(sp, A=A, B=B, levels=[135, 136, 137], output=("full"))
p_full

array([[ 94572.33402093, 102005.12626943],
       [ 94829.60971572, 102283.52257037],
       [ 95065.55729093, 102538.16444369]])

The level kwarg respects the specified ordering. In the example below the output goes upwards from the lowest level.

p_full = vertical.pressure_on_hybrid_levels(sp, A=A, B=B, levels=[137, 136, 135], output=("full"))
p_full

array([[ 95065.55729093, 102538.16444369],
       [ 94829.60971572, 102283.52257037],
       [ 94572.33402093, 102005.12626943]])

Computing geopotential thickness

relative_geopotential_thickness_on_hybrid_levels() compute the geopotential thickness between the surface and the full hybrid levels. It requires temperature and specific humidity profiles.

z_thickness = vertical.relative_geopotential_thickness_on_hybrid_levels(t, q, sp, A, B)

Note that t and q must contain the same model levels in ascending order with respect to the model level number. The model level range must be contiguous and must include the bottom-most level, but not all the levels must be present. E.g. for our data we can compute the relative geopotential thickness for the lowest 3 levels as shown below. Please note that even in this case we use the full A and B arrays.

z_thickness = vertical.relative_geopotential_thickness_on_hybrid_levels(t[-3:], q[-3:], sp, A, B)
z_thickness

array([[493.46642874, 540.17339902],
       [283.62077829, 310.1849329 ],
       [ 91.62752892, 100.19125537]])

Computing geopotential

We can compute the geopotential on full hybrid levels with geopotential_on_hybrid_levels().

z = vertical.geopotential_on_hybrid_levels(t, q, zs, sp, A, B)
z.shape, z[-3:]

((137, 2),
 array([[5774.59167288, -289.70135684],
        [5564.74602243, -519.68982296],
        [5372.75277306, -729.68350049]]))

Just like for the relative geopotential thickness a subset of levels can be used. E.g. to compute the geopotential for the lowest 3 levels we can write:

z = vertical.geopotential_on_hybrid_levels(t[-3:], q[-3:], zs, sp, A, B)
z.shape, z

((3, 2),
 array([[5774.59167288, -289.70135684],
        [5564.74602243, -519.68982296],
        [5372.75277306, -729.68350049]]))

Computing height

The height on full hybrid levels can be computed with height_on_hybrid_levels(). We can specify the height type via h_type and the reference level via h_reference. The following options are available:

• h_type: “geometric” or “geopotential”

• h_reference: “ground” or “sea”

h = vertical.height_on_hybrid_levels(t, q, zs, sp, A, B, h_type="geometric", h_reference="ground")
h.shape, h[-3:]

((137, 2),
 array([[50.32847687, 55.08137028],
        [28.92629165, 31.62937732],
        [ 9.34500108, 10.21640974]]))

h = vertical.height_on_hybrid_levels(t, q, zs, sp, A, B, h_type="geometric", h_reference="sea")
h.shape, h[-3:]

((137, 2),
 array([[588.89890268, -29.54118008],
        [567.49671746, -52.99317304],
        [547.91542689, -74.40614062]]))

Just like for the previous methods a subset of levels can be used. E.g. to compute the geometric height above ground for the lowest 3 levels we can write:

h = vertical.height_on_hybrid_levels(t[-3:], q[-3:], 0, sp, A, B, h_type="geometric", h_reference="ground")
h.shape, h

((3, 2),
 array([[50.31996922, 55.0828335 ],
        [28.92140188, 31.63021754],
        [ 9.34342138, 10.21668113]]))