flux_subs.py

import numpy as np
from util_subs import (kappa, visc_air)

# ---------------------------------------------------------------------

def cdn_calc(u10n, usr, Ta, grav, meth):
    """
    Calculates neutral drag coefficient

    Parameters
    ----------
    u10n : float
        neutral 10m wind speed [m/s]
    usr : float
        friction velocity      [m/s]
    Ta   : float
        air temperature        [K]
    grav : float
        gravity               [m/s^2]
    meth : str

    Returns
    -------
    cdn : float
    zo  : float
    """
    cdn = np.zeros(Ta.shape)*np.nan
    if (meth == "S80"): # eq. 14 Smith 1980
        cdn = (0.61+0.063*u10n)*0.001
    elif (meth == "LP82"):
        #  Large & Pond 1981 u10n <11m/s & eq. 21 Large & Pond 1982
        cdn = np.where(u10n < 11, 1.2*0.001, (0.49+0.065*u10n)*0.001)
    elif (meth == "S88" or meth == "UA" or meth == "ecmwf" or meth == "C30" or
          meth == "C35" or meth == "Beljaars"): #  or meth == "C40"
        cdn = cdn_from_roughness(u10n, usr, Ta, grav, meth)
    elif (meth == "YT96"):
        # convert usr in eq. 21 to cdn to expand for low wind speeds
        cdn = np.power((0.10038+u10n*2.17e-3+np.power(u10n, 2)*2.78e-3 -
                        np.power(u10n, 3)*4.4e-5)/u10n, 2)
    elif (meth == "NCAR"): # eq. 11 Large and Yeager 2009
        cdn = np.where(u10n > 0.5, (0.142+2.7/u10n+u10n/13.09 -
                                    3.14807e-10*np.power(u10n, 6))*1e-3,
                       (0.142+2.7/0.5+0.5/13.09 -
                        3.14807e-10*np.power(0.5, 6))*1e-3)
        cdn = np.where(u10n > 33, 2.34e-3, np.copy(cdn))
        cdn = np.maximum(np.copy(cdn), 0.1e-3)
    else:
        raise ValueError("unknown method cdn: "+meth)

    zo = 10/np.exp(kappa/np.sqrt(cdn))
    
    return cdn, zo
# ---------------------------------------------------------------------


def cdn_from_roughness(u10n, usr, Ta, grav, meth):
    """
    Calculates neutral drag coefficient from roughness length

    Parameters
    ----------
    u10n : float
        neutral 10m wind speed [m/s]
    usr : float
        friction velocity      [m/s]
    Ta   : float
        air temperature        [K]
    grav : float                [m/s]
        gravity
    meth : str

    Returns
    -------
    cdn : float
    """
    #cdn = (0.61+0.063*u10n)*0.001
    zo, zc, zs = np.zeros(Ta.shape), np.zeros(Ta.shape), np.zeros(Ta.shape)
    for it in range(5):
        if (meth == "S88"):
            # Charnock roughness length (eq. 4 in Smith 88)
            zc = 0.011*np.power(usr, 2)/grav
            #  smooth surface roughness length (eq. 6 in Smith 88)
            zs = 0.11*visc_air(Ta)/usr
            zo = zc + zs  #  eq. 7 & 8 in Smith 88
        elif (meth == "UA"):
            # valid for 0<u<18m/s # Zeng et al. 1998 (24)
            zo = 0.013*np.power(usr, 2)/grav+0.11*visc_air(Ta)/usr
        elif (meth == "C30"): # eq. 25 Fairall et al. 1996a
            a = 0.011*np.ones(Ta.shape)
            a = np.where(u10n > 10, 0.011+(u10n-10)*(0.018-0.011)/(18-10),
                         np.where(u10n > 18, 0.018, a))
            zo = a*np.power(usr, 2)/grav+0.11*visc_air(Ta)/usr
        elif (meth == "C35"): # eq.6-11 Edson et al. (2013)
            zo = (0.11*visc_air(Ta)/usr +
                  np.minimum(0.0017*19-0.0050, 0.0017*u10n-0.0050) *
                  np.power(usr, 2)/grav)
        elif ((meth == "ecmwf" or meth == "Beljaars")):
            # eq. (3.26) p.38 over sea IFS Documentation cy46r1
            zo = 0.018*np.power(usr, 2)/grav+0.11*visc_air(Ta)/usr
        else:
            raise ValueError("unknown method for cdn_from_roughness "+meth)
            
        cdn = np.power(kappa/np.log(10/zo), 2)
    return cdn
