diff --git a/environment.yml b/environment.yml
index cacd62710b..53d0599f3d 100644
--- a/environment.yml
+++ b/environment.yml
@@ -25,6 +25,7 @@ dependencies:
   - jinja2
   - libnetcdf !=4.9.1  # to avoid hdf5 warnings
   - nc-time-axis
+  - ncdata
   - nested-lookup
   - netcdf4
   - numpy !=1.24.3,<2.0.0  # avoid pulling 2.0.0rcX
diff --git a/esmvalcore/cmor/_fixes/cmip6/cesm2.py b/esmvalcore/cmor/_fixes/cmip6/cesm2.py
index 0c5c0eed94..a414bf6d95 100644
--- a/esmvalcore/cmor/_fixes/cmip6/cesm2.py
+++ b/esmvalcore/cmor/_fixes/cmip6/cesm2.py
@@ -1,9 +1,9 @@
 """Fixes for CESM2 model."""
 
-from shutil import copyfile
-
+import ncdata
+import ncdata.iris
+import ncdata.netcdf4
 import numpy as np
-from netCDF4 import Dataset
 
 from ..common import SiconcFixScalarCoord
 from ..fix import Fix
@@ -19,24 +19,16 @@
 class Cl(Fix):
     """Fixes for ``cl``."""
 
-    def _fix_formula_terms(
-        self,
-        filepath,
-        output_dir,
-        add_unique_suffix=False,
-    ):
+    @staticmethod
+    def _fix_formula_terms(dataset: ncdata.NcData) -> None:
         """Fix ``formula_terms`` attribute."""
