diff --git a/doc/quickstart/find_data.rst b/doc/quickstart/find_data.rst
index f33a4741b1..3d2a09a55a 100644
--- a/doc/quickstart/find_data.rst
+++ b/doc/quickstart/find_data.rst
@@ -363,6 +363,32 @@ Key           Description                   Default value if not specified
               file                          default DRS is used)
 ============= ============================= =================================
 
+ESMValCore automatically makes native ICON data `UGRID
+<https://ugrid-conventions.github.io/ugrid-conventions/>`__-compliant when
+loading the data.
+The UGRID conventions provide a standardized format to store data on
+unstructured grids, which is required by many software packages or tools to
+work correctly.
+An example is the horizontal regridding of native ICON data to a regular grid.
+While the built-in :ref:`unstructured_nearest scheme <built-in regridding
+schemes>` can handle unstructured grids not in UGRID format, using more complex
+regridding algorithms (for example provided by the
+:doc:`iris-esmf-regrid:index` package through :ref:`generic regridding
+schemes`) requires the input data in UGRID format.
+The following code snippet provides a preprocessor that regrids native ICON
+data to a 1°x1° grid using `ESMF's first-order conservative regridding
+algorithm <https://earthsystemmodeling.org/regrid/#regridding-methods>`__:
+
+.. code-block:: yaml
+
+   preprocessors:
+     regrid_icon:
+       regrid:
+         target_grid: 1x1
+         scheme:
+           reference: esmf_regrid.experimental.unstructured_scheme:regrid_unstructured_to_rectilinear
+           method: conservative
+
 .. hint::
 
    To use the :func:`~esmvalcore.preprocessor.extract_levels` preprocessor on
diff --git a/doc/recipe/preprocessor.rst b/doc/recipe/preprocessor.rst
index dc947514cd..c64dba23fe 100644
--- a/doc/recipe/preprocessor.rst
+++ b/doc/recipe/preprocessor.rst
@@ -932,17 +932,17 @@ arguments is also supported. The call function for such schemes must be defined
 The `regrid` module will automatically pass the source and grid cubes as inputs
 of the scheme. An example of this usage is
 the :func:`~esmf_regrid.schemes.regrid_rectilinear_to_rectilinear`
-scheme available in :doc:`iris-esmf-regrid:index`:git 
+scheme available in :doc:`iris-esmf-regrid:index`:
 
 .. code-block:: yaml
 
-    preprocessors:
-      regrid_preprocessor:
-        regrid:
-          target_grid: 2.5x2.5
-          scheme:
-            reference: esmf_regrid.schemes:regrid_rectilinear_to_rectilinear
-            mdtol: 0.7
+  preprocessors:
+    regrid_preprocessor:
+      regrid:
+        target_grid: 2.5x2.5
+        scheme:
+          reference: esmf_regrid.schemes:regrid_rectilinear_to_rectilinear
+          mdtol: 0.7
 
 .. _ensemble statistics:
 
diff --git a/environment.yml b/environment.yml
index 94d156f0e3..cfff47377f 100644
--- a/environment.yml
+++ b/environment.yml
@@ -17,7 +17,7 @@ dependencies:
   - geopy
   - humanfriendly
   - importlib_resources
-  - iris>=3.2.1
+  - iris>=3.4.0
   - iris-esmf-regrid
   - isodate
   - jinja2
diff --git a/esmvalcore/cmor/_fixes/icon/_base_fixes.py b/esmvalcore/cmor/_fixes/icon/_base_fixes.py
index ecb722109e..7fe0ab1c1e 100644
--- a/esmvalcore/cmor/_fixes/icon/_base_fixes.py
+++ b/esmvalcore/cmor/_fixes/icon/_base_fixes.py
@@ -88,6 +88,7 @@ def get_horizontal_grid(self, cube):
         # Download file if necessary
         logger.debug("Attempting to download ICON grid file from '%s' to '%s'",
                      grid_url, grid_path)
+        self.CACHE_DIR.mkdir(parents=True, exist_ok=True)
         with requests.get(grid_url, stream=True,
                           timeout=self.TIMEOUT) as response:
             response.raise_for_status()
diff --git a/esmvalcore/cmor/_fixes/icon/icon.py b/esmvalcore/cmor/_fixes/icon/icon.py
index 8d36eac222..b65bb8b5ae 100644
--- a/esmvalcore/cmor/_fixes/icon/icon.py
+++ b/esmvalcore/cmor/_fixes/icon/icon.py
@@ -11,6 +11,7 @@
 from iris import NameConstraint
 from iris.coords import AuxCoord, DimCoord
 from iris.cube import CubeList
+from iris.experimental.ugrid import Connectivity, Mesh
 
 from esmvalcore.iris_helpers import add_leading_dim_to_cube, date2num
 
@@ -63,9 +64,9 @@ def fix_metadata(self, cubes):
         else:
             lon_idx = None
 
-        # Fix cell index for unstructured grid if necessary
-        if self._cell_index_needs_fixing(lat_idx, lon_idx):
-            self._fix_unstructured_cell_index(cube, lat_idx)
+        # Fix unstructured mesh of unstructured grid if present
+        if self._is_unstructured_grid(lat_idx, lon_idx):
+            self._fix_mesh(cube, lat_idx)
 
         # Fix scalar coordinates
         self.fix_scalar_coords(cube)
@@ -75,8 +76,7 @@ def fix_metadata(self, cubes):
 
         return CubeList([cube])
 
-    def _add_coord_from_grid_file(self, cube, coord_name,
-                                  target_coord_long_name):
+    def _add_coord_from_grid_file(self, cube, coord_name):
         """Add coordinate from grid file to cube.
 
         Note
@@ -89,10 +89,8 @@ def _add_coord_from_grid_file(self, cube, coord_name,
         cube: iris.cube.Cube
             ICON data to which the coordinate from the grid file is added.
         coord_name: str
-            Name of the coordinate in the grid file. Must be one of
-            ``'grid_latitude'``, ``'grid_longitude'``.
-        target_coord_long_name: str
-            Long name that is assigned to the newly added coordinate.
+            Name of the coordinate to add from the grid file. Must be one of
+            ``'latitude'``, ``'longitude'``.
 
         Raises
         ------
@@ -102,35 +100,43 @@ def _add_coord_from_grid_file(self, cube, coord_name,
             coordinate.
 
         """
