diff --git a/.github/workflows/test_suite.yml b/.github/workflows/test_suite.yml
index 7d51c64a9..40d5a6520 100644
--- a/.github/workflows/test_suite.yml
+++ b/.github/workflows/test_suite.yml
@@ -86,7 +86,6 @@ jobs:
         ./run_test_jan2010_integration.sh python3
         ./run_test_advection.sh python3
         ./run_test_error_handling_test.sh python3
-        cd -
 
   # MPI tests on latest Ubuntu OS without XIOS support
   test-ubuntu-mpi-noxios:
@@ -246,4 +245,3 @@ jobs:
         ./run_test_jan2010_integration.sh python
         ./run_test_advection.sh python
         ./run_test_error_handling_test.sh python
-        cd -
diff --git a/core/src/ParaGridIO.cpp b/core/src/ParaGridIO.cpp
index 8db1204fb..7ea559ea5 100644
--- a/core/src/ParaGridIO.cpp
+++ b/core/src/ParaGridIO.cpp
@@ -33,14 +33,15 @@ ParaGridIO::ParaGridIO(ParametricGrid& grid)
     , dimensionKeys({
           // clang-format off
           // Accept post-May 2024 (xdim, ydim, zdim) dimension names and pre-May 2024 (x, y, z)
-        { "yx", ModelArray::Type::H },
         { "ydimxdim", ModelArray::Type::H },
-        { "yxdg_comp", ModelArray::Type::DG },
+        { "y_dimx_dim", ModelArray::Type::H },
         { "ydimxdimdg_comp", ModelArray::Type::DG },
-        { "yxdgstress_comp", ModelArray::Type::DGSTRESS },
+        { "y_dimx_dimdg_comp", ModelArray::Type::DG },
         { "ydimxdimdgstress_comp", ModelArray::Type::DGSTRESS },
-        { "ycgxcg", ModelArray::Type::CG },
+        { "y_dimx_dimdgstress_comp", ModelArray::Type::DGSTRESS },
+        { "y_cgx_cg", ModelArray::Type::CG },
         { "yvertexxvertexncoords", ModelArray::Type::VERTEX },
+        { "y_vertexx_vertexncoords", ModelArray::Type::VERTEX },
           // clang-format on
       })
     , isDG({
diff --git a/core/src/Xios.cpp b/core/src/Xios.cpp
index 6332e654d..63fa0106c 100644
--- a/core/src/Xios.cpp
+++ b/core/src/Xios.cpp
@@ -501,83 +501,112 @@ void Xios::setupDomains()
         // Set domain extents based on model metadata
         size_t counter = 0;
         for (ModelArray::Dimension& dim : ModelArray::typeDimensions[type]) {
-            const std::string& domainName = ModelArray::definedDimensions[dim].name;
             if (counter == 0) {
+                const std::string dimName = "x";
+                int ni_glo;
+                int ni;
+                int ibegin;
                 if (dim == ModelArray::Dimension::X) {
-                    cxios_set_domain_ni_glo(domain, metadata.getGlobalExtentX());
-                    cxios_set_domain_ni(domain, metadata.getLocalExtentX());
-                    cxios_set_domain_ibegin(domain, metadata.getLocalCornerX());
+                    ni_glo = metadata.getGlobalExtentX();
+                    ni = metadata.getLocalExtentX();
+                    ibegin = metadata.getLocalCornerX();
                 } else if (dim == ModelArray::Dimension::XVERTEX) {
-                    cxios_set_domain_ni_glo(domain, metadata.getGlobalExtentX() + 1);
-                    cxios_set_domain_ni(domain, metadata.getLocalExtentX() + 1);
-                    cxios_set_domain_ibegin(domain, metadata.getLocalCornerX());
+                    ni_glo = metadata.getGlobalExtentX() + 1;
+                    ni = metadata.getLocalExtentX() + 1;
+                    ibegin = metadata.getLocalCornerX();
                 } else if (dim == ModelArray::Dimension::XCG) {
-                    cxios_set_domain_ni_glo(domain, CGDEGREE * metadata.getGlobalExtentX() + 1);
-                    cxios_set_domain_ni(domain, CGDEGREE * metadata.getLocalExtentX() + 1);
-                    cxios_set_domain_ibegin(domain, CGDEGREE * metadata.getLocalCornerX());
+                    ni_glo = CGDEGREE * metadata.getGlobalExtentX() + 1;
+                    ni = CGDEGREE * metadata.getLocalExtentX() + 1;
+                    ibegin = CGDEGREE * metadata.getLocalCornerX();
                 } else {
                     throw std::runtime_error(
-                        "Xios: Could not set domain extents based on dimension '"
-                        + ModelArray::definedDimensions.at(dim).name + "'");
+                        "Xios: Could not set domain extents based on dimension '" + dimName + "'");
                 }
+                cxios_set_domain_ni_glo(domain, ni_glo);
                 if (!cxios_is_defined_domain_ni_glo(domain)) {
                     throw std::runtime_error(
                         "Xios: Failed to set global x-size for domain '" + domainId + "'");
                 }
+                cxios_set_domain_ni(domain, ni);
                 if (!cxios_is_defined_domain_ni(domain)) {
                     throw std::runtime_error(
                         "Xios: Failed to set local x-size for domain '" + domainId + "'");
                 }
+                cxios_set_domain_ibegin(domain, ibegin);
                 if (!cxios_is_defined_domain_ibegin(domain)) {
                     throw std::runtime_error(
                         "Xios: Failed to set local starting x-index for domain '" + domainId + "'");
                 }
-                cxios_set_domain_dim_i_name(domain, domainName.c_str(), domainName.length());
+                std::vector<double> lonvalue;
+                for (int i = 0; i < ni; i++) {
+                    lonvalue.push_back(ibegin + i);
+                }
+                cxios_set_domain_lonvalue_1d(domain, lonvalue.data(), &ni);
+                if (!cxios_is_defined_domain_lonvalue_1d(domain)) {
+                    throw std::runtime_error(
+                        "Xios: Failed to set local x-indices for domain '" + domainId + "'");
+                }
+                cxios_set_domain_dim_i_name(domain, dimName.c_str(), dimName.length());
                 if (!cxios_is_defined_domain_dim_i_name(domain)) {
                     throw std::runtime_error(
                         "Xios: Failed to set x-coordinate name for domain '" + domainId + "'");
                 }
-                cxios_set_domain_lon_name(domain, domainName.c_str(), domainName.length());
+                cxios_set_domain_lon_name(domain, dimName.c_str(), dimName.length());
                 if (!cxios_is_defined_domain_lon_name(domain)) {
                     throw std::runtime_error(
                         "Xios: Failed to set longitude name for domain '" + domainId + "'");
                 }
             } else if (counter == 1) {
+                const std::string dimName = "y";
+                int nj_glo;
+                int nj;
+                int jbegin;
                 if (dim == ModelArray::Dimension::Y) {
-                    cxios_set_domain_nj_glo(domain, metadata.getGlobalExtentY());
-                    cxios_set_domain_nj(domain, metadata.getLocalExtentY());
-                    cxios_set_domain_jbegin(domain, metadata.getLocalCornerY());
+                    nj_glo = metadata.getGlobalExtentY();
+                    nj = metadata.getLocalExtentY();
+                    jbegin = metadata.getLocalCornerY();
                 } else if (dim == ModelArray::Dimension::YVERTEX) {
-                    cxios_set_domain_nj_glo(domain, metadata.getGlobalExtentY() + 1);
-                    cxios_set_domain_nj(domain, metadata.getLocalExtentY() + 1);
-                    cxios_set_domain_jbegin(domain, metadata.getLocalCornerY());
+                    nj_glo = metadata.getGlobalExtentY() + 1;
+                    nj = metadata.getLocalExtentY() + 1;
+                    jbegin = metadata.getLocalCornerY();
                 } else if (dim == ModelArray::Dimension::YCG) {
-                    cxios_set_domain_nj_glo(domain, CGDEGREE * metadata.getGlobalExtentY() + 1);
-                    cxios_set_domain_nj(domain, CGDEGREE * metadata.getLocalExtentY() + 1);
-                    cxios_set_domain_jbegin(domain, CGDEGREE * metadata.getLocalCornerY());
+                    nj_glo = CGDEGREE * metadata.getGlobalExtentY() + 1;
+                    nj = CGDEGREE * metadata.getLocalExtentY() + 1;
+                    jbegin = CGDEGREE * metadata.getLocalCornerY();
                 } else {
                     throw std::runtime_error(
-                        "Xios: Could not set domain extents based on dimension '"
-                        + ModelArray::definedDimensions.at(dim).name + "'");
+                        "Xios: Could not set domain extents based on dimension '" + dimName + "'");
                 }
+                cxios_set_domain_nj_glo(domain, nj_glo);
                 if (!cxios_is_defined_domain_nj_glo(domain)) {
                     throw std::runtime_error(
                         "Xios: Failed to set global y-size for domain '" + domainId + "'");
                 }
+                cxios_set_domain_nj(domain, nj);
                 if (!cxios_is_defined_domain_nj(domain)) {
                     throw std::runtime_error(
                         "Xios: Failed to set local y-size for domain '" + domainId + "'");
                 }
+                cxios_set_domain_jbegin(domain, jbegin);
                 if (!cxios_is_defined_domain_jbegin(domain)) {
                     throw std::runtime_error(
                         "Xios: Failed to set local starting y-index for domain '" + domainId + "'");
                 }
-                cxios_set_domain_dim_j_name(domain, domainName.c_str(), domainName.length());
+                std::vector<double> latvalue;
+                for (int j = 0; j < nj; j++) {
+                    latvalue.push_back(jbegin + j);
+                }
+                cxios_set_domain_latvalue_1d(domain, latvalue.data(), &nj);
+                if (!cxios_is_defined_domain_latvalue_1d(domain)) {
+                    throw std::runtime_error(
+                        "Xios: Failed to set local y-indices for domain '" + domainId + "'");
+                }
+                cxios_set_domain_dim_j_name(domain, dimName.c_str(), dimName.length());
                 if (!cxios_is_defined_domain_dim_j_name(domain)) {
                     throw std::runtime_error(
                         "Xios: Failed to set y-coordinate name for domain '" + domainId + "'");
                 }
-                cxios_set_domain_lat_name(domain, domainName.c_str(), domainName.length());
+                cxios_set_domain_lat_name(domain, dimName.c_str(), dimName.length());
                 if (!cxios_is_defined_domain_lat_name(domain)) {
                     throw std::runtime_error(
                         "Xios: Failed to set latitude name for domain '" + domainId + "'");
@@ -928,12 +957,12 @@ void Xios::setupFields()
 
         // Create map for field types
         const std::map<std::string, ModelArray::Type> dimensionKeys = {
-            { "yx", ModelArray::Type::H },
             { "ydimxdim", ModelArray::Type::H },
-            { "yxdg_comp", ModelArray::Type::DG },
+            { "y_dimx_dim", ModelArray::Type::H },
             { "ydimxdimdg_comp", ModelArray::Type::DG },
-            { "yxdgstress_comp", ModelArray::Type::DGSTRESS },
+            { "y_dimx_dimdg_comp", ModelArray::Type::DG },
             { "ydimxdimdgstress_comp", ModelArray::Type::DGSTRESS },
+            { "y_dimx_dimdgstress_comp", ModelArray::Type::DGSTRESS },
             { "y_cgx_cg", ModelArray::Type::CG },
             { "yvertexxvertexncoords", ModelArray::Type::VERTEX },
             { "y_vertexx_vertexncoords", ModelArray::Type::VERTEX },
diff --git a/core/src/discontinuousgalerkin/ModelArrayDetails.cpp b/core/src/discontinuousgalerkin/ModelArrayDetails.cpp
index 8e6132983..4897a499d 100644
--- a/core/src/discontinuousgalerkin/ModelArrayDetails.cpp
+++ b/core/src/discontinuousgalerkin/ModelArrayDetails.cpp
@@ -27,10 +27,10 @@ namespace Nextsim {
 // clang-format off
 std::map<ModelArray::Dimension, ModelArray::DimensionSpec> ModelArray::definedDimensions = {
     // set default global size, local size and start position for each dimension
-    { ModelArray::Dimension::X, { "xdim", "x", 0, 0, 0 } },
-    { ModelArray::Dimension::Y, { "ydim", "y", 0, 0, 0 } },
-    { ModelArray::Dimension::XVERTEX, { "xvertex", "xvertex", 1, 1, 0 } }, // defined as x + 1
-    { ModelArray::Dimension::YVERTEX, { "yvertex", "yvertex", 1, 1, 0 } }, // defined as y + 1
+    { ModelArray::Dimension::X, { "x_dim", "xdim", 0, 0, 0 } },
+    { ModelArray::Dimension::Y, { "y_dim", "ydim", 0, 0, 0 } },
+    { ModelArray::Dimension::XVERTEX, { "x_vertex", "xvertex", 1, 1, 0 } }, // defined as x + 1
+    { ModelArray::Dimension::YVERTEX, { "y_vertex", "yvertex", 1, 1, 0 } }, // defined as y + 1
     { ModelArray::Dimension::XCG, { "x_cg", "x_cg", CGDEGREE, CGDEGREE, 0 } },
     { ModelArray::Dimension::YCG, { "y_cg", "y_cg", CGDEGREE, CGDEGREE, 0 } },
     // The DG components are also included here to store the names
diff --git a/core/src/finitevolume/ModelArrayDetails.cpp b/core/src/finitevolume/ModelArrayDetails.cpp
index 4c5a83338..e80745bff 100644
--- a/core/src/finitevolume/ModelArrayDetails.cpp
+++ b/core/src/finitevolume/ModelArrayDetails.cpp
@@ -13,10 +13,10 @@
 
 namespace Nextsim {
 std::map<ModelArray::Dimension, ModelArray::DimensionSpec> ModelArray::definedDimensions = {
-    { ModelArray::Dimension::X, { "xdim", "x", 0, 0, 0 } },
-    { ModelArray::Dimension::Y, { "ydim", "y", 0, 0, 0 } },
-    { ModelArray::Dimension::XVERTEX, { "xvertex", "xvertex", 1, 1, 0 } }, // defined as x + 1
-    { ModelArray::Dimension::YVERTEX, { "yvertex", "yvertex", 1, 1, 0 } }, // defined as y + 1
+    { ModelArray::Dimension::X, { "x_dim", "xdim", 0, 0, 0 } },
+    { ModelArray::Dimension::Y, { "y_dim", "ydim", 0, 0, 0 } },
+    { ModelArray::Dimension::XVERTEX, { "x_vertex", "xvertex", 1, 1, 0 } }, // defined as x + 1
+    { ModelArray::Dimension::YVERTEX, { "y_vertex", "yvertex", 1, 1, 0 } }, // defined as y + 1
 };
 
 ModelArray::TypeDimensions ModelArray::typeDimensions = {
diff --git a/core/src/include/Xios.hpp b/core/src/include/Xios.hpp
index f6b1861fe..71cef9f4d 100644
--- a/core/src/include/Xios.hpp
+++ b/core/src/include/Xios.hpp
@@ -151,12 +151,17 @@ class Xios : public Configured<Xios> {
     xios::CAxis* getAxis(const std::string& axisId);
 
     /* Domain */
+    // NOTE: Dimension names get processed as <dim>_<domainId> by XIOS, so we define the domainIds
+    //       so that these coincide with the altnames when applied to dimensions x and y.
     std::map<ModelArray::Type, std::string> domainIds = {
-        { ModelArray::Type::H, "HDomain" },
-        { ModelArray::Type::VERTEX, "VertexDomain" },
-        { ModelArray::Type::DG, "HDomain" },
-        { ModelArray::Type::DGSTRESS, "HDomain" },
-        { ModelArray::Type::CG, "CGDomain" },
+        // Standard cell-based x- and y-dimensions (alt. names x_dim and y_dim)
+        { ModelArray::Type::H, "dim" },
+        { ModelArray::Type::DG, "dim" },
+        { ModelArray::Type::DGSTRESS, "dim" },
+        // Vertex-based x- and y-dimensions (alt. names x_vertex and y_vertex)
+        { ModelArray::Type::VERTEX, "vertex" },
+        // CG-based x- and y-dimensions (alt. names x_cg and y_cg)
+        { ModelArray::Type::CG, "cg" },
     };
     xios::CDomainGroup* getDomainGroup();
     xios::CDomain* getDomain(const std::string& domainId);
diff --git a/core/src/include/gridNames.hpp b/core/src/include/gridNames.hpp
index e3313f89a..ca95c2b9f 100644
--- a/core/src/include/gridNames.hpp
+++ b/core/src/include/gridNames.hpp
@@ -40,8 +40,8 @@ static const std::string vIOStressName = "viostress";
 static const std::string coordsName = "coords";
 static const std::string latitudeName = "latitude";
 static const std::string longitudeName = "longitude";
-static const std::string xName = "x";
-static const std::string yName = "y";
+static const std::string xName = "x_dim";
+static const std::string yName = "y_dim";
 static const std::string gridAzimuthName = "grid_azimuth";
 
 static const std::string mdiName = "missing_value";
diff --git a/core/src/include/xios_c_interface.hpp b/core/src/include/xios_c_interface.hpp
index 845e3c294..e9a84a5c9 100644
--- a/core/src/include/xios_c_interface.hpp
+++ b/core/src/include/xios_c_interface.hpp
@@ -90,10 +90,14 @@ bool cxios_is_defined_domain_lat_name(xios::CDomain* domain_hdl);
 bool cxios_is_defined_domain_lon_name(xios::CDomain* domain_hdl);
 bool cxios_is_defined_domain_ni_glo(xios::CDomain* domain_hdl);
 bool cxios_is_defined_domain_nj_glo(xios::CDomain* domain_hdl);
+void cxios_set_domain_lonvalue_1d(xios::CDomain* domain_hdl, double* data, const int* ni);
+void cxios_set_domain_latvalue_1d(xios::CDomain* domain_hdl, double* data, const int* nj);
 bool cxios_is_defined_domain_ni(xios::CDomain* domain_hdl);
 bool cxios_is_defined_domain_nj(xios::CDomain* domain_hdl);
 bool cxios_is_defined_domain_ibegin(xios::CDomain* domain_hdl);
 bool cxios_is_defined_domain_jbegin(xios::CDomain* domain_hdl);
+bool cxios_is_defined_domain_lonvalue_1d(xios::CDomain* domain_hdl);
+bool cxios_is_defined_domain_latvalue_1d(xios::CDomain* domain_hdl);
 
 // grid group methods
 void cxios_gridgroup_handle_create(xios::CGridGroup** _ret, const char* _id, int _id_len);
diff --git a/core/src/iodef.xml.jinja b/core/src/iodef.xml.jinja
index 18b3db77f..3b93c85cb 100644
--- a/core/src/iodef.xml.jinja
+++ b/core/src/iodef.xml.jinja
@@ -9,9 +9,9 @@
 
     <!-- Axes defining extrusions for certain fields based on degrees of freedom per node -->
     <axis_definition>
-      <axis id="DGAxis" n_glo="{{ DGCOMP }}" name="dg_comp" dim_name="DGAxis"/>
-      <axis id="DGSAxis" n_glo="{{ DGSTRESSCOMP }}" name="dgstress_comp" dim_name="DGSAxis"/>
-      <axis id="VertexAxis" n_glo="2" name="ncoords" dim_name="VertexAxis"/>
+      <axis id="DGAxis" n_glo="{{ DGCOMP }}" name="dg_comp" dim_name="dg_comp"/>
+      <axis id="DGSAxis" n_glo="{{ DGSTRESSCOMP }}" name="dgstress_comp" dim_name="dgstress_comp"/>
+      <axis id="VertexAxis" n_glo="2" name="ncoords" dim_name="ncoords"/>
     </axis_definition>
 
     <grid_definition>
diff --git a/core/test/CMakeLists.txt b/core/test/CMakeLists.txt
index e84ca786d..08d3d6dce 100644
--- a/core/test/CMakeLists.txt
+++ b/core/test/CMakeLists.txt
@@ -12,7 +12,7 @@ set(ModulesRoot "../src/modules")
 
 # Generate netCDF files for tests
 execute_process(
-    COMMAND "${Python_EXECUTABLE}" "${CMAKE_CURRENT_SOURCE_DIR}/old_names.py"
+    COMMAND "${Python_EXECUTABLE}" "${CMAKE_CURRENT_SOURCE_DIR}/alt_names.py"
     WORKING_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}"
 )
 execute_process(
@@ -271,6 +271,9 @@ if(ENABLE_MPI)
         )
 
         add_mpi_test("${XIOS_TESTS}")
+
+        # Ensure XiosWrite_test is run before XiosRead_test
+        set_tests_properties(testXiosRead_MPI2 PROPERTIES DEPENDS testXiosWrite_MPI2)
     else()
         add_executable(testHaloExchangeCB_MPI3
             "HaloExchangeCB_test.cpp"
diff --git a/core/test/ParaGridIO_input_test.py b/core/test/ParaGridIO_input_test.py
index 379a8fd7b..b3241d2f1 100644
--- a/core/test/ParaGridIO_input_test.py
+++ b/core/test/ParaGridIO_input_test.py
@@ -37,10 +37,10 @@
     time_meta[:] = the_time
     time_meta.units = "seconds since 1970-01-01T00:00:00Z"
 
-    xDim = ncFile.createDimension("x", nx)
-    yDim = ncFile.createDimension("y", ny)
-    xVertexDim = ncFile.createDimension("xvertex", nx + 1)
-    yVertexDim = ncFile.createDimension("yvertex", ny + 1)
+    xDim = ncFile.createDimension("x_dim", nx)
+    yDim = ncFile.createDimension("y_dim", ny)
+    xVertexDim = ncFile.createDimension("x_vertex", nx + 1)
+    yVertexDim = ncFile.createDimension("y_vertex", ny + 1)
     xcg_dim = ncFile.createDimension("x_cg", nx * ncg + 1)
     ycg_dim = ncFile.createDimension("y_cg", ny * ncg + 1)
     dg_comp = ncFile.createDimension("dg_comp", n_dg)
@@ -49,9 +49,9 @@
     tDim = ncFile.createDimension("time", None)
 
     # Position and time variables
-    nc_lons = ncFile.createVariable("longitude", "f8", ("y", "x"))
+    nc_lons = ncFile.createVariable("longitude", "f8", ("y_dim", "x_dim"))
     nc_lons[:, :] = element_x
-    nc_lats = ncFile.createVariable("latitude", "f8", ("y", "x"))
+    nc_lats = ncFile.createVariable("latitude", "f8", ("y_dim", "x_dim"))
     nc_lats[:, :] = element_y
 
     nc_times = ncFile.createVariable("time", "f8", ("time"))
@@ -59,7 +59,7 @@
     # (unix_times_e, era5_times) = create_era5_times(start_time, stop_time)
     # For each field and time, get the corresponding file name for each dataset
     # get the source data
-    data = ncFile.createVariable(field_name, "f8", ("time", "y", "x"))
+    data = ncFile.createVariable(field_name, "f8", ("time", "y_dim", "x_dim"))
     data[0, :, :] = time_data
     # 'fill' the time axis
     nc_times[0] = the_time
diff --git a/core/test/ParaGrid_test.cpp b/core/test/ParaGrid_test.cpp
index 06dd30c1a..fbdd3ea9c 100644
--- a/core/test/ParaGrid_test.cpp
+++ b/core/test/ParaGrid_test.cpp
@@ -489,12 +489,12 @@ TEST_CASE("Check an exception is thrown for an invalid file name")
 }
 
 #ifdef USE_MPI
-MPI_TEST_CASE("Check if a file with the old dimension names can be read", 2)
+MPI_TEST_CASE("Check if a file with the alternative dimension names can be read", 2)
 #else
-TEST_CASE("Check if a file with the old dimension names can be read")
+TEST_CASE("Check if a file with the alternative dimension names can be read")
 #endif
 {
-    std::string inputFilename = std::string(TEST_FILE_SOURCE) + "/old_names.nc";
+    std::string inputFilename = std::string(TEST_FILE_SOURCE) + "/alt_names.nc";
 
     Module::setImplementation<IStructure>("Nextsim::ParametricGrid");
 
diff --git a/core/test/XiosRead_test.cpp b/core/test/XiosRead_test.cpp
index a71c17661..f7385ac1c 100644
--- a/core/test/XiosRead_test.cpp
+++ b/core/test/XiosRead_test.cpp
@@ -20,7 +20,7 @@
 #include <filesystem>
 
 const std::string testFilesDir = TEST_FILES_DIR;
-const std::string restartFilename = testFilesDir + "/xios_test_input.nc";
+const std::string restartFilename = testFilesDir + "/xios_test_output.nc";
 const std::string forcingFilename = testFilesDir + "/xios_test_forcing.nc";
 
 static const int DGCOMP = 6;
@@ -79,8 +79,11 @@ MPI_TEST_CASE("TestXiosRead", 2)
 
     xiosHandler.close_context_definition();
 
-    // Check the input files exists
-    REQUIRE(std::filesystem::exists(restartFilename));
+    // Check the input files exist
+    if (!std::filesystem::exists(restartFilename)) {
+        throw std::runtime_error(
+            "XiosRead_test: Input file not found. Did you run XiosWrite_test?");
+    }
     REQUIRE(std::filesystem::exists(forcingFilename));
 
     // Check calendar step is zero initially
@@ -169,6 +172,14 @@ MPI_TEST_CASE("TestXiosRead", 2)
         REQUIRE(xiosHandler.getCalendarStep() == ts + 1);
     }
 
+    if (rank == 0) {
+        // TODO: Re-enable file splitting (#898)
+        // std::filesystem::remove("xios_test_output_20230317171100-20230317201059.nc");
+        // std::filesystem::remove("xios_test_output_20230317201100-20230317231059.nc");
+        std::filesystem::remove("xios_test_output.nc");
+        std::filesystem::remove("xios_test_diagnostic.nc");
+    }
+
     xiosHandler.context_finalize();
     Finalizer::finalize();
 }
diff --git a/core/test/XiosWrite_test.cpp b/core/test/XiosWrite_test.cpp
index 33f4697e7..5a6069a47 100644
--- a/core/test/XiosWrite_test.cpp
+++ b/core/test/XiosWrite_test.cpp
@@ -51,7 +51,8 @@ MPI_TEST_CASE("TestXiosWrite", 2)
     config << "field_names = " << maskName << "," << coordsName << "," << hiceName << ","
            << ticeName << "," << uName << std::endl;
     config << "period = P0-0T01:30:00" << std::endl;
-    config << "split_period = P0-0T03:00:00" << std::endl;
+    // TODO: Re-enable file splitting (#898)
+    // config << "split_period = P0-0T03:00:00" << std::endl;
     config << "[XiosDiagnostic]" << std::endl;
     config << "filename = " << diagnosticFilename << std::endl;
     config << "field_names = " << hsnowName << std::endl;
@@ -194,15 +195,13 @@ MPI_TEST_CASE("TestXiosWrite", 2)
         grid.dumpModelState(restarts, restartOutputFilename, true);
     }
 
-    // Check the files have indeed been created then remove it
-    REQUIRE(std::filesystem::exists("xios_test_output_20230317171100-20230317201059.nc"));
-    REQUIRE(std::filesystem::exists("xios_test_output_20230317201100-20230317231059.nc"));
+    // Check the files have indeed been created
+    // NOTE: We don't remove them because they are used in XiosRead_test
+    // TODO: Re-enable file splitting (#898)
+    // REQUIRE(std::filesystem::exists("xios_test_output_20230317171100-20230317201059.nc"));
+    // REQUIRE(std::filesystem::exists("xios_test_output_20230317201100-20230317231059.nc"));
+    REQUIRE(std::filesystem::exists("xios_test_output.nc"));
     REQUIRE(std::filesystem::exists("xios_test_diagnostic.nc"));
-    if (rank == 0) {
-        std::filesystem::remove("xios_test_output_20230317171100-20230317201059.nc");
-        std::filesystem::remove("xios_test_output_20230317201100-20230317231059.nc");
-        std::filesystem::remove("xios_test_diagnostic.nc");
-    }
 
     xiosHandler.context_finalize();
     Finalizer::finalize();
diff --git a/core/test/old_names.py b/core/test/alt_names.py
similarity index 80%
rename from core/test/old_names.py
rename to core/test/alt_names.py
index 1fe6226d4..a8919e37a 100644
--- a/core/test/old_names.py
+++ b/core/test/alt_names.py
@@ -1,4 +1,4 @@
-"""Generate the old_names.nc file used by core/test/ParaGrid_test.cpp."""
+"""Generate the alt_names.nc file used by core/test/ParaGrid_test.cpp."""
 
 import time
 
@@ -16,7 +16,7 @@
     n_dgstress = 1
     n_coords = 2
 
-    ncFile = netCDF4.Dataset("old_names.nc", "w", format="NETCDF4")
+    ncFile = netCDF4.Dataset("alt_names.nc", "w", format="NETCDF4")
 
     ncFile.structure_name = "parametric_rectangular"
 
@@ -28,8 +28,8 @@
     formatted.format = "%Y-%m-%dT%H:%M:%SZ"
     formatted[0] = "2010-01-01T00:00:00Z"
 
-    yDim = ncFile.createDimension("y", nfirst)
-    xDim = ncFile.createDimension("x", nsecond)
+    yDim = ncFile.createDimension("ydim", nfirst)
+    xDim = ncFile.createDimension("xdim", nsecond)
     yVertexDim = ncFile.createDimension("yvertex", nfirst + 1)
     xVertexDim = ncFile.createDimension("xvertex", nsecond + 1)
     ycg_dim = ncFile.createDimension("y_cg", nfirst * ncg + 1)
@@ -38,21 +38,21 @@
     dgs_comp = ncFile.createDimension("dgstress_comp", n_dgstress)
     n_coords_comp = ncFile.createDimension("ncoords", n_coords)
 
-    field_dims = ("y", "x")
+    field_dims = ("ydim", "xdim")
     coord_dims = ("yvertex", "xvertex", "ncoords")
 
     # Array coordinates
     space = 25000.0  # 25 km in metres
 
-    x_vertex1d = np.arange(0, space * (nsecond + 1), space)
-    y_vertex1d = np.arange(0, space * (nfirst + 1), space)
+    xvertex1d = np.arange(0, space * (nsecond + 1), space)
+    yvertex1d = np.arange(0, space * (nfirst + 1), space)
 
-    x_vertex = np.resize(x_vertex1d, (nfirst + 1, nsecond + 1))
-    y_vertex = np.resize(y_vertex1d, (nsecond + 1, nfirst + 1)).transpose()
+    xvertex = np.resize(xvertex1d, (nfirst + 1, nsecond + 1))
+    yvertex = np.resize(yvertex1d, (nsecond + 1, nfirst + 1)).transpose()
 
     coords = ncFile.createVariable("coords", "f8", coord_dims)
-    coords[:, :, 0] = x_vertex
-    coords[:, :, 1] = y_vertex
+    coords[:, :, 0] = xvertex
+    coords[:, :, 1] = yvertex
 
     x_cell1d = np.arange(space * 0.5, space * (nsecond + 0.5), space)
     y_cell1d = np.arange(space * 0.5, space * (nfirst + 0.5), space)
@@ -60,9 +60,9 @@
     x_cell = np.resize(x_cell1d, (nfirst, nsecond))
     y_cell = np.resize(y_cell1d, (nsecond, nfirst)).transpose()
 
-    elem_x = ncFile.createVariable("x", "f8", field_dims)
+    elem_x = ncFile.createVariable("xdim", "f8", field_dims)
     elem_x[:, :] = x_cell
-    elem_y = ncFile.createVariable("y", "f8", field_dims)
+    elem_y = ncFile.createVariable("ydim", "f8", field_dims)
     elem_y[:, :] = y_cell
 
     grid_azimuth = ncFile.createVariable("grid_azimuth", "f8", field_dims)
diff --git a/core/test/testmodelarraydetails/ModelArrayDetails.cpp b/core/test/testmodelarraydetails/ModelArrayDetails.cpp
index 4e61f331e..06d984246 100644
--- a/core/test/testmodelarraydetails/ModelArrayDetails.cpp
+++ b/core/test/testmodelarraydetails/ModelArrayDetails.cpp
@@ -1,5 +1,4 @@
 /*!
- *
  * @author  Tim Spain <timothy.spain@nersc.no>
  */
 
diff --git a/core/test/xios_test_forcing.cdl b/core/test/xios_test_forcing.cdl
index 3a0c048d1..4f6f1f754 100644
--- a/core/test/xios_test_forcing.cdl
+++ b/core/test/xios_test_forcing.cdl
@@ -1,10 +1,10 @@
 netcdf xios_test_forcing {
 dimensions:
 	axis_nbounds = 2 ;
-	x = 4 ;
-	y = 2 ;
-  xvertex = 5 ;
-  yvertex = 3 ;
+	x_dim = 4 ;
+	y_dim = 2 ;
+  x_vertex = 5 ;
+  y_vertex = 3 ;
   x_cg = 9 ;
   y_cg = 5 ;
   dg_comp = 3 ;
@@ -12,16 +12,16 @@ dimensions:
   ncoords = 2 ;
 	time_counter = UNLIMITED ; // (5 currently)
 variables:
-	float x(x) ;
-		x:axis = "X" ;
-		x:standard_name = "longitude" ;
-		x:long_name = "Longitude" ;
-		x:units = "degrees_east" ;
-	float y(y) ;
-		y:axis = "Y" ;
-		y:standard_name = "latitude" ;
-		y:long_name = "Latitude" ;
-		y:units = "degrees_north" ;
+	float x_dim(x_dim) ;
+		x_dim:axis = "X" ;
+		x_dim:standard_name = "longitude" ;
+		x_dim:long_name = "Longitude" ;
+		x_dim:units = "degrees_east" ;
+	float y_dim(y_dim) ;
+		y_dim:axis = "Y" ;
+		y_dim:standard_name = "latitude" ;
+		y_dim:long_name = "Latitude" ;
+		y_dim:units = "degrees_north" ;
 	double time_instant(time_counter) ;
 		time_instant:standard_name = "time" ;
 		time_instant:long_name = "Time axis" ;
@@ -39,7 +39,7 @@ variables:
 		time_counter:time_origin = "2020-01-23 00:08:15" ;
 		time_counter:bounds = "time_counter_bounds" ;
 	double time_counter_bounds(time_counter, axis_nbounds) ;
-	float hsnow(time_counter, y, x) ;
+	float hsnow(time_counter, y_dim, x_dim) ;
 		hsnow:online_operation = "instant" ;
 		hsnow:interval_operation = "5400 s" ;
 		hsnow:interval_write = "5400 s" ;
@@ -55,10 +55,6 @@ variables:
 
 data:
 
- y = -1, 0 ;
-
- x = -1, -0.5, 0, 0.5 ;
-
  time_instant = 99334965, 99340365, 99345765, 99351165, 99356565 ;
 
  time_instant_bounds =
diff --git a/core/test/xios_test_input.cdl b/core/test/xios_test_input.cdl
index 36673bdfc..49720846d 100644
--- a/core/test/xios_test_input.cdl
+++ b/core/test/xios_test_input.cdl
@@ -1,36 +1,36 @@
 netcdf xios_test_input {
 dimensions:
 	axis_nbounds = 2 ;
-	x = 4 ;
-	y = 2 ;
-  xvertex = 5 ;
-  yvertex = 3 ;
+	x_dim = 4 ;
+	y_dim = 2 ;
+  x_vertex = 5 ;
+  y_vertex = 3 ;
   x_cg = 9 ;
   y_cg = 5 ;
   dg_comp = 6 ;
   dgstress_comp = 8 ;
   ncoords = 2 ;
 variables:
-	float x(x) ;
-		x:axis = "X" ;
-		x:standard_name = "longitude" ;
-		x:long_name = "Longitude" ;
-		x:units = "degrees_east" ;
-	float y(y) ;
-		y:axis = "Y" ;
-		y:standard_name = "latitude" ;
-		y:long_name = "Latitude" ;
-		y:units = "degrees_north" ;
-	float xvertex(xvertex) ;
-		xvertex:axis = "XVERTEX" ;
-		xvertex:standard_name = "vertex longitude" ;
-		xvertex:long_name = "Vertex longitude" ;
-		xvertex:units = "degrees_east" ;
-	float yvertex(yvertex) ;
-		yvertex:axis = "YVERTEX" ;
-		yvertex:standard_name = "vertex latitude" ;
-		yvertex:long_name = "Vertex latitude" ;
-		yvertex:units = "degrees_north" ;
+	float x_dim(x_dim) ;
+		x_dim:axis = "X" ;
+		x_dim:standard_name = "longitude" ;
+		x_dim:long_name = "Longitude" ;
+		x_dim:units = "degrees_east" ;
+	float y_dim(y_dim) ;
+		y_dim:axis = "Y" ;
+		y_dim:standard_name = "latitude" ;
+		y_dim:long_name = "Latitude" ;
+		y_dim:units = "degrees_north" ;
+	float x_vertex(x_vertex) ;
+		x_vertex:axis = "x_vertex" ;
+		x_vertex:standard_name = "vertex longitude" ;
+		x_vertex:long_name = "Vertex longitude" ;
+		x_vertex:units = "degrees_east" ;
+	float y_vertex(y_vertex) ;
+		y_vertex:axis = "y_vertex" ;
+		y_vertex:standard_name = "vertex latitude" ;
+		y_vertex:long_name = "Vertex latitude" ;
+		y_vertex:units = "degrees_north" ;
 	float x_cg(x_cg) ;
 		x_cg:axis = "XCG" ;
 		x_cg:standard_name = "CG longitude" ;
@@ -56,16 +56,16 @@ variables:
 		VertexAxis:standard_name = "Dimension" ;
     VertexAxis:long_name = "Spatial dimension" ;
 		VertexAxis:units = "-" ;
-	float coords(yvertex, xvertex, ncoords) ;
+	float coords(y_vertex, x_vertex, ncoords) ;
 		coords:online_operation = "once" ;
 		coords:coordinates = "lat lon VertexAxis" ;
-  float hice(y, x, dg_comp) ;
+  float hice(y_dim, x_dim, dg_comp) ;
 		hice:online_operation = "once" ;
 		hice:coordinates = "lat lon DGAxis" ;
-	float mask(y, x) ;
+	float mask(y_dim, x_dim) ;
 		mask:online_operation = "once" ;
 		mask:coordinates = "lat lon" ;
-	float tice(y, x, dgstress_comp) ;
+	float tice(y_dim, x_dim, dgstress_comp) ;
 		tice:online_operation = "once" ;
 		tice:coordinates = "lat lon DGSAxis" ;
 	float u(y_cg, x_cg) ;
@@ -81,14 +81,6 @@ variables:
 
 data:
 
- y = -1, 0 ;
-
- x = -1, -0.5, 0, 0.5 ;
-
- yvertex = -1, 0, 1 ;
-
- xvertex = -1, -0.5, 0, 0.5, 1 ;
-
  coords =
   0, 0,
   1, 0,
diff --git a/physics/test/era5_test_data.py b/physics/test/era5_test_data.py
index c4b369cb6..3c36418a2 100644
--- a/physics/test/era5_test_data.py
+++ b/physics/test/era5_test_data.py
@@ -15,10 +15,10 @@
 
 ncFile.structure_name = "parametric_rectangular"
 
-xDim = ncFile.createDimension("x", nx)
-yDim = ncFile.createDimension("y", ny)
-xVertexDim = ncFile.createDimension("xvertex", nx + 1)
-yVertexDim = ncFile.createDimension("yvertex", ny + 1)
+xDim = ncFile.createDimension("x_dim", nx)
+yDim = ncFile.createDimension("y_dim", ny)
+xVertexDim = ncFile.createDimension("x_vertex", nx + 1)
+yVertexDim = ncFile.createDimension("y_vertex", ny + 1)
 n_coords_comp = ncFile.createDimension("ncoords", n_coords)
 time_dim = ncFile.createDimension("time", None)
 
@@ -29,8 +29,8 @@
 formatted.format = "%Y-%m-%dT%H:%M:%SZ"
 formatted[0] = "2000-01-01T00:00:00Z"
 
-hfield_dims = ("y", "x")
-timefield_dims = ("time", "y", "x")
+hfield_dims = ("y_dim", "x_dim")
+timefield_dims = ("time", "y_dim", "x_dim")
 
 # fmt: off
 mask33 = np.array(
@@ -98,7 +98,7 @@
 lat = 90 - (x_coords**2 + y_coords**2) ** 0.5
 lon = np.rad2deg(np.arctan2(y_coords, x_coords))
 
-coords = ncFile.createVariable("coords", "f8", ("yvertex", "xvertex", "ncoords"))
+coords = ncFile.createVariable("coords", "f8", ("y_vertex", "x_vertex", "ncoords"))
 coords[:, :, 0] = lon
 coords[:, :, 1] = lat
 
diff --git a/physics/test/topaz_test_data.py b/physics/test/topaz_test_data.py
index d3002e27a..1e7ae4f55 100644
--- a/physics/test/topaz_test_data.py
+++ b/physics/test/topaz_test_data.py
@@ -15,10 +15,10 @@
 
 ncFile.structure_name = "parametric_rectangular"
 
-xDim = ncFile.createDimension("x", nx)
-yDim = ncFile.createDimension("y", ny)
-xVertexDim = ncFile.createDimension("xvertex", nx + 1)
-yVertexDim = ncFile.createDimension("yvertex", ny + 1)
+xDim = ncFile.createDimension("x_dim", nx)
+yDim = ncFile.createDimension("y_dim", ny)
+xVertexDim = ncFile.createDimension("x_vertex", nx + 1)
+yVertexDim = ncFile.createDimension("y_vertex", ny + 1)
 n_coords_comp = ncFile.createDimension("ncoords", n_coords)
 time_dim = ncFile.createDimension("time", None)
 
@@ -29,8 +29,8 @@
 formatted.format = "%Y-%m-%dT%H:%M:%SZ"
 formatted[0] = "2000-01-01T00:00:00Z"
 
-hfield_dims = ("y", "x")
-timefield_dims = ("time", "y", "x")
+hfield_dims = ("y_dim", "x_dim")
+timefield_dims = ("time", "y_dim", "x_dim")
 
 # fmt: off
 mask33 = np.array(
@@ -98,7 +98,7 @@
 lat = 90 - (x_coords**2 + y_coords**2) ** 0.5
 lon = np.rad2deg(np.arctan2(y_coords, x_coords))
 
-coords = ncFile.createVariable("coords", "f8", ("yvertex", "xvertex", "ncoords"))
+coords = ncFile.createVariable("coords", "f8", ("y_vertex", "x_vertex", "ncoords"))
 coords[:, :, 0] = lon
 coords[:, :, 1] = lat
 
diff --git a/run/era5_topaz4_maker.py b/run/era5_topaz4_maker.py
index 93ae53f73..86c38f5eb 100644
--- a/run/era5_topaz4_maker.py
+++ b/run/era5_topaz4_maker.py
@@ -238,12 +238,12 @@ def topaz4_source_file_name(unix_time, path):
     era5_ncFile = netCDF4.Dataset(era5_out_file, "w", format="NETCDF4")
     era5_ncFile.structure_name = target_structure
 
-    xDim = era5_ncFile.createDimension("xdim", nx)
-    yDim = era5_ncFile.createDimension("ydim", ny)
+    xDim = era5_ncFile.createDimension("x_dim", nx)
+    yDim = era5_ncFile.createDimension("y_dim", ny)
     tDim = era5_ncFile.createDimension("time", None)
 
-    hfield_dims = ("ydim", "xdim")
-    timefield_dims = ("time", "ydim", "xdim")
+    hfield_dims = ("y_dim", "x_dim")
+    timefield_dims = ("time", "y_dim", "x_dim")
 
     # Position and time variables
     nc_lons = era5_ncFile.createVariable("longitude", "f8", hfield_dims)
@@ -352,8 +352,8 @@ def topaz4_source_file_name(unix_time, path):
     print(f"Writing TOPAZ4 data to {topaz_out_file}")
     topaz_ncFile.structure_name = target_structure
 
-    xDim = topaz_ncFile.createDimension("xdim", nx)
-    yDim = topaz_ncFile.createDimension("ydim", ny)
+    xDim = topaz_ncFile.createDimension("x_dim", nx)
+    yDim = topaz_ncFile.createDimension("y_dim", ny)
     tDim = topaz_ncFile.createDimension("time", None)
 
     (unix_times_t, topaz4_times) = create_topaz_times(start_time, stop_time)
@@ -401,8 +401,8 @@ def topaz4_source_file_name(unix_time, path):
                 # Now interpolate the source data to the target grid
                 time_data = np.zeros((nx, ny))
                 proj_string = source_file["stereographic"].proj4
-                source_x = source_file["x"][:]
-                source_y = source_file["y"][:]
+                source_x = source_file["x_dim"][:]
+                source_y = source_file["y_dim"][:]
                 time_data = topaz4_interpolate(
                     element_lon,
                     element_lat,
@@ -441,8 +441,8 @@ def topaz4_source_file_name(unix_time, path):
         v_source_data_tgrid = np.zeros((nx, ny))
         # Interpolate the current components on the TOPAZ basis on to the new grid
         proj_string = source_file["stereographic"].proj4
-        source_x = source_file["x"][:]
-        source_y = source_file["y"][:]
+        source_x = source_file["x_dim"][:]
+        source_y = source_file["y_dim"][:]
         u_source_data_tgrid = topaz4_interpolate(
             element_lon, element_lat, u_source_data, source_x, source_y, proj_string
         )
diff --git a/run/init_para128x128.py b/run/init_para128x128.py
index 85b849b4f..1312e4d00 100644
--- a/run/init_para128x128.py
+++ b/run/init_para128x128.py
@@ -23,18 +23,18 @@
 formatted.format = "%Y-%m-%dT%H:%M:%SZ"
 formatted[0] = "2000-01-01T00:00:00Z"
 
-xDim = ncFile.createDimension("x", nx)
-yDim = ncFile.createDimension("y", ny)
+xDim = ncFile.createDimension("x_dim", nx)
+yDim = ncFile.createDimension("y_dim", ny)
 nLay = ncFile.createDimension("nLayers", nLayers)
-xVertexDim = ncFile.createDimension("xvertex", nx + 1)
-yVertexDim = ncFile.createDimension("yvertex", ny + 1)
+xVertexDim = ncFile.createDimension("x_vertex", nx + 1)
+yVertexDim = ncFile.createDimension("y_vertex", ny + 1)
 xcg_dim = ncFile.createDimension("x_cg", nx * ncg + 1)
 ycg_dim = ncFile.createDimension("y_cg", ny * ncg + 1)
 dg_comp = ncFile.createDimension("dg_comp", n_dg)
 dgs_comp = ncFile.createDimension("dgstress_comp", n_dgstress)
 n_coords_comp = ncFile.createDimension("ncoords", n_coords)
 
-hfield_dims = ("y", "x")
+hfield_dims = ("y_dim", "x_dim")
 
 # fmt: off
 mask33 = np.array(
@@ -108,7 +108,7 @@
 lon %= 360.0
 lon -= 180.0
 
-coords = ncFile.createVariable("coords", "f8", ("yvertex", "xvertex", "ncoords"))
+coords = ncFile.createVariable("coords", "f8", ("y_vertex", "x_vertex", "ncoords"))
 coords[:, :, 0] = lon
 coords[:, :, 1] = lat
 
@@ -165,7 +165,7 @@
 hice[:, :] = cice[:, :] * 2
 hsnow = ncFile.createVariable("hsnow", "f8", hfield_dims)
 hsnow[:, :] = cice[:, :] / 2
-tice = ncFile.createVariable("tice", "f8", ("nLayers", "y", "x"))
+tice = ncFile.createVariable("tice", "f8", ("nLayers", "y_dim", "x_dim"))
 tice[0, :, :] = -0.5 - cice[:, :]
 tice[1, :, :] = -1.5 - cice[:, :]
 tice[2, :, :] = -2.5 - cice[:, :]
diff --git a/run/init_topaz128x128.py b/run/init_topaz128x128.py
index c7653724d..a1eb1e783 100644
--- a/run/init_topaz128x128.py
+++ b/run/init_topaz128x128.py
@@ -92,18 +92,18 @@ def topaz4_interpolate(target_lon_deg, target_lat_deg, data, lat_array):
     formatted.format = "%Y-%m-%dT%H:%M:%SZ"
     formatted[0] = "2010-01-01T00:00:00Z"
 
-    xDim = ncFile.createDimension("x", nx)
-    yDim = ncFile.createDimension("y", ny)
-    nLay = ncFile.createDimension("z", nLayers)
-    xVertexDim = ncFile.createDimension("xvertex", nx + 1)
-    yVertexDim = ncFile.createDimension("yvertex", ny + 1)
+    xDim = ncFile.createDimension("x_dim", nx)
+    yDim = ncFile.createDimension("y_dim", ny)
+    nLay = ncFile.createDimension("nLayers", nLayers)
+    xVertexDim = ncFile.createDimension("x_vertex", nx + 1)
+    yVertexDim = ncFile.createDimension("y_vertex", ny + 1)
     xcg_dim = ncFile.createDimension("x_cg", nx * ncg + 1)
     ycg_dim = ncFile.createDimension("y_cg", ny * ncg + 1)
     dg_comp = ncFile.createDimension("dg_comp", n_dg)
     dgs_comp = ncFile.createDimension("dgstress_comp", n_dgstress)
     n_coords_comp = ncFile.createDimension("ncoords", n_coords)
 
-    hfield_dims = ("y", "x")
+    hfield_dims = ("y_dim", "x_dim")
 
     # Array coordinates
     array_size1d = 20.0  # About 20˚ square
@@ -156,7 +156,7 @@ def topaz4_interpolate(target_lon_deg, target_lat_deg, data, lat_array):
     nc_lats = ncFile.createVariable("latitude", "f8", hfield_dims)
     nc_lats[:, :] = element_lat
 
-    coords = ncFile.createVariable("coords", "f8", ("yvertex", "xvertex", "ncoords"))
+    coords = ncFile.createVariable("coords", "f8", ("y_vertex", "x_vertex", "ncoords"))
     coords[:, :, 0] = lon
     coords[:, :, 1] = lat
 
@@ -206,7 +206,7 @@ def topaz4_interpolate(target_lon_deg, target_lat_deg, data, lat_array):
     mu = -0.055
 
     # Ice temperature
-    tice = ncFile.createVariable("tice", "f8", ("z", "y", "x"))
+    tice = ncFile.createVariable("tice", "f8", ("nLayers", "y_dim", "x_dim"))
     ice_melt = mu * 5  # Melting point of sea ice (salinity = 5) in ˚C
     # Tice outside the ice pack is the melting point of pure water ice, which is conveniently 0˚C
     ice_temp2d = np.fmin(sst_data, ice_melt) * isice
diff --git a/run/make_init.py b/run/make_init.py
index f0724af54..6ca64c9e9 100644
--- a/run/make_init.py
+++ b/run/make_init.py
@@ -19,13 +19,13 @@
 formatted.format = "%Y-%m-%dT%H:%M:%SZ"
 formatted[0] = "2000-01-01T00:00:00Z"
 
-x_dim = ncFile.createDimension("xdim", nx)
-y_dim = ncFile.createDimension("ydim", ny)
-xvertex_dim = ncFile.createDimension("xvertex", nx + 1)
-yvertex_dim = ncFile.createDimension("yvertex", ny + 1)
+x_dim = ncFile.createDimension("x_dim", nx)
+y_dim = ncFile.createDimension("y_dim", ny)
+xvertex_dim = ncFile.createDimension("x_vertex", nx + 1)
+yvertex_dim = ncFile.createDimension("y_vertex", ny + 1)
 coords_dim = ncFile.createDimension("ncoords", n_coords)
 
-hfield_dims = ("ydim", "xdim")
+hfield_dims = ("y_dim", "x_dim")
 
 mask = ncFile.createVariable("mask", "f8", hfield_dims)
 # fmt: off
@@ -113,8 +113,8 @@
 
 # coordinates
 # element centres
-x_var = ncFile.createVariable("x", "f8", hfield_dims)
-y_var = ncFile.createVariable("y", "f8", hfield_dims)
+x_var = ncFile.createVariable("x_dim", "f8", hfield_dims)
+y_var = ncFile.createVariable("y_dim", "f8", hfield_dims)
 
 d_distance = 150000  # 150 km element spacing
 
@@ -125,7 +125,7 @@
         x_var[j, i] = x
         y_var[j, i] = y
 
-coords = ncFile.createVariable("coords", "f8", ("yvertex", "xvertex", "ncoords"))
+coords = ncFile.createVariable("coords", "f8", ("y_vertex", "x_vertex", "ncoords"))
 
 for j in range(ny + 1):
     y = j * d_distance
diff --git a/run/make_init25kmNH_test.py b/run/make_init25kmNH_test.py
index 3dbf52aa8..d5112475e 100644
--- a/run/make_init25kmNH_test.py
+++ b/run/make_init25kmNH_test.py
@@ -5,8 +5,8 @@
 
 grid = netCDF4.Dataset("25km_NH.nc", "r")
 
-nx = grid.dimensions["x"].size
-ny = grid.dimensions["y"].size
+nx = grid.dimensions["x_dim"].size
+ny = grid.dimensions["y_dim"].size
 ncg = 1
 n_dg = 1
 n_dgstress = 1
@@ -24,10 +24,10 @@
 formatted.format = "%Y-%m-%dT%H:%M:%SZ"
 formatted[0] = "2000-01-01T00:00:00Z"
 
-xDim = ncFile.createDimension("xdim", nx)
-yDim = ncFile.createDimension("ydim", ny)
-xVertexDim = ncFile.createDimension("xvertex", nx + 1)
-yVertexDim = ncFile.createDimension("yvertex", ny + 1)
+xDim = ncFile.createDimension("x_dim", nx)
+yDim = ncFile.createDimension("y_dim", ny)
+xVertexDim = ncFile.createDimension("x_vertex", nx + 1)
+yVertexDim = ncFile.createDimension("y_vertex", ny + 1)
 xcg_dim = ncFile.createDimension("x_cg", nx * ncg + 1)
 ycg_dim = ncFile.createDimension("y_cg", ny * ncg + 1)
 dg_comp = ncFile.createDimension("dg_comp", n_dg)
@@ -36,7 +36,7 @@
 
 grid_mask = grid["mask"]
 
-mask = ncFile.createVariable("mask", "f8", ("xdim", "ydim"))
+mask = ncFile.createVariable("mask", "f8", ("x_dim", "y_dim"))
 mask[:, :] = grid_mask[:, :]
 antimask = 1 - mask[:, :]
 
@@ -53,7 +53,7 @@
 node_lat[-1, -1] = grid["lat_corners"][-1, -1, 2]
 node_lat[-1, 0:-1] = grid["lat_corners"][-1, :, 3]
 
-coords = ncFile.createVariable("coords", "f8", ("xvertex", "yvertex", "ncoords"))
+coords = ncFile.createVariable("coords", "f8", ("x_vertex", "y_vertex", "ncoords"))
 coords[:, :, 0] = node_lon
 coords[:, :, 1] = node_lat
 
@@ -61,8 +61,8 @@
     "longitude",
     "f8",
     (
-        "xdim",
-        "ydim",
+        "x_dim",
+        "y_dim",
     ),
 )
 elem_lon[:, :] = grid["plon"][:, :]
@@ -70,8 +70,8 @@
     "latitude",
     "f8",
     (
-        "xdim",
-        "ydim",
+        "x_dim",
+        "y_dim",
     ),
 )
 elem_lat[:, :] = grid["plat"][:, :]
@@ -81,8 +81,8 @@
     "cice",
     "f8",
     (
-        "xdim",
-        "ydim",
+        "x_dim",
+        "y_dim",
     ),
 )
 cice[:, :] = mask[:, :] * 0.95
@@ -90,8 +90,8 @@
     "hice",
     "f8",
     (
-        "xdim",
-        "ydim",
+        "x_dim",
+        "y_dim",
     ),
 )
 hice[:, :] = cice[:, :] * 2
@@ -99,24 +99,24 @@
     "hsnow",
     "f8",
     (
-        "xdim",
-        "ydim",
+        "x_dim",
+        "y_dim",
     ),
 )
 hsnow[:, :] = cice[:, :] / 2
-tsurf = ncFile.createVariable("tsurf", "f8", ("xdim", "ydim"))
+tsurf = ncFile.createVariable("tsurf", "f8", ("x_dim", "y_dim"))
 tsurf[:, :] = -0.5 - cice[:, :]
-tbott = ncFile.createVariable("tbottom", "f8", ("xdim", "ydim"))
+tbott = ncFile.createVariable("tbottom", "f8", ("x_dim", "y_dim"))
 tbott[:, :] = -1.8
-tintr = ncFile.createVariable("tinterior", "f8", ("xdim", "ydim"))
+tintr = ncFile.createVariable("tinterior", "f8", ("x_dim", "y_dim"))
 tintr[:, :] = 0.5 * (tsurf[:, :] + tbott[:, :])
 
 sst = ncFile.createVariable(
     "sst",
     "f8",
     (
-        "xdim",
-        "ydim",
+        "x_dim",
+        "y_dim",
     ),
 )
 sst[:, :] = -cice[:, :]
@@ -124,8 +124,8 @@
     "sss",
     "f8",
     (
-        "xdim",
-        "ydim",
+        "x_dim",
+        "y_dim",
     ),
 )
 sss[:, :] = cice[:, :] * 33.68
@@ -133,8 +133,8 @@
     "u",
     "f8",
     (
-        "xdim",
-        "ydim",
+        "x_dim",
+        "y_dim",
     ),
 )
 u[:, :] = 0.0
@@ -142,8 +142,8 @@
     "v",
     "f8",
     (
-        "xdim",
-        "ydim",
+        "x_dim",
+        "y_dim",
     ),
 )
 v[:, :] = 0.0
diff --git a/run/make_init_base.py b/run/make_init_base.py
index e1cc61aa2..af6a727d4 100644
--- a/run/make_init_base.py
+++ b/run/make_init_base.py
@@ -86,18 +86,18 @@ def _init_vars_and_file(self, nFirst, nSecond):
         structure_name = "parametric_rectangular"
         self._ncFile.structure_name = structure_name
 
-        self._ncFile.createDimension("ydim", nFirst)
-        self._ncFile.createDimension("xdim", nSecond)
-        self._ncFile.createDimension("yvertex", nFirst + 1)
-        self._ncFile.createDimension("xvertex", nSecond + 1)
+        self._ncFile.createDimension("y_dim", nFirst)
+        self._ncFile.createDimension("x_dim", nSecond)
+        self._ncFile.createDimension("y_vertex", nFirst + 1)
+        self._ncFile.createDimension("x_vertex", nSecond + 1)
         self._ncFile.createDimension("y_cg", nFirst * self._nCg + 1)
         self._ncFile.createDimension("x_cg", nSecond * self._nCg + 1)
         self._ncFile.createDimension("dg_comp", self._nDg)
         self._ncFile.createDimension("dgstress_comp", self._nDgStress)
         self._ncFile.createDimension("ncoords", self._nCoords)
 
-        self._field_dims = ("ydim", "xdim")
-        self._coord_dims = ("yvertex", "xvertex", "ncoords")
+        self._field_dims = ("y_dim", "x_dim")
+        self._coord_dims = ("y_vertex", "x_vertex", "ncoords")
 
     def make_cartesian_grid(self, nFirst, nSecond, res):
         """
@@ -129,9 +129,9 @@ def make_cartesian_grid(self, nFirst, nSecond, res):
                 px[j, i] = (j + 0.5) * res
                 py[j, i] = (i + 0.5) * res
 
-        elem_x = self._ncFile.createVariable("x", "f8", self._field_dims)
+        elem_x = self._ncFile.createVariable("x_dim", "f8", self._field_dims)
         elem_x[:, :] = px
-        elem_y = self._ncFile.createVariable("y", "f8", self._field_dims)
+        elem_y = self._ncFile.createVariable("y_dim", "f8", self._field_dims)
         elem_y[:, :] = py
 
         grid_azimuth = self._ncFile.createVariable(
@@ -162,8 +162,8 @@ def make_geographic_grid(
         grid = netCDF4.Dataset(f"{grid_file}", "r")
 
         # Grid dimensions. We're dealing with files written in and for FORTRAN, so y is the first dimension and x the second.
-        nfirst = grid.dimensions["y"].size
-        nsecond = grid.dimensions["x"].size
+        nfirst = grid.dimensions["y_dim"].size
+        nsecond = grid.dimensions["x_dim"].size
 
         # Initialise the arrays and the output file
         self._init_vars_and_file(nfirst, nsecond)
diff --git a/run/make_init_geographic.py b/run/make_init_geographic.py
index efebdc40d..5998cd341 100644
--- a/run/make_init_geographic.py
+++ b/run/make_init_geographic.py
@@ -117,8 +117,8 @@ def topaz4_source_file_name(time_struct):
     source_file_name = topaz4_source_file_name(data_time)
     source_file = netCDF4.Dataset(source_file_name, "r")
     proj_string = source_file["stereographic"].proj4
-    source_x = source_file["x"][:]
-    source_y = source_file["y"][:]
+    source_x = source_file["x_dim"][:]
+    source_y = source_file["y_dim"][:]
 
     element_lon = init_base.get_element_longitude()
     element_lat = init_base.get_element_latitude()
diff --git a/run/make_init_para24x30.py b/run/make_init_para24x30.py
index f7103b29d..c0aff27a8 100644
--- a/run/make_init_para24x30.py
+++ b/run/make_init_para24x30.py
@@ -25,16 +25,16 @@
     formatted.format = "%Y-%m-%dT%H:%M:%SZ"
     formatted[0] = "2010-01-01T00:00:00Z"
 
-    yDim = ncFile.createDimension("ydim", nfirst)
-    xDim = ncFile.createDimension("xdim", nsecond)
-    yVertexDim = ncFile.createDimension("yvertex", nfirst + 1)
-    xVertexDim = ncFile.createDimension("xvertex", nsecond + 1)
+    yDim = ncFile.createDimension("y_dim", nfirst)
+    xDim = ncFile.createDimension("x_dim", nsecond)
+    yVertexDim = ncFile.createDimension("y_vertex", nfirst + 1)
+    xVertexDim = ncFile.createDimension("x_vertex", nsecond + 1)
     dg_comp = ncFile.createDimension("dg_comp", n_dg)
     dgs_comp = ncFile.createDimension("dgstress_comp", n_dgstress)
     n_coords_comp = ncFile.createDimension("ncoords", n_coords)
-    field_dims = ("ydim", "xdim")
-    dg_dims = ("ydim", "xdim", "dg_comp")
-    coord_dims = ("yvertex", "xvertex", "ncoords")
+    field_dims = ("y_dim", "x_dim")
+    dg_dims = ("y_dim", "x_dim", "dg_comp")
+    coord_dims = ("y_vertex", "x_vertex", "ncoords")
 
     # Array coordinates
     space = 25000.0  # 25 km in metres
@@ -57,9 +57,9 @@
     x_cell = np.resize(x_cell1d, (nfirst, nsecond))
     y_cell = np.resize(y_cell1d, (nsecond, nfirst)).transpose()
 
-    elem_x = ncFile.createVariable("x", "f8", field_dims)
+    elem_x = ncFile.createVariable("x_dim", "f8", field_dims)
     elem_x[:, :] = x_cell
-    elem_y = ncFile.createVariable("y", "f8", field_dims)
+    elem_y = ncFile.createVariable("y_dim", "f8", field_dims)
     elem_y[:, :] = y_cell
 
     grid_azimuth = ncFile.createVariable("grid_azimuth", "f8", field_dims)
diff --git a/run/make_init_rect20x30.py b/run/make_init_rect20x30.py
index aa26e7f92..466f228ea 100644
--- a/run/make_init_rect20x30.py
+++ b/run/make_init_rect20x30.py
@@ -11,10 +11,10 @@
 
 ncFile.structure_name = "simple_rectangular"
 
-xDim = ncFile.createDimension("xdim", nx)
-yDim = ncFile.createDimension("ydim", ny)
+xDim = ncFile.createDimension("x_dim", nx)
+yDim = ncFile.createDimension("y_dim", ny)
 
-hfield_dims = ("ydim", "xdim")
+hfield_dims = ("y_dim", "x_dim")
 
 mask = ncFile.createVariable("mask", "f8", hfield_dims)
 # fmt: off
@@ -90,7 +90,7 @@
 hice[:, :] = cice[:, :] * 2
 hsnow = ncFile.createVariable("hsnow", "f8", hfield_dims)
 hsnow[:, :] = cice[:, :] / 2
-tsurf = ncFile.createVariable("tsurf", "f8", ("ydim", "xdim"))
+tsurf = ncFile.createVariable("tsurf", "f8", ("y_dim", "x_dim"))
 tsurf[:, :] = -0.5 - cice[:, :]
 
 mdi = -(2.0**300)
diff --git a/test/run_test_advection.sh b/test/run_test_advection.sh
index 8853f90fc..657e29dff 100755
--- a/test/run_test_advection.sh
+++ b/test/run_test_advection.sh
@@ -1,25 +1,30 @@
 #!/bin/sh
+# Prepares, executes and test the nextSIM-DG advection test
+
+# Log commands
+set -x
 
-# Get the name of the python executable
-if [ $# -lt 1 ]
-then
-    PYTHON=python
+# Get the name of the Python executable
+if [ $# -lt 1 ]; then
+  PYTHON=python
 else
-    PYTHON=$1
+  PYTHON=$1
 fi
 
-
-# Prepares, executes and test the nextSIM-DG advection test
-
+# Set filenames
 restart_file=advection_test_init.nc
 final_file=advection_test.restart.nc
 diag_file=advection_test.diagnostic.nc
 
-$PYTHON make_advection_test_data.py
-echo run_advection_test.sh
+# Setup
+${PYTHON} make_advection_test_data.py
+
+# Run the advection test and check it corrects the expected files
 ./run_advection_test.sh
 $PYTHON test_advection.py
 test_return_value=$?
+
+# Cleanup
 rm $restart_file $out_file $diag_file
 rm nextsim.*.log
 
diff --git a/test/run_test_error_handling_test.sh b/test/run_test_error_handling_test.sh
index b7a91ef10..cca529471 100755
--- a/test/run_test_error_handling_test.sh
+++ b/test/run_test_error_handling_test.sh
@@ -1,26 +1,27 @@
 #!/bin/sh
+# Prepares and executes nextSIM-DG integration test based on the January 2010
+# Arctic simulation, but with erroneous forcing data. Tests if the return code
+# is non-zero.
+
+# Log commands
+set -x
 
-# Get the name of the python executable
-if [ $# -lt 1 ]
-then
-    PYTHON=python
+# Get the name of the Python executable
+if [ $# -lt 1 ]; then
+  PYTHON=python
 else
-    PYTHON=$1
+  PYTHON=$1
 fi
 
+# Setup
+${PYTHON} make_init25kmNH_test_data.py
+${PYTHON} era5_topaz4_test_data.py --errors
 
-# Prepares and executes nextSIM-DG integration test based on the January 2010
-# Arctic simulation, but with erroneous forcing data. Tests if the return code
-# is non-zero.
-
-$PYTHON make_init25kmNH_test_data.py
-$PYTHON era5_topaz4_test_data.py --errors
-echo run_integration_test.sh
+# Run the integration test and check it raises the expected errors
 ./run_integration_test.sh
 
 # The test passes if the model returns a non-zero exit code
-if [ $? -ne 0 ]
-then
+if [ $? -ne 0 ]; then
   echo "Success: Error caught and model exits with a non-zero exit code"
   exit 0
 else
diff --git a/test/run_test_jan2010_integration.sh b/test/run_test_jan2010_integration.sh
index 88a11e7d1..beae5dc29 100755
--- a/test/run_test_jan2010_integration.sh
+++ b/test/run_test_jan2010_integration.sh
@@ -1,28 +1,33 @@
 #!/bin/sh
+# Prepares, executes and test the nextSIM-DG integration test based on the
+# January 2010 Arctic simulation
+
+# Log commands
+set -x
 
-# Get the name of the python executable
-if [ $# -lt 1 ]
-then
-    PYTHON=python
+# Get the name of the Python executable
+if [ $# -lt 1 ]; then
+  PYTHON=python
 else
-    PYTHON=$1
+  PYTHON=$1
 fi
 
-
-# Prepares, executes and test the nextSIM-DG integration test based on the
-# January 2010 Arctic simulation
-
+# Set filenames
 restart_file=init_25km_NH.nc
 era5_file=25km_NH.ERA5_2010-01-01T000000_test_data.nc
 topaz_file=25km_NH.TOPAZ4_2010-01-01T000000_test_data.nc
 out_file=out.integration_test.nc
 
-$PYTHON make_init25kmNH_test_data.py
-$PYTHON era5_topaz4_test_data.py
-echo run_integration_test.sh
+# Setup
+${PYTHON} make_init25kmNH_test_data.py
+${PYTHON} era5_topaz4_test_data.py
+
+# Run the integration test and check it corrects the expected files
 ./run_integration_test.sh
-$PYTHON test_integration.py
+${PYTHON} test_integration.py
 test_return_value=$?
+
+# Cleanup
 rm $restart_file $era5_file $topaz_file $out_file
 rm nextsim.*.log