# ---------------------------------------------------------------------


def cd_calc(cdn, hin, hout, psim):
    """
    Calculates drag coefficient at reference height

    Parameters
    ----------
    cdn : float
        neutral drag coefficient
    hin : float
        wind speed height       [m]
    hout : float
        reference height        [m]
    psim : float
        momentum stability function

    Returns
    -------
    cd : float
    """
    cd = (cdn/np.power(1+(np.sqrt(cdn)*(np.log(hin/hout)-psim))/kappa, 2))
    return cd
# ---------------------------------------------------------------------


def ctqn_calc(corq, zol, cdn, usr, zo, Ta, meth):
    """
    Calculates neutral heat and moisture exchange coefficients

    Parameters
    ----------
    corq : flag to select
           "ct" or "cq"
    zol  : float
        height over MO length
    cdn  : float
        neutral drag coefficient
    usr : float
        friction velocity      [m/s]
    zo   : float
        surface roughness       [m]
    Ta   : float
        air temperature         [K]
    meth : str

    Returns
    -------
    ctqn : float
        neutral heat exchange coefficient
    zotq : float
        roughness length for t or q
    """
    if (meth == "S80" or meth == "S88" or meth == "YT96"):
        cqn = np.ones(Ta.shape)*1.20*0.001  # from S88
        ctn = np.ones(Ta.shape)*1.00*0.001
        zot = 10/(np.exp(np.power(kappa, 2) / (ctn*np.log(10/zo))))
        zoq = 10/(np.exp(np.power(kappa, 2) / (cqn*np.log(10/zo))))
    elif (meth == "LP82"):
        cqn = np.where((zol <= 0), 1.15*0.001, 1*0.001)
        ctn = np.where((zol <= 0), 1.13*0.001, 0.66*0.001)
        zot = 10/(np.exp(np.power(kappa, 2) / (ctn*np.log(10/zo))))
        zoq = 10/(np.exp(np.power(kappa, 2) / (cqn*np.log(10/zo))))
    elif (meth == "NCAR"):
        cqn = np.maximum(34.6*0.001*np.sqrt(cdn), 0.1e-3)
        ctn = np.maximum(np.where(zol <= 0, 32.7*0.001*np.sqrt(cdn),
                                  18*0.001*np.sqrt(cdn)), 0.1e-3)
        zot = 10/(np.exp(np.power(kappa, 2) / (ctn*np.log(10/zo))))
        zoq = 10/(np.exp(np.power(kappa, 2) / (cqn*np.log(10/zo))))
    elif (meth == "UA"):
        # Zeng et al. 1998 (25)
        rr = usr*zo/visc_air(Ta)
        zoq = zo/np.exp(2.67*np.power(rr, 1/4)-2.57)
        zot = np.copy(zoq)
        cqn = np.power(kappa, 2)/(np.log(10/zo)*np.log(10/zoq))
        ctn = np.power(kappa, 2)/(np.log(10/zo)*np.log(10/zoq))
    elif (meth == "C30"):
        rr = zo*usr/visc_air(Ta)
        zoq = np.minimum(5e-5/np.power(rr, 0.6), 1.15e-4)  # moisture roughness
        zot = np.copy(zoq)  # temperature roughness
        cqn = np.power(kappa, 2)/np.log(10/zo)/np.log(10/zoq)
        ctn = np.power(kappa, 2)/np.log(10/zo)/np.log(10/zot)
    elif (meth == "C35"):
        rr = zo*usr/visc_air(Ta)
        zoq = np.minimum(5.8e-5/np.power(rr, 0.72), 1.6e-4) # moisture roughness
        zot = np.copy(zoq)  # temperature roughness
        cqn = np.power(kappa, 2)/np.log(10/zo)/np.log(10/zoq)
        ctn = np.power(kappa, 2)/np.log(10/zo)/np.log(10/zot)
    elif (meth == "ecmwf" or meth == "Beljaars"):
        # eq. (3.26) p.38 over sea IFS Documentation cy46r1
        zot = 0.40*visc_air(Ta)/usr
        zoq = 0.62*visc_air(Ta)/usr
        cqn = np.power(kappa, 2)/np.log(10/zo)/np.log(10/zoq)
        ctn = np.power(kappa, 2)/np.log(10/zo)/np.log(10/zot)
    else:
        raise ValueError("unknown method ctqn: "+meth)

    if corq == "ct":
        ctqn = ctn
        zotq = zot
    elif corq == "cq":
        ctqn = cqn
        zotq = zoq
    else:
        raise ValueError("unknown flag - should be ct or cq: "+corq)
        
    return ctqn, zotq