-        allowed_coord_names = ('grid_latitude', 'grid_longitude')
-        if coord_name not in allowed_coord_names:
+        # The following dict maps from desired coordinate name in output file
+        # (dict keys) to coordinate name in grid file (dict values)
+        coord_names_mapping = {
+            'latitude': 'grid_latitude',
+            'longitude': 'grid_longitude',
+        }
+        if coord_name not in coord_names_mapping:
             raise ValueError(
-                f"coord_name must be one of {allowed_coord_names}, got "
+                f"coord_name must be one of {list(coord_names_mapping)}, got "
                 f"'{coord_name}'")
-        horizontal_grid = self.get_horizontal_grid(cube)
+        coord_name_in_grid = coord_names_mapping[coord_name]
 
-        # Use 'cell_area' as dummy cube to extract coordinates
+        # Use 'cell_area' as dummy cube to extract desired coordinates
         # Note: it might be necessary to expand this when more coord_names are
         # supported
+        horizontal_grid = self.get_horizontal_grid(cube)
         grid_cube = horizontal_grid.extract_cube(
             NameConstraint(var_name='cell_area'))
-        coord = grid_cube.coord(coord_name)
+        coord = grid_cube.coord(coord_name_in_grid)
 
-        # Find index of horizontal coordinate (= single unnamed dimension)
+        # Find index of mesh dimension (= single unnamed dimension)
         n_unnamed_dimensions = cube.ndim - len(cube.dim_coords)
         if n_unnamed_dimensions != 1:
             raise ValueError(
                 f"Cannot determine coordinate dimension for coordinate "
-                f"'{target_coord_long_name}', cube does not contain a single "
-                f"unnamed dimension:\n{cube}")
+                f"'{coord_name}', cube does not contain a single unnamed "
+                f"dimension:\n{cube}")
         coord_dims = ()
         for idx in range(cube.ndim):
             if not cube.coords(dimensions=idx, dim_coords=True):
                 coord_dims = (idx,)
                 break
 
+        # Adapt coordinate names so that the coordinate can be referenced with
+        # 'cube.coord(coord_name)'; the exact name will be set at a later stage
         coord.standard_name = None
-        coord.long_name = target_coord_long_name
+        coord.long_name = coord_name
         cube.add_aux_coord(coord, coord_dims)
 
     def _add_time(self, cube, cubes):
@@ -210,7 +216,7 @@ def _fix_lat(self, cube):
         # Add latitude coordinate if not already present
         if not cube.coords(lat_name):
             try:
-                self._add_coord_from_grid_file(cube, 'grid_latitude', lat_name)
+                self._add_coord_from_grid_file(cube, 'latitude')
             except Exception as exc:
                 msg = "Failed to add missing latitude coordinate to cube"
                 raise ValueError(msg) from exc
@@ -227,14 +233,14 @@ def _fix_lon(self, cube):
         # Add longitude coordinate if not already present
         if not cube.coords(lon_name):
             try:
-                self._add_coord_from_grid_file(
-                    cube, 'grid_longitude', lon_name)
+                self._add_coord_from_grid_file(cube, 'longitude')
             except Exception as exc:
                 msg = "Failed to add missing longitude coordinate to cube"
                 raise ValueError(msg) from exc
 
-        # Fix metadata
+        # Fix metadata and convert to [0, 360]
         lon = self.fix_lon_metadata(cube, lon_name)
+        self._set_range_in_0_360(lon)
 
         return cube.coord_dims(lon)
 
@@ -274,42 +280,197 @@ def _fix_time(self, cube, cubes):
 
         return cube
 
