Interpolating from hybrid to pressure levels

import earthkit.meteo.vertical.array as vertical

Hybrid levels

Hybrid levels are terrain following at the bottom and isobaric at the top of the atmosphere with a transition in between. This is the vertical coordinate system in the IFS model.

For details about the hybrid levels see IFS Documentation CY47R3 - Part IV Physical processes, Chapter 2, Section 2.2.1.. Also see the /how-tos/hybrid_levels.ipynb notebook for the related computations in earthkit-meteo.

Getting the data

First, we get some sample data containing 137 IFS model levels (hybrid full-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]
t = DATA.t  # temperature [K]
q = DATA.q  # specific humidity [kg/kg]

sp.shape, t.shape, q.shape

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

t[::40]

array([[197.06175232, 198.54808044],
       [230.29292297, 219.00386047],
       [234.56352234, 229.85160828],
       [264.58778381, 292.04481506]])

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,))

Using interpolate_hybrid_to_pressure_levels

We can use interpolate_hybrid_to_pressure_levels() for hybrid to pressure level interpolation. The example below interpolates temperature to pressure levels.

target_p = [85000.0, 50000.0]  # Pa

t_res = vertical.interpolate_hybrid_to_pressure_levels(
    t,  # data to interpolate
    target_p,
    sp,
    A,
    B,
    interpolation="linear",
)

for pv, rv in zip(target_p, t_res):
    print(pv, rv)

85000.0 [263.50982741 287.70299692]
50000.0 [238.00383748 259.50822691]

Using a subset of levels

It is possible to use only a subset of the model levels in the input data. This can significantly speed up the computations and reduce memory usage.

In this case the model level range in the input must be contiguous and include the bottom-most level. The following example only uses the lowest 50 model levels above the surface.

# using the 50 lowest model levels as input
# Please note we still need to use the full A and B arrays
t_res = vertical.interpolate_hybrid_to_pressure_levels(
    t[-50:],  # data to interpolate
    target_p,
    sp,
    A,
    B,
    interpolation="linear",
)

for pv, rv in zip(target_p, t_res):
    print(pv, rv)

85000.0 [263.50982741 287.70299692]
50000.0 [238.00383748 259.50822691]

Using aux levels

Since the lowest (full) model level is always above the surface the interpolation fails if the target pressure is between the surface and the lowest model level. We can check the pressure on the surface and the lowest model level in our data:

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

(array([ 95178.337944  , 102659.81019512]),
 array([[ 95065.55729093, 102538.16444369]]))

As a consequence, interpolation interpolation to 95100 Pa results in nan for the first point:

target_p = [95100.0]  # Pa

t_res = vertical.interpolate_hybrid_to_pressure_levels(
    t,  # data to interpolate
    target_p,
    sp,
    A,
    B,
    interpolation="linear",
)

for pv, rv in zip(target_p, t_res):
    print(pv, rv)

95100.0 [         nan 292.08454145]

To overcome this problem we can define “aux” levels with a prescribed input data. As a demonstration, we set the temperature values on the surface in all the input points with the aux_bottom_* kwargs.

target_p = [95100.0]  # Pa

t_res = vertical.interpolate_hybrid_to_pressure_levels(
    t,  # data to interpolate
    target_p,
    sp,
    A,
    B,
    aux_bottom_p=[95178.337944, 102659.81019512],
    aux_bottom_data=[270.0, 293.0],
    interpolation="linear",
)

for pv, rv in zip(target_p, t_res):
    print(pv, rv)

95100.0 [269.21926951 292.08454145]

Using interpolate_monotonic

We can also use interpolate_monotonic() to carry out the interpolation. This is a generic method and we need an extra step to compute the pressure on (full) model levels using pressure_on_hybrid_levels(). This pressure data then can be used to interpolate multiple input data arrays.

First, compute the pressure all the model levels in the input data.

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

(137, 2)

Next, interpolate the temperature to the target pressure levels.

target_p = [85000.0, 50000.0]  # Pa

t_res = vertical.interpolate_monotonic(t, p, target_p, interpolation="linear")

for pv, rv in zip(target_p, t_res):
    print(pv, rv)

85000.0 [263.50982741 287.70299692]
50000.0 [238.00383748 259.50822691]

Repeat the same computation for specific humidity.

q_res = vertical.interpolate_monotonic(q, p, target_p, interpolation="linear")

for pv, rv in zip(target_p, q_res):
    print(pv, rv)

85000.0 [0.00119029 0.00621813]
50000.0 [0.00010948 0.00046553]

Using a subset of levels

We can also use a subset of levels. This can significantly speed up the computations and reduce memory usage.

In this case the model level range in the input data must be contiguous and include the bottom-most level. The example below only uses the lowest 50 model levels above the surface.

# compute pressure on the 50 lowest model levels
# Please note we still need to use the full A and B arrays.
levels = list(range(138 - 50, 138))  # 97-137
p_sub = vertical.pressure_on_hybrid_levels(sp, A=A, B=B, levels=levels, output="full")


# interpolate to pressure levels
t_res = vertical.interpolate_monotonic(t[-50:], p_sub, target_p, interpolation="linear")

for pv, rv in zip(target_p, t_res):
    print(pv, rv)

85000.0 [263.50982741 287.70299692]
50000.0 [238.00383748 259.50822691]

Using aux levels

As described above, the lowest (full) model level is always above the surface so the interpolation does not work if the target pressure is between the surface and the lowest model level.

# this pressure is larger than the pressure on the lowest model level in point 1
target_p = [95100.0]  # Pa

t_res = vertical.interpolate_monotonic(t, p, target_p, interpolation="linear")

for pv, rv in zip(target_p, t_res):
    print(pv, rv)

95100.0 [         nan 292.08454145]

To overcome this problem we can define “aux” levels with a prescribed input data. As a demonstration, we set the temperature values on the surface in all the input points with the aux_max_level_* kwargs.

target_p = [95100.0]  # Pa

t_res = vertical.interpolate_monotonic(
    t,
    p,
    target_p,
    aux_max_level_coord=[95178.337944, 102659.81019512],
    aux_max_level_data=[270.0, 293.0],
    interpolation="linear",
)

for pv, rv in zip(target_p, t_res):
    print(pv, rv)

95100.0 [269.21926951 292.08454145]