diff --git a/lib/lattice/Lattice_transfer.h b/lib/lattice/Lattice_transfer.h
index 44f0337d..f988f310 100644
--- a/lib/lattice/Lattice_transfer.h
+++ b/lib/lattice/Lattice_transfer.h
@@ -599,6 +599,51 @@ unvectorizeToLexOrdArray(std::vector<sobj> &out, const Lattice<vobj> &in)
     extract1(in_vobj, out_ptrs, 0);
   }
 }
+
+template<typename vobj, typename sobj>
+typename std::enable_if<isSIMDvectorized<vobj>::value && !isSIMDvectorized<sobj>::value, void>::type 
+unvectorizeToRevLexOrdArray(std::vector<sobj> &out, const Lattice<vobj> &in)
+{
+
+  typedef typename vobj::vector_type vtype;
+  
+  GridBase* in_grid = in._grid;
+  out.resize(in_grid->lSites());
+  
+  int ndim = in_grid->Nd();
+  int in_nsimd = vtype::Nsimd();
+
+  std::vector<std::vector<int> > in_icoor(in_nsimd);
+      
+  for(int lane=0; lane < in_nsimd; lane++){
+    in_icoor[lane].resize(ndim);
+    in_grid->iCoorFromIindex(in_icoor[lane], lane);
+  }
+  
+  parallel_for(int in_oidx = 0; in_oidx < in_grid->oSites(); in_oidx++){ //loop over outer index
+    //Assemble vector of pointers to output elements
+    std::vector<sobj*> out_ptrs(in_nsimd);
+
+    std::vector<int> in_ocoor(ndim);
+    in_grid->oCoorFromOindex(in_ocoor, in_oidx);
+
+    std::vector<int> lcoor(in_grid->Nd());
+      
+    for(int lane=0; lane < in_nsimd; lane++){
+      for(int mu=0;mu<ndim;mu++)
+	lcoor[mu] = in_ocoor[mu] + in_grid->_rdimensions[mu]*in_icoor[lane][mu];
+
+      int lex;
+      Lexicographic::IndexFromCoorReversed(lcoor, lex, in_grid->_ldimensions);
+      out_ptrs[lane] = &out[lex];
+    }
+    
+    //Unpack into those ptrs
+    const vobj & in_vobj = in._odata[in_oidx];
+    extract1(in_vobj, out_ptrs, 0);
+  }
+}
+
 //Copy SIMD-vectorized lattice to array of scalar objects in lexicographic order
 template<typename vobj, typename sobj>
 typename std::enable_if<isSIMDvectorized<vobj>::value 
@@ -648,6 +693,54 @@ vectorizeFromLexOrdArray( std::vector<sobj> &in, Lattice<vobj> &out)
   }
 }
 