+    def _get_mesh(self, cube):
+        """Create mesh from horizontal grid file.
+
+        Note
+        ----
+        This functions creates a new :class:`iris.experimental.ugrid.Mesh` from
+        the ``clat`` (already present in the cube), ``clon`` (already present
+        in the cube), ``vertex_index``, ``vertex_of_cell``, ``vlat``, and
+        ``vlon`` variables of the horizontal grid file.
+
+        We do not use :func:`iris.experimental.ugrid.Mesh.from_coords` with the
+        existing latitude and longitude coordinates here because this would
+        produce lots of duplicated entries for the node coordinates. The reason
+        for this is that the node coordinates are constructed from the bounds;
+        since each node is contained 6 times in the bounds array (each node is
+        shared by 6 neighboring cells) the number of nodes is 6 times higher
+        with :func:`iris.experimental.ugrid.Mesh.from_coords` compared to using
+        the information already present in the horizontal grid file.
+
+        """
+        horizontal_grid = self.get_horizontal_grid(cube)
+
+        # Extract connectivity (i.e., the mapping cell faces -> cell nodes)
+        # from the the horizontal grid file (in ICON jargon called
+        # 'vertex_of_cell'; since UGRID expects a different dimension ordering
+        # we transpose the cube here)
+        vertex_of_cell = horizontal_grid.extract_cube(
+            NameConstraint(var_name='vertex_of_cell'))
+        vertex_of_cell.transpose()
+
+        # Extract start index used to name nodes from the the horizontal grid
+        # file
+        start_index = self._get_start_index(horizontal_grid)
+
+        # Extract face coordinates from cube (in ICON jargon called 'cell
+        # latitude' and 'cell longitude')
+        face_lat = cube.coord('latitude')
+        face_lon = cube.coord('longitude')
+
+        # Extract node coordinates from horizontal grid
+        (node_lat, node_lon) = self._get_node_coords(horizontal_grid)
+
+        # The bounds given by the face coordinates slightly differ from the
+        # bounds determined by the connectivity. We arbitrarily assume here
+        # that the information given by the connectivity is correct.
+        conn_node_inds = vertex_of_cell.data - start_index
+
+        # Latitude: there might be slight numerical differences (-> check that
+        # the differences are very small before fixing it)
+        if not np.allclose(face_lat.bounds, node_lat.points[conn_node_inds]):
+            raise ValueError(
+                "Cannot create mesh from horizontal grid file: latitude "
+                "bounds of the face coordinate ('clat_vertices' in the grid "
+                "file) differ from the corresponding values calculated from "
+                "the connectivity ('vertex_of_cell') and the node coordinate "
+                "('vlat')")
+        face_lat.bounds = node_lat.points[conn_node_inds]
+
+        # Longitude: there might be differences at the poles, where the
+        # longitude information does not matter (-> check that the only large
+        # differences are located at the poles). In addition, values might
+        # differ by 360°, which is also okay.
+        face_lon_bounds_to_check = face_lon.bounds % 360
+        node_lon_conn_to_check = node_lon.points[conn_node_inds] % 360
+        idx_notclose = ~np.isclose(face_lon_bounds_to_check,
+                                   node_lon_conn_to_check)
+        if not np.allclose(np.abs(face_lat.bounds[idx_notclose]), 90.0):
+            raise ValueError(
+                "Cannot create mesh from horizontal grid file: longitude "
+                "bounds of the face coordinate ('clon_vertices' in the grid "
+                "file) differ from the corresponding values calculated from "
+                "the connectivity ('vertex_of_cell') and the node coordinate "
+                "('vlon'). Note that these values are allowed to differ by "
+                "360° or at the poles of the grid.")
+        face_lon.bounds = node_lon.points[conn_node_inds]
+
+        # Create mesh
+        connectivity = Connectivity(
+            indices=vertex_of_cell.data,
+            cf_role='face_node_connectivity',
+            start_index=start_index,
+            location_axis=0,
+        )
+        mesh = Mesh(
+            topology_dimension=2,
+            node_coords_and_axes=[(node_lat, 'y'), (node_lon, 'x')],
+            connectivities=[connectivity],
+            face_coords_and_axes=[(face_lat, 'y'), (face_lon, 'x')],
+        )
+        return mesh
+
+    def _get_node_coords(self, horizontal_grid):
+        """Get node coordinates from horizontal grid.
+
+        Extract node coordinates from dummy variable 'dual_area' in horizontal
+        grid file (in ICON jargon called 'vertex latitude' and 'vertex
+        longitude'), remove their bounds (not accepted by UGRID), and adapt
+        metadata.
+
+        """
+        dual_area_cube = horizontal_grid.extract_cube(
+            NameConstraint(var_name='dual_area'))
+        node_lat = dual_area_cube.coord(var_name='vlat')
+        node_lon = dual_area_cube.coord(var_name='vlon')
+
+        # Fix metadata
+        node_lat.bounds = None
+        node_lon.bounds = None
+        node_lat.var_name = 'nlat'
+        node_lon.var_name = 'nlon'
+        node_lat.standard_name = 'latitude'
+        node_lon.standard_name = 'longitude'
+        node_lat.long_name = 'node latitude'
+        node_lon.long_name = 'node longitude'
+        node_lat.convert_units('degrees_north')
+        node_lon.convert_units('degrees_east')
+
+        # Convert longitude to [0, 360]
+        self._set_range_in_0_360(node_lon)
+
+        return (node_lat, node_lon)
+
+    def _fix_mesh(self, cube, mesh_idx):
+        """Fix mesh."""
+        # Remove any already-present dimensional coordinate describing the mesh
+        # dimension
+        if cube.coords(dimensions=mesh_idx, dim_coords=True):
+            cube.remove_coord(cube.coord(dimensions=mesh_idx, dim_coords=True))
+
+        # Add dimensional coordinate that describes the mesh dimension
+        index_coord = DimCoord(
+            np.arange(cube.shape[mesh_idx[0]]),
+            var_name='i',
+            long_name=('first spatial index for variables stored on an '
+                       'unstructured grid'),
+            units='1',
+        )
+        cube.add_dim_coord(index_coord, mesh_idx)
+
+        # Create mesh and replace the original latitude and longitude
+        # coordinates with their new mesh versions
+        mesh = self._get_mesh(cube)
+        cube.remove_coord('latitude')
+        cube.remove_coord('longitude')
+        for mesh_coord in mesh.to_MeshCoords('face'):
+            cube.add_aux_coord(mesh_coord, mesh_idx)
+
     @staticmethod
-    def _cell_index_needs_fixing(lat_idx, lon_idx):
-        """Check if cell index coordinate of unstructured grid needs fixing."""
+    def _get_start_index(horizontal_grid):
+        """Get start index used to name nodes from horizontal grid.
+
+        Extract start index used to name nodes from the the horizontal grid
+        file (in ICON jargon called 'vertex_index').
+
+        Note
+        ----
+        UGRID expects this to be a int32.
+
+        """
+        vertex_index = horizontal_grid.extract_cube(
+            NameConstraint(var_name='vertex_index'))
+        return np.int32(np.min(vertex_index.data))
+
+    @staticmethod
+    def _is_unstructured_grid(lat_idx, lon_idx):
+        """Check if data is defined on an unstructured grid."""
         # If either latitude or longitude are not present (i.e., the
-        # corresponding index is None), no fix is necessary
+        # corresponding index is None), no unstructured grid is present
         if lat_idx is None:
             return False
         if lon_idx is None:
             return False
 
-        # If latitude and longitude do not share their dimensions, no fix is
-        # necessary
+        # If latitude and longitude do not share their dimensions, no
+        # unstructured grid is present
         if lat_idx != lon_idx:
             return False
 
-        # If latitude and longitude are multi-dimensional (i.e., curvilinear
-        # instead of unstructured grid is given), no fix is necessary
+        # If latitude and longitude are multi-dimensional (e.g., curvilinear
+        # grid), no unstructured grid is present
         if len(lat_idx) != 1:
             return False
 
         return True
 
     @staticmethod
