Repository: ESMValGroup/ESMValCore

Issue:
Numba causing installation issues with python=3.9
conda=4.9.2 env builds well but installing doesn't:
```
    Running setup.py install for llvmlite ... error
    ERROR: Command errored out with exit status 1:
     command: /home/valeriu/miniconda3/envs/esmvalcore-492/bin/python3.9 -u -c 'import sys, setuptools, tokenize; sys.argv[0] = '"'"'/tmp/pip-install-elou0qm8/llvmlite/setup.py'"'"'; __file__='"'"'/tmp/pip-install-elou0qm8/llvmlite/setup.py'"'"';f=getattr(tokenize, '"'"'open'"'"', open)(__file__);code=f.read().replace('"'"'\r\n'"'"', '"'"'\n'"'"');f.close();exec(compile(code, __file__, '"'"'exec'"'"'))' install --record /tmp/pip-record-woct46sr/install-record.txt --single-version-externally-managed --compile --install-headers /home/valeriu/miniconda3/envs/esmvalcore-492/include/python3.9/llvmlite
         cwd: /tmp/pip-install-elou0qm8/llvmlite/
    Complete output (29 lines):
    running install
    running build
    got version from file /tmp/pip-install-elou0qm8/llvmlite/llvmlite/_version.py {'version': '0.34.0', 'full': 'c5889c9e98c6b19d5d85ebdd982d64a03931f8e2'}
    running build_ext
    /home/valeriu/miniconda3/envs/esmvalcore-492/bin/python3.9 /tmp/pip-install-elou0qm8/llvmlite/ffi/build.py
    LLVM version... Traceback (most recent call last):
      File "/tmp/pip-install-elou0qm8/llvmlite/ffi/build.py", line 105, in main_posix
        out = subprocess.check_output([llvm_config, '--version'])
      File "/home/valeriu/miniconda3/envs/esmvalcore-492/lib/python3.9/subprocess.py", line 420, in check_output
        return run(*popenargs, stdout=PIPE, timeout=timeout, check=True,
      File "/home/valeriu/miniconda3/envs/esmvalcore-492/lib/python3.9/subprocess.py", line 501, in run
        with Popen(*popenargs, **kwargs) as process:
      File "/home/valeriu/miniconda3/envs/esmvalcore-492/lib/python3.9/subprocess.py", line 947, in __init__
        self._execute_child(args, executable, preexec_fn, close_fds,
      File "/home/valeriu/miniconda3/envs/esmvalcore-492/lib/python3.9/subprocess.py", line 1819, in _execute_child
        raise child_exception_type(errno_num, err_msg, err_filename)
    FileNotFoundError: [Errno 2] No such file or directory: 'llvm-config'
    
    During handling of the above exception, another exception occurred:
    
    Traceback (most recent call last):
      File "/tmp/pip-install-elou0qm8/llvmlite/ffi/build.py", line 191, in <module>
        main()
      File "/tmp/pip-install-elou0qm8/llvmlite/ffi/build.py", line 181, in main
        main_posix('linux', '.so')
      File "/tmp/pip-install-elou0qm8/llvmlite/ffi/build.py", line 107, in main_posix
        raise RuntimeError("%s failed executing, please point LLVM_CONFIG "
    RuntimeError: llvm-config failed executing, please point LLVM_CONFIG to the path for llvm-config
    error: command '/home/valeriu/miniconda3/envs/esmvalcore-492/bin/python3.9' failed with exit code 1
    ----------------------------------------
ERROR: Command errored out with exit status 1: /home/valeriu/miniconda3/envs/esmvalcore-492/bin/python3.9 -u -c 'import sys, setuptools, tokenize; sys.argv[0] = '"'"'/tmp/pip-install-elou0qm8/llvmlite/setup.py'"'"'; __file__='"'"'/tmp/pip-install-elou0qm8/llvmlite/setup.py'"'"';f=getattr(tokenize, '"'"'open'"'"', open)(__file__);code=f.read().replace('"'"'\r\n'"'"', '"'"'\n'"'"');f.close();exec(compile(code, __file__, '"'"'exec'"'"'))' install --record /tmp/pip-record-woct46sr/install-record.txt --single-version-externally-managed --compile --install-headers /home/valeriu/miniconda3/envs/esmvalcore-492/include/python3.9/llvmlite Check the logs for full command output.
```

Discussion:
valeriupredoi: having a look right now! :beer: 

