Skip to content

GitLab

Commit c0ee1bf6 authored by sbiri's avatar sbiri
Browse files

changed gustiness application; added function for output variables in... 

changed gustiness application; added function for output variables in util_subs; changed get_gust in flux_subs to input also stable ug
parent 06469c60
Hide whitespace changes
Inline Side-by-side
Showing with 135 additions and 95 deletions
Code/AirSeaFluxCode.py
View file @ c0ee1bf6
......@@ -12,14 +12,12 @@ class S88:
def _wind_iterate(self, ind):
if self.gust[0] in range(1, 6):
self.wind[ind] = np.sqrt(np.power(np.copy(self.spd[ind]), 2) +
np.power(get_gust(self.gust[1],
self.theta[ind],
self.usr[ind],
self.tsrv[ind],
self.gust[2],
self.grav[ind]), 2))
# ratio of gusty to horizontal wind
self.GustFact = apply_GF(self.gust, self.spd, self.wind, "u")
np.power(get_gust(
self.gust[1], self.gust[2],
self.gust[3], self.theta[ind],
self.usr[ind], self.tsrv[ind],
self.grav[ind]), 2))
self.GustFact = np.sqrt(self.wind/self.spd)
self.u10n[ind] = self.spd[ind]-self.usr[ind]/kappa / \
self.GustFact[ind]*(np.log(self.h_in[0, ind]/self.ref10) -
self.psim[ind])
......@@ -140,7 +138,6 @@ class S88:
Rb = self.grav*self.h_in[1]*self.dtv/(self.T*np.power(self.wind, 2))
# eq. 12 Grachev & Fairall 1997 # DO.THIS
self.monob = self.h_in[1]/12.0/Rb
# where does 12.0 come from??
# ------------
......@@ -195,7 +192,6 @@ class S88:
"ref": ["u10n", "t10n", "q10n"]}
new_vars["all"] = new_vars["ref"] + new_vars["flux"]
new_vars = new_vars[tol[0]]
# I'm sure there are better ways of doing this
# extract tolerance values by deleting flag from tol
tvals = np.delete(np.copy(tol), 0)
tol_vals = list([float(tt) for tt in tvals])
......@@ -324,14 +320,9 @@ class S88:
self.u10n = np.where(self.u10n < 0, 0.5, self.u10n)
self.itera[ind] = np.full(1, it)
# remove effect of gustiness following Fairall et al. (2003)
# usr is divided by (GustFact)^0.5 (here applied to sensible and
# latent as well as tau)
# GustFact should be 1 if gust is OFF or gust[0]=2
self.GustFact = apply_GF(self.gust, self.spd, self.wind, "TSF")
self.tau = self.rho*np.power(self.usr, 2)/self.GustFact[0]
self.sensible = self.rho*self.cp*(self.usr/self.GustFact[1])*self.tsr
self.latent = self.rho*self.lv*(self.usr/self.GustFact[2])*self.qsr
self.tau = self.rho*np.power(self.usr, 2)
self.sensible = self.rho*self.cp*self.usr*self.tsr
self.latent = self.rho*self.lv*self.usr*self.qsr
# Set the new variables (for comparison against "old")
new = np.array([np.copy(getattr(self, i)) for i in new_vars])
......@@ -388,9 +379,15 @@ class S88:
self._get_humidity() # Get the Relative humidity
self._flag(out=out) # Get flags
self.GustFact = apply_GF(self.gust, self.spd, self.wind, "TSF")
self.tau = self.rho*np.power(self.usr, 2)/self.GustFact[0]
self.sensible = self.rho*self.cp*(self.usr/self.GustFact[1])*self.tsr
self.latent = self.rho*self.lv*(self.usr/self.GustFact[2])*self.qsr
self.GustFact = apply_GF(self.gust, self.spd, self.wind, "u")
# remove effect of gustiness following Fairall et al. (2003)
# usr is divided by (GustFact)^0.5
self.u10n = self.spd-self.usr/kappa/self.GustFact*(
np.log(self.h_in[0]/self.ref10)-self.psim) # C.4-7
self.usrGF = self.usr/self.GustFact
self.uref = self.spd-self.usr/kappa/self.GustFact * \
(np.log(self.h_in[0]/self.h_out[0])-self.psim +
psim_calc(self.h_out[0]/self.monob, self.meth))
......@@ -401,13 +398,16 @@ class S88:
self.qref = self.qair-self.qsr/kappa * \
(np.log(self.h_in[2]/self.h_out[2])-self.psiq +
psit_calc(self.h_out[2]/self.monob, self.meth))
self.psim_ref = psim_calc(self.h_out[0]/self.monob, self.meth)
self.psit_ref = psit_calc(self.h_out[1]/self.monob, self.meth)
self.psiq_ref = psit_calc(self.h_out[2]/self.monob, self.meth)
if self.wl == 0:
self.dtwl = np.zeros(self.T.shape)*self.msk
# reset to zero if not used
# Do not calculate lhf if a measure of humidity is not input
# This gets filled into a pd dataframe and so no need to specify y
# This gets filled into a pd dataframe and so no need to specify
# dimension of array
if self.hum[0] == 'no':
self.latent, self.qsr, self.q10n = np.empty(3)
......@@ -425,15 +425,8 @@ class S88:
pass
# Combine all output variables into a pandas array
if out_var == "limited":
out_var = ("tau", "sensible", "latent", "uref", "tref", "qref")
res_vars = ("tau", "sensible", "latent", "monob", "cd", "cd10n", "ct",
"ct10n", "cq", "cq10n", "tsrv", "tsr", "qsr", "usr",
"psim", "psit", "psiq", "u10n", "t10n", "q10n", "zo",
"zot", "zoq", "uref", "tref", "qref", "dter", "dqer",
"dtwl", "qair", "qsea", "Rl", "Rs", "Rnl", "ug", "Rb",
"rh", "rho", "cp", "tkt", "lv",
"itera") if out_var is None else out_var
res_vars = get_outvars(out_var, self.cskin, self.gust)
res = np.zeros((len(res_vars), len(self.spd)))
for i, value in enumerate(res_vars):
......@@ -443,14 +436,18 @@ class S88:
res[:, self.ind] = np.nan
# set missing values where data have non acceptable values
if self.hum[0] != 'no':
res[:-1] = np.asarray([
res = np.asarray([
np.where(self.q10n < 0, np.nan, res[i][:])
for i in range(len(res_vars)-1)])
for i in range(len(res_vars))])
# len(res_vars)-1 instead of len(res_vars) in order to keep
# itera= -1 for no convergence
res[:-1] = np.asarray([
res = np.asarray([
np.where(self.u10n < 0, np.nan, res[i][:])
for i in range(len(res_vars)-1)])
for i in range(len(res_vars))])
res = np.asarray([
np.where(((self.t10n < 173) | (self.t10n > 373)), np.nan,
res[i][:])
for i in range(len(res_vars))])
else:
warnings.warn("Warning: the output will contain values for points"
" that have not converged and negative values "
......@@ -458,6 +455,8 @@ class S88:
resAll = pd.DataFrame(data=res.T, index=range(self.nlen),
columns=res_vars)
if "itera" in res_vars:
resAll["itera"] = self.itera # restore itera
resAll["flag"] = self.flag
......@@ -490,11 +489,12 @@ class S88:
try:
gust = self.default_gust
except AttributeError:
gust = [0, 0, 0] # gustiness OFF
gust = [0, 0, 0, 0] # gustiness OFF
# gust = [1, 1.2, 800]
elif ((np.size(gust) < 3) and (gust == 0)):
gust = [0, 0, 0]
gust = [0, 0, 0, 0]
assert np.size(gust) == 3, "gust input must be a 3x1 array"
assert np.size(gust) == 4, "gust input must be a 4x1 array"
assert gust[0] in range(6), "gust at position 0 must be 0 to 5"
self.gust = gust
......@@ -559,7 +559,7 @@ class UA(S88):
def __init__(self):
self.meth = "UA"
self.default_gust = [3, 1, 1000]
self.default_gust = [1, 1.2, 600]
self.u_lo = [-999, -999]
self.u_hi = [18, 18]
......@@ -570,14 +570,14 @@ class C30(S88):
def __init__(self):
self.meth = "C30"
self.default_gust = [3, 1.2, 600]
self.default_gust = [1, 1.2, 600]
self.skin = "C35"
class C35(C30):
def __init__(self):
self.meth = "C35"
self.default_gust = [3, 1.2, 600]
self.default_gust = [1, 1.2, 600]
self.skin = "C35"
......@@ -593,7 +593,7 @@ class ecmwf(C30):
def __init__(self):
self.meth = "ecmwf"
self.default_gust = [3, 1, 1000]
self.default_gust = [1, 1.2, 600]
self.skin = "ecmwf"
......@@ -604,7 +604,7 @@ class Beljaars(C30):
def __init__(self):
self.meth = "Beljaars"
self.default_gust = [3, 1, 1000]
self.default_gust = [1, 1.2, 600]
self.skin = "Beljaars"
......@@ -654,13 +654,14 @@ def AirSeaFluxCode(spd, T, SST, meth, lat=None, hum=None, P=None, hin=18,
wl : int
warm layer correction default is 0, to switch on set to 1
gust : int
3x1 [x, beta, zi] x=0 gustiness is OFF, x=1-5 gustiness is ON and
use gustiness factor: 1. Fairall et al. 2003, 2. GF is removed
4x1 [x, beta, zi, ustb] x=0 gustiness is OFF, x=1-5 gustiness is ON
and use gustiness factor: 1. Fairall et al. 2003, 2. GF is removed
from TSFs u10n, uref, 3. GF=1, 4. following Zeng et al. 1998 or
Brodeau et al. 2006, 5. following C35 matlab code;
beta gustiness parameter, beta=1 for UA, ecmwf & Beljaars,
beta=1.2 for COARE, zi PBL height (m) 600 for COARE,
1000 for UA, ecmwf & Beljaars
beta gustiness parameter, default is 1.2,
zi PBL height (m) default is 600,
min is the value for gust speed in stable conditions,
default is 0.01ms^{-1}
qmeth : str
is the saturation evaporation method to use amongst
"HylandWexler","Hardy","Preining","Wexler","GoffGratch","WMO",
......@@ -678,18 +679,19 @@ def AirSeaFluxCode(spd, T, SST, meth, lat=None, hum=None, P=None, hin=18,
number of iterations (default = 10)
out : int
set 0 to set points that have not converged, negative values of
u10n and q10n to missing (default)
u10n, q10n or T10n out of limits to missing (default)
set 1 to keep points
out_var : str
optional. user can define pandas array of variables to be output.
the default full pandas array is :
out_var = ("tau", "sensible", "latent", "monob", "cd", "cd10n",
"ct", "ct10n", "cq", "cq10n", "tsrv", "tsr", "qsr",
"usr", "psim", "psit", "psiq", "u10n", "t10n",
"q10n", "zo", "zot", "zoq", "uref", "tref", "qref",
"dter", "dqer", "dtwl", "qair", "qsea", "Rl", "Rs",
"Rnl", "ug", "Rb", "rh", "rho", "cp", "tkt", "lv",
"itera")
"ct", "ct10n", "cq", "cq10n", "tsrv", "tsr", "qsr",
"usr", "psim", "psit", "psiq", "psim_ref", "psit_ref",
"psiq_ref", "u10n", "t10n", "q10n", "zo", "zot", "zoq",
"uref", "tref", "qref", "dter", "dqer", "dtwl", "tkt",
"qair", "qsea", "Rl", "Rs", "Rnl", "ug", "usrGF",
"GustFact", "Rb", "rh", "rho", "cp", "lv", "theta",
"itera")
the "limited" pandas array is:
out_var = ("tau", "sensible", "latent", "uref", "tref", "qref")
the user can define a custom pandas array of variables to output.
......@@ -718,32 +720,38 @@ def AirSeaFluxCode(spd, T, SST, meth, lat=None, hum=None, P=None, hin=18,
15. momentum stability function (psim)
16. heat stability function (psit)
17. moisture stability function (psiq)
18. 10m neutral wind speed (u10n)
19. 10m neutral temperature (t10n)
20. 10m neutral specific humidity (q10n)
21. surface roughness length (zo)
22. heat roughness length (zot)
23. moisture roughness length (zoq)
24. wind speed at reference height (uref)
25. temperature at reference height (tref)
26. specific humidity at reference height (qref)
27. cool-skin temperature depression (dter)
28. cool-skin humidity depression (dqer)
29. warm layer correction (dtwl)
30. specific humidity of air (qair)
31. specific humidity at sea surface (qsea)
32. downward longwave radiation (Rl)
33. downward shortwave radiation (Rs)
34. downward net longwave radiation (Rnl)
35. gust wind speed (ug)
36. Bulk Richardson number (Rib)
37. relative humidity (rh)
38. air density (rho)
39 specific heat of moist air (cp)
40. thickness of the viscous layer (delta)
41. lv latent heat of vaporization (Jkg−1)
42. number of iterations until convergence
43. flag ("n": normal, "o": out of nominal range,
18. momentum stability function at hout (psim_ref)
19. heat stability function at hout (psit_ref)
20. moisture stability function at hout (psiq_ref)
21. 10m neutral wind speed (u10n)
22. 10m neutral temperature (t10n)
23. 10m neutral specific humidity (q10n)
24. surface roughness length (zo)
25. heat roughness length (zot)
26. moisture roughness length (zoq)
27. wind speed at reference height (uref)
28. temperature at reference height (tref)
29. specific humidity at reference height (qref)
30. cool-skin temperature depression (dter)
31. cool-skin humidity depression (dqer)
32. warm layer correction (dtwl)
33. thickness of the viscous layer (delta)
34. specific humidity of air (qair)
35. specific humidity at sea surface (qsea)
36. downward longwave radiation (Rl)
37. downward shortwave radiation (Rs)
38. downward net longwave radiation (Rnl)
39. gust wind speed (ug)
40. star wind speed/GustFact (usrGF)
41. Gustiness Factor (GustFact)
42. Bulk Richardson number (Rb)
43. relative humidity (rh)
44. air density (rho)
45. specific heat of moist air (cp)
46. lv latent heat of vaporization (Jkg−1)
47. potential temperature (theta)
48. number of iterations until convergence
49. flag ("n": normal, "o": out of nominal range,
"u": u10n<0, "q":q10n<0 or q>0.04
"m": missing,
"l": Rib<-0.5 or Rib>0.2 or z/L>1000,
......
Code/flux_subs.py
View file @ c0ee1bf6
......@@ -560,7 +560,7 @@ def psim_stab(zol, meth):
# ---------------------------------------------------------------------
def get_gust(beta, Ta, usr, tsrv, zi, grav):
def get_gust(beta, zi, ustb, Ta, usr, tsrv, grav):
"""
Compute gustiness.
......@@ -568,14 +568,16 @@ def get_gust(beta, Ta, usr, tsrv, zi, grav):
----------
beta : float
constant
zi : int
scale height of the boundary layer depth [m]
ustb : float
gust wind in stable conditions [m/s]
Ta : float
air temperature [K]
usr : float
friction velocity [m/s]
tsrv : float
star virtual temperature of air [K]
zi : int
scale height of the boundary layer depth [m]
grav : float
gravity
......@@ -587,8 +589,8 @@ def get_gust(beta, Ta, usr, tsrv, zi, grav):
Ta = Ta+273.16
# minus sign to allow cube root
Bf = (-grav/Ta)*usr*tsrv
ug = np.ones(np.shape(Ta))*0.2
ug = np.where(Bf > 0, beta*np.power(Bf*zi, 1/3), 0.2)
ug = np.ones(np.shape(Ta))*ustb
ug = np.where(Bf > 0, np.maximum(beta*np.power(Bf*zi, 1/3), ustb), ustb)
return ug
# ---------------------------------------------------------------------
......@@ -609,7 +611,7 @@ def apply_GF(gust, spd, wind, step):
2: gustiness is switched ON and GF is removed from TSFs u10n, uref
3: gustiness is switched ON and GF=1
4: gustiness is switched ON following Zeng et al. 1998 or
Brodeau et al. 2006 for ECMWF
Brodeau et al. 2006
5: gustiness is switched ON following C35 matlab code
spd : float
wind speed [ms^{-1}]
......@@ -637,7 +639,7 @@ def apply_GF(gust, spd, wind, step):
elif step == "TSF":
# remove effect of gustiness from TSFs
# here it is a 3xspd.shape array
GustFact = np.empty([3, spd.shape[0]], dtype=float)*np.nan
GustFact = np.ones([3, spd.shape[0]], dtype=float)
GustFact[0, :] = wind/spd
GustFact[1:3, :] = wind*0+1
# following Fairall et al. (2003)
......@@ -673,6 +675,8 @@ def get_strs(hin, monob, wind, zo, zot, zoq, dt, dq, cd, ct, cq, meth):
temperature difference [K]
dq : float
specific humidity difference [g/kg]
cd : float
drag coefficient
ct : float
temperature exchange coefficient
cq : float
......@@ -767,7 +771,6 @@ def get_tsrv(tsr, qsr, Ta, qair):
virtual star temperature (K)
"""
# as in aerobulk One_on_L in mod_phymbl.f90
tsrv = tsr*(1+0.6077*qair)+0.6077*Ta*qsr
return tsrv
......@@ -796,8 +799,6 @@ def get_Rb(grav, usr, hin_u, hin_t, tv, dtv, wind, monob, meth):
wind speed (m/s)
monob : float
Monin-Obukhov length from previous iteration step (m)
psim : float
momentum stability function
meth : str
bulk parameterisation method option: "S80", "S88", "LP82", "YT96",
"UA", "NCAR", "C30", "C35", "ecmwf", "Beljaars"
......@@ -880,10 +881,6 @@ def get_Ltsrv(tsrv, grav, tv, usr):
M-O length (m)
"""
# tmp = (g*kappa*tsrv /
# np.maximum(np.power(usr, 2)*Ta*(1+0.6077*qair), 1e-9))
# tmp = np.minimum(np.abs(tmp), 200)*np.sign(tmp)
# monob = 1/np.copy(tmp)
tsrv = np.maximum(np.abs(tsrv), 1e-9)*np.sign(tsrv)
monob = (np.power(usr, 2)*tv)/(grav*kappa*tsrv)
return monob
......
Code/util_subs.py
View file @ c0ee1bf6
......@@ -174,4 +174,39 @@ def set_flag(miss, rh, u10n, q10n, t10n, Rb, hin, monob, itera, out=0):
(np.char.find(flag.astype(str), 'u') == -1)),
flag+[","]+["i"], flag))
return flag
\ No newline at end of file
return flag
# ---------------------------------------------------------------------
def get_outvars(out_var, cskin, gust):
if out_var is None: # full output
if cskin == 1 and gust[0] == 0: # skin ON and gust OFF
res_vars = ("tau", "sensible", "latent", "monob", "cd", "cd10n",
"ct", "ct10n", "cq", "cq10n", "tsrv", "tsr", "qsr",
"usr", "psim", "psit", "psiq", "psim_ref", "psit_ref",
"psiq_ref", "u10n", "t10n", "q10n", "zo", "zot", "zoq",
"uref", "tref", "qref", "dter", "dqer", "dtwl", "tkt",
"qair", "qsea", "Rl", "Rs", "Rnl", "Rb", "rh", "rho",
"cp", "lv", "theta", "itera")
elif cskin == 0 and gust[0] != 0: # skin OFF and gust ON
res_vars = ("tau", "sensible", "latent", "monob", "cd", "cd10n",
"ct", "ct10n", "cq", "cq10n", "tsrv", "tsr", "qsr",
"usr", "psim", "psit", "psiq", "psim_ref", "psit_ref",
"psiq_ref", "u10n", "t10n", "q10n", "zo", "zot", "zoq",
"uref", "tref", "qref", "qair", "qsea", "ug", "usrGF",
"GustFact", "Rb", "rh", "rho", "cp", "lv", "theta",
"itera")
else:
res_vars = ("tau", "sensible", "latent", "monob", "cd", "cd10n",
"ct", "ct10n", "cq", "cq10n", "tsrv", "tsr", "qsr",
"usr", "psim", "psit", "psiq", "psim_ref", "psit_ref",
"psiq_ref", "u10n", "t10n", "q10n", "zo", "zot", "zoq",
"uref", "tref", "qref", "dter", "dqer", "dtwl", "tkt",
"qair", "qsea", "Rl", "Rs", "Rnl", "ug", "usrGF",
"GustFact", "Rb", "rh", "rho", "cp", "lv", "theta",
"itera")
elif out_var == "limited":
res_vars = ("tau", "sensible", "latent", "uref", "tref", "qref")
else:
res_vars = out_var
return res_vars
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