diff --git a/include/bout/adios_object.hxx b/include/bout/adios_object.hxx
index 63a2d4f85f..ee7e82e45a 100755
--- a/include/bout/adios_object.hxx
+++ b/include/bout/adios_object.hxx
@@ -16,11 +16,23 @@
 
 #if BOUT_HAS_ADIOS2
 
+#include "bout/assert.hxx"
+#include "bout/boutcomm.hxx"
 #include "bout/boutexception.hxx"
+#include "bout/field2d.hxx"
+#include "bout/field3d.hxx"
+#include "bout/fieldperp.hxx"
+#include "bout/globals.hxx"
+#include "bout/mesh.hxx"
+#include "bout/utils.hxx"
 
 #include <adios2.h>
 #include <memory>
 #include <mpi.h>
+#include <string>
+#include <tuple>
+#include <utility>
+#include <vector>
 
 namespace bout {
 
@@ -29,21 +41,17 @@ void ADIOSInit(const std::string configFile, MPI_Comm comm);
 void ADIOSFinalize();
 
 using ADIOSPtr = std::shared_ptr<adios2::ADIOS>;
-using EnginePtr = std::shared_ptr<adios2::Engine>;
-using IOPtr = std::shared_ptr<adios2::IO>;
 
 ADIOSPtr GetADIOSPtr();
-IOPtr GetIOPtr(const std::string IOName);
 
 class ADIOSStream {
 public:
   adios2::IO io;
-  adios2::Variable<double> vTime;
-  adios2::Variable<int> vStep;
   int adiosStep = 0;
 
   /** create or return the ADIOSStream based on the target file name */
-  static ADIOSStream& ADIOSGetStream(const std::string& fname, adios2::Mode mode);
+  static ADIOSStream& ADIOSGetStream(const std::string& fname, adios2::Mode mode,
+                                     const std::string& engineType = "BP5");
 
   ~ADIOSStream();
 
@@ -56,6 +64,91 @@ public:
     return v;
   }
 
+  template <class T>
+  void Get(const std::string& varname, T& value, adios2::Mode mode = adios2::Mode::Sync) {
+    auto variable = io.InquireVariable<T>(varname);
+    ASSERT1(variable);
+    ASSERT1(variable.ShapeID() == adios2::ShapeID::GlobalValue);
+    engine().Get(variable, &value, mode);
+  }
+
+  template <class T>
+  void Get(const std::string& varname, Array<T>& value,
+           adios2::Mode mode = adios2::Mode::Sync) {
+    GetArrayLike(varname, value, mode);
+  }
+
+  template <class T>
+  void Get(const std::string& varname, Matrix<T>& value,
+           adios2::Mode mode = adios2::Mode::Sync) {
+    GetArrayLike(varname, value, mode);
+  }
+
+  template <class T>
+  void Get(const std::string& varname, Tensor<T>& value,
+           adios2::Mode mode = adios2::Mode::Sync) {
+    GetArrayLike(varname, value, mode);
+  }
+
+  void Get(const std::string& varname, Field2D& value,
+           adios2::Mode mode = adios2::Mode::Sync) {
+    value.allocate();
+    GetField(varname, {"x", "y"}, *value.getMesh(), &value(0, 0), mode);
+  }
+
+  void Get(const std::string& varname, Field3D& value,
+           adios2::Mode mode = adios2::Mode::Sync) {
+    value.allocate();
+    GetField(varname, {"x", "y", "z"}, *value.getMesh(), &value(0, 0, 0), mode);
+  }
+
+  void Get(const std::string& varname, FieldPerp& value,
+           adios2::Mode mode = adios2::Mode::Sync) {
+    value.allocate();
+    GetField(varname, {"x", "z"}, *value.getMesh(), &value(0, 0), mode);
+  }
+
+  template <class T>
+  void Put(const std::string& varname, T value, adios2::Mode mode = adios2::Mode::Sync) {
+    if (BoutComm::rank() != 0) {
+      return;
+    }
+    engine().Put(GetValueVariable<T>(varname), value, mode);
+  }
+
+  template <class T>
+  void Put(const std::string& varname, const Array<T>& value,
+           adios2::Mode mode = adios2::Mode::Sync) {
+    PutArrayLike(varname, value, mode);
+  }
+
+  template <class T>
+  void Put(const std::string& varname, const Matrix<T>& value,
+           adios2::Mode mode = adios2::Mode::Sync) {
+    PutArrayLike(varname, value, mode);
+  }
+
+  template <class T>
+  void Put(const std::string& varname, const Tensor<T>& value,
+           adios2::Mode mode = adios2::Mode::Sync) {
+    PutArrayLike(varname, value, mode);
+  }
+
+  void Put(const std::string& varname, const Field2D& value,
+           adios2::Mode mode = adios2::Mode::Sync) {
+    PutField(varname, {"x", "y"}, *value.getMesh(), &value(0, 0), mode);
+  }
+
+  void Put(const std::string& varname, const Field3D& value,
+           adios2::Mode mode = adios2::Mode::Sync) {
+    PutField(varname, {"x", "y", "z"}, *value.getMesh(), &value(0, 0, 0), mode);
+  }
+
+  void Put(const std::string& varname, const FieldPerp& value,
+           adios2::Mode mode = adios2::Mode::Sync) {
+    PutField(varname, {"x", "z"}, *value.getMesh(), &value(0, 0), mode);
+  }
+
   template <class T>
   adios2::Variable<T>
   GetArrayVariable(const std::string& varname, const adios2::Dims& shape,
@@ -78,7 +171,7 @@ public:
     if (not engine_) {
       engine_ = io.Open(fname, file_mode);
       if (not engine_) {
-        throw BoutException("Could not open ADIOS file '{:s}' for writing", fname);
+        throw BoutException("Could not open ADIOS file '{:s}'", fname);
       }
     }
     return engine_;
@@ -101,25 +194,184 @@ public:
 
   void finish() {
     if (engine_) {
-      engine().EndStep();
+      if (isInStep) {
+        engine().EndStep();
+        isInStep = false;
+      }
       engine().Close();
+      engine_ = adios2::Engine();
     }
   }
 
 private:
-  ADIOSStream(const std::string& fname, adios2::Mode mode)
-      : fname(fname), file_mode(mode) {
-
-    ADIOSPtr adiosp = GetADIOSPtr();
-    std::string ioname = "write_" + fname;
-    try {
-      io = adiosp->AtIO(ioname);
-    } catch (const std::invalid_argument& e) {
-      io = adiosp->DeclareIO(ioname);
-      io.SetEngine("BP5");
+  ADIOSStream(const std::string& fname, adios2::Mode mode, const std::string& engineType);
+
+  struct FieldSelection {
+    adios2::Dims shape;
+    adios2::Dims start;
+    adios2::Dims count;
+    adios2::Dims mem_start;
+    adios2::Dims mem_count;
+
+    auto selection() const { return adios2::Box<adios2::Dims>{start, count}; }
+    auto memorySelection() const {
+      return adios2::Box<adios2::Dims>{mem_start, mem_count};
     }
   };
 
+  void GetField(const std::string& varname, const std::vector<std::string>& dim_names,
+                const Mesh& mesh, BoutReal* data, adios2::Mode mode) {
+    auto variable = io.InquireVariable<BoutReal>(varname);
+    ASSERT1(variable);
+    ASSERT1(variable.ShapeID() == adios2::ShapeID::GlobalArray);
+
+    auto selection = makeFieldSelection(dim_names, mesh);
+    variable.SetSelection(selection.selection());
+    variable.SetMemorySelection(selection.memorySelection());
+    engine().Get(variable, data, mode);
+  }
+
+  void PutField(const std::string& varname, const std::vector<std::string>& dim_names,
+                const Mesh& mesh, const BoutReal* data, adios2::Mode mode) {
+    auto selection = makeFieldSelection(dim_names, mesh);
+    auto variable =
+        GetArrayVariable<BoutReal>(varname, selection.shape, dim_names, BoutComm::rank());
+    variable.SetSelection(selection.selection());
+    variable.SetMemorySelection(selection.memorySelection());
+    engine().Put(variable, data, mode);
+  }
+
+  FieldSelection makeFieldSelection(const std::vector<std::string>& dim_names,
+                                    const Mesh& mesh) const {
+    ASSERT1(!dim_names.empty());
+    ASSERT1(dim_names.size() <= 3);
+    ASSERT1(dim_names[0] == "x");
+
+    FieldSelection selection;
+    selection.shape.push_back(static_cast<std::size_t>(mesh.GlobalNx));
+    selection.start.push_back(static_cast<std::size_t>(mesh.MapGlobalX));
+    selection.count.push_back(static_cast<std::size_t>(mesh.MapCountX));
+    selection.mem_start.push_back(static_cast<std::size_t>(mesh.MapLocalX));
+    selection.mem_count.push_back(static_cast<std::size_t>(mesh.LocalNx));
+
+    if (dim_names.size() > 1) {
+      if (dim_names[1] == "y") {
+        selection.shape.push_back(static_cast<std::size_t>(mesh.GlobalNy));
+        selection.start.push_back(static_cast<std::size_t>(mesh.MapGlobalY));
+        selection.count.push_back(static_cast<std::size_t>(mesh.MapCountY));
+        selection.mem_start.push_back(static_cast<std::size_t>(mesh.MapLocalY));
+        selection.mem_count.push_back(static_cast<std::size_t>(mesh.LocalNy));
+      } else if (dim_names[1] == "z") {
+        selection.shape.push_back(static_cast<std::size_t>(mesh.GlobalNz));
+        selection.start.push_back(static_cast<std::size_t>(mesh.MapGlobalZ));
+        selection.count.push_back(static_cast<std::size_t>(mesh.MapCountZ));
+        selection.mem_start.push_back(static_cast<std::size_t>(mesh.MapLocalZ));
+        selection.mem_count.push_back(static_cast<std::size_t>(mesh.LocalNz));
+      } else {
+        ASSERT1(false);
+      }
+    }
+
+    if (dim_names.size() > 2) {
+      ASSERT1(dim_names[1] == "y");
+      ASSERT1(dim_names[2] == "z");
+      selection.shape.push_back(static_cast<std::size_t>(mesh.GlobalNz));
+      selection.start.push_back(static_cast<std::size_t>(mesh.MapGlobalZ));
+      selection.count.push_back(static_cast<std::size_t>(mesh.MapCountZ));
+      selection.mem_start.push_back(static_cast<std::size_t>(mesh.MapLocalZ));
+      selection.mem_count.push_back(static_cast<std::size_t>(mesh.LocalNz));
+    }
+
+    return selection;
+  }
+
+  template <class Container>
+  void GetArrayLike(const std::string& varname, Container& value, adios2::Mode mode) {
+    using T = typename Container::data_type;
+    auto variable = io.InquireVariable<T>(varname);
+    ASSERT1(variable);
+    ASSERT1(variable.ShapeID() == adios2::ShapeID::GlobalArray);
+
+    const auto shape = variable.Shape();
+    auto dims_attr = io.InquireAttribute<std::string>(varname + "/__xarray_dimensions__");
+    auto read_shape = shape;
+
+    if (dims_attr) {
+      const auto dim_names = dims_attr.Data();
+      if (!dim_names.empty() && dim_names[0] == "rank") {
+        ASSERT1(!shape.empty());
+        ASSERT1(static_cast<std::size_t>(BoutComm::rank()) < shape[0]);
+
+        adios2::Dims start{static_cast<std::size_t>(BoutComm::rank())};
+        adios2::Dims count{1};
+        for (std::size_t i = 1; i < shape.size(); i++) {
+          start.push_back(0);
+          count.push_back(shape[i]);
+        }
+
+        variable.SetSelection(adios2::Box<adios2::Dims>{start, count});
+        variable.SetMemorySelection(adios2::Box<adios2::Dims>{start, count});
+        read_shape = adios2::Dims(shape.begin() + 1, shape.end());
+      } else if (!dim_names.empty() && dim_names[0] == "x") {
+        ASSERT1(globals::mesh);
+        auto selection = makeFieldSelection(dim_names, *globals::mesh);
+        variable.SetSelection(selection.selection());
+        variable.SetMemorySelection(selection.memorySelection());
+        read_shape = selection.mem_count;
+      }
+    }
+
+    constexpr auto ndims = std::tuple_size_v<decltype(value.shape())>;
+    ASSERT1(read_shape.size() == ndims);
+
+    resizeForShape(value, read_shape, std::make_index_sequence<ndims>{});
+    engine().Get(variable, value.begin(), mode);
+  }
+
+  template <class Container>
+  void PutArrayLike(const std::string& varname, const Container& value,
+                    adios2::Mode mode) {
+    using T = typename Container::data_type;
+    auto var = GetArrayVariable<T>(varname, makeShape(BoutComm::size(), value),
+                                   makeDimNames(value), BoutComm::rank());
+    var.SetSelection(adios2::Box<adios2::Dims>{makeStart(value), makeShape(1, value)});
+    engine().Put<T>(var, value.begin(), mode);
+  }
+
+  template <class Container>
+  adios2::Dims makeShape(std::size_t first, const Container& value) const {
+    return std::apply(
+        [first](auto... sizes) {
+          return adios2::Dims{first, static_cast<std::size_t>(sizes)...};
+        },
+        value.shape());
+  }
+
+  template <class Container>
+  adios2::Dims makeStart(const Container& value) const {
+    constexpr auto ndims = std::tuple_size_v<decltype(value.shape())>;
+    adios2::Dims start(ndims + 1, 0);
+    start[0] = static_cast<std::size_t>(BoutComm::rank());
+    return start;
+  }
+
+  template <class Container>
+  std::vector<std::string> makeDimNames(const Container& value) const {
+    constexpr auto ndims = std::tuple_size_v<decltype(value.shape())>;
+    std::vector<std::string> dim_names{"rank"};
+    dim_names.reserve(ndims + 1);
+    for (std::size_t i = 0; i < ndims; i++) {
+      dim_names.push_back("dim_" + std::to_string(i));
+    }
+    return dim_names;
+  }
+
+  template <class Container, std::size_t... I>
+  void resizeForShape(Container& value, const adios2::Dims& shape,
+                      std::index_sequence<I...> /*indices*/) {
+    value.reallocate(static_cast<typename Container::size_type>(shape[I])...);
+  }
+
   std::string fname;
   adios2::Mode file_mode;
   adios2::Engine engine_;
@@ -128,10 +380,6 @@ private:
   bool isInStep = false;
 };
 
-/** Set user parameters for an IO group */
-void ADIOSSetParameters(const std::string& input, char delimKeyValue, char delimItem,
-                        adios2::IO& io);
-
 } // namespace bout
 
 #endif //BOUT_HAS_ADIOS2
diff --git a/src/sys/adios_object.cxx b/src/sys/adios_object.cxx
index 1b3e0a12e2..596be66b93 100644
--- a/src/sys/adios_object.cxx
+++ b/src/sys/adios_object.cxx
@@ -7,6 +7,7 @@
 
 #include <adios2.h>
 
+#include <string>
 #include <unordered_map>
 
 namespace bout {
@@ -14,6 +15,30 @@ namespace bout {
 static ADIOSPtr adios = nullptr;
 static std::unordered_map<std::string, ADIOSStream> adiosStreams;
 
+namespace {
+std::string GetIOName(const std::string& fname, adios2::Mode mode) {
+  switch (mode) {
+  case adios2::Mode::Read:
+    return "read_" + fname;
+  case adios2::Mode::ReadRandomAccess:
+    return "read_random_access_" + fname;
+  case adios2::Mode::Append:
+    return "append_" + fname;
+  case adios2::Mode::Write:
+    return "write_" + fname;
+  default:
+    return fname;
+  }
+}
+
+adios2::IO GetIO(const std::string& fname, adios2::Mode mode,
+                 const std::string& engineType) {
+  auto io = GetADIOSPtr()->DeclareIO(GetIOName(fname, mode));
+  io.SetEngine(engineType);
+  return io;
+}
+} // namespace
+
 void ADIOSInit(MPI_Comm comm) { adios = std::make_shared<adios2::ADIOS>(comm); }
 
 void ADIOSInit(const std::string configFile, MPI_Comm comm) {
@@ -37,16 +62,6 @@ ADIOSPtr GetADIOSPtr() {
   return adios;
 }
 
-IOPtr GetIOPtr(const std::string IOName) {
-  auto adios = GetADIOSPtr();
-  IOPtr io = nullptr;
-  try {
-    io = std::make_shared<adios2::IO>(adios->AtIO(IOName));
-  } catch (std::invalid_argument& e) {
-  }
-  return io;
-}
-
 ADIOSStream::~ADIOSStream() {
   if (engine_) {
     if (isInStep) {
@@ -57,42 +72,19 @@ ADIOSStream::~ADIOSStream() {
   }
 }
 
-ADIOSStream& ADIOSStream::ADIOSGetStream(const std::string& fname, adios2::Mode mode) {
-  auto it = adiosStreams.find(fname);
+ADIOSStream::ADIOSStream(const std::string& fname, adios2::Mode mode,
+                         const std::string& engineType)
+    : io(GetIO(fname, mode, engineType)), fname(fname), file_mode(mode) {}
+
+ADIOSStream& ADIOSStream::ADIOSGetStream(const std::string& fname, adios2::Mode mode,
+                                         const std::string& engineType) {
+  const auto key = GetIOName(fname, mode);
+  auto it = adiosStreams.find(key);
   if (it == adiosStreams.end()) {
-    it = adiosStreams.emplace(fname, ADIOSStream(fname, mode)).first;
+    it = adiosStreams.emplace(key, ADIOSStream(fname, mode, engineType)).first;
   }
   return it->second;
 }
 
-void ADIOSSetParameters(const std::string& input, char delimKeyValue, char delimItem,
-                        adios2::IO& io) {
-  auto lf_Trim = [](std::string& input) {
-    input.erase(0, input.find_first_not_of(" \n\r\t")); // prefixing spaces
-    input.erase(input.find_last_not_of(" \n\r\t") + 1); // suffixing spaces
-  };
-
-  std::istringstream inputSS(input);
-  std::string parameter;
-  while (std::getline(inputSS, parameter, delimItem)) {
-    const size_t position = parameter.find(delimKeyValue);
-    if (position == std::string::npos) {
-      throw BoutException("ADIOSSetParameters(): wrong format for IO parameter "
-                          + parameter + ", format must be key" + delimKeyValue
-                          + "value for each entry");
-    }
-
-    std::string key = parameter.substr(0, position);
-    lf_Trim(key);
-    std::string value = parameter.substr(position + 1);
-    lf_Trim(value);
-    if (value.length() == 0) {
-      throw BoutException("ADIOS2SetParameters: empty value in IO parameter " + parameter
-                          + ", format must be key" + delimKeyValue + "value");
-    }
-    io.SetParameter(key, value);
-  }
-}
-
 } // namespace bout
 #endif //BOUT_HAS_ADIOS2
diff --git a/src/sys/options/options_adios.cxx b/src/sys/options/options_adios.cxx
index 844284a400..4abf55bcd5 100644
--- a/src/sys/options/options_adios.cxx
+++ b/src/sys/options/options_adios.cxx
@@ -10,157 +10,41 @@
 #include "bout/bout_types.hxx"
 #include "bout/boutexception.hxx"
 #include "bout/field2d.hxx"
-#include "bout/globals.hxx"
-#include "bout/mesh.hxx"
 #include "bout/options.hxx"
 #include "bout/options_io.hxx"
 #include "bout/output.hxx"
 #include "bout/sys/timer.hxx"
 #include "bout/sys/variant.hxx"
-#include "bout/traits.hxx"
 #include "bout/utils.hxx"
 
 #include <adios2.h> // IWYU pragma: keep
 #include <fmt/format.h>
 #include <fmt/ranges.h>
 
-#include <algorithm>
 #include <cstddef>
 #include <cstdint>
 #include <exception>
-#include <iterator>
 #include <stdexcept>
 #include <string>
-#include <utility>
-#include <vector>
 
 namespace {
-auto to_int_dims(const adios2::Dims& dims) {
-  std::vector<int> int_dims;
-  int_dims.reserve(dims.size());
-  std::transform(dims.begin(), dims.end(), std::back_inserter(int_dims),
-                 [](auto dim) { return static_cast<int>(dim); });
-  return int_dims;
-}
-
-// Helper class to construct Adios hyperslices
-struct Selection {
-  // Offset of this processor's data into the global array
-  adios2::Dims start;
-  // The size of the mapped region
-  adios2::Dims count;
-  // Where the actual data starts in data pointer (to exclude ghost cells)
-  adios2::Dims mem_start;
-  // The actual size of data pointer in memory (including ghost cells)
-  adios2::Dims mem_count;
-  // Global shape, including boundaries but not guard cells
-  adios2::Dims shape;
-  // Shape of the local variable to read into
-  std::vector<int> dims;
-
-  // Distributed Field/Array/Matrix/Tensor
-  bool should_set_selection{false};
-
-  Selection(const std::vector<std::string>& dim_names, const std::vector<int>& dim_sizes,
-            const Mesh& mesh) {
-    const auto ndims = dim_names.size();
-    const bool dim0_is_rank = ndims > 0 ? (dim_names[0] == "rank") : false;
-    const bool dim0_is_x = ndims > 0 ? (dim_names[0] == "x") : false;
-    const bool dim1_is_y = ndims > 1 ? (dim_names[1] == "y") : false;
-    const bool dim1_is_z = ndims > 1 ? (dim_names[1] == "z") : false;
-    const bool dim2_is_z = ndims > 2 ? (dim_names[2] == "z") : false;
-
-    should_set_selection = dim0_is_rank or (ndims == 2 and dim0_is_x and dim1_is_y)
-                           or (ndims == 2 and dim0_is_x and dim1_is_z)
-                           or (ndims == 3 and dim0_is_x and dim1_is_y and dim2_is_z);
-
-    if (dim0_is_rank) {
-      const auto ndim_sizes = dim_sizes.size();
-      ASSERT3(ndim_sizes > 1);
-
-      // This is a distributed array, so the local variable is going
-      // to be shape (dim_sizes[1]...) (that is, drop the rank)
-      dims.push_back(dim_sizes[1]);
-      // but we tell adios to read our rank's bit with the full ndims
-      start = {static_cast<std::size_t>(BoutComm::rank()), 0};
-      count = {std::size_t{1}, static_cast<std::size_t>(dim_sizes[0])};
-
-      if (ndim_sizes > 2) {
-        dims.push_back(dim_sizes[2]);
-        start.push_back(0);
-        count.push_back(dim_sizes[1]);
-      }
-      if (ndim_sizes > 3) {
-        dims.push_back(dim_sizes[3]);
-        start.push_back(0);
-        count.push_back(dim_sizes[2]);
-      }
-      mem_count = count;
-      mem_start = start;
-      return;
-    }
-
-    if (ndims > 0) {
-      dims.push_back(dim0_is_x ? mesh.LocalNx : dim_sizes[0]);
-    }
-    if (ndims > 1) {
-      if (dim1_is_y) {
-        dims.push_back(mesh.LocalNy);
-      } else if (dim1_is_z) {
-        dims.push_back(mesh.LocalNz);
-      } else {
-        dims.push_back(dim_sizes[1]);
-      }
-    }
-    if (ndims > 2) {
-      dims.push_back(dim2_is_z ? mesh.LocalNz : dim_sizes[2]);
-    }
-
-    shape.push_back(static_cast<std::size_t>(mesh.GlobalNx));
-    start.push_back(static_cast<std::size_t>(mesh.MapGlobalX));
-    count.push_back(static_cast<std::size_t>(mesh.MapCountX));
-    mem_start.push_back(static_cast<std::size_t>(mesh.MapLocalX));
-    mem_count.push_back(static_cast<std::size_t>(mesh.LocalNx));
-
-    if (dim1_is_y) {
-      shape.push_back(static_cast<std::size_t>(mesh.GlobalNy));
-      start.push_back(static_cast<std::size_t>(mesh.MapGlobalY));
-      count.push_back(static_cast<std::size_t>(mesh.MapCountY));
-      mem_start.push_back(static_cast<std::size_t>(mesh.MapLocalY));
-      mem_count.push_back(static_cast<std::size_t>(mesh.LocalNy));
-    } else if (dim1_is_z) {
-      shape.push_back(static_cast<std::size_t>(mesh.GlobalNz));
-      start.push_back(static_cast<std::size_t>(mesh.MapGlobalZ));
-      count.push_back(static_cast<std::size_t>(mesh.MapCountZ));
-      mem_start.push_back(static_cast<std::size_t>(mesh.MapLocalZ));
-      mem_count.push_back(static_cast<std::size_t>(mesh.LocalNz));
-    }
-
-    if (dim2_is_z) {
-      shape.push_back(static_cast<std::size_t>(mesh.GlobalNz));
-      start.push_back(static_cast<std::size_t>(mesh.MapGlobalZ));
-      count.push_back(static_cast<std::size_t>(mesh.MapCountZ));
-      mem_start.push_back(static_cast<std::size_t>(mesh.MapLocalZ));
-      mem_count.push_back(static_cast<std::size_t>(mesh.LocalNz));
+std::size_t getLocalNDims(adios2::IO& io, const std::string& name,
+                          const adios2::Dims& shape) {
+  auto dims_attr = io.InquireAttribute<std::string>(name + "/__xarray_dimensions__");
+  if (dims_attr) {
+    const auto dim_names = dims_attr.Data();
+    if (!dim_names.empty() && dim_names[0] == "rank") {
+      ASSERT1(!shape.empty());
+      return shape.size() - 1;
     }
   }
+  return shape.size();
+}
 
-  auto selection() const { return adios2::Box<adios2::Dims>{start, count}; }
-  auto memorySelection() const { return adios2::Box<adios2::Dims>{mem_start, mem_count}; }
-};
-
-template <template <class> class T, class U>
-auto read_variable(adios2::IO& io, adios2::Engine& reader, const std::string& name,
-                   const Selection& selection, T<U>& value) {
-  auto variable = io.InquireVariable<U>(name);
-
-  if (selection.should_set_selection) {
-    variable.SetSelection(selection.selection());
-    variable.SetMemorySelection(selection.memorySelection());
-  }
-
+template <class T>
+Options readArrayVariable(bout::ADIOSStream& stream, const std::string& name, T& value) {
   try {
-    reader.Get<U>(variable, value.begin(), adios2::Mode::Sync);
+    stream.Get(name, value);
   } catch (const std::exception& e) {
     output_warn.write("ADIOS exception while reading '{}': {}\n", name, e.what());
     return Options{};
@@ -169,84 +53,77 @@ auto read_variable(adios2::IO& io, adios2::Engine& reader, const std::string& na
 }
 
 template <class T>
-Options readVariable(adios2::Engine& reader, adios2::IO& io, const std::string& name,
+Options readVariable(bout::ADIOSStream& stream, const std::string& name,
                      const std::string& type) {
-  auto variable = io.InquireVariable<T>(name);
-
-  using bout::globals::mesh;
+  auto variable = stream.io.InquireVariable<T>(name);
 
   if (variable.ShapeID() == adios2::ShapeID::GlobalValue) {
     T value;
-    reader.Get<T>(variable, &value, adios2::Mode::Sync);
+    stream.Get(name, value);
     return Options(value);
   }
 
   if (variable.ShapeID() == adios2::ShapeID::LocalArray) {
     throw BoutException("ADIOS reader did not implement reading local arrays like `{}` "
                         "'{}' in file '{}'\n",
-                        type, name, reader.Name());
+                        type, name, stream.engine().Name());
   }
 
   if (type != "double" && type != "float" && type != "int32_t" && type != "int64_t") {
     throw BoutException("ADIOS reader did not implement reading arrays that are not "
                         "`double`/`float`/`int32_t`/`int64_t` type. "
                         "Found `{}` '{}' in file '{}'\n",
-                        type, name, reader.Name());
+                        type, name, stream.engine().Name());
   }
 
   if (type == "double" && sizeof(BoutReal) != sizeof(double)) {
     throw BoutException(
         "ADIOS does not allow for implicit type conversions. BoutReal type is "
         "float but found `{}` '{}' in file '{}'\n",
-        type, name, reader.Name());
+        type, name, stream.engine().Name());
   }
 
   if (type == "float" && sizeof(BoutReal) != sizeof(float)) {
     throw BoutException(
         "ADIOS reader does not allow for implicit type conversions. BoutReal type is "
         "double but found `{}` '{}' in file '{}'",
-        type, name, reader.Name());
+        type, name, stream.engine().Name());
   }
 
-  const auto dims = to_int_dims(variable.Shape());
-  auto dims_attr =
-      io.InquireAttribute<std::string>(fmt::format("{}/__xarray_dimensions__", name))
-          .Data();
-
-  const auto selection = Selection(dims_attr, dims, *mesh);
-  const auto ndims = selection.dims.size();
+  const auto shape = variable.Shape();
+  const auto ndims = getLocalNDims(stream.io, name, shape);
 
   switch (ndims) {
   case 1: {
     if (type == "float" or type == "double") {
-      Array<BoutReal> value(dims[0]);
-      return read_variable(io, reader, name, selection, value);
+      Array<BoutReal> value;
+      return readArrayVariable(stream, name, value);
     }
     if (type == "int32_t" or type == "int64_t") {
-      Array<int> value(dims[0]);
-      return read_variable(io, reader, name, selection, value);
+      Array<int> value;
+      return readArrayVariable(stream, name, value);
     }
     break;
   }
   case 2: {
     if (type == "float" or type == "double") {
-      Matrix<BoutReal> value(selection.dims[0], selection.dims[1]);
-      return read_variable(io, reader, name, selection, value);
+      Matrix<BoutReal> value;
+      return readArrayVariable(stream, name, value);
     }
     if (type == "int32_t" or type == "int64_t") {
-      Matrix<int> value(selection.dims[0], selection.dims[1]);
-      return read_variable(io, reader, name, selection, value);
+      Matrix<int> value;
+      return readArrayVariable(stream, name, value);
     }
     break;
   }
   case 3: {
     if (type == "float" or type == "double") {
-      Tensor<BoutReal> value(selection.dims[0], selection.dims[1], selection.dims[2]);
-      return read_variable(io, reader, name, selection, value);
+      Tensor<BoutReal> value;
+      return readArrayVariable(stream, name, value);
     }
     if (type == "int32_t" or type == "int64_t") {
-      Tensor<int> value(selection.dims[0], selection.dims[1], selection.dims[2]);
-      return read_variable(io, reader, name, selection, value);
+      Tensor<int> value;
+      return readArrayVariable(stream, name, value);
     }
     break;
   }
@@ -254,52 +131,52 @@ Options readVariable(adios2::Engine& reader, adios2::IO& io, const std::string&
   default:
     break;
   }
-  auto dims_str = fmt::format("[{}]", fmt::join(dims, ", "));
+  auto dims_str = fmt::format("[{}]", fmt::join(shape, ", "));
   throw BoutException(
       "ADIOS reader failed to read '{}' (shape: {}, type: '{}') in file '{}'\n", name,
-      dims_str, type, reader.Name());
+      dims_str, type, stream.engine().Name());
 }
 
-Options readVariable(adios2::Engine& reader, adios2::IO& io, const std::string& name,
+Options readVariable(bout::ADIOSStream& stream, const std::string& name,
                      const std::string& type) {
   if (type == adios2::GetType<std::string>()) {
-    return readVariable<std::string>(reader, io, name, type);
+    return readVariable<std::string>(stream, name, type);
   }
   if (type == adios2::GetType<char>()) {
-    return readVariable<char>(reader, io, name, type);
+    return readVariable<char>(stream, name, type);
   }
   if (type == adios2::GetType<int8_t>()) {
-    return readVariable<int8_t>(reader, io, name, type);
+    return readVariable<int8_t>(stream, name, type);
   }
   if (type == adios2::GetType<int16_t>()) {
-    return readVariable<int16_t>(reader, io, name, type);
+    return readVariable<int16_t>(stream, name, type);
   }
   if (type == adios2::GetType<int32_t>()) {
-    return readVariable<int32_t>(reader, io, name, type);
+    return readVariable<int32_t>(stream, name, type);
   }
   if (type == adios2::GetType<int64_t>()) {
-    return readVariable<int64_t>(reader, io, name, type);
+    return readVariable<int64_t>(stream, name, type);
   }
   if (type == adios2::GetType<uint8_t>()) {
-    return readVariable<uint8_t>(reader, io, name, type);
+    return readVariable<uint8_t>(stream, name, type);
   }
   if (type == adios2::GetType<uint16_t>()) {
-    return readVariable<uint16_t>(reader, io, name, type);
+    return readVariable<uint16_t>(stream, name, type);
   }
   if (type == adios2::GetType<uint32_t>()) {
-    return readVariable<uint32_t>(reader, io, name, type);
+    return readVariable<uint32_t>(stream, name, type);
   }
   if (type == adios2::GetType<uint64_t>()) {
-    return readVariable<uint64_t>(reader, io, name, type);
+    return readVariable<uint64_t>(stream, name, type);
   }
   if (type == adios2::GetType<float>()) {
-    return readVariable<float>(reader, io, name, type);
+    return readVariable<float>(stream, name, type);
   }
   if (type == adios2::GetType<double>()) {
-    return readVariable<double>(reader, io, name, type);
+    return readVariable<double>(stream, name, type);
   }
   if (type == adios2::GetType<long double>()) {
-    return readVariable<long double>(reader, io, name, type);
+    return readVariable<long double>(stream, name, type);
   }
 
   output_warn.write("ADIOS readVariable can't read type '{}' (variable '{}')\n", type,
@@ -356,16 +233,9 @@ Options OptionsADIOS::read([[maybe_unused]] bool lazy) {
   const Timer timer("io");
 
   // Open file
-  const ADIOSPtr adiosp = GetADIOSPtr();
-  adios2::IO io;
-  const std::string ioname = "read_" + filename;
-  try {
-    io = adiosp->AtIO(ioname);
-  } catch (const std::invalid_argument& e) {
-    io = adiosp->DeclareIO(ioname);
-  }
-
-  adios2::Engine reader = io.Open(filename, adios2::Mode::ReadRandomAccess);
+  ADIOSStream& stream =
+      ADIOSStream::ADIOSGetStream(filename, adios2::Mode::ReadRandomAccess);
+  adios2::Engine& reader = stream.engine();
   if (!reader) {
     throw BoutException("Could not open ADIOS file '{:s}' for reading\n", filename);
   }
@@ -373,7 +243,7 @@ Options OptionsADIOS::read([[maybe_unused]] bool lazy) {
   Options result;
 
   // Iterate over all variables
-  for (const auto& varpair : io.AvailableVariables()) {
+  for (const auto& varpair : stream.io.AvailableVariables()) {
     const auto& var_name = varpair.first; // Name of the variable
 
     auto it = varpair.second.find("Type");
@@ -388,21 +258,21 @@ Options OptionsADIOS::read([[maybe_unused]] bool lazy) {
     // Note: Copying the value rather than simple assignment is used
     // because the Options assignment operator overwrites full_name.
     var = 0; // Setting is_section to false
-    var.value = readVariable(reader, io, var_name, var_type).value;
+    var.value = readVariable(stream, var_name, var_type).value;
     var.attributes["source"] = filename;
 
     // Get variable attributes
-    for (const auto& attpair : io.AvailableAttributes(var_name, "/", true)) {
+    for (const auto& attpair : stream.io.AvailableAttributes(var_name, "/", true)) {
       const auto& att_name = attpair.first; // Attribute name
       const auto& att = attpair.second;     // attribute params
 
       auto it = att.find("Type");
       const std::string& att_type = it->second;
-      readAttribute(io, att_name, att_type, var);
+      readAttribute(stream.io, att_name, att_type, var);
     }
   }
 
-  reader.Close();
+  stream.finish();
 
   return result;
 }
@@ -421,142 +291,22 @@ void OptionsADIOS::verifyTimesteps() const {
 }
 } // namespace bout
 
-const std::vector<std::string> DIMS_X = {"x"};
-const std::vector<std::string> DIMS_XY = {"x", "y"};
-const std::vector<std::string> DIMS_XZ = {"x", "z"};
-const std::vector<std::string> DIMS_XYZ = {"x", "y", "z"};
-
 namespace {
-using bout::utils::tuple_index_sequence;
-
-template <class Tuple, std::size_t... I>
-auto make_shape_impl(std::size_t first, const Tuple& t,
-                     std::index_sequence<I...> /* index */) {
-  return adios2::Dims{first, static_cast<std::size_t>(std::get<I>(t))...};
-}
-// Return an `adios2::Dims` with value ``{first, value.shape()[0]...}``
-template <class T>
-auto make_shape(std::size_t first, const T& value) {
-  const auto shape = value.shape();
-  return make_shape_impl(first, shape, tuple_index_sequence<decltype(shape)>{});
-}
-
-template <class Tuple, std::size_t... I>
-auto make_start_impl(const Tuple& t, std::index_sequence<I...> /* index */) {
-  // Hey look, a legitimate use of the comma operator to get a bunch
-  // of zeros the length of the index_sequence!
-  return adios2::Dims{static_cast<std::size_t>(BoutComm::rank()),
-                      (std::get<I>(t), std::size_t{0})...};
-}
-// Return an `adios2::Dims` with value ``{rank, 0...}``, with as many zeros as the dimension of ``T``
-template <class T>
-auto make_start(const T& value) {
-  const auto shape = value.shape();
-  return make_start_impl(shape, tuple_index_sequence<decltype(shape)>{});
-}
-
-template <std::size_t... I>
-auto make_dims_impl(std::index_sequence<I...> /*index*/) {
-  using namespace std::string_literals;
-  return std::vector{"rank"s, (fmt::format("dim_{}", I))...};
-}
-// Return vector of dimension names: `{"rank", "dim_0", ...}`
-template <class T>
-auto make_dims(const T& value) {
-  return make_dims_impl(tuple_index_sequence<decltype(value.shape())>{});
-}
-
 /// Visit a variant type, and put the data into a NcVar
 struct ADIOSPutVarVisitor {
   ADIOSPutVarVisitor(const std::string& name, bout::ADIOSStream& stream)
       : varname(name), stream(stream) {}
   template <typename T>
   void operator()(const T& value) {
-    adios2::Variable<T> var = stream.GetValueVariable<T>(varname);
-    stream.engine().Put(var, value);
-  }
-
-  void operator()(const Array<int>& value) { amt_put_helper(value); }
-  void operator()(const Array<BoutReal>& value) { amt_put_helper(value); }
-  void operator()(const Matrix<int>& value) { amt_put_helper(value); }
-  void operator()(const Matrix<BoutReal>& value) { amt_put_helper(value); }
-  void operator()(const Tensor<int>& value) { amt_put_helper(value); }
-  void operator()(const Tensor<BoutReal>& value) { amt_put_helper(value); }
-  void operator()(bool value) {
-    // Scalars are only written from processor 0
-    if (BoutComm::rank() != 0) {
-      return;
-    }
-    stream.engine().Put(stream.GetValueVariable<int>(varname), static_cast<int>(value));
-  }
-  void operator()(int value) {
-    // Scalars are only written from processor 0
-    if (BoutComm::rank() != 0) {
-      return;
-    }
-    stream.engine().Put(stream.GetValueVariable<int>(varname), value);
-  }
-  void operator()(BoutReal value) {
-    // Scalars are only written from processor 0
-    if (BoutComm::rank() != 0) {
-      return;
-    }
-    stream.engine().Put(stream.GetValueVariable<BoutReal>(varname), value);
-  }
-
-  void operator()(const std::string& value) {
-    // Scalars are only written from processor 0
-    if (BoutComm::rank() != 0) {
-      return;
-    }
-    stream.engine().Put<std::string>(stream.GetValueVariable<std::string>(varname), value,
-                                     adios2::Mode::Sync);
+    stream.Put(varname, value);
   }
 
-  void operator()(const Field2D& value) {
-    // Empty dim_sizes is fine because we provide full list of dim_names
-    auto selection = Selection(DIMS_XY, {}, *value.getMesh());
-    auto var = stream.GetArrayVariable<BoutReal>(varname, selection.shape, DIMS_XY,
-                                                 BoutComm::rank());
-    var.SetSelection(selection.selection());
-    var.SetMemorySelection(selection.memorySelection());
-    stream.engine().Put(var, &value(0, 0));
-  }
-
-  void operator()(const Field3D& value) {
-    // Empty dim_sizes is fine because we provide full list of dim_names
-    auto selection = Selection(DIMS_XYZ, {}, *value.getMesh());
-    auto var = stream.GetArrayVariable<BoutReal>(varname, selection.shape, DIMS_XYZ,
-                                                 BoutComm::rank());
-    var.SetSelection(selection.selection());
-    var.SetMemorySelection(selection.memorySelection());
-    stream.engine().Put(var, &value(0, 0, 0));
-  }
-
-  void operator()(const FieldPerp& value) {
-    // Empty dim_sizes is fine because we provide full list of dim_names
-    auto selection = Selection(DIMS_XZ, {}, *value.getMesh());
-    auto var = stream.GetArrayVariable<BoutReal>(varname, selection.shape, DIMS_XZ,
-                                                 BoutComm::rank());
-    var.SetSelection(selection.selection());
-    var.SetMemorySelection(selection.memorySelection());
-    stream.engine().Put<BoutReal>(var, &value(0, 0));
-  }
+  void operator()(bool value) { stream.Put(varname, static_cast<int>(value)); }
 
 private:
   // Ok to keep const/refs here as visitor is only a temporary
   const std::string& varname; // NOLINT(*-avoid-const-or-ref-data-members)
   bout::ADIOSStream& stream;  // NOLINT(*-avoid-const-or-ref-data-members)
-
-  // helper for `Array`, `Matrix`, `Tensor`
-  template <template <class> class T, class U>
-  void amt_put_helper(const T<U>& value) {
-    const auto shape = make_shape(static_cast<std::size_t>(BoutComm::size()), value);
-    auto var =
-        stream.GetArrayVariable<U>(varname, shape, make_dims(value), BoutComm::rank());
-    var.SetSelection({make_start(value), make_shape(1, value)});
-    stream.engine().Put<U>(var, value.begin());
-  }
 };
 
 /// Visit a variant type, and put the data into a NcVar
diff --git a/tests/unit/sys/test_options_adios2.cxx b/tests/unit/sys/test_options_adios2.cxx
index 13bbe3db54..a0fcf20b2e 100644
--- a/tests/unit/sys/test_options_adios2.cxx
+++ b/tests/unit/sys/test_options_adios2.cxx
@@ -27,7 +27,8 @@ class OptionsAdios2Test : public FakeMeshFixture {
   static void TearDownTestSuite() { bout::ADIOSFinalize(); }
 
   ~OptionsAdios2Test() override {
-    bout::ADIOSStream::ADIOSGetStream(filename, adios2::Mode::Deferred).engine().Close();
+    bout::ADIOSStream::ADIOSGetStream(filename, adios2::Mode::Write).finish();
+    bout::ADIOSStream::ADIOSGetStream(filename, adios2::Mode::Append).finish();
   }
 
   // A temporary filename