diff --git a/AirSeaFluxCode.py b/AirSeaFluxCode.py
index d0df1d52231e9fca260c4bc42fd38a19c1efe875..de19904d60c57b49b3c2924a7fe728e786638f71 100644
--- a/AirSeaFluxCode.py
+++ b/AirSeaFluxCode.py
@@ -161,9 +161,12 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
np.copy(T)*np.power(1000/P,287.1/1004.67)) # potential T
sst = np.where(SST < 200, np.copy(SST)+CtoK, np.copy(SST))
qair, qsea = get_hum(hum, T, sst, P, qmeth)
- Rb = np.empty(sst.shape)
+ # mask to preserve missing values when initialising variables
+ msk = np.empty(sst.shape)
+ msk = np.where(np.isnan(spd+T+SST), np.nan, 1)
+ Rb = np.empty(sst.shape)*msk
# specific heat
- cp = 1004.67*(1+0.00084*qsea) #1005*np.ones(spd.shape)#
+ cp = 1004.67*(1+0.00084*qsea)
#lapse rate
tlapse = gamma("dry", SST, T, qair/1000, cp)
Ta = np.where(T < 200, np.copy(T)+CtoK+tlapse*h_in[1],
@@ -177,7 +180,6 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
rh = np.ones(sst.shape)*80
elif (hum[0] == 'rh'):
rh = hum[1]
- # rh = np.where(rh > 100, np.nan, rh)
elif (hum[0] == 'Td'):
Td = hum[1] # dew point temperature (K)
Td = np.where(Td < 200, np.copy(Td)+CtoK, np.copy(Td))
@@ -185,7 +187,6 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
esd = 611.21*np.exp(17.502*((Td-CtoK)/(Td-32.19)))
es = 611.21*np.exp(17.502*((T-CtoK)/(T-32.19)))
rh = 100*esd/es
- # rh = np.where(rh > 100, np.nan, rh)
flag = np.empty(spd.shape, dtype="object")
flag[:] = "n"
flag = np.where(np.isnan(spd+T+SST+hum[1]+P+Rs+Rl), "m", flag)
@@ -206,18 +207,18 @@ 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)
- psim, psit, psiq = (np.zeros(spd.shape), np.zeros(spd.shape),
- np.zeros(spd.shape))
+ psim, psit, psiq = (np.zeros(spd.shape)*msk, np.zeros(spd.shape)*msk,
+ np.zeros(spd.shape)*msk)
cd = cd_calc(cd10n, h_in[0], ref_ht, psim)
- tsr, tsrv = np.zeros(spd.shape), np.zeros(spd.shape)
- qsr = np.zeros(spd.shape)
+ tsr, tsrv, qsr = (np.zeros(spd.shape)*msk, np.zeros(spd.shape)*msk,
+ np.zeros(spd.shape)*msk)
# cskin parameters
- tkt = 0.001*np.ones(T.shape)
- dter = np.ones(T.shape)*0.3
+ tkt = 0.001*np.ones(T.shape)*msk
+ dter = np.ones(T.shape)*0.3*msk
dqer = dter*0.622*lv*qsea/(287.1*np.power(sst, 2))
Rnl = 0.97*(5.67e-8*np.power(sst-0.3*cskin, 4)-Rl)
Qs = 0.945*Rs
- dtwl = np.ones(T.shape)*0.3
+ dtwl = np.ones(T.shape)*0.3*msk
skt = np.copy(sst)
# gustiness adjustment
if (gust[0] == 1 and meth == "UA"):
@@ -229,25 +230,25 @@ 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 = 1e-4*np.ones(spd.shape)
- zot, zoq = 1e-4*np.ones(spd.shape), 1e-4*np.ones(spd.shape)
+ zo, zot, zoq = (1e-4*np.ones(spd.shape)*msk, 1e-4*np.ones(spd.shape)*msk,
+ 1e-4*np.ones(spd.shape)*msk)
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
+ monob = -100*np.ones(spd.shape)*msk # Monin-Obukhov length
tsr = (dt+dter*cskin-dtwl*wl)*kappa/(np.log(h_in[1]/zot) -
psit_calc(h_in[1]/monob, meth))
qsr = (dq+dqer*cskin)*kappa/(np.log(h_in[2]/zoq) -
psit_calc(h_in[2]/monob, meth))
# set-up to feed into iteration loop
it, ind = 0, np.where(spd > 0)
- 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)
+ ii, itera = True, -1*np.ones(spd.shape)*msk
+ tau = 0.05*np.ones(spd.shape)*msk
+ sensible = 5*np.ones(spd.shape)*msk
+ latent = 65*np.ones(spd.shape)*msk
# iteration loop
while np.any(ii):
it += 1
@@ -371,8 +372,8 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
qsr[ind], h_in[:, ind], Ta[ind],
sst[ind]-dter[ind]*cskin+dtwl[ind]*wl,
qair[ind], qsea[ind], wind[ind],
- np.copy(monob[ind]), psim[ind],
- meth)
+ np.copy(monob[ind]), zo[ind],
+ zot[ind], psim[ind], meth)
# sst[ind]-dter[ind]*cskin+dtwl[ind]*wl
psim[ind] = psim_calc(h_in[0, ind]/monob[ind], meth)
psit[ind] = psit_calc(h_in[1, ind]/monob[ind], meth)
@@ -396,14 +397,12 @@ 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",
np.where((u10n < 0) &
(np.char.find(flag.astype(str), 'u') == -1),
flag+[","]+["u"], flag))
- # u10n = np.where(u10n < 0, np.nan, u10n)
utmp = np.copy(u10n)
utmp = np.where(utmp < 0, np.nan, utmp)
itera[ind] = np.ones(1)*it
@@ -464,7 +463,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
qref = (qair-qsr/kappa*(np.log(h_in[2]/h_out[2]) -
psit+psit_calc(h_out[2]/monob, meth)))
if (wl == 0):
- dtwl = np.zeros(T.shape) # reset to zero if not used
+ dtwl = np.zeros(T.shape)*msk # reset to zero if not used
flag = np.where((q10n < 0) & (flag == "n"), "q",
np.where((q10n < 0) & (flag != "n"), flag+[","]+["q"],
flag))
diff --git a/flux_subs.py b/flux_subs.py
index 3482b3387c9f70898633b60fd6a0c2cbd2618d09..771cf071c17a4e08bd276c8e719249c5da4b0baf 100755
--- a/flux_subs.py
+++ b/flux_subs.py
@@ -878,8 +878,8 @@ def get_gust(beta, Ta, usr, tsrv, zi, lat):
# ---------------------------------------------------------------------
-def get_L(L, lat, usr, tsr, qsr, hin, Ta, sst, qair, qsea, wind, monob, psim,
- meth):
+def get_L(L, lat, usr, tsr, qsr, hin, Ta, sst, qair, qsea, wind, monob, zo,
+ zot, psim, meth):
"""
calculates Monin-Obukhov length and virtual star temperature
@@ -915,6 +915,10 @@ def get_L(L, lat, usr, tsr, qsr, hin, Ta, sst, qair, qsea, wind, monob, psim,
wind speed (m/s)
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", "LY04", "C30", "C35", "ecmwf", "Beljaars"
@@ -929,7 +933,7 @@ def get_L(L, lat, usr, tsr, qsr, hin, Ta, sst, qair, qsea, wind, monob, psim,
"""
g = gc(lat)
Rb = np.empty(sst.shape)
- # as in NCAR, LY04
+ # as in aerobulk One_on_L in mod_phymbl.f90
tsrv = tsr*(1+0.6077*qair)+0.6077*Ta*qsr
# from eq. 3.24 ifs Cy46r1 pp. 37
thvs = sst*(1+0.6077*qsea) # virtual SST
@@ -937,7 +941,7 @@ def get_L(L, lat, usr, tsr, qsr, hin, Ta, sst, qair, qsea, wind, monob, psim,
tv = 0.5*(thvs+Ta*(1+0.6077*qair)) # estimate tv within surface layer
# adjust wind to T sensor's height
uz = (wind-usr/kappa*(np.log(hin[0]/hin[1])-psim +
- psim_calc(hin[1]/monob, meth)))
+ psim_calc(hin[1]/monob, meth)))
Rb = g*dthv*hin[1]/(tv*uz*uz)
if (L == "S80"):
temp = (g*kappa*tsrv /
@@ -946,8 +950,6 @@ def get_L(L, lat, usr, tsr, qsr, hin, Ta, sst, qair, qsea, wind, monob, psim,
temp = np.minimum(np.abs(temp), 10/hin[0])*np.sign(temp)
monob = 1/np.copy(temp)
elif (L == "ecmwf"):
- zo = (0.11*visc_air(Ta)/usr+0.018*np.power(usr, 2)/g)
- zot = 0.40*visc_air(Ta)/usr
zol = (Rb*(np.power(np.log((hin[1]+zo)/zo)-psim_calc((hin[1]+zo) /
monob, meth) +
psim_calc(zo/monob, meth), 2) /