Repository: ESMValGroup/ESMValCore

Issue:
ICON: specifiy vertical grid file in recipe
ICON 3D output does not always have the vertical grid information (i.e., the altitude) available. In this case, this information is stored in a separate file. It would be helpful to specify this file (absolute or relative to `auxiliary_data_dir` in the recipe).

Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2067?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup) Report
> Merging [#2067](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2067?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup) (f4b5248) into [main](https://app.codecov.io/gh/ESMValGroup/ESMValCore/commit/83b12a233b9bd3ae3dc37540db9ac07d14e3b613?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup) (83b12a2) will **increase** coverage by `0.02%`.
> The diff coverage is `100.00%`.

```diff
@@            Coverage Diff             @@
##             main    #2067      +/-   ##
==========================================
+ Coverage   92.79%   92.82%   +0.02%     
==========================================
  Files         236      236              
  Lines       12485    12534      +49     
==========================================
+ Hits        11586    11635      +49     
  Misses        899      899              
```


| [Impacted Files](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2067?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup) | Coverage Δ | |
|---|---|---|
| [esmvalcore/cmor/\_fixes/icon/\_base\_fixes.py](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2067?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup#diff-ZXNtdmFsY29yZS9jbW9yL19maXhlcy9pY29uL19iYXNlX2ZpeGVzLnB5) | `100.00% <100.00%> (ø)` | |
| [esmvalcore/cmor/\_fixes/icon/icon.py](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2067?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup#diff-ZXNtdmFsY29yZS9jbW9yL19maXhlcy9pY29uL2ljb24ucHk=) | `100.00% <100.00%> (ø)` | |
| [esmvalcore/cmor/\_fixes/native\_datasets.py](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2067?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup#diff-ZXNtdmFsY29yZS9jbW9yL19maXhlcy9uYXRpdmVfZGF0YXNldHMucHk=) | `100.00% <100.00%> (ø)` | |



Patch:
diff --git a/doc/quickstart/find_data.rst b/doc/quickstart/find_data.rst
index fcd27a66a..9f9963928 100644
--- a/doc/quickstart/find_data.rst
+++ b/doc/quickstart/find_data.rst
@@ -370,33 +370,15 @@ This is useful for diagnostics that do not support input data in UGRID format
 if you want to use the built-in :ref:`unstructured_nearest scheme <built-in
 regridding schemes>` regridding scheme.
 
-Similar to any other fix, the ICON fix allows the use of :ref:`extra
-facets<extra_facets>`.
-By default, the file :download:`icon-mappings.yml
-</../esmvalcore/config/extra_facets/icon-mappings.yml>` is used for that
-purpose.
-For some variables, extra facets are necessary; otherwise ESMValTool cannot
-read them properly.
-Supported keys for extra facets are:
-
-============= ============================= =================================
-Key           Description                   Default value if not specified
-============= ============================= =================================
-``latitude``  Standard name of the latitude ``latitude``
-              coordinate in the raw input
-              file
-``longitude`` Standard name of the          ``longitude``
-              longitude coordinate in the
-              raw input file
-``raw_name``  Variable name of the          CMOR variable name of the
-              variable in the raw input     corresponding variable
-              file
-``ugrid``     Automatic UGRIDization of     ``True``
-              the input data
-``var_type``  Variable type of the          No default (needs to be specified
-              variable in the raw input     in extra facets or recipe if
-              file                          default DRS is used)
-============= ============================= =================================
+For 3D ICON variables, ESMValCore tries to add the pressure level information
+(from the variables `pfull` and `phalf`) and/or altitude information (from the
+variables `zg` and `zghalf`) to the preprocessed output files.
+If neither of these variables are available in the input files, it is possible
+to specify the location of files that include the corresponding `zg` or
+`zghalf` variables with the facets ``zg_file`` and/or ``zghalf_file`` in the
+recipe or the extra facets.
+The paths to these files can be specified absolute or relative (to
+``auxiliary_data_dir`` as defined in the :ref:`user configuration file`).
 
 .. hint::
 
@@ -404,6 +386,9 @@ Key           Description                   Default value if not specified
    native ICON data, you need to specify the name of the vertical coordinate
    (e.g., ``coordinate: air_pressure``) since native ICON output usually
    provides a 3D air pressure field instead of a simple 1D vertical coordinate.
+   This also works if your files only contain altitude information (in this
+   case, the US standard atmosphere is used to convert between altitude and
+   pressure levels; see :ref:`Vertical interpolation` for details).
    Example:
 
    .. code-block:: yaml
@@ -415,6 +400,46 @@ Key           Description                   Default value if not specified
           scheme: linear
           coordinate: air_pressure
 
+Similar to any other fix, the ICON fix allows the use of :ref:`extra
+facets<extra_facets>`.
+By default, the file :download:`icon-mappings.yml
+</../esmvalcore/config/extra_facets/icon-mappings.yml>` is used for that
+purpose.
+For some variables, extra facets are necessary; otherwise ESMValTool cannot
+read them properly.
+Supported keys for extra facets are:
+
+=================== ================================ ===================================
+Key                 Description                      Default value if not specified
+=================== ================================ ===================================
+``latitude``        Standard name of the latitude    ``latitude``
+                    coordinate in the raw input
+                    file
+``longitude``       Standard name of the             ``longitude``
+                    longitude coordinate in the
+                    raw input file
+``raw_name``        Variable name of the             CMOR variable name of the
+                    variable in the raw input        corresponding variable
+                    file
+``ugrid``           Automatic UGRIDization of        ``True``
+                    the input data
+``var_type``        Variable type of the             No default (needs to be specified
+                    variable in the raw input        in extra facets or recipe if
+                    file                             default DRS is used)
+``zg_file``         Absolute or relative (to         If possible, use `zg` variable
+                    ``auxiliary_data_dir`` defined   provided by the raw input file
+                    in the
+                    :ref:`user configuration file`)
+                    path to the input file that
+                    contains `zg`
+``zghalf_file``     Absolute or relative (to         If possible, use `zghalf` variable
+                    ``auxiliary_data_dir`` defined   provided by the raw input file
+                    in the
+                    :ref:`user configuration file`)
+                    path to the input file that
+                    contains `zghalf`
+=================== ================================ ===================================
+
 .. hint::
 
    In order to read cell area files (``areacella`` and ``areacello``), one
diff --git a/environment.yml b/environment.yml
index 26809991e..66c4782c8 100644
--- a/environment.yml
+++ b/environment.yml
@@ -35,7 +35,7 @@ dependencies:
   - py-cordex
   - pybtex
   - python>=3.9
-  - python-stratify
+  - python-stratify>=0.3
   - pyyaml
   - requests
   - scipy>=1.6
diff --git a/esmvalcore/cmor/_fixes/icon/_base_fixes.py b/esmvalcore/cmor/_fixes/icon/_base_fixes.py
index f2502ac0a..0590c2026 100644
--- a/esmvalcore/cmor/_fixes/icon/_base_fixes.py
+++ b/esmvalcore/cmor/_fixes/icon/_base_fixes.py
@@ -1,6 +1,7 @@
 """Fix base classes for ICON on-the-fly CMORizer."""
 
 import logging
+import os
 import warnings
 from datetime import datetime
 from pathlib import Path
@@ -178,6 +179,75 @@ class IconFix(NativeDatasetFix):
 
         return (node_lat, node_lon)
 
+    def _get_path_from_facet(self, facet, description=None):
+        """Try to get path from facet."""
+        if description is None:
+            description = 'File'
+        path = Path(os.path.expandvars(self.extra_facets[facet])).expanduser()
+        if not path.is_file():
+            new_path = self.session['auxiliary_data_dir'] / path
+            if not new_path.is_file():
+                raise FileNotFoundError(
+                    f"{description} '{path}' given by facet '{facet}' does "
+                    f"not exist (specify a valid absolute path or a path "
+                    f"relative to the auxiliary_data_dir "
+                    f"'{self.session['auxiliary_data_dir']}')"
+                )
+            path = new_path
+        return path
+
+    def add_additional_cubes(self, cubes):
+        """Add additional user-defined cubes to list of cubes (in-place).
+
+        An example use case is adding a vertical coordinate (e.g., `zg`) to the
+        dataset if the vertical coordinate data is stored in a separate ICON
+        output file.
+
+        Currently, the following cubes can be added:
+
+        - `zg` (`geometric_height_at_full_level_center`) from facet `zg_file`.
+          This can be used as vertical coordinate.
+        - `zghalf` (`geometric_height_at_half_level_center`) from facet
+          `zghalf_file`. This can be used as bounds for the vertical
+          coordinate.
+
+        Note
+        ----
+        Files can be specified as absolute or relative (to
+        ``auxiliary_data_dir`` as defined in the :ref:`user configuration
+        file`) paths.
+
+        Parameters
+        ----------
+        cubes: iris.cube.CubeList
+            Input cubes which will be modified in place.
+
+        Returns
+        -------
+        iris.cube.CubeList
+            Modified cubes. The cubes are modified in place; they are just
+            returned out of convenience for easy access.
+
+        Raises
+        ------
+        InputFilesNotFound
+            A specified file does not exist.
+
+        """
+        facets_to_consider = [
+            'zg_file',
+            'zghalf_file',
+        ]
+        for facet in facets_to_consider:
+            if facet not in self.extra_facets:
+                continue
+            path_to_add = self._get_path_from_facet(facet)
+            logger.debug("Adding cubes from %s", path_to_add)
+            new_cubes = self._load_cubes(path_to_add)
+            cubes.extend(new_cubes)
+
+        return cubes
+
     def get_horizontal_grid(self, cube):
         """Get copy of ICON horizontal grid from global attribute of cube.
 
@@ -318,6 +388,14 @@ class IconFix(NativeDatasetFix):
                 category=UserWarning,
                 module='iris',
             )
+            warnings.filterwarnings(
+                'ignore',
+                message="Failed to create 'height' dimension coordinate: The "
+                        "'height' DimCoord bounds array must be strictly "
+                        "monotonic.",
+                category=UserWarning,
+                module='iris',
+            )
             cubes = iris.load(str(path))
         return cubes
 
diff --git a/esmvalcore/cmor/_fixes/icon/icon.py b/esmvalcore/cmor/_fixes/icon/icon.py
index bef17d5a9..aef9c9c14 100644
--- a/esmvalcore/cmor/_fixes/icon/icon.py
+++ b/esmvalcore/cmor/_fixes/icon/icon.py
@@ -24,6 +24,7 @@ class AllVars(IconFix):
 
     def fix_metadata(self, cubes):
         """Fix metadata."""
+        cubes = self.add_additional_cubes(cubes)
         cube = self.get_cube(cubes)
 
         # Fix time
@@ -152,46 +153,66 @@ class AllVars(IconFix):
             f"'{self.vardef.short_name}', cube and other cubes in file do not "
             f"contain it")
 
+    def _get_z_coord(self, cubes, points_name, bounds_name=None):
+        """Get z-coordinate without metadata (reversed)."""
+        points_cube = iris.util.reverse(
+            cubes.extract_cube(NameConstraint(var_name=points_name)),
+            'height',
+        )
+        points = points_cube.core_data()
+
+        # Get bounds if possible
+        if bounds_name is not None:
+            bounds_cube = iris.util.reverse(
+                cubes.extract_cube(NameConstraint(var_name=bounds_name)),
+                'height',
+            )
+            bounds = bounds_cube.core_data()
+            bounds = da.stack(
+                (bounds[..., :-1, :], bounds[..., 1:, :]), axis=-1
+            )
+        else:
+            bounds = None
+
+        z_coord = AuxCoord(
+            points,
+            bounds=bounds,
+            units=points_cube.units,
+        )
+        return z_coord
+
     def _fix_height(self, cube, cubes):
         """Fix height coordinate of cube."""
         # Reverse entire cube along height axis so that index 0 is surface
         # level
         cube = iris.util.reverse(cube, 'height')
 
-        # Add air_pressure coordinate if possible
-        # (make sure to also reverse pressure cubes)
+        # If possible, extract reversed air_pressure coordinate from list of
+        # cubes and add it to cube
+        # Note: pfull/phalf have dimensions (time, height, spatial_dim)
         if cubes.extract(NameConstraint(var_name='pfull')):
-            plev_points_cube = iris.util.reverse(
-                cubes.extract_cube(NameConstraint(var_name='pfull')),
-                'height',
-            )
-            air_pressure_points = plev_points_cube.core_data()
-
-            # Get bounds from half levels and reshape array
             if cubes.extract(NameConstraint(var_name='phalf')):
-                plev_bounds_cube = iris.util.reverse(
-                    cubes.extract_cube(NameConstraint(var_name='phalf')),
-                    'height',
-                )
-                air_pressure_bounds = plev_bounds_cube.core_data()
-                air_pressure_bounds = da.stack(
-                    (air_pressure_bounds[:, :-1], air_pressure_bounds[:, 1:]),
-                    axis=-1)
+                phalf = 'phalf'
             else:
-                air_pressure_bounds = None
-
-            # Setup air pressure coordinate with correct metadata and add to
-            # cube
-            air_pressure_coord = AuxCoord(
-                air_pressure_points,
-                bounds=air_pressure_bounds,
-                var_name='plev',
-                standard_name='air_pressure',
-                long_name='pressure',
-                units=plev_points_cube.units,
-                attributes={'positive': 'down'},
-            )
-            cube.add_aux_coord(air_pressure_coord, np.arange(cube.ndim))
+                phalf = None
+            plev_coord = self._get_z_coord(cubes, 'pfull', bounds_name=phalf)
+            self.fix_plev_metadata(cube, plev_coord)
+            cube.add_aux_coord(plev_coord, np.arange(cube.ndim))
+
+        # If possible, extract reversed altitude coordinate from list of cubes
+        # and add it to cube
+        # Note: zg/zghalf have dimensions (height, spatial_dim)
+        if cubes.extract(NameConstraint(var_name='zg')):
+            if cubes.extract(NameConstraint(var_name='zghalf')):
+                zghalf = 'zghalf'
+            else:
+                zghalf = None
+            alt_coord = self._get_z_coord(cubes, 'zg', bounds_name=zghalf)
+            self.fix_alt16_metadata(cube, alt_coord)
+
+            # Altitude coordinate only spans height and spatial dimensions (no
+            # time) -> these are always the last two dimensions in the cube
+            cube.add_aux_coord(alt_coord, np.arange(cube.ndim)[-2:])
 
         # Fix metadata
         z_coord = cube.coord('height')
diff --git a/esmvalcore/cmor/_fixes/native_datasets.py b/esmvalcore/cmor/_fixes/native_datasets.py
index 72e9561c6..1e14e5c95 100644
--- a/esmvalcore/cmor/_fixes/native_datasets.py
+++ b/esmvalcore/cmor/_fixes/native_datasets.py
@@ -231,6 +231,36 @@ class NativeDatasetFix(Fix):
         coord.long_name = 'time'
         return coord
 
+    @staticmethod
+    def fix_alt16_metadata(cube, coord=None):
+        """Fix metadata of alt16 coordinate (in-place).
+
+        Parameters
+        ----------
+        cube: iris.cube.Cube
+            Input cube.
+        coord: str or iris.coords.Coord or None, optional (default: None)
+            Coordinate for which metadata will be fixed in-place. If ``None``,
+            assume the coordinate's name is `altitude`.
+
+        Returns
+        -------
+        iris.coords.AuxCoord or iris.coords.DimCoord
+            Fixed altitude coordinate. The coordinate is altered in-place; it
+            is just returned out of convenience for easy access.
+
+        """
+        if coord is None:
+            coord = cube.coord('altitude')
+        elif isinstance(coord, str):
+            coord = cube.coord(coord)
+        coord.var_name = 'alt16'
+        coord.standard_name = 'altitude'
+        coord.long_name = 'altitude'
+        coord.convert_units('m')
+        coord.attributes['positive'] = 'up'
+        return coord
+
     @staticmethod
     def fix_height_metadata(cube, coord=None):
         """Fix metadata of height coordinate (in-place).
diff --git a/esmvalcore/preprocessor/_other.py b/esmvalcore/preprocessor/_other.py
index 36bbccafb..14594ba7c 100644
--- a/esmvalcore/preprocessor/_other.py
+++ b/esmvalcore/preprocessor/_other.py
@@ -4,6 +4,7 @@ import logging
 from collections import defaultdict
 
 import dask.array as da
+import numpy as np
 
 logger = logging.getLogger(__name__)
 
@@ -71,3 +72,16 @@ def _group_products(products, by_key):
 
     grouped = _groupby(products, keyfunc=grouper)
     return grouped.items()
+
+
+def get_array_module(*args):
+    """Return the best matching array module.
+
+    If at least one of the arguments is a :class:`dask.array.Array` object,
+    the :mod:`dask.array` module is returned. In all other cases the
+    :mod:`numpy` module is returned.
+    """
+    for arg in args:
+        if isinstance(arg, da.Array):
+            return da
+    return np
diff --git a/esmvalcore/preprocessor/_regrid.py b/esmvalcore/preprocessor/_regrid.py
index a14abfbe6..4972a7ce8 100644
--- a/esmvalcore/preprocessor/_regrid.py
+++ b/esmvalcore/preprocessor/_regrid.py
@@ -14,13 +14,13 @@ from typing import Dict
 import iris
 import numpy as np
 import stratify
-from dask import array as da
 from geopy.geocoders import Nominatim
 from iris.analysis import AreaWeighted, Linear, Nearest, UnstructuredNearest
 from iris.util import broadcast_to_shape
 
 from ..cmor._fixes.shared import add_altitude_from_plev, add_plev_from_altitude
 from ..cmor.table import CMOR_TABLES
+from ._other import get_array_module
 from ._regrid_esmpy import ESMF_REGRID_METHODS
 from ._regrid_esmpy import regrid as esmpy_regrid
 from ._supplementary_vars import add_ancillary_variable, add_cell_measure
@@ -797,22 +797,19 @@ def _vertical_interpolate(cube, src_levels, levels, interpolation,
                                               cube.coord_dims(src_levels))
 
     # force mask onto data as nan's
-    cube.data = da.ma.filled(cube.core_data(), np.nan)
+    npx = get_array_module(cube.core_data())
+    data = npx.ma.filled(cube.core_data(), np.nan)
 
     # Now perform the actual vertical interpolation.
     new_data = stratify.interpolate(levels,
                                     src_levels_broadcast,
-                                    cube.core_data(),
+                                    data,
                                     axis=z_axis,
                                     interpolation=interpolation,
                                     extrapolation=extrapolation)
 
     # Calculate the mask based on the any NaN values in the interpolated data.
-    mask = np.isnan(new_data)
-
-    if np.any(mask):
-        # Ensure that the data is masked appropriately.
-        new_data = np.ma.array(new_data, mask=mask, fill_value=_MDI)
+    new_data = npx.ma.masked_where(npx.isnan(new_data), new_data)
 
     # Construct the resulting cube with the interpolated data.
     return _create_cube(cube, new_data, src_levels, levels.astype(float))
diff --git a/setup.py b/setup.py
index 05d5a9ac5..5743b3f60 100755
--- a/setup.py
+++ b/setup.py
@@ -58,7 +58,7 @@ REQUIREMENTS = {
         'scipy>=1.6',
         'scitools-iris>=3.4.0',
         'shapely[vectorized]',
-        'stratify',
+        'stratify>=0.3',
         'yamale',
     ],
     # Test dependencies


Test patch:
diff --git a/tests/integration/cmor/_fixes/icon/test_icon.py b/tests/integration/cmor/_fixes/icon/test_icon.py
index ada031f39..686cf80a5 100644
--- a/tests/integration/cmor/_fixes/icon/test_icon.py
+++ b/tests/integration/cmor/_fixes/icon/test_icon.py
@@ -1,5 +1,6 @@
 """Tests for the ICON on-the-fly CMORizer."""
 from datetime import datetime
+from pathlib import Path
 from unittest import mock
 
 import iris
@@ -15,6 +16,7 @@ from esmvalcore.cmor._fixes.icon._base_fixes import IconFix
 from esmvalcore.cmor._fixes.icon.icon import AllVars, Clwvi, Siconc, Siconca
 from esmvalcore.cmor.fix import Fix
 from esmvalcore.cmor.table import CoordinateInfo, get_var_info
+from esmvalcore.config import CFG
 from esmvalcore.config._config import get_extra_facets
 from esmvalcore.dataset import Dataset
 
@@ -90,7 +92,25 @@ def cubes_2d_lat_lon_grid():
     return CubeList([cube])
 
 
-def _get_fix(mip, short_name, fix_name):
+@pytest.fixture
+def simple_unstructured_cube():
+    """Simple cube with unstructured grid."""
+    time_coord = DimCoord([0], var_name='time', standard_name='time',
+                          units='days since 1850-01-01')
+    height_coord = DimCoord([0, 1, 2], var_name='height')
+    lat_coord = AuxCoord([0.0, 1.0], var_name='lat', standard_name='latitude',
+                         long_name='latitude', units='degrees_north')
+    lon_coord = AuxCoord([0.0, 1.0], var_name='lon',
+                         standard_name='longitude', long_name='longitude',
+                         units='degrees_east')
+    cube = Cube([[[0.0, 1.0], [2.0, 3.0], [4.0, 5.0]]], var_name='ta',
+                units='K',
+                dim_coords_and_dims=[(time_coord, 0), (height_coord, 1)],
+                aux_coords_and_dims=[(lat_coord, 2), (lon_coord, 2)])
+    return cube
+
+
+def _get_fix(mip, short_name, fix_name, session=None):
     """Load a fix from esmvalcore.cmor._fixes.icon.icon."""
     dataset = Dataset(
         project='ICON',
@@ -101,26 +121,26 @@ def _get_fix(mip, short_name, fix_name):
     extra_facets = get_extra_facets(dataset, ())
     vardef = get_var_info(project='ICON', mip=mip, short_name=short_name)
     cls = getattr(esmvalcore.cmor._fixes.icon.icon, fix_name)
-    fix = cls(vardef, extra_facets=extra_facets)
+    fix = cls(vardef, extra_facets=extra_facets, session=session)
     return fix
 
 
-def get_fix(mip, short_name):
+def get_fix(mip, short_name, session=None):
     """Load a variable fix from esmvalcore.cmor._fixes.icon.icon."""
     fix_name = short_name[0].upper() + short_name[1:]
-    return _get_fix(mip, short_name, fix_name)
+    return _get_fix(mip, short_name, fix_name, session=session)
 
 
-def get_allvars_fix(mip, short_name):
+def get_allvars_fix(mip, short_name, session=None):
     """Load the AllVars fix from esmvalcore.cmor._fixes.icon.icon."""
-    return _get_fix(mip, short_name, 'AllVars')
+    return _get_fix(mip, short_name, 'AllVars', session=session)
 
 
-def fix_metadata(cubes, mip, short_name):
+def fix_metadata(cubes, mip, short_name, session=None):
     """Fix metadata of cubes."""
-    fix = get_fix(mip, short_name)
+    fix = get_fix(mip, short_name, session=session)
     cubes = fix.fix_metadata(cubes)
-    fix = get_allvars_fix(mip, short_name)
+    fix = get_allvars_fix(mip, short_name, session=session)
     cubes = fix.fix_metadata(cubes)
     return cubes
 
@@ -175,18 +195,20 @@ def check_time(cube):
     assert time.attributes == {}
 
 
-def check_height(cube, plev_has_bounds=True):
-    """Check height coordinate of cube."""
+def check_model_level_metadata(cube):
+    """Check metadata of model_level coordinate."""
     assert cube.coords('model level number', dim_coords=True)
     height = cube.coord('model level number', dim_coords=True)
     assert height.var_name == 'model_level'
     assert height.standard_name is None
     assert height.long_name == 'model level number'
     assert height.units == 'no unit'
-    np.testing.assert_array_equal(height.points, np.arange(47))
-    assert height.bounds is None
     assert height.attributes == {'positive': 'up'}
+    return height
 
+
+def check_air_pressure_metadata(cube):
+    """Check metadata of air_pressure coordinate."""
     assert cube.coords('air_pressure', dim_coords=False)
     plev = cube.coord('air_pressure', dim_coords=False)
     assert plev.var_name == 'plev'
@@ -194,6 +216,16 @@ def check_height(cube, plev_has_bounds=True):
     assert plev.long_name == 'pressure'
     assert plev.units == 'Pa'
     assert plev.attributes == {'positive': 'down'}
+    return plev
+
+
+def check_height(cube, plev_has_bounds=True):
+    """Check height coordinate of cube."""
+    height = check_model_level_metadata(cube)
+    np.testing.assert_array_equal(height.points, np.arange(47))
+    assert height.bounds is None
+
+    plev = check_air_pressure_metadata(cube)
     assert cube.coord_dims('air_pressure') == (0, 1, 2)
 
     np.testing.assert_allclose(
@@ -1553,3 +1585,168 @@ def test_get_mesh_not_cached(monkeypatch):
     monkeypatch.setattr(fix, '_create_mesh', mock.Mock())
     fix.get_mesh(cube)
     fix._create_mesh.assert_called_once_with(cube)
+
+
+# Test _get_path_from_facet
+
+
+@pytest.mark.parametrize(
+    'path,description,output',
+    [
+        ('{tmp_path}/a.nc', None, '{tmp_path}/a.nc'),
+        ('b.nc', 'Grid file', '{tmp_path}/b.nc'),
+    ],
+)
+def test_get_path_from_facet(path, description, output, tmp_path):
+    """Test fix."""
+    session = CFG.start_session('my session')
+    session['auxiliary_data_dir'] = tmp_path
+    path = path.format(tmp_path=tmp_path)
+    fix = get_allvars_fix('Amon', 'tas', session=session)
+    fix.extra_facets['test_path'] = path
+
+    # Create empty dummy file
+    output = output.format(tmp_path=tmp_path)
+    with open(output, 'w'):
+        pass
+
+    out_path = fix._get_path_from_facet('test_path', description=description)
+
+    assert isinstance(out_path, Path)
+    assert out_path == Path(output.format(tmp_path=tmp_path))
+
+
+@pytest.mark.parametrize(
+    'path,description',
+    [
+        ('{tmp_path}/a.nc', None),
+        ('b.nc', 'Grid file'),
+    ],
+)
+def test_get_path_from_facet_fail(path, description, tmp_path):
+    """Test fix."""
+    session = CFG.start_session('my session')
+    session['auxiliary_data_dir'] = tmp_path
+    path = path.format(tmp_path=tmp_path)
+    fix = get_allvars_fix('Amon', 'tas', session=session)
+    fix.extra_facets['test_path'] = path
+
+    with pytest.raises(FileNotFoundError, match=description):
+        fix._get_path_from_facet('test_path', description=description)
+
+
+# Test add_additional_cubes
+
+
+@pytest.mark.parametrize('facet', ['zg_file', 'zghalf_file'])
+@pytest.mark.parametrize('path', ['{tmp_path}/a.nc', 'a.nc'])
+def test_add_additional_cubes(path, facet, tmp_path):
+    """Test fix."""
+    session = CFG.start_session('my session')
+    session['auxiliary_data_dir'] = tmp_path
+    path = path.format(tmp_path=tmp_path)
+    fix = get_allvars_fix('Amon', 'tas', session=session)
+    fix.extra_facets[facet] = path
+
+    # Save temporary cube
+    cube = Cube(0, var_name=facet)
+    iris.save(cube, tmp_path / 'a.nc')
+
+    cubes = CubeList([])
+    new_cubes = fix.add_additional_cubes(cubes)
+
+    assert new_cubes is cubes
+    assert len(cubes) == 1
+    assert cubes[0].var_name == facet
+
+
+@pytest.mark.parametrize('facet', ['zg_file', 'zghalf_file'])
+@pytest.mark.parametrize('path', ['{tmp_path}/a.nc', 'a.nc'])
+def test_add_additional_cubes_fail(path, facet, tmp_path):
+    """Test fix."""
+    session = CFG.start_session('my session')
+    session['auxiliary_data_dir'] = tmp_path
+    path = path.format(tmp_path=tmp_path)
+    fix = get_allvars_fix('Amon', 'tas', session=session)
+    fix.extra_facets[facet] = path
+
+    cubes = CubeList([])
+    with pytest.raises(FileNotFoundError, match='File'):
+        fix.add_additional_cubes(cubes)
+
+
+# Test _fix_height
+
+
+@pytest.mark.parametrize('bounds', [True, False])
+def test_fix_height_plev(bounds, simple_unstructured_cube):
+    """Test fix."""
+    cube = simple_unstructured_cube[:, 1:, :]
+    pfull_cube = simple_unstructured_cube[:, 1:, :]
+    pfull_cube.var_name = 'pfull'
+    pfull_cube.units = 'Pa'
+    cubes = CubeList([cube, pfull_cube])
+    if bounds:
+        phalf_cube = simple_unstructured_cube.copy()
+        phalf_cube.var_name = 'phalf'
+        phalf_cube.units = 'Pa'
+        cubes.append(phalf_cube)
+    fix = get_allvars_fix('Amon', 'ta')
+
+    fixed_cube = fix._fix_height(cube, cubes)
+
+    expected_data = [[[4.0, 5.0], [2.0, 3.0]]]
+    np.testing.assert_allclose(fixed_cube.data, expected_data)
+
+    height = check_model_level_metadata(fixed_cube)
+    np.testing.assert_array_equal(height.points, [0, 1])
+    assert height.bounds is None
+
+    plev = check_air_pressure_metadata(fixed_cube)
+    assert fixed_cube.coord_dims('air_pressure') == (0, 1, 2)
+    np.testing.assert_allclose(plev.points, expected_data)
+    if bounds:
+        expected_bnds = [[[[4.0, 2.0], [5.0, 3.0]], [[2.0, 0.0], [3.0, 1.0]]]]
+        np.testing.assert_allclose(plev.bounds, expected_bnds)
+    else:
+        assert plev.bounds is None
+
+
+@pytest.mark.parametrize('bounds', [True, False])
+def test_fix_height_alt16(bounds, simple_unstructured_cube):
+    """Test fix."""
+    cube = simple_unstructured_cube[:, 1:, :]
+    zg_cube = simple_unstructured_cube[0, 1:, :]
+    zg_cube.var_name = 'zg'
+    zg_cube.units = 'm'
+    cubes = CubeList([cube, zg_cube])
+    if bounds:
+        zghalf_cube = simple_unstructured_cube[0, :, :]
+        zghalf_cube.var_name = 'zghalf'
+        zghalf_cube.units = 'm'
+        cubes.append(zghalf_cube)
+    fix = get_allvars_fix('Amon', 'ta')
+
+    fixed_cube = fix._fix_height(cube, cubes)
+
+    expected_data = [[[4.0, 5.0], [2.0, 3.0]]]
+    np.testing.assert_allclose(fixed_cube.data, expected_data)
+
+    height = check_model_level_metadata(fixed_cube)
+    np.testing.assert_array_equal(height.points, [0, 1])
+    assert height.bounds is None
+
+    assert fixed_cube.coords('altitude', dim_coords=False)
+    alt16 = fixed_cube.coord('altitude', dim_coords=False)
+    assert alt16.var_name == 'alt16'
+    assert alt16.standard_name == 'altitude'
+    assert alt16.long_name == 'altitude'
+    assert alt16.units == 'm'
+    assert alt16.attributes == {'positive': 'up'}
+    assert fixed_cube.coord_dims('altitude') == (1, 2)
+    np.testing.assert_allclose(alt16.points, expected_data[0])
+    if bounds:
+        expected_bnds = [[[4.0, 2.0], [5.0, 3.0]], [[2.0, 0.0], [3.0, 1.0]]]
+        np.testing.assert_allclose(alt16.bounds, expected_bnds)
+    else:
+        assert alt16.bounds is None
diff --git a/tests/integration/cmor/_fixes/test_native_datasets.py b/tests/integration/cmor/_fixes/test_native_datasets.py
index a48b37e95..371c2c5a2 100644
--- a/tests/integration/cmor/_fixes/test_native_datasets.py
+++ b/tests/integration/cmor/_fixes/test_native_datasets.py
@@ -41,6 +41,11 @@ def sample_cube():
         long_name='air_pressure',
         units='hPa',
     )
+    alt16_coord = AuxCoord(
+        [2.0, 4.0],
+        long_name='altitude',
+        units='km',
+    )
     height_coord = AuxCoord(
         [2.0, 4.0],
         long_name='height',
@@ -68,7 +73,8 @@ def sample_cube():
                              (plev_coord, 1),
                              (lat_coord, 2),
                              (lon_coord, 3)],
-        aux_coords_and_dims=[(height_coord, 1),
+        aux_coords_and_dims=[(alt16_coord, 1),
+                             (height_coord, 1),
                              (coord_with_bounds, 1)],
     )
     return cube
@@ -375,6 +381,46 @@ def test_fix_time_metadata_from_coord(sample_cube, fix):
     assert out_coord.bounds is None
 
 
+def test_fix_alt16_metadata(sample_cube, fix):
+    """Test ``fix_alt16_metadata``."""
+    out_coord = fix.fix_alt16_metadata(sample_cube)
+    assert out_coord is sample_cube.coord('altitude')
+    assert out_coord.standard_name == 'altitude'
+    assert out_coord.var_name == 'alt16'
+    assert out_coord.long_name == 'altitude'
+    assert out_coord.units == 'm'
+    assert out_coord.attributes['positive'] == 'up'
+    np.testing.assert_allclose(out_coord.points, [2000.0, 4000.0])
+    assert out_coord.bounds is None
+
+
+def test_fix_alt16_metadata_from_str(sample_cube, fix):
+    """Test ``fix_alt16_metadata`` from string."""
+    out_coord = fix.fix_alt16_metadata(sample_cube, coord='altitude')
+    assert out_coord is sample_cube.coord('altitude')
+    assert out_coord.standard_name == 'altitude'
+    assert out_coord.var_name == 'alt16'
+    assert out_coord.long_name == 'altitude'
+    assert out_coord.units == 'm'
+    assert out_coord.attributes['positive'] == 'up'
+    np.testing.assert_allclose(out_coord.points, [2000.0, 4000.0])
+    assert out_coord.bounds is None
+
+
+def test_fix_alt16_metadata_from_coord(sample_cube, fix):
+    """Test ``fix_alt16_metadata`` from string."""
+    coord = AuxCoord([2.0], units='m')
+    out_coord = fix.fix_alt16_metadata(sample_cube, coord=coord)
+    assert out_coord is coord
+    assert out_coord.standard_name == 'altitude'
+    assert out_coord.var_name == 'alt16'
+    assert out_coord.long_name == 'altitude'
+    assert out_coord.units == 'm'
+    assert out_coord.attributes['positive'] == 'up'
+    np.testing.assert_allclose(out_coord.points, [2.0])
+    assert out_coord.bounds is None
+
+
 def test_fix_height_metadata(sample_cube, fix):
     """Test ``fix_height_metadata``."""
     out_coord = fix.fix_height_metadata(sample_cube)
diff --git a/tests/integration/preprocessor/_regrid/test_extract_levels.py b/tests/integration/preprocessor/_regrid/test_extract_levels.py
index 8687f1922..b0869244a 100644
--- a/tests/integration/preprocessor/_regrid/test_extract_levels.py
+++ b/tests/integration/preprocessor/_regrid/test_extract_levels.py
@@ -15,6 +15,7 @@ from tests.unit.preprocessor._regrid import _make_cube, _make_vcoord
 
 
 class Test(tests.Test):
+
     def setUp(self):
         """Prepare tests."""
         shape = (3, 2, 2)
@@ -57,20 +58,52 @@ class Test(tests.Test):
         levels = [0.5, 1.5]
         scheme = 'linear'
         result = extract_levels(self.cube, levels, scheme)
-        expected = np.array([[[[2., 3.], [4., 5.]], [[6., 7.], [8., 9.]]],
-                             [[[14., 15.], [16., 17.]], [[18., 19.],
-                                                         [20., 21.]]]])
+        expected = np.ma.array([
+            [
+                [[2., 3.], [4., 5.]],
+                [[6., 7.], [8., 9.]],
+            ],
+            [
+                [[14., 15.], [16., 17.]],
+                [[18., 19.], [20., 21.]],
+            ],
+        ])
         self.assert_array_equal(result.data, expected)
         self.shape[self.z_dim] = len(levels)
         self.assertEqual(result.shape, tuple(self.shape))
 
+    def test_interpolation__linear_lazy(self):
+        levels = [0.5, 1.5]
+        scheme = 'linear'
+        cube = self.cube.copy(self.cube.lazy_data())
+        result = extract_levels(cube, levels, scheme)
+        self.assertTrue(result.has_lazy_data())
+        expected = np.ma.array([
+            [
+                [[2., 3.], [4., 5.]],
+                [[6., 7.], [8., 9.]],
+            ],
+            [
+                [[14., 15.], [16., 17.]],
+                [[18., 19.], [20., 21.]],
+            ],
+        ])
+        self.assert_array_equal(result.data, expected)
+
     def test_interpolation__nearest(self):
         levels = [0.49, 1.51]
         scheme = 'nearest'
         result = extract_levels(self.cube, levels, scheme)
-        expected = np.array([[[[0., 1.], [2., 3.]], [[8., 9.], [10., 11.]]],
-                             [[[12., 13.], [14., 15.]], [[20., 21.],
-                                                         [22., 23.]]]])
+        expected = np.ma.array([
+            [
+                [[0., 1.], [2., 3.]],
+                [[8., 9.], [10., 11.]],
+            ],
+            [
+                [[12., 13.], [14., 15.]],
+                [[20., 21.], [22., 23.]],
+            ],
+        ])
         self.assert_array_equal(result.data, expected)
         self.shape[self.z_dim] = len(levels)
         self.assertEqual(result.shape, tuple(self.shape))
diff --git a/tests/unit/preprocessor/_other/test_other.py b/tests/unit/preprocessor/_other/test_other.py
index ab98895c0..803b66d99 100644
--- a/tests/unit/preprocessor/_other/test_other.py
+++ b/tests/unit/preprocessor/_other/test_other.py
@@ -2,6 +2,7 @@
 
 import unittest
 
+import dask.array as da
 import iris.coord_categorisation
 import iris.coords
 import numpy as np
@@ -10,11 +11,16 @@ from iris.cube import Cube
 from numpy.testing import assert_array_equal
 
 from esmvalcore.preprocessor import PreprocessorFile
-from esmvalcore.preprocessor._other import _group_products, clip
+from esmvalcore.preprocessor._other import (
+    _group_products,
+    clip,
+    get_array_module,
+)
 
 
 class TestOther(unittest.TestCase):
     """Test class for _other."""
+
     def test_clip(self):
         """Test clip function."""
         cube = Cube(np.array([-10, 0, 10]))
@@ -67,5 +73,23 @@ def test_group_products_string_list():
     assert grouped_by_list == grouped_by_string
 
 
+def test_get_array_module_da():
+
+    npx = get_array_module(da.array([1, 2]))
+    assert npx is da
+
+
+def test_get_array_module_np():
+
+    npx = get_array_module(np.array([1, 2]))
+    assert npx is np
+
+
+def test_get_array_module_mixed():
+
+    npx = get_array_module(da.array([1]), np.array([1]))
+    assert npx is da
+
+
 if __name__ == '__main__':
     unittest.main()
diff --git a/tests/unit/preprocessor/_regrid/test_extract_levels.py b/tests/unit/preprocessor/_regrid/test_extract_levels.py
index 638d61fbc..f53fbac70 100644
--- a/tests/unit/preprocessor/_regrid/test_extract_levels.py
+++ b/tests/unit/preprocessor/_regrid/test_extract_levels.py
@@ -1,6 +1,4 @@
 """Unit tests for :func:`esmvalcore.preprocessor.regrid.extract_levels`."""
-import os
-import tempfile
 import unittest
 from unittest import mock
 
@@ -35,8 +33,6 @@ class Test(tests.Test):
         self.schemes = [
             'linear', 'nearest', 'linear_extrapolate', 'nearest_extrapolate',
         ]
-        descriptor, self.filename = tempfile.mkstemp('.nc')
-        os.close(descriptor)
 
     def test_invalid_scheme__unknown(self):
         levels = mock.sentinel.levels
@@ -160,7 +156,7 @@ class Test(tests.Test):
         # Check the _create_cube args ...
         self.assertEqual(len(args), 4)
         self.assertEqual(args[0], self.cube)
-        self.assert_array_equal(args[1], new_data)
+        self.assert_array_equal(args[1], np.ma.array(new_data))
         self.assert_array_equal(args[2],
                                 self.cube.coord(axis='z', dim_coords=True))
         self.assert_array_equal(args[3], levels)
@@ -281,15 +277,8 @@ class Test(tests.Test):
         masked = ma.empty(self.shape)
         masked.mask = mask
         cube = _make_cube(masked, dtype=self.dtype)
-        # save cube to test the lazy data interpolation too
-        iris.save(cube, self.filename)
         with mock.patch('stratify.interpolate',
                         return_value=new_data) as mocker:
-            # first test lazy
-            loaded_cube = iris.load_cube(self.filename)
-            result_from_lazy = extract_levels(loaded_cube, levels, scheme)
-            self.assertEqual(result_from_lazy, self.created_cube)
-            # then test realized
             result = extract_levels(cube, levels, scheme)
             self.assertEqual(result, self.created_cube)
             args, kwargs = mocker.call_args
@@ -300,7 +289,7 @@ class Test(tests.Test):
             src_levels_broadcast = np.broadcast_to(pts.reshape(self.z, 1, 1),
                                                    cube.shape)
             self.assert_array_equal(args[1], src_levels_broadcast)
-            self.assert_array_equal(args[2], cube.data)
+            self.assert_array_equal(args[2], np.ma.filled(masked, np.nan))
             # Check the stratify.interpolate kwargs ...
             self.assertEqual(
                 kwargs, dict(axis=0, interpolation=scheme,