From 7cea3c2717a426f70df4d8005adf2130fcef049b Mon Sep 17 00:00:00 2001
From: Joe Wallwork <jw2423@cam.ac.uk>
Date: Tue, 24 Mar 2026 15:12:59 +0000
Subject: [PATCH 1/9] Give halos in XIOS case

---
 core/src/ModelMetadata.cpp | 7 +------
 1 file changed, 1 insertion(+), 6 deletions(-)

diff --git a/core/src/ModelMetadata.cpp b/core/src/ModelMetadata.cpp
index 648c86423..7ace2952a 100644
--- a/core/src/ModelMetadata.cpp
+++ b/core/src/ModelMetadata.cpp
@@ -15,9 +15,8 @@
 #include "include/ParallelNetcdfFile.hpp"
 #ifdef USE_XIOS
 #include "include/Xios.hpp"
-#else
-#include "include/Halo.hpp"
 #endif
+#include "include/Halo.hpp"
 #endif
 #include "include/gridNames.hpp"
 #include <array>
@@ -311,12 +310,8 @@ void ModelMetadata::setDimensionsFromFile(const std::string& filename)
                 localLength = dim.getSize();
                 start = 0;
             }
-#if USE_XIOS
-            ModelArray::setDimension(dimType, dim.getSize(), localLength, start);
-#else
             ModelArray::setDimension(
                 dimType, dim.getSize(), localLength + 2 * Halo::haloWidth, start);
-#endif
 #else
             ModelArray::setDimension(dimType, dim.getSize());
 #endif

From 43499b77281c2a23cb3ffa776291c4a9ae81de30 Mon Sep 17 00:00:00 2001
From: Joe Wallwork <jw2423@cam.ac.uk>
Date: Tue, 24 Mar 2026 15:42:09 +0000
Subject: [PATCH 2/9] Attempt at only reading inner block

---
 core/src/Xios.cpp | 54 ++++++++++++++++++++++++++++++++++++-----------
 1 file changed, 42 insertions(+), 12 deletions(-)

diff --git a/core/src/Xios.cpp b/core/src/Xios.cpp
index 8a7279c0d..3e81672c8 100644
--- a/core/src/Xios.cpp
+++ b/core/src/Xios.cpp
@@ -19,6 +19,7 @@
 
 #include "StructureModule/include/ParametricGrid.hpp"
 #include "include/Finalizer.hpp"
+#include "include/Halo.hpp"
 #include "include/ModelMPI.hpp"
 #include "include/ModelMetadata.hpp"
 #include "include/ParallelNetcdfFile.hpp"
@@ -1681,6 +1682,8 @@ void Xios::read(const std::string& fieldId, ModelArray& modelarray)
     const ModelArray::Type& type = modelarray.getType();
     const ModelArray::Type& expectedType = getFieldType(fieldId);
 
+    Halo halo(modelarray);
+
     // Account for fields to be read in as HField but converted to DGField
     if (inputFieldsToConvert.count(fieldId)) {
         if (expectedType != ModelArray::Type::H) {
@@ -1692,9 +1695,11 @@ void Xios::read(const std::string& fieldId, ModelArray& modelarray)
                 "Xios::read: field '" + fieldId + "' was expected to be converted to a DGField");
         }
         HField inputarray(ModelArray::Type::H);
+        ModelArray::DataType& inputData = inputarray.getDataRef();
+        halo.getInnerBlock(modelarray.getDataRef(), inputData);
+        inputData.resize(halo.getInnerSize(), modelarray.nComponents());
         auto& dims = inputarray.dimensions();
-        cxios_read_data_k82(
-            fieldId.c_str(), fieldId.length(), inputarray.getData(), dims[0], dims[1]);
+        cxios_read_data_k82(fieldId.c_str(), fieldId.length(), inputData.data(), dims[0], dims[1]);
         modelarray = 0;
         // FIXME: Conversion with overloaded '=' operator is known to be problematic
         modelarray = inputarray;
@@ -1702,23 +1707,48 @@ void Xios::read(const std::string& fieldId, ModelArray& modelarray)
     }
 
     // Other field types should not need converting
-    auto& dims = modelarray.dimensions();
     if (type != expectedType) {
         throw std::runtime_error(
             "Xios::read: field '" + fieldId + "' does not have the expected type");
     }