# ---------------------------------------------------------------------
# ---------------------------------------------------------------------

def ctq_calc(cdn, cd, ctqn, hin, hout, psitq):
    """
    Calculates heat and moisture exchange coefficients at reference height

    Parameters
    ----------
    cdn : float
        neutral drag coefficient
    cd  : float
        drag coefficient at reference height
    ctqn : float
        neutral heat or moisture exchange coefficient
    hin : float
        original temperature or humidity sensor height [m]
    hout : float
        reference height                   [m]
    psitq : float
        heat or moisture stability function

    Returns
    -------
    ctq : float
       heat or moisture exchange coefficient
    """
    ctq = (ctqn*np.sqrt(cd/cdn) /
          (1+ctqn*((np.log(hin/hout)-psitq)/(kappa*np.sqrt(cdn)))))

    return ctq
# ---------------------------------------------------------------------

def get_stabco(meth):
    """
    Gives the coefficients \\alpha, \\beta, \\gamma for stability functions

    Parameters
    ----------
    meth : str

    Returns
    -------
    coeffs : float
    """
    alpha, beta, gamma = 0, 0, 0
    if (meth == "S80" or meth == "S88" or meth == "NCAR" or
        meth == "UA" or meth == "ecmwf" or meth == "C30" or
        meth == "C35" or meth == "Beljaars"):
        alpha, beta, gamma = 16, 0.25, 5  # Smith 1980, from Dyer (1974)
    elif (meth == "LP82"):
        alpha, beta, gamma = 16, 0.25, 7
    elif (meth == "YT96"):
        alpha, beta, gamma = 20, 0.25, 5
    else:
        raise ValueError("unknown method stabco: "+meth)
    coeffs = np.zeros(3)
    coeffs[0] = alpha
    coeffs[1] = beta
    coeffs[2] = gamma
    return coeffs
# ---------------------------------------------------------------------


def psim_calc(zol, meth):
    """
    Calculates momentum stability function

    Parameters
    ----------
    zol : float
        height over MO length
    meth : str

    Returns
    -------
    psim : float
    """
    if (meth == "ecmwf"):
        psim = psim_ecmwf(zol)
    elif (meth == "C30" or meth == "C35"):
        psim = psiu_26(zol, meth)
    elif (meth == "Beljaars"): # Beljaars (1997) eq. 16, 17
        psim = np.where(zol < 0, psim_conv(zol, meth), psi_Bel(zol))
    else:
        psim = np.where(zol < 0, psim_conv(zol, meth),
                        psim_stab(zol, meth))
    return psim
# ---------------------------------------------------------------------


