diff --git a/VERSION b/VERSION
index 9084fa2..26aaba0 100644
--- a/VERSION
+++ b/VERSION
@@ -1 +1 @@
-1.1.0
+1.2.0
diff --git a/src/orca-jedi/CMakeLists.txt b/src/orca-jedi/CMakeLists.txt
index 582ec3f..d4aad4e 100644
--- a/src/orca-jedi/CMakeLists.txt
+++ b/src/orca-jedi/CMakeLists.txt
@@ -27,7 +27,7 @@ nemo_io/NemoFieldReader.h
 nemo_io/NemoFieldReader.cc
 nemo_io/NemoFieldWriter.h
 nemo_io/NemoFieldWriter.cc
-nemo_io/OrcaIndex.h
+nemo_io/AtlasIndex.h
 nemo_io/ReadServer.h
 nemo_io/ReadServer.cc
 nemo_io/WriteServer.h
diff --git a/src/orca-jedi/geometry/Geometry.cc b/src/orca-jedi/geometry/Geometry.cc
index a2d5d2f..59bb96a 100644
--- a/src/orca-jedi/geometry/Geometry.cc
+++ b/src/orca-jedi/geometry/Geometry.cc
@@ -46,12 +46,19 @@ Geometry::Geometry(const eckit::Configuration & config,
                    const eckit::mpi::Comm & comm) :
                       comm_(comm), vars_(orcaVariableFactory(config)),
                       n_levels_(config.getInt("number levels")),