+template<typename vobj, typename sobj>
+typename std::enable_if<isSIMDvectorized<vobj>::value 
+                    && !isSIMDvectorized<sobj>::value, void>::type 
+vectorizeFromRevLexOrdArray( std::vector<sobj> &in, Lattice<vobj> &out)
+{
+
+  typedef typename vobj::vector_type vtype;
+  
+  GridBase* grid = out._grid;
+  assert(in.size()==grid->lSites());
+  
+  int ndim     = grid->Nd();
+  int nsimd    = vtype::Nsimd();
+
+  std::vector<std::vector<int> > icoor(nsimd);
+      
+  for(int lane=0; lane < nsimd; lane++){
+    icoor[lane].resize(ndim);
+    grid->iCoorFromIindex(icoor[lane],lane);
+  }
+  
+  parallel_for(uint64_t oidx = 0; oidx < grid->oSites(); oidx++){ //loop over outer index
+    //Assemble vector of pointers to output elements
+    std::vector<sobj*> ptrs(nsimd);
+
+    std::vector<int> ocoor(ndim);
+    grid->oCoorFromOindex(ocoor, oidx);
+
+    std::vector<int> lcoor(grid->Nd());
+      
+    for(int lane=0; lane < nsimd; lane++){
+
+      for(int mu=0;mu<ndim;mu++){
+	lcoor[mu] = ocoor[mu] + grid->_rdimensions[mu]*icoor[lane][mu];
+      }
+
+      int lex;
+      Lexicographic::IndexFromCoorReversed(lcoor, lex, grid->_ldimensions);
+      ptrs[lane] = &in[lex];
+    }
+    
+    //pack from those ptrs
+    vobj vecobj;
+    merge1(vecobj, ptrs, 0);
+    out._odata[oidx] = vecobj; 
+  }
+}
+
 //Convert a Lattice from one precision to another
 template<class VobjOut, class VobjIn>
 void precisionChange(Lattice<VobjOut> &out, const Lattice<VobjIn> &in){
diff --git a/lib/serialisation/Hdf5IO.h b/lib/serialisation/Hdf5IO.h
index 12625ab8..940bb11a 100644
--- a/lib/serialisation/Hdf5IO.h
+++ b/lib/serialisation/Hdf5IO.h
@@ -22,6 +22,8 @@
 
 namespace Grid
 {
+  template<class vobj> class Lattice;
+
   class Hdf5Writer: public Writer<Hdf5Writer>
   {
   public:
@@ -33,6 +35,8 @@ namespace Grid
     template <typename U>
     void writeDefault(const std::string &s, const U &x);
     template <typename U>
+    void writeDefault(const std::string &s, const Lattice<U> &field);
+    template <typename U>
     typename std::enable_if<element<std::vector<U>>::is_number, void>::type
     writeDefault(const std::string &s, const std::vector<U> &x);
     template <typename U>
@@ -60,6 +64,8 @@ namespace Grid
     template <typename U>
     void readDefault(const std::string &s, U &output);
     template <typename U>
+    void readDefault(const std::string &s, Lattice<U> &field);
+    template <typename U>
     typename std::enable_if<element<std::vector<U>>::is_number, void>::type
     readDefault(const std::string &s, std::vector<U> &x);
     template <typename U>
@@ -98,7 +104,40 @@ namespace Grid
   
   template <>
   void Hdf5Writer::writeDefault(const std::string &s, const std::string &x);
-  
+
+
+  template <typename U>
+  void Hdf5Writer::writeDefault(const std::string &s, const Lattice<U> &field)
+  {
+    // alias scalar types
+    typedef std::vector<typename U::scalar_object> ScalarLattice;
+    typedef typename U::scalar_type                ScalarType;
+
+    ScalarLattice scalField;
+
+    unvectorizeToRevLexOrdArray(scalField, field);
+
+    std::vector<hsize_t> dim;
+    std::vector<size_t>  tDim;
+
+    tensorDim(tDim, scalField[0]);
+    for (auto &d: field._grid->GlobalDimensions())
+    {
+      dim.push_back(d);
+    }
+    for (auto &d: tDim)
+    {
+      dim.push_back(d);
+    }
+
+    // write to file
+    H5NS::DataSpace dataSpace(dim.size(), dim.data());
+    H5NS::DataSet dataSet;
+      
+    dataSet = group_.createDataSet(s, Hdf5Type<ScalarType>::type(), dataSpace);
+    dataSet.write(scalField.data(), Hdf5Type<ScalarType>::type());
+  }
+
   template <typename U>
   typename std::enable_if<element<std::vector<U>>::is_number, void>::type
   Hdf5Writer::writeDefault(const std::string &s, const std::vector<U> &x)
@@ -169,6 +208,29 @@ namespace Grid
   template <>
   void Hdf5Reader::readDefault(const std::string &s, std::string &x);
   
+  template <typename U>
+  void Hdf5Reader::readDefault(const std::string &s, Lattice<U> &field)
+  {
+    // alias scalar types
+    typedef std::vector<typename U::scalar_object> ScalarLattice;
+    typedef typename U::scalar_type                ScalarType;
+
+    ScalarLattice        scalField;
+    H5NS::DataSet        dataSet;
+    std::vector<hsize_t> dim;
+    hsize_t              size = 1;
+   
+    dataSet = group_.openDataSet(s);
+    for (auto &d: field._grid->GlobalDimensions())
+    {
+      dim.push_back(d);
+      size *= d;
+    }
+    scalField.resize(size);
+    dataSet.read(scalField.data(), Hdf5Type<ScalarType>::type());
+    vectorizeFromRevLexOrdArray(scalField, field);
+  }
+
   template <typename U>
   typename std::enable_if<element<std::vector<U>>::is_number, void>::type
   Hdf5Reader::readDefault(const std::string &s, std::vector<U> &x)
diff --git a/lib/serialisation/VectorUtils.h b/lib/serialisation/VectorUtils.h
index 6df9416d..e1095488 100644
--- a/lib/serialisation/VectorUtils.h
+++ b/lib/serialisation/VectorUtils.h
@@ -2,7 +2,7 @@
 #define GRID_SERIALISATION_VECTORUTILS_H
 
 #include <type_traits>
-#include <Grid/tensors/Tensors.h>
+#include <Grid/tensors/Tensors.h>         
 
 namespace Grid {
   // Pair IO utilities /////////////////////////////////////////////////////////
@@ -78,6 +78,48 @@ namespace Grid {
     typedef typename std::vector<std::vector<typename TensorToVec<T>::type>> type;
   };
 
+  template <typename T>
+  void tensorDim(std::vector<size_t> &dim, const T &t, const bool wipe = true)
+  {
+    if (wipe)
+    {
+      dim.clear();
+    }
+  }
+
+  template <typename T>
+  void tensorDim(std::vector<size_t> &dim, const iScalar<T> &t, const bool wipe = true)
+  {
+    if (wipe)
+    {
+      dim.clear();
+    }
+    tensorDim(dim, t._internal, false);
+  }
+
+  template <typename T, int N>
+  void tensorDim(std::vector<size_t> &dim, const iVector<T, N> &t, const bool wipe = true)
+  {
+    if (wipe)
+    {
+      dim.clear();
+    }
+    dim.push_back(N);
+    tensorDim(dim, t._internal[0], false);
+  }
+
+  template <typename T, int N>
+  void tensorDim(std::vector<size_t> &dim, const iMatrix<T, N> &t, const bool wipe = true)
+  {
+    if (wipe)
+    {
+      dim.clear();
+    }
+    dim.push_back(N);
+    dim.push_back(N);
+    tensorDim(dim, t._internal[0][0], false);
+  }
+
   template <typename T>
   typename TensorToVec<T>::type tensorToVec(const T &t)
   {