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.