AirSeaFluxCode.py

import numpy as np
import sys
import logging
#import metpy.constants as mpcon
from flux_subs import (kappa, CtoK, get_heights, cdn_calc, cd_calc, get_skin,
                       psim_calc, psit_calc, ctcq_calc, ctcqn_calc, get_gust,
                       gc, q_calc, qsea_calc, qsat26sea, qsat26air,
                       visc_air, psit_26, psiu_26)


def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi=600,
                   Rl=None, Rs=None, jcool=1, meth="S88", n=10):
    """ Calculates momentum and heat fluxes using different parameterizations

    Parameters
    ----------
        meth : str
            "S80","S88","LP82","YT96","UA","LY04","C30","C35","C40","ERA5"
        spd : float
            relative wind speed in m/s (is assumed as magnitude difference
            between wind and surface current vectors)
        T : float
            air temperature in K (will convert if < 200)
        SST : float
            sea surface temperature in K (will convert if < 200)
        lat : float
            latitude (deg)
        RH : float
            relative humidity in %
        P : float
            air pressure
        hin : float
            sensor heights in m (array of 1->3 values: 1 -> u=t=q; /
            2 -> u,t=q; 3 -> u,t,q ) default 10m
        hout : float
            output height, default is 10m
        zi : int
            PBL height (m) called in C35
        Rl : float
            downward longwave radiation (W/m^2)
        Rs : float
            downward shortwave radiation (W/m^2)
        jcool : bool
            0 if sst is true ocean skin temperature called in COARE
        n : int
            number of iterations

    Returns
    -------
        res : array that contains
                       1. momentum flux (W/m^2)
                       2. sensible heat (W/m^2)
                       3. latent heat (W/m^2)
                       4. Monin-Obhukov length (mb)
                       5. drag coefficient (cd)
                       6. neutral drag coefficient (cdn)
                       7. heat exhange coefficient (ct)
                       8. neutral heat exhange coefficient (ctn)
                       9. moisture exhange coefficient (cq)
                       10. neutral moisture exhange coefficient (cqn)
                       11. star virtual temperature (tsrv)
                       12. star temperature (tsr)
                       13. star humidity (qsr)
                       14. star velocity (usr)
                       15. momentum stability function (psim)
                       16. heat stability funciton (psit)
                       17. 10m neutral velocity (u10n)
                       18. 10m neutral temperature (t10n)
                       19. 10m neutral virtual temperature (tv10n)
                       20. 10m neutral specific humidity (q10n)
                       21. surface roughness length (zo)
                       22. heat roughness length (zot)
                       23. moisture roughness length (zoq)
                       24. velocity at reference height (urefs)
                       25. temperature at reference height (trefs)
                       26. specific humidity at reference height (qrefs)
                       27. number of iterations until convergence
        ind : int
            the indices in the matrix for the points that did not converge
            after the maximum number of iterations
    The code is based on bform.f and flux_calc.R modified by S. Biri
    """
    logging.basicConfig(filename='flux_calc.log',
                        format='%(asctime)s %(message)s',level=logging.INFO)
    ref_ht, tlapse = 10, 0.0098   # reference height, lapse rate
    hh_in = get_heights(hin)      # heights of input measurements/fields
    hh_out = get_heights(hout)    # desired height of output variables
    if np.all(np.isnan(lat)):     # set latitude to 45deg if empty
        lat=45*np.ones(spd.shape)
    g = gc(lat, None)             # acceleration due to gravity
    ctn, ct, cqn, cq = (np.zeros(spd.shape)*np.nan, np.zeros(spd.shape)*np.nan,
                        np.zeros(spd.shape)*np.nan, np.zeros(spd.shape)*np.nan)
    # if input values are nan break
    if (np.all(np.isnan(spd)) or np.all(np.isnan(T)) or np.all(np.isnan(SST))):
        sys.exit("input wind, T or SST is empty")
        logging.debug('all input is nan')
    if (np.all(np.isnan(RH)) and meth == "C35"):
        RH = np.ones(spd.shape)*80  # if empty set to default for COARE3.5
    elif (np.all(np.isnan(RH))):
        sys.exit("input RH is empty")
        logging.debug('input RH is empty')
    if (np.all(np.isnan(P)) and (meth == "C30" or meth == "C40")):
        P = np.ones(spd.shape)*1015  # if empty set to default for COARE3.0
    elif ((np.all(np.isnan(P))) and meth == "C35"):
        P = np.ones(spd.shape)*1013  # if empty set to default for COARE3.5
    elif (np.all(np.isnan(P))):
        sys.exit("input P is empty")
        logging.debug('input P is empty')
    if (np.all(np.isnan(Rl)) and meth == "C30"):
        Rl = np.ones(spd.shape)*150    # set to default for COARE3.0
    elif ((np.all(np.isnan(Rl)) and meth == "C35") or
          (np.all(np.isnan(Rl)) and meth == "C40")):
        Rl = np.ones(spd.shape)*370    # set to default for COARE3.5
    if (np.all(np.isnan(Rs)) and meth == "C30"):
        Rs = np.ones(spd.shape)*370  # set to default for COARE3.0
    elif ((np.all(np.isnan(Rs))) and (meth == "C35" or meth == "C40")):
        Rs = np.ones(spd.shape)*150  # set to default for COARE3.5
    if ((np.all(np.isnan(zi))) and (meth == "C30" or meth == "C35" or
        meth == "C40")):
        zi = 600  # set to default for COARE3.5
    elif ((np.all(np.isnan(zi))) and (meth == "ERA5" or meth == "UA")):
        zi = 1000
    if (np.all(np.isnan(lat)) and (meth == "C30" or meth == "C35" or
        meth == "C40")):
        lat=45*np.ones(np.shape(spd))
    ####
    th = np.where(np.nanmax(T) < 200, (np.copy(T)+CtoK) *
                  np.power(1000/P,287.1/1004.67),
                  np.copy(T)*np.power(1000/P,287.1/1004.67))  # potential T
    Ta = np.where(np.nanmax(T) < 200, np.copy(T)+CtoK+tlapse*hh_in[1],
                  np.copy(T)+tlapse*hh_in[1])  # convert to Kelvin if needed
    sst = np.where(np.nanmax(SST) < 200, np.copy(SST)+CtoK, np.copy(SST))
    if (meth == "C30" or meth == "C35" or meth == "C40"  or meth == "UA" or
        meth == "ERA5"):
        qsea = qsat26sea(sst, P)/1000  # surface water specific humidity (g/kg)
        Q, _ = qsat26air(T, P, RH)     # specific humidity of air (g/kg)
        qair = Q/1000
        del Q
        logging.info('method %s | qsea:%s, qair:%s', meth,
                     np.ma.median(qsea[~np.isnan(qsea)]),
                     np.ma.median(qair[~np.isnan(qair)]))
    else:
        qsea = qsea_calc(sst, P)
        qair = q_calc(Ta, RH, P)
        logging.info('method %s | qsea:%s, qair:%s', meth,
                     np.ma.median(qsea[~np.isnan(qsea)]),
                     np.ma.median(qair[~np.isnan(qair)]))
    if (np.all(np.isnan(qsea)) or np.all(np.isnan(qair))):
        print("qsea and qair cannot be nan")
        logging.info('method %s qsea and qair cannot be nan | sst:%s, Ta:%s,'
                      'P:%s, RH:%s', meth, np.ma.median(sst[~np.isnan(sst)]),
                      np.ma.median(Ta[~np.isnan(Ta)]),
                      np.ma.median(P[~np.isnan(P)]),
                      np.ma.median(RH[~np.isnan(RH)]))
    # first guesses
    dt = Ta - sst
    dq = qair - qsea
    t10n, q10n = np.copy(Ta), np.copy(qair)
    tv10n = t10n*(1 + 0.61*q10n)
    #  Zeng et al. 1998
    tv=th*(1.+0.61*qair)   # virtual potential T
    dtv=dt*(1.+0.61*qair)+0.61*th*dq
    # ------------
