Source code for pycif.plugins.transforms.complex.diagmet.utils.boundary_layer_height

import numpy as np




def boundary_layer_height(transf, inout_datastore, ddi, mapper):
    """Interpolate the virtual potential temperature to a fixed 25 m reference height.

    .. note::

        Despite its name (inherited from the historical ``diagmet.f90``), this
        step does not compute the PBL height: CHIMERE reuses the ECMWF boundary
        layer height (``hght``) directly, passed through unchanged elsewhere in
        the transform. It only linearly interpolates the virtual potential
        temperature ``po`` (computed by
        :func:`~pycif.plugins.transforms.complex.diagmet.utils.friction_velocity.friction_velocity`)
        to :math:`z_\\mathrm{therm}=25` m, storing the result as ``potts`` for
        :func:`~pycif.plugins.transforms.complex.diagmet.utils.obukov_length.obukov_length`.

    See :doc:`/documentation/doc-models/chimere/diagmet` (section 7) for the
    full derivation.

    Args:
        transf (Plugin): diagmet transform instance.
        inout_datastore (dict): mutable datastore.
        ddi (datetime): current sub-simulation date.
        mapper (dict): transform mapper.
    """
    
    # Parameters
    ztherm = 25.0  # Height of thermals start
    
    # Variables from other functions
    po = transf.diag_misc["po"]
    alti = transf.diag_misc["alti"]
    
    ntime, nlev, nlat, nlon = alti.shape
    
    # Compute
    ind_therm = np.maximum(0, np.argmin(alti <= ztherm, axis=1) - 1)
    t, i, j = np.meshgrid(np.arange(ntime), np.arange(nlat), np.arange(nlon),
                          indexing="ij")
    po0 = po[t, ind_therm, i, j]
    po1 = po[t, ind_therm + 1, i, j]
    alt0 = alti[t, ind_therm, i, j]
    alt1 = alti[t, ind_therm + 1, i, j]
    
    potts = po0 + (po1 - po0) * (ztherm - alt0) / (alt1 - alt0)
    
    transf.diag_misc["potts"] = potts