-    def _fix_unstructured_cell_index(cube, horizontal_idx):
-        """Fix unstructured cell index coordinate."""
-        if cube.coords(dimensions=horizontal_idx, dim_coords=True):
-            cube.remove_coord(cube.coord(dimensions=horizontal_idx,
-                                         dim_coords=True))
-        index_coord = DimCoord(
-            np.arange(cube.shape[horizontal_idx[0]]),
-            var_name='i',
-            long_name=('first spatial index for variables stored on an '
-                       'unstructured grid'),
-            units='1',
-        )
-        cube.add_dim_coord(index_coord, horizontal_idx)
+    def _set_range_in_0_360(lon_coord):
+        """Convert longitude coordinate to [0, 360]."""
+        lon_coord.points = (lon_coord.points + 360.0) % 360.0
+        if lon_coord.bounds is not None:
+            lon_coord.bounds = (lon_coord.bounds + 360.0) % 360.0
 
 
 class Clwvi(IconFix):
diff --git a/esmvalcore/preprocessor/_regrid.py b/esmvalcore/preprocessor/_regrid.py
index e4edbce9fe..c85f2d6434 100644
--- a/esmvalcore/preprocessor/_regrid.py
+++ b/esmvalcore/preprocessor/_regrid.py
@@ -478,14 +478,11 @@ def regrid(cube, target_grid, scheme, lat_offset=True, lon_offset=True):
     target_grid : Cube or str or dict
         The (location of a) cube that specifies the target or reference grid
         for the regridding operation.
-
         Alternatively, a string cell specification may be provided,
         of the form ``MxN``, which specifies the extent of the cell, longitude
         by latitude (degrees) for a global, regular target grid.
-
         Alternatively, a dictionary with a regional target grid may
         be specified (see above).
-
     scheme : str or dict
         The regridding scheme to perform. If both source and target grid are
         structured (regular or irregular), can be one of the built-in schemes
@@ -574,6 +571,7 @@ def regrid(cube, target_grid, scheme, lat_offset=True, lon_offset=True):
         raise ValueError(f'Expecting a cube, got {target_grid}.')
 
     if isinstance(scheme, dict):
+        scheme = dict(scheme)  # do not overwrite original scheme
         try:
             object_ref = scheme.pop("reference")
         except KeyError as key_err:
diff --git a/setup.py b/setup.py
index 7c55a74194..b36a60f248 100755
--- a/setup.py
+++ b/setup.py
@@ -55,7 +55,7 @@
         'pyyaml',
         'requests',
         'scipy>=1.6',
-        'scitools-iris>=3.2.1',
+        'scitools-iris>=3.4.0',
         'shapely[vectorized]',
         'stratify',
         'yamale',
diff --git a/tests/integration/cmor/_fixes/icon/test_icon.py b/tests/integration/cmor/_fixes/icon/test_icon.py
index e831a3a952..706fad7e06 100644
--- a/tests/integration/cmor/_fixes/icon/test_icon.py
+++ b/tests/integration/cmor/_fixes/icon/test_icon.py
@@ -9,11 +9,25 @@
 from iris.coords import AuxCoord, DimCoord
 from iris.cube import Cube, CubeList
 
+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 get_var_info
 from esmvalcore.config._config import get_extra_facets
 
+# TODO: ADAPT ONCE THE FILE IS AVAILABLE ON GITHUB
+TEST_GRID_FILE_URI = (
+    'https://seafile.zfn.uni-bremen.de/f/f66b000792434878bcbf/?dl=1'
+)
+TEST_GRID_FILE_NAME = 'f66b000792434878bcbf'
+
+
+@pytest.fixture(autouse=True)
+def tmp_cache_dir(monkeypatch, tmp_path):
+    """Use temporary path as cache directory for all tests in this module."""
+    monkeypatch.setattr(IconFix, 'CACHE_DIR', tmp_path)
+
+
 # Note: test_data_path is defined in tests/integration/cmor/_fixes/conftest.py
 
 
@@ -189,6 +203,7 @@ def check_lat(cube):
     assert lat.standard_name == 'latitude'
     assert lat.long_name == 'latitude'
     assert lat.units == 'degrees_north'