#    rho = P*100/(287.1*tv10n)
    rho = P*100/(287.1*(T+CtoK)*(1+0.61*qair))
    # rho=P*100/(287.1*sst*(1+0.61*qsea))  # Zeng et al. 1998
    lv = (2.501-0.00237*SST)*1e6
    cp = 1004.67*(1 + 0.00084*qsea)
#    cp = 1004.67  # Zeng et al. 1998, C3.0, C3.5
    u10n = np.copy(spd)
    monob = -100*np.ones(u10n.shape)  # Monin-Obukhov length
    cdn = cdn_calc(u10n, Ta, None, lat, meth)
    psim, psit, psiq = (np.zeros(u10n.shape), np.zeros(u10n.shape),
                        np.zeros(u10n.shape))
    cd = cd_calc(cdn, hh_in[0], ref_ht, psim)
    tsr, tsrv = np.zeros(u10n.shape), np.zeros(u10n.shape)
    qsr = np.zeros(u10n.shape)
    if (meth == "UA"):
        wind = np.where(dtv >= 0,np.where(spd > 0.1, spd, 0.1),
                        np.sqrt(np.power(np.copy(spd),2)+np.power(0.5,2)))
        usr = 0.06
        for i in range(5):
            zo = 0.013*np.power(usr,2)/g+0.11*visc_air(Ta)/usr
            usr=kappa*wind/np.log(hh_in[0]/zo)
        Rb = g*hh_in[0]*dtv/(tv*wind**2)
        zol = np.where(Rb >= 0, Rb*np.log(hh_in[0]/zo) /
                       (1-5*np.where(Rb < 0.19, Rb, 0.19)),
                       Rb*np.log(hh_in[0]/zo))
        monob = hh_in[0]/zol
        zo = 0.013*np.power(usr, 2)/g + 0.11*visc_air(Ta)/usr
        zot = zo/np.exp(2.67*np.power(usr*zo/visc_air(Ta), 0.25)-2.57)
        zoq = zot
        logging.info('method %s | wind:%s, usr:%s,'
                     'zo:%s, zot:%s, zoq:%s, Rb:%s, monob:%s', meth,
                     np.ma.median(wind[~np.isnan(wind)]),
                     np.ma.median(usr[~np.isnan(usr)]),
                     np.ma.median(zo[~np.isnan(zo)]),
                     np.ma.median(zot[~np.isnan(zot)]),
                     np.ma.median(zoq[~np.isnan(zoq)]),
                     np.ma.median(Rb[~np.isnan(Rb)]),
                     np.ma.median(monob[~np.isnan(monob)]))
    elif (meth == "ERA5"):
        wind = np.sqrt(np.power(np.copy(spd), 2)+np.power(0.5, 2))
        usr = np.sqrt(cd*wind**2)
        Rb = ((g*hh_in[0]*((2*dt)/(Ta+sst-g*hh_in[0]) +
                0.61*dq))/np.power(wind, 2))
        zo = 0.11*visc_air(Ta)/usr+0.018*np.power(usr, 2)/g
        zot = 0.40*visc_air(Ta)/usr
        zol = (Rb*((np.log((hh_in[0]+zo)/zo)-psim_calc((hh_in[0]+zo) /
               monob, meth)+psim_calc(zo/monob, meth)) /
               (np.log((hh_in[0]+zo)/zot) -
               psit_calc((hh_in[0]+zo)/monob, meth) +
               psit_calc(zot/monob, meth))))
        monob = hh_in[0]/zol
        logging.info('method %s | wind:%s, usr:%s,'
                     'zo:%s, zot:%s, zoq:%s, Rb:%s, monob:%s', meth,
                     np.ma.median(wind[~np.isnan(wind)]),
                     np.ma.median(usr[~np.isnan(usr)]),
                     np.ma.median(zo[~np.isnan(zo)]),
                     np.ma.median(zot[~np.isnan(zot)]),
                     np.ma.median(Rb[~np.isnan(Rb)]),
                     np.ma.median(monob[~np.isnan(monob)]))
    elif (meth == "C30" or meth == "C35" or meth == "C40"):
        wind = np.sqrt(np.power(np.copy(spd), 2)+np.power(0.5, 2))
        usr = 0.035*wind*np.log(10/1e-4)/np.log(hh_in[0]/1e-4)
        a = 0.011*np.ones(T.shape)
        a = np.where(wind > 10, 0.011+(wind-10)/(18-10)*(0.018-0.011),
                         np.where(wind > 18, 0.018, a))
        zo = a*np.power(usr, 2)/g+0.11*visc_air(T)/usr
        rr = zo*usr/visc_air(T)
        zoq = np.minimum(5e-5/np.power(rr, 0.6), 1.15e-4)
        zot=zoq
        Rb = g*hh_in[0]*dtv/((T+CtoK)*np.power(wind, 2))
        zol =  (Rb*((np.log((hh_in[0]+zo)/zo)-psim_calc((hh_in[0]+zo) /
               monob, meth)+psim_calc(zo/monob, meth)) /
               (np.log((hh_in[0]+zo)/zot) -
               psit_calc((hh_in[0]+zo)/monob, meth) +
               psit_calc(zot/monob, meth))))
        monob = hh_in[0]/zol