def psit_calc(zol, meth):
    """
    Calculates heat stability function

    Parameters
    ----------
    zol : float
        height over MO length
    meth : str
        parameterisation method

    Returns
    -------
    psit : float
    """
    if (meth == "ecmwf"):
        psit = np.where(zol < 0, psi_conv(zol, meth),
                        psi_ecmwf(zol))
    elif (meth == "C30" or meth == "C35"):
        psit = psit_26(zol)
    elif (meth == "Beljaars"): # Beljaars (1997) eq. 16, 17
        psit = np.where(zol < 0, psi_conv(zol, meth), psi_Bel(zol))
    else:
        psit = np.where(zol < 0, psi_conv(zol, meth),
                        psi_stab(zol, meth))
    return psit
# ---------------------------------------------------------------------


def psi_Bel(zol):
    """
    Calculates momentum/heat stability function

    Parameters
    ----------
    zol : float
        height over MO length
    meth : str
        parameterisation method

    Returns
    -------
    psit : float
    """
    a, b, c, d = 0.7, 0.75, 5, 0.35
    psi = -(a*zol+b*(zol-c/d)*np.exp(-d*zol)+b*c/d)
    return psi
# ---------------------------------------------------------------------


def psi_ecmwf(zol):
    """
    Calculates heat stability function for stable conditions
    for method ecmwf

    Parameters
    ----------
    zol : float
        height over MO length

    Returns
    -------
    psit : float
    """
    # eq (3.22) p. 37 IFS Documentation cy46r1
    a, b, c, d = 1, 2/3, 5, 0.35
    psit = -b*(zol-c/d)*np.exp(-d*zol)-np.power(1+(2/3)*a*zol, 1.5)-(b*c)/d+1
    return psit
# ---------------------------------------------------------------------


def psit_26(zol):
    """
    Computes temperature structure function as in C35

    Parameters
    ----------
    zol : float
        height over MO length

    Returns
    -------
    psi : float
    """
    b, d = 2/3, 0.35
    dzol = np.minimum(d*zol, 50)
    psi = -1*((1+b*zol)**1.5+b*(zol-14.28)*np.exp(-dzol)+8.525)
    k = np.where(zol < 0)
    x = np.sqrt(1-15*zol[k])
    psik = 2*np.log((1+x)/2)
    x = np.power(1-34.15*zol[k], 1/3)
    psic = (1.5*np.log((1+x+np.power(x, 2))/3)-np.sqrt(3) *
            np.arctan((1+2*x)/np.sqrt(3))+4*np.arctan(1)/np.sqrt(3))
    f = np.power(zol[k], 2)/(1+np.power(zol[k], 2))
    psi[k] = (1-f)*psik+f*psic
    return psi
# ---------------------------------------------------------------------


def psi_conv(zol, meth):
    """
    Calculates heat stability function for unstable conditions

    Parameters
    ----------
    zol : float
        height over MO length
    meth : str
        parameterisation method

    Returns
    -------
    psit : float
    """
    coeffs = get_stabco(meth)
    alpha, beta = coeffs[0], coeffs[1]
    xtmp = np.power(1-alpha*zol, beta)
    psit = 2*np.log((1+np.power(xtmp, 2))*0.5)
    return psit
# ---------------------------------------------------------------------


def psi_stab(zol, meth):
    """
    Calculates heat stability function for stable conditions

    Parameters
    ----------
    zol : float
        height over MO length
    meth : str
        parameterisation method

    Returns
    -------
    psit : float
    """
    coeffs = get_stabco(meth)
    gamma = coeffs[2]
    psit = -gamma*zol
    return psit
# ---------------------------------------------------------------------


def psim_ecmwf(zol):
    """
    Calculates momentum stability function for method ecmwf

    Parameters
    ----------
    zol : float
        height over MO length

    Returns
    -------
    psim : float
    """
    # eq (3.20, 3.22) p. 37 IFS Documentation cy46r1
    coeffs = get_stabco("ecmwf")
    alpha, beta = coeffs[0], coeffs[1]
    xtmp = np.power(1-alpha*zol, beta)
    a, b, c, d = 1, 2/3, 5, 0.35
    psim = np.where(zol < 0, np.pi/2-2*np.arctan(xtmp) +
                    np.log((np.power(1+xtmp, 2)*(1+np.power(xtmp, 2)))/8),
                    -b*(zol-c/d)*np.exp(-d*zol)-a*zol-(b*c)/d)
    return psim