-                      grid_(config.getString("grid name")),
                       eckit_timer_(new eckit::Timer("Geometry(ORCA): ", oops::Log::trace()))
 {
     eckit_timer_->start();
     log_status();
     params_.validateAndDeserialize(config);
+    {
+      eckit::PathName grid_spec_path(params_.gridName.value());
+      if (grid_spec_path.exists()) {
+        grid_ = atlas::Grid{atlas::Grid::Spec{grid_spec_path}};
+      } else {
+        grid_ = atlas::Grid{params_.gridName.value()};
+      }
+    }
     int64_t halo = params_.sourceMeshHalo.value();
     if ( ( (params_.partitioner.value() == "serial") || (comm.size() == 1) )
          && (halo > 0) ) {
diff --git a/src/orca-jedi/nemo_io/AtlasIndex.h b/src/orca-jedi/nemo_io/AtlasIndex.h
new file mode 100644
index 0000000..dbd6d9f
--- /dev/null
+++ b/src/orca-jedi/nemo_io/AtlasIndex.h
@@ -0,0 +1,220 @@
+/*
+ * (C) British Crown Copyright 2024 Met Office
+ */
+
+#pragma once
+
+#include <algorithm>
+#include <sstream>
+#include <limits>
+#include <map>
+#include <utility>
+#include <string>
+#include <memory>
+#include <vector>
+
+#include "eckit/exception/Exceptions.h"
+
+#include "oops/util/Logger.h"
+
+#include "atlas/mesh.h"
+#include "atlas/grid.h"
+#include "atlas/parallel/omp/omp.h"
+#include "atlas-orca/grid/OrcaGrid.h"
+
+namespace orcamodel {
+
+/// \brief Interface from netcdf i, j coordinates in a field to 1D index of atlas
+///        field data.
+class AtlasIndexToBufferIndex {
+ public:
+  virtual ~AtlasIndexToBufferIndex() {}
+  virtual int64_t operator()(const int i, const int j) const = 0;
+  virtual int64_t operator()(const size_t inode) const = 0;
+  virtual std::pair<int, int> ij(const size_t inode) const  = 0;
+  virtual size_t nx() const = 0;
+  virtual size_t ny() const = 0;
+};
+
+/// \brief Indexer from the orca netcdf i, j field to an index of a 1D buffer
+///        of orca data.
+class OrcaIndexToBufferIndex : public AtlasIndexToBufferIndex {
+ private:
+  atlas::OrcaGrid orcaGrid_;
+  int32_t ix_glb_max;
+  int32_t ix_glb_min;
+  int32_t iy_glb_max;
+  int32_t iy_glb_min;
+  int32_t glbarray_offset;
+  int32_t glbarray_jstride;
+  size_t nx_;
+  size_t ny_;
+  const atlas::Mesh mesh_;
+
+ public:
+  static std::string name() {return "ORCA";}
+  size_t nx() const {return nx_;}
+  size_t ny() const {return ny_;}
+
+  explicit OrcaIndexToBufferIndex(const atlas::Mesh& mesh) : orcaGrid_(mesh.grid()), mesh_(mesh) {
+    iy_glb_max = orcaGrid_.ny() + orcaGrid_.haloNorth() - 1;
+    ix_glb_max = orcaGrid_.nx() + orcaGrid_.haloEast() - 1;
+
+    nx_ = orcaGrid_.nx() + orcaGrid_.haloWest() + orcaGrid_.haloEast();
+    ny_ = orcaGrid_.ny() + orcaGrid_.haloSouth() + orcaGrid_.haloNorth();
+
+    // vector of local indices: necessary for remote indices of ghost nodes
+    iy_glb_min = -orcaGrid_.haloSouth();
+    ix_glb_min = -orcaGrid_.haloWest();
+    glbarray_offset  = -(nx_ * iy_glb_min) - ix_glb_min;
+    glbarray_jstride = nx_;
+  }
+
+  /// \brief Index of a 1D array corresponding to point i, j
+  /// \param i
+  /// \param j
+  /// \return index of a matching 1D array
+  int64_t operator()(const int i, const int j) const {
+      ATLAS_ASSERT(i <= ix_glb_max,
+          std::to_string(i) + " > " + std::to_string(ix_glb_max));
+      ATLAS_ASSERT(j <= iy_glb_max,
+          std::to_string(j) + " > " + std::to_string(iy_glb_max));
+      ATLAS_ASSERT(i >= ix_glb_min,
+          std::to_string(i) + " < " + std::to_string(ix_glb_min));
+      ATLAS_ASSERT(j >= iy_glb_min,
+          std::to_string(j) + " < " + std::to_string(iy_glb_min));
+      return glbarray_offset + j * glbarray_jstride + i;
+  }
+
+  /// \brief Index of a 1D array corresponding to a mesh node index
+  /// \param inode
+  /// \return index of a matching 1D array
+  int64_t operator()(const size_t inode) const {
+      auto ij = atlas::array::make_view<int32_t, 2>(mesh_.nodes().field("ij"));
+      return (*this)(ij(inode, 0), ij(inode, 1));
+  }
+
+  /// \brief i, j pair corresponding to a node number. Only use this for diagnostic purposes as
+  ///        it is not robust at the orca halo.
+  /// \param inode
+  /// \return std::pair of the 2D indices corresponding to the node.
+  std::pair<int, int> ij(const size_t inode) const {
+    auto ij = atlas::array::make_view<int32_t, 2>(mesh_.nodes().field("ij"));
+    int i = ij(inode, 0) >= 0 ? ij(inode, 0) : nx_ + ij(inode, 0);
+    const int ci = i >= nx_ ? i - nx_ : i;
+    int j = ij(inode, 1) + 1;
+    const int cj = j >= ny_ ? ny_ - 1 : j;
+    return std::pair<int, int>{ci, cj};
+  }
+};
+
+/// \brief Indexer from the regular lon lat structured grid netcdf i, j field to an index of a
+///        1D buffer of regular lon lat field data.
+class RegLonLatIndexToBufferIndex : public AtlasIndexToBufferIndex {
+ public:
+  static std::string name() {return "structured";}
+  size_t nx() const {return nx_;}
+  size_t ny() const {return ny_;}
+
+  explicit RegLonLatIndexToBufferIndex(const atlas::Mesh& mesh) : grid_(mesh.grid()) {
+    ATLAS_ASSERT(grid_.regular(), "RegLonLatIndexToBufferIndex only works with regular grids");
+    nx_ = grid_.nx(0);
+    ny_ = grid_.ny();
+
+    iy_glb_max = ny_ - 1;
+    ix_glb_max = nx_ - 1;
+
+    const size_t num_nodes = mesh.nodes().size();
+    inode2ij.resize(num_nodes);
+    auto gidx = atlas::array::make_view<atlas::gidx_t, 1>(mesh.nodes().global_index());
+    atlas_omp_parallel_for(size_t inode = 0; inode < num_nodes; ++inode) {
+      const int global_index = gidx(inode) - 1;
+      const int i = global_index % nx_;
+      const int j = std::floor(global_index / nx_);
+      inode2ij[inode] = std::pair<int, int>(i, j);
+    }
+    ATLAS_ASSERT(inode2ij.size() <= nx_*ny_,
+        std::to_string(inode2ij.size()) + " > " + std::to_string(nx_*ny_));
+  }
+
+  /// \brief Index of a 1D array corresponding to point i, j
+  /// \param i
+  /// \param j
+  /// \return index of a matching 1D array
+  int64_t operator()(const int i, const int j) const {
+    ATLAS_ASSERT(i >= 0, std::to_string(i) + " < 0");
+    ATLAS_ASSERT(i <= ix_glb_max,
+        std::to_string(i) + " > " + std::to_string(ix_glb_max));
+    ATLAS_ASSERT(j >= 0, std::to_string(j) + " < 0");
+    ATLAS_ASSERT(j <= iy_glb_max,
+        std::to_string(j) + " > " + std::to_string(iy_glb_max));
+    return j * nx_ + i;
+  }
+  /// \brief Index of a 1D array corresponding to a mesh node index
+  /// \param inode
+  /// \return index of a matching 1D array
+  int64_t operator()(const size_t inode) const {
+    auto[i, j] = inode2ij[inode];
+    return (*this)(i, j);
+  }
+
+  /// \brief i, j pair corresponding to a node number.
+  /// \param inode
+  /// \return std::pair of the 2D indices corresponding to the node.
+  std::pair<int, int> ij(const size_t inode) const {
+    return inode2ij[inode];
+  }
+
+ private:
+  atlas::StructuredGrid grid_;
+  size_t nx_;
+  size_t ny_;
+  int32_t ix_glb_max;
+  int32_t iy_glb_max;
+  std::vector<std::pair<int, int>> inode2ij;
+};
+
+/// \brief Factory for creating AtlasIndexToBufferIndex objects.
+class AtlasIndexToBufferIndexCreator {
+ public:
+  static void register_type(std::string name, AtlasIndexToBufferIndexCreator *factory) {
+    get_factory()[name] = factory;
+  }
+  virtual std::unique_ptr<AtlasIndexToBufferIndex> create_unique(const atlas::Mesh& mesh) = 0;
+  static std::unique_ptr<AtlasIndexToBufferIndex> create_unique(std::string name,
+      const atlas::Mesh& mesh) {
+    std::unique_ptr<AtlasIndexToBufferIndex> AtlasIndexToBufferIndex =
+      std::move(get_factory()[name]->create_unique(mesh));
+    return AtlasIndexToBufferIndex;
+  }
+  static std::map<std::string, AtlasIndexToBufferIndexCreator*> &get_factory() {
+    static std::map<std::string, AtlasIndexToBufferIndexCreator*> creator_map;
+    return creator_map;
+  }
+};
+
+/// \brief Factory for creating OrcaIndexToBufferIndex objects.
+class OrcaIndexToBufferIndexCreator : public AtlasIndexToBufferIndexCreator {
+ public:
+  OrcaIndexToBufferIndexCreator()
+    { AtlasIndexToBufferIndexCreator::register_type(OrcaIndexToBufferIndex::name(), this); }
+  std::unique_ptr<AtlasIndexToBufferIndex> create_unique(const atlas::Mesh& mesh) {
+    std::unique_ptr<AtlasIndexToBufferIndex> o2b_index(new OrcaIndexToBufferIndex(mesh));
+    return o2b_index;
+  }
+};
+static OrcaIndexToBufferIndexCreator orcaIndexToBufferIndexCreator;
+
+/// \brief Factory for creating RegLonLatIndexToBufferIndex objects.
+class RegLonLatIndexToBufferIndexCreator : public AtlasIndexToBufferIndexCreator {
+ public:
+  RegLonLatIndexToBufferIndexCreator()
+    { AtlasIndexToBufferIndexCreator::register_type(RegLonLatIndexToBufferIndex::name(), this); }
+  std::unique_ptr<AtlasIndexToBufferIndex> create_unique(const atlas::Mesh& mesh) {
+    std::unique_ptr<AtlasIndexToBufferIndex> r2b_index(new RegLonLatIndexToBufferIndex(mesh));
+    return r2b_index;
+  }
+};
+static RegLonLatIndexToBufferIndexCreator regLonLatIndexToBufferIndexCreator;
+
+}  // namespace orcamodel
diff --git a/src/orca-jedi/nemo_io/NemoFieldReader.cc b/src/orca-jedi/nemo_io/NemoFieldReader.cc
index 58d4a60..3249128 100644
--- a/src/orca-jedi/nemo_io/NemoFieldReader.cc
+++ b/src/orca-jedi/nemo_io/NemoFieldReader.cc
@@ -20,8 +20,7 @@
 #include "oops/util/Logger.h"
 #include "oops/util/Duration.h"
 
-#include "atlas-orca/grid/OrcaGrid.h"
-#include "orca-jedi/nemo_io/OrcaIndex.h"
+#include "orca-jedi/nemo_io/AtlasIndex.h"
 
 namespace orcamodel {
 
diff --git a/src/orca-jedi/nemo_io/OrcaIndex.h b/src/orca-jedi/nemo_io/OrcaIndex.h
deleted file mode 100644
index 96719d1..0000000
--- a/src/orca-jedi/nemo_io/OrcaIndex.h
+++ /dev/null
@@ -1,67 +0,0 @@
-/*
- * (C) British Crown Copyright 2024 Met Office
- */
-
-#pragma once
-
-#include <algorithm>
-#include <sstream>
-#include <limits>
-#include <map>
-#include <utility>
-#include <string>
-#include <memory>
-
-#include "eckit/exception/Exceptions.h"
-
-#include "oops/util/Logger.h"
-
-#include "atlas/mesh.h"
-#include "atlas-orca/grid/OrcaGrid.h"
-
-namespace orcamodel {
-
-/// \brief Indexer from the orca netcdf i, j field to an index of a 1D buffer
-///        of orca data.
-struct OrcaIndexToBufferIndex {
- private:
-  atlas::OrcaGrid orcaGrid_;
-  int32_t ix_glb_max;
-  int32_t iy_glb_max;
-  int32_t glbarray_offset;
-  int32_t glbarray_jstride;
-  size_t nx_;
-  size_t ny_;
-
- public:
-  size_t nx() const {return nx_;}
-  size_t ny() const {return ny_;}
-
-  explicit OrcaIndexToBufferIndex(const atlas::Mesh& mesh) : orcaGrid_(mesh.grid()) {
-    iy_glb_max = orcaGrid_.ny() + orcaGrid_.haloNorth() - 1;
-    ix_glb_max = orcaGrid_.nx() + orcaGrid_.haloEast() - 1;
-
-    nx_ = orcaGrid_.nx() + orcaGrid_.haloWest() + orcaGrid_.haloEast();
-    ny_ = orcaGrid_.ny() + orcaGrid_.haloSouth() + orcaGrid_.haloNorth();
-
-    // vector of local indices: necessary for remote indices of ghost nodes
-    int iy_glb_min = -orcaGrid_.haloSouth();
-    int ix_glb_min = -orcaGrid_.haloWest();
-    glbarray_offset  = -(nx_ * iy_glb_min) - ix_glb_min;
-    glbarray_jstride = nx_;
-  }
-
-  /// \brief Index of a 1D array corresponding to point i, j
-  /// \param i
-  /// \param j
-  /// \return index of a matching 1D array
-  int64_t operator()(const int i, const int j) const {
-      ATLAS_ASSERT(i <= ix_glb_max,
-          std::to_string(i) + " > " + std::to_string(ix_glb_max));
-      ATLAS_ASSERT(j <= iy_glb_max,
-          std::to_string(j) + " > " + std::to_string(iy_glb_max));
-      return glbarray_offset + j * glbarray_jstride + i;
-  }
-};
-
-}  // namespace orcamodel
diff --git a/src/orca-jedi/nemo_io/ReadServer.cc b/src/orca-jedi/nemo_io/ReadServer.cc
index d901793..552598a 100644
--- a/src/orca-jedi/nemo_io/ReadServer.cc
+++ b/src/orca-jedi/nemo_io/ReadServer.cc
@@ -9,6 +9,7 @@
 #include <string>
 #include <vector>
 #include <memory>
+#include <utility>
 
 #include "oops/util/Logger.h"
 #include "atlas/parallel/omp/omp.h"
@@ -19,8 +20,10 @@ namespace orcamodel {
 ReadServer::ReadServer(std::shared_ptr<eckit::Timer> eckit_timer,
   const eckit::PathName& file_path, const atlas::Mesh& mesh) :
   mesh_(mesh),
-  orca_buffer_indices_(mesh),
   eckit_timer_(eckit_timer) {
+  buffer_indices_ = std::move(AtlasIndexToBufferIndexCreator::create_unique(
+          mesh.grid().type(), mesh));
+
   if (myrank == mpiroot) {
     reader_ = std::make_unique<NemoFieldReader>(file_path);
   }
@@ -37,7 +40,7 @@ template<class T> void ReadServer::read_var_on_root(const std::string& var_name,
               std::vector<T>& buffer) const {
   oops::Log::trace() << "State(ORCA)::nemo_io::ReadServer::read_var_on_root "
     << var_name << std::endl;
-  size_t size = orca_buffer_indices_.nx() * orca_buffer_indices_.ny();
+  size_t size = buffer_indices_->nx() * buffer_indices_->ny();
   if (myrank == mpiroot) {
     buffer = reader_->read_var_slice<T>(var_name, t_index, z_index);
   } else {
@@ -60,7 +63,7 @@ template void ReadServer::read_var_on_root<float>(const std::string& var_name,
 template<class T> void ReadServer::read_vertical_var_on_root(const std::string& var_name,
               const size_t n_levels,
               std::vector<T>& buffer) const {
-  size_t size = orca_buffer_indices_.nx() * orca_buffer_indices_.ny();
+  size_t size = buffer_indices_->nx() * buffer_indices_->ny();
   if (myrank == mpiroot) {
     buffer = reader_->read_vertical_var<T>(var_name, n_levels);
   } else {
@@ -83,13 +86,12 @@ template<class T> void ReadServer::fill_field(const std::vector<T>& buffer,
       atlas::array::ArrayView<T, 2>& field_view) const {
     oops::Log::trace() << "State(ORCA)::nemo_io::ReadServer::fill_field" << std::endl;
     auto ghost = atlas::array::make_view<int32_t, 1>(this->mesh_.nodes().ghost());
-    auto ij = atlas::array::make_view<int32_t, 2>(this->mesh_.nodes().field("ij"));
     const size_t num_nodes = field_view.shape(0);
     // "ReadServer buffer size does not equal the number of horizontal nodes in the field_view"
     ASSERT(num_nodes <= buffer.size());
     atlas_omp_parallel_for(size_t inode = 0; inode < num_nodes; ++inode) {
       if (ghost(inode)) continue;
-      const int64_t ibuf = orca_buffer_indices_(ij(inode, 0), ij(inode, 1));
+      const int64_t ibuf = (*buffer_indices_)(inode);
       field_view(inode, z_index) = buffer[ibuf];
     }
 }
@@ -146,7 +148,7 @@ template<class T> void ReadServer::read_var(const std::string& var_name,
     << var_name << std::endl;
 
   size_t n_levels = field_view.shape(1);
-  size_t size = orca_buffer_indices_.nx() * orca_buffer_indices_.ny();
+  size_t size = buffer_indices_->nx() * buffer_indices_->ny();
 
   std::vector<T> buffer;
   // For each level
diff --git a/src/orca-jedi/nemo_io/ReadServer.h b/src/orca-jedi/nemo_io/ReadServer.h
index d2044c2..84ba25e 100644
--- a/src/orca-jedi/nemo_io/ReadServer.h
+++ b/src/orca-jedi/nemo_io/ReadServer.h
@@ -15,7 +15,7 @@
 #include "atlas/grid.h"
 #include "atlas-orca/grid/OrcaGrid.h"
 
-#include "orca-jedi/nemo_io/OrcaIndex.h"
+#include "orca-jedi/nemo_io/AtlasIndex.h"
 #include "orca-jedi/nemo_io/NemoFieldReader.h"
 
 #include "eckit/exception/Exceptions.h"
@@ -61,7 +61,7 @@ void log_status() const {
   const size_t mpiroot = 0;
   const size_t myrank = atlas::mpi::rank();
   const atlas::Mesh& mesh_;
-  const OrcaIndexToBufferIndex orca_buffer_indices_;
+  std::unique_ptr<AtlasIndexToBufferIndex> buffer_indices_;
   std::unique_ptr<NemoFieldReader> reader_;
   std::shared_ptr<eckit::Timer> eckit_timer_;
 };
diff --git a/src/orca-jedi/nemo_io/WriteServer.cc b/src/orca-jedi/nemo_io/WriteServer.cc
index 47c5d02..957c123 100644
--- a/src/orca-jedi/nemo_io/WriteServer.cc
+++ b/src/orca-jedi/nemo_io/WriteServer.cc
@@ -5,6 +5,7 @@
 #include "orca-jedi/nemo_io/WriteServer.h"
 
 #include<algorithm>
+#include<utility>
 
 #include "atlas-orca/grid/OrcaGrid.h"
 
@@ -54,7 +55,7 @@ template<class T> std::vector<T> WriteServer::sort_buffer(
   const std::vector<T> & buffer) const {
   oops::Log::trace() << "State(ORCA)::nemo_io::WriteServer::sort_buffer "
                      << std::endl;
-  std::vector<T> sorted_buffer(orca_indices_.nx() * orca_indices_.ny());
+  std::vector<T> sorted_buffer(buffer_indices_->nx() * buffer_indices_->ny());
   for (size_t i_buf = 0; i_buf < buffer.size(); ++i_buf) {
     const size_t source_index = unsorted_buffer_indices_[i_buf];
     ASSERT(source_index < sorted_buffer.size());
@@ -79,7 +80,7 @@ template<class T> std::vector<T> WriteServer::gather_field_on_root(
                      << std::endl;
 
   std::vector<T> local_buffer;
-  size_t size = orca_indices_.nx() * orca_indices_.ny();
+  size_t size = buffer_indices_->nx() * buffer_indices_->ny();
   const bool has_mv = static_cast<bool>(missingValue);
 
   // handle the serial data case
@@ -136,7 +137,7 @@ template<class T> void WriteServer::write_vol_var(const std::string& var_name,
   oops::Log::trace() << "State(ORCA)::nemo_io::WriteServer::write_vol_var "
                      << var_name << std::endl;
   std::vector<T> buffer;
-  size_t size = orca_indices_.nx() * orca_indices_.ny();
+  size_t size = buffer_indices_->nx() * buffer_indices_->ny();
   if (myrank == mpiroot) {
     buffer.resize(n_levels_*size);
   }
@@ -197,7 +198,7 @@ void WriteServer::write_dimensions() {
   atlas::array::ArrayView<double, 2> lonlat{atlas::array::make_view<double, 2>(
                                               mesh_.nodes().lonlat())};
   atlas::field::MissingValue missingValue(mesh_.nodes().lonlat());
-  size_t size = orca_indices_.nx() * orca_indices_.ny();
+  size_t size = buffer_indices_->nx() * buffer_indices_->ny();
 
   const std::vector<double> lon_buffer = this->gather_field_on_root(lonlat, 0, missingValue);
   const std::vector<double> lat_buffer = this->gather_field_on_root(lonlat, 1, missingValue);
@@ -214,29 +215,28 @@ WriteServer::WriteServer(std::shared_ptr<eckit::Timer> eckit_timer,
     const atlas::Mesh& mesh,
     const std::vector<util::DateTime> datetimes,
     const std::vector<double> depths,
-    bool is_serial) : mesh_(mesh),
-  orca_indices_(mesh), eckit_timer_(eckit_timer), is_serial_(is_serial),
+    bool is_serial) : mesh_(mesh), eckit_timer_(eckit_timer), is_serial_(is_serial),
   n_levels_(depths.size()) {
   oops::Log::trace() << "State(ORCA)::nemo_io::WriteServer::WriteServer" << std::endl;
+  buffer_indices_ = std::move(AtlasIndexToBufferIndexCreator::create_unique(
+          mesh.grid().type(), mesh));
+
   std::vector<size_t> local_buf_indices;
 
   if (is_serial_) {
     if (myrank == mpiroot) {
       writer_ = std::make_unique<NemoFieldWriter>(file_path, datetimes,
-                                                  orca_indices_.nx(), orca_indices_.ny(),
+                                                  buffer_indices_->nx(), buffer_indices_->ny(),
                                                   depths);
     }
     this->write_dimensions();
-    recvcnt_ = orca_indices_.nx() * orca_indices_.ny();
+    recvcnt_ = buffer_indices_->nx() * buffer_indices_->ny();
     return;
   }
 
-  auto ij = atlas::array::make_view<int32_t, 2>(mesh_.nodes().field("ij"));
   auto orcaGrid = atlas::OrcaGrid(mesh_.grid());
   for (atlas::idx_t i_node = 0; i_node < mesh_.nodes().size(); ++i_node) {
-    const size_t i = ij(i_node, 0);
-    const size_t j = ij(i_node, 1);
-    local_buf_indices.emplace_back(orca_indices_(i, j));
+    local_buf_indices.emplace_back((*buffer_indices_)(i_node));
   }
 
   // gather all remote buffer indices
@@ -253,11 +253,11 @@ WriteServer::WriteServer(std::shared_ptr<eckit::Timer> eckit_timer,
       recvcnt_ += recvcounts_[jproc];
   }
 
-  ASSERT(recvcnt_ >= static_cast<int>(orca_indices_.nx()*orca_indices_.ny()));
+  ASSERT(recvcnt_ >= static_cast<int>(buffer_indices_->nx()*buffer_indices_->ny()));
 
   if (myrank == mpiroot) {
     writer_ = std::make_unique<NemoFieldWriter>(file_path, datetimes,
-                                                orca_indices_.nx(), orca_indices_.ny(),
+                                                buffer_indices_->nx(), buffer_indices_->ny(),
                                                 depths);
     unsorted_buffer_indices_.resize(recvcnt_);
   }
diff --git a/src/orca-jedi/nemo_io/WriteServer.h b/src/orca-jedi/nemo_io/WriteServer.h
index f739fea..2f63f15 100644
--- a/src/orca-jedi/nemo_io/WriteServer.h
+++ b/src/orca-jedi/nemo_io/WriteServer.h
@@ -15,7 +15,7 @@
 #include "atlas/field/MissingValue.h"
 #include "atlas-orca/grid/OrcaGrid.h"
 
-#include "orca-jedi/nemo_io/OrcaIndex.h"
+#include "orca-jedi/nemo_io/AtlasIndex.h"
 #include "orca-jedi/nemo_io/NemoFieldWriter.h"
 
 #include "eckit/exception/Exceptions.h"
@@ -63,7 +63,7 @@ class WriteServer {
   const size_t mpiroot = 0;
   const size_t myrank = atlas::mpi::rank();
   const atlas::Mesh& mesh_;
-  const OrcaIndexToBufferIndex orca_indices_;
+  std::unique_ptr<AtlasIndexToBufferIndex> buffer_indices_;
   std::vector<size_t> unsorted_buffer_indices_;
   std::vector<int> recvcounts_;
   std::vector<int> recvdispls_;
diff --git a/src/orca-jedi/state/State.cc b/src/orca-jedi/state/State.cc
index 3bed699..b7d53b8 100644
--- a/src/orca-jedi/state/State.cc
+++ b/src/orca-jedi/state/State.cc
@@ -125,6 +125,7 @@ State::~State() {
 }
 
 void State::subsetFieldSet(const oops::Variables & variables) {
+  oops::Log::trace() << "State(ORCA)::subsetFieldSet subsetting." << std::endl;
   atlas::FieldSet subset;
   for (int iVar = 0; iVar < variables.size(); iVar++) {
     auto variable = variables[iVar].name();
@@ -142,6 +143,7 @@ void State::subsetFieldSet(const oops::Variables & variables) {
     stateFields_.add(subset[variable]);
   }
   vars_ = variables;
+  oops::Log::trace() << "State(ORCA)::subsetFieldSet complete" << std::endl;
 }
 
 
diff --git a/src/tests/CMakeLists.txt b/src/tests/CMakeLists.txt
index 615505a..47da01a 100644
--- a/src/tests/CMakeLists.txt
+++ b/src/tests/CMakeLists.txt
@@ -34,6 +34,7 @@ ecbuild_add_test( TARGET test_orcamodel_hofx_ssh_parallel_checkerboard
                   ARGS testinput/hofx_nc_ssh_checkerboard.yaml
                   COMMAND orcamodel_hofx.x )
 
+# Disabled until halos are available from atlas-orca
 #ecbuild_add_test( TARGET test_orcamodel_hofx_ssh_parallel_eorca025
 #                  OMP 1
 #                  MPI 2
@@ -46,6 +47,18 @@ ecbuild_add_test( TARGET test_orcamodel_hofx_ssh_eorca025
                   ARGS testinput/hofx_nc_ssh_eorca025.yaml
                   COMMAND orcamodel_hofx.x )
 
+ecbuild_add_test( TARGET test_orcamodel_hofx_ssh_amm1
+                  OMP 1
+                  MPI 1
+                  ARGS testinput/hofx_nc_ssh_amm1.yaml
+                  COMMAND orcamodel_hofx.x )
+
+ecbuild_add_test( TARGET test_orcamodel_hofx3D_ssh_amm1
+                  OMP 1
+                  MPI 1
+                  ARGS testinput/hofx3d_nc_ssh_amm1.yaml
+                  COMMAND orcamodel_hofx3D.x )
+
 ecbuild_add_test( TARGET test_orcamodel_hofx_ssh_parallel_serial
                   OMP 1
                   MPI 2
diff --git a/src/tests/Data/CMakeLists.txt b/src/tests/Data/CMakeLists.txt
index 8c0e2a9..5cc20f9 100644
--- a/src/tests/Data/CMakeLists.txt
+++ b/src/tests/Data/CMakeLists.txt
@@ -20,6 +20,10 @@ hofx_two_vars_obs.nc
 hofx_potm_obs.nc
 hofx_prof_2var_obs.nc
 hofx_ssh_obs.nc
+amm1_coords.nc
+amm1_nemo.nc
+amm1_atlas_grid_spec.yaml
+amm_ssh_obs.nc
 )
 
 foreach(FILENAME ${orcajedi_test_data})
diff --git a/src/tests/Data/amm1_atlas_grid_spec.yaml b/src/tests/Data/amm1_atlas_grid_spec.yaml
new file mode 100644
index 0000000..8d05a7b
--- /dev/null
+++ b/src/tests/Data/amm1_atlas_grid_spec.yaml
@@ -0,0 +1,16 @@
+type : "regional"
+nx : 15
+ny : 19
+north : 67
+south : 38
+east : -5
+west : -20
+y_numbering: +1
+
+
+check :
+  size : 285
+  lonlat(first)     : [-20.0,38.0]
+  lonlat(last)      : [-5.0,67.0]
+  uid : 43ff9de7f0cbd13663eca41a22dc5dcb
+  bounding_box(n,w,s,e) : [67,-20,38,-5]
diff --git a/src/tests/Data/amm1_coords.nc b/src/tests/Data/amm1_coords.nc
new file mode 100644
index 0000000..01ed05a
--- /dev/null
+++ b/src/tests/Data/amm1_coords.nc
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:cc66af3234fa0396e90dc1842cf2f9103a82a37ccac07c23da8d3bb490f6718b
+size 9500
diff --git a/src/tests/Data/amm1_nemo.nc b/src/tests/Data/amm1_nemo.nc
new file mode 100644
index 0000000..9faceea
--- /dev/null
+++ b/src/tests/Data/amm1_nemo.nc
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:1ea2a156cde3a9737e11e0c3a69a58673ec41c7a0bffc14783dd35859c807e6f
+size 89301
diff --git a/src/tests/Data/amm_ssh_obs.nc b/src/tests/Data/amm_ssh_obs.nc
new file mode 100644
index 0000000..1955dda
--- /dev/null
+++ b/src/tests/Data/amm_ssh_obs.nc
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:d5239a49f946cf573daed8d2f38bcb584cee40adf6af5632536be372b72bb43e
+size 87319
diff --git a/src/tests/orca-jedi/CMakeLists.txt b/src/tests/orca-jedi/CMakeLists.txt
index 6396284..1d2c64f 100644
--- a/src/tests/orca-jedi/CMakeLists.txt
+++ b/src/tests/orca-jedi/CMakeLists.txt
@@ -13,7 +13,7 @@ ecbuild_add_test( TARGET  test_orcajedi_nemo_io_field_writer.x
                   CONDITION  eckit_HAVE_MPI
                   LIBS    orcamodel )
 
-ecbuild_add_test( TARGET  test_orcajedi_nemo_io_read_server.x
+ecbuild_add_test( TARGET  test_orcajedi_nemo_io_read_server_MPI2.x
                   SOURCES test_nemo_io_read_server.cc
                   ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT}
                   MPI        2
diff --git a/src/tests/orca-jedi/test_interpolator.cc b/src/tests/orca-jedi/test_interpolator.cc
index 3218deb..f5eeffe 100644
--- a/src/tests/orca-jedi/test_interpolator.cc
+++ b/src/tests/orca-jedi/test_interpolator.cc
@@ -4,6 +4,7 @@
 
 #include<sstream>
 #include<cmath>
+#include <map>
 
 #include "eckit/log/Bytes.h"
 #include "eckit/config/LocalConfiguration.h"
@@ -29,105 +30,184 @@ namespace test {
 
 const double ATOL = 1e-6;
 
+struct InterpTestSettingsFixture {
+ public:
+  eckit::LocalConfiguration geometry_config;
+  eckit::LocalConfiguration state_config;
+  eckit::LocalConfiguration interpolator_config;
+  size_t nlocs, nlevs;
+  std::vector<double> lons;
+  std::vector<double> lats;
+  oops::Variables surf_vars, vol_vars;
+  std::vector<double> surf_values;
+  std::vector<double> vol_values;
+};
+
+const double missing_value = util::missingValue<double>();
+
 //-----------------------------------------------------------------------------
 
-CASE("test basic interpolator") {
-  int nlevs = 3;
-  eckit::LocalConfiguration config;
-  config.set("grid name", "ORCA2_T");
-  config.set("number levels", nlevs);
-
-  std::vector<eckit::LocalConfiguration> nemo_var_mappings(4);
-  nemo_var_mappings[0].set("name", "sea_ice_area_fraction")
-    .set("nemo field name", "iiceconc")
-    .set("model space", "surface");
-  nemo_var_mappings[1].set("name", "sea_ice_area_fraction_error")
-    .set("nemo field name", "sic_tot_var")
-    .set("model space", "surface")
-    .set("variable type", "background variance");
-  nemo_var_mappings[2].set("name", "sea_surface_foundation_temperature")
-    .set("nemo field name", "votemper")
-    .set("model space", "surface");
-  nemo_var_mappings[3].set("name", "sea_water_potential_temperature")
-    .set("nemo field name", "votemper")
-    .set("model space", "volume");
-  config.set("nemo variables", nemo_var_mappings);
-  Geometry geometry(config, eckit::mpi::comm());
-
-  eckit::LocalConfiguration interp_conf;
-  interp_conf.set("type", "unstructured-bilinear-lonlat");
-  interp_conf.set("non_linear", "missing-if-all-missing-real32");
-  eckit::LocalConfiguration interpolator_conf;
-  interpolator_conf.set("atlas-interpolator", interp_conf);
-
-  OrcaInterpolatorParameters params;
-  params.validateAndDeserialize(interpolator_conf);
-
-  std::vector<double> lons({0, 120, 270});
-  std::vector<double> lats({88, 0, 30});
-  const int nlocs = 3;
-
-  // create a state from the test data
-  eckit::LocalConfiguration state_config;
-  std::vector<std::string> state_variables {
-      "sea_ice_area_fraction",
-      "sea_surface_foundation_temperature",
-      "sea_water_potential_temperature"};
-  state_config.set("state variables", state_variables);
-  state_config.set("date", "2021-06-30T00:00:00Z");
-  state_config.set("nemo field file", "../Data/orca2_t_nemo.nc");
-  state_config.set("variance field file", "../Data/orca2_t_bkg_var.nc");
-  OrcaStateParameters stateParams;
-  stateParams.validateAndDeserialize(state_config);
-  State state(geometry, stateParams);
-
-  SECTION("test interpolator succeeds even with no locations") {
-    Interpolator interpolator(interpolator_conf, geometry, {}, {});
+CASE("test  interpolator") {
+  std::map<std::string, InterpTestSettingsFixture> settings_map{
+      {"ORCA2_T", InterpTestSettingsFixture()}, {"AMM1", InterpTestSettingsFixture()}};
+
+  // ORCA2_T settings
+  {
+    settings_map["ORCA2_T"].nlocs = 3;
+    settings_map["ORCA2_T"].nlevs = 3;
+    settings_map["ORCA2_T"].geometry_config.set("grid name", "ORCA2_T");
+    settings_map["ORCA2_T"].geometry_config.set("number levels", settings_map["ORCA2_T"].nlevs);
+
+    std::vector<eckit::LocalConfiguration> nemo_var_mappings(4);
+    nemo_var_mappings[0].set("name", "sea_ice_area_fraction")
+      .set("nemo field name", "iiceconc")
+      .set("model space", "surface");
+    nemo_var_mappings[1].set("name", "sea_ice_area_fraction_error")
+      .set("nemo field name", "sic_tot_var")
+      .set("model space", "surface")
+      .set("variable type", "background variance");
+    nemo_var_mappings[2].set("name", "sea_surface_foundation_temperature")
+      .set("nemo field name", "votemper")
+      .set("model space", "surface");
+    nemo_var_mappings[3].set("name", "sea_water_potential_temperature")
+      .set("nemo field name", "votemper")
+      .set("model space", "volume");
+    settings_map["ORCA2_T"].geometry_config.set("nemo variables", nemo_var_mappings);
+    eckit::LocalConfiguration interp_conf;
+    interp_conf.set("type", "unstructured-bilinear-lonlat");
+    interp_conf.set("non_linear", "missing-if-all-missing-real32");
+    settings_map["ORCA2_T"].interpolator_config.set("atlas-interpolator", interp_conf);
+
+    std::vector<std::string> state_variables {
+        "sea_ice_area_fraction",
+        "sea_surface_foundation_temperature",
+        "sea_water_potential_temperature"};
+    settings_map["ORCA2_T"].state_config.set("state variables", state_variables);
+    settings_map["ORCA2_T"].state_config.set("date", "2021-06-30T00:00:00Z");
+    settings_map["ORCA2_T"].state_config.set("nemo field file", "../Data/orca2_t_nemo.nc");
+    settings_map["ORCA2_T"].state_config.set("variance field file", "../Data/orca2_t_bkg_var.nc");
+
+    settings_map["ORCA2_T"].lons = std::vector<double>{0, 120, 270};
+    settings_map["ORCA2_T"].lats = std::vector<double>{88, 0, 30};
+
+    settings_map["ORCA2_T"].surf_vars = oops::Variables{{
+        oops::Variable{"sea_ice_area_fraction"},
+        oops::Variable{"sea_surface_foundation_temperature"}}};
+    settings_map["ORCA2_T"].surf_values = std::vector<double>{
+        1,             missing_value, 0,
+        18.4888916016, missing_value, 18.1592999503};
+
+    settings_map["ORCA2_T"].vol_vars = oops::Variables{
+        {oops::Variable{"sea_water_potential_temperature"}}};
+    settings_map["ORCA2_T"].vol_values = std::vector<double>{
+        18.4888916016, missing_value, 18.1592999503,
+        17.9419364929, missing_value, 17.75000288,
+        missing_value, missing_value, missing_value};
   }
 
-  Interpolator interpolator(interpolator_conf, geometry, lats, lons);
-
-  SECTION("test interpolator.apply fails missing variable") {
-    oops::Variables variables{{oops::Variable{"NOTAVARIABLE"}}};
-    std::vector<double> vals(3);
-    std::vector<bool> mask(3, true);
-    EXPECT_THROWS_AS(interpolator.apply(variables, state, mask, vals),
-        eckit::BadParameter);
+  // AMM1 settings
+  {
+    settings_map["AMM1"].nlocs = 3;
+    settings_map["AMM1"].nlevs = 3;
+    settings_map["AMM1"].geometry_config.set("grid name", "../Data/amm1_atlas_grid_spec.yaml");
+    settings_map["AMM1"].geometry_config.set("number levels", settings_map["AMM1"].nlevs);
+
+    std::vector<eckit::LocalConfiguration> nemo_var_mappings(3);
+    nemo_var_mappings[0].set("name", "sea_surface_height_anomaly")
+      .set("nemo field name", "sossheig")
+      .set("model space", "surface");
+    nemo_var_mappings[1].set("name", "mass_concentration_of_chlorophyll_in_sea_water")
+      .set("nemo field name", "CHL")
+      .set("model space", "surface");
+    nemo_var_mappings[2].set("name", "sea_water_potential_temperature")
+      .set("nemo field name", "votemper")
+      .set("model space", "volume");
+    settings_map["AMM1"].geometry_config.set("nemo variables", nemo_var_mappings);
+    eckit::LocalConfiguration interp_conf;
+    interp_conf.set("type", "unstructured-bilinear-lonlat");
+    interp_conf.set("non_linear", "missing-if-all-missing-real32");
+    settings_map["AMM1"].interpolator_config.set("atlas-interpolator", interp_conf);
+
+    std::vector<std::string> state_variables {
+        "sea_surface_height_anomaly",
+        "mass_concentration_of_chlorophyll_in_sea_water",
+        "sea_water_potential_temperature"};
+    settings_map["AMM1"].state_config.set("state variables", state_variables);
+    settings_map["AMM1"].state_config.set("date", "2021-06-30T00:00:00Z");
+    settings_map["AMM1"].state_config.set("nemo field file", "../Data/amm1_nemo.nc");
+    settings_map["AMM1"].state_config.set("variance field file", "../Data/amm1_nemo.nc");
+
+    settings_map["AMM1"].lons = std::vector<double>{-17.5, -6.78, -16.1};
+    settings_map["AMM1"].lats = std::vector<double>{58.16, 58.91, 63.55};
+
+    settings_map["AMM1"].surf_vars = oops::Variables{{
+        oops::Variable{"sea_surface_height_anomaly"},
+        oops::Variable{"mass_concentration_of_chlorophyll_in_sea_water"}}};
+    settings_map["AMM1"].surf_values = std::vector<double>{
+        -0.470139563084, -0.286416769028, -0.433749824762,
+        0.345775008202,  0.83959633112,  2.09327220917};
+
+    settings_map["AMM1"].vol_vars = oops::Variables{
+        {oops::Variable{"sea_water_potential_temperature"}}};
+    settings_map["AMM1"].vol_values = std::vector<double>{
+        11.9501609802, 13.9884538651, 10.1916904449,
+        11.9500904083, 13.8567619324, 10.1912517548,
+        11.9499549866, 13.7289009094, 10.1908035278};
   }
 
-  SECTION("test interpolator.apply") {
-    // two variables at n locations
-    std::vector<double> vals(2*nlocs);
-    std::vector<bool> mask(nlocs, true);
-    interpolator.apply(oops::Variables{{oops::Variable{"sea_ice_area_fraction"},
-        oops::Variable{"sea_surface_foundation_temperature"}}}, state, mask, vals);
-
-    double missing_value = util::missingValue<double>();
-    std::vector<double> testvals = {1, missing_value, 0, 18.4888916016,
-                                    missing_value, 18.1592999503};
-
-    for (size_t i=0; i < testvals.size(); ++i) {
-      std::cout << "vals[" << i << "] " << std::setprecision(12) << vals[i]
-                << " testvals[" << i << "] " << testvals[i] << std::endl;
-      EXPECT(std::abs(vals[i] - testvals[i]) < ATOL);
+  for (const auto& [key, settings] : settings_map) {
+    Geometry geometry(settings.geometry_config, eckit::mpi::comm());
+
+    OrcaInterpolatorParameters params;
+    params.validateAndDeserialize(settings.interpolator_config);
+
+    // create a state from the test data
+    OrcaStateParameters stateParams;
+    stateParams.validateAndDeserialize(settings.state_config);
+    State state(geometry, stateParams);
+
+    SECTION("test " + key + " interpolator succeeds even with no locations") {
+      Interpolator interpolator(settings.interpolator_config, geometry, {}, {});
     }
-  }
-  SECTION("test interpolator.apply multiple levels") {
-    std::vector<double> vals(nlevs*nlocs);
-    std::vector<bool> mask(nlocs, true);
-    interpolator.apply(oops::Variables{{oops::Variable{"sea_water_potential_temperature"}}},
-                                       state, mask, vals);
-
-    double missing_value = util::missingValue<double>();
-    std::vector<double> testvals = {
-      18.4888916016, missing_value, 18.1592999503,
-      17.9419364929, missing_value, 17.75000288,
-      missing_value, missing_value, missing_value};
-
-    for (size_t i=0; i < testvals.size(); ++i) {
-      std::cout << "vals[" << i << "] " << std::setprecision(12) << vals[i]
-                << " testvals[" << i << "] " << testvals[i] << std::endl;
-      EXPECT(std::abs(vals[i] - testvals[i]) < ATOL);
+
+    Interpolator interpolator(settings.interpolator_config, geometry, settings.lats, settings.lons);
+
+    SECTION("test " + key + " interpolator.apply fails missing variable") {
+      oops::Variables variables{{oops::Variable{"NOTAVARIABLE"}}};
+      std::vector<double> vals(3);
+      std::vector<bool> mask(3, true);
+      EXPECT_THROWS_AS(interpolator.apply(variables, state, mask, vals),
+          eckit::BadParameter);
+    }
+
+    SECTION("test " + key + " interpolator.apply") {
+      // two variables at n locations
+      std::vector<double> vals(2*settings.nlocs);
+      std::vector<bool> mask(settings.nlocs, true);
+
+      interpolator.apply(settings.surf_vars, state, mask, vals);
+
+      for (size_t i=0; i < settings.surf_values.size(); ++i) {
+        std::cout << "vals[" << i << "] " << std::setprecision(12) << vals[i]
+                  << " kgo_values[" << i << "] " << settings.surf_values[i] << std::endl;
+      }
+      for (size_t i=0; i < settings.surf_values.size(); ++i) {
+        EXPECT(std::abs(vals[i] - settings.surf_values[i]) < ATOL);
+      }
+    }
+    SECTION("test " + key + " interpolator.apply multiple levels") {
+      std::vector<double> vals(settings.nlevs*settings.nlocs);
+      std::vector<bool> mask(settings.nlocs, true);
+
+      interpolator.apply(settings.vol_vars, state, mask, vals);
+
+      for (size_t i=0; i < settings.vol_values.size(); ++i) {
+        std::cout << "vals[" << i << "] " << std::setprecision(12) << vals[i]
+                  << " kgo_values[" << i << "] " << settings.vol_values[i] << std::endl;
+      }
+      for (size_t i=0; i < settings.vol_values.size(); ++i) {
+        EXPECT(std::abs(vals[i] - settings.vol_values[i]) < ATOL);
+      }
     }
   }
 }
diff --git a/src/tests/orca-jedi/test_nemo_io_read_server.cc b/src/tests/orca-jedi/test_nemo_io_read_server.cc
index d3d1d0c..e6baa33 100644
--- a/src/tests/orca-jedi/test_nemo_io_read_server.cc
+++ b/src/tests/orca-jedi/test_nemo_io_read_server.cc
@@ -2,6 +2,8 @@
  * (C) British Crown Copyright 2024 Met Office
  */
 
+#include <fstream>
+
 #include "eckit/log/Bytes.h"
 #include "eckit/testing/Test.h"
 
@@ -22,7 +24,7 @@
 
 #include "orca-jedi/nemo_io/ReadServer.h"
 #include "orca-jedi/nemo_io/NemoFieldReader.h"
-#include "orca-jedi/nemo_io/OrcaIndex.h"
+#include "orca-jedi/nemo_io/AtlasIndex.h"
 
 #include "tests/orca-jedi/OrcaModelTestEnvironment.h"
 
@@ -31,8 +33,111 @@ namespace test {
 
 
 //-----------------------------------------------------------------------------
+//
+CASE("test MPI distributed reads structured grid field array view") {
+  eckit::PathName test_data_path("../Data/amm1_nemo.nc");
+  eckit::PathName grid_spec_path("../Data/amm1_atlas_grid_spec.yaml");
+  auto partitioner_names = std::vector<std::string>{"serial"};
+  for (std::string partitioner_name : partitioner_names) {
+    atlas::Grid grid{atlas::Grid::Spec{grid_spec_path}};
+
+    auto meshgen_config = grid.meshgenerator();
+    std::cout << "meshgen_config: " << meshgen_config << std::endl;
+    atlas::MeshGenerator meshgen(meshgen_config);
+
+    auto partitioner_config = grid.partitioner();
+    partitioner_config.set("type", partitioner_name);
+    auto partitioner = atlas::grid::Partitioner(partitioner_config);
+
+    auto mesh = meshgen.generate(grid, partitioner);
+    std::unique_ptr<AtlasIndexToBufferIndex> atlas_index(
+      AtlasIndexToBufferIndexCreator::create_unique(grid.type(), mesh));
+
+    auto funcSpace = atlas::functionspace::NodeColumns(mesh);
+
+    auto eckit_timer = std::make_shared<eckit::Timer>(
+      "Geometry(ORCA): ", oops::Log::trace());
+    eckit_timer->start();
+    ReadServer read_server(eckit_timer, test_data_path, mesh);
+
+    constexpr float kgo_missing_value = -32768.;
+    SECTION(partitioner_name + " get field _FillValue") {
+      auto missing_value = read_server.read_fillvalue<float>("sossheig");
+      if (std::abs(missing_value - kgo_missing_value) >= 1e-6) {
+        std::cout << "missing_value: " << missing_value << " == " << kgo_missing_value
+            << " diff " << std::abs(missing_value - kgo_missing_value) << std::endl;
+      }
+      EXPECT(std::abs(missing_value - kgo_missing_value) < 1e-6);
+
+      auto default_missing_value = read_server.read_fillvalue<float>("nav_lev");
+      if (std::abs(default_missing_value - std::numeric_limits<float>::lowest()) >= 1e-6) {
+        std::cout << "default_missing_value: " << missing_value << " == "
+            << std::numeric_limits<float>::lowest()
+            << " diff " << std::abs(default_missing_value - std::numeric_limits<float>::lowest())
+            << std::endl;
+      }
+      EXPECT(std::abs(default_missing_value - std::numeric_limits<float>::lowest()) < 1e-6);
+    }
+
+    SECTION(partitioner_name + " surface field") {
+      std::string nemo_name{"sossheig"};
+      atlas::Field field(funcSpace.createField<float>(
+            atlas::option::name(nemo_name) |
+            atlas::option::levels(1)));
+      auto field_view = atlas::array::make_view<float, 2>(field);
+      field.metadata().set("missing_value", kgo_missing_value);
+      field.metadata().set("missing_value_type", "approximately-equals");
+      field.metadata().set("missing_value_epsilon", 1e-6);
+
+      read_server.read_var<float>(nemo_name, 0, field_view);
+
+      std::vector<float> raw_data;
+      {
+        NemoFieldReader field_reader(test_data_path);
+        raw_data = field_reader.read_var_slice<float>(nemo_name, 0, 0);
+      }
+
+      auto ghost = atlas::array::make_view<int32_t, 1>(mesh.nodes().ghost());
+      for (size_t i = 0; i < field_view.size(); ++i) {
+        if (ghost(i)) continue;
+        float raw_value = raw_data[(*atlas_index)(i)];
+        if (raw_value != field_view(i, 0)) {
+          std::cout << "mismatch: i " << i
+            << " " << raw_value << " != " << field_view(i, 0) << std::endl;
+        }
+        EXPECT_EQUAL(raw_value, field_view(i, 0));
+      }
+    }
+
+    SECTION(partitioner_name + " volume field") {
+      atlas::Field field(funcSpace.createField<double>(
+            atlas::option::name("votemper") |
+            atlas::option::levels(3)));
+      auto field_view = atlas::array::make_view<double, 2>(field);
+
+      NemoFieldReader field_reader(test_data_path);
+      read_server.read_var<double>("votemper", 0, field_view);
+
+      auto ghost = atlas::array::make_view<int32_t, 1>(mesh.nodes().ghost());
+      std::vector<double> raw_data;
+      for (int k =0; k <3; k++) {
+        raw_data = field_reader.read_var_slice<double>("votemper", 0, k);
+        for (int i = 0; i < field_view.shape(0); ++i) {
+          double raw_value = raw_data[(*atlas_index)(i)];
+          if (ghost(i)) continue;
+          if (raw_value != field_view(i, k)) {
+            std::cout << "mismatch: "
+              << raw_value << " "
+              << " with mesh " << field_view(i, k) << std::endl;
+          }
+          EXPECT_EQUAL(raw_value, field_view(i, k));
+        }
+      }
+    }
+  }
+}
 
-CASE("test MPI distributed reads field array view") {
+CASE("test MPI distributed reads orca grid field array view") {
   eckit::PathName test_data_path("../Data/orca2_t_nemo.nc");
 
   auto partitioner_names = std::vector<std::string>{"serial", "checkerboard"};
@@ -46,7 +151,8 @@ CASE("test MPI distributed reads field array view") {
     auto partitioner = atlas::grid::Partitioner(partitioner_config);
     auto mesh = meshgen.generate(grid, partitioner);
     auto funcSpace = atlas::functionspace::NodeColumns(mesh);
-    auto orca_index = OrcaIndexToBufferIndex(mesh.grid());
+    std::unique_ptr<AtlasIndexToBufferIndex> orca_index(
+      AtlasIndexToBufferIndexCreator::create_unique(grid.type(), mesh));
 
     std::vector<double> raw_data;
     {
@@ -68,17 +174,19 @@ CASE("test MPI distributed reads field array view") {
     auto ij = atlas::array::make_view<int32_t, 2>(mesh.nodes().field("ij"));
 
     auto ghost = atlas::array::make_view<int32_t, 1>(mesh.nodes().ghost());
-    for (size_t i = 0; i < field_view.size(); ++i) {
-      if (ghost(i)) continue;
-      double raw_value = raw_data[orca_index(ij(i, 0), ij(i, 1))];
-      if (raw_value != field_view(i, 0)) {
-        std::cout << "mismatch: "
-          << " ij(" << i << ", 0) " << ij(i, 0)
-          << " ij(" << i << ", 1) " << ij(i, 1)
-          << " 1 proc " << raw_value << " "
-          << " with mesh " << field_view(i, 0) << std::endl;
+    SECTION(partitioner_name + " surface field") {
+      for (size_t i = 0; i < field_view.size(); ++i) {
+        if (ghost(i)) continue;
+        double raw_value = raw_data[(*orca_index)(ij(i, 0), ij(i, 1))];
+        if (raw_value != field_view(i, 0)) {
+          std::cout << "mismatch: "
+            << " ij(" << i << ", 0) " << ij(i, 0)
+            << " ij(" << i << ", 1) " << ij(i, 1)
+            << " 1 proc " << raw_value << " "
+            << " with mesh " << field_view(i, 0) << std::endl;
+        }
+        EXPECT_EQUAL(raw_value, field_view(i, 0));
       }
-      EXPECT_EQUAL(raw_value, field_view(i, 0));
     }
 
     SECTION(partitioner_name + " volume field") {
@@ -97,7 +205,7 @@ CASE("test MPI distributed reads field array view") {
       for (int k =0; k <3; k++) {
         raw_data = field_reader.read_var_slice<double>("votemper", 0, k);
         for (int i = 0; i < field_view.shape(0); ++i) {
-          double raw_value = raw_data[orca_index(ij(i, 0), ij(i, 1))];
+          double raw_value = raw_data[(*orca_index)(ij(i, 0), ij(i, 1))];
           if (ghost(i)) continue;
           if (raw_value != field_view(i, k)) {
             std::cout << "mismatch: "
diff --git a/src/tests/orca-jedi/test_nemo_io_write_server.cc b/src/tests/orca-jedi/test_nemo_io_write_server.cc
index bad7e72..7da4c92 100644
--- a/src/tests/orca-jedi/test_nemo_io_write_server.cc
+++ b/src/tests/orca-jedi/test_nemo_io_write_server.cc
@@ -22,7 +22,7 @@
 
 #include "orca-jedi/nemo_io/WriteServer.h"
 #include "orca-jedi/nemo_io/NemoFieldReader.h"
-#include "orca-jedi/nemo_io/OrcaIndex.h"
+#include "orca-jedi/nemo_io/AtlasIndex.h"
 
 #include "tests/orca-jedi/OrcaModelTestEnvironment.h"
 
@@ -56,181 +56,191 @@ void applyMPISerialised(const Functor& functor, const std::string& partitioner_n
 
 CASE("test MPI distributed field array view write to disk") {
   auto partitioner_names = std::vector<std::string>{"serial", "checkerboard"};
-  std::string grid_name("ORCA2_T");
+  auto grid_names = std::vector<std::string>{"ORCA2_T",
+     "../Data/amm1_atlas_grid_spec.yaml"};
   size_t nparts = atlas::mpi::size();
   if (nparts == 1)
     partitioner_names = std::vector<std::string>{"serial"};
-  for (std::string partitioner_name : partitioner_names) {
-    // setup atlas data before write
-    eckit::PathName test_data_path(std::string("../testoutput/write_orca2_t_")
-        + partitioner_name + "_" + std::to_string(nparts) + ".nc");
-    atlas::OrcaGrid grid(grid_name);
-
-    auto meshgen_config = grid.meshgenerator();
-    atlas::MeshGenerator meshgen(meshgen_config);
-    auto partitioner_config = grid.partitioner();
-    partitioner_config.set("type", partitioner_name);
-    auto partitioner = atlas::grid::Partitioner(partitioner_config);
-    auto mesh = meshgen.generate(grid, partitioner);
-
-    auto funcSpace = atlas::functionspace::NodeColumns(mesh);
-    auto orca2buffer = OrcaIndexToBufferIndex(mesh);
-
-    if (grid_name == "ORCA2_T") {
-      SECTION(partitioner_name + "_" + std::to_string(nparts) + " test ORCA indexing") {
-        EXPECT(orca2buffer.nx() == 182);
-        EXPECT(orca2buffer.ny() == 149);
+  for (std::string grid_name : grid_names) {
+    for (std::string partitioner_name : partitioner_names) {
+      // setup atlas data before write
+      eckit::PathName test_data_path;
+      eckit::PathName grid_spec_path(grid_name);
+      atlas::Grid grid;
+      if (grid_spec_path.exists()) {
+        grid = atlas::Grid{atlas::Grid::Spec{grid_spec_path}};
+        test_data_path = eckit::PathName{std::string("../testoutput/write_" + grid.type() + "_")
+          + partitioner_name + "_" + std::to_string(nparts) + ".nc"};
+      } else {
+        grid = atlas::OrcaGrid{grid_name};
+        test_data_path = eckit::PathName{std::string("../testoutput/write_" + grid_name + "_")
+            + partitioner_name + "_" + std::to_string(nparts) + ".nc"};
       }
-    }
 
-    auto ghost = atlas::array::make_view<int32_t, 1>(mesh.nodes().ghost());
-
-    auto field_ice = funcSpace.createField<double>(
-      atlas::option::name("iiceconc") | atlas::option::levels(1));
-    auto field_t0_temp = funcSpace.createField<float>(
-      atlas::option::name("votemper0") | atlas::option::levels(3));
-    auto field_t1_temp = funcSpace.createField<float>(
-      atlas::option::name("votemper1") | atlas::option::levels(3));
-
-    field_ice.metadata().set("missing_value", 0);
-    field_ice.metadata().set("missing_value_type", "approximately-equals");
-    field_ice.metadata().set("missing_value_epsilon", 1e-6);
-    field_t0_temp.metadata().set("missing_value", 0);
-    field_t0_temp.metadata().set("missing_value_type", "approximately-equals");
-    field_t0_temp.metadata().set("missing_value_epsilon", 1e-6);
-    field_t1_temp.metadata().set("missing_value", 0);
-    field_t1_temp.metadata().set("missing_value_type", "approximately-equals");
-    field_t1_temp.metadata().set("missing_value_epsilon", 1e-6);
-
-    atlas::field::MissingValue ice_mv(field_ice);
-    atlas::field::MissingValue temp0_mv(field_t0_temp);
-    atlas::field::MissingValue temp1_mv(field_t1_temp);
-
-    auto field_view_ice = atlas::array::make_view<double, 2>(field_ice);
-    auto field_view_t0_temp = atlas::array::make_view<float, 2>(field_t0_temp);
-    auto field_view_t1_temp = atlas::array::make_view<float, 2>(field_t1_temp);
-
-    // fill fields
-    {
-      const size_t num_nodes = field_view_ice.shape(0);
-      for (size_t inode = 0; inode < num_nodes; ++inode) {
-        if (ghost(inode)) continue;
-        field_view_ice(inode, 0) = inode;
-        for (size_t ilevel = 0; ilevel < 3; ++ilevel) {
-          field_view_t0_temp(inode, ilevel) = inode + ilevel*5;
-          field_view_t1_temp(inode, ilevel) = inode + ilevel*10;
+      auto meshgen_config = grid.meshgenerator();
+      atlas::MeshGenerator meshgen(meshgen_config);
+      auto partitioner_config = grid.partitioner();
+      partitioner_config.set("type", partitioner_name);
+      auto partitioner = atlas::grid::Partitioner(partitioner_config);
+      auto mesh = meshgen.generate(grid, partitioner);
+
+      auto funcSpace = atlas::functionspace::NodeColumns(mesh);
+      std::unique_ptr<AtlasIndexToBufferIndex> atlas2buffer(
+        AtlasIndexToBufferIndexCreator::create_unique(grid.type(), mesh));
+
+      if (grid_name == "ORCA2_T") {
+        SECTION(partitioner_name + "_" + std::to_string(nparts) + " test ORCA indexing") {
+          EXPECT(atlas2buffer->nx() == 182);
+          EXPECT(atlas2buffer->ny() == 149);
         }
       }
-    }
-
-    funcSpace.haloExchange(field_ice);
-    funcSpace.haloExchange(field_t0_temp);
-    funcSpace.haloExchange(field_t1_temp);
-
-    SECTION(partitioner_name + "_" + std::to_string(nparts)
-             + " write data to file via write server") {
-      std::vector<util::DateTime> datetimes{util::DateTime("1970-01-01T00:00:00Z"),
-                                            util::DateTime("1970-01-02T00:00:00Z")};
-      std::shared_ptr<eckit::Timer> eckit_timer =
-          std::make_shared<eckit::Timer>("write_server tests: ", oops::Log::debug());
-      const bool serial_distribution = (atlas::mpi::size() == 1 || partitioner_name == "serial");
-      WriteServer writer(eckit_timer, test_data_path, mesh,
-                         datetimes, {1, 2, 3}, serial_distribution);
-
-      writer.write_surf_var<double>("iiceconc", 0, ice_mv, field_view_ice);
-      writer.write_vol_var<float>("votemper", 0, temp0_mv, field_view_t0_temp);
-      writer.write_vol_var<float>("votemper", 1, temp1_mv, field_view_t1_temp);
-    }
-
-    SECTION(partitioner_name + "_" + std::to_string(nparts) + " File exists") {
-      auto check_file_exists = [&](){
-        // wait up to 20 seconds for the file system...
-        for (int wait_count=0; wait_count < 10; ++wait_count) {
-          if (test_data_path.exists()) break;
-          sleep(2);
-        }
-        EXPECT(test_data_path.exists());
-      };
-      applyMPISerialised(check_file_exists, partitioner_name);
-    }
-
-    SECTION(partitioner_name + "_" + std::to_string(nparts) + " check latitude/longitude data") {
-      auto check_lon_lat = [&](){
-        NemoFieldReader field_reader(test_data_path);
-        std::vector<atlas::PointXY> data = field_reader.read_locs();
 
-        auto lonlat = atlas::array::make_view<double, 2>(mesh.nodes().lonlat());
-        auto ij = atlas::array::make_view<int32_t, 2>(mesh.nodes().field("ij"));
-        const size_t num_nodes = ij.shape(0);
-        EXPECT(num_nodes <= data.size());
+      auto ghost = atlas::array::make_view<int32_t, 1>(mesh.nodes().ghost());
+
+      auto field_ice = funcSpace.createField<double>(
+        atlas::option::name("iiceconc") | atlas::option::levels(1));
+      auto field_t0_temp = funcSpace.createField<float>(
+        atlas::option::name("votemper0") | atlas::option::levels(3));
+      auto field_t1_temp = funcSpace.createField<float>(
+        atlas::option::name("votemper1") | atlas::option::levels(3));
+
+      field_ice.metadata().set("missing_value", 0);
+      field_ice.metadata().set("missing_value_type", "approximately-equals");
+      field_ice.metadata().set("missing_value_epsilon", 1e-6);
+      field_t0_temp.metadata().set("missing_value", 0);
+      field_t0_temp.metadata().set("missing_value_type", "approximately-equals");
+      field_t0_temp.metadata().set("missing_value_epsilon", 1e-6);
+      field_t1_temp.metadata().set("missing_value", 0);
+      field_t1_temp.metadata().set("missing_value_type", "approximately-equals");
+      field_t1_temp.metadata().set("missing_value_epsilon", 1e-6);
+
+      atlas::field::MissingValue ice_mv(field_ice);
+      atlas::field::MissingValue temp0_mv(field_t0_temp);
+      atlas::field::MissingValue temp1_mv(field_t1_temp);
+
+      auto field_view_ice = atlas::array::make_view<double, 2>(field_ice);
+      auto field_view_t0_temp = atlas::array::make_view<float, 2>(field_t0_temp);
+      auto field_view_t1_temp = atlas::array::make_view<float, 2>(field_t1_temp);
+
+      // fill fields
+      {
+        const size_t num_nodes = field_view_ice.shape(0);
         for (size_t inode = 0; inode < num_nodes; ++inode) {
           if (ghost(inode)) continue;
-          const int64_t ibuf = orca2buffer(ij(inode, 0), ij(inode, 1));
-          EXPECT(data[ibuf](0) == lonlat(inode, 0));
-          EXPECT(data[ibuf](1) == lonlat(inode, 1));
+          field_view_ice(inode, 0) = inode;
+          for (size_t ilevel = 0; ilevel < 3; ++ilevel) {
+            field_view_t0_temp(inode, ilevel) = inode + ilevel*5;
+            field_view_t1_temp(inode, ilevel) = inode + ilevel*10;
+          }
         }
-      };
-      applyMPISerialised(check_lon_lat, partitioner_name);
-    }
+      }
 
-    SECTION(partitioner_name + "_" + std::to_string(nparts) + " check surface data") {
-      auto check_surface_data = [&](){
-        NemoFieldReader field_reader(test_data_path);
-        std::vector<double> data = field_reader.read_var_slice<double>("iiceconc", 0, 0);
+      funcSpace.haloExchange(field_ice);
+      funcSpace.haloExchange(field_t0_temp);
+      funcSpace.haloExchange(field_t1_temp);
+
+      SECTION(partitioner_name + "_" + std::to_string(nparts)
+               + " write data to file via write server") {
+        std::vector<util::DateTime> datetimes{util::DateTime("1970-01-01T00:00:00Z"),
+                                              util::DateTime("1970-01-02T00:00:00Z")};
+        std::shared_ptr<eckit::Timer> eckit_timer =
+            std::make_shared<eckit::Timer>("write_server tests: ", oops::Log::debug());
+        const bool serial_distribution = (atlas::mpi::size() == 1 || partitioner_name == "serial");
+        WriteServer writer(eckit_timer, test_data_path, mesh,
+                           datetimes, {1, 2, 3}, serial_distribution);
+
+        writer.write_surf_var<double>("iiceconc", 0, ice_mv, field_view_ice);
+        writer.write_vol_var<float>("votemper", 0, temp0_mv, field_view_t0_temp);
+        writer.write_vol_var<float>("votemper", 1, temp1_mv, field_view_t1_temp);
+      }
 
-        auto ij = atlas::array::make_view<int32_t, 2>(mesh.nodes().field("ij"));
-        const size_t num_nodes = ij.shape(0);
-        EXPECT(num_nodes <= data.size());
-        for (size_t inode = 0; inode < num_nodes; ++inode) {
-          if (ghost(inode)) continue;
-          const int64_t ibuf = orca2buffer(ij(inode, 0), ij(inode, 1));
-          if (ice_mv(field_view_ice(inode, 0))) {
-            EXPECT(std::abs(data[ibuf] - 1e+20) < 1e-6);
-          } else {
-            EXPECT(data[ibuf] == field_view_ice(inode, 0));
+      SECTION(partitioner_name + "_" + std::to_string(nparts) + " File exists") {
+        auto check_file_exists = [&](){
+          // wait up to 20 seconds for the file system...
+          for (int wait_count=0; wait_count < 10; ++wait_count) {
+            if (test_data_path.exists()) break;
+            sleep(2);
           }
-        }
-      };
-      applyMPISerialised(check_surface_data, partitioner_name);
-    }
+          EXPECT(test_data_path.exists());
+        };
+        applyMPISerialised(check_file_exists, partitioner_name);
+      }
+
+      SECTION(partitioner_name + "_" + std::to_string(nparts) + " check latitude/longitude data") {
+        auto check_lon_lat = [&](){
+          NemoFieldReader field_reader(test_data_path);
+          std::vector<atlas::PointXY> data = field_reader.read_locs();
 
-    SECTION(partitioner_name + "_" + std::to_string(nparts) + " check volume data") {
-      auto check_volume_data = [&](){
-        NemoFieldReader field_reader(test_data_path);
-        std::vector<float> data_t0_l0 = field_reader.read_var_slice<float>("votemper", 0, 0);
-        std::vector<float> data_t0_l1 = field_reader.read_var_slice<float>("votemper", 0, 1);
-        std::vector<float> data_t0_l2 = field_reader.read_var_slice<float>("votemper", 0, 2);
-        std::vector<float> data_t1_l0 = field_reader.read_var_slice<float>("votemper", 1, 0);
-        std::vector<float> data_t1_l1 = field_reader.read_var_slice<float>("votemper", 1, 1);
-        std::vector<float> data_t1_l2 = field_reader.read_var_slice<float>("votemper", 1, 2);
-
-        auto ij = atlas::array::make_view<int32_t, 2>(mesh.nodes().field("ij"));
-
-        auto check_slice = [ghost, ij, orca2buffer](const std::vector<float>& test_data,
-            const atlas::array::ArrayView<float, 2>& original_fv,
-            const atlas::field::MissingValue mv,
-            const size_t ilev) {
-          const size_t num_nodes = ij.shape(0);
-          EXPECT(num_nodes <= test_data.size());
+          auto lonlat = atlas::array::make_view<double, 2>(mesh.nodes().lonlat());
+          const size_t num_nodes = mesh.nodes().size();
+          EXPECT(num_nodes <= data.size());
           for (size_t inode = 0; inode < num_nodes; ++inode) {
             if (ghost(inode)) continue;
-            const int64_t ibuf = orca2buffer(ij(inode, 0), ij(inode, 1));
-            if (mv(original_fv(inode, ilev))) continue;
-            if (std::abs(test_data[ibuf] - original_fv(inode, ilev)) >= 1e-7)
-              std::cout << "test_data[" << ibuf << "] " << test_data[ibuf]
-                        << " original_fv(" << inode << ", " << ilev << ") "
-                        << original_fv(inode, ilev) << std::endl;
-            EXPECT(std::abs(test_data[ibuf] - original_fv(inode, ilev)) < 1e-7);
+            const int64_t ibuf = (*atlas2buffer)(inode);
+            EXPECT(data[ibuf](0) == lonlat(inode, 0));
+            EXPECT(data[ibuf](1) == lonlat(inode, 1));
           }
         };
-        check_slice(data_t0_l0, field_view_t0_temp, temp0_mv, 0);
-        check_slice(data_t0_l1, field_view_t0_temp, temp0_mv, 1);
-        check_slice(data_t0_l2, field_view_t0_temp, temp0_mv, 2);
-        check_slice(data_t1_l0, field_view_t1_temp, temp1_mv, 0);
-        check_slice(data_t1_l1, field_view_t1_temp, temp1_mv, 1);
-        check_slice(data_t1_l2, field_view_t1_temp, temp1_mv, 2);
-      };
-      applyMPISerialised(check_volume_data, partitioner_name);
+        applyMPISerialised(check_lon_lat, partitioner_name);
+      }
+
+      SECTION(partitioner_name + "_" + std::to_string(nparts) + " check surface data") {
+        auto check_surface_data = [&](){
+          NemoFieldReader field_reader(test_data_path);
+          std::vector<double> data = field_reader.read_var_slice<double>("iiceconc", 0, 0);
+
+          const size_t num_nodes = mesh.nodes().size();
+          EXPECT(num_nodes <= data.size());
+          for (size_t inode = 0; inode < num_nodes; ++inode) {
+            if (ghost(inode)) continue;
+            const int64_t ibuf = (*atlas2buffer)(inode);
+            if (ice_mv(field_view_ice(inode, 0))) {
+              EXPECT(std::abs(data[ibuf] - 1e+20) < 1e-6);
+            } else {
+              EXPECT(data[ibuf] == field_view_ice(inode, 0));
+            }
+          }
+        };
+        applyMPISerialised(check_surface_data, partitioner_name);
+      }
+
+      SECTION(partitioner_name + "_" + std::to_string(nparts) + " check volume data") {
+        auto check_volume_data = [&](){
+          NemoFieldReader field_reader(test_data_path);
+          std::vector<float> data_t0_l0 = field_reader.read_var_slice<float>("votemper", 0, 0);
+          std::vector<float> data_t0_l1 = field_reader.read_var_slice<float>("votemper", 0, 1);
+          std::vector<float> data_t0_l2 = field_reader.read_var_slice<float>("votemper", 0, 2);
+          std::vector<float> data_t1_l0 = field_reader.read_var_slice<float>("votemper", 1, 0);
+          std::vector<float> data_t1_l1 = field_reader.read_var_slice<float>("votemper", 1, 1);
+          std::vector<float> data_t1_l2 = field_reader.read_var_slice<float>("votemper", 1, 2);
+
+          auto check_slice = [ghost, mesh](const std::vector<float>& test_data,
+              const std::unique_ptr<AtlasIndexToBufferIndex>& o2b,
+              const atlas::array::ArrayView<float, 2>& original_fv,
+              const atlas::field::MissingValue mv,
+              const size_t ilev) {
+            const size_t num_nodes = mesh.nodes().size();
+            EXPECT(num_nodes <= test_data.size());
+            for (size_t inode = 0; inode < num_nodes; ++inode) {
+              if (ghost(inode)) continue;
+              const int64_t ibuf = (*o2b)(inode);
+              if (mv(original_fv(inode, ilev))) continue;
+              if (std::abs(test_data[ibuf] - original_fv(inode, ilev)) >= 1e-7)
+                std::cout << "test_data[" << ibuf << "] " << test_data[ibuf]
+                          << " original_fv(" << inode << ", " << ilev << ") "
+                          << original_fv(inode, ilev) << std::endl;
+              EXPECT(std::abs(test_data[ibuf] - original_fv(inode, ilev)) < 1e-7);
+            }
+          };
+          check_slice(data_t0_l0, atlas2buffer, field_view_t0_temp, temp0_mv, 0);
+          check_slice(data_t0_l1, atlas2buffer, field_view_t0_temp, temp0_mv, 1);
+          check_slice(data_t0_l2, atlas2buffer, field_view_t0_temp, temp0_mv, 2);
+          check_slice(data_t1_l0, atlas2buffer, field_view_t1_temp, temp1_mv, 0);
+          check_slice(data_t1_l1, atlas2buffer, field_view_t1_temp, temp1_mv, 1);
+          check_slice(data_t1_l2, atlas2buffer, field_view_t1_temp, temp1_mv, 2);
+        };
+        applyMPISerialised(check_volume_data, partitioner_name);
+      }
     }
   }
 }
diff --git a/src/tests/testinput/CMakeLists.txt b/src/tests/testinput/CMakeLists.txt
index e7c0964..ca1ee29 100644
--- a/src/tests/testinput/CMakeLists.txt
+++ b/src/tests/testinput/CMakeLists.txt
@@ -9,9 +9,11 @@ list( APPEND orcajedi_test_input
   hofx_nc_ssh.yaml
   hofx_nc_ssh_checkerboard.yaml
   hofx_nc_ssh_eorca025.yaml
+  hofx_nc_ssh_amm1.yaml
   hofx3d_nc_sst.yaml
   hofx3d_nc_potm.yaml
   hofx3d_nc_prof_2vars.yaml
+  hofx3d_nc_ssh_amm1.yaml
   odb_ice_query.yaml
   odb_ice_mapping.yaml
   test_name_map.yaml
diff --git a/src/tests/testinput/hofx3d_nc_ssh_amm1.yaml b/src/tests/testinput/hofx3d_nc_ssh_amm1.yaml
new file mode 100644
index 0000000..ef10910
--- /dev/null
+++ b/src/tests/testinput/hofx3d_nc_ssh_amm1.yaml
@@ -0,0 +1,48 @@
+# (C) British Crown Copyright 2024 Met Office
+
+time window:
+  begin: 2021-06-28T23:00:00Z
+  length: P2D
+geometry :
+  nemo variables:
+    - name: sea_surface_height_anomaly
+      nemo field name: sossheig
+      model space: surface
+    - name: sea_surface_height_anomaly_background_error
+      nemo field name: sossheig
+      model space: surface
+      variable type: background variance
+  grid name: Data/amm1_atlas_grid_spec.yaml
+  number levels: 1
+  partitioner: serial
+  source mesh halo: 1
+state :
+  date: 2021-06-29T23:00:00Z
+  state variables:
+    - sea_surface_height_anomaly
+    - sea_surface_height_anomaly_background_error
+  nemo field file: &field_file Data/amm1_nemo.nc
+  variance field file: *field_file
+observations:
+  observers:
+  - obs space:
+      name: Sea Surface Height
+      obsdatain:
+        engine:
+          type: H5File
+          obsfile: Data/amm_ssh_obs.nc
+      simulated variables: [seaSurfaceHeightAnomaly]
+    get values:
+      atlas-interpolator:
+        type: unstructured-bilinear-lonlat
+        non_linear: missing-if-all-missing-real32
+    obs operator:
+      name: Composite
+      components:
+      - name: Identity
+        observation alias file: testinput/test_name_map.yaml
+      - name: BackgroundErrorIdentity
+        observation alias file: testinput/test_name_map.yaml
+test:
+  reference filename: testoutput/test_hofx3d_nc_ssh_amm1.ref
+  float absolute tolerance: 5.0e-7 # tolerance for MPI sum-order related differences
diff --git a/src/tests/testinput/hofx_nc_ssh_amm1.yaml b/src/tests/testinput/hofx_nc_ssh_amm1.yaml
new file mode 100644
index 0000000..7d73d6e
--- /dev/null
+++ b/src/tests/testinput/hofx_nc_ssh_amm1.yaml
@@ -0,0 +1,71 @@
+# (C) British Crown Copyright 2024 Met Office
+
+forecast length : P2D
+time window:
+  begin: 2021-06-28T23:00:00Z
+  length: P2D
+geometry :
+  nemo variables:
+    - name: sea_surface_height_anomaly
+      nemo field name: sossheig
+      model space: surface
+    - name: sea_surface_height_anomaly_background_error
+      nemo field name: sossheig
+      model space: surface
+      variable type: background variance
+  grid name: Data/amm1_atlas_grid_spec.yaml
+  number levels: 1
+  partitioner: serial
+  source mesh halo: 1
+initial condition :
+  date: 2021-06-28T23:00:00Z
+  state variables:
+    - sea_surface_height_anomaly
+    - sea_surface_height_anomaly_background_error
+  nemo field file: &field_file Data/amm1_nemo.nc
+  variance field file: *field_file
+model :
+  name: PseudoModel
+  tstep: P1D
+  states:
+    - date: 2021-06-29T23:00:00Z
+      nemo field file: *field_file
+      variance field file: *field_file
+      state variables: &state_variables [ sea_surface_height_anomaly ]
+    - date: 2021-06-30T23:00:00Z
+      nemo field file: *field_file
+      variance field file: *field_file
+      state variables: *state_variables
+observations:
+  observers:
+  - obs space:
+      name: Sea Surface Height
+      obsdatain:
+        engine:
+          type: H5File
+          obsfile: Data/amm_ssh_obs.nc
+      obsdataout:
+        engine:
+          type: H5File
+          obsfile: testoutput/hofx_nc_ssh_jopa.nc
+      simulated variables: [seaSurfaceHeightAnomaly]
+    get values:
+      atlas-interpolator:
+        type: unstructured-bilinear-lonlat
+        non_linear: missing-if-all-missing-real32
+    obs operator:
+      name: Composite
+      components:
+      - name: Identity
+        observation alias file: testinput/test_name_map.yaml
+      - name: BackgroundErrorIdentity
+        observation alias file: testinput/test_name_map.yaml
+    obs filters:
+    - filter: Variable Assignment
+      assignments:
+      - name: GrossErrorProbability/seaSurfaceHeightAnomaly
+        type: float
+        value: 0.04
+test:
+  reference filename: testoutput/test_hofx_nc_ssh_amm1.ref
+  float absolute tolerance: 5.0e-7 # tolerance for MPI sum-order related differences
diff --git a/src/tests/testoutput/CMakeLists.txt b/src/tests/testoutput/CMakeLists.txt
index 17e60a4..f0dc157 100644
--- a/src/tests/testoutput/CMakeLists.txt
+++ b/src/tests/testoutput/CMakeLists.txt
@@ -5,9 +5,11 @@ list( APPEND orcajedi_test_output
   test_hofx_nc_sst.ref
   test_hofx_nc_ssh.ref
   test_hofx_nc_ssh_eorca025.ref
+  test_hofx_nc_ssh_amm1.ref
   test_hofx3d_nc_sst.ref
   test_hofx3d_nc_potm.ref
   test_hofx3d_nc_prof_2vars.ref
+  test_hofx3d_nc_ssh_amm1.ref
 )
 
 foreach(FILENAME ${orcajedi_test_output})
diff --git a/src/tests/testoutput/test_hofx3d_nc_ssh_amm1.ref b/src/tests/testoutput/test_hofx3d_nc_ssh_amm1.ref
new file mode 100644
index 0000000..df9f439
--- /dev/null
+++ b/src/tests/testoutput/test_hofx3d_nc_ssh_amm1.ref
@@ -0,0 +1,12 @@
+Test     : State: 
+Test     :  Model state valid at time: 2021-06-29T23:00:00Z
+Test     :     2 variables: sea_surface_height_anomaly, sea_surface_height_anomaly_background_error
+Test     :     atlas field norms:
+Test     :         sea_surface_height_anomaly: 3.35877e-02
+Test     :         sea_surface_height_anomaly_background_error: 3.35877e-02
+
+Test     : H(x): 
+Test     : Sea Surface Height nobs= 19 Min=-4.86354e-01, Max=-1.43070e-01, RMS=4.07939e-01
+
+Test     : End H(x)
+
diff --git a/src/tests/testoutput/test_hofx_nc_ssh_amm1.ref b/src/tests/testoutput/test_hofx_nc_ssh_amm1.ref
new file mode 100644
index 0000000..c07a05b
--- /dev/null
+++ b/src/tests/testoutput/test_hofx_nc_ssh_amm1.ref
@@ -0,0 +1,18 @@
+Test     : Initial state:
+Test     :  Model state valid at time: 2021-06-28T23:00:00Z
+Test     :     2 variables: sea_surface_height_anomaly, sea_surface_height_anomaly_background_error
+Test     :     atlas field norms:
+Test     :         sea_surface_height_anomaly: 3.35877e-02
+Test     :         sea_surface_height_anomaly_background_error: 3.35877e-02
+
+Test     : Final state:
+Test     :  Model state valid at time: 2021-06-30T23:00:00Z
+Test     :     2 variables: sea_surface_height_anomaly, sea_surface_height_anomaly_background_error
+Test     :     atlas field norms:
+Test     :         sea_surface_height_anomaly: 3.35877e-02
+Test     :         sea_surface_height_anomaly_background_error: 3.35877e-02
+
+Test     : H(x):
+Test     : Sea Surface Height nobs= 19 Min=-4.86354e-01, Max=-1.43070e-01, RMS=4.07939e-01
+
+Test     : End H(x)