diff --git a/MANIFEST.in b/MANIFEST.in
index 99b801e827..c64ba3dd7b 100644
--- a/MANIFEST.in
+++ b/MANIFEST.in
@@ -5,7 +5,7 @@ include CHANGES COPYING COPYING.LESSER
 recursive-include lib/iris/tests/results *.cml *.cdl *.txt *.xml *.json
 recursive-include lib/iris/etc *
 include lib/iris/fileformats/_pyke_rules/*.k?b
-include lib/iris/tests/stock*.npz
+include lib/iris/tests/stock/file_headers/*
 
 include requirements/*.txt
 
diff --git a/lib/iris/tests/stock/__init__.py b/lib/iris/tests/stock/__init__.py
index 42ddfac161..d1941a7369 100644
--- a/lib/iris/tests/stock/__init__.py
+++ b/lib/iris/tests/stock/__init__.py
@@ -4,7 +4,7 @@
 # See COPYING and COPYING.LESSER in the root of the repository for full
 # licensing details.
 """
-A collection of routines which create standard Cubes for test purposes.
+A collection of routines which create standard Cubes/files for test purposes.
 
 """
 
diff --git a/lib/iris/tests/stock/file_headers/README.md b/lib/iris/tests/stock/file_headers/README.md
new file mode 100644
index 0000000000..e725b17027
--- /dev/null
+++ b/lib/iris/tests/stock/file_headers/README.md
@@ -0,0 +1,6 @@
+A directory of text files containing file header strings that include
+placeholders - designed to be interpreted by 
+[Python's string.Template](https://docs.python.org/3/library/string.html#string.Template).
+
+* `.cdl` files: used by [`tests.stock.netcdf`](../netcdf.py) to create
+synthetic NetCDF files of the required dimensions.
diff --git a/lib/iris/tests/stock/file_headers/xios_2D_face_half_levels.cdl b/lib/iris/tests/stock/file_headers/xios_2D_face_half_levels.cdl
new file mode 100644
index 0000000000..b135546f2d
--- /dev/null
+++ b/lib/iris/tests/stock/file_headers/xios_2D_face_half_levels.cdl
@@ -0,0 +1,58 @@
+// For now: have omitted all optional information (on edges + connectivity),
+// *except* for face locations.
+netcdf ${DATASET_NAME} {
+dimensions:
+	axis_nbounds = 2 ;
+	Two = 2 ;
+	nMesh2d_half_levels_node = ${NUM_NODES} ;
+	nMesh2d_half_levels_face = ${NUM_FACES} ;
+	nMesh2d_half_levels_vertex = 4 ;
+	time_counter = UNLIMITED ; // (1 currently)
+variables:
+	int Mesh2d_half_levels ;
+		Mesh2d_half_levels:cf_role = "mesh_topology" ;
+		Mesh2d_half_levels:long_name = "Topology data of 2D unstructured mesh" ;
+		Mesh2d_half_levels:topology_dimension = 2 ;
+		Mesh2d_half_levels:node_coordinates = "Mesh2d_half_levels_node_x Mesh2d_half_levels_node_y" ;
+		Mesh2d_half_levels:face_coordinates = "Mesh2d_half_levels_face_x Mesh2d_half_levels_face_y" ;
+		Mesh2d_half_levels:face_node_connectivity = "Mesh2d_half_levels_face_nodes" ;
+	float Mesh2d_half_levels_node_x(nMesh2d_half_levels_node) ;
+		Mesh2d_half_levels_node_x:standard_name = "longitude" ;
+		Mesh2d_half_levels_node_x:long_name = "Longitude of mesh nodes." ;
+		Mesh2d_half_levels_node_x:units = "degrees_east" ;
+	float Mesh2d_half_levels_node_y(nMesh2d_half_levels_node) ;
+		Mesh2d_half_levels_node_y:standard_name = "latitude" ;
+		Mesh2d_half_levels_node_y:long_name = "Latitude of mesh nodes." ;
+		Mesh2d_half_levels_node_y:units = "degrees_north" ;
+	float Mesh2d_half_levels_face_x(nMesh2d_half_levels_face) ;
+		Mesh2d_half_levels_face_x:standard_name = "longitude" ;
+		Mesh2d_half_levels_face_x:long_name = "Characteristic longitude of mesh faces." ;
+		Mesh2d_half_levels_face_x:units = "degrees_east" ;
+	float Mesh2d_half_levels_face_y(nMesh2d_half_levels_face) ;
+		Mesh2d_half_levels_face_y:standard_name = "latitude" ;
+		Mesh2d_half_levels_face_y:long_name = "Characteristic latitude of mesh faces." ;
+		Mesh2d_half_levels_face_y:units = "degrees_north" ;
+	int Mesh2d_half_levels_face_nodes(nMesh2d_half_levels_face, nMesh2d_half_levels_vertex) ;
+		Mesh2d_half_levels_face_nodes:cf_role = "face_node_connectivity" ;
+		Mesh2d_half_levels_face_nodes:long_name = "Maps every face to its corner nodes." ;
+		Mesh2d_half_levels_face_nodes:start_index = 0 ;
+	double time_instant(time_counter) ;
+		time_instant:standard_name = "time" ;
+		time_instant:long_name = "Time axis" ;
+		time_instant:calendar = "gregorian" ;
+		time_instant:units = "seconds since 2016-01-01 15:00:00" ;
+		time_instant:time_origin = "2016-01-01 15:00:00" ;
+		time_instant:bounds = "time_instant_bounds" ;
+	double time_instant_bounds(time_counter, axis_nbounds) ;
+	double thing(time_counter, nMesh2d_half_levels_face) ;
+//      Fictional phenomenon to simplify and avoid un-realistic data/units in the required file.
+		thing:long_name = "thingness" ;
+		thing:mesh = "Mesh2d_half_levels" ;
+		thing:location = "face" ;
+		thing:coordinates = "time_instant Mesh2d_half_levels_face_y Mesh2d_half_levels_face_x" ;
+
+// global attributes:
+		:name = "${DATASET_NAME}" ;
+//		original name = "lfric_ngvat_2D_1t_face_half_levels_main_conv_rain"
+		:Conventions = "UGRID" ;
+}
diff --git a/lib/iris/tests/stock/file_headers/xios_3D_face_full_levels.cdl b/lib/iris/tests/stock/file_headers/xios_3D_face_full_levels.cdl
new file mode 100644
index 0000000000..e4f32de7b7
--- /dev/null
+++ b/lib/iris/tests/stock/file_headers/xios_3D_face_full_levels.cdl
@@ -0,0 +1,61 @@
+// For now: have omitted all optional information (on edges + connectivity),
+// *except* for face locations.
+netcdf ${DATASET_NAME} {
+dimensions:
+	axis_nbounds = 2 ;
+	Two = 2 ;
+	nMesh2d_full_levels_node = ${NUM_NODES} ;
+	nMesh2d_full_levels_face = ${NUM_FACES} ;
+	nMesh2d_full_levels_vertex = 4 ;
+	full_levels = ${NUM_LEVELS} ;
+	time_counter = UNLIMITED ; // (1 currently)
+variables:
+	int Mesh2d_full_levels ;
+		Mesh2d_full_levels:cf_role = "mesh_topology" ;
+		Mesh2d_full_levels:long_name = "Topology data of 2D unstructured mesh" ;
+		Mesh2d_full_levels:topology_dimension = 2 ;
+		Mesh2d_full_levels:node_coordinates = "Mesh2d_full_levels_node_x Mesh2d_full_levels_node_y" ;
+		Mesh2d_full_levels:face_coordinates = "Mesh2d_full_levels_face_x Mesh2d_full_levels_face_y" ;
+		Mesh2d_full_levels:face_node_connectivity = "Mesh2d_full_levels_face_nodes" ;
+	float Mesh2d_full_levels_node_x(nMesh2d_full_levels_node) ;
+		Mesh2d_full_levels_node_x:standard_name = "longitude" ;
+		Mesh2d_full_levels_node_x:long_name = "Longitude of mesh nodes." ;
+		Mesh2d_full_levels_node_x:units = "degrees_east" ;
+	float Mesh2d_full_levels_node_y(nMesh2d_full_levels_node) ;
+		Mesh2d_full_levels_node_y:standard_name = "latitude" ;
+		Mesh2d_full_levels_node_y:long_name = "Latitude of mesh nodes." ;
+		Mesh2d_full_levels_node_y:units = "degrees_north" ;
+	float Mesh2d_full_levels_face_x(nMesh2d_full_levels_face) ;
+		Mesh2d_full_levels_face_x:standard_name = "longitude" ;
+		Mesh2d_full_levels_face_x:long_name = "Characteristic longitude of mesh faces." ;
+		Mesh2d_full_levels_face_x:units = "degrees_east" ;
+	float Mesh2d_full_levels_face_y(nMesh2d_full_levels_face) ;
+		Mesh2d_full_levels_face_y:standard_name = "latitude" ;
+		Mesh2d_full_levels_face_y:long_name = "Characteristic latitude of mesh faces." ;
+		Mesh2d_full_levels_face_y:units = "degrees_north" ;
+	int Mesh2d_full_levels_face_nodes(nMesh2d_full_levels_face, nMesh2d_full_levels_vertex) ;
+		Mesh2d_full_levels_face_nodes:cf_role = "face_node_connectivity" ;
+		Mesh2d_full_levels_face_nodes:long_name = "Maps every face to its corner nodes." ;
+		Mesh2d_full_levels_face_nodes:start_index = 0 ;
+	float full_levels(full_levels) ;
+		full_levels:name = "full_levels" ;
+	double time_instant(time_counter) ;
+		time_instant:standard_name = "time" ;
+		time_instant:long_name = "Time axis" ;
+		time_instant:calendar = "gregorian" ;
+		time_instant:units = "seconds since 2016-01-01 15:00:00" ;
+		time_instant:time_origin = "2016-01-01 15:00:00" ;
+		time_instant:bounds = "time_instant_bounds" ;
+	double time_instant_bounds(time_counter, axis_nbounds) ;
+	double thing(time_counter, full_levels, nMesh2d_full_levels_face) ;
+//      Fictional phenomenon to simplify and avoid un-realistic data/units in the required file.
+		thing:long_name = "thingness" ;
+		thing:mesh = "Mesh2d_full_levels" ;
+		thing:location = "face" ;
+		thing:coordinates = "time_instant Mesh2d_full_levels_face_y Mesh2d_full_levels_face_x" ;
+
+// global attributes:
+		:name = "${DATASET_NAME}" ;
+//		original name = "lfric_ngvat_3D_1t_full_level_face_grid_main_u3"
+		:Conventions = "UGRID" ;
+}
diff --git a/lib/iris/tests/stock/file_headers/xios_3D_face_half_levels.cdl b/lib/iris/tests/stock/file_headers/xios_3D_face_half_levels.cdl
new file mode 100644
index 0000000000..a193dbe451
--- /dev/null
+++ b/lib/iris/tests/stock/file_headers/xios_3D_face_half_levels.cdl
@@ -0,0 +1,61 @@
+// For now: have omitted all optional information (on edges + connectivity),
+// *except* for face locations.
+netcdf ${DATASET_NAME} {
+dimensions:
+	axis_nbounds = 2 ;
+	Two = 2 ;
+	nMesh2d_half_levels_node = ${NUM_NODES} ;
+	nMesh2d_half_levels_face = ${NUM_FACES} ;
+	nMesh2d_half_levels_vertex = 4 ;
+	half_levels = ${NUM_LEVELS} ;
+	time_counter = UNLIMITED ; // (1 currently)
+variables:
+	int Mesh2d_half_levels ;
+		Mesh2d_half_levels:cf_role = "mesh_topology" ;
+		Mesh2d_half_levels:long_name = "Topology data of 2D unstructured mesh" ;
+		Mesh2d_half_levels:topology_dimension = 2 ;
+		Mesh2d_half_levels:node_coordinates = "Mesh2d_half_levels_node_x Mesh2d_half_levels_node_y" ;
+		Mesh2d_half_levels:face_coordinates = "Mesh2d_half_levels_face_x Mesh2d_half_levels_face_y" ;
+		Mesh2d_half_levels:face_node_connectivity = "Mesh2d_half_levels_face_nodes" ;
+	float Mesh2d_half_levels_node_x(nMesh2d_half_levels_node) ;
+		Mesh2d_half_levels_node_x:standard_name = "longitude" ;
+		Mesh2d_half_levels_node_x:long_name = "Longitude of mesh nodes." ;
+		Mesh2d_half_levels_node_x:units = "degrees_east" ;
+	float Mesh2d_half_levels_node_y(nMesh2d_half_levels_node) ;
+		Mesh2d_half_levels_node_y:standard_name = "latitude" ;
+		Mesh2d_half_levels_node_y:long_name = "Latitude of mesh nodes." ;
+		Mesh2d_half_levels_node_y:units = "degrees_north" ;
+	float Mesh2d_half_levels_face_x(nMesh2d_half_levels_face) ;
+		Mesh2d_half_levels_face_x:standard_name = "longitude" ;
+		Mesh2d_half_levels_face_x:long_name = "Characteristic longitude of mesh faces." ;
+		Mesh2d_half_levels_face_x:units = "degrees_east" ;
+	float Mesh2d_half_levels_face_y(nMesh2d_half_levels_face) ;
+		Mesh2d_half_levels_face_y:standard_name = "latitude" ;
+		Mesh2d_half_levels_face_y:long_name = "Characteristic latitude of mesh faces." ;
+		Mesh2d_half_levels_face_y:units = "degrees_north" ;
+	int Mesh2d_half_levels_face_nodes(nMesh2d_half_levels_face, nMesh2d_half_levels_vertex) ;
+		Mesh2d_half_levels_face_nodes:cf_role = "face_node_connectivity" ;
+		Mesh2d_half_levels_face_nodes:long_name = "Maps every face to its corner nodes." ;
+		Mesh2d_half_levels_face_nodes:start_index = 0 ;
+	float half_levels(half_levels) ;
+		half_levels:name = "half_levels" ;
+	double time_instant(time_counter) ;
+		time_instant:standard_name = "time" ;
+		time_instant:long_name = "Time axis" ;
+		time_instant:calendar = "gregorian" ;
+		time_instant:units = "seconds since 2016-01-01 15:00:00" ;
+		time_instant:time_origin = "2016-01-01 15:00:00" ;
+		time_instant:bounds = "time_instant_bounds" ;
+	double time_instant_bounds(time_counter, axis_nbounds) ;
+	double thing(time_counter, half_levels, nMesh2d_half_levels_face) ;
+//      Fictional phenomenon to simplify and avoid un-realistic data/units in the required file.
+		thing:long_name = "thingness" ;
+		thing:mesh = "Mesh2d_half_levels" ;
+		thing:location = "face" ;
+		thing:coordinates = "time_instant Mesh2d_half_levels_face_y Mesh2d_half_levels_face_x" ;
+
+// global attributes:
+		:name = "${DATASET_NAME}" ;
+//		original name = "lfric_ngvat_3D_1t_half_level_face_grid_derived_theta_in_w3"
+		:Conventions = "UGRID" ;
+}
diff --git a/lib/iris/tests/stock/netcdf.py b/lib/iris/tests/stock/netcdf.py
new file mode 100644
index 0000000000..78dc08eafd
--- /dev/null
+++ b/lib/iris/tests/stock/netcdf.py
@@ -0,0 +1,215 @@
+# Copyright Iris contributors
+#
+# This file is part of Iris and is released under the LGPL license.
+# See COPYING and COPYING.LESSER in the root of the repository for full
+# licensing details.
+"""Routines for generating synthetic NetCDF files from template headers."""
+
+from pathlib import Path
+from string import Template
+import subprocess
+
+import netCDF4
+import numpy as np
+
+
+def _file_from_cdl_template(
+    temp_file_dir, dataset_name, dataset_type, template_subs
+):
+    """Shared template filling behaviour.
+
+    Substitutes placeholders in the appropriate CDL template, saves to a
+    NetCDF file.
+
+    """
+    nc_write_path = (
+        Path(temp_file_dir).joinpath(dataset_name).with_suffix(".nc")
+    )
+    # Fetch the specified CDL template type.
+    templates_dir = Path(__file__).parent / "file_headers"
+    template_filepath = templates_dir.joinpath(dataset_type).with_suffix(
+        ".cdl"
+    )
+    # Substitute placeholders.
+    with open(template_filepath) as file:
+        template_string = Template(file.read())
+    cdl = template_string.substitute(template_subs)
+
+    # Spawn an "ncgen" command to create an actual NetCDF file from the
+    # CDL string.
+    subprocess.run(
+        ["ncgen", "-o" + str(nc_write_path)],
+        input=cdl,
+        encoding="ascii",
+        check=True,
+    )
+
+    return nc_write_path
+
+
+def _add_standard_data(nc_path, unlimited_dim_size=0):
+    """Shared data populating behaviour.
+
+    Adds placeholder data to the variables in a NetCDF file, accounting for
+    dimension size, 'dimension coordinates' and a possible unlimited dimension.
+
+    """
+
+    ds = netCDF4.Dataset(nc_path, "r+")
+
+    unlimited_dim_names = [
+        dim for dim in ds.dimensions if ds.dimensions[dim].size == 0
+    ]
+    # Data addition dependent on this assumption:
+    assert len(unlimited_dim_names) < 2
+
+    # Fill variables data with placeholder numbers.
+    for var in ds.variables.values():
+        shape = list(var.shape)
+        dims = var.dimensions
+        # Fill the unlimited dimension with the input size.
+        shape = [
+            unlimited_dim_size if dim == unlimited_dim_names[0] else size
+            for dim, size in zip(dims, shape)
+        ]
+        data = np.zeros(shape, dtype=var.dtype)
+        if len(var.dimensions) == 1 and var.dimensions[0] == var.name:
+            # Fill the var with ascending values (not all zeroes),
+            # so it can be a dim-coord.
+            data = np.arange(data.size, dtype=data.dtype).reshape(data.shape)
+        var[:] = data
+
+    ds.close()
+
+
+def create_file__xios_2d_face_half_levels(
+    temp_file_dir, dataset_name, n_faces=866, n_times=1
+):
+    """Create a synthetic NetCDF file with XIOS-like content.
+
+    Starts from a template CDL headers string, modifies to the input
+    dimensions then adds data of the correct size.
+
+    Parameters
+    ----------
+    temp_file_dir : str or pathlib.Path
+        The directory in which to place the created file.
+    dataset_name : str
+        The name for the NetCDF dataset and also the created file.
+    n_faces, n_times: int
+        Dimension sizes for the dataset.
+
+    Returns
+    -------
+    str
+        Path of the created NetCDF file.
+
+    """
+    dataset_type = "xios_2D_face_half_levels"
+
+    # Set the placeholder substitutions.
+    template_subs = {
+        "DATASET_NAME": dataset_name,
+        "NUM_NODES": n_faces + 2,
+        "NUM_FACES": n_faces,
+    }
+
+    # Create a NetCDF file based on the dataset type template and substitutions.
+    nc_path = _file_from_cdl_template(
+        temp_file_dir, dataset_name, dataset_type, template_subs
+    )
+
+    # Populate with the standard set of data, sized correctly.
+    _add_standard_data(nc_path, unlimited_dim_size=n_times)
+
+    return str(nc_path)
+
+
+def create_file__xios_3d_face_half_levels(
+    temp_file_dir, dataset_name, n_faces=866, n_times=1, n_levels=38
+):
+    """Create a synthetic NetCDF file with XIOS-like content.
+
+    Starts from a template CDL headers string, modifies to the input
+    dimensions then adds data of the correct size.
+
+    Parameters
+    ----------
+    temp_file_dir : str or pathlib.Path
+        The directory in which to place the created file.
+    dataset_name : str
+        The name for the NetCDF dataset and also the created file.
+    n_faces, n_times, n_levels: int
+        Dimension sizes for the dataset.
+
+    Returns
+    -------
+    str
+        Path of the created NetCDF file.
+
+    """
+    dataset_type = "xios_3D_face_half_levels"
+
+    # Set the placeholder substitutions.
+    template_subs = {
+        "DATASET_NAME": dataset_name,
+        "NUM_NODES": n_faces + 2,
+        "NUM_FACES": n_faces,
+        "NUM_LEVELS": n_levels,
+    }
+
+    # Create a NetCDF file based on the dataset type template and
+    # substitutions.
+    nc_path = _file_from_cdl_template(
+        temp_file_dir, dataset_name, dataset_type, template_subs
+    )
+
+    # Populate with the standard set of data, sized correctly.
+    _add_standard_data(nc_path, unlimited_dim_size=n_times)
+
+    return str(nc_path)
+
+
+def create_file__xios_3d_face_full_levels(
+    temp_file_dir, dataset_name, n_faces=866, n_times=1, n_levels=39
+):
+    """Create a synthetic NetCDF file with XIOS-like content.
+
+    Starts from a template CDL headers string, modifies to the input
+    dimensions then adds data of the correct size.
+
+    Parameters
+    ----------
+    temp_file_dir : str or pathlib.Path
+        The directory in which to place the created file.
+    dataset_name : str
+        The name for the NetCDF dataset and also the created file.
+    n_faces, n_times, n_levels: int
+        Dimension sizes for the dataset.
+
+    Returns
+    -------
+    str
+        Path of the created NetCDF file.
+
+    """
+    dataset_type = "xios_3D_face_full_levels"
+
+    # Set the placeholder substitutions.
+    template_subs = {
+        "DATASET_NAME": dataset_name,
+        "NUM_NODES": n_faces + 2,
+        "NUM_FACES": n_faces,
+        "NUM_LEVELS": n_levels,
+    }
+
+    # Create a NetCDF file based on the dataset type template and
+    # substitutions.
+    nc_path = _file_from_cdl_template(
+        temp_file_dir, dataset_name, dataset_type, template_subs
+    )
+
+    # Populate with the standard set of data, sized correctly.
+    _add_standard_data(nc_path, unlimited_dim_size=n_times)
+
+    return str(nc_path)
diff --git a/lib/iris/tests/unit/tests/stock/__init__.py b/lib/iris/tests/unit/tests/stock/__init__.py
new file mode 100644
index 0000000000..f91390c2b3
--- /dev/null
+++ b/lib/iris/tests/unit/tests/stock/__init__.py
@@ -0,0 +1,6 @@
+# Copyright Iris contributors
+#
+# This file is part of Iris and is released under the LGPL license.
+# See COPYING and COPYING.LESSER in the root of the repository for full
+# licensing details.
+"""Unit tests for the :mod:`iris.tests.stock` module."""
diff --git a/lib/iris/tests/unit/tests/stock/test_netcdf.py b/lib/iris/tests/unit/tests/stock/test_netcdf.py
new file mode 100644
index 0000000000..afb94e695f
--- /dev/null
+++ b/lib/iris/tests/unit/tests/stock/test_netcdf.py
@@ -0,0 +1,141 @@
+# Copyright Iris contributors
+#
+# This file is part of Iris and is released under the LGPL license.
+# See COPYING and COPYING.LESSER in the root of the repository for full
+# licensing details.
+"""Unit tests for the :mod:`iris.tests.stock.netcdf` module."""
+
+# Import iris.tests first so that some things can be initialised before
+# importing anything else.
+import iris.tests as tests
+
+import shutil
+import tempfile
+
+from iris import load_cube
+from iris.tests.stock import netcdf
+from iris.experimental.ugrid import Mesh, MeshCoord, PARSE_UGRID_ON_LOAD
+
+
+class XIOSFileMixin(tests.IrisTest):
+    @classmethod
+    def setUpClass(cls):
+        # Create a temp directory for transient test files.
+        cls.temp_dir = tempfile.mkdtemp()
+
+    @classmethod
+    def tearDownClass(cls):
+        # Destroy the temp directory.
+        shutil.rmtree(cls.temp_dir)
+
+    def create_synthetic_file(self, **create_kwargs):
+        # Should be overridden to invoke one of the create_file_ functions.
+        # E.g.
+        # return netcdf.create_file__xios_2d_face_half_levels(
+        #     temp_file_dir=self.temp_dir, dataset_name="mesh", **create_kwargs
+        # )
+        raise NotImplementedError
+
+    def create_synthetic_test_cube(self, **create_kwargs):
+        file_path = self.create_synthetic_file(**create_kwargs)
+        with PARSE_UGRID_ON_LOAD.context():
+            cube = load_cube(file_path)
+        return cube
+
+    def check_cube(self, cube, shape, location, level):
+        # Basic checks on the primary data cube.
+        self.assertEqual(cube.var_name, "thing")
+        self.assertEqual(cube.long_name, "thingness")
+        self.assertEqual(cube.shape, shape)
+
+        # Also a few checks on the attached mesh-related information.
+        last_dim = cube.ndim - 1
+        self.assertIsInstance(cube.mesh, Mesh)
+        self.assertEqual(cube.mesh_dim(), last_dim)
+        self.assertEqual(cube.location, location)
+        for coord_name in ("longitude", "latitude"):
+            coord = cube.coord(coord_name)
+            self.assertIsInstance(coord, MeshCoord)
+            self.assertEqual(coord.shape, (shape[last_dim],))
+        self.assertTrue(cube.mesh.var_name.endswith(f"{level}_levels"))
+
+
+class Test_create_file__xios_2d_face_half_levels(XIOSFileMixin):
+    def create_synthetic_file(self, **create_kwargs):
+        return netcdf.create_file__xios_2d_face_half_levels(
+            temp_file_dir=self.temp_dir, dataset_name="mesh", **create_kwargs
+        )
+
+    def test_basic_load(self):
+        cube = self.create_synthetic_test_cube()
+        self.check_cube(cube, shape=(1, 866), location="face", level="half")
+
+    def test_scale_mesh(self):
+        cube = self.create_synthetic_test_cube(n_faces=10)
+        self.check_cube(cube, shape=(1, 10), location="face", level="half")
+
+    def test_scale_time(self):
+        cube = self.create_synthetic_test_cube(n_times=3)
+        self.check_cube(cube, shape=(3, 866), location="face", level="half")
+
+
+class Test_create_file__xios_3d_face_half_levels(XIOSFileMixin):
+    def create_synthetic_file(self, **create_kwargs):
+        return netcdf.create_file__xios_3d_face_half_levels(
+            temp_file_dir=self.temp_dir, dataset_name="mesh", **create_kwargs
+        )
+
+    def test_basic_load(self):
+        cube = self.create_synthetic_test_cube()
+        self.check_cube(
+            cube, shape=(1, 38, 866), location="face", level="half"
+        )
+
+    def test_scale_mesh(self):
+        cube = self.create_synthetic_test_cube(n_faces=10)
+        self.check_cube(cube, shape=(1, 38, 10), location="face", level="half")
+
+    def test_scale_time(self):
+        cube = self.create_synthetic_test_cube(n_times=3)
+        self.check_cube(
+            cube, shape=(3, 38, 866), location="face", level="half"
+        )
+
+    def test_scale_levels(self):
+        cube = self.create_synthetic_test_cube(n_levels=10)
+        self.check_cube(
+            cube, shape=(1, 10, 866), location="face", level="half"
+        )
+
+
+class Test_create_file__xios_3d_face_full_levels(XIOSFileMixin):
+    def create_synthetic_file(self, **create_kwargs):
+        return netcdf.create_file__xios_3d_face_full_levels(
+            temp_file_dir=self.temp_dir, dataset_name="mesh", **create_kwargs
+        )
+
+    def test_basic_load(self):
+        cube = self.create_synthetic_test_cube()
+        self.check_cube(
+            cube, shape=(1, 39, 866), location="face", level="full"
+        )
+
+    def test_scale_mesh(self):
+        cube = self.create_synthetic_test_cube(n_faces=10)
+        self.check_cube(cube, shape=(1, 39, 10), location="face", level="full")
+
+    def test_scale_time(self):
+        cube = self.create_synthetic_test_cube(n_times=3)
+        self.check_cube(
+            cube, shape=(3, 39, 866), location="face", level="full"
+        )
+
+    def test_scale_levels(self):
+        cube = self.create_synthetic_test_cube(n_levels=10)
+        self.check_cube(
+            cube, shape=(1, 10, 866), location="face", level="full"
+        )
+
+
+if __name__ == "__main__":
+    tests.main()