# ---------------------------------------------------------------------


def psiu_26(zol, meth):
    """
    Computes velocity structure function C35

    Parameters
    ----------
    zol : float
        height over MO length

    Returns
    -------
    psi : float
    """
    if (meth == "C30"):
        dzol = np.minimum(0.35*zol, 50) # stable
        psi = -1*((1+zol)+0.6667*(zol-14.28)*np.exp(-dzol)+8.525)
        k = np.where(zol < 0) # unstable
        x = (1-15*zol[k])**0.25
        psik = (2*np.log((1+x)/2)+np.log((1+x*x)/2)-2*np.arctan(x) +
                2*np.arctan(1))
        x = (1-10.15*zol[k])**(1/3)
        psic = (1.5*np.log((1+x+x*x)/3) -
                np.sqrt(3)*np.arctan((1+2*x)/np.sqrt(3)) +
                4*np.arctan(1)/np.sqrt(3))
        f = zol[k]**2/(1+zol[k]**2)
        psi[k] = (1-f)*psik+f*psic
    elif (meth == "C35"): #  or meth == "C40"
        dzol = np.minimum(50, 0.35*zol)  # stable
        a, b, c, d = 0.7, 3/4, 5, 0.35
        psi = -1*(a*zol+b*(zol-c/d)*np.exp(-dzol)+b*c/d)
        k = np.where(zol < 0)  # unstable
        x = np.power(1-15*zol[k], 1/4)
        psik = 2*np.log((1+x)/2)+np.log((1+x*x)/2)-2*np.arctan(x)+2*np.arctan(1)
        x = np.power(1-10.15*zol[k], 1/3)
        psic = (1.5*np.log((1+x+np.power(x, 2))/3)-np.sqrt(3) *
                np.arctan((1+2*x)/np.sqrt(3))+4*np.arctan(1)/np.sqrt(3))
        f = np.power(zol[k], 2)/(1+np.power(zol[k], 2))
        psi[k] = (1-f)*psik+f*psic
    return psi
#------------------------------------------------------------------------------



def psim_conv(zol, meth):
    """
    Calculates momentum stability function for unstable conditions

    Parameters
    ----------
    zol : float
        height over MO length
    meth : str
        parameterisation method

    Returns
    -------
    psim : float
    """
    coeffs = get_stabco(meth)
    alpha, beta = coeffs[0], coeffs[1]
    xtmp = np.power(1-alpha*zol, beta)
    psim = (2*np.log((1+xtmp)*0.5)+np.log((1+np.power(xtmp, 2))*0.5) -
            2*np.arctan(xtmp)+np.pi/2)
    return psim
# ---------------------------------------------------------------------


def psim_stab(zol, meth):
    """
    Calculates momentum stability function for stable conditions

    Parameters
    ----------
    zol : float
        height over MO length
    meth : str
        parameterisation method

    Returns
    -------
    psim : float
    """
    coeffs = get_stabco(meth)
    gamma = coeffs[2]
    psim = -gamma*zol
    return psim
# ---------------------------------------------------------------------


def get_gust(beta, Ta, usr, tsrv, zi, grav):
    """
    Computes gustiness

    Parameters
    ----------
    beta : float
        constant
    Ta : float
        air temperature   [K]
    usr : float
        friction velocity [m/s]
    tsrv : float
        star virtual temperature of air [K]
    zi : int
        scale height of the boundary layer depth [m]
    grav : float
        gravity

    Returns
    -------
    ug : float        [m/s]
    """
    if (np.nanmax(Ta) < 200):  # convert to K if in Celsius
        Ta = Ta+273.16
    # minus sign to allow cube root
    Bf = (-grav/Ta)*usr*tsrv
    ug = np.ones(np.shape(Ta))*0.2
    ug = np.where(Bf > 0, beta*np.power(Bf*zi, 1/3), 0.2)
    return ug