#        wetc = 0.622*lv*qsea/(287.1*np.power(sst, 2))
        tkt = 0.001*np.ones(T.shape)
        Rnl = 0.97*(5.67e-8*np.power(sst-0.3*jcool+CtoK, 4)-Rl)
        dter, dqer = np.ones(T.shape)*0.3, np.zeros(T.shape)*np.nan
        logging.info('method %s | wind:%s, usr:%s,'
                     'zo:%s, zot:%s, zoq:%s, Rb:%s, monob:%s', meth,
                     np.ma.median(wind[~np.isnan(wind)]),
                     np.ma.median(usr[~np.isnan(usr)]),
                     np.ma.median(zo[~np.isnan(zo)]),
                     np.ma.median(zot[~np.isnan(zot)]),
                     np.ma.median(Rb[~np.isnan(Rb)]),
                     np.ma.median(monob[~np.isnan(monob)]))
    else:
        wind = np.copy(spd)
        zo, zot = 0.0001*np.ones(u10n.shape), 0.0001*np.ones(u10n.shape)
        usr = np.sqrt(cd * wind**2)
        logging.info('method %s | wind:%s, usr:%s,'
                     'zo:%s, zot:%s, zoq:%s, Rb:%s, monob:%s', meth,
                     np.ma.median(wind[~np.isnan(wind)]),
                     np.ma.median(usr[~np.isnan(usr)]),
                     np.ma.median(zo[~np.isnan(zo)]),
                     np.ma.median(zot[~np.isnan(zot)]),
                     np.ma.median(monob[~np.isnan(monob)]))
    # tolerance for u,t,q,usr,tsr,qsr
    tol = np.array([1e-06, 0.01, 5e-05, 1e-06, 0.001, 5e-07])
    it, ind, ii, itera = 0, np.where(spd > 0), True, np.zeros(spd.shape)*np.nan
    while np.any(ii):
        it += 1
        if it > n:
            break
        old = np.array([np.copy(u10n[ind]), np.copy(t10n[ind]),
                       np.copy(q10n[ind]), np.copy(usr[ind]),
                       np.copy(tsr[ind]), np.copy(qsr[ind])])
        cdn[ind] = cdn_calc(u10n[ind], Ta[ind], None, lat[ind], meth)
        if (np.all(np.isnan(cdn))):
            break  # sys.exit("cdn cannot be nan")
            logging.info('break %s at iteration %s cdn<0', meth, it)
        zo[ind] = ref_ht/np.exp(kappa/np.sqrt(cdn[ind]))
        psim[ind] = psim_calc(hh_in[0]/monob[ind], meth)
        cd[ind] = cd_calc(cdn[ind], hh_in[0], ref_ht, psim[ind])
        ctn[ind], cqn[ind] = ctcqn_calc(hh_in[1]/monob[ind], cdn[ind],
                                        u10n[ind], zo[ind], Ta[ind], meth)
        psit[ind] = psit_calc(hh_in[1]/monob[ind], meth)
        psiq[ind] = psit_calc(hh_in[2]/monob[ind], meth)
        ct[ind], cq[ind] = ctcq_calc(cdn[ind], cd[ind], ctn[ind], cqn[ind],
                                     hh_in[1], hh_in[2], ref_ht,
                                     psit[ind], psiq[ind])
        if (meth == "UA"):
            usr[ind] = np.where(hh_in[0]/monob[ind] < -1.574, kappa*wind[ind] /
                                (np.log(-1.574*monob[ind]/zo[ind]) -
                                psim_calc(-1.574, meth) +
                                psim_calc(zo[ind]/monob[ind], meth) +
                                1.14*(np.power(-hh_in[0]/monob[ind],1/3) -
                                np.power(1.574,1/3))),
                                np.where((hh_in[0]/monob[ind] > -1.574) &
                                (hh_in[0]/monob[ind] < 0),
                                kappa*wind[ind]/(np.log(hh_in[0]/zo[ind]) -
                                psim_calc(hh_in[0]/monob[ind], meth) +
                                psim_calc(zo[ind]/monob[ind], meth)),
                                np.where((hh_in[0]/monob[ind] > 0) &
                                (hh_in[0]/monob[ind] < 1),
                                kappa*wind[ind]/(np.log(hh_in[0]/zo[ind]) +
                                5*hh_in[0]/monob[ind]-5*zo[ind]/monob[ind]),
                                kappa*wind[ind]/(np.log(monob[ind]/zo[ind]) +
                                5-5*zo[ind]/monob[ind] +
                                5*np.log(hh_in[0]/monob[ind]) +
                                hh_in[0]/monob[ind]-1))))
                                # Zeng et al. 1998 (7-10)
            tsr[ind] = np.where(hh_in[1]/monob[ind] < -0.465, kappa*dt[ind] /
                                (np.log((-0.465*monob[ind])/zot[ind]) -
                                psit_calc(-0.465, meth)+0.8 *
                                (np.power(0.465,-1/3) -
                                np.power(-hh_in[1]/monob[ind],-1/3))),
                                np.where((hh_in[1]/monob[ind]>-0.465) &
                                (hh_in[1]/monob[ind]<0),
                                kappa*dt[ind]/(np.log(hh_in[1]/zot[ind]) -
                                psit_calc(hh_in[1]/monob[ind], meth) +
                                psit_calc(zot[ind]/monob[ind], meth)),
                                np.where((hh_in[1]/monob[ind]>0) &
                                (hh_in[1]/monob[ind]<1),
                                kappa*dt[ind]/(np.log(hh_in[1]/zot[ind]) +
                                5*hh_in[1]/monob[ind]-5*zot[ind]/monob[ind]),
                                kappa*dt[ind]/(np.log(monob[ind]/zot[ind])+5 -
                                5**zot[ind]/monob[ind] +
                                5*np.log(hh_in[1]/monob[ind]) +
                                hh_in[1]/monob[ind]-1))))
                                # Zeng et al. 1998 (11-14)
            qsr[ind] = np.where(hh_in[2]/monob[ind] < -0.465, kappa*dq[ind] /
                                (np.log((-0.465*monob[ind])/zoq[ind]) -
                                psit_calc(-0.465, meth) +
                                psit_calc(zoq[ind]/monob[ind], meth) +
                                0.8*(np.power(0.465,-1/3) -
                                np.power(-hh_in[2]/monob[ind],-1/3))),
                                np.where((hh_in[2]/monob[ind]>-0.465) &
                                (hh_in[2]/monob[ind]<0),
                                kappa*dq[ind]/(np.log(hh_in[1]/zot[ind]) -
                                psit_calc(hh_in[2]/monob[ind], meth) +
                                psit_calc(zoq[ind]/monob[ind], meth)),
                                np.where((hh_in[2]/monob[ind]>0) &
                                (hh_in[2]/monob[ind]<1), kappa*dq[ind] /
                                (np.log(hh_in[1]/zoq[ind]) +
                                5*hh_in[2]/monob[ind]-5*zoq[ind]/monob[ind]),
                                kappa*dq[ind]/(np.log(monob[ind]/zoq[ind])+5 -
                                5*zoq[ind]/monob[ind] +
                                5*np.log(hh_in[2]/monob[ind]) +
                                hh_in[2]/monob[ind]-1))))
        elif (meth == "C30" or meth == "C35" or meth == "C40"):
            usr[ind] = (wind[ind]*kappa/(np.log(hh_in[0]/zo[ind]) -
                        psiu_26(hh_in[0]/monob[ind], meth)))
            dter[ind], dqer[ind], tkt[ind] = get_skin(sst[ind], qsea[ind],
                                                      rho[ind], Rl[ind],
                                                      Rs[ind], Rnl[ind],
                                                      cp[ind], lv[ind],
                                                      np.copy(tkt[ind]),
                                                      usr[ind], tsr[ind],
                                                      qsr[ind], lat[ind])
            logging.info('method %s | dter = %s | Rnl = %s '
                         '| usr = %s | tsr = %s | qsr = %s', meth,
                         np.ma.median(dter[~np.isnan(dter)]),
                         np.ma.median(Rnl[~np.isnan(Rnl)]),
                         np.ma.median(usr[~np.isnan(usr)]),
                         np.ma.median(tsr[~np.isnan(tsr)]),
                         np.ma.median(qsr[~np.isnan(qsr)]))
            qsr[ind] = ((dq[ind]+dqer[ind]*jcool)*(kappa /
                        (np.log(hin[2]/zoq[ind])-psit_26(hin[2]/monob[ind]))))
            tsr[ind] = ((dt[ind]+dter[ind]*jcool)*(kappa /
                        (np.log(hin[1]/zot[ind])-psit_26(hin[1]/monob[ind]))))
            Rnl[ind] = 0.97*(5.67e-8*np.power(SST[ind] -
                             dter[ind]*jcool+CtoK, 4)-Rl[ind])
        else:
            usr[ind] = np.sqrt(cd[ind]*wind[ind]**2)
            tsr[ind] = ct[ind]*wind[ind]*dt[ind]/usr[ind]
            qsr[ind] = cq[ind]*wind[ind]*dq[ind]/usr[ind]
        fact = (np.log(hh_in[1]/ref_ht)-psit[ind])/kappa
        t10n[ind] = Ta[ind] - (tsr[ind]*fact)
        fact = (np.log(hh_in[2]/ref_ht)-psiq[ind])/kappa
        q10n[ind] = qair[ind] - (qsr[ind]*fact)
        tv10n[ind] = t10n[ind]*(1+0.61*q10n[ind])
        if (meth == "UA"):
            tsrv[ind] = tsr[ind]*(1.+0.61*qair[ind])+0.61*th[ind]*qsr[ind]
            monob[ind] = (tv[ind]*np.power(usr[ind], 2))/(kappa*g[ind]*tsrv[ind])
        elif (meth == "C30" or meth == "C35" or meth == "C40"):
            tsrv[ind] = tsr[ind]+0.61*(T[ind]+CtoK)*qsr[ind]
            zol[ind] = (kappa*g[ind]*hh_in[0]/(T[ind]+CtoK)*(tsr[ind] +
                        0.61*(T[ind]+CtoK)*qsr[ind])/np.power(usr[ind], 2))
            monob[ind] = hh_in[0]/zol[ind]
        elif (meth == "ERA5"):
            tsrv[ind] = tsr[ind]+0.61*t10n[ind]*qsr[ind]
            Rb[ind] = ((g[ind]*hh_in[0]*((2*dt[ind])/(Ta[ind]+sst[ind]-g[ind] *
                       hh_in[0])+0.61*dq[ind]))/np.power(wind[ind], 2))
            zo[ind] = (0.11*visc_air(Ta[ind])/usr[ind]+0.018 *
                       np.power(usr[ind], 2)/g[ind])
            zot[ind] = 0.40*visc_air(Ta[ind])/usr[ind]
            zol[ind] = (Rb[ind]*((np.log((hh_in[0]+zo[ind])/zo[ind]) -
                        psim_calc((hh_in[0]+zo[ind])/monob[ind], meth) +
                        psim_calc(zo[ind]/monob[ind], meth)) /
                        (np.log((hh_in[0]+zo[ind])/zot[ind]) -
                        psit_calc((hh_in[0]+zo[ind])/monob[ind], meth) +
                        psit_calc(zot[ind]/monob[ind], meth))))
            monob[ind] = hh_in[0]/zol[ind]
        else:
            tsrv[ind] = tsr[ind]+0.61*t10n[ind]*qsr[ind]
            monob[ind] = (tv10n[ind]*usr[ind]**2)/(g[ind]*kappa*tsrv[ind])
        psim[ind] = psim_calc(hh_in[0]/monob[ind], meth)
        psit[ind] = psit_calc(hh_in[1]/monob[ind], meth)
        psiq[ind] = psit_calc(hh_in[2]/monob[ind], meth)
        if (meth == "UA"):
            wind[ind] = np.where(dtv[ind] >= 0, np.where(spd[ind] > 0.1,
                                 spd[ind], 0.1),
                                 np.sqrt(np.power(np.copy(spd[ind]), 2) +
                                 np.power(get_gust(1,tv[ind], usr[ind],
                                 tsrv[ind], zi, lat[ind]), 2)))
                                 # Zeng et al. 1998 (20)
            u10n[ind] = (wind[ind]-(usr[ind]/kappa)*(np.log(hh_in[0]/10) -
                         psim[ind]))
            u10n[u10n < 0] = np.nan
        elif (meth == "C30" or meth == "C35" or meth == "C40"):
            wind[ind] = np.sqrt(np.power(np.copy(spd[ind]), 2) +
                                np.power(get_gust(1.2, Ta[ind], usr[ind],
                                tsrv[ind], zi, lat[ind]), 2))
            u10n[ind] = ((wind[ind] + usr[ind]/kappa*(np.log(10/hh_in[0])-
                         psiu_26(10/monob[ind], meth) +
                         psiu_26(hh_in[0]/monob[ind], meth)) +
                         psiu_26(10/monob[ind], meth)*usr[ind]/kappa /
                         (wind[ind]/spd[ind])))
            u10n[u10n < 0] = np.nan
        elif (meth == "ERA5"):
            wind[ind] = np.sqrt(np.power(np.copy(spd[ind]), 2) +
                                np.power(get_gust(1, Ta[ind], usr[ind],
                                tsrv[ind], zi, lat[ind]), 2))
            u10n[ind] = (spd[ind]+(usr[ind]/kappa)*(np.log(hh_in[0] /
                         ref_ht)-psim[ind]))
            u10n[u10n < 0] = np.nan
        else:
            u10n[ind] = (wind[ind]-(usr[ind]/kappa)*(np.log(hh_in[0]/10) -
                         psim[ind]))
            u10n[u10n < 0] = np.nan
        new = np.array([np.copy(u10n[ind]), np.copy(t10n[ind]),
                       np.copy(q10n[ind]), np.copy(usr[ind]),
                       np.copy(tsr[ind]), np.copy(qsr[ind])])
        d = np.abs(new-old)
        ind = np.where((d[0, :] > tol[0])+(d[1, :] > tol[1]) +
                       (d[2, :] > tol[2])+(d[3, :] > tol[3]) +
                       (d[4, :] > tol[4])+(d[5, :] > tol[5]))
        itera[ind] = np.ones(1)*it
        if np.shape(ind)[0] == 0:
            break
        else:
            ii = ((d[0, ind] > tol[0])+(d[1, ind] > tol[1]) +
                  (d[2, ind] > tol[2])+(d[3, ind] > tol[3]) +
                  (d[4, ind] > tol[4])+(d[5, ind] > tol[5]))
    logging.info('method %s | # of iterations:%s', meth, np.ma.median(it))
    logging.info('method %s | # of points that did not converge :%s', meth,
                  np.shape(ind))
    # calculate output parameters
