import numpy as np
import sys
import logging
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, qsat_sea, qsat_air, visc_air, psit_26, psiu_26)
def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi, Rl, Rs, jcool,
gust, meth, qmeth, n):
""" Calculates momentum and heat fluxes using different parameterizations
Parameters
----------
meth : str
"S80","S88","LP82","YT96","UA","LY04","C30","C35","C40","ERA5"
qmeth : str
is the saturation evaporation method to use amongst
"HylandWexler","Hardy","Preining","Wexler","GoffGratch","CIMO",
"MagnusTetens","Buck","Buck2","WMO","WMO2000","Sonntag","Bolton",
"IAPWS","MurphyKoop"]
default is Buck2
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 3x1 or 3xn) default 10m ([10, 10, 10])
hout : float
output height, default is 10m
zi : int
PBL height (m)
Rl : float
downward longwave radiation (W/m^2)
Rs : float
downward shortwave radiation (W/m^2)
jcool : int
0 if sst is true ocean skin temperature, else 1
gust : int
1 to include the effect of gustiness, else 0
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 function (psit)
17. moisture stability function (psiq)
18. 10m neutral velocity (u10n)
19. 10m neutral temperature (t10n)
20. 10m neutral virtual temperature (tv10n)
21. 10m neutral specific humidity (q10n)
22. surface roughness length (zo)
23. heat roughness length (zot)
24. moisture roughness length (zoq)
25. velocity at reference height (urefs)
26. temperature at reference height (trefs)
27. specific humidity at reference height (qrefs)
28. 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
h_in = get_heights(hin, len(spd)) # heights of input measurements/fields
h_out = get_heights(hout, 1) # 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
# 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
elif (np.all(np.isnan(Rl))):
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
elif (np.all(np.isnan(zi))):
zi = 800
####
th = np.where(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(T < 200, np.copy(T)+CtoK+tlapse*h_in[1],
np.copy(T)+tlapse*h_in[1]) # convert to Kelvin if needed
sst = np.where(SST < 200, np.copy(SST)+CtoK, np.copy(SST))
if qmeth:
qsea = qsat_sea(sst, P, qmeth)/1000 # surface water q (g/kg)
qair = qsat_air(T, P, RH, qmeth)/1000 # q of air (g/kg)
logging.info('method %s and q method %s | qsea:%s, qair:%s', meth,
qmeth, np.nanmedian(qsea), np.nanmedian(qair))
else:
qsea = qsat_sea(sst, P, "Buck2")/1000 # surface water q (g/kg)
qair = qsat_air(T, P, RH, "Buck2")/1000 # q of air (g/kg)
logging.info('method %s | qsea:%s, qair:%s', meth,
np.nanmedian(qsea[~np.isnan(qsea)]),
np.nanmedian(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.nanmedian(sst), np.nanmedian(Ta),
np.nanmedian(P), np.nanmedian(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)
lv = (2.501-0.00237*SST)*1e6
cp = 1004.67*(1 + 0.00084*qsea)
u10n = np.copy(spd)
monob = -100*np.ones(spd.shape) # Monin-Obukhov length
cdn = cdn_calc(u10n, Ta, None, lat, meth)
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)
psim, psit, psiq = (np.zeros(spd.shape), np.zeros(spd.shape),
np.zeros(spd.shape))
cd = cd_calc(cdn, h_in[0], ref_ht, psim)
tsr, tsrv = np.zeros(spd.shape), np.zeros(spd.shape)
qsr = np.zeros(spd.shape)
# jcool parameters
tkt = 0.001*np.ones(T.shape)
Rnl = 0.97*(5.67e-8*np.power(sst-0.3*jcool+CtoK, 4)-Rl)
dter = np.ones(T.shape)*0.3
dqer = dter*0.622*lv*qsea/(287.1*np.power(sst, 2))
if (gust == 1 and 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)))
elif (gust == 1):
wind = np.sqrt(np.power(np.copy(spd), 2)+np.power(0.5, 2))
elif (gust == 0):
wind = np.copy(spd)
if (meth == "UA"):
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(h_in[0]/zo)
Rb = g*h_in[0]*dtv/(tv*wind**2)
zol = np.where(Rb >= 0, Rb*np.log(h_in[0]/zo) /
(1-5*np.where(Rb < 0.19, Rb, 0.19)),
Rb*np.log(h_in[0]/zo))
monob = h_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.nanmedian(wind), np.nanmedian(usr), np.nanmedian(zo),
np.nanmedian(zot), np.nanmedian(zoq), np.nanmedian(Rb),
np.nanmedian(monob))
elif (meth == "ERA5"):
usr = np.sqrt(cd*np.power(wind, 2))
Rb = ((g*h_in[0]*((2*dt)/(Ta+sst-g*h_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((h_in[0]+zo)/zo)-psim_calc((h_in[0]+zo) /
monob, meth)+psim_calc(zo/monob, meth)) /
(np.log((h_in[0]+zo)/zot) -
psit_calc((h_in[0]+zo)/monob, meth) +
psit_calc(zot/monob, meth))))
monob = h_in[0]/zol
logging.info('method %s | wind:%s, usr:%s, '
'zo:%s, zot:%s, Rb:%s, monob:%s', meth,
np.nanmedian(wind), np.nanmedian(usr), np.nanmedian(zo),
np.nanmedian(zot), np.nanmedian(Rb), np.nanmedian(monob))
elif (meth == "C30" or meth == "C35" or meth == "C40"):
usr = np.sqrt(cd*np.power(wind, 2))
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*h_in[0]*dtv/((T+CtoK)*np.power(wind, 2))
zol = (Rb*((np.log((h_in[0]+zo)/zo)-psim_calc((h_in[0]+zo) /
monob, meth)+psim_calc(zo/monob, meth)) /
(np.log((h_in[0]+zo)/zot) -
psit_calc((h_in[0]+zo)/monob, meth) +
psit_calc(zot/monob, meth))))
monob = h_in[0]/zol
logging.info('method %s | wind:%s, usr:%s, '
'zo:%s, zot:%s, Rb:%s, monob:%s', meth,
np.nanmedian(wind), np.nanmedian(usr), np.nanmedian(zo),
np.nanmedian(zot), np.nanmedian(Rb), np.nanmedian(monob))
else:
zo, zot = 0.0001*np.ones(spd.shape), 0.0001*np.ones(spd.shape)
usr = np.sqrt(cd*np.power(wind, 2))
logging.info('method %s | wind:%s, usr:%s, '
'zo:%s, zot:%s, monob:%s', meth,
np.nanmedian(wind), np.nanmedian(usr), np.nanmedian(zo),
np.nanmedian(zot), np.nanmedian(monob))
tsr = (dt+dter*jcool)*kappa/(np.log(hin[1]/zot) -
psit_calc(h_in[1]/monob, meth))
qsr = (dq+dqer*jcool)*kappa/(np.log(hin[2]/zot) -
psit_calc(h_in[2]/monob, meth))
# tolerance for u,t,q,usr,tsr,qsr
tol = np.array([0.01, 0.01, 5e-05, 0.005, 0.001, 5e-07])
it, ind = 0, np.where(spd > 0)
ii, itera = True, np.zeros(spd.shape)*np.nan
while np.any(ii):
it += 1
if it > n:
break
old = np.array([np.copy(u10n), np.copy(t10n), np.copy(q10n),
np.copy(usr), np.copy(tsr), np.copy(qsr)])
cdn[ind] = cdn_calc(u10n[ind], Ta[ind], None, lat[ind], meth)
if (np.all(np.isnan(cdn))):
break
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(h_in[0, ind]/monob[ind], meth)
cd[ind] = cd_calc(cdn[ind], h_in[0, ind], ref_ht, psim[ind])
ctn[ind], cqn[ind] = ctcqn_calc(h_in[1, ind]/monob[ind], cdn[ind],
u10n[ind], zo[ind], Ta[ind], meth)
psit[ind] = psit_calc(h_in[1, ind]/monob[ind], meth)
psiq[ind] = psit_calc(h_in[2, ind]/monob[ind], meth)
ct[ind], cq[ind] = ctcq_calc(cdn[ind], cd[ind], ctn[ind], cqn[ind],
h_in[1, ind], h_in[2, ind], ref_ht,
psit[ind], psiq[ind])
if (meth == "UA"):
usr[ind] = np.where(h_in[0, ind]/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(-h_in[0, ind]/monob[ind], 1/3) -
np.power(1.574, 1/3))),
np.where((h_in[0, ind]/monob[ind] > -1.574) &
(h_in[0, ind]/monob[ind] < 0),
kappa*wind[ind]/(np.log(h_in[0, ind]/zo[ind]) -
psim_calc(h_in[0, ind]/monob[ind], meth) +
psim_calc(zo[ind]/monob[ind], meth)),
np.where((h_in[0, ind]/monob[ind] > 0) &
(h_in[0, ind]/monob[ind] < 1),
kappa*wind[ind]/(np.log(h_in[0, ind]/zo[ind]) +
5*h_in[0, ind]/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(h_in[0, ind]/monob[ind]) +
h_in[0, ind]/monob[ind]-1))))
# Zeng et al. 1998 (7-10)
tsr[ind] = np.where(h_in[1, ind]/monob[ind] < -0.465,
kappa*(dt[ind] +
dter[ind]*jcool) /
(np.log((-0.465*monob[ind])/zot[ind]) -
psit_calc(-0.465, meth)+0.8 *
(np.power(0.465, -1/3) -
np.power(-h_in[1, ind]/monob[ind], -1/3))),
np.where((h_in[1, ind]/monob[ind] > -0.465) &
(h_in[1, ind]/monob[ind] < 0),
kappa*(dt[ind]+dter[ind]*jcool) /
(np.log(h_in[1, ind]/zot[ind]) -
psit_calc(h_in[1, ind]/monob[ind], meth) +
psit_calc(zot[ind]/monob[ind], meth)),
np.where((h_in[1, ind]/monob[ind] > 0) &
(h_in[1, ind]/monob[ind] < 1),
kappa*(dt[ind]+dter[ind]*jcool) /
(np.log(h_in[1, ind]/zot[ind]) +
5*h_in[1, ind]/monob[ind]-5*zot[ind] /
monob[ind]),
kappa*(dt[ind]+dter[ind]*jcool) /
(np.log(monob[ind]/zot[ind])+5 -
5*zot[ind]/monob[ind] +
5*np.log(h_in[1, ind]/monob[ind]) +
h_in[1, ind]/monob[ind]-1))))
# Zeng et al. 1998 (11-14)
qsr[ind] = np.where(h_in[2, ind]/monob[ind] < -0.465,
kappa*(dq[ind] +
dqer[ind]*jcool) /
(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(-h_in[2, ind]/monob[ind], -1/3))),
np.where((h_in[2, ind]/monob[ind] > -0.465) &
(h_in[2, ind]/monob[ind] < 0),
kappa*(dq[ind]+dqer[ind]*jcool) /
(np.log(h_in[1, ind]/zot[ind]) -
psit_calc(h_in[2, ind]/monob[ind], meth) +
psit_calc(zoq[ind]/monob[ind], meth)),
np.where((h_in[2, ind]/monob[ind] > 0) &
(h_in[2, ind]/monob[ind]<1), kappa*(dq[ind] +
dqer[ind]*jcool) /
(np.log(h_in[1, ind]/zoq[ind]) +
5*h_in[2, ind]/monob[ind]-5*zoq[ind]/monob[ind]),
kappa*(dq[ind]+dqer[ind]*jcool) /
(np.log(monob[ind]/zoq[ind])+5 -
5*zoq[ind]/monob[ind] +
5*np.log(h_in[2, ind]/monob[ind]) +
h_in[2, ind]/monob[ind]-1))))
elif (meth == "C30" or meth == "C35" or meth == "C40"):
usr[ind] = (wind[ind]*kappa/(np.log(h_in[0, ind]/zo[ind]) -
psiu_26(h_in[0, ind]/monob[ind], meth)))
logging.info('method %s | dter = %s | Rnl = %s '
'| usr = %s | tsr = %s | qsr = %s', meth,
np.nanmedian(dter), np.nanmedian(Rnl),
np.nanmedian(usr), np.nanmedian(tsr),
np.nanmedian(qsr))
qsr[ind] = ((dq[ind]+dqer[ind]*jcool)*(kappa/(np.log(hin[2, ind] /
zoq[ind])-psit_26(hin[2, ind]/monob[ind]))))
tsr[ind] = ((dt[ind]+dter[ind]*jcool)*(kappa/(np.log(hin[1, ind] /
zot[ind])-psit_26(hin[1, ind]/monob[ind]))))
else:
usr[ind] = (wind[ind]*kappa/(np.log(h_in[0, ind]/zo[ind]) -
psim_calc(h_in[0, ind]/monob[ind], meth)))
tsr[ind] = ct[ind]*wind[ind]*(dt[ind]+dter[ind]*jcool)/usr[ind]
qsr[ind] = cq[ind]*wind[ind]*(dq[ind]+dqer[ind]*jcool)/usr[ind]
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])
Rnl[ind] = 0.97*(5.67e-8*np.power(SST[ind] -
dter[ind]*jcool+CtoK, 4)-Rl[ind])
t10n[ind] = (Ta[ind] +
(tsr[ind]*(np.log(h_in[1, ind]/ref_ht)-psit[ind])/kappa))
q10n[ind] = (qair[ind] +
(qsr[ind]*(np.log(h_in[2, ind]/ref_ht)-psiq[ind])/kappa))
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]*h_in[0, ind]/(T[ind]+CtoK)*(tsr[ind] +
0.61*(T[ind]+CtoK)*qsr[ind])/np.power(usr[ind], 2))
monob[ind] = h_in[0, ind]/zol[ind]
elif (meth == "ERA5"):
tsrv[ind] = tsr[ind]+0.61*t10n[ind]*qsr[ind]
Rb[ind] = ((g[ind]*h_in[0, ind]*((2*dt[ind])/(Ta[ind] +
sst[ind]-g[ind]*h_in[0, ind])+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((h_in[0, ind]+zo[ind])/zo[ind]) -
psim_calc((h_in[0, ind]+zo[ind])/monob[ind], meth) +
psim_calc(zo[ind]/monob[ind], meth)) /
(np.log((h_in[0, ind]+zo[ind])/zot[ind]) -
psit_calc((h_in[0, ind]+zo[ind])/monob[ind], meth) +
psit_calc(zot[ind]/monob[ind], meth))))
monob[ind] = h_in[0, ind]/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(h_in[0, ind]/monob[ind], meth)
psit[ind] = psit_calc(h_in[1, ind]/monob[ind], meth)
psiq[ind] = psit_calc(h_in[2, ind]/monob[ind], meth)
if (gust == 1 and 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)
elif (gust == 1 and (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))
elif (gust == 1):
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))
elif (gust == 0):
wind[ind] = np.copy(spd[ind])
if (meth == "UA"):
u10n[ind] = (wind[ind]+(usr[ind]/kappa)*(np.log(h_in[0, ind]/10) -
psim[ind]))
u10n[u10n < 0] = np.nan
elif (meth == "C30" or meth == "C35" or meth == "C40"):
u10n[ind] = ((wind[ind] + usr[ind]/kappa*(np.log(10/h_in[0, ind]) -
psiu_26(10/monob[ind], meth) +
psiu_26(h_in[0, ind]/monob[ind], meth)) +
psiu_26(10/monob[ind], meth)*usr[ind]/kappa /
(wind[ind]/spd[ind])))
u10n[u10n < 0] = np.nan
elif (meth == "ERA5"):
u10n[ind] = (spd[ind]+(usr[ind]/kappa)*(np.log(h_in[0, ind] /
ref_ht)-psim[ind]))
u10n[u10n < 0] = np.nan
else:
u10n[ind] = (wind[ind]+(usr[ind]/kappa)*(np.log(h_in[0, ind]/10) -
psim[ind]))
u10n[u10n < 0] = np.nan
itera[ind] = np.ones(1)*it
new = np.array([np.copy(u10n), np.copy(t10n), np.copy(q10n),
np.copy(usr), np.copy(tsr), np.copy(qsr)])
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]))
if (ind[0].size == 0):
ii = False
else:
ii = True
logging.info('method %s | # of iterations:%s', meth, it)
logging.info('method %s | # of points that did not converge :%s', meth,
ind[0].size)
# calculate output parameters
rho = (0.34838*P)/(tv10n)
t10n = t10n-(273.16+tlapse*ref_ht)
sensible = -rho*cp*usr*tsr
latent = -rho*lv*usr*qsr
if (gust == 1):
tau = rho*np.power(usr, 2)*(spd/wind)
elif (gust == 0):
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(h_in[0]/h_out[0])-psim +
psim_calc(h_out[0]/monob, meth)))
trefs = (Ta-(tsr/kappa)*(np.log(h_in[1]/h_out[1])-psit +
psit_calc(h_out[0]/monob, meth)))
trefs = trefs-(273.16+tlapse*h_out[1])
qrefs = (qair-(qsr/kappa)*(np.log(h_in[2]/h_out[2]) -
psit+psit_calc(h_out[2]/monob, meth)))
res = np.zeros((28, 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][:] = psiq
res[17][:] = u10n
res[18][:] = t10n
res[19][:] = tv10n
res[20][:] = q10n
res[21][:] = zo
res[22][:] = zot
res[23][:] = zoq
res[24][:] = urefs
res[25][:] = trefs
res[26][:] = qrefs
res[27][:] = itera
return res, ind