# ---------------------------------------------------------------------

def get_strs(hin, monob, wind, zo, zot, zoq, dt, dq, cd, ct, cq, meth):
    """
    calculates star wind speed, temperature and specific humidity

    Parameters
    ----------
    hin : float
        sensor heights [m]
    monob : float
        M-O length     [m]
    wind : float
        wind speed     [m/s]
    zo : float
        momentum roughness length    [m]
    zot : float
        temperature roughness length [m]
    zoq : float
        moisture roughness length    [m]
    dt : float
        temperature difference       [K]
    dq : float
        specific humidity difference [g/kg]
    ct : float
        temperature exchange coefficient
    cq : float
        moisture exchange coefficient
    meth : str
        bulk parameterisation method option: "S80", "S88", "LP82", "YT96", "UA",
        "NCAR", "C30", "C35", "ecmwf", "Beljaars"

    Returns
    -------
    usr : float
        friction wind speed [m/s]
    tsr : float
        star temperature    [K]
    qsr : float
        star specific humidity [g/kg]

    """
    if (meth == "UA"):

        # Zeng et al. 1998
        # away from extremes UA follows e.g. S80
        usr = wind*np.power(cd, 1/2)
        tsr = ct*wind*dt/usr
        qsr = cq*wind*dq/usr

        # momentum
        hol0 = hin[0]/np.copy(monob)
        # very unstable (Zeng et al. 1998 eq 7)
        usr = np.where(hol0 <= -1.574, wind*kappa / (np.log(-1.574*monob/zo)-psim_calc(-1.574, meth) + psim_calc(zo/monob, meth) + 1.14*(np.power(-hin[0]/monob, 1/3) - np.power(1.574, 1/3))), usr)
        # very stable (Zeng et al. 1998 eq 10)
        usr = np.where(hol0 > 1, wind*kappa/(np.log(monob/zo)+5 - 5*zo/monob + 5*np.log(hin[0]/monob) + hin[0]/monob-1), usr)

        # temperature
        hol1 = hin[1]/np.copy(monob)
        # very unstable (Zeng et al. 1998 eq 11)
        tsr = np.where(hol1 < -0.465, kappa*dt / (np.log((-0.465*monob)/zot) - psit_calc(-0.465, meth)+0.8*(np.power(0.465, -1/3) - np.power(-hin[1]/monob, -1/3))), tsr)
        # very stable (Zeng et al. 1998 eq 14)
        tsr = np.where(hol1 > 1, kappa*(dt) / (np.log(monob/zot)+5 - 5*zot/monob+5*np.log(hin[1]/monob) + hin[1]/monob-1), tsr)

        # humidity
        hol2 = hin[2]/monob
        # very unstable (Zeng et al. 1998 eq 11)
        qsr = np.where(hol2 < -0.465, kappa*dq / (np.log((-0.465*monob)/zoq) - psit_calc(-0.465, meth)+psit_calc(zoq/monob, meth) + 0.8*(np.power(0.465, -1/3) - np.power(-hin[2]/monob, -1/3))), qsr)
        # very stable (Zeng et al. 1998 eq 14)
        qsr = np.where(hol2 > 1, kappa*dq / (np.log(monob/zoq)+5-5*zoq/monob + 5*np.log(hin[2]/monob) + hin[2]/monob-1), qsr)

    else:
        usr = wind*np.power(cd, 1/2)
        tsr = ct*wind*dt/usr
        qsr = cq*wind*dq/usr
    return usr, tsr, qsr
# ---------------------------------------------------------------------



# --------------------------------------------------------------------------------
# --------------------------------------------------------------------------------
# ---------------------------------------------------------------------
 