-        new_path = self.get_fixed_filepath(
-            output_dir, filepath, add_unique_suffix=add_unique_suffix
+        lev = dataset.variables["lev"]
+        lev.set_attrval("formula_terms", "p0: p0 a: a b: b ps: ps")
+        lev.set_attrval(
+            "standard_name", "atmosphere_hybrid_sigma_pressure_coordinate"
         )
-        copyfile(filepath, new_path)
-        dataset = Dataset(new_path, mode="a")
-        dataset.variables["lev"].formula_terms = "p0: p0 a: a b: b ps: ps"
-        dataset.variables[
-            "lev"
-        ].standard_name = "atmosphere_hybrid_sigma_pressure_coordinate"
-        dataset.close()
-        return new_path
+        lev.set_attrval("units", "1")
+        dataset.variables["lev_bnds"].attributes.pop("units")
 
     def fix_file(self, filepath, output_dir, add_unique_suffix=False):
         """Fix hybrid pressure coordinate.
@@ -63,45 +55,33 @@ def fix_file(self, filepath, output_dir, add_unique_suffix=False):
         -------
         str
             Path to the fixed file.
-
         """
-        new_path = self._fix_formula_terms(
-            filepath, output_dir, add_unique_suffix=add_unique_suffix
+        dataset = ncdata.netcdf4.from_nc4(filepath)
+        self._fix_formula_terms(dataset)
+
+        # Correct order of bounds data
+        a_bnds = dataset.variables["a_bnds"]
+        a_bnds.data = a_bnds.data[::-1, :]
+        b_bnds = dataset.variables["b_bnds"]
+        b_bnds.data = b_bnds.data[::-1, :]
+
+        # Correct lev and lev_bnds data
+        lev = dataset.variables["lev"]
+        lev.data = dataset.variables["a"].data + dataset.variables["b"].data
+        lev_bnds = dataset.variables["lev_bnds"]
+        lev_bnds.data = (
+            dataset.variables["a_bnds"].data + dataset.variables["b_bnds"].data
         )
-        dataset = Dataset(new_path, mode="a")
-        dataset.variables["a_bnds"][:] = dataset.variables["a_bnds"][::-1, :]
-        dataset.variables["b_bnds"][:] = dataset.variables["b_bnds"][::-1, :]
-        dataset.close()
-        return new_path
-
-    def fix_metadata(self, cubes):
-        """Fix ``atmosphere_hybrid_sigma_pressure_coordinate``.
-
-        See discussion in #882 for more details on that.
-
-        Parameters
-        ----------
-        cubes : iris.cube.CubeList
-            Input cubes.
 
-        Returns
-        -------
-        iris.cube.CubeList
+        # Remove 'title' attribute that duplicates long name
+        for var_name in dataset.variables:
+            dataset.variables[var_name].attributes.pop("title", None)
 
-        """
-        cube = self.get_cube_from_list(cubes)
-        lev_coord = cube.coord(var_name="lev")
-        a_coord = cube.coord(var_name="a")
-        b_coord = cube.coord(var_name="b")
-        lev_coord.points = a_coord.core_points() + b_coord.core_points()
-        lev_coord.bounds = a_coord.core_bounds() + b_coord.core_bounds()
-        lev_coord.units = "1"
-        return cubes
+        return self.ncdata_to_iris(dataset, filepath)
 
 
 Cli = Cl
 
-
 Clw = Cl
 
 
@@ -119,7 +99,6 @@ def fix_metadata(self, cubes):
         Returns
         -------
         iris.cube.CubeList
-
         """
         cube = self.get_cube_from_list(cubes)
         add_scalar_depth_coord(cube)
@@ -130,8 +109,7 @@ class Prw(Fix):
     """Fixes for tas."""
 
     def fix_metadata(self, cubes):
-        """
-        Fix latitude_bounds and longitude_bounds data type and round to 4 d.p.
+        """Fix latitude_bounds and longitude_bounds dtype and round to 4 d.p.
 
         Parameters
         ----------
@@ -141,7 +119,6 @@ def fix_metadata(self, cubes):
         Returns
         -------
         iris.cube.CubeList
-
         """
         for cube in cubes:
             for coord_name in ["latitude", "longitude"]:
@@ -159,8 +136,7 @@ class Tas(Prw):
     """Fixes for tas."""
 
     def fix_metadata(self, cubes):
-        """
-        Add height (2m) coordinate.
+        """Add height (2m) coordinate.
 
         Fix also done for prw.
         Fix latitude_bounds and longitude_bounds data type and round to 4 d.p.
@@ -173,7 +149,6 @@ def fix_metadata(self, cubes):
         Returns
         -------
         iris.cube.CubeList
-
         """
         super().fix_metadata(cubes)
         # Specific code for tas
@@ -197,7 +172,6 @@ def fix_metadata(self, cubes):
         Returns
         -------
         iris.cube.CubeList
-
         """
         cube = self.get_cube_from_list(cubes)
         add_scalar_typeland_coord(cube)
@@ -218,7 +192,6 @@ def fix_metadata(self, cubes):
         Returns
         -------
         iris.cube.CubeList
-
         """
         cube = self.get_cube_from_list(cubes)
         add_scalar_typesea_coord(cube)
@@ -232,8 +205,7 @@ class Tos(Fix):
     """Fixes for tos."""
 
     def fix_metadata(self, cubes):
-        """
-        Round times to 1 d.p. for monthly means.
+        """Round times to 1 d.p. for monthly means.
 
         Required to get hist-GHG and ssp245-GHG Omon tos to concatenate.
 
@@ -245,7 +217,6 @@ def fix_metadata(self, cubes):
         Returns
         -------
         iris.cube.CubeList
-
         """
         cube = self.get_cube_from_list(cubes)
 
@@ -271,7 +242,6 @@ def fix_metadata(self, cubes):
         Returns
         -------
         iris.cube.CubeList
-
         """
         for cube in cubes:
             if cube.coords(axis="Z"):
diff --git a/esmvalcore/cmor/_fixes/cmip6/cesm2_waccm.py b/esmvalcore/cmor/_fixes/cmip6/cesm2_waccm.py
index 156a656b5f..a26f61f058 100644
--- a/esmvalcore/cmor/_fixes/cmip6/cesm2_waccm.py
+++ b/esmvalcore/cmor/_fixes/cmip6/cesm2_waccm.py
@@ -1,6 +1,7 @@
 """Fixes for CESM2-WACCM model."""
 
-from netCDF4 import Dataset
+import ncdata.iris
+import ncdata.netcdf4
 
 from ..common import SiconcFixScalarCoord
 from .cesm2 import Cl as BaseCl
@@ -10,7 +11,7 @@
 
 
 class Cl(BaseCl):
-    """Fixes for cl."""
+    """Fixes for ``cl``."""
 
     def fix_file(self, filepath, output_dir, add_unique_suffix=False):
         """Fix hybrid pressure coordinate.
@@ -37,31 +38,31 @@ def fix_file(self, filepath, output_dir, add_unique_suffix=False):
         -------
         str
             Path to the fixed file.
-
         """
-        new_path = self._fix_formula_terms(
-            filepath, output_dir, add_unique_suffix=add_unique_suffix
-        )
-        dataset = Dataset(new_path, mode="a")
-        dataset.variables["a_bnds"][:] = dataset.variables["a_bnds"][:, ::-1]
-        dataset.variables["b_bnds"][:] = dataset.variables["b_bnds"][:, ::-1]
-        dataset.close()
-        return new_path
+        dataset = ncdata.netcdf4.from_nc4(filepath)
+        self._fix_formula_terms(dataset)
 
+        # Correct order of bounds data
+        a_bnds = dataset.variables["a_bnds"]
+        a_bnds.data = a_bnds.data[:, ::-1]
+        b_bnds = dataset.variables["b_bnds"]
+        b_bnds.data = b_bnds.data[:, ::-1]
 
-Cli = Cl
+        # Remove 'title' attribute that duplicates long name
+        for var_name in dataset.variables:
+            dataset.variables[var_name].attributes.pop("title", None)
 
+        return self.ncdata_to_iris(dataset, filepath)
 
-Clw = Cl
 
+Cli = Cl
 
-Fgco2 = BaseFgco2
+Clw = Cl
 
+Fgco2 = BaseFgco2
 
 Omon = BaseOmon
 
-
 Siconc = SiconcFixScalarCoord
 
-
 Tas = BaseTas
diff --git a/esmvalcore/cmor/_fixes/fix.py b/esmvalcore/cmor/_fixes/fix.py
index 973ac57d0b..ff5af625b6 100644
--- a/esmvalcore/cmor/_fixes/fix.py
+++ b/esmvalcore/cmor/_fixes/fix.py
@@ -9,8 +9,12 @@
 from collections.abc import Sequence
 from pathlib import Path
 from typing import TYPE_CHECKING, Any, Optional
+from warnings import catch_warnings, filterwarnings
 
 import dask
+import iris
+import ncdata.iris
+import ncdata.threadlock_sharing
 import numpy as np
 from cf_units import Unit
 from iris.coords import Coord, CoordExtent
@@ -36,6 +40,9 @@
 logger = logging.getLogger(__name__)
 generic_fix_logger = logging.getLogger(f"{__name__}.genericfix")
 
+# Enable lock sharing between Iris and ncdata
+ncdata.threadlock_sharing.enable_lockshare(iris=True)
+
 
 class Fix:
     """Base class for dataset fixes."""
@@ -157,6 +164,60 @@ def get_cube_from_list(
                 return cube
         raise ValueError(f'Cube for variable "{short_name}" not found')
 
+    def ncdata_to_iris(
+        self,
+        dataset: ncdata.NcData,
+        filepath: Path,
+    ) -> CubeList:
+        """Convert an :obj:`~ncdata.NcData` object to an Iris cubelist.
+
+        This function mimics the behaviour of
+        :func:`esmvalcore.preprocessor.load`.
+
+        Parameters
+        ----------
+        dataset:
+            The :obj:`~ncdata.NcData` object to convert.
+        filepath:
+            The path that the dataset was loaded from.
+
+        Returns
+        -------
+        iris.cube.CubeList
+            :obj:`iris.cube.CubeList` containing the requested cube.
+
+        """
+        # Filter warnings
+        with catch_warnings():
+            # Ignore warnings about missing cell measures that are stored in
+            # a separate file for CMIP data.
+            filterwarnings(
+                message="Missing CF-netCDF measure variable .*",
+                category=UserWarning,
+                module="iris",
+                action="ignore",
+            )
+            cubes = ncdata.iris.to_iris(dataset)
+
+        cube = self.get_cube_from_list(cubes)
+
+        # Restore the lat/lon coordinate units that iris changes to degrees
+        for coord_name in ["latitude", "longitude"]:
+            try:
+                coord = cube.coord(coord_name)
+            except iris.exceptions.CoordinateNotFoundError:
+                pass
+            else:
+                if coord.var_name in dataset.variables:
+                    nc_coord = dataset.variables[coord.var_name]
+                    coord.units = nc_coord.attributes["units"].value
+
+        # Add the source file as an attribute to support grouping by file
+        # when calling fix_metadata.
+        cube.attributes["source_file"] = str(filepath)
+
+        return iris.cube.CubeList([cube])
+
     def fix_data(self, cube: Cube) -> Cube:
         """Apply fixes to the data of the cube.
 
diff --git a/esmvalcore/preprocessor/_io.py b/esmvalcore/preprocessor/_io.py
index 5f83b1946c..ea8fa90a8d 100644
--- a/esmvalcore/preprocessor/_io.py
+++ b/esmvalcore/preprocessor/_io.py
@@ -49,21 +49,15 @@ def _get_attr_from_field_coord(ncfield, coord_name, attr):
     return None
 
 
-def _load_callback(raw_cube, field, _):
-    """Use this callback to fix anything Iris tries to break."""
-    # Remove attributes that cause issues with merging and concatenation
-    _delete_attributes(
-        raw_cube, ("creation_date", "tracking_id", "history", "comment")
-    )
-    for coord in raw_cube.coords():
-        # Iris chooses to change longitude and latitude units to degrees
-        # regardless of value in file, so reinstating file value
+def _restore_lat_lon_units(cube, field, filename):  # pylint: disable=unused-argument
+    """Use this callback to restore the original lat/lon units."""
+    # Iris chooses to change longitude and latitude units to degrees
+    # regardless of value in file, so reinstating file value
+    for coord in cube.coords():
         if coord.standard_name in ["longitude", "latitude"]:
             units = _get_attr_from_field_coord(field, coord.var_name, "units")
             if units is not None:
                 coord.units = units
-        # CMOR sometimes adds a history to the coordinates.
-        _delete_attributes(coord, ("history",))
 
 
 def _delete_attributes(iris_object, atts):
@@ -73,7 +67,7 @@ def _delete_attributes(iris_object, atts):
 
 
 def load(
-    file: str | Path,
+    file: str | Path | iris.cube.Cube,
     ignore_warnings: Optional[list[dict]] = None,
 ) -> CubeList:
     """Load iris cubes from string or Path objects.
@@ -81,7 +75,8 @@ def load(
     Parameters
     ----------
     file:
-        File to be loaded. Could be string or POSIX Path object.
+        File to be loaded. If ``file`` is already an Iris Cube, it will be
+        put in a CubeList and returned.
     ignore_warnings:
         Keyword arguments passed to :func:`warnings.filterwarnings` used to
         ignore warnings issued by :func:`iris.load_raw`. Each list element
@@ -97,6 +92,9 @@ def load(
     ValueError
         Cubes are empty.
     """
+    if isinstance(file, iris.cube.Cube):
+        return iris.cube.CubeList([file])
+
     file = Path(file)
     logger.debug("Loading:\n%s", file)
 
@@ -139,7 +137,7 @@ def load(
         # warnings.filterwarnings
         # (see https://github.com/SciTools/cf-units/issues/240)
         with suppress_errors():
-            raw_cubes = iris.load_raw(file, callback=_load_callback)
+            raw_cubes = iris.load_raw(file, callback=_restore_lat_lon_units)
     logger.debug("Done with loading %s", file)
 
     if not raw_cubes:
@@ -386,6 +384,15 @@ def concatenate(cubes, check_level=CheckLevels.DEFAULT):
     if len(cubes) == 1:
         return cubes[0]
 
+    for cube in cubes:
+        # Remove attributes that cause issues with merging and concatenation
+        _delete_attributes(
+            cube, ("creation_date", "tracking_id", "history", "comment")
+        )
+        for coord in cube.coords():
+            # CMOR sometimes adds a history to the coordinates.
+            _delete_attributes(coord, ("history",))
+
     cubes = _concatenate_cubes_by_experiment(cubes)
 
     merge_cube_attributes(cubes)
diff --git a/pyproject.toml b/pyproject.toml
index 797e9bc81b..6f4ed9b8da 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -48,6 +48,7 @@ dependencies = [
     "isodate>=0.7.0",
     "jinja2",
     "nc-time-axis",  # needed by iris.plot
+    "ncdata",
     "nested-lookup",
     "netCDF4",
     "numpy!=1.24.3,<2.0.0",  # avoid pulling 2.0.0rc1
diff --git a/tests/integration/cmor/_fixes/cmip6/test_cesm2.py b/tests/integration/cmor/_fixes/cmip6/test_cesm2.py
index 5d504f6084..cd7ccdd046 100644
--- a/tests/integration/cmor/_fixes/cmip6/test_cesm2.py
+++ b/tests/integration/cmor/_fixes/cmip6/test_cesm2.py
@@ -1,8 +1,5 @@
 """Tests for the fixes of CESM2."""
 
-import os
-import unittest.mock
-
 import iris
 import numpy as np
 import pytest
@@ -60,10 +57,7 @@ def test_get_cl_fix():
 )
 
 
-@unittest.mock.patch(
-    "esmvalcore.cmor._fixes.cmip6.cesm2.Fix.get_fixed_filepath", autospec=True
-)
-def test_cl_fix_file(mock_get_filepath, tmp_path, test_data_path):
+def test_cl_fix_file(tmp_path, test_data_path):
     """Test ``fix_file`` for ``cl``."""
     nc_path = test_data_path / "cesm2_cl.nc"
     cubes = iris.load(str(nc_path))
@@ -82,22 +76,10 @@ def test_cl_fix_file(mock_get_filepath, tmp_path, test_data_path):
     assert not raw_cube.coords("air_pressure")
 
     # Apply fix
-    mock_get_filepath.return_value = os.path.join(
-        tmp_path, "fixed_cesm2_cl.nc"
-    )
-    fix = Cl(None)
-    fixed_file = fix.fix_file(nc_path, tmp_path)
-    mock_get_filepath.assert_called_once_with(
-        tmp_path, nc_path, add_unique_suffix=False
-    )
-    fixed_cubes = iris.load(fixed_file)
-    assert len(fixed_cubes) == 2
-    var_names = [cube.var_name for cube in fixed_cubes]
-    assert "cl" in var_names
-    assert "ps" in var_names
-    fixed_cl_cube = fixed_cubes.extract_cube(
-        "cloud_area_fraction_in_atmosphere_layer"
-    )
+    vardef = get_var_info("CMIP6", "Amon", "cl")
+    fix = Cl(vardef)
+    (fixed_cl_cube,) = fix.fix_file(nc_path, tmp_path)
+    assert fixed_cl_cube.var_name == "cl"
     fixed_air_pressure_coord = fixed_cl_cube.coord("air_pressure")
     assert fixed_air_pressure_coord.points is not None
     assert fixed_air_pressure_coord.bounds is not None
@@ -107,73 +89,10 @@ def test_cl_fix_file(mock_get_filepath, tmp_path, test_data_path):
     np.testing.assert_allclose(
         fixed_air_pressure_coord.bounds, AIR_PRESSURE_BOUNDS
     )
-
-
-@pytest.fixture
-def cl_cubes():
-    """``cl`` cube."""
-    time_coord = iris.coords.DimCoord(
-        [0.0, 1.0],
-        var_name="time",
-        standard_name="time",
-        units="days since 1850-01-01 00:00:00",
-    )
-    a_coord = iris.coords.AuxCoord(
-        [0.1, 0.2, 0.1],
-        bounds=[[0.0, 0.15], [0.15, 0.25], [0.25, 0.0]],
-        var_name="a",
-        units="1",
-    )
-    b_coord = iris.coords.AuxCoord(
-        [0.9, 0.3, 0.1],
-        bounds=[[1.0, 0.8], [0.8, 0.25], [0.25, 0.0]],
-        var_name="b",
-        units="1",
-    )
-    lev_coord = iris.coords.DimCoord(
-        [999.0, 99.0, 9.0],
-        var_name="lev",
-        standard_name="atmosphere_hybrid_sigma_pressure_coordinate",
-        units="hPa",
-        attributes={"positive": "up"},
-    )
-    lat_coord = iris.coords.DimCoord(
-        [0.0, 1.0], var_name="lat", standard_name="latitude", units="degrees"
-    )
-    lon_coord = iris.coords.DimCoord(
-        [0.0, 1.0], var_name="lon", standard_name="longitude", units="degrees"
-    )
-    coord_specs = [
-        (time_coord, 0),
-        (lev_coord, 1),
-        (lat_coord, 2),
-        (lon_coord, 3),
-    ]
-    cube = iris.cube.Cube(
-        np.arange(2 * 3 * 2 * 2).reshape(2, 3, 2, 2),
-        var_name="cl",
-        standard_name="cloud_area_fraction_in_atmosphere_layer",
-        units="%",
-        dim_coords_and_dims=coord_specs,
-        aux_coords_and_dims=[(a_coord, 1), (b_coord, 1)],
-    )
-    return iris.cube.CubeList([cube])
-
-
-def test_cl_fix_metadata(cl_cubes):
-    """Test ``fix_metadata`` for ``cl``."""
-    vardef = get_var_info("CMIP6", "Amon", "cl")
-    fix = Cl(vardef)
-    out_cubes = fix.fix_metadata(cl_cubes)
-    out_cube = out_cubes.extract_cube(
-        "cloud_area_fraction_in_atmosphere_layer"
-    )
-    lev_coord = out_cube.coord(var_name="lev")
+    lev_coord = fixed_cl_cube.coord(var_name="lev")
     assert lev_coord.units == "1"
-    np.testing.assert_allclose(lev_coord.points, [1.0, 0.5, 0.2])
-    np.testing.assert_allclose(
-        lev_coord.bounds, [[1.0, 0.95], [0.95, 0.5], [0.5, 0.0]]
-    )
+    np.testing.assert_allclose(lev_coord.points, [1, 3])
+    np.testing.assert_allclose(lev_coord.bounds, [[-1, 2], [2, 5]])
 
 
 def test_get_cli_fix():
diff --git a/tests/integration/cmor/_fixes/cmip6/test_cesm2_waccm.py b/tests/integration/cmor/_fixes/cmip6/test_cesm2_waccm.py
index 363bf0d80c..036e57bd7d 100644
--- a/tests/integration/cmor/_fixes/cmip6/test_cesm2_waccm.py
+++ b/tests/integration/cmor/_fixes/cmip6/test_cesm2_waccm.py
@@ -1,8 +1,5 @@
 """Tests for the fixes of CESM2-WACCM."""
 
-import os
-import unittest.mock
-
 import iris
 import numpy as np
 import pytest
@@ -22,6 +19,7 @@
 from esmvalcore.cmor._fixes.common import SiconcFixScalarCoord
 from esmvalcore.cmor._fixes.fix import GenericFix
 from esmvalcore.cmor.fix import Fix
+from esmvalcore.cmor.table import get_var_info
 
 
 def test_get_cl_fix():
@@ -35,21 +33,12 @@ def test_cl_fix():
     assert issubclass(Cl, BaseCl)
 
 
-@unittest.mock.patch(
-    "esmvalcore.cmor._fixes.cmip6.cesm2.Fix.get_fixed_filepath", autospec=True
-)
-def test_cl_fix_file(mock_get_filepath, tmp_path, test_data_path):
+def test_cl_fix_file(tmp_path, test_data_path):
     """Test ``fix_file`` for ``cl``."""
     nc_path = test_data_path / "cesm2_waccm_cl.nc"
-    mock_get_filepath.return_value = os.path.join(
-        tmp_path, "fixed_cesm2_waccm_cl.nc"
-    )
-    fix = Cl(None)
-    fixed_file = fix.fix_file(nc_path, tmp_path)
-    mock_get_filepath.assert_called_once_with(
-        tmp_path, nc_path, add_unique_suffix=False
-    )
-    fixed_cube = iris.load_cube(fixed_file)
+    vardef = get_var_info("CMIP6", "Amon", "cl")
+    fix = Cl(vardef)
+    (fixed_cube,) = fix.fix_file(nc_path, tmp_path)
     lev_coord = fixed_cube.coord(var_name="lev")
     a_coord = fixed_cube.coord(var_name="a")
     b_coord = fixed_cube.coord(var_name="b")
diff --git a/tests/integration/preprocessor/_io/test_concatenate.py b/tests/integration/preprocessor/_io/test_concatenate.py
index 1fef5f9693..f600e1d095 100644
--- a/tests/integration/preprocessor/_io/test_concatenate.py
+++ b/tests/integration/preprocessor/_io/test_concatenate.py
@@ -348,6 +348,26 @@ def test_concatenate_by_experiment_first(self):
         assert_array_equal(result.coord("time").points, np.arange(7))
         assert_array_equal(result.data, np.array([0, 0, 0, 1, 1, 1, 1]))
 
+    def test_concatenate_remove_unwanted_attributes(self):
+        """Test concatenate removes unwanted attributes."""
+        attributes = ("history", "creation_date", "tracking_id", "comment")
+        for i, cube in enumerate(self.raw_cubes):
+            for attr in attributes:
+                cube.attributes[attr] = f"{attr}-{i}"
+        concatenated = _io.concatenate(self.raw_cubes)
+        assert not set(attributes) & set(concatenated.attributes)
+
+    def test_concatenate_remove_unwanted_attributes_from_coords(self):
+        """Test concatenate removes unwanted attributes from coords."""
+        attributes = ("history",)
+        for i, cube in enumerate(self.raw_cubes):
+            for coord in cube.coords():
+                for attr in attributes:
+                    coord.attributes[attr] = f"{attr}-{i}"
+        concatenated = _io.concatenate(self.raw_cubes)
+        for coord in concatenated.coords():
+            assert not set(attributes) & set(coord.attributes)
+
     def test_concatenate_differing_attributes(self):
         """Test concatenation of cubes with different attributes."""
         cubes = CubeList(self.raw_cubes)
diff --git a/tests/integration/preprocessor/_io/test_load.py b/tests/integration/preprocessor/_io/test_load.py
index 4c76ba2651..4baffe7cec 100644
--- a/tests/integration/preprocessor/_io/test_load.py
+++ b/tests/integration/preprocessor/_io/test_load.py
@@ -52,45 +52,11 @@ def test_load(self):
             (cube.coord("latitude").points == np.array([1, 2])).all()
         )
 
-    def test_callback_remove_attributes(self):
-        """Test callback remove unwanted attributes."""
-        attributes = ("history", "creation_date", "tracking_id", "comment")
-        for _ in range(2):
-            cube = _create_sample_cube()
-            for attr in attributes:
-                cube.attributes[attr] = attr
-            self._save_cube(cube)
-        for temp_file in self.temp_files:
-            cubes = load(temp_file)
-            cube = cubes[0]
-            self.assertEqual(1, len(cubes))
-            self.assertTrue((cube.data == np.array([1, 2])).all())
-            self.assertTrue(
-                (cube.coord("latitude").points == np.array([1, 2])).all()
-            )
-            for attr in attributes:
-                self.assertTrue(attr not in cube.attributes)
-
-    def test_callback_remove_attributes_from_coords(self):
-        """Test callback remove unwanted attributes from coords."""
-        attributes = ("history",)
-        for _ in range(2):
-            cube = _create_sample_cube()
-            for coord in cube.coords():
-                for attr in attributes:
-                    coord.attributes[attr] = attr
-            self._save_cube(cube)
-        for temp_file in self.temp_files:
-            cubes = load(temp_file)
-            cube = cubes[0]
-            self.assertEqual(1, len(cubes))
-            self.assertTrue((cube.data == np.array([1, 2])).all())
-            self.assertTrue(
-                (cube.coord("latitude").points == np.array([1, 2])).all()
-            )
-            for coord in cube.coords():
-                for attr in attributes:
-                    self.assertTrue(attr not in coord.attributes)
+    def test_load_noop(self):
+        """Test loading an Iris cube."""
+        cube = _create_sample_cube()
+        cubes = load(cube)
+        assert [cube] == cubes
 
     def test_callback_fix_lat_units(self):
         """Test callback for fixing units."""