Introduction
Obukhov Length is a well-known length scale in boundary layer meteorology. It characterizes stability effects (buoyancy versus shear turbulence) in the turbulent layer of the atmosphere near the surface. It is used in air pollution modelling and to characterize the shape of wind profiles in wind energy applications (among many other applications).
The following procedure describes how the Obukhov length (L) can be calculated from ERA5 data.
The formula used in the code, computes the inverse of the Obukhov length (1/L [m-1]) as eq (5.7c) from Stull (1988):
where vk = 0.4 is the Von Karman constant, g is the gravitational acceleration, tvst is the turbulent temperature scale, tv2 is the virtual temperature and ust is the friction velocity.
The output is a 2D field from input data taken from the Climate Data Store (CDS) ERA5 hourly data on single levels from 1959 to present dataset.
You will need:
- Python3
The CDS API installed. For more details, please follow the instructions here.
- The ecCodes library to read and write GRIB data.
Step 1: Download input data from the CDS
The variables that are needed from the ERA5 data, are:
- 'surface_pressure',
- '2m_dewpoint_temperature',
- '2m_temperature',
- 'instantaneous_eastward_turbulent_surface_stress',
- 'instantaneous_moisture_flux',
- 'instantaneous_northward_turbulent_surface_stress',
- 'instantaneous_surface_sensible_heat_flux',
- 'standard_deviation_of_filtered_subgrid_orography'
Please be careful about the ECMWF ERA5 sign convention for fluxes. All downward fluxes are positive, so during daytime the surface sensible heat flux (sshf) is typically negative.
The Python script used the CDS API to download the ERA5 data from the CDS. The suggested procedure is:
- Copy the script below to a text editor on your computer
- Edit the 'year', 'month', 'day', 'time', and output filename ('data_2020010112.grib' in this example) in the script to meet your requirements. Note: 'grid' and 'area' keywords can also be included.
- Save the script (for example with the filename as 'get_data.py')
Run the script:
python3 get_data.py
# **************************** LICENSE START ***********************************
#
# Copyright 2022 ECMWF. This software is distributed under the terms
# of the Apache License version 2.0. In applying this license, ECMWF does not
# waive the privileges and immunities granted to it by virtue of its status as
# an Intergovernmental Organization or submit itself to any jurisdiction.
#
# ***************************** LICENSE END ************************************
import cdsapi
c = cdsapi.Client()
c.retrieve('reanalysis-era5-single-levels', {
'product_type': 'reanalysis',
'format': 'grib',
'variable': [
'surface_pressure',
'2m_dewpoint_temperature',
'2m_temperature',
'instantaneous_eastward_turbulent_surface_stress',
'instantaneous_moisture_flux',
'instantaneous_northward_turbulent_surface_stress',
'instantaneous_surface_sensible_heat_flux',
'standard_deviation_of_filtered_subgrid_orography'
],
'year': '2020',
'month': '01',
'day': '01',
'time': '12:00',
}, 'data_2020010112.grib')
Running the script produces a file in the current working directory called 'data_2020010112.grib' . This is a GRIB file containing the ERA5 variables needed.
Step 2: Compute inverse Obukhov length 1/L
The suggested procedure to run the Python script to compute the inverse Obukhov length (1/L) is :
- Copy the script below to a text editor
- Save the script as 'calculate-Obukhov-length.py'
- Run the script 'calculate-Obukhov-length.py' with the correct arguments:
python3 calculate-Obukhov-length.py data_2020010112.grib output.grib
# **************************** LICENSE START ***********************************
#
# Copyright 2022 ECMWF. This software is distributed under the terms
# of the Apache License version 2.0. In applying this license, ECMWF does not
# waive the privileges and immunities granted to it by virtue of its status as
# an Intergovernmental Organization or submit itself to any jurisdiction.
#
# ***************************** LICENSE END ************************************
import sys
import numpy as np
from eccodes import codes_count_in_file, codes_get, codes_get_values
from eccodes import codes_clone, codes_set, codes_set_values
from eccodes import codes_set_key_vals, codes_write
from eccodes import codes_grib_new_from_file, codes_release
from eccodes import CodesInternalError
def qswat(t, p):
"""
Computes saturation q (with respect to water)
In t : Temperature (K)
p : Pressure (Pa)
Out qswt: Saturation specific humidity (kg/kg)
"""
rkbol = 1.380658e-23
rnavo = 6.0221367e+23
r = rnavo * rkbol
rmd = 28.9644
rmv = 18.0153
rd = 1000 * r / rmd # Gas constant for air
rv = 1000 * r / rmv # Gas constant for water vapour
restt = 611.21
r2es = restt * rd / rv
r3les = 17.502
r4les = 32.19
retv = rv / rd - 1 # Tv = T(1+retv*q)
rtt = 273.16 # Melting point (0 Celcius)
foeew = r2es * np.exp((r3les * (t - rtt)) / (t - r4les))
qs = foeew / p
zcor = 1 / (1 - retv * qs)
qs = qs * zcor
return qs
def calculate_Obukhov_length(data_file, output_file):
try:
gid_list = None
with open(data_file) as f_in:
mcount = codes_count_in_file(f_in)
gid_list = [codes_grib_new_from_file(f_in) for i in range(mcount)]
# Input data: sp, 2d, 2t, ie, ishf, iews, inss
data = {}
for i in range(mcount):
gid = gid_list[i]
sn = codes_get(gid, 'shortName')
vals = codes_get_values(gid) # Type: class 'numpy.ndarray'
data[sn] = vals
rd = 287.06 # Gas constant
retv = 0.6078 # Tv = T*(1+retv*q)
cp = 1004.7 # Air heat capacity at constant pressure
lmelt = 2500800 # Latent heat of evaporation
vk = 0.4 # VonKarman constant
g = 9.81 # Gravity acceleration
q2 = qswat(data['2d'], data['sp']) # q2 is saturation value at 2d
tv2 = data['2t'] * (1 + retv * q2) # Virtual temperature
rho = data['sp'] / (rd * tv2) # Air density
tau = np.sqrt(data['iews']**2 + data['inss']**2)
# Turb. surface stress
ust = np.maximum(np.sqrt(tau / rho), 0.001) # Friction velocity
wt = -data['ishf'] / (rho * cp) # Turbulent heat flux
wq = -data['ie'] / rho # Turbulent moisture flux
wtv = wt + retv * data['2t'] * wq # Virtual turbulent heat flux
tvst = -wtv / ust # Turbulent temperature scale
qst = -wq / ust # Turbulent moisture scale
Linv = vk * g * tvst / (tv2 * ust**2) # Inverse Obukhov length
# Stull eq. (5.7c)
Linvrange = 100. # Range of useful Linv values to control
# grib-coding
Linv = np.maximum(np.minimum(Linv, Linvrange), -Linvrange)
print('Linv max/min/mean: %e %e %e' %
(Linv.max(), Linv.min(), np.abs(Linv).mean()))
# Write to a file
with open(output_file, 'wb') as f_out:
gid = codes_clone(gid_list[0])
codes_set(gid, 'paramId', 80) # Experimental parameter
codes_set_values(gid, Linv)
codes_write(gid, f_out)
codes_release(gid)
except CodesInternalError as ex:
sys.stderr.write('Error: %s\n' % ex)
sys.exit(-100)
finally:
if gid_list != None:
for gid in gid_list:
codes_release(gid)
def print_usage():
print('python calculate-Obukhov-length.py data.grib output.grib')
def main(argv):
if len(argv) != 2:
print_usage()
sys.exit(100)
calculate_Obukhov_length(argv[0], argv[1])
if __name__ == '__main__':
main(sys.argv[1:])
Running the script produces a file in the current working directory called 'output.grib' . This is a GRIB file containing the inverse of the Obukhov length.
If you require the output data in netCDF, the GRIB file with the results can be converted to netCDF on your local system by using the ecCodes command grib_to_netcdf.
Limitations
In areas with significant orography, the inverse Obukhov length (1/L), as computed above, may be less usefull for two reasons: (a) Monin-Obkhov similarity becomes questionable in areas with orography, and (b) the ERA5 momentum fluxes, as needed for the computation of 1/L, contain the sum of turbulent flux and momentum flux due to unresolved small scale orography. The latter is parametrized by the so-called Turbulent Orographic Form Drag scheme (TOFD). To identify such areas, the retrieval job also includes parameter sdfor (standard deviation of filtered subgrid orography). It is recommended to use 1/L only in areas where sdfor < 50 m. In that case the contribution of TOFD to the total turbulent stress is small (Beljaars et al. 2004).
References
Stull, Roland B., 1988: An Introduction to Boundary Layer Meteorology. Springer, 670 pp.
Beljaars, A., Brown, A. and Wood, N., 2004: A new parametrization of turbulent orographic form drag. QJRMS, 130, pp.1327-1347.
