Skip to content

GitLab

Commit 65e4f22c authored by sbiri's avatar sbiri
Browse files

Update AirSeaFluxCode.py, flux_subs.py files

parent 0deb28e4
Hide whitespace changes
Inline Side-by-side
Showing with 25 additions and 20 deletions
AirSeaFluxCode.py
View file @ 65e4f22c
......@@ -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
flux_subs.py
View file @ 65e4f22c
......@@ -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
# ---------------------------------------------------------------------
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