+    assert lat.attributes == {}
     np.testing.assert_allclose(
         lat.points,
         [-45.0, -45.0, -45.0, -45.0, 45.0, 45.0, 45.0, 45.0],
@@ -219,22 +234,23 @@ def check_lon(cube):
     assert lon.standard_name == 'longitude'
     assert lon.long_name == 'longitude'
     assert lon.units == 'degrees_east'
+    assert lon.attributes == {}
     np.testing.assert_allclose(
         lon.points,
-        [-135.0, -45.0, 45.0, 135.0, -135.0, -45.0, 45.0, 135.0],
+        [225.0, 315.0, 45.0, 135.0, 225.0, 315.0, 45.0, 135.0],
         rtol=1e-5
     )
     np.testing.assert_allclose(
         lon.bounds,
         [
-            [-135.0, -90.0, -180.0],
-            [-45.0, 0.0, -90.0],
-            [45.0, 90.0, 0.0],
-            [135.0, 180.0, 90.0],
-            [-180.0, -90.0, -135.0],
-            [-90.0, 0.0, -45.0],
-            [0.0, 90.0, 45.0],
-            [90.0, 180.0, 135.0],
+            [0.0, 270.0, 180.0],
+            [0.0, 0.0, 270.0],
+            [0.0, 90.0, 0.0],
+            [0.0, 180.0, 90.0],
+            [180.0, 270.0, 0.0],
+            [270.0, 0.0, 0.0],
+            [0.0, 90.0, 0.0],
+            [90.0, 180.0, 0.0],
         ],
         rtol=1e-5
     )
@@ -242,11 +258,15 @@ def check_lon(cube):
 
 
 def check_lat_lon(cube):
-    """Check latitude, longitude and spatial index coordinates of cube."""
+    """Check latitude, longitude and mesh of cube."""
     lat = check_lat(cube)
     lon = check_lon(cube)
 
-    # Check spatial index coordinate
+    # Check that latitude and longitude are mesh coordinates
+    assert cube.coords('latitude', mesh_coords=True)
+    assert cube.coords('longitude', mesh_coords=True)
+
+    # Check dimensional coordinate describing the mesh
     assert cube.coords('first spatial index for variables stored on an '
                        'unstructured grid', dim_coords=True)
     i_coord = cube.coord('first spatial index for variables stored on an '
@@ -263,6 +283,131 @@ def check_lat_lon(cube):
     assert cube.coord_dims(lat) == cube.coord_dims(lon)
     assert cube.coord_dims(lat) == cube.coord_dims(i_coord)
 
+    # Check the mesh itself
+    assert cube.location == 'face'
+    mesh = cube.mesh
+    check_mesh(mesh)
+
+
+def check_mesh(mesh):
+    """Check the mesh."""
+    assert mesh is not None
+    assert mesh.var_name is None
+    assert mesh.standard_name is None
+    assert mesh.long_name is None
+    assert mesh.units == 'unknown'
+    assert mesh.attributes == {}
+    assert mesh.cf_role == 'mesh_topology'
+    assert mesh.topology_dimension == 2
+
+    # Check face coordinates
+    assert len(mesh.coords(include_faces=True)) == 2
+
+    mesh_face_lat = mesh.coord(include_faces=True, axis='y')
+    assert mesh_face_lat.var_name == 'lat'
+    assert mesh_face_lat.standard_name == 'latitude'
+    assert mesh_face_lat.long_name == 'latitude'
+    assert mesh_face_lat.units == 'degrees_north'
+    assert mesh_face_lat.attributes == {}
+    np.testing.assert_allclose(
+        mesh_face_lat.points,
+        [-45.0, -45.0, -45.0, -45.0, 45.0, 45.0, 45.0, 45.0],
+        rtol=1e-5
+    )
+    np.testing.assert_allclose(
+        mesh_face_lat.bounds,
+        [
+            [-90.0, 0.0, 0.0],
+            [-90.0, 0.0, 0.0],
+            [-90.0, 0.0, 0.0],
+            [-90.0, 0.0, 0.0],
+            [0.0, 0.0, 90.0],
+            [0.0, 0.0, 90.0],
+            [0.0, 0.0, 90.0],
+            [0.0, 0.0, 90.0],
+        ],
+        rtol=1e-5
+    )
+
+    mesh_face_lon = mesh.coord(include_faces=True, axis='x')
+    assert mesh_face_lon.var_name == 'lon'
+    assert mesh_face_lon.standard_name == 'longitude'
+    assert mesh_face_lon.long_name == 'longitude'
+    assert mesh_face_lon.units == 'degrees_east'
+    assert mesh_face_lon.attributes == {}
+    np.testing.assert_allclose(
+        mesh_face_lon.points,
+        [225.0, 315.0, 45.0, 135.0, 225.0, 315.0, 45.0, 135.0],
+        rtol=1e-5
+    )
+    np.testing.assert_allclose(
+        mesh_face_lon.bounds,
+        [
+            [0.0, 270.0, 180.0],
+            [0.0, 0.0, 270.0],
+            [0.0, 90.0, 0.0],
+            [0.0, 180.0, 90.0],
+            [180.0, 270.0, 0.0],
+            [270.0, 0.0, 0.0],
+            [0.0, 90.0, 0.0],
+            [90.0, 180.0, 0.0],
+        ],
+        rtol=1e-5
+    )
+
+    # Check node coordinates
+    assert len(mesh.coords(include_nodes=True)) == 2
+
+    mesh_node_lat = mesh.coord(include_nodes=True, axis='y')
+    assert mesh_node_lat.var_name == 'nlat'
+    assert mesh_node_lat.standard_name == 'latitude'
+    assert mesh_node_lat.long_name == 'node latitude'
+    assert mesh_node_lat.units == 'degrees_north'
+    assert mesh_node_lat.attributes == {}
+    np.testing.assert_allclose(
+        mesh_node_lat.points,
+        [-90.0, 0.0, 0.0, 0.0, 0.0, 90.0],
+        rtol=1e-5
+    )
+    assert mesh_node_lat.bounds is None
+
+    mesh_node_lon = mesh.coord(include_nodes=True, axis='x')
+    assert mesh_node_lon.var_name == 'nlon'
+    assert mesh_node_lon.standard_name == 'longitude'
+    assert mesh_node_lon.long_name == 'node longitude'
+    assert mesh_node_lon.units == 'degrees_east'
+    assert mesh_node_lon.attributes == {}
+    np.testing.assert_allclose(
+        mesh_node_lon.points,
+        [0.0, 180.0, 270.0, 0.0, 90, 0.0],
+        rtol=1e-5
+    )
+    assert mesh_node_lon.bounds is None
+
+    # Check connectivity
+    assert len(mesh.connectivities()) == 1
+    conn = mesh.connectivity()
+    assert conn.var_name is None
+    assert conn.standard_name is None
+    assert conn.long_name is None
+    assert conn.units == 'unknown'
+    assert conn.attributes == {}
+    assert conn.cf_role == 'face_node_connectivity'
+    assert conn.start_index == 1
+    assert conn.location_axis == 0
+    assert conn.shape == (8, 3)
+    np.testing.assert_array_equal(
+        conn.indices,
+        [[1, 3, 2],
+         [1, 4, 3],
+         [1, 5, 4],
+         [1, 2, 5],
+         [2, 3, 6],
+         [3, 4, 6],
+         [4, 5, 6],
+         [5, 2, 6]],
+    )
+
 
 def check_typesi(cube):
     """Check scalar typesi coordinate of cube."""
@@ -761,22 +906,14 @@ def test_add_time_fail():
         fix._add_time(cube, cubes)
 
 
-def test_add_latitude(cubes_2d, tmp_path):
+def test_add_latitude(cubes_2d):
     """Test fix."""
     # Remove latitude from tas cube to test automatic addition
     tas_cube = cubes_2d.extract_cube(NameConstraint(var_name='tas'))
     tas_cube.remove_coord('latitude')
-    tas_cube.attributes['grid_file_uri'] = (
-        'https://github.com/ESMValGroup/ESMValCore/raw/main/tests/'
-        'integration/cmor/_fixes/test_data/icon_grid.nc'
-    )
     cubes = CubeList([tas_cube])
     fix = get_allvars_fix('Amon', 'tas')
 
-    # Temporary overwrite default cache location for downloads
-    original_cache_dir = fix.CACHE_DIR
-    fix.CACHE_DIR = tmp_path
-
     assert len(fix._horizontal_grids) == 0
     fixed_cubes = fix.fix_metadata(cubes)
 
@@ -784,28 +921,17 @@ def test_add_latitude(cubes_2d, tmp_path):
     assert cube.shape == (1, 8)
     check_lat_lon(cube)
     assert len(fix._horizontal_grids) == 1
-    assert 'icon_grid.nc' in fix._horizontal_grids
-
-    # Restore cache location
-    fix.CACHE_DIR = original_cache_dir
+    assert TEST_GRID_FILE_NAME in fix._horizontal_grids
 
 
-def test_add_longitude(cubes_2d, tmp_path):
+def test_add_longitude(cubes_2d):
     """Test fix."""
     # Remove longitude from tas cube to test automatic addition
     tas_cube = cubes_2d.extract_cube(NameConstraint(var_name='tas'))
     tas_cube.remove_coord('longitude')
-    tas_cube.attributes['grid_file_uri'] = (
-        'https://github.com/ESMValGroup/ESMValCore/raw/main/tests/'
-        'integration/cmor/_fixes/test_data/icon_grid.nc'
-    )
     cubes = CubeList([tas_cube])
     fix = get_allvars_fix('Amon', 'tas')
 
-    # Temporary overwrite default cache location for downloads
-    original_cache_dir = fix.CACHE_DIR
-    fix.CACHE_DIR = tmp_path
-
     assert len(fix._horizontal_grids) == 0
     fixed_cubes = fix.fix_metadata(cubes)
 
@@ -813,29 +939,18 @@ def test_add_longitude(cubes_2d, tmp_path):
     assert cube.shape == (1, 8)
     check_lat_lon(cube)
     assert len(fix._horizontal_grids) == 1
-    assert 'icon_grid.nc' in fix._horizontal_grids
-
-    # Restore cache location
-    fix.CACHE_DIR = original_cache_dir
+    assert TEST_GRID_FILE_NAME in fix._horizontal_grids
 
 
-def test_add_latitude_longitude(cubes_2d, tmp_path):
+def test_add_latitude_longitude(cubes_2d):
     """Test fix."""
     # Remove latitude and longitude from tas cube to test automatic addition
     tas_cube = cubes_2d.extract_cube(NameConstraint(var_name='tas'))
     tas_cube.remove_coord('latitude')
     tas_cube.remove_coord('longitude')
-    tas_cube.attributes['grid_file_uri'] = (
-        'https://github.com/ESMValGroup/ESMValCore/raw/main/tests/'
-        'integration/cmor/_fixes/test_data/icon_grid.nc'
-    )
     cubes = CubeList([tas_cube])
     fix = get_allvars_fix('Amon', 'tas')
 
-    # Temporary overwrite default cache location for downloads
-    original_cache_dir = fix.CACHE_DIR
-    fix.CACHE_DIR = tmp_path
-
     assert len(fix._horizontal_grids) == 0
     fixed_cubes = fix.fix_metadata(cubes)
 
@@ -843,17 +958,16 @@ def test_add_latitude_longitude(cubes_2d, tmp_path):
     assert cube.shape == (1, 8)
     check_lat_lon(cube)
     assert len(fix._horizontal_grids) == 1
-    assert 'icon_grid.nc' in fix._horizontal_grids
-
-    # Restore cache location
-    fix.CACHE_DIR = original_cache_dir
+    assert TEST_GRID_FILE_NAME in fix._horizontal_grids
 
 
 def test_add_latitude_fail(cubes_2d):
     """Test fix."""
-    # Remove latitude from tas cube to test automatic addition
+    # Remove latitude and grid file attribute from tas cube to test automatic
+    # addition
     tas_cube = cubes_2d.extract_cube(NameConstraint(var_name='tas'))
     tas_cube.remove_coord('latitude')
+    tas_cube.attributes = {}
     cubes = CubeList([tas_cube])
     fix = get_allvars_fix('Amon', 'tas')
 
@@ -864,9 +978,11 @@ def test_add_latitude_fail(cubes_2d):
 
 def test_add_longitude_fail(cubes_2d):
     """Test fix."""
-    # Remove longitude from tas cube to test automatic addition
+    # Remove longitude and grid file attribute from tas cube to test automatic
+    # addition
     tas_cube = cubes_2d.extract_cube(NameConstraint(var_name='tas'))
     tas_cube.remove_coord('longitude')
+    tas_cube.attributes = {}
     cubes = CubeList([tas_cube])
     fix = get_allvars_fix('Amon', 'tas')
 
@@ -881,8 +997,7 @@ def test_add_coord_from_grid_file_fail_invalid_coord():
 
     msg = r"coord_name must be one of .* got 'invalid_coord_name'"
     with pytest.raises(ValueError, match=msg):
-        fix._add_coord_from_grid_file(mock.sentinel.cube, 'invalid_coord_name',
-                                      'invalid_target_name')
+        fix._add_coord_from_grid_file(mock.sentinel.cube, 'invalid_coord_name')
 
 
 def test_add_coord_from_grid_file_fail_no_url():
@@ -892,58 +1007,36 @@ def test_add_coord_from_grid_file_fail_no_url():
     msg = ("Cube does not contain the attribute 'grid_file_uri' necessary to "
            "download the ICON horizontal grid file")
     with pytest.raises(ValueError, match=msg):
-        fix._add_coord_from_grid_file(Cube(0), 'grid_latitude', 'latitude')
+        fix._add_coord_from_grid_file(Cube(0), 'latitude')
 
 
-def test_add_coord_from_grid_fail_no_unnamed_dim(cubes_2d, tmp_path):
+def test_add_coord_from_grid_fail_no_unnamed_dim(cubes_2d):
     """Test fix."""
     # Remove latitude from tas cube to test automatic addition
     tas_cube = cubes_2d.extract_cube(NameConstraint(var_name='tas'))
     tas_cube.remove_coord('latitude')
-    tas_cube.attributes['grid_file_uri'] = (
-        'https://github.com/ESMValGroup/ESMValCore/raw/main/tests/'
-        'integration/cmor/_fixes/test_data/icon_grid.nc'
-    )
     index_coord = DimCoord(np.arange(8), var_name='ncells')
     tas_cube.add_dim_coord(index_coord, 1)
     fix = get_allvars_fix('Amon', 'tas')
 
-    # Temporary overwrite default cache location for downloads
-    original_cache_dir = fix.CACHE_DIR
-    fix.CACHE_DIR = tmp_path
-
     msg = ("Cannot determine coordinate dimension for coordinate 'latitude', "
            "cube does not contain a single unnamed dimension")
     with pytest.raises(ValueError, match=msg):
-        fix._add_coord_from_grid_file(tas_cube, 'grid_latitude', 'latitude')
-
-    # Restore cache location
-    fix.CACHE_DIR = original_cache_dir
+        fix._add_coord_from_grid_file(tas_cube, 'latitude')
 
 
-def test_add_coord_from_grid_fail_two_unnamed_dims(cubes_2d, tmp_path):
+def test_add_coord_from_grid_fail_two_unnamed_dims(cubes_2d):
     """Test fix."""
     # Remove latitude from tas cube to test automatic addition
     tas_cube = cubes_2d.extract_cube(NameConstraint(var_name='tas'))
     tas_cube.remove_coord('latitude')
-    tas_cube.attributes['grid_file_uri'] = (
-        'https://github.com/ESMValGroup/ESMValCore/raw/main/tests/'
-        'integration/cmor/_fixes/test_data/icon_grid.nc'
-    )
     tas_cube = iris.util.new_axis(tas_cube)
     fix = get_allvars_fix('Amon', 'tas')
 
-    # Temporary overwrite default cache location for downloads
-    original_cache_dir = fix.CACHE_DIR
-    fix.CACHE_DIR = tmp_path
-
     msg = ("Cannot determine coordinate dimension for coordinate 'latitude', "
            "cube does not contain a single unnamed dimension")
     with pytest.raises(ValueError, match=msg):
-        fix._add_coord_from_grid_file(tas_cube, 'grid_latitude', 'latitude')
-
-    # Restore cache location
-    fix.CACHE_DIR = original_cache_dir
+        fix._add_coord_from_grid_file(tas_cube, 'latitude')
 
 
 @mock.patch('esmvalcore.cmor._fixes.icon._base_fixes.requests', autospec=True)
@@ -970,10 +1063,6 @@ def test_get_horizontal_grid_cached_in_file(mock_requests, tmp_path):
     grid_cube = Cube(0, var_name='grid')
     iris.save(grid_cube, str(tmp_path / 'grid_file.nc'))
 
-    # Temporary overwrite default cache location for downloads
-    original_cache_dir = fix.CACHE_DIR
-    fix.CACHE_DIR = tmp_path
-
     grid = fix.get_horizontal_grid(cube)
     assert isinstance(grid, CubeList)
     assert len(grid) == 1
@@ -982,16 +1071,10 @@ def test_get_horizontal_grid_cached_in_file(mock_requests, tmp_path):
     assert 'grid_file.nc' in fix._horizontal_grids
     assert mock_requests.mock_calls == []
 
-    # Restore cache location
-    fix.CACHE_DIR = original_cache_dir
-
 
-def test_get_horizontal_grid_cache_file_too_old(tmp_path):
+def test_get_horizontal_grid_cache_file_too_old(tmp_path, monkeypatch):
     """Test fix."""
-    cube = Cube(0, attributes={
-        'grid_file_uri': 'https://github.com/ESMValGroup/ESMValCore/raw/main/'
-                         'tests/integration/cmor/_fixes/test_data/'
-                         'icon_grid.nc'})
+    cube = Cube(0, attributes={'grid_file_uri': TEST_GRID_FILE_URI})
     fix = get_allvars_fix('Amon', 'tas')
     assert len(fix._horizontal_grids) == 0
 
@@ -1001,35 +1084,31 @@ def test_get_horizontal_grid_cache_file_too_old(tmp_path):
 
     # Temporary overwrite default cache location for downloads and cache
     # validity duration
-    original_cache_dir = fix.CACHE_DIR
-    original_cache_validity = fix.CACHE_VALIDITY
-    fix.CACHE_DIR = tmp_path
-    fix.CACHE_VALIDITY = -1
+    monkeypatch.setattr(fix, 'CACHE_VALIDITY', -1)
 
     grid = fix.get_horizontal_grid(cube)
     assert isinstance(grid, CubeList)
-    assert len(grid) == 1
-    assert grid[0].var_name == 'cell_area'
+    assert len(grid) == 4
+    var_names = [cube.var_name for cube in grid]
+    assert 'cell_area' in var_names
+    assert 'dual_area' in var_names
+    assert 'vertex_index' in var_names
+    assert 'vertex_of_cell' in var_names
     assert len(fix._horizontal_grids) == 1
-    assert 'icon_grid.nc' in fix._horizontal_grids
-
-    # Restore cache location
-    fix.CACHE_DIR = original_cache_dir
-    fix.CACHE_VALIDITY = original_cache_validity
+    assert TEST_GRID_FILE_NAME in fix._horizontal_grids
 
 
 # Test with single-dimension cubes
 
 
-def test_only_time():
+def test_only_time(monkeypatch):
     """Test fix."""
     # We know that ta has dimensions time, plev19, latitude, longitude, but the
     # ICON CMORizer is designed to check for the presence of each dimension
     # individually. To test this, remove all but one dimension of ta to create
     # an artificial, but realistic test case.
     vardef = get_var_info('ICON', 'Amon', 'ta')
-    original_dimensions = vardef.dimensions
-    vardef.dimensions = ['time']
+    monkeypatch.setattr(vardef, 'dimensions', ['time'])
     extra_facets = get_extra_facets('ICON', 'ICON', 'Amon', 'ta', ())
     fix = AllVars(vardef, extra_facets=extra_facets)
 
@@ -1062,19 +1141,18 @@ def test_only_time():
     np.testing.assert_allclose(new_time_coord.bounds,
                                [[-0.5, 0.5], [0.5, 1.5]])
 
-    # Restore original dimensions of ta
-    vardef.dimensions = original_dimensions
+    # Check that no mesh has been created
+    assert cube.mesh is None
 
 
-def test_only_height():
+def test_only_height(monkeypatch):
     """Test fix."""
     # We know that ta has dimensions time, plev19, latitude, longitude, but the
     # ICON CMORizer is designed to check for the presence of each dimension
     # individually. To test this, remove all but one dimension of ta to create
     # an artificial, but realistic test case.
     vardef = get_var_info('ICON', 'Amon', 'ta')
-    original_dimensions = vardef.dimensions
-    vardef.dimensions = ['plev19']
+    monkeypatch.setattr(vardef, 'dimensions', ['plev19'])
     extra_facets = get_extra_facets('ICON', 'ICON', 'Amon', 'ta', ())
     fix = AllVars(vardef, extra_facets=extra_facets)
 
@@ -1107,19 +1185,18 @@ def test_only_height():
     np.testing.assert_allclose(new_height_coord.points, [1.0, 10.0])
     assert new_height_coord.bounds is None
 
-    # Restore original dimensions of ta
-    vardef.dimensions = original_dimensions
+    # Check that no mesh has been created
+    assert cube.mesh is None
 
 
-def test_only_latitude():
+def test_only_latitude(monkeypatch):
     """Test fix."""
     # We know that ta has dimensions time, plev19, latitude, longitude, but the
     # ICON CMORizer is designed to check for the presence of each dimension
     # individually. To test this, remove all but one dimension of ta to create
     # an artificial, but realistic test case.
     vardef = get_var_info('ICON', 'Amon', 'ta')
-    original_dimensions = vardef.dimensions
-    vardef.dimensions = ['latitude']
+    monkeypatch.setattr(vardef, 'dimensions', ['latitude'])
     extra_facets = get_extra_facets('ICON', 'ICON', 'Amon', 'ta', ())
     fix = AllVars(vardef, extra_facets=extra_facets)
 
@@ -1151,19 +1228,18 @@ def test_only_latitude():
     np.testing.assert_allclose(new_lat_coord.points, [0.0, 10.0])
     assert new_lat_coord.bounds is None
 
-    # Restore original dimensions of ta
-    vardef.dimensions = original_dimensions
+    # Check that no mesh has been created
+    assert cube.mesh is None
 
 
-def test_only_longitude():
+def test_only_longitude(monkeypatch):
     """Test fix."""
     # We know that ta has dimensions time, plev19, latitude, longitude, but the
     # ICON CMORizer is designed to check for the presence of each dimension
     # individually. To test this, remove all but one dimension of ta to create
     # an artificial, but realistic test case.
     vardef = get_var_info('ICON', 'Amon', 'ta')
-    original_dimensions = vardef.dimensions
-    vardef.dimensions = ['longitude']
+    monkeypatch.setattr(vardef, 'dimensions', ['longitude'])
     extra_facets = get_extra_facets('ICON', 'ICON', 'Amon', 'ta', ())
     fix = AllVars(vardef, extra_facets=extra_facets)
 
@@ -1195,8 +1271,8 @@ def test_only_longitude():
     np.testing.assert_allclose(new_lon_coord.points, [0.0, 180.0])
     assert new_lon_coord.bounds is None
 
-    # Restore original dimensions of ta
-    vardef.dimensions = original_dimensions
+    # Check that no mesh has been created
+    assert cube.mesh is None
 
 
 # Test variable not available in file
@@ -1221,3 +1297,38 @@ def test_invalid_time_units(cubes_2d):
     msg = "Expected time units"
     with pytest.raises(ValueError, match=msg):
         fix.fix_metadata(cubes_2d)
+
+
+# Test mesh creation fails because bounds do not match vertices
+
+
+def test_get_mesh_fail_invalid_clat_bounds(cubes_2d):
+    """Test fix."""
+    # Slightly modify latitude bounds from tas cube to make mesh creation fail
+    tas_cube = cubes_2d.extract_cube(NameConstraint(var_name='tas'))
+    lat_bnds = tas_cube.coord('latitude').bounds.copy()
+    lat_bnds[0, 0] = 40.0
+    tas_cube.coord('latitude').bounds = lat_bnds
+    cubes = CubeList([tas_cube])
+    fix = get_allvars_fix('Amon', 'tas')
+
+    msg = ("Cannot create mesh from horizontal grid file: latitude bounds of "
+           "the face coordinate")
+    with pytest.raises(ValueError, match=msg):
+        fix.fix_metadata(cubes)
+
+
+def test_get_mesh_fail_invalid_clon_bounds(cubes_2d):
+    """Test fix."""
+    # Slightly modify longitude bounds from tas cube to make mesh creation fail
+    tas_cube = cubes_2d.extract_cube(NameConstraint(var_name='tas'))
+    lon_bnds = tas_cube.coord('longitude').bounds.copy()
+    lon_bnds[0, 1] = 40.0
+    tas_cube.coord('longitude').bounds = lon_bnds
+    cubes = CubeList([tas_cube])
+    fix = get_allvars_fix('Amon', 'tas')
+
+    msg = ("Cannot create mesh from horizontal grid file: longitude bounds "
+           "of the face coordinate")
+    with pytest.raises(ValueError, match=msg):
+        fix.fix_metadata(cubes)
diff --git a/tests/integration/cmor/_fixes/test_data/icon_2d.nc b/tests/integration/cmor/_fixes/test_data/icon_2d.nc
index b6b7d66fe7..f650492ebd 100644
Binary files a/tests/integration/cmor/_fixes/test_data/icon_2d.nc and b/tests/integration/cmor/_fixes/test_data/icon_2d.nc differ
diff --git a/tests/integration/cmor/_fixes/test_data/icon_3d.nc b/tests/integration/cmor/_fixes/test_data/icon_3d.nc
index 4c4e5fd911..a3f03badff 100644
Binary files a/tests/integration/cmor/_fixes/test_data/icon_3d.nc and b/tests/integration/cmor/_fixes/test_data/icon_3d.nc differ
diff --git a/tests/integration/cmor/_fixes/test_data/icon_grid.nc b/tests/integration/cmor/_fixes/test_data/icon_grid.nc
index 882f63cedc..0e4a0634eb 100644
Binary files a/tests/integration/cmor/_fixes/test_data/icon_grid.nc and b/tests/integration/cmor/_fixes/test_data/icon_grid.nc differ