Source code for meteodatalab.operators.hzerocl

"""Algorithm for computation of height of zero degree isotherm."""

# Standard library
from typing import cast

# Third-party
import numpy as np
import xarray as xr

# Local
from .. import metadata
from .destagger import destagger


[docs] def fhzerocl( t: xr.DataArray, hhl: xr.DataArray, extrapolate: bool = False ) -> xr.DataArray: """Height of the zero deg C isotherm in m amsl. The algorithm searches from the top of the atmosphere downwards. When extrapolation is enabled, the search may extend past the lowest layers of the model by means of a linear model. The resulting height may be below the surface of the earth. Parameters ---------- t : xr.DataArray Air temperature in K. hhl : xr.DataArray Heights of the interfaces between vertical layers in m amsl. extrapolate : bool, optional Allow the extrapolation of the search below the lowest model layer. Defaults to False. Returns ------- xr.DataArray Height of the zero deg C isotherm in m amsl. """ # Physical constants t0 = 273.15 # Heights of layer mid surfaces (where t is defined) hfl = destagger(hhl, "z") tkm1 = t.shift(z=1) # 3d field with values of height for those levels where temperature # is > 0 and it was < 0 on the level below. Otherwise values are NaN. height0 = hfl.where((t >= t0) & (tkm1 < t0)) # The previous condition can be satisfied on multiple levels. # Take the k indices of the maximum height value where the condition is satisfied maxind = cast(xr.DataArray, height0.fillna(-1).argmax(dim="z")) if extrapolate: # The full column is below freezing below_ground = (t < t0).all("z") # Temperature is increasing with decreasing altitude positive_dt = t[{"z": -1}] - t[{"z": -2}] > 1e-10 # Allow t0 to be outside of the [height1, height2] range cond = np.logical_and(below_ground, positive_dt) maxind = maxind.where(~cond, -1) # Compute the 2D fields with height values where T is > 0 and < 0 on level below height2 = hfl[{"z": maxind}] # Compute the 2D fields with height values where T is < 0 and > 0 on level above height1 = hfl[{"z": maxind - 1}] # The height level where T == t0 must be between [height1, height2] t2 = t[{"z": maxind}] t1 = tkm1[{"z": maxind}] hzerocl = height1 + (height2 - height1) * (t0 - t1) / (t2 - t1) return xr.DataArray( data=hzerocl.where(hzerocl > 0), attrs=metadata.override(t.message, shortName="HZEROCL"), )