diff --git a/Code/AirSeaFluxCode.py b/Code/AirSeaFluxCode.py
index 74f80e6cbe918f75faebbd3d45bee11389974996..b0839adf6d13ae6947ee830d538a57c21eadf49d 100644
--- a/Code/AirSeaFluxCode.py
+++ b/Code/AirSeaFluxCode.py
@@ -10,7 +10,7 @@ from cs_wl_subs import *
class S88:
def _wind_iterate(self, ind):
- if self.gust[0] in [1, 2, 3]:
+ if self.gust[0] in range(1, 5):
self.wind[ind] = np.sqrt(np.power(np.copy(self.spd[ind]), 2) +
np.power(get_gust(self.gust[1],
self.theta[ind],
@@ -19,16 +19,14 @@ class S88:
self.gust[2],
self.grav[ind]), 2))
# ratio of gusty to horizontal wind
- self.GustFact[ind] = self.wind[ind]/self.spd[ind]
- # remove effect of gustiness following Fairall et al. (2003)
- self.u10n[ind] = self.spd[ind]-self.usr[ind]/kappa/np.sqrt(
- self.GustFact[ind])*(np.log(self.h_in[0, ind]/self.ref10) -
- self.psim[ind])
- if self.gust[0] == 2:
- self.GustFact[ind] = 1
+ self.GustFact = apply_GF(self.gust, self.spd, self.wind, "u")
+ 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])
+ if self.gust[0] == 3:
# option to not remove GustFact
self.u10n[ind] = self.wind[ind]-self.usr[ind]/kappa*(
- np.log(self.h_in[0, ind]/self.ref10)-self.psim[ind])
+ np.log(self.h_in[0, ind]/self.ref10)-self.psim[ind])
else:
# initalisation of wind
self.wind[ind] = np.copy(self.spd[ind])
@@ -330,15 +328,10 @@ class S88:
# 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.tau = self.rho*np.power(self.usr, 2)/self.GustFact
- if self.gust[0] == 3:
- self.sensible = self.rho*self.cp*self.usr*self.tsr
- self.latent = self.rho*self.lv*self.usr*self.qsr
- else:
- self.sensible = self.rho*self.cp*self.usr / \
- np.sqrt(self.GustFact)*self.tsr
- self.latent = self.rho*self.lv*self.usr / \
- np.sqrt(self.GustFact)*self.qsr
+ 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
# Set the new variables (for comparison against "old")
new = np.array([np.copy(getattr(self, i)) for i in new_vars])
@@ -395,9 +388,10 @@ 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, "u")
# remove effect of gustiness following Fairall et al. (2003)
# usr is divided by (GustFact)^0.5
- self.uref = self.spd-self.usr/kappa/np.sqrt(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))
# include lapse rate adjustment as theta is well-mixed
@@ -501,8 +495,7 @@ class S88:
gust = [0, 0, 0]
assert np.size(gust) == 3, "gust input must be a 3x1 array"
- assert gust[0] in [
- 0, 1, 2, 3], "gust at position 0 must be 0, 1, 2 or 3"
+ assert gust[0] in range(6), "gust at position 0 must be 0 to 5"
self.gust = gust
def _class_flag(self):
diff --git a/Code/flux_subs.py b/Code/flux_subs.py
index 171a9c8f7d3ffd4cfa489ffb7322612c9c572450..b24aa88361242f4246e51ec8391ff86c0ed54e7f 100755
--- a/Code/flux_subs.py
+++ b/Code/flux_subs.py
@@ -593,6 +593,64 @@ def get_gust(beta, Ta, usr, tsrv, zi, grav):
# ---------------------------------------------------------------------
+def apply_GF(gust, spd, wind, step):
+ """
+ Apply gustiness factor according if gustiness ON.
+
+ There are different ways to remove the effect of gustiness according to
+ the user's choice.
+
+ Parameters
+ ----------
+ gust : int
+ option on how to apply gustiness
+ 0: gustiness is switched OFF
+ 1: gustiness is switched ON following Fairall et al.
+ 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
+ 5: gustiness is switched ON following C35 matlab code
+ spd : float
+ wind speed [ms^{-1}]
+ wind : float
+ wind speed including gust [ms^{-1}]
+ step : str
+ step during AirSeaFluxCode the GF is applied: "u", "TSF"
+
+ Returns
+ -------
+ GustFact : float
+ gustiness factor.
+
+ """
+ # 1. following C35 documentation, 2. use GF to TSF, u10n uzout,
+ # 3. GF=1, 4. UA, 5. C35 code 6. ecmwf aerobulk)
+ # ratio of gusty to horizontal wind; gustiness factor
+ if step in ["u"]:
+ GustFact = wind*0+1
+ if gust[0] in [1, 2]:
+ GustFact = np.sqrt(wind/spd)
+ elif gust[0] == 5:
+ # as in C35 matlab code
+ GustFact = wind/spd
+ 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[0, :] = wind/spd
+ GustFact[1:3, :] = wind*0+1
+ # following Fairall et al. (2003)
+ if gust[0] == 2:
+ # usr is divided by (GustFact)^0.5 (here applied to sensible and
+ # latent as well as tau)
+ GustFact[1:3, :] = np.sqrt(wind/spd)
+ elif gust[0] == 3:
+ GustFact[0, :] = wind*0+1
+ return GustFact
+# ---------------------------------------------------------------------
+
+
def get_strs(hin, monob, wind, zo, zot, zoq, dt, dq, cd, ct, cq, meth):
"""
Calculate star wind speed, temperature and specific humidity.