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

import numpy as np




def sv_heat_flux(transf, inout_datastore, ddi, mapper):
    """Compute sensible and virtual heat fluxes from surface temperature and humidity fields.

    Derives the surface virtual potential-temperature flux ``potf`` from the
    ECMWF accumulated sensible/latent heat fluxes, and a convective velocity
    scale ``wstar0`` from ``potf`` using a fixed convective height scale of
    1500 m. Both are stored in ``transf.diag_misc`` for use by
    :func:`~pycif.plugins.transforms.complex.diagmet.utils.friction_velocity.friction_velocity`
    and :func:`~pycif.plugins.transforms.complex.diagmet.utils.obukov_length.obukov_length`.
    See :doc:`/documentation/doc-models/chimere/diagmet` (section 5) 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
    p0 = 1e5  # Reference pressure(Pa)
    xkappa = 0.2857  # Kappa
    R = 287.04  # R constant
    Cp = 1005.  # Cp
    Lv = 2.45e6  # Lv
    Rv = 461.5  # Rv
    rcp = R / Cp  # R / Cp
    rl = R / Lv  # R / Lv
    g = 9.81  # Gravity constant
    chs0 = 1500. # Convective height scale(m)
    
    # Inputs
    slhf = inout_datastore["inputs"][("meteo", "slhf")][ddi]["spec"].values
    sshf = inout_datastore["inputs"][("meteo", "sshf")][ddi]["spec"].values
    temp2 = inout_datastore["inputs"][("meteo", "tem2")][ddi]["spec"].values
    temp = transf.diag_misc["temp"]
    sphu = transf.diag_misc["sphu"]
    pres = transf.diag_misc["pres"]
    auf = transf.diag_misc["auf"]
    
    # Compute sensible and virtual heat fluxes in K.m/s
    th = temp * (p0 / pres) ** xkappa
    
    heat = rcp * temp2 * sshf / pres[:, [0]]  # XXX humf <0 => potf < 0.
    humf = rl * temp2 * slhf / pres[:, [0]]  # XXX humf <0 => potf < 0.
    potf = heat * (1 + 0.61 * sphu[:, [0]]) \
           + humf * 0.61 * th[:, [0]] + rcp * temp2 * auf / pres[:, [0]]

    potfm = np.maximum(potf, 1.e-6)
    wstar0 = (g * potfm * chs0 / temp2) ** (1. / 3.)
    
    # Save variables for later
    transf.diag_misc["wstar0"] = wstar0
    transf.diag_misc["potf"] = potf   #notzero(potf,1.e-30) ?
    transf.diag_misc["potfm"] = potfm
    transf.diag_misc["th"] = th