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

import numpy as np
from .mean_z0_shf_extra_urban_temp import mean_z0_shf_extra_urban_temp
from .sv_heat_flux import sv_heat_flux




def friction_velocity(transf, inout_datastore, ddi, mapper):
    """Compute friction velocity u* using the Louis (1982) stability function.

    Only run when ``usta: recompute`` (the default); otherwise :math:`u_*` is
    read directly from ECMWF. Builds a bulk Richardson number between the
    surface and 10 m from virtual potential temperature and wind speed, applies
    the Louis (1982) stability correction to the neutral drag coefficient, and
    combines the mechanical and convective (``wstar0``, see
    :func:`~pycif.plugins.transforms.complex.diagmet.utils.sv_heat_flux.sv_heat_flux`)
    contributions. See :doc:`/documentation/doc-models/chimere/diagmet`
    (section 6) 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.
    """
    # Wind speed at 10m
    w10m = np.sqrt(
        inout_datastore["inputs"][("meteo", "u10m")][ddi]["spec"] ** 2
        + inout_datastore["inputs"][("meteo", "v10m")][ddi]["spec"] ** 2
    ).values
    
    # Inputs
    pres = transf.diag_misc["pres"]
    temp = transf.diag_misc["temp"]
    sphu = transf.diag_misc["sphu"]
    water = transf.diag_misc["water"]
    
    # Correction by urban heat
    awf = transf.diag_misc["awf"]
    az0 = transf.diag_misc["az0"]
    auf = transf.diag_misc["auf"]
    pm = transf.diag_misc["pm"]
    
    w10s = w10m * awf
    
    # SV heat flux
    wstar0 = transf.diag_misc["wstar0"]
    
    # Fixed parameters
    woff = 0.5  # Wind offset to smooth Richardson numbers (m / s)
    g = 9.81  # Gravity constant
    xkappa = 0.2857  # Kappa
    vkarm = 0.4  # von Karman constant
    p0 = 1e5  # Reference pressure(Pa)
    
    # Compute auxiliary parameters
    po = (1.0 + 0.61 * sphu - water) * temp * (p0 / pres) ** xkappa
    th = temp * (p0 / pres) ** xkappa
    transf.diag_misc["po"] = po
    
    
    if transf.usta == 'recompute':
        # Compute ustar
        zustar = 10.
        vustar = np.maximum(w10s, woff)
    
        dtheta = po[:, [1]] - po[:, [0]]
        rich = dtheta * g * zustar / (th[:, [0]] * vustar ** 2)
        cdnm = vkarm / np.log(zustar / az0)
        cdn2 = cdnm * cdnm
        facm = 75.0 * cdn2 * np.sqrt(zustar / az0)

        fm = np.where(
          rich < 0.,
          1.0 - 10.0 * rich / (1.0 + facm * np.sqrt(-rich)),
          1.0 / (1.0 + 10.0 * rich / np.sqrt(1. + 5. * rich))
         )
    
        cd = cdnm * np.sqrt(fm)
        ustar = cd * np.sqrt(vustar ** 2 + (1.2 * wstar0) ** 2)
    
        inout_datastore["outputs"][("meteo", "usta")][ddi]["spec"] = \
          0 * inout_datastore["outputs"][("meteo", "hght")][ddi]["spec"] + ustar