Patch:
diff --git a/doc/requirements.txt b/doc/requirements.txt
index a055635a3..7d67ccad6 100644
--- a/doc/requirements.txt
+++ b/doc/requirements.txt
@@ -2,7 +2,6 @@ autodocsumm
 dask[array]
 fiona
 netCDF4
-numba
 numpy
 pandas
 pillow
diff --git a/esmvalcore/preprocessor/_derive/_shared.py b/esmvalcore/preprocessor/_derive/_shared.py
index 260242ec9..17556fc9b 100644
--- a/esmvalcore/preprocessor/_derive/_shared.py
+++ b/esmvalcore/preprocessor/_derive/_shared.py
@@ -1,8 +1,10 @@
 """Auxiliary derivation functions used for multiple variables."""
 
 import logging
+from copy import deepcopy
 
 import iris
+import numpy as np
 
 from esmvalcore.iris_helpers import var_name_constraint
 
@@ -25,3 +27,104 @@ def cloud_area_fraction(cubes, tau_constraint, plev_constraint):
         new_cube = new_cube.collapsed('air_pressure', iris.analysis.SUM)
 
     return new_cube
+
+
+def pressure_level_widths(cube, ps_cube, top_limit=0.0):
+    """Create a cube with pressure level widths.
+
+    This is done by taking a 2D surface pressure field as lower bound.
+
+    Parameters
+    ----------
+    cube : iris.cube.Cube
+        Cube containing ``air_pressure`` coordiante.
+    ps_cube : iris.cube.Cube
+        Cube containing ``surface_air_pressure``.
+    top_limit : float
+        Pressure in Pa.
+
+    Returns
+    -------
+    iris.cube.Cube
+        Cube of same shape as ``cube`` containing pressure level widths.
+
+    """
+    pressure_array = _create_pressure_array(cube, ps_cube, top_limit)
+
+    data = _get_pressure_level_widths(pressure_array)
+    p_level_widths_cube = cube.copy(data=data)
+    p_level_widths_cube.rename('pressure level widths')
+    p_level_widths_cube.units = ps_cube.units
+
+    return p_level_widths_cube
+
+
+def _create_pressure_array(cube, ps_cube, top_limit):
+    """Create an array filled with the ``air_pressure`` coord values.
+
+    The array is created from ``cube`` with the same dimensions as ``cube``.
+    This array is then sandwiched with a 2D array containing the surface
+    pressure and a 2D array containing the top pressure limit.
+
+    """
+    # Create 4D array filled with pressure level values
+    p_levels = cube.coord('air_pressure').points.astype(np.float32)
+    p_4d_array = iris.util.broadcast_to_shape(p_levels, cube.shape, [1])
+
+    # Create 4d array filled with surface pressure values
+    shape = cube.shape
+    ps_4d_array = iris.util.broadcast_to_shape(ps_cube.data, shape, [0, 2, 3])
+
+    # Set pressure levels below the surface pressure to NaN
+    pressure_4d = np.where((ps_4d_array - p_4d_array) < 0, np.NaN, p_4d_array)
+
+    # Make top_limit last pressure level
+    top_limit_array = np.full(ps_cube.shape, top_limit, dtype=np.float32)
+    data = top_limit_array[:, np.newaxis, :, :]
+    pressure_4d = np.concatenate((pressure_4d, data), axis=1)
+
+    # Make surface pressure the first pressure level
+    data = ps_cube.data[:, np.newaxis, :, :]
+    pressure_4d = np.concatenate((data, pressure_4d), axis=1)
+
+    return pressure_4d
+
+
+def _get_pressure_level_widths(array, air_pressure_axis=1):
+    """Compute pressure level widths.
+
+    For a 1D array with pressure level columns, return a 1D array with
+    pressure level widths.
+
+    """
+    array = np.copy(array)
+    if np.any(np.diff(array, axis=air_pressure_axis) > 0.0):
+        raise ValueError("Pressure level value increased with height")
+
+    # Calculate centers
+    indices = [slice(None)] * array.ndim
+    array_shifted = np.roll(array, -1, axis=air_pressure_axis)
+    index_0 = deepcopy(indices)
+    index_0[air_pressure_axis] = 0
+    array_shifted[tuple(index_0)] = array[tuple(index_0)]
+    index_neg1 = deepcopy(indices)
+    index_neg1[air_pressure_axis] = -1
+    index_neg2 = deepcopy(indices)
+    index_neg2[air_pressure_axis] = -2
+    array[tuple(index_neg2)] = array[tuple(index_neg1)]
+    array_centers = (array + array_shifted) / 2.0
+    index_range_neg1 = deepcopy(indices)
+    index_range_neg1[air_pressure_axis] = slice(None, -1)
+    array_centers = array_centers[tuple(index_range_neg1)]
+
+    # Remove NaNs (replace them with surface pressure)
+    dim_map = np.arange(array_centers.ndim)
+    dim_map = np.delete(dim_map, air_pressure_axis)
+    array_centers_surface = iris.util.broadcast_to_shape(
+        array_centers[tuple(index_0)], array_centers.shape, dim_map)
+    array_centers = np.where(np.isnan(array_centers), array_centers_surface,
+                             array_centers)
+
+    # Calculate level widths
+    p_level_widths = -np.diff(array_centers, axis=air_pressure_axis)
+    return p_level_widths
diff --git a/esmvalcore/preprocessor/_derive/toz.py b/esmvalcore/preprocessor/_derive/toz.py
index e478250d5..b6de06a7c 100644
--- a/esmvalcore/preprocessor/_derive/toz.py
+++ b/esmvalcore/preprocessor/_derive/toz.py
@@ -2,17 +2,16 @@
 
 import cf_units
 import iris