-    if ((type == ModelArray::Type::H) || (type == ModelArray::Type::CG)) {
-        cxios_read_data_k82(
-            fieldId.c_str(), fieldId.length(), modelarray.getData(), dims[0], dims[1]);
+    if (type == ModelArray::Type::H) {
+        HField inputarray(ModelArray::Type::H);
+        ModelArray::DataType& inputData = inputarray.getDataRef();
+        halo.getInnerBlock(modelarray.getDataRef(), inputData);
+        inputData.resize(halo.getInnerSize(), modelarray.nComponents());
+        auto& dims = inputarray.dimensions();
+        cxios_read_data_k82(fieldId.c_str(), fieldId.length(), inputData.data(), dims[0], dims[1]);
+    } else if (type == ModelArray::Type::CG) {
+        CGField inputarray(ModelArray::Type::CG);
+        ModelArray::DataType& inputData = inputarray.getDataRef();
+        halo.getInnerBlock(modelarray.getDataRef(), inputData);
+        inputData.resize(halo.getInnerSize(), modelarray.nComponents());
+        auto& dims = inputarray.dimensions();
+        cxios_read_data_k82(fieldId.c_str(), fieldId.length(), inputData.data(), dims[0], dims[1]);
     } else if (type == ModelArray::Type::VERTEX) {
-        cxios_read_data_k83(fieldId.c_str(), fieldId.length(), modelarray.getData(), dims[0],
-            dims[1], ModelArray::size(ModelArray::Dimension::NCOORDS));
+        VertexField inputarray(ModelArray::Type::VERTEX);
+        ModelArray::DataType& inputData = inputarray.getDataRef();
+        halo.getInnerBlock(modelarray.getDataRef(), inputData);
+        inputData.resize(halo.getInnerSize(), modelarray.nComponents());
+        auto& dims = inputarray.dimensions();
+        cxios_read_data_k83(fieldId.c_str(), fieldId.length(), inputData.data(), dims[0], dims[1],
+            ModelArray::size(ModelArray::Dimension::NCOORDS));
     } else if (type == ModelArray::Type::DG) {
-        cxios_read_data_k83(fieldId.c_str(), fieldId.length(), modelarray.getData(), dims[0],
-            dims[1], ModelArray::size(ModelArray::Dimension::DG));
+        DGField inputarray(ModelArray::Type::DG);
+        ModelArray::DataType& inputData = inputarray.getDataRef();
+        halo.getInnerBlock(modelarray.getDataRef(), inputData);
+        inputData.resize(halo.getInnerSize(), modelarray.nComponents());
+        auto& dims = inputarray.dimensions();
+        cxios_read_data_k83(fieldId.c_str(), fieldId.length(), inputData.data(), dims[0], dims[1],
+            ModelArray::size(ModelArray::Dimension::DG));
     } else if (type == ModelArray::Type::DGSTRESS) {
-        cxios_read_data_k83(fieldId.c_str(), fieldId.length(), modelarray.getData(), dims[0],
-            dims[1], ModelArray::size(ModelArray::Dimension::DGSTRESS));
+        DGSField inputarray(ModelArray::Type::DGSTRESS);
+        ModelArray::DataType& inputData = inputarray.getDataRef();
+        halo.getInnerBlock(modelarray.getDataRef(), inputData);
+        inputData.resize(halo.getInnerSize(), modelarray.nComponents());
+        auto& dims = inputarray.dimensions();
+        cxios_read_data_k83(fieldId.c_str(), fieldId.length(), inputData.data(), dims[0], dims[1],
+            ModelArray::size(ModelArray::Dimension::DGSTRESS));
     } else {
         throw std::invalid_argument(
             "Only HFields, VertexFields, DGFields, DGSFields, and CGFields are supported");

From 470eb8a32bf5aa113c189ecd18338a3dac662494 Mon Sep 17 00:00:00 2001
From: Joe Wallwork <jw2423@cam.ac.uk>
Date: Tue, 24 Mar 2026 15:52:08 +0000
Subject: [PATCH 3/9] Define and use getInnerShape

---
 core/src/Xios.cpp         | 10 +++++-----
 core/src/include/Halo.hpp |  5 +++++
 2 files changed, 10 insertions(+), 5 deletions(-)

diff --git a/core/src/Xios.cpp b/core/src/Xios.cpp
index 3e81672c8..10762f213 100644
--- a/core/src/Xios.cpp
+++ b/core/src/Xios.cpp
@@ -1698,7 +1698,7 @@ void Xios::read(const std::string& fieldId, ModelArray& modelarray)
         ModelArray::DataType& inputData = inputarray.getDataRef();
         halo.getInnerBlock(modelarray.getDataRef(), inputData);
         inputData.resize(halo.getInnerSize(), modelarray.nComponents());
-        auto& dims = inputarray.dimensions();
+        auto& dims = halo.getInnerShape();
         cxios_read_data_k82(fieldId.c_str(), fieldId.length(), inputData.data(), dims[0], dims[1]);
         modelarray = 0;
         // FIXME: Conversion with overloaded '=' operator is known to be problematic
@@ -1716,21 +1716,21 @@ void Xios::read(const std::string& fieldId, ModelArray& modelarray)
         ModelArray::DataType& inputData = inputarray.getDataRef();
         halo.getInnerBlock(modelarray.getDataRef(), inputData);
         inputData.resize(halo.getInnerSize(), modelarray.nComponents());
-        auto& dims = inputarray.dimensions();
+        auto& dims = halo.getInnerShape();
         cxios_read_data_k82(fieldId.c_str(), fieldId.length(), inputData.data(), dims[0], dims[1]);
     } else if (type == ModelArray::Type::CG) {
         CGField inputarray(ModelArray::Type::CG);
         ModelArray::DataType& inputData = inputarray.getDataRef();
         halo.getInnerBlock(modelarray.getDataRef(), inputData);
         inputData.resize(halo.getInnerSize(), modelarray.nComponents());
-        auto& dims = inputarray.dimensions();
+        auto& dims = halo.getInnerShape();
         cxios_read_data_k82(fieldId.c_str(), fieldId.length(), inputData.data(), dims[0], dims[1]);
     } else if (type == ModelArray::Type::VERTEX) {
         VertexField inputarray(ModelArray::Type::VERTEX);
         ModelArray::DataType& inputData = inputarray.getDataRef();
         halo.getInnerBlock(modelarray.getDataRef(), inputData);
         inputData.resize(halo.getInnerSize(), modelarray.nComponents());
-        auto& dims = inputarray.dimensions();
+        auto& dims = halo.getInnerShape();
         cxios_read_data_k83(fieldId.c_str(), fieldId.length(), inputData.data(), dims[0], dims[1],
             ModelArray::size(ModelArray::Dimension::NCOORDS));
     } else if (type == ModelArray::Type::DG) {
@@ -1738,7 +1738,7 @@ void Xios::read(const std::string& fieldId, ModelArray& modelarray)
         ModelArray::DataType& inputData = inputarray.getDataRef();
         halo.getInnerBlock(modelarray.getDataRef(), inputData);
         inputData.resize(halo.getInnerSize(), modelarray.nComponents());
-        auto& dims = inputarray.dimensions();
+        auto& dims = halo.getInnerShape();
         cxios_read_data_k83(fieldId.c_str(), fieldId.length(), inputData.data(), dims[0], dims[1],
             ModelArray::size(ModelArray::Dimension::DG));
     } else if (type == ModelArray::Type::DGSTRESS) {
diff --git a/core/src/include/Halo.hpp b/core/src/include/Halo.hpp
index db24fa33e..923ccc97b 100644
--- a/core/src/include/Halo.hpp
+++ b/core/src/include/Halo.hpp
@@ -434,6 +434,11 @@ class Halo {
      */
     size_t getInnerSize() { return m_innerNx * m_innerNy; }
 
+    /*!
+     * @brief Returns shape of the inner array
+     */
+    const MultiDim getInnerShape() const { return std::vector({ m_innerNx, m_innerNy }); }
+
     /*!
      * @brief Get inner block from source array and copy into target
      *

From d1c40449c3aa002db3da2c83d00cf55c7958e22e Mon Sep 17 00:00:00 2001
From: Joe Wallwork <jw2423@cam.ac.uk>
Date: Tue, 24 Mar 2026 16:18:52 +0000
Subject: [PATCH 4/9] Use tempData approach

---
 core/src/Xios.cpp | 62 +++++++++++++++++++++--------------------------
 1 file changed, 28 insertions(+), 34 deletions(-)

diff --git a/core/src/Xios.cpp b/core/src/Xios.cpp
index 10762f213..f3b4fc8b4 100644
--- a/core/src/Xios.cpp
+++ b/core/src/Xios.cpp
@@ -1682,6 +1682,7 @@ void Xios::read(const std::string& fieldId, ModelArray& modelarray)
     const ModelArray::Type& type = modelarray.getType();
     const ModelArray::Type& expectedType = getFieldType(fieldId);
 
+    // Create a halo object based on the input array
     Halo halo(modelarray);
 
     // Account for fields to be read in as HField but converted to DGField
@@ -1694,12 +1695,21 @@ void Xios::read(const std::string& fieldId, ModelArray& modelarray)
             throw std::runtime_error(
                 "Xios::read: field '" + fieldId + "' was expected to be converted to a DGField");
         }
+
+        // Create a HField to read into
         HField inputarray(ModelArray::Type::H);
-        ModelArray::DataType& inputData = inputarray.getDataRef();
-        halo.getInnerBlock(modelarray.getDataRef(), inputData);
-        inputData.resize(halo.getInnerSize(), modelarray.nComponents());
+
+        // Read into a temporary array without the halo
+        ModelArray::DataType tempData;
+        tempData.resize(halo.getInnerSize(), inputarray.nComponents());
         auto& dims = halo.getInnerShape();
-        cxios_read_data_k82(fieldId.c_str(), fieldId.length(), inputData.data(), dims[0], dims[1]);
+        cxios_read_data_k82(fieldId.c_str(), fieldId.length(), tempData.data(), dims[0], dims[1]);
+
+        // Set the inner block of the array based on what was read and exchange halos
+        halo.setInnerBlock(tempData, inputarray.getDataRef());
+        halo.exchangeHalos(inputarray.getDataRef());
+
+        // Set the DGField based on the contents of the HField
         modelarray = 0;
         // FIXME: Conversion with overloaded '=' operator is known to be problematic
         modelarray = inputarray;
@@ -1711,48 +1721,32 @@ void Xios::read(const std::string& fieldId, ModelArray& modelarray)
         throw std::runtime_error(
             "Xios::read: field '" + fieldId + "' does not have the expected type");
     }
+
+    // Read into a temporary array without the halo
+    ModelArray::DataType tempData;
+    tempData.resize(halo.getInnerSize(), modelarray.nComponents());
+    auto& dims = halo.getInnerShape();
     if (type == ModelArray::Type::H) {
-        HField inputarray(ModelArray::Type::H);
-        ModelArray::DataType& inputData = inputarray.getDataRef();
-        halo.getInnerBlock(modelarray.getDataRef(), inputData);
-        inputData.resize(halo.getInnerSize(), modelarray.nComponents());
-        auto& dims = halo.getInnerShape();
-        cxios_read_data_k82(fieldId.c_str(), fieldId.length(), inputData.data(), dims[0], dims[1]);
+        cxios_read_data_k82(fieldId.c_str(), fieldId.length(), tempData.data(), dims[0], dims[1]);
     } else if (type == ModelArray::Type::CG) {
-        CGField inputarray(ModelArray::Type::CG);
-        ModelArray::DataType& inputData = inputarray.getDataRef();
-        halo.getInnerBlock(modelarray.getDataRef(), inputData);
-        inputData.resize(halo.getInnerSize(), modelarray.nComponents());
-        auto& dims = halo.getInnerShape();
-        cxios_read_data_k82(fieldId.c_str(), fieldId.length(), inputData.data(), dims[0], dims[1]);
+        cxios_read_data_k82(fieldId.c_str(), fieldId.length(), tempData.data(), dims[0], dims[1]);
     } else if (type == ModelArray::Type::VERTEX) {
-        VertexField inputarray(ModelArray::Type::VERTEX);
-        ModelArray::DataType& inputData = inputarray.getDataRef();
-        halo.getInnerBlock(modelarray.getDataRef(), inputData);
-        inputData.resize(halo.getInnerSize(), modelarray.nComponents());
-        auto& dims = halo.getInnerShape();
-        cxios_read_data_k83(fieldId.c_str(), fieldId.length(), inputData.data(), dims[0], dims[1],
+        cxios_read_data_k83(fieldId.c_str(), fieldId.length(), tempData.data(), dims[0], dims[1],
             ModelArray::size(ModelArray::Dimension::NCOORDS));
     } else if (type == ModelArray::Type::DG) {
-        DGField inputarray(ModelArray::Type::DG);
-        ModelArray::DataType& inputData = inputarray.getDataRef();
-        halo.getInnerBlock(modelarray.getDataRef(), inputData);
-        inputData.resize(halo.getInnerSize(), modelarray.nComponents());
-        auto& dims = halo.getInnerShape();
-        cxios_read_data_k83(fieldId.c_str(), fieldId.length(), inputData.data(), dims[0], dims[1],
+        cxios_read_data_k83(fieldId.c_str(), fieldId.length(), tempData.data(), dims[0], dims[1],
             ModelArray::size(ModelArray::Dimension::DG));
     } else if (type == ModelArray::Type::DGSTRESS) {
-        DGSField inputarray(ModelArray::Type::DGSTRESS);
-        ModelArray::DataType& inputData = inputarray.getDataRef();
-        halo.getInnerBlock(modelarray.getDataRef(), inputData);
-        inputData.resize(halo.getInnerSize(), modelarray.nComponents());
-        auto& dims = inputarray.dimensions();
-        cxios_read_data_k83(fieldId.c_str(), fieldId.length(), inputData.data(), dims[0], dims[1],
+        cxios_read_data_k83(fieldId.c_str(), fieldId.length(), tempData.data(), dims[0], dims[1],
             ModelArray::size(ModelArray::Dimension::DGSTRESS));
     } else {
         throw std::invalid_argument(
             "Only HFields, VertexFields, DGFields, DGSFields, and CGFields are supported");
     }
+
+    // Set the inner block of the array based on what was read and exchange halos
+    halo.setInnerBlock(tempData, modelarray.getDataRef());
+    halo.exchangeHalos(modelarray.getDataRef());
 }
 }
 

From a87f0ae37dc54324003338d7f211a9b31b42e9c2 Mon Sep 17 00:00:00 2001
From: Joe Wallwork <jw2423@cam.ac.uk>
Date: Tue, 24 Mar 2026 16:42:13 +0000
Subject: [PATCH 5/9] Allow Halos based on constant ModelArrays

---
 core/src/include/Halo.hpp | 12 ++++++++++++
 1 file changed, 12 insertions(+)

diff --git a/core/src/include/Halo.hpp b/core/src/include/Halo.hpp
index 923ccc97b..de7a34afa 100644
--- a/core/src/include/Halo.hpp
+++ b/core/src/include/Halo.hpp
@@ -53,6 +53,18 @@ class Halo {
         initializeHaloMetadata();
     }
 
+    /*!
+     * @brief Constructs a halo object from constant ModelArray
+     * @param ma ModelArray object to create halo from
+     */
+    Halo(const ModelArray& ma)
+    {
+        m_numComps = ma.nComponents();
+        isVertex = ma.getType() == ModelArray::Type::VERTEX;
+        setSpatialDims();
+        initializeHaloMetadata();
+    }
+
     /*!
      * @brief Constructs a halo object from DGVectorHolder
      * @param dgvh DGVectorHolder object to create halo from

From 3ee279857d1367caef0a02a8762dd64591432b94 Mon Sep 17 00:00:00 2001
From: Joe Wallwork <jw2423@cam.ac.uk>
Date: Tue, 24 Mar 2026 16:42:27 +0000
Subject: [PATCH 6/9] Handle halos when writing with XIOS

---
 core/src/Xios.cpp | 24 ++++++++++++++++--------
 1 file changed, 16 insertions(+), 8 deletions(-)

diff --git a/core/src/Xios.cpp b/core/src/Xios.cpp
index f3b4fc8b4..e32f3f715 100644
--- a/core/src/Xios.cpp
+++ b/core/src/Xios.cpp
@@ -1634,7 +1634,6 @@ void Xios::write(const std::string& fieldId, const ModelArray& modelarray)
     if (modelarray.nDimensions() != 2) {
         throw std::invalid_argument("Only ModelArrays of dimension 2 are supported");
     }
-    auto& dims = modelarray.dimensions();
     const ModelArray::Type& type = modelarray.getType();
     domainWritten[domainIds[type]] = true;
 
@@ -1645,19 +1644,28 @@ void Xios::write(const std::string& fieldId, const ModelArray& modelarray)
             "Xios::write: field '" + fieldId + "' does not have the expected type");
     }
 
+    // Create a halo object based on the input array
+    Halo halo(modelarray);
+
+    // Copy data into a temporary array without the halo
+    ModelArray::DataType tempData;
+    tempData.resize(halo.getInnerSize(), modelarray.nComponents());
+    halo.getInnerBlock(modelarray.data(), tempData);
+    auto& dims = halo.getInnerShape();
+
     // Write out according to field type
     if ((type == ModelArray::Type::H) || (type == ModelArray::Type::CG)) {
         cxios_write_data_k82(
-            fieldId.c_str(), fieldId.length(), modelarray.getData(), dims[0], dims[1], -1);
+            fieldId.c_str(), fieldId.length(), tempData.data(), dims[0], dims[1], -1);
     } else if (type == ModelArray::Type::VERTEX) {
-        cxios_write_data_k83(fieldId.c_str(), fieldId.length(), modelarray.getData(), dims[0],
-            dims[1], ModelArray::size(ModelArray::Dimension::NCOORDS), -1);
+        cxios_write_data_k83(fieldId.c_str(), fieldId.length(), tempData.data(), dims[0], dims[1],
+            ModelArray::size(ModelArray::Dimension::NCOORDS), -1);
     } else if (type == ModelArray::Type::DG) {
-        cxios_write_data_k83(fieldId.c_str(), fieldId.length(), modelarray.getData(), dims[0],
-            dims[1], ModelArray::size(ModelArray::Dimension::DG), -1);
+        cxios_write_data_k83(fieldId.c_str(), fieldId.length(), tempData.data(), dims[0], dims[1],
+            ModelArray::size(ModelArray::Dimension::DG), -1);
     } else if (type == ModelArray::Type::DGSTRESS) {
-        cxios_write_data_k83(fieldId.c_str(), fieldId.length(), modelarray.getData(), dims[0],
-            dims[1], ModelArray::size(ModelArray::Dimension::DGSTRESS), -1);
+        cxios_write_data_k83(fieldId.c_str(), fieldId.length(), tempData.data(), dims[0], dims[1],
+            ModelArray::size(ModelArray::Dimension::DGSTRESS), -1);
     } else {
         throw std::invalid_argument(
             "Only HFields, VertexFields, DGFields, DGSFields, and CGFields are supported");

From bf8f9dd476c224308dd285b705f306ad11a0f1cb Mon Sep 17 00:00:00 2001
From: Joe Wallwork <jw2423@cam.ac.uk>
Date: Thu, 26 Mar 2026 13:45:59 +0000
Subject: [PATCH 7/9] Implement getFieldLocalDoFs

---
 core/src/Xios.cpp                     | 18 ++++++++++++++++++
 core/src/include/Xios.hpp             |  1 +
 core/src/include/xios_c_interface.hpp |  6 ++++--
 3 files changed, 23 insertions(+), 2 deletions(-)

diff --git a/core/src/Xios.cpp b/core/src/Xios.cpp
index e32f3f715..b5f6f3e78 100644
--- a/core/src/Xios.cpp
+++ b/core/src/Xios.cpp
@@ -1123,6 +1123,24 @@ void Xios::setDiagnosticFieldType(const std::string& fieldId, const ModelArray::
     setFieldType(fieldId, fieldType, DIAGNOSTIC);
 }
 
+/*!
+ * Get the rank-local degrees of freedom for a given ModelArray.
+ *
+ * @param modelarray The ModelArray to get DoF information for
+ * @return MultiDim vector containing the number of DoFs in each direction
+ */
+const ModelArray::MultiDim Xios::getFieldLocalDoFs(const ModelArray& modelarray)
+{
+    const ModelArray::Type& type = modelarray.getType();
+    const std::string& domainId = domainIds[type];
+    xios::CDomain* domain = getDomain(domainId);
+    int ni;
+    cxios_get_domain_ni(domain, &ni);
+    int nj;
+    cxios_get_domain_nj(domain, &nj);
+    return std::vector<size_t> { (size_t)ni, (size_t)nj };
+}
+
 /*!
  * @brief   Do an initial read of input files to deduce field dimensions.
  *
diff --git a/core/src/include/Xios.hpp b/core/src/include/Xios.hpp
index 8854dd3ce..268218147 100644
--- a/core/src/include/Xios.hpp
+++ b/core/src/include/Xios.hpp
@@ -80,6 +80,7 @@ class Xios : public Configured<Xios> {
         = { sstName, sssName, mldName, uName, vName, sshName };
     void setPrognosticFieldType(const std::string& fieldId, const ModelArray::Type& type);
     void setDiagnosticFieldType(const std::string& fieldId, const ModelArray::Type& type);
+    const ModelArray::MultiDim getFieldLocalDoFs(const ModelArray& modelarray);
 
     enum {
         ENABLED_KEY,
diff --git a/core/src/include/xios_c_interface.hpp b/core/src/include/xios_c_interface.hpp
index b90ff70ab..fc3a020d8 100644
--- a/core/src/include/xios_c_interface.hpp
+++ b/core/src/include/xios_c_interface.hpp
@@ -83,6 +83,10 @@ void cxios_set_domain_ni(xios::CDomain* domain_hdl, int ni);
 void cxios_set_domain_nj(xios::CDomain* domain_hdl, int nj);
 void cxios_set_domain_ibegin(xios::CDomain* domain_hdl, int ibegin);
 void cxios_set_domain_jbegin(xios::CDomain* domain_hdl, int jbegin);
+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);
+int cxios_get_domain_ni(xios::CDomain* domain_hdl, int* ni);
+int cxios_get_domain_nj(xios::CDomain* domain_hdl, int* nj);
 bool cxios_is_defined_domain_type(xios::CDomain* domain_hdl);
 bool cxios_is_defined_domain_dim_i_name(xios::CDomain* domain_hdl);
 bool cxios_is_defined_domain_dim_j_name(xios::CDomain* domain_hdl);
@@ -90,8 +94,6 @@ 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);

From 2042c124d4e51a2efc5aaee6598de37ed79341f0 Mon Sep 17 00:00:00 2001
From: Joe Wallwork <jw2423@cam.ac.uk>
Date: Thu, 26 Mar 2026 13:53:32 +0000
Subject: [PATCH 8/9] Account for halos in forcing test

---
 core/test/CMakeLists.txt           |  2 ++
 core/test/XiosReadForcing_test.cpp | 26 ++++++++++++++++++++++----
 2 files changed, 24 insertions(+), 4 deletions(-)

diff --git a/core/test/CMakeLists.txt b/core/test/CMakeLists.txt
index 5cceb9ab8..f18263dc4 100644
--- a/core/test/CMakeLists.txt
+++ b/core/test/CMakeLists.txt
@@ -223,6 +223,8 @@ if(ENABLE_MPI)
                 "${MODEL_INCLUDE_DIR}"
                 "${XIOS_INCLUDE_LIST}"
                 "${ModulesRoot}/StructureModule"
+                "../../dynamics/src"
+                "../../dynamics/src/include"
         )
         target_link_libraries(testXiosReadForcing_MPI2 PRIVATE nextsimlib doctest::doctest)
 
diff --git a/core/test/XiosReadForcing_test.cpp b/core/test/XiosReadForcing_test.cpp
index 5b1ec02ea..59b3f4003 100644
--- a/core/test/XiosReadForcing_test.cpp
+++ b/core/test/XiosReadForcing_test.cpp
@@ -7,7 +7,9 @@
 #undef INFO
 
 #include "StructureModule/include/ParametricGrid.hpp"
+#include "cgVector.hpp"
 #include "include/Finalizer.hpp"
+#include "include/Halo.hpp"
 #include "include/Model.hpp"
 #include "include/ModelMPI.hpp"
 #include "include/NextsimModule.hpp"
@@ -83,8 +85,6 @@ MPI_TEST_CASE("TestXiosReadForcing", 2)
     REQUIRE(xiosHandler.getCalendarStep() == 0);
 
     // Deduce the local lengths of the two dimensions
-    const size_t nx = ModelArray::size(ModelArray::Dimension::X);
-    const size_t ny = ModelArray::size(ModelArray::Dimension::Y);
     REQUIRE(ModelArray::nComponents(ModelArray::Type::DG) == DGCOMP);
     REQUIRE(ModelArray::size(ModelArray::Dimension::DG) == DGCOMP);
 
@@ -101,9 +101,18 @@ MPI_TEST_CASE("TestXiosReadForcing", 2)
         // Read ERA5 forcings from file and check they take the expected values
         for (const auto& [fieldName, modelarray] :
             pio->readForcingTimeStatic(era5ForcingFieldNames, time, era5ForcingFilename).data) {
+            // Extract the inner block of the ModelArray
+            Halo halo(modelarray);
+            ModelArray::DataType tempData;
+            tempData.resize(halo.getInnerSize(), modelarray.nComponents());
+            halo.getInnerBlock(modelarray.data(), tempData);
+            // Check that the inner block has the expected values
+            auto& ndofs = xiosHandler.getFieldLocalDoFs(modelarray);
+            auto nx = ndofs[0];
+            auto ny = ndofs[1];
             for (size_t j = 0; j < ny; ++j) {
                 for (size_t i = 0; i < nx; ++i) {
-                    REQUIRE(modelarray(i, j) == doctest::Approx(0.1 * ts));
+                    REQUIRE(tempData(i, j) == doctest::Approx(0.1 * ts));
                 }
             }
         }
@@ -114,9 +123,18 @@ MPI_TEST_CASE("TestXiosReadForcing", 2)
             for (const auto& [fieldName, modelarray] :
                 pio->readForcingTimeStatic(topazForcingFieldNames, time, topazForcingFilename)
                     .data) {
+                // Extract the inner block of the ModelArray
+                Halo halo(modelarray);
+                ModelArray::DataType tempData;
+                tempData.resize(halo.getInnerSize(), modelarray.nComponents());
+                halo.getInnerBlock(modelarray.data(), tempData);
+                // Check that the inner block has the expected values
+                auto& ndofs = xiosHandler.getFieldLocalDoFs(modelarray);
+                auto nx = ndofs[0];
+                auto ny = ndofs[1];
                 for (size_t j = 0; j < ny; ++j) {
                     for (size_t i = 0; i < nx; ++i) {
-                        REQUIRE(modelarray(i, j) == doctest::Approx(ts + 1));
+                        REQUIRE(tempData(i, j) == doctest::Approx(ts + 1));
                     }
                 }
             }

From 8c56ff2334cc1be8dd67bb7067e295b61d573228 Mon Sep 17 00:00:00 2001
From: Joe Wallwork <jw2423@cam.ac.uk>
Date: Thu, 26 Mar 2026 16:25:03 +0000
Subject: [PATCH 9/9] Account for halos in diagnostics tests

---
 core/test/CMakeLists.txt               |  3 +++
 core/test/XiosReadDiagnostic_test.cpp  | 17 ++++++++++++-----
 core/test/XiosWriteDiagnostic_test.cpp |  6 ++++--
 3 files changed, 19 insertions(+), 7 deletions(-)

diff --git a/core/test/CMakeLists.txt b/core/test/CMakeLists.txt
index f18263dc4..4fce8d260 100644
--- a/core/test/CMakeLists.txt
+++ b/core/test/CMakeLists.txt
@@ -241,6 +241,9 @@ if(ENABLE_MPI)
                 "${MODEL_INCLUDE_DIR}"
                 "${XIOS_INCLUDE_LIST}"
                 "${ModulesRoot}/StructureModule"
+                "${ModulesRoot}/StructureModule"
+                "../../dynamics/src"
+                "../../dynamics/src/include"
         )
         target_link_libraries(testXiosReadDiagnostic_MPI3 PRIVATE nextsimlib doctest::doctest)
 
diff --git a/core/test/XiosReadDiagnostic_test.cpp b/core/test/XiosReadDiagnostic_test.cpp
index 7ab24a1a4..72f3fcba9 100644
--- a/core/test/XiosReadDiagnostic_test.cpp
+++ b/core/test/XiosReadDiagnostic_test.cpp
@@ -6,7 +6,9 @@
 #include <doctest/extensions/doctest_mpi.h>
 #undef INFO
 
+#include "cgVector.hpp"
 #include "include/Finalizer.hpp"
+#include "include/Halo.hpp"
 #include "include/Model.hpp"
 #include "include/ModelMPI.hpp"
 #include "include/NextsimModule.hpp"
@@ -70,19 +72,24 @@ MPI_TEST_CASE("TestXiosReadDiagnostic", 3)
     Xios& xiosHandler = Xios::getInstance();
     REQUIRE(xiosHandler.getCalendarStep() == 0);
 
-    // Deduce the local lengths of the two dimensions
-    const size_t nx = ModelArray::size(ModelArray::Dimension::X);
-    const size_t ny = ModelArray::size(ModelArray::Dimension::Y);
-
     // Read restarts from file and check they take the expected values
     // NOTE: The ParametricGrid is created and the XIOS context definition is closed in the call to
     //       StructureFactory::stateFromFile()
     for (const auto [fieldName, modelarray] :
         StructureFactory::stateFromFile(diagnosticFilename).data) {
         REQUIRE(fieldName == hsnowName);
+        // Extract the inner block of the ModelArray
+        Halo halo(modelarray);
+        ModelArray::DataType tempData;
+        tempData.resize(halo.getInnerSize(), modelarray.nComponents());
+        halo.getInnerBlock(modelarray.data(), tempData);
+        // Check that the inner block has the expected values
+        auto& ndofs = xiosHandler.getFieldLocalDoFs(modelarray);
+        auto nx = ndofs[0];
+        auto ny = ndofs[1];
         for (size_t j = 0; j < ny; ++j) {
             for (size_t i = 0; i < nx; ++i) {
-                REQUIRE(modelarray(i, j) == doctest::Approx(0.15));
+                REQUIRE(tempData(i, j) == doctest::Approx(0.15));
             }
         }
     }
diff --git a/core/test/XiosWriteDiagnostic_test.cpp b/core/test/XiosWriteDiagnostic_test.cpp
index f937d2625..4a08b520b 100644
--- a/core/test/XiosWriteDiagnostic_test.cpp
+++ b/core/test/XiosWriteDiagnostic_test.cpp
@@ -64,8 +64,6 @@ MPI_TEST_CASE("TestXiosWriteDiagnostic", 2)
     Xios& xiosHandler = Xios::getInstance();
 
     // Set ModelArray dimensions
-    const size_t nx = ModelArray::size(ModelArray::Dimension::X);
-    const size_t ny = ModelArray::size(ModelArray::Dimension::Y);
     REQUIRE(ModelArray::nComponents(ModelArray::Type::DG) == DGCOMP);
     REQUIRE(ModelArray::size(ModelArray::Dimension::DG) == DGCOMP);
 
@@ -92,11 +90,15 @@ MPI_TEST_CASE("TestXiosWriteDiagnostic", 2)
     // Simulate 4 iterations (timesteps)
     ModelMetadata& metadata = ModelMetadata::getInstance();
     const Duration& timestep = metadata.stepLength();
+    const size_t nx = ModelArray::size(ModelArray::Dimension::X);
+    const size_t ny = ModelArray::size(ModelArray::Dimension::Y);
     int rank;
     MPI_Comm_rank(test_comm, &rank);
     for (int ts = 0; ts <= 4; ts++) {
 
         // Update diagnostics
+        // NOTE: This example writes the (spatially constant) value into the halo region as well as
+        //       the interior. It doesn't matter because the halo values get discarded.
         for (size_t j = 0; j < ny; ++j) {
             for (size_t i = 0; i < nx; ++i) {
                 hsnow(i, j) = 0.1 * ts;