diff --git a/include/bout/adios_object.hxx b/include/bout/adios_object.hxx
index 12da8ce4c7..07482a0dcb 100755
--- a/include/bout/adios_object.hxx
+++ b/include/bout/adios_object.hxx
@@ -57,9 +57,9 @@ public:
   }
 
   template <class T>
-  adios2::Variable<T> GetArrayVariable(const std::string& varname, adios2::Dims& shape,
-                                       const std::vector<std::string>& dimNames,
-                                       int rank) {
+  adios2::Variable<T>
+  GetArrayVariable(const std::string& varname, const adios2::Dims& shape,
+                   const std::vector<std::string>& dimNames, int rank) {
     adios2::Variable<T> v = io.InquireVariable<T>(varname);
     if (!v) {
       adios2::Dims start(shape.size());
diff --git a/include/bout/array.hxx b/include/bout/array.hxx
index 2c42f15aad..2305fe3b22 100644
--- a/include/bout/array.hxx
+++ b/include/bout/array.hxx
@@ -29,6 +29,7 @@
 #include <algorithm>
 #include <map>
 #include <memory>
+#include <tuple>
 #include <vector>
 
 #if BOUT_USE_OPENMP
@@ -283,6 +284,9 @@ public:
     return ptr->size();
   }
 
+  /// Return shape of the array (the `size()` in a length-1 tuple)
+  std::tuple<size_type> shape() const { return std::make_tuple(size()); };
+
   /*!
    * Returns true if the data is unique to this Array.
    *
diff --git a/include/bout/options.hxx b/include/bout/options.hxx
index dbb926838a..1c552cd4f9 100644
--- a/include/bout/options.hxx
+++ b/include/bout/options.hxx
@@ -12,7 +12,7 @@
 * options and allows access to all sub-sections
 *
 **************************************************************************
-* Copyright 2010-2024 BOUT++ contributors
+* Copyright 2010-2025 BOUT++ contributors
 *
 * Contact: Ben Dudson, dudson2@llnl.gov
 *
@@ -39,6 +39,7 @@ class Options;
 #ifndef OPTIONS_H
 #define OPTIONS_H
 
+#include "bout/array.hxx"
 #include "bout/bout_types.hxx"
 #include "bout/field2d.hxx"
 #include "bout/field3d.hxx"
@@ -203,7 +204,8 @@ public:
   /// The type used to store values
   using ValueType =
       bout::utils::variant<bool, int, BoutReal, std::string, Field2D, Field3D, FieldPerp,
-                           Array<BoutReal>, Matrix<BoutReal>, Tensor<BoutReal>>;
+                           Array<BoutReal>, Array<int>, Matrix<BoutReal>, Matrix<int>,
+                           Tensor<BoutReal>, Tensor<int>>;
 
   /// Methods to iterate over `Options`
   auto begin() { return std::begin(children); }
@@ -963,9 +965,15 @@ Options& Options::assign<>(FieldPerp val, std::string source);
 template <>
 Options& Options::assign<>(Array<BoutReal> val, std::string source);
 template <>
+Options& Options::assign<>(Array<int> val, std::string source);
+template <>
 Options& Options::assign<>(Matrix<BoutReal> val, std::string source);
 template <>
+Options& Options::assign<>(Matrix<int> val, std::string source);
+template <>
 Options& Options::assign<>(Tensor<BoutReal> val, std::string source);
+template <>
+Options& Options::assign<>(Tensor<int> val, std::string source);
 
 /// Specialised similar comparison methods
 template <>
@@ -975,25 +983,31 @@ inline bool Options::similar<BoutReal>(BoutReal lhs, BoutReal rhs) const {
 
 /// Specialised as routines
 template <>
-std::string Options::as<std::string>(const std::string& similar_to) const;
+auto Options::as(const std::string& similar_to) const -> std::string;
+template <>
+auto Options::as(const int& similar_to) const -> int;
+template <>
+auto Options::as(const BoutReal& similar_to) const -> BoutReal;
+template <>
+auto Options::as(const bool& similar_to) const -> bool;
 template <>
-int Options::as<int>(const int& similar_to) const;
+auto Options::as(const Field2D& similar_to) const -> Field2D;
 template <>
-BoutReal Options::as<BoutReal>(const BoutReal& similar_to) const;
+auto Options::as(const Field3D& similar_to) const -> Field3D;
 template <>
-bool Options::as<bool>(const bool& similar_to) const;
+auto Options::as(const FieldPerp& similar_to) const -> FieldPerp;
 template <>
-Field2D Options::as<Field2D>(const Field2D& similar_to) const;
+auto Options::as(const Array<BoutReal>& similar_to) const -> Array<BoutReal>;
 template <>
-Field3D Options::as<Field3D>(const Field3D& similar_to) const;
+auto Options::as(const Array<int>& similar_to) const -> Array<int>;
 template <>
-FieldPerp Options::as<FieldPerp>(const FieldPerp& similar_to) const;
+auto Options::as(const Matrix<BoutReal>& similar_to) const -> Matrix<BoutReal>;
 template <>
-Array<BoutReal> Options::as<Array<BoutReal>>(const Array<BoutReal>& similar_to) const;
+auto Options::as(const Matrix<int>& similar_to) const -> Matrix<int>;
 template <>
-Matrix<BoutReal> Options::as<Matrix<BoutReal>>(const Matrix<BoutReal>& similar_to) const;
+auto Options::as(const Tensor<BoutReal>& similar_to) const -> Tensor<BoutReal>;
 template <>
-Tensor<BoutReal> Options::as<Tensor<BoutReal>>(const Tensor<BoutReal>& similar_to) const;
+auto Options::as(const Tensor<int>& similar_to) const -> Tensor<int>;
 
 /// Convert \p value to string
 std::string toString(const Options& value);
diff --git a/include/bout/sys/type_name.hxx b/include/bout/sys/type_name.hxx
index be4d096b4d..fcd251bf7a 100644
--- a/include/bout/sys/type_name.hxx
+++ b/include/bout/sys/type_name.hxx
@@ -4,13 +4,19 @@
 #ifndef TYPE_NAME_HXX
 #define TYPE_NAME_HXX
 
+#include "bout/array.hxx"
 #include "bout/bout_types.hxx"
+
 #include <string>
 #include <typeinfo>
 
 class Field2D;
 class Field3D;
 class FieldPerp;
+template <class T>
+class Matrix;
+template <class T>
+class Tensor;
 
 namespace bout {
 namespace utils {
@@ -41,6 +47,19 @@ std::string typeName<Field3D>();
 
 template <>
 std::string typeName<FieldPerp>();
+
+template <>
+std::string typeName<Array<int>>();
+template <>
+std::string typeName<Array<BoutReal>>();
+template <>
+std::string typeName<Matrix<int>>();
+template <>
+std::string typeName<Matrix<BoutReal>>();
+template <>
+std::string typeName<Tensor<int>>();
+template <>
+std::string typeName<Tensor<BoutReal>>();
 } // namespace utils
 } // namespace bout
 
diff --git a/include/bout/traits.hxx b/include/bout/traits.hxx
index 4e03d6b526..5796c65879 100644
--- a/include/bout/traits.hxx
+++ b/include/bout/traits.hxx
@@ -1,6 +1,7 @@
 #ifndef BOUT_TRAITS_H
 #define BOUT_TRAITS_H
 
+#include <tuple>
 #include <type_traits>
 
 class Field;
@@ -112,6 +113,10 @@ using EnableIfFieldPerp =
 /// Enable a function if T is a subclass of Options
 template <class T>
 using EnableIfOptions = std::enable_if_t<std::is_base_of_v<Options, T>>;
+
+/// Create a `std::index_sequence` with the length of a tuple ``T``
+template <typename T>
+using tuple_index_sequence = std::make_index_sequence<std::tuple_size_v<T>>;
 } // namespace utils
 } // namespace bout
 
diff --git a/include/bout/utils.hxx b/include/bout/utils.hxx
index b8c5a5c944..f8d6f08856 100644
--- a/include/bout/utils.hxx
+++ b/include/bout/utils.hxx
@@ -5,9 +5,9 @@
  * simple but common calculations
  *
  **************************************************************************
- * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu
+ * Copyright 2010 - 2025 BOUT++ contributors
  *
- * Contact: Ben Dudson, bd512@york.ac.uk
+ * Contact: Ben Dudson, dudson2@llnl.gov
  * 
  * This file is part of BOUT++.
  *
@@ -517,18 +517,18 @@ std::string toString(const T& val) {
 /// where the type may be std::string.
 inline std::string toString(const std::string& val) { return val; }
 
-template <>
-inline std::string toString<>(const Array<BoutReal>& UNUSED(val)) {
+template <typename T>
+inline std::string toString(const Array<T>& UNUSED(val)) {
   return "<Array>";
 }
 
-template <>
-inline std::string toString<>(const Matrix<BoutReal>& UNUSED(val)) {
+template <typename T>
+inline std::string toString(const Matrix<T>& UNUSED(val)) {
   return "<Matrix>";
 }
 
-template <>
-inline std::string toString<>(const Tensor<BoutReal>& UNUSED(val)) {
+template <typename T>
+inline std::string toString(const Tensor<T>& UNUSED(val)) {
   return "<Tensor>";
 }
 
diff --git a/src/mesh/data/gridfromfile.cxx b/src/mesh/data/gridfromfile.cxx
index 5625522795..5078882136 100644
--- a/src/mesh/data/gridfromfile.cxx
+++ b/src/mesh/data/gridfromfile.cxx
@@ -2,6 +2,7 @@
 #include "bout/traits.hxx"
 #include <bout/griddata.hxx>
 
+#include <bout/array.hxx>
 #include <bout/boutexception.hxx>
 #include <bout/constants.hxx>
 #include <bout/fft.hxx>
@@ -11,6 +12,9 @@
 #include <bout/sys/timer.hxx>
 #include <bout/unused.hxx>
 #include <bout/utils.hxx>
+
+#include <algorithm>
+#include <iterator>
 #include <utility>
 
 GridFile::GridFile(std::string gridfilename)
@@ -471,7 +475,26 @@ bool GridFile::get([[maybe_unused]] Mesh* m, [[maybe_unused]] std::vector<int>&
                    [[maybe_unused]] GridDataSource::Direction dir) {
   TRACE("GridFile::get(vector<int>)");
 
-  return false;
+  if (not data.isSet(name)) {
+    return false;
+  }
+
+  const auto full_var = data[name].as<Array<int>>();
+
+  // Check size
+  if (full_var.size() < len + offset) {
+    throw BoutException("{} has length {}. Expected {} elements + {} offset", name,
+                        full_var.size(), len, offset);
+  }
+
+  // Ensure that output variable has the correct size
+  var.resize(len);
+
+  const auto* it = std::begin(full_var);
+  std::advance(it, offset);
+  std::copy_n(it, len, std::begin(var));
+
+  return true;
 }
 
 bool GridFile::get(Mesh* UNUSED(m), std::vector<BoutReal>& var, const std::string& name,
diff --git a/src/mesh/impls/bout/boutmesh.cxx b/src/mesh/impls/bout/boutmesh.cxx
index 574902ea7b..9f366d9456 100644
--- a/src/mesh/impls/bout/boutmesh.cxx
+++ b/src/mesh/impls/bout/boutmesh.cxx
@@ -2141,6 +2141,37 @@ void BoutMesh::topology() {
     add_target(ny_inner - 1, 0, nx);
   }
 
+  // Additional limiters
+  // Each limiter needs 3 indices: A Y index, start and end X indices
+  int limiter_count = 0;
+  Mesh::get(limiter_count, "limiter_count", 0);
+  if (limiter_count > 0) {
+    std::vector<int> limiter_yinds;
+    if (!source->get(this, limiter_yinds, "limiter_yinds", limiter_count)) {
+      throw BoutException("Couldn't read limiter_yinds vector of length {} from mesh",
+                          limiter_count);
+    }
+    std::vector<int> limiter_xstarts;
+    if (!source->get(this, limiter_xstarts, "limiter_xstarts", limiter_count)) {
+      throw BoutException("Couldn't read limiter_xstarts vector of length {} from mesh",
+                          limiter_count);
+    }
+    std::vector<int> limiter_xends;
+    if (!source->get(this, limiter_xends, "limiter_xends", limiter_count)) {
+      throw BoutException("Couldn't read limiter_xend vector of length {} from mesh",
+                          limiter_count);
+    }
+
+    for (int i = 0; i < limiter_count; ++i) {
+      int const yind = limiter_yinds[i];
+      int const xstart = limiter_xstarts[i];
+      int const xend = limiter_xends[i];
+      output_info.write("Adding a limiter between y={} and {}. X indices {} to {}\n",
+                        yind, yind + 1, xstart, xend);
+      add_target(yind, xstart, xend);
+    }
+  }
+
   if ((ixseps_inner > 0)
       && (((PE_YIND * MYSUB > jyseps1_1) && (PE_YIND * MYSUB <= jyseps2_1))
           || ((PE_YIND * MYSUB > jyseps1_2) && (PE_YIND * MYSUB <= jyseps2_2)))) {
diff --git a/src/sys/options.cxx b/src/sys/options.cxx
index 24c9e933c8..077b2a6d85 100644
--- a/src/sys/options.cxx
+++ b/src/sys/options.cxx
@@ -333,15 +333,30 @@ Options& Options::assign<>(Array<BoutReal> val, std::string source) {
   return *this;
 }
 template <>
+Options& Options::assign<>(Array<int> val, std::string source) {
+  _set_no_check(std::move(val), std::move(source));
+  return *this;
+}
+template <>
 Options& Options::assign<>(Matrix<BoutReal> val, std::string source) {
   _set_no_check(std::move(val), std::move(source));
   return *this;
 }
 template <>
+Options& Options::assign<>(Matrix<int> val, std::string source) {
+  _set_no_check(std::move(val), std::move(source));
+  return *this;
+}
+template <>
 Options& Options::assign<>(Tensor<BoutReal> val, std::string source) {
   _set_no_check(std::move(val), std::move(source));
   return *this;
 }
+template <>
+Options& Options::assign<>(Tensor<int> val, std::string source) {
+  _set_no_check(std::move(val), std::move(source));
+  return *this;
+}
 
 namespace {
 /// Use FieldFactory to evaluate expression
@@ -739,70 +754,69 @@ struct ConvertContainer {
 };
 } // namespace
 
-template <>
-Array<BoutReal> Options::as<Array<BoutReal>>(const Array<BoutReal>& similar_to) const {
-  if (is_section) {
-    throw BoutException(_("Option {:s} has no value"), full_name);
+namespace {
+// `Options::as` helper for `Array`, `Matrix`, `Tensor` (amt)
+template <class T>
+T as_amt(const Options& self, const T& similar_to) {
+  if (self.isSection()) {
+    throw BoutException(_("Option {:s} has no value"), self.str());
   }
 
-  Array<BoutReal> result = bout::utils::visit(
-      ConvertContainer<Array<BoutReal>>{
-          fmt::format(
-              _("Value for option {:s} cannot be converted to an Array<BoutReal>"),
-              full_name),
+  const T result = bout::utils::visit(
+      ConvertContainer<T>{
+          fmt::format(_("Value for option {:s} cannot be converted to an {}"), self.str(),
+                      bout::utils::typeName<T>()),
           similar_to},
-      value);
+      self.value);
 
-  // Mark this option as used
-  value_used = true;
+  printNameValueSourceLine(self, bout::utils::typeName<T>());
 
-  printNameValueSourceLine(*this, "Array<BoutReal>");
+  output_info.write("{} {} = {}", _("\tOption "), self.str(), bout::utils::typeName<T>());
+  if (self.hasAttribute("source")) {
+    // Specify the source of the setting
+    output_info.write(" ({})",
+                      bout::utils::variantToString(self.attributes.at("source")));
+  }
+  output_info.write("\n");
 
   return result;
 }
+} // namespace
 
 template <>
-Matrix<BoutReal> Options::as<Matrix<BoutReal>>(const Matrix<BoutReal>& similar_to) const {
-  if (is_section) {
-    throw BoutException(_("Option {:s} has no value"), full_name);
-  }
-
-  auto result = bout::utils::visit(
-      ConvertContainer<Matrix<BoutReal>>{
-          fmt::format(
-              _("Value for option {:s} cannot be converted to an Matrix<BoutReal>"),
-              full_name),
-          similar_to},
-      value);
-
-  // Mark this option as used
+auto Options::as(const Array<BoutReal>& similar_to) const -> Array<BoutReal> {
   value_used = true;
-
-  printNameValueSourceLine(*this, "Matrix<BoutReal>");
-
-  return result;
+  return as_amt(*this, similar_to);
 }
 
 template <>
-Tensor<BoutReal> Options::as<Tensor<BoutReal>>(const Tensor<BoutReal>& similar_to) const {
-  if (is_section) {
-    throw BoutException(_("Option {:s} has no value"), full_name);
-  }
+auto Options::as(const Array<int>& similar_to) const -> Array<int> {
+  value_used = true;
+  return as_amt(*this, similar_to);
+}
 
-  auto result = bout::utils::visit(
-      ConvertContainer<Tensor<BoutReal>>{
-          fmt::format(
-              _("Value for option {:s} cannot be converted to an Tensor<BoutReal>"),
-              full_name),
-          similar_to},
-      value);
+template <>
+auto Options::as(const Matrix<BoutReal>& similar_to) const -> Matrix<BoutReal> {
+  value_used = true;
+  return as_amt(*this, similar_to);
+}
 
-  // Mark this option as used
+template <>
+auto Options::as(const Matrix<int>& similar_to) const -> Matrix<int> {
   value_used = true;
+  return as_amt(*this, similar_to);
+}
 
-  printNameValueSourceLine(*this, "Tensor<BoutReal>");
+template <>
+auto Options::as(const Tensor<BoutReal>& similar_to) const -> Tensor<BoutReal> {
+  value_used = true;
+  return as_amt(*this, similar_to);
+}
 
-  return result;
+template <>
+auto Options::as(const Tensor<int>& similar_to) const -> Tensor<int> {
+  value_used = true;
+  return as_amt(*this, similar_to);
 }
 
 // Note: This is defined here rather than in the header
@@ -946,12 +960,17 @@ struct GetDimensions {
   std::vector<int> operator()([[maybe_unused]] int value) { return {1}; }
   std::vector<int> operator()([[maybe_unused]] BoutReal value) { return {1}; }
   std::vector<int> operator()([[maybe_unused]] const std::string& value) { return {1}; }
-  std::vector<int> operator()(const Array<BoutReal>& array) { return {array.size()}; }
-  std::vector<int> operator()(const Matrix<BoutReal>& array) {
+  template <typename T>
+  std::vector<int> operator()(const Array<T>& array) {
+    return {array.size()};
+  }
+  template <typename T>
+  std::vector<int> operator()(const Matrix<T>& array) {
     const auto shape = array.shape();
     return {std::get<0>(shape), std::get<1>(shape)};
   }
-  std::vector<int> operator()(const Tensor<BoutReal>& array) {
+  template <typename T>
+  std::vector<int> operator()(const Tensor<T>& array) {
     const auto shape = array.shape();
     return {std::get<0>(shape), std::get<1>(shape), std::get<2>(shape)};
   }
diff --git a/src/sys/options/options_adios.cxx b/src/sys/options/options_adios.cxx
index 09797647ed..4a41a014eb 100644
--- a/src/sys/options/options_adios.cxx
+++ b/src/sys/options/options_adios.cxx
@@ -1,4 +1,5 @@
 #include "bout/build_defines.hxx"
+#include "bout/traits.hxx"
 
 #if BOUT_HAS_ADIOS2
 
@@ -12,8 +13,11 @@
 #include "bout/sys/timer.hxx"
 
 #include "adios2.h"
+
+#include <cstddef>
 #include <exception>
-#include <iostream>
+#include <tuple>
+#include <utility>
 #include <vector>
 
 namespace bout {
@@ -313,6 +317,37 @@ 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, Tuple&& t,
+                     std::index_sequence<I...> /* index */) {
+  return adios2::Dims{first,
+                      static_cast<std::size_t>(std::get<I>(std::forward<Tuple>(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(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>(std::forward<Tuple>(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)>{});
+}
+} // namespace
+
 /// Visit a variant type, and put the data into a NcVar
 struct ADIOSPutVarVisitor {
   ADIOSPutVarVisitor(const std::string& name, ADIOSStream& stream)
@@ -323,9 +358,25 @@ struct ADIOSPutVarVisitor {
     stream.engine.Put(var, value);
   }
 
+  void operator()(const Array<int>& value) { put_helper(value); }
+  void operator()(const Array<BoutReal>& value) { put_helper(value); }
+  void operator()(const Matrix<int>& value) { put_helper(value); }
+  void operator()(const Matrix<BoutReal>& value) { put_helper(value); }
+  void operator()(const Tensor<int>& value) { put_helper(value); }
+  void operator()(const Tensor<BoutReal>& value) { put_helper(value); }
+
 private:
   const std::string& varname;
   ADIOSStream& stream;
+
+  // helper for `Array`, `Matrix`, `Tensor`
+  template <template <class> class T, class U>
+  void 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, DIMS_NONE, BoutComm::rank());
+    var.SetSelection({make_start(value), make_shape(1, value)});
+    stream.engine.Put<U>(var, value.begin());
+  }
 };
 
 template <>
@@ -472,46 +523,6 @@ void ADIOSPutVarVisitor::operator()<FieldPerp>(const FieldPerp& value) {
   stream.engine.Put<BoutReal>(var, &value(0, 0));
 }
 
-template <>
-void ADIOSPutVarVisitor::operator()<Array<BoutReal>>(const Array<BoutReal>& value) {
-  // Pointer to data. Assumed to be contiguous array
-  adios2::Dims shape = {(size_t)BoutComm::size(), (size_t)value.size()};
-  adios2::Dims start = {(size_t)BoutComm::rank(), 0};
-  adios2::Dims count = {1, shape[1]};
-  adios2::Variable<BoutReal> var =
-      stream.GetArrayVariable<BoutReal>(varname, shape, DIMS_NONE, BoutComm::rank());
-  var.SetSelection({start, count});
-  stream.engine.Put<BoutReal>(var, value.begin());
-}
-
-template <>
-void ADIOSPutVarVisitor::operator()<Matrix<BoutReal>>(const Matrix<BoutReal>& value) {
-  // Pointer to data. Assumed to be contiguous array
-  auto s = value.shape();
-  adios2::Dims shape = {(size_t)BoutComm::size(), (size_t)std::get<0>(s),
-                        (size_t)std::get<1>(s)};
-  adios2::Dims start = {(size_t)BoutComm::rank(), 0, 0};
-  adios2::Dims count = {1, shape[1], shape[2]};
-  adios2::Variable<BoutReal> var =
-      stream.GetArrayVariable<BoutReal>(varname, shape, DIMS_NONE, BoutComm::rank());
-  var.SetSelection({start, count});
-  stream.engine.Put<BoutReal>(var, value.begin());
-}
-
-template <>
-void ADIOSPutVarVisitor::operator()<Tensor<BoutReal>>(const Tensor<BoutReal>& value) {
-  // Pointer to data. Assumed to be contiguous array
-  auto s = value.shape();
-  adios2::Dims shape = {(size_t)BoutComm::size(), (size_t)std::get<0>(s),
-                        (size_t)std::get<1>(s), (size_t)std::get<2>(s)};
-  adios2::Dims start = {(size_t)BoutComm::rank(), 0, 0, 0};
-  adios2::Dims count = {1, shape[1], shape[2], shape[3]};
-  adios2::Variable<BoutReal> var =
-      stream.GetArrayVariable<BoutReal>(varname, shape, DIMS_NONE, BoutComm::rank());
-  var.SetSelection({start, count});
-  stream.engine.Put<BoutReal>(var, value.begin());
-}
-
 /// Visit a variant type, and put the data into a NcVar
 struct ADIOSPutAttVisitor {
   ADIOSPutAttVisitor(const std::string& varname, const std::string& attrname,
diff --git a/src/sys/options/options_netcdf.cxx b/src/sys/options/options_netcdf.cxx
index 94b62f0f02..f6edb26dbc 100644
--- a/src/sys/options/options_netcdf.cxx
+++ b/src/sys/options/options_netcdf.cxx
@@ -5,13 +5,15 @@
 #include "options_netcdf.hxx"
 
 #include "bout/bout.hxx"
-#include "bout/globals.hxx"
+#include "bout/bout_types.hxx"
 #include "bout/mesh.hxx"
 #include "bout/sys/timer.hxx"
+#include "bout/traits.hxx"
+
+#include <fmt/format.h>
 
 #include <climits>
 #include <exception>
-#include <iostream>
 #include <netcdf>
 #include <vector>
 
@@ -90,6 +92,10 @@ void readGroup(const std::string& filename, const NcGroup& group, Options& resul
         Array<double> value(static_cast<int>(dims[0].getSize()));
         var.getVar(value.begin());
         result[var_name] = value;
+      } else if (var_type == ncInt or var_type == ncShort) {
+        Array<int> value(static_cast<int>(dims[0].getSize()));
+        var.getVar(value.begin());
+        result[var_name] = value;
       } else if ((var_type == ncString) or (var_type == ncChar)) {
         std::string value;
         value.resize(dims[0].getSize());
@@ -104,6 +110,11 @@ void readGroup(const std::string& filename, const NcGroup& group, Options& resul
                              static_cast<int>(dims[1].getSize()));
         var.getVar(value.begin());
         result[var_name] = value;
+      } else if (var_type == ncInt) {
+        Matrix<int> value(static_cast<int>(dims[0].getSize()),
+                          static_cast<int>(dims[1].getSize()));
+        var.getVar(value.begin());
+        result[var_name] = value;
       }
       break;
     }
@@ -213,6 +224,13 @@ struct NcTypeVisitor {
   NcType operator()(const T& UNUSED(t)) {
     return {}; // Null object by default
   }
+
+  NcType operator()([[maybe_unused]] const Array<int>& t) { return ncInt; }
+  NcType operator()([[maybe_unused]] const Array<BoutReal>& t) { return ncDouble; }
+  NcType operator()([[maybe_unused]] const Matrix<int>& t) { return ncInt; }
+  NcType operator()([[maybe_unused]] const Matrix<BoutReal>& t) { return ncDouble; }
+  NcType operator()([[maybe_unused]] const Tensor<int>& t) { return ncInt; }
+  NcType operator()([[maybe_unused]] const Tensor<BoutReal>& t) { return ncDouble; }
 };
 
 template <>
@@ -255,20 +273,9 @@ NcType NcTypeVisitor::operator()<FieldPerp>(const FieldPerp& UNUSED(t)) {
   return operator()<BoutReal>(0.0);
 }
 
-/// Visit a variant type, returning dimensions
-struct NcDimVisitor {
-  NcDimVisitor(NcGroup& group) : group(group) {}
-  template <typename T>
-  std::vector<NcDim> operator()(const T& UNUSED(value)) {
-    return {};
-  }
-
-private:
-  NcGroup& group;
-};
-
+namespace {
+// Find dimension with name and size in group, creating it if it doesn't exist
 NcDim findDimension(NcGroup& group, const std::string& name, unsigned int size) {
-  // Get the dimension
   try {
     auto dim = group.getDim(name, NcGroup::ParentsAndCurrent);
     if (dim.isNull()) {
@@ -293,6 +300,40 @@ NcDim findDimension(NcGroup& group, const std::string& name, unsigned int size)
   }
 }
 
+using bout::utils::tuple_index_sequence;
+
+template <class Tuple, std::size_t... I>
+auto make_dims_impl(NcGroup& group, Tuple&& t, std::index_sequence<I...> /* index */) {
+  return std::vector{findDimension(group, fmt::format("dim_{}", I),
+                                   std::get<I>(std::forward<Tuple>(t)))...};
+}
+// Help get the dimensions for `Array`, `Matrix`, `Tensor` (amt)
+template <class T>
+auto make_amt_dims(NcGroup& group, const T& value) {
+  const auto shape = value.shape();
+  return make_dims_impl(group, shape, tuple_index_sequence<decltype(shape)>{});
+}
+} // namespace
+
+/// Visit a variant type, returning dimensions
+struct NcDimVisitor {
+  NcDimVisitor(NcGroup& group) : group(group) {}
+  template <typename T>
+  std::vector<NcDim> operator()(const T& UNUSED(value)) {
+    return {};
+  }
+
+  auto operator()(const Array<int>& value) { return make_amt_dims(group, value); }
+  auto operator()(const Array<BoutReal>& value) { return make_amt_dims(group, value); }
+  auto operator()(const Matrix<int>& value) { return make_amt_dims(group, value); }
+  auto operator()(const Matrix<BoutReal>& value) { return make_amt_dims(group, value); }
+  auto operator()(const Tensor<int>& value) { return make_amt_dims(group, value); }
+  auto operator()(const Tensor<BoutReal>& value) { return make_amt_dims(group, value); }
+
+private:
+  NcGroup& group;
+};
+
 template <>
 std::vector<NcDim> NcDimVisitor::operator()<Field2D>(const Field2D& value) {
   auto xdim = findDimension(group, "x", value.getNx());
@@ -337,6 +378,13 @@ struct NcPutVarVisitor {
     var.putVar(&value);
   }
 
+  void operator()(const Array<int>& value) { var.putVar(value.begin()); }
+  void operator()(const Array<BoutReal>& value) { var.putVar(value.begin()); }
+  void operator()(const Matrix<int>& value) { var.putVar(value.begin()); }
+  void operator()(const Matrix<BoutReal>& value) { var.putVar(value.begin()); }
+  void operator()(const Tensor<int>& value) { var.putVar(value.begin()); }
+  void operator()(const Tensor<BoutReal>& value) { var.putVar(value.begin()); }
+
 private:
   NcVar& var;
 };
@@ -383,6 +431,19 @@ struct NcPutVarCountVisitor {
     var.putVar(start, &value);
   }
 
+  void operator()(const Array<int>& value) { var.putVar(start, count, value.begin()); }
+  void operator()(const Array<BoutReal>& value) {
+    var.putVar(start, count, value.begin());
+  }
+  void operator()(const Matrix<int>& value) { var.putVar(start, count, value.begin()); }
+  void operator()(const Matrix<BoutReal>& value) {
+    var.putVar(start, count, value.begin());
+  }
+  void operator()(const Tensor<int>& value) { var.putVar(start, count, value.begin()); }
+  void operator()(const Tensor<BoutReal>& value) {
+    var.putVar(start, count, value.begin());
+  }
+
 private:
   NcVar& var;
   const std::vector<size_t>& start; ///< Starting (corner) index
diff --git a/src/sys/type_name.cxx b/src/sys/type_name.cxx
index c1082c499a..68272bd473 100644
--- a/src/sys/type_name.cxx
+++ b/src/sys/type_name.cxx
@@ -41,5 +41,31 @@ template <>
 std::string typeName<FieldPerp>() {
   return "FieldPerp";
 }
+
+template <>
+std::string typeName<Array<int>>() {
+  return "Array<int>";
+}
+template <>
+std::string typeName<Array<BoutReal>>() {
+  return "Array<BoutReal>";
+}
+template <>
+std::string typeName<Matrix<int>>() {
+  return "Matrix<int>";
+}
+template <>
+std::string typeName<Matrix<BoutReal>>() {
+  return "Matrix<BoutReal>";
+}
+template <>
+std::string typeName<Tensor<int>>() {
+  return "Tensor<int>";
+}
+template <>
+std::string typeName<Tensor<BoutReal>>() {
+  return "Tensor<BoutReal>";
+}
+
 } // namespace utils
 } // namespace bout
diff --git a/tests/integrated/test-options-adios/test-options-adios.cxx b/tests/integrated/test-options-adios/test-options-adios.cxx
index 60604e1aa3..ba0bbe3898 100644
--- a/tests/integrated/test-options-adios/test-options-adios.cxx
+++ b/tests/integrated/test-options-adios/test-options-adios.cxx
@@ -3,6 +3,7 @@
 
 #include "bout/options_io.hxx"
 #include "bout/optionsreader.hxx"
+#include "bout/utils.hxx"
 
 using bout::OptionsIO;
 
@@ -45,18 +46,31 @@ int main(int argc, char** argv) {
   ///////////////////////////
   // Write fields
 
-  Options fields;
-  fields["f2d"] = Field2D(1.0);
-  fields["f3d"] = Field3D(2.0);
-  fields["fperp"] = FieldPerp(3.0);
-  auto f = OptionsIO::create({{"file", "fields.bp"}, {"type", "adios"}});
-  /*
+  {
+    constexpr int nx = 2;
+    constexpr int ny = 3;
+    Matrix<int> matrix_in(nx, ny);
+    int count = 0;
+
+    for (int i = 0; i < nx; ++i) {
+      for (int j = 0; j < ny; ++j) {
+        matrix_in(i, j) = ++count;
+      }
+    }
+
+    Options fields{{"f2d", Field2D(1.0)},
+                   {"f3d", Field3D(2.0)},
+                   {"fperp", FieldPerp(3.0)},
+                   {"matrix", matrix_in}};
+    auto f = OptionsIO::create({{"file", "fields.bp"}, {"type", "adios"}});
+    /*
      write() for adios only buffers data but does not guarantee writing to disk
      unless singleWriteFile is set to true 
      */
-  f->write(fields);
-  // indicate completion of step, required to get data on disk
-  f->verifyTimesteps();
+    f->write(fields);
+    // indicate completion of step, required to get data on disk
+    f->verifyTimesteps();
+  }
 
   ///////////////////////////
   // Read fields
diff --git a/tests/unit/sys/test_options.cxx b/tests/unit/sys/test_options.cxx
index 67e405953d..b74cfdcb9f 100644
--- a/tests/unit/sys/test_options.cxx
+++ b/tests/unit/sys/test_options.cxx
@@ -1450,3 +1450,21 @@ TEST_F(OptionsTest, Iterate) {
   ASSERT_EQ(option["value1"], 2);
   ASSERT_EQ(option["value2"], 3);
 }
+
+TEST_F(OptionsTest, MatrixInt) {
+  constexpr int nx = 2;
+  constexpr int ny = 3;
+  Matrix<int> matrix_in(nx, ny);
+  int count = 0;
+
+  for (int i = 0; i < nx; ++i) {
+    for (int j = 0; j < ny; ++j) {
+      matrix_in(i, j) = ++count;
+    }
+  }
+
+  Options options{{"int_matrix", matrix_in}};
+
+  const auto matrix_out = options["int_matrix"].as<Matrix<int>>();
+  ASSERT_EQ(matrix_out(nx - 1, ny - 1), matrix_in(nx - 1, ny - 1));
+}
diff --git a/tests/unit/sys/test_options_netcdf.cxx b/tests/unit/sys/test_options_netcdf.cxx
index 127874b2ad..b16d838b54 100644
--- a/tests/unit/sys/test_options_netcdf.cxx
+++ b/tests/unit/sys/test_options_netcdf.cxx
@@ -1,6 +1,7 @@
 // Test reading and writing to NetCDF
 
 #include "bout/build_defines.hxx"
+#include "bout/utils.hxx"
 
 #if BOUT_HAS_NETCDF && !BOUT_HAS_LEGACY_NETCDF
 
@@ -93,6 +94,29 @@ TEST_F(OptionsNetCDFTest, ReadWriteField3D) {
   EXPECT_DOUBLE_EQ(value(1, 1, 1), 2.4);
 }
 
+TEST_F(OptionsNetCDFTest, ReadWriteMatrixInt) {
+  constexpr int nx = 2;
+  constexpr int ny = 3;
+  Matrix<int> matrix_in(nx, ny);
+  int count = 0;
+
+  for (int i = 0; i < nx; ++i) {
+    for (int j = 0; j < ny; ++j) {
+      matrix_in(i, j) = ++count;
+    }
+  }
+
+  {
+    Options options{{"int_matrix", matrix_in}};
+    OptionsIO::create(filename)->write(options);
+  }
+
+  Options data = OptionsIO::create(filename)->read();
+
+  const auto matrix_out = data["int_matrix"].as<Matrix<int>>();
+  ASSERT_EQ(matrix_out(nx - 1, ny - 1), matrix_in(nx - 1, ny - 1));
+}
+
 TEST_F(OptionsNetCDFTest, Groups) {
   {
     Options options;