diff --git a/TurTLE/TEST.py b/TurTLE/TEST.py
index 0b7d8df5db4f05804d3ed0091f24dff2ec6ecff3..72fc9609ea962787550a6082d3b1e65f0bb0f69e 100644
--- a/TurTLE/TEST.py
+++ b/TurTLE/TEST.py
@@ -125,6 +125,8 @@ class TEST(_code):
             self,
             dns_type = None):
         pars = {}
+        if dns_type in ['field_output_test']:
+            pars['number_operations'] = int(1)
         if dns_type == 'test_interpolation':
             pars['nparticles'] = 3
             pars['tracers0_integration_steps'] = int(4)
diff --git a/cpp/field.cpp b/cpp/field.cpp
index 513c4b9c99d467faca85d399bdcd620728dfcff3..16091a6bdc3ae2dfb3627f467fdf2bff0f0206cf 100644
--- a/cpp/field.cpp
+++ b/cpp/field.cpp
@@ -24,16 +24,18 @@
 
 
 
+#include "field.hpp"
+#include "fftw_tools.hpp"
+#include "scope_timer.hpp"
+#include "shared_array.hpp"
+#include "field_binary_IO.hpp"
+
 #include <sys/stat.h>
 #include <cmath>
 #include <cstdlib>
 #include <algorithm>
 #include <cassert>
 #include <random>
-#include "fftw_tools.hpp"
-#include "field.hpp"
-#include "scope_timer.hpp"
-#include "shared_array.hpp"
 
 template <typename rnumber,
           field_backend be,
@@ -178,6 +180,56 @@ void field<rnumber, be, fc>::dft()
     this->real_space_representation = false;
 }
 
+template <typename rnumber,
+          field_backend be,
+          field_components fc>
+int field<rnumber, be, fc>::io_binary(
+        const std::string fname,
+        const int iteration,
+        const bool read)
+{
+    const std::string full_fname = (
+            fname +
+            "_i" +
+            std::to_string(iteration) +
+            ".bin");
+    if (this->real_space_representation)
+    {
+        field_binary_IO<rnumber, REAL, fc> *bin_IO = new field_binary_IO <rnumber, REAL, fc>(
+                this->rlayout->sizes,
+                this->rlayout->subsizes,
+                this->rlayout->starts,
+                this->rlayout->comm);
+        if(read)
+            bin_IO->read(
+                    full_fname,
+                    this->get_rdata());
+        else
+            bin_IO->write(
+                    full_fname,
+                    this->get_rdata());
+        delete bin_IO;
+    }
+    else
+    {
+        field_binary_IO<rnumber, COMPLEX, fc> *bin_IO = new field_binary_IO <rnumber, COMPLEX, fc>(
+                this->clayout->sizes,
+                this->clayout->subsizes,
+                this->clayout->starts,
+                this->clayout->comm);
+        if(read)
+            bin_IO->read(
+                    full_fname,
+                    this->get_cdata());
+        else
+            bin_IO->write(
+                    full_fname,
+                    this->get_cdata());
+        delete bin_IO;
+    }
+    return EXIT_SUCCESS;
+}
+
 template <typename rnumber,
           field_backend be,
           field_components fc>
diff --git a/cpp/field_binary_IO.cpp b/cpp/field_binary_IO.cpp
index 288ef26a9cdb490895e908c5c3f377ed13c756de..0d11599d5180c73ee590d156f1d82e0200ae0f48 100644
--- a/cpp/field_binary_IO.cpp
+++ b/cpp/field_binary_IO.cpp
@@ -166,6 +166,7 @@ int field_binary_IO<rnumber, fr, fc>::read(
                     mpi_type<rnumber>(fr),
                     MPI_STATUS_IGNORE);
         MPI_File_close(&f);
+        MPI_Info_free(&info);
     }
     return EXIT_SUCCESS;
 }
@@ -205,6 +206,7 @@ int field_binary_IO<rnumber, fr, fc>::write(
                     mpi_type<rnumber>(fr),
                     MPI_STATUS_IGNORE);
         MPI_File_close(&f);
+        MPI_Info_free(&info);
     }
 
     return EXIT_SUCCESS;
diff --git a/cpp/full_code/field_output_test.cpp b/cpp/full_code/field_output_test.cpp
index 89c2883a7cc48d84755004a6ba8b2a486b09dcde..ec90948dab64280ad99cee472a297ba80fd1d97e 100644
--- a/cpp/full_code/field_output_test.cpp
+++ b/cpp/full_code/field_output_test.cpp
@@ -23,12 +23,14 @@
 
 
 
-#include <string>
-#include <cmath>
-#include <random>
 #include "field_output_test.hpp"
 #include "scope_timer.hpp"
+#include "hdf5_tools.hpp"
+#include "shared_array.hpp"
 
+#include <string>
+#include <cmath>
+#include <random>
 
 template <typename rnumber>
 int field_output_test<rnumber>::initialize(void)