-import numba
-import numpy as np
 from scipy import constants
 
 from ._baseclass import DerivedVariableBase
+from ._shared import pressure_level_widths
 
 # Constants
 AVOGADRO_CONST = constants.value('Avogadro constant')
 AVOGADRO_CONST_UNIT = constants.unit('Avogadro constant')
-STANDARD_GRAVITY = 9.81
-STANDARD_GRAVITY_UNIT = cf_units.Unit('m s^-2')
+STANDARD_GRAVITY = constants.value('standard acceleration of gravity')
+STANDARD_GRAVITY_UNIT = constants.unit('standard acceleration of gravity')
 MW_AIR = 29
 MW_AIR_UNIT = cf_units.Unit('g mol^-1')
 MW_O3 = 48
@@ -47,9 +46,9 @@ class DerivedVariable(DerivedVariableBase):
         ps_cube = cubes.extract_strict(
             iris.Constraint(name='surface_air_pressure'))
 
-        p_layer_widths = _pressure_level_widths(tro3_cube,
-                                                ps_cube,
-                                                top_limit=0.0)
+        p_layer_widths = pressure_level_widths(tro3_cube,
+                                               ps_cube,
+                                               top_limit=0.0)
         toz_cube = (tro3_cube * p_layer_widths / STANDARD_GRAVITY * MW_O3 /
                     MW_AIR)
         toz_cube = toz_cube.collapsed('air_pressure', iris.analysis.SUM)
@@ -63,130 +62,3 @@ class DerivedVariable(DerivedVariableBase):
         toz_cube.units = 'DU'
 
         return toz_cube