def get_tsrv(tsr, qsr, Ta, qair):
    """
    calculates virtual star temperature
 
    Parameters
    ----------
    tsr : float
        star temperature (K)
    qsr : float
        star specific humidity (g/kg)
    Ta : float
        air temperature (K)
    qair : float
        air specific humidity (g/kg)
 
    Returns
    -------
    tsrv : float
        virtual star temperature (K)
 
    """
    # as in aerobulk One_on_L in mod_phymbl.f90
    tsrv = tsr*(1+0.6077*qair)+0.6077*Ta*qsr
    return tsrv
 
# ---------------------------------------------------------------------
 
def get_Rb(grav, usr, hin_u, hin_t, tv, dtv, wind, monob, meth):
    """
    calculates bulk Richardson number
 
    Parameters
    ----------
    grav : float
        acceleration due to gravity (m/s2)
    usr : float
        friction wind speed (m/s)
    hin_u : float
        u sensor height (m)
    hin_t : float
        t sensor height (m)
    tv : float
        virtual temperature (K)
    dtv : float
        virtual temperature difference, air and sea (K)
    wind : float
        wind speed (m/s)
    monob : float
        Monin-Obukhov length from previous iteration step (m)
    psim : float
        momentum stability function
    meth : str
        bulk parameterisation method option: "S80", "S88", "LP82", "YT96",
        "UA", "NCAR", "C30", "C35", "ecmwf", "Beljaars"
 
    Returns
    -------
    Rb  : float
       Richardson number
 
    """
    # now input dtv
    #tvs = sst*(1+0.6077*qsea) # virtual SST
    #dtv = tv - tvs          # virtual air - sea temp. diff
    # adjust wind to t measurement height
    uz = (wind-usr/kappa*(np.log(hin_u/hin_t)-psim_calc(hin_u/monob, meth) +
                          psim_calc(hin_t/monob, meth)))
    Rb = grav*dtv*hin_t/(tv*uz*uz)

    return Rb
 
# ---------------------------------------------------------------------
 
def get_LRb(Rb, hin_t, monob, zo, zot, meth):
    """
    calculates Monin-Obukhov length following ecmwf (IFS Documentation cy46r1)
    default for methods ecmwf and Beljaars
 
    Parameters
    ----------
    Rb  : float
       Richardson number
    hin_t : float
        t sensor height (m)
    monob : float
        Monin-Obukhov length from previous iteration step (m)
    zo   : float
        surface roughness       (m)
    zot   : float
        temperature roughness length       (m)
    meth : str
        bulk parameterisation method option: "S80", "S88", "LP82", "YT96",
        "UA", "NCAR", "C30", "C35", "ecmwf", "Beljaars"
 
    Returns
    -------
    monob : float
        M-O length (m)
 
    """
    zol = (Rb*(np.power(np.log((hin_t+zo)/zo)-psim_calc((hin_t+zo) /
                monob, meth) + psim_calc(zo/monob, meth), 2) /
               (np.log((hin_t+zo)/zot) -
               psit_calc((hin_t+zo)/monob, meth) +
               psit_calc(zot/monob, meth))))
    monob = hin_t/zol
    return monob
 
# ---------------------------------------------------------------------
 
def get_Ltsrv(tsrv, grav, tv, usr):
    """
    calculates Monin-Obukhov length from tsrv
 
    Parameters
    ----------
    tsrv : float
        virtual star temperature (K)
    grav : float
        acceleration due to gravity (m/s2)
    tv : float
        virtual temperature (K)
    usr : float
        friction wind speed (m/s)
 
    Returns
    -------
    monob : float
        M-O length (m)
 
    """
    #tmp = (g*kappa*tsrv /
    #        np.maximum(np.power(usr, 2)*Ta*(1+0.6077*qair), 1e-9))
    #tmp = np.minimum(np.abs(tmp), 200)*np.sign(tmp)
    #monob = 1/np.copy(tmp)

    tsrv = np.maximum(np.abs(tsrv), 1e-9)*np.sign(tsrv)
    monob = (np.power(usr, 2)*tv)/(grav*kappa*tsrv)
        
    return monob
 
# ---------------------------------------------------------------------