diff --git a/AirSeaFluxCode.py b/AirSeaFluxCode.py
index f3c1f7662b82052267aa0b2e7f9dda9926110b44..1862e6341da9f7509017396010a7104faff8761e 100644
--- a/AirSeaFluxCode.py
+++ b/AirSeaFluxCode.py
@@ -3,18 +3,23 @@ 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, q_calc, qsea_calc, qsat26sea, qsat26air,
- visc_air, psit_26, psiu_26)
+ gc, qsat_sea, qsat_air, 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):
+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)
@@ -29,8 +34,7 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi=600,
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
+ sensor heights in m (array 3x1 or 3xn) default 10m ([10, 10, 10])
hout : float
output height, default is 10m
zi : int
@@ -41,6 +45,8 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi=600,
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
@@ -62,18 +68,19 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi=600,
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
+ 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
@@ -81,15 +88,12 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi=600,
"""
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
+ 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
- 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)
- cdn = 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")
@@ -121,34 +125,28 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi=600,
elif ((np.all(np.isnan(zi))) and (meth == "ERA5" or meth == "UA")):
zi = 1000
####
- th = np.where(np.nanmax(T) < 200, (np.copy(T)+CtoK) *
+ 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(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)]))
+ 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 = qsea_calc(sst, P)
- qair = q_calc(Ta, RH, P)
+ 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.ma.median(qsea[~np.isnan(qsea)]),
- np.ma.median(qair[~np.isnan(qair)]))
+ 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.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)]))
+ '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
@@ -159,70 +157,66 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi=600,
dtv=dt*(1.+0.61*qair)+0.61*th*dq
# ------------
rho = P*100/(287.1*tv10n)
-# 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*np.ones(Ta.shape) # Zeng et al. 1998, C3.0, C3.5
u10n = np.copy(spd)
- monob = -100*np.ones(u10n.shape) # Monin-Obukhov length
+ monob = -100*np.ones(spd.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 jcool == 1:
+ 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 (meth == "UA"):
+ 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(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) /
+ 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(hh_in[0]/zo))
- monob = hh_in[0]/zol
+ 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.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)]))
+ np.nanmedian(wind), np.nanmedian(usr), np.nanmedian(zo),
+ np.nanmedian(zot), np.nanmedian(zoq), np.nanmedian(Rb),
+ np.nanmedian(monob))
elif (meth == "ERA5"):
- wind = np.sqrt(np.power(np.copy(spd), 2)+np.power(0.5, 2))
usr = np.sqrt(cd*np.power(wind, 2))
- Rb = ((g*hh_in[0]*((2*dt)/(Ta+sst-g*hh_in[0]) +
+ 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((hh_in[0]+zo)/zo)-psim_calc((hh_in[0]+zo) /
+ zol = (Rb*((np.log((h_in[0]+zo)/zo)-psim_calc((h_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) +
+ (np.log((h_in[0]+zo)/zot) -
+ psit_calc((h_in[0]+zo)/monob, meth) +
psit_calc(zot/monob, meth))))
- monob = hh_in[0]/zol
+ monob = h_in[0]/zol
logging.info('method %s | wind:%s, usr:%s, '
'zo:%s, zot:%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)]))
+ 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"):
- wind = np.sqrt(np.power(np.copy(spd), 2)+np.power(0.5, 2))
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),
@@ -231,38 +225,30 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi=600,
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) /
+ 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((hh_in[0]+zo)/zot) -
- psit_calc((hh_in[0]+zo)/monob, meth) +
+ (np.log((h_in[0]+zo)/zot) -
+ psit_calc((h_in[0]+zo)/monob, meth) +
psit_calc(zot/monob, meth))))
- monob = hh_in[0]/zol
+ monob = h_in[0]/zol
logging.info('method %s | wind:%s, usr:%s, '
'zo:%s, zot:%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)]))
+ np.nanmedian(wind), np.nanmedian(usr), np.nanmedian(zo),
+ np.nanmedian(zot), np.nanmedian(Rb), np.nanmedian(monob))
else:
- wind = np.copy(spd)
- zo, zot = 0.0001*np.ones(u10n.shape), 0.0001*np.ones(u10n.shape)
+ 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, 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)]))
+ '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(hh_in[1]/monob, meth))
+ psit_calc(h_in[1]/monob, meth))
qsr = (dq+dqer*jcool)*kappa/(np.log(hin[2]/zot) -
- psit_calc(hh_in[2]/monob, meth))
+ 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]) #[1e-06, 0.01, 5e-07, 1e-06, 0.001, 5e-07]
+ 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):
@@ -276,180 +262,183 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi=600,
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(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],
+ 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(hh_in[1]/monob[ind], meth)
- psiq[ind] = psit_calc(hh_in[2]/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)
ct[ind], cq[ind] = ctcq_calc(cdn[ind], cd[ind], ctn[ind], cqn[ind],
- hh_in[1], hh_in[2], ref_ht,
+ h_in[1, ind], h_in[2, ind], ref_ht,
psit[ind], psiq[ind])
if (meth == "UA"):
- usr[ind] = np.where(hh_in[0]/monob[ind] < -1.574, kappa*wind[ind] /
+ 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(-hh_in[0]/monob[ind], 1/3) -
+ 1.14*(np.power(-h_in[0, ind]/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]) +
+ 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, 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(hh_in[0]/monob[ind]) +
- hh_in[0]/monob[ind]-1))))
+ 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(hh_in[1]/monob[ind] < -0.465, kappa*(dt[ind] +
+ 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(-hh_in[1]/monob[ind], -1/3))),
- np.where((hh_in[1]/monob[ind] > -0.465) &
- (hh_in[1]/monob[ind] < 0),
+ 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(hh_in[1]/zot[ind]) -
- psit_calc(hh_in[1]/monob[ind], meth) +
+ (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((hh_in[1]/monob[ind] > 0) &
- (hh_in[1]/monob[ind] < 1),
+ np.where((h_in[1, ind]/monob[ind] > 0) &
+ (h_in[1, ind]/monob[ind] < 1),
kappa*(dt[ind]+dter[ind]*jcool) /
- (np.log(hh_in[1]/zot[ind]) +
- 5*hh_in[1]/monob[ind]-5*zot[ind]/monob[ind]),
+ (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(hh_in[1]/monob[ind]) +
- hh_in[1]/monob[ind]-1))))
+ 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(hh_in[2]/monob[ind] < -0.465, kappa*(dq[ind] +
+ 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(-hh_in[2]/monob[ind], -1/3))),
- np.where((hh_in[2]/monob[ind] > -0.465) &
- (hh_in[2]/monob[ind] < 0),
+ 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(hh_in[1]/zot[ind]) -
- psit_calc(hh_in[2]/monob[ind], meth) +
+ (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((hh_in[2]/monob[ind] > 0) &
- (hh_in[2]/monob[ind]<1), kappa*(dq[ind] +
+ np.where((h_in[2, ind]/monob[ind] > 0) &
+ (h_in[2, ind]/monob[ind]<1), kappa*(dq[ind] +
dqer[ind]*jcool) /
- (np.log(hh_in[1]/zoq[ind]) +
- 5*hh_in[2]/monob[ind]-5*zoq[ind]/monob[ind]),
+ (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(hh_in[2]/monob[ind]) +
- hh_in[2]/monob[ind]-1))))
+ 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(hh_in[0]/zo[ind]) -
- psiu_26(hh_in[0]/monob[ind], meth)))
+ 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.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]))))
+ 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] = np.sqrt(cd[ind]*np.power(wind[ind], 2))
- # = (wind[ind]*kappa/(np.log(hh_in[0]/zo[ind]) -
- # psim_calc(hh_in[0]/monob[ind], meth)))
+ 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]
- if (jcool == 1):
- 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])
- else:
- dter[ind] = np.zeros(spd[ind].shape)
- dqer[ind] = np.zeros(spd[ind].shape)
- tkt[ind] = np.zeros(spd[ind].shape)
+ 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(hh_in[1]/ref_ht)-psit[ind])/kappa))
- q10n[ind] = (qair[ind] -
- (qsr[ind]*(np.log(hh_in[2]/ref_ht)-psiq[ind])/kappa))
+ 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]*hh_in[0]/(T[ind]+CtoK)*(tsr[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] = hh_in[0]/zol[ind]
+ 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]*hh_in[0]*((2*dt[ind])/(Ta[ind]+sst[ind]-g[ind] *
- hh_in[0])+0.61*dq[ind]))/np.power(wind[ind], 2))
+ 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((hh_in[0]+zo[ind])/zo[ind]) -
- psim_calc((hh_in[0]+zo[ind])/monob[ind], meth) +
+ 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((hh_in[0]+zo[ind])/zot[ind]) -
- psit_calc((hh_in[0]+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] = hh_in[0]/zol[ind]
+ 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(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"):
+ 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)
- 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"):
+ 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))
- u10n[ind] = ((wind[ind] + usr[ind]/kappa*(np.log(10/hh_in[0])-
+ 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(hh_in[0]/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
+ # 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] /
+ u10n[ind] = (spd[ind]+(usr[ind]/kappa)*(np.log(h_in[0, ind] /
ref_ht)-psim[ind]))
- u10n[u10n < 0] = np.nan
+ # u10n[u10n < 0] = np.nan
else:
- u10n[ind] = (wind[ind]+(usr[ind]/kappa)*(np.log(hh_in[0]/10) -
+ u10n[ind] = (wind[ind]+(usr[ind]/kappa)*(np.log(h_in[0, ind]/10) -
psim[ind]))
- u10n[u10n < 0] = np.nan
+ # 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)])
@@ -461,41 +450,29 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi=600,
ii = False
else:
ii = True
- # if ((it == 3) or (it == 6) or (it == 10)):
- # print('Method {}, # {} tol u10n: {} | t10n: {} | q10n: {} | '
- # 'u*: {} | t*: {} | q*: {}'.format(meth, it,
- # np.ma.max(d[0, ~np.isnan(d[0,:])]),
- # np.ma.max(d[1, ~np.isnan(d[1,:])]),
- # np.ma.max(d[2, ~np.isnan(d[2,:])]),
- # np.ma.max(d[3, ~np.isnan(d[3,:])]),
- # np.ma.max(d[4, ~np.isnan(d[4,:])]),
- # np.ma.max(d[5, ~np.isnan(d[5,:])])),
- # file=open('tol_mid.txt','a'))
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)
-# rho = P*100./(287.1*(T+CtoK)*(1+0.61*qair)) # C35
t10n = t10n-(273.16+tlapse*ref_ht)
sensible = -rho*cp*usr*tsr
latent = -rho*lv*usr*qsr
- if (meth == "C30" or meth == "C35" or meth == "C40" or meth == "UA" or
- meth == "ERA5"):
+ if (gust == 1):
tau = rho*np.power(usr, 2)*(spd/wind)
- else:
+ 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(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)))
+ 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
@@ -512,15 +489,16 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi=600,
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
+ 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