-
-
-# Helper functions
-def _pressure_level_widths(tro3_cube, ps_cube, top_limit=0.0):
-    """Create a cube with pressure level widths.
-
-    This is done by taking a 2D surface pressure field as lower bound.
-
-    Parameters
-    ----------
-        tro3_cube : iris.cube.Cube
-            `Cube` containing `mole_fraction_of_ozone_in_air`.
-        ps_cube : iris.cube.Cube
-            `Cube` containing `surface_air_pressure`.
-        top_limit : float
-            Pressure in Pa.
-
-    Returns
-    -------
-    iris.cube.Cube
-    `Cube` of same shape as `tro3_cube` containing pressure level widths.
-
-    """
-    pressure_array = _create_pressure_array(tro3_cube, ps_cube, top_limit)
-
-    data = _apply_pressure_level_widths(pressure_array)
-    p_level_widths_cube = tro3_cube.copy(data=data)
-    p_level_widths_cube.rename('pressure level widths')
-    p_level_widths_cube.units = ps_cube.units
-
-    return p_level_widths_cube
-
-
-def _create_pressure_array(tro3_cube, ps_cube, top_limit):
-    """Create an array filled with the `air_pressure` coord values.
-
-    The array is created from the `tro3_cube` with the same dimensions
-    as `tro3_cube`. This array is then sandwiched with a 2D array
-    containing the surface pressure and a 2D array containing the top
-    pressure limit.
-    """
-    # Create 4D array filled with pressure level values
-    p_levels = tro3_cube.coord('air_pressure').points.astype(np.float32)
-    p_4d_array = iris.util.broadcast_to_shape(p_levels, tro3_cube.shape, [1])
-
-    # Create 4d array filled with surface pressure values
-    shape = tro3_cube.shape
-    ps_4d_array = iris.util.broadcast_to_shape(ps_cube.data, shape, [0, 2, 3])
-
-    # Set pressure levels below the surface pressure to NaN
-    pressure_4d = np.where((ps_4d_array - p_4d_array) < 0, np.NaN, p_4d_array)
-
-    # Make top_limit last pressure level
-    top_limit_array = np.full(ps_cube.shape, top_limit, dtype=np.float32)
-    data = top_limit_array[:, np.newaxis, :, :]
-    pressure_4d = np.concatenate((pressure_4d, data), axis=1)
-
-    # Make surface pressure the first pressure level
-    data = ps_cube.data[:, np.newaxis, :, :]
-    pressure_4d = np.concatenate((data, pressure_4d), axis=1)
-
-    return pressure_4d
-
-
-def _apply_pressure_level_widths(array, air_pressure_axis=1):
-    """Compute pressure level widths.
-
-    For a 1D array with pressure level columns, return a 1D array with
-    pressure level widths.
-    """
-    return np.apply_along_axis(_p_level_widths, air_pressure_axis, array)
-
-
-@numba.jit()  # ~10x faster
-def _p_level_widths(array):
-    """Create pressure level widths.
-
-    The array with pressure levels is assumed to be monotonic and the
-    values are decreasing.
-
-    The first element is the lower boundary (surface pressure), the last
-    value is the upper boundary. Thicknesses are only calculated for the
-    values between these boundaries, the returned array, therefore,
-    contains two elements less.
-
-    >>> _p_level_widths(np.array([1020, 1000, 700, 500, 5]))
-    array([170., 250., 595.], dtype=float32)
-
-    >>> _p_level_widths(np.array([990, np.NaN, 700, 500, 5]))
-    array([  0., 390., 595.], dtype=float32)
-    """
-    surface_pressure = array[0]
-    top_limit = array[-1]
-    array = array[1:-1]
-
-    p_level_widths = np.full(array.shape, np.NAN, dtype=np.float32)
-
-    last_pressure_level = len(array) - 1
-    for i, val in enumerate(array):
-        # numba would otherwise initialize it to 0 and
-        # hide bugs that would occur in raw Python
-        bounds_width = np.NAN
-        if np.isnan(val):
-            bounds_width = 0
-        else:
-            # Distance to lower bound
-            if i == 0 or np.isnan(array[i - 1]):
-                # First pressure level with value
-                dist_to_lower_bound = surface_pressure - val
-            else:
-                dist_to_lower_bound = 0.5 * (array[i - 1] - val)
-
-            # Distance to upper bound
-            if i == last_pressure_level:  # last pressure level
-                dist_to_upper_bound = val - top_limit
-            else:
-                dist_to_upper_bound = 0.5 * (val - array[i + 1])
-
-            # Check monotonicity - all distances must be >= 0
-            if dist_to_lower_bound < 0.0 or dist_to_upper_bound < 0.0:
-                raise ValueError("Pressure level value increased with "
-                                 "height.")
-
-            bounds_width = dist_to_lower_bound + dist_to_upper_bound
-
-        p_level_widths[i] = bounds_width
-    return p_level_widths
diff --git a/package/meta.yaml b/package/meta.yaml
index e0963b043..f5ee7471e 100644
--- a/package/meta.yaml
+++ b/package/meta.yaml
@@ -48,7 +48,6 @@ requirements:
     - fire
     - nc-time-axis
     - netCDF4
-    - numba
     - numpy
     - prov
     - psutil
