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

import numpy as np




def obukov_length(transf, inout_datastore, ddi, mapper):
    """Compute the Obukhov length L from sensible heat flux and friction velocity.

    Derives the Obukhov length ``obuk``, aerodynamic resistance ``aerr``
    (integrated Monin-Obukhov flux-profile relationships, stable or unstable),
    and convective velocity scale ``wsta`` from the near-surface potential
    temperature (``potts``), virtual heat flux (``potf``/``potfm``), friction
    velocity, and roughness length. See
    :doc:`/documentation/doc-models/chimere/diagmet` (section 10) 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
    vkarm = 0.4  # von Karman constant
    g = 9.81  # Gravity constant
    
    # Inputs
    potts = transf.diag_misc["potts"][:, np.newaxis]
    potf = transf.diag_misc["potf"]
    potfm = transf.diag_misc["potfm"]
    az0 = transf.diag_misc["az0"]
    ustar = inout_datastore["outputs"][("meteo", "usta")][ddi]["spec"].values
    hght = inout_datastore["outputs"][("meteo", "hght")][ddi]["spec"].values
    
    # Compute Obukov
    zchim2 = \
        inout_datastore["outputs"][("meteo", "alti")][ddi]["spec"][:, [0]].values * 0.5
    obuklen = -potts * ustar ** 3 / (vkarm * g * potf) # XXX potf < 0 => obuklen < 0
    
    zeta = zchim2 / obuklen
    zeta0 = az0 / obuklen # XXX problem: obuklen < 0. 
    zln = np.log(zchim2 / az0)
    eta0 = (1. - 15. * zeta0) ** 0.25
    eta = (1. - 15. * zeta) ** 0.25
    
    ra = np.where(
        obuklen >= 0,
        (zln + 4.7 * (zeta - zeta0)) / (vkarm * ustar),
        (zln + np.log((eta0 ** 2 + 1) * (eta0 + 1) ** 2
                   / (eta ** 2 + 1) / (eta + 1) ** 2)
         + 2 * (np.arctan(eta) - np.arctan(eta0))) / (vkarm * ustar)
    )
    wstar = np.where(obuklen >= 0, 0, (g * potfm * hght / potts) ** (1. / 3.))
    
    # Save outputs
    ref_dataarray = inout_datastore["inputs"][("meteo", "hght")][ddi]["spec"] * 0
    inout_datastore["outputs"][("meteo", "obuk")][ddi]["spec"] = ref_dataarray + obuklen
    inout_datastore["outputs"][("meteo", "aerr")][ddi]["spec"] = ref_dataarray + ra
    inout_datastore["outputs"][("meteo", "wsta")][ddi]["spec"] = ref_dataarray + wstar