#    rho = (0.34838*P)/(tv10n)
    rho = P*100./(287.1*(T+CtoK)*(1+0.61*qair))  # C35
    t10n = t10n-(273.16+tlapse*ref_ht)
    sensible = -1*tsr*usr*cp*rho
    latent = -1*qsr*usr*lv*rho
    if (meth == "C30" or meth == "C35" or meth == "C40" or meth == "UA" or
        meth == "ERA5"):
        tau = rho*np.power(usr, 2)*(spd/wind)
    else:
        tau = rho*np.power(usr, 2)
    zo = ref_ht/np.exp(kappa/cdn**0.5)
    zot = ref_ht/(np.exp(kappa**2/(ctn*np.log(ref_ht/zo))))
    zoq = ref_ht/(np.exp(kappa**2/(cqn*np.log(ref_ht/zo))))
    urefs = (spd-(usr/kappa)*(np.log(hh_in[0]/hh_out[0])-psim +
             psim_calc(hh_out[0]/monob, meth)))
    trefs = (Ta-(tsr/kappa)*(np.log(hh_in[1]/hh_out[1])-psit +
             psit_calc(hh_out[0]/monob, meth)))
    trefs = trefs-(273.16+tlapse*hh_out[1])
    qrefs = (qair-(qsr/kappa)*(np.log(hh_in[2]/hh_out[2]) -
             psit+psit_calc(hh_out[2]/monob, meth)))
    res = np.zeros((27, len(spd)))
    res[0][:] = tau
    res[1][:] = sensible
    res[2][:] = latent
    res[3][:] = monob
    res[4][:] = cd
    res[5][:] = cdn
    res[6][:] = ct
    res[7][:] = ctn
    res[8][:] = cq
    res[9][:] = cqn
    res[10][:] = tsrv
    res[11][:] = tsr
    res[12][:] = qsr
    res[13][:] = usr
    res[14][:] = psim
    res[15][:] = psit
    res[16][:] = u10n
    res[17][:] = t10n
    res[18][:] = tv10n
    res[19][:] = q10n
    res[20][:] = zo
    res[21][:] = zot
    res[22][:] = zoq
    res[23][:] = urefs
    res[24][:] = trefs
    res[25][:] = qrefs
    res[26][:] = itera
    return res, ind