diff --git a/setup.py b/setup.py
index b9da19c27..d777e0739 100755
--- a/setup.py
+++ b/setup.py
@@ -36,7 +36,6 @@ REQUIREMENTS = {
         'fire',
         'nc-time-axis',  # needed by iris.plot
         'netCDF4',
-        'numba',
         'numpy',
         'prov[dot]',
         'psutil',


Test patch:
diff --git a/tests/unit/preprocessor/_derive/test_p_level_widths.py b/tests/unit/preprocessor/_derive/test_p_level_widths.py
deleted file mode 100644
index 8bba96fa2..000000000
--- a/tests/unit/preprocessor/_derive/test_p_level_widths.py
+++ /dev/null
@@ -1,88 +0,0 @@
-"""Tests for toz variable derivation functions."""
-
-import numpy as np
-import pytest
-
-from esmvalcore.preprocessor._derive.toz import _p_level_widths
-
-
-def test_col_is_not_monotonic():
-    """Test for non-monotonic column."""
-    plev = 1000
-    top_limit = 5
-    col = np.array([1, 2, 3, 2, 1])
-    col = np.insert(col, 0, plev)
-    col = np.append(col, top_limit)
-    with pytest.raises(ValueError):
-        _p_level_widths(col)
-
-
-def test_keeping_column_length():
-    """Test for level widths keeping column length."""
-    plev = 1000
-    top_limit = 5
-    col = np.array([1000, 900, 800])
-    col = np.insert(col, 0, plev)
-    col = np.append(col, top_limit)
-    assert len(_p_level_widths(col)) == len(col) - 2
-
-
-def test_low_lev_surf_press():
-    """Test for lowest level equal to surface pressure."""
-    plev = 1000
-    top_limit = 5
-    col = np.array([1000, 900, 800])
-    col = np.insert(col, 0, plev)
-    col = np.append(col, top_limit)
-    result = np.array([50, 100, 845])
-    assert all(_p_level_widths(col) == result)
-
-
-def test_low_lev_above_surf_press():
-    """Test for lowest level above surface pressure."""
-    plev = 1020
-    top_limit = 5
-    col = np.array([1000, 900, 800])
-    col = np.insert(col, 0, plev)
-    col = np.append(col, top_limit)
-    result = np.array([70, 100, 845])
-    assert all(_p_level_widths(col) == result)
-
-
-def test_low_lev_below_surf_press():
-    """Test for lowest level below surface pressure."""
-    plev = 970
-    top_limit = 5
-    col = np.array([np.NaN, 900, 800])
-    col = np.insert(col, 0, plev)
-    col = np.append(col, top_limit)
-    result = np.array([0, 120, 845])
-    assert all(_p_level_widths(col) == result)
-
-    col = np.array([np.NaN, np.NaN, 900, 800])
-    col = np.insert(col, 0, plev)
-    col = np.append(col, top_limit)
-    result = np.array([0, 0, 120, 845])
-    assert all(_p_level_widths(col) == result)
-
-
-def test_high_level_top_limit():
-    """Test for highest level equal to top limit."""
-    plev = 1020
-    top_limit = 5
-    col = np.array([1000, 900, 5])
-    col = np.insert(col, 0, plev)
-    col = np.append(col, top_limit)
-    result = np.array([70, 50 + 895 / 2, 895 / 2])
-    assert all(_p_level_widths(col) == result)
-
-
-def test_high_level_above_top_limit():
-    """Test for highest level above top limit."""
-    plev = 1020
-    top_limit = 5
-    col = np.array([1000, 900, 3])
-    col = np.insert(col, 0, plev)
-    col = np.append(col, top_limit)
-    with pytest.raises(ValueError):
-        _p_level_widths(col)
diff --git a/tests/unit/preprocessor/_derive/test_shared.py b/tests/unit/preprocessor/_derive/test_shared.py
new file mode 100644
index 000000000..77e4b84f6
--- /dev/null
+++ b/tests/unit/preprocessor/_derive/test_shared.py
@@ -0,0 +1,138 @@
+"""Unit tests for shared variable derivation functions."""
+
+import numpy as np
+import pytest
+
+from esmvalcore.preprocessor._derive._shared import _get_pressure_level_widths
+
+
+def test_col_is_not_monotonic():
+    """Test for non-monotonic column."""
+    plev = 1000
+    top_limit = 5
+    col = np.array([1, 2, 3, 2, 1])
+    col = np.insert(col, 0, plev)
+    col = np.append(col, top_limit)
+    with pytest.raises(ValueError):
+        _get_pressure_level_widths(col, air_pressure_axis=0)
+    with pytest.raises(ValueError):
+        _get_pressure_level_widths(np.atleast_2d(col), air_pressure_axis=1)
+    with pytest.raises(ValueError):
+        _get_pressure_level_widths(np.atleast_3d(col), air_pressure_axis=1)
+
+
+def test_keeping_column_length():
+    """Test for level widths keeping column length."""
+    plev = 1000
+    top_limit = 5
+    col = np.array([1000, 900, 800])
+    col = np.insert(col, 0, plev)
+    col = np.append(col, top_limit)
+    assert (len(_get_pressure_level_widths(col, air_pressure_axis=0)) ==
+            len(col) - 2)
+    col = np.atleast_2d(col)
+    assert (_get_pressure_level_widths(col, air_pressure_axis=1).shape ==
+            (1, 3))
+    col = np.atleast_3d(col)
+    assert (_get_pressure_level_widths(col, air_pressure_axis=1).shape ==
+            (1, 3, 1))
+
+
+def test_low_lev_surf_press():
+    """Test for lowest level equal to surface pressure."""
+    plev = 1000
+    top_limit = 5
+    col = np.array([1000, 900, 800])
+    col = np.insert(col, 0, plev)
+    col = np.append(col, top_limit)
+    result = np.array([50, 100, 845])
+    assert np.array_equal(_get_pressure_level_widths(col, air_pressure_axis=0),
+                          result)
+    col = np.atleast_2d(col)
+    assert np.array_equal(_get_pressure_level_widths(col, air_pressure_axis=1),
+                          np.atleast_2d(result))
+    col = np.atleast_3d(col)
+    assert np.array_equal(_get_pressure_level_widths(col, air_pressure_axis=1),
+                          np.atleast_3d(result))
+
+
+def test_low_lev_above_surf_press():
+    """Test for lowest level above surface pressure."""
+    plev = 1020
+    top_limit = 5
+    col = np.array([1000, 900, 800])
+    col = np.insert(col, 0, plev)
+    col = np.append(col, top_limit)
+    result = np.array([70, 100, 845])
+    assert np.array_equal(_get_pressure_level_widths(col, air_pressure_axis=0),
+                          result)
+    col = np.atleast_2d(col)
+    assert np.array_equal(_get_pressure_level_widths(col, air_pressure_axis=1),
+                          np.atleast_2d(result))
+    col = np.atleast_3d(col)
+    assert np.array_equal(_get_pressure_level_widths(col, air_pressure_axis=1),
+                          np.atleast_3d(result))
+
+
+def test_low_lev_below_surf_press():
+    """Test for lowest level below surface pressure."""
+    plev = 970
+    top_limit = 5
+    col = np.array([np.NaN, 900, 800])
+    col = np.insert(col, 0, plev)
+    col = np.append(col, top_limit)
+    result = np.array([0, 120, 845])
+    assert np.array_equal(_get_pressure_level_widths(col, air_pressure_axis=0),
+                          result)
+    col = np.atleast_2d(col)
+    assert np.array_equal(_get_pressure_level_widths(col, air_pressure_axis=1),
+                          np.atleast_2d(result))
+    col = np.atleast_3d(col)
+    assert np.array_equal(_get_pressure_level_widths(col, air_pressure_axis=1),
+                          np.atleast_3d(result))
+
+    col = np.array([np.NaN, np.NaN, 900, 800])
+    col = np.insert(col, 0, plev)
+    col = np.append(col, top_limit)
+    result = np.array([0, 0, 120, 845])
+    assert np.array_equal(_get_pressure_level_widths(col, air_pressure_axis=0),
+                          result)
+    col = np.atleast_2d(col)
+    assert np.array_equal(_get_pressure_level_widths(col, air_pressure_axis=1),
+                          np.atleast_2d(result))
+    col = np.atleast_3d(col)
+    assert np.array_equal(_get_pressure_level_widths(col, air_pressure_axis=1),
+                          np.atleast_3d(result))
+
+
+def test_high_level_top_limit():
+    """Test for highest level equal to top limit."""
+    plev = 1020
+    top_limit = 5
+    col = np.array([1000, 900, 5])
+    col = np.insert(col, 0, plev)
+    col = np.append(col, top_limit)
+    result = np.array([70, 50 + 895 / 2, 895 / 2])
+    assert np.array_equal(_get_pressure_level_widths(col, air_pressure_axis=0),
+                          result)
+    col = np.atleast_2d(col)
+    assert np.array_equal(_get_pressure_level_widths(col, air_pressure_axis=1),
+                          np.atleast_2d(result))
+    col = np.atleast_3d(col)
+    assert np.array_equal(_get_pressure_level_widths(col, air_pressure_axis=1),
+                          np.atleast_3d(result))
+
+
+def test_high_level_above_top_limit():
+    """Test for highest level above top limit."""
+    plev = 1020
+    top_limit = 5
+    col = np.array([1000, 900, 3])
+    col = np.insert(col, 0, plev)
+    col = np.append(col, top_limit)
+    with pytest.raises(ValueError):
+        _get_pressure_level_widths(col, air_pressure_axis=0)
+    with pytest.raises(ValueError):
+        _get_pressure_level_widths(np.atleast_2d(col), air_pressure_axis=1)
+    with pytest.raises(ValueError):
+        _get_pressure_level_widths(np.atleast_3d(col), air_pressure_axis=1)
diff --git a/tests/unit/preprocessor/_derive/test_toz.py b/tests/unit/preprocessor/_derive/test_toz.py
new file mode 100644
index 000000000..262deab1e
--- /dev/null
+++ b/tests/unit/preprocessor/_derive/test_toz.py
@@ -0,0 +1,84 @@
+"""Test derivation of ``toz``."""
+import dask.array as da
+import iris
+import numpy as np
+import pytest
+
+import esmvalcore.preprocessor._derive.toz as toz
+from .test_co2s import get_coord_spec, get_ps_cube
+
+
+@pytest.fixture
+def masked_cubes():
+    """Masked O3 cube."""
+    coord_spec = get_coord_spec()
+    o3_data = da.ma.masked_less([[[[0.0, -1.0],
+                                   [-1.0, -1.0]],
+                                  [[1.0, 2.0],
+                                   [3.0, -1.0]],
+                                  [[2.0, 2.0],
+                                   [2.0, 2.0]]]], 0.0)
+    o3_cube = iris.cube.Cube(
+        o3_data,
+        var_name='o3',
+        standard_name='mole_fraction_of_ozone_in_air',
+        units='1e-9',
+        dim_coords_and_dims=coord_spec,
+    )
+    ps_cube = get_ps_cube()
+    return iris.cube.CubeList([o3_cube, ps_cube])
+
+
+@pytest.fixture
+def unmasked_cubes():
+    """Unmasked O3 cube."""
+    coord_spec = get_coord_spec()
+    o3_data = da.array([[[[2.0, 1.0],
+                          [0.8, 1.0]],
+                         [[1.5, 0.8],
+                          [2.0, 3.0]],
+                         [[4.0, 1.0],
+                          [3.0, 2.0]]]])
+    o3_cube = iris.cube.Cube(
+        o3_data,
+        var_name='o3',
+        standard_name='mole_fraction_of_ozone_in_air',
+        units='1e-9',
+        dim_coords_and_dims=coord_spec,
+    )
+    ps_cube = get_ps_cube()
+    return iris.cube.CubeList([o3_cube, ps_cube])
+
+
+def test_toz_calculate_masked_cubes(masked_cubes):
+    """Test function ``calculate`` with masked cube."""
+    derived_var = toz.DerivedVariable()
+    out_cube = derived_var.calculate(masked_cubes)
+    assert not np.ma.is_masked(out_cube.data)
+    np.testing.assert_allclose(out_cube.data,
+                               [[[1.2988646378902597, 0.7871906896304607],
+                                 [1.6924599827054907, 0.9446288275565529]]])
+    assert out_cube.units == 'DU'
+
+
+def test_toz_calculate_unmasked_cubes(unmasked_cubes):
+    """Test function ``calculate`` with unmasked cube."""
+    derived_var = toz.DerivedVariable()
+    out_cube = derived_var.calculate(unmasked_cubes)
+    assert not np.ma.is_masked(out_cube.data)
+    np.testing.assert_allclose(out_cube.data,
+                               [[[2.65676858, 0.39359534],
+                                 [2.04669579, 0.94462883]]])
+    assert out_cube.units == 'DU'
+
+
+@pytest.mark.parametrize('project,out', [
+    ('CMIP5', [{'short_name': 'tro3'}, {'short_name': 'ps'}]),
+    ('TEST', [{'short_name': 'tro3'}, {'short_name': 'ps'}]),
+    ('CMIP6', [{'short_name': 'o3'}, {'short_name': 'ps'}]),
+])
+def test_toz_required(project, out):
+    """Test function ``required``."""
+    derived_var = toz.DerivedVariable()
+    output = derived_var.required(project)
+    assert output == out