@@ -50,6 +52,13 @@ int field_output_test<rnumber>::read_parameters()
 {
     TIMEZONE("field_output_test::read_parameters");
     this->test::read_parameters();
+    hid_t parameter_file = H5Fopen(
+            (this->simname + std::string(".h5")).c_str(),
+            H5F_ACC_RDONLY,
+            H5P_DEFAULT);
+    this->random_seed = hdf5_tools::read_value<int>(parameter_file, "/parameters/field_random_seed");
+    this->number_operations = hdf5_tools::read_value<int>(parameter_file, "/parameters/number_operations");
+    H5Fclose(parameter_file);
     MPI_Barrier(this->comm);
     return EXIT_SUCCESS;
 }
@@ -63,27 +72,105 @@ int field_output_test<rnumber>::do_work(void)
             this->nx, this->ny, this->nz,
             this->comm,
             FFTW_ESTIMATE);
-    std::default_random_engine rgen;
-    std::normal_distribution<rnumber> rdist;
-    rgen.seed(1);
-
-    // fill up scal_field
-    scal_field->real_space_representation = true;
-    scal_field->RLOOP(
-            [&](ptrdiff_t rindex,
-                ptrdiff_t xindex,
-                ptrdiff_t yindex,
-                ptrdiff_t zindex){
-            scal_field->rval(rindex) = rdist(rgen);
-            });
-
-    scal_field->io(
+    field<rnumber, FFTW, ONE> *backup_field = new field<rnumber, FFTW, ONE>(
+            this->nx, this->ny, this->nz,
+            this->comm,
+            FFTW_ESTIMATE);
+
+
+    kspace<FFTW, SMOOTH> *kk = new kspace<FFTW, SMOOTH>(
+            scal_field->clayout, this->dkx, this->dky, this->dkz);
+
+    make_gaussian_random_field(
+        kk,
+        scal_field,
+        random_seed,
+        0.4,
+        1.0,
+        6.0 / this->nx);
+    scal_field->ift();
+
+    *backup_field = *scal_field;
+
+    for (int i=0; i<this->number_operations; i++)
+    {
+        DEBUG_MSG("IO operation %d\n", i);
+        // write data
+        scal_field->io(
             this->simname + std::string("_fields.h5"),
-            "scal_field",
-            0,
+            "scal_field_",
+            i,
             false);
+        scal_field->io_binary(
+            this->simname + std::string("_field.bin"),
+            i,
+            false);
+
+        // read HDF5
+        scal_field->io(
+            this->simname + std::string("_fields.h5"),
+            "scal_field_",
+            i,
+            true);
+        // RLOOP uses a "#pragma openmp parallel"
+        // so we need a shared array
+        shared_array<bool> comparison_0_ok(
+                1,
+                [&](bool* local_comparison){
+                    std::fill_n(local_comparison, 1, true);
+                });
+        scal_field->RLOOP(
+                [&](const ptrdiff_t rindex,
+                    const ptrdiff_t xindex,
+                    const ptrdiff_t yindex,
+                    const ptrdiff_t zindex){
+                comparison_0_ok.getMine()[0] = (
+                        comparison_0_ok.getMine()[0] &&
+                        (scal_field->rval(rindex) ==
+                         backup_field->rval(rindex)));
+                });
+        comparison_0_ok.mergeParallel(
+                [&](const int idx,
+                    const bool& v1,
+                    const bool& v2) -> bool {
+                return v1 && v2;});
+        if (!comparison_0_ok.getMasterData()[0])
+            DEBUG_MSG("error with HDF5 comparison\n");
+
+        // read binary
+        scal_field->io_binary(
+            this->simname + std::string("_field"),
+            i,
+            true);
+        // RLOOP uses a "#pragma openmp parallel"
+        // so we need a shared array
+        shared_array<bool> comparison_1_ok(
+                1,
+                [&](bool* local_comparison){
+                    std::fill_n(local_comparison, 1, true);
+                });
+        scal_field->RLOOP(
+                [&](const ptrdiff_t rindex,
+                    const ptrdiff_t xindex,
+                    const ptrdiff_t yindex,
+                    const ptrdiff_t zindex){
+                comparison_1_ok.getMine()[0] = (
+                        comparison_1_ok.getMine()[0] &&
+                        (scal_field->rval(rindex) ==
+                         backup_field->rval(rindex)));
+                });
+        comparison_1_ok.mergeParallel(
+                [&](const int idx,
+                    const bool& v1,
+                    const bool& v2) -> bool {
+                return v1 && v2;});
+        if (!comparison_1_ok.getMasterData()[0])
+            DEBUG_MSG("error with binary comparison\n");
+    }
 
     // deallocate
+    delete kk;
+    delete backup_field;
     delete scal_field;
     return EXIT_SUCCESS;
 }
diff --git a/cpp/full_code/field_output_test.hpp b/cpp/full_code/field_output_test.hpp
index 72c832450de7d921700cf11df36670f857ba365d..75ecafa391395c228d634844cd433ece99a27ed8 100644
--- a/cpp/full_code/field_output_test.hpp
+++ b/cpp/full_code/field_output_test.hpp
@@ -42,6 +42,10 @@ template <typename rnumber>
 class field_output_test: public test
 {
     public:
+
+        int random_seed;
+        int number_operations;
+
         field_output_test(
                 const MPI_Comm COMMUNICATOR,
                 const std::string &simulation_name):