From d7915eece590d68807fea2ce5f11ec2837a554b6 Mon Sep 17 00:00:00 2001
From: sbiri <sbiri@noc.ac.uk>
Date: Wed, 17 Mar 2021 09:47:53 +0000
Subject: [PATCH] changed get_L to always output Rib following ecmwf's eq.
3.24. changed flag "l" in AirSeaFluxCode to be set when Rib<-0.5 or Rib>0.2
---
AirSeaFluxCode.py | 16 ++++++++--------
flux_subs.py | 22 ++++++++++------------
2 files changed, 18 insertions(+), 20 deletions(-)
diff --git a/AirSeaFluxCode.py b/AirSeaFluxCode.py
index 1c6287f..046061d 100644
--- a/AirSeaFluxCode.py
+++ b/AirSeaFluxCode.py
@@ -122,13 +122,13 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
34. downward longwave radiation (Rl)
35. downward shortwave radiation (Rs)
36. downward net longwave radiation (Rnl)
- 37. flag ("n": normal, "o": out of nominal range,
+ 37. gust wind speed (ug)
+ 38. Bulk Richardson number (Rib)
+ 39. relative humidity (rh)
+ 40. flag ("n": normal, "o": out of nominal range,
"u": u10n<0, "q":q10n<0
- "m": missing, "l": z/L<0.01,
+ "m": missing, "l": Rib<-0.5 or Rib>0.2,
"i": convergence fail at n)
- 38. gust wind speed (ug)
- 39. Bulk Richardson number (Rib)
- 40. relative humidity (rh)
2021 / Author S. Biri
"""
@@ -349,7 +349,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")):
+ elif (gust[0] == 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(gust[1], Ta[ind], usr[ind],
tsrv[ind], gust[2], lat[ind]), 2))
@@ -427,8 +427,8 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
flag = np.where((q10n < 0) & (flag == "n"), "q",
np.where((q10n < 0) & (flag != "n"), flag+[","]+["q"],
flag))
- flag = np.where((np.abs(hin[0]/monob) < 0.01) & (flag == "n"), "l",
- np.where((np.abs(hin[0]/monob) < 0.01) & (flag != "n"),
+ flag = np.where((Rb < -0.5 | Rb > 0.2) & (flag == "n"), "l",
+ np.where((Rb < -0.5 | Rb > 0.2) & (flag != "n"),
flag+[","]+["l"], flag))
flag = np.where((itera == -1) & (flag == "n"), "i",
np.where((itera == -1) & (flag != "n"), flag+[","]+["i"],
diff --git a/flux_subs.py b/flux_subs.py
index 3175abd..2f6f740 100755
--- a/flux_subs.py
+++ b/flux_subs.py
@@ -947,25 +947,23 @@ 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
+ 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
+ dthv = (Ta-sst)*(1+0.6077*qair)+0.6077*Ta*(qair-qsea)
+ 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)))
+ Rb = g*dthv*hin[1]/(tv*uz*uz)
if (L == "S80"):
- # as in NCAR, LY04
- tsrv = tsr*(1+0.6077*qair)+0.6077*Ta*qsr
temp = (g*kappa*tsrv /
np.maximum(np.power(usr, 2)*Ta*(1+0.6077*qair), 1e-9))
temp = np.minimum(np.abs(temp), 200)*np.sign(temp)
temp = np.minimum(np.abs(temp), 10/hin[0])*np.sign(temp)
monob = 1/np.copy(temp)
elif (L == "ecmwf"):
- Rb = np.empty(sst.shape)
- 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
- dthv = (Ta-sst)*(1+0.6077*qair)+0.6077*Ta*(qair-qsea)
- 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)))
- Rb = g*dthv*hin[1]/(tv*uz*uz)
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) /
--
GitLab