Skip to content

GitLab

Commit 01d7c95d authored by sbiri's avatar sbiri
Browse files

LY04 updated to aerobulk's version

parent ea7a0324
Hide whitespace changes
Inline Side-by-side
Showing with 27 additions and 13 deletions
AirSeaFluxCode.py
View file @ 01d7c95d
......@@ -189,8 +189,6 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
u10n = np.copy(spd)
usr = 0.035*u10n
cd10n = cdn_calc(u10n, usr, Ta, lat, meth)
ct10n, ct, cq10n, 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(cd10n, h_in[0], ref_ht, psim)
......@@ -214,8 +212,14 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
wind = np.copy(spd)
# stars and roughness lengths
usr = np.sqrt(cd*np.power(wind, 2))
zo = 0.0001*np.ones(spd.shape)
zot, zoq = 0.0001*np.ones(spd.shape), 0.0001*np.ones(spd.shape)
zo = 1e-4*np.ones(spd.shape)
zot, zoq = 1e-4*np.ones(spd.shape), 1e-4*np.ones(spd.shape)
ct10n = np.power(kappa, 2)/(np.log(h_in[0]/zo)*np.log(h_in[1]/zot))
cq10n = np.power(kappa, 2)/(np.log(h_in[0]/zo)*np.log(h_in[2]/zoq))
ct = np.power(kappa, 2)/((np.log(h_in[0]/zo)-psim) *
(np.log(h_in[1]/zot)-psit))
cq = np.power(kappa, 2)/((np.log(h_in[0]/zo)-psim) *
(np.log(h_in[2]/zoq)-psiq))
monob = -100*np.ones(spd.shape) # Monin-Obukhov length
tsr = (dt+dter*cskin-dtwl*wl)*kappa/(np.log(h_in[1]/zot) -
psit_calc(h_in[1]/monob, meth))
......@@ -223,7 +227,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
psit_calc(h_in[2]/monob, meth))
# set-up to feed into iteration loop
it, ind = 0, np.where(spd > 0)
ii, itera = True, np.zeros(spd.shape)*np.nan
ii, itera = True, -1*np.ones(spd.shape)
tau = 0.05*np.ones(spd.shape)
sensible = 5*np.ones(spd.shape)
latent = 65*np.ones(spd.shape)
......@@ -258,6 +262,11 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
ct[ind], cq[ind] = ctcq_calc(cd10n[ind], cd[ind], ct10n[ind], cq10n[ind],
h_in[:, ind], [ref_ht, ref_ht, ref_ht],
psit[ind], psiq[ind])
if (meth == "LY04"):
cd = np.maximum(np.copy(cd), 1e-4)
ct = np.maximum(np.copy(ct), 1e-4)
cq = np.maximum(np.copy(cq), 1e-4)
zo = np.minimum(np.copy(zo), 0.0025)
usr[ind], tsr[ind], qsr[ind] = get_strs(h_in[:, ind], monob[ind],
wind[ind], zo[ind], zot[ind],
zoq[ind], dt[ind], dq[ind],
......@@ -358,7 +367,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
np.power(get_gust(gust[1], tv[ind], usr[ind],
tsrv[ind], gust[2], lat[ind]), 2)))
# Zeng et al. 1998 (20)
elif (gust[0] == 1 and (meth == "C30" or meth == "C35")): # or meth == "C40"
elif (gust[0] == 1 and (meth == "C30" or meth == "C35")):
wind[ind] = np.sqrt(np.power(np.copy(spd[ind]), 2) +
np.power(get_gust(gust[1], Ta[ind], usr[ind],
tsrv[ind], gust[2], lat[ind]), 2))
......@@ -370,6 +379,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
wind[ind] = np.copy(spd[ind])
u10n[ind] = wind[ind]-usr[ind]/kappa*(np.log(h_in[0, ind]/10) -
psim[ind])
# usr[ind]/np.sqrt(cd10n[ind])
if (it < 4): # make sure you allow small negative values convergence
u10n = np.where(u10n < 0, 0.5, u10n)
flag = np.where((u10n < 0) & (flag == "n"), "u",
......@@ -406,7 +416,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
else:
ii = True
itera[ind] = -1
# itera = np.where(itera > n, -1, itera)
itera = np.where(itera > n, -1, itera)
logging.info('method %s | # of iterations:%s', meth, it)
logging.info('method %s | # of points that did not converge :%s \n', meth,
ind[0].size)
......
flux_subs.py
View file @ 01d7c95d
......@@ -39,9 +39,12 @@ def cdn_calc(u10n, usr, Ta, lat, meth="S80"):
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 == "LY04"):
cdn = np.where(u10n >= 0.5,
(0.142+(2.7/u10n)+(u10n/13.09))*0.001,
(0.142+(2.7/0.5)+(0.5/13.09))*0.001)
cdn = np.where(u10n > 0.25, (0.142+2.7/u10n+u10n/13.09 -
3.14807e-10*np.power(u10n, 6))*1e-3,
(0.142+2.7/0.25+0.25/13.09 -
3.14807e-10*np.power(0.25, 6))*1e-3)
cdn = np.where(u10n > 33, 2.34e-3, np.copy(cdn))
cdn = np.maximum(np.copy(cdn), 0.1e-3)
else:
print("unknown method cdn: "+meth)
return cdn
......@@ -156,8 +159,9 @@ def ctcqn_calc(zol, cdn, usr, zo, Ta, meth="S80"):
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)
elif (meth == "LY04"):
cqn = 34.6*0.001*np.sqrt(cdn)
ctn = np.where(zol <= 0, 32.7*0.001*np.sqrt(cdn), 18*0.001*np.sqrt(cdn))
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)
elif (meth == "UA"):
# Zeng et al. 1998 (25)
rr = usr*zo/visc_air(Ta)
......@@ -274,7 +278,7 @@ def psim_calc(zol, meth="S80"):
"""
if (meth == "ecmwf"):
psim = psim_ecmwf(zol)
elif (meth == "C30" or meth == "C35"): # or meth == "C40"
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))
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment