diff --git a/AirSeaFluxCode.py b/AirSeaFluxCode.py
index d322648bde4331ed8275ce2fc7d83d4ff8d130c0..8f1afec3f70222605b8c585e25fb0ff4be684fd7 100644
--- a/AirSeaFluxCode.py
+++ b/AirSeaFluxCode.py
@@ -126,6 +126,8 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
37. gust wind speed (ug)
38. Bulk Richardson number (Rib)
39. relative humidity (rh)
+ 40. thickness of the viscous layer (delta)
+ 41. lv latent heat of vaporization (Jkg−1)
40. flag ("n": normal, "o": out of nominal range,
"u": u10n<0, "q":q10n<0
"m": missing, "l": Rib<-0.5 or Rib>0.2,
@@ -294,7 +296,8 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
meth)
if ((cskin == 1) and (wl == 0)):
if (skin == "C35"):
- dter[ind], dqer[ind], tkt[ind] = cs_C35(sst[ind], qsea[ind],
+ dter[ind], dqer[ind], tkt[ind] = cs_C35(np.copy(sst[ind]),
+ qsea[ind],
rho[ind], Rs[ind],
Rnl[ind],
cp[ind], lv[ind],
@@ -304,7 +307,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
elif (skin == "ecmwf"):
dter[ind] = cs_ecmwf(rho[ind], Rs[ind], Rnl[ind], cp[ind],
lv[ind], usr[ind], tsr[ind], qsr[ind],
- sst[ind], lat[ind])
+ np.copy(sst[ind]), lat[ind])
dqer[ind] = (dter[ind]*0.622*lv[ind]*qsea[ind] /
(287.1*np.power(sst[ind], 2)))
elif (skin == "Beljaars"):
@@ -313,6 +316,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
tsr[ind], qsr[ind], lat[ind],
np.copy(Qs[ind]))
dqer = dter*0.622*lv*qsea/(287.1*np.power(sst, 2))
+ skt = np.copy(sst)+dter
elif ((cskin == 1) and (wl == 1)):
if (skin == "C35"):
dter[ind], dqer[ind], tkt[ind] = cs_C35(sst[ind], qsea[ind],
@@ -520,7 +524,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
flag+[","]+["o"], flag))
- res = np.zeros((40, len(spd)))
+ res = np.zeros((41, len(spd)))
res[0][:] = tau
res[1][:] = sensible
res[2][:] = latent
@@ -561,14 +565,15 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
res[37][:] = Rb
res[38][:] = rh
res[39][:] = tkt
+ res[40][:] = lv
if (out == 0):
res[:, ind] = np.nan
# set missing values where data have non acceptable values
res = np.asarray([np.where(q10n < 0, np.nan,
- res[i][:]) for i in range(40)])
+ res[i][:]) for i in range(41)])
res = np.asarray([np.where(u10n < 0, np.nan,
- res[i][:]) for i in range(40)])
+ res[i][:]) for i in range(41)])
# output with pandas
resAll = pd.DataFrame(data=res.T, index=range(len(spd)),
columns=["tau", "shf", "lhf", "L", "cd", "cdn", "ct",
@@ -577,7 +582,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
"t10n", "tv10n", "q10n", "zo", "zot", "zoq",
"uref", "tref", "qref", "iteration", "dter",
"dqer", "dtwl", "qair", "qsea", "Rl", "Rs",
- "Rnl", "ug", "Rib", "rh", "delta"])
+ "Rnl", "ug", "Rib", "rh", "delta", "lv"])
resAll["flag"] = flag
return resAll
diff --git a/flux_subs.py b/flux_subs.py
index 14182fdd965f8810d1fab702ee4aa43f8f997726..e1573f6830917c1959dc0c3e6ca335e09ee38b20 100755
--- a/flux_subs.py
+++ b/flux_subs.py
@@ -687,7 +687,7 @@ def cs_ecmwf(rho, Rs, Rnl, cp, lv, usr, tsr, qsr, sst, lat):
if (np.nanmin(sst) < 200): # if sst in Celsius convert to Kelvin
sst = sst+CtoK
aw = np.maximum(1e-5, 1e-5*(sst-CtoK))
- Rns = 0.945*Rs # net solar radiation (albedo correction)
+ Rns = 0.945*Rs # (net solar radiation (albedo correction)
shf = rho*cp*usr*tsr
lhf = rho*lv*usr*qsr
Qnsol = shf+lhf+Rnl # eq. 8.152
@@ -1029,49 +1029,49 @@ def get_strs(hin, monob, wind, zo, zot, zoq, dt, dq, dter, dqer, dtwl, ct, cq,
5*np.log(hin[0]/monob) +
hin[0]/monob-1))))
# Zeng et al. 1998 (7-10)
- tsr = np.where(hin[1]/monob < -0.465, kappa*(dt+dter*cskin-dtwl*wl) /
+ tsr = np.where(hin[1]/monob < -0.465, kappa*(dt-dter*cskin-dtwl*wl) /
(np.log((-0.465*monob)/zot) -
psit_calc(-0.465, meth)+0.8*(np.power(0.465, -1/3) -
np.power(-hin[1]/monob, -1/3))),
- np.where(hin[1]/monob < 0, kappa*(dt+dter*cskin-dtwl*wl) /
+ np.where(hin[1]/monob < 0, kappa*(dt-dter*cskin-dtwl*wl) /
(np.log(hin[1]/zot) -
psit_calc(hin[1]/monob, meth) +
psit_calc(zot/monob, meth)),
np.where(hin[1]/monob <= 1,
- kappa*(dt+dter*cskin-dtwl*wl) /
+ kappa*(dt-dter*cskin-dtwl*wl) /
(np.log(hin[1]/zot) +
5*hin[1]/monob-5*zot/monob),
- kappa*(dt+dter*cskin-dtwl*wl) /
+ kappa*(dt-dter*cskin-dtwl*wl) /
(np.log(monob/zot)+5 -
5*zot/monob+5*np.log(hin[1]/monob) +
hin[1]/monob-1))))
# Zeng et al. 1998 (11-14)
- qsr = np.where(hin[2]/monob < -0.465, kappa*(dq+dqer*cskin) /
+ qsr = np.where(hin[2]/monob < -0.465, kappa*(dq-dqer*cskin) /
(np.log((-0.465*monob)/zoq) -
psit_calc(-0.465, meth)+psit_calc(zoq/monob, meth) +
0.8*(np.power(0.465, -1/3) -
np.power(-hin[2]/monob, -1/3))),
np.where(hin[2]/monob < 0,
- kappa*(dq+dqer*cskin)/(np.log(hin[1]/zot) -
+ kappa*(dq-dqer*cskin)/(np.log(hin[1]/zot) -
psit_calc(hin[2]/monob, meth) +
psit_calc(zoq/monob, meth)),
np.where(hin[2]/monob <= 1,
- kappa*(dq+dqer*cskin) /
+ kappa*(dq-dqer*cskin) /
(np.log(hin[1]/zoq)+5*hin[2]/monob -
5*zoq/monob),
- kappa*(dq+dqer*cskin)/
+ kappa*(dq-dqer*cskin)/
(np.log(monob/zoq)+5-5*zoq/monob +
5*np.log(hin[2]/monob) +
hin[2]/monob-1))))
- elif (meth == "C30" or meth == "C35"):
+ elif (meth == "C30" or meth == "C35"): # or meth == "C40"
usr = (wind*kappa/(np.log(hin[0]/zo)-psiu_26(hin[0]/monob, meth)))
- tsr = ((dt+dter*cskin-dtwl*wl)*(kappa/(np.log(hin[1]/zot) -
+ tsr = ((dt-dter*cskin-dtwl*wl)*(kappa/(np.log(hin[1]/zot) -
psit_26(hin[1]/monob))))
- qsr = ((dq+dqer*cskin)*(kappa/(np.log(hin[2]/zoq) -
+ qsr = ((dq-dqer*cskin)*(kappa/(np.log(hin[2]/zoq) -
psit_26(hin[2]/monob))))
else:
usr = (wind*kappa/(np.log(hin[0]/zo)-psim_calc(hin[0]/monob, meth)))
- tsr = ct*wind*(dt+dter*cskin-dtwl*wl)/usr
- qsr = cq*wind*(dq+dqer*cskin)/usr
+ tsr = ct*wind*(dt-dter*cskin-dtwl*wl)/usr
+ qsr = cq*wind*(dq-dqer*cskin)/usr
return usr, tsr, qsr
# ---------------------------------------------------------------------