diff --git a/TurTLE/DNS.py b/TurTLE/DNS.py
index 606014b09a87b490f1dda781fc6bc6e501d12eb8..94a2091f8a7113fc8a13946e3f63795f4e1679db 100644
--- a/TurTLE/DNS.py
+++ b/TurTLE/DNS.py
@@ -997,22 +997,7 @@ class DNS(_code):
                         f,
                         checkpoint_field + '/complex/{0}'.format(0))
             else:
-                if self.dns_type == 'NSVE_Stokes_particles':
-                    data = self.generate_vector_field(
-                       write_to_file = False,
-                       spectra_slope = 2.0,
-                       rseed = self.parameters['field_random_seed'],
-                       amplitude = self.parameters['initial_field_amplitude'])
-                    f[checkpoint_field + '/complex/{0}'.format(0)] = data
-                elif self.dns_type != 'NSVE' or self.parameters['field_random_seed'] == 0:
-                    data = self.generate_vector_field(
-                       write_to_file = False,
-                       spectra_slope = 2.0,
-                       rseed = self.parameters['field_random_seed'],
-                       amplitude = 0.05)
-                    f[checkpoint_field + '/complex/{0}'.format(0)] = data
-                else:
-                    f.require_group(checkpoint_field + '/complex')
+                f.require_group(checkpoint_field + '/complex')
             f.close()
         if self.dns_type == 'kraichnan_field':
             if not os.path.exists(self.get_checkpoint_0_fname()):
@@ -1028,17 +1013,21 @@ class DNS(_code):
             opt = None,
             iter0 = 0):
         """Launches jobs. Prepares initial condition and writes parameters to file if necessary."""
+        if (iter0 == 0):
+            self.check_current_vorticity_exists = False
+            self.check_current_velocity_exists = False
+        else:
+            if (self.dns_type in ['NSVE']):
+                self.check_current_vorticity_exists = True
+                self.check_current_velocity_exists = False
+            if self.dns_type in ['NSE', 'NSE_alt_dealias']:
+                self.check_current_vorticity_exists = False
+                self.check_current_velocity_exists = True
         if not os.path.exists(self.get_data_file_name()):
             if not os.path.isdir(self.work_dir):
                 os.makedirs(self.work_dir)
             self.generate_initial_condition(opt = opt)
             self.write_par(iter0 = iter0)
-        if (('test' in self.dns_type) or
-            (self.dns_type in ['kraichnan_field', 'NSVE'])):
-            self.check_current_vorticity_exists = False
-        if self.dns_type in ['NSE', 'NSE_alt_dealias']:
-            self.check_current_vorticity_exists = False
-            self.check_current_velocity_exists = True
         self.run(
                 nb_processes = opt.nb_processes,
                 nb_threads_per_process = opt.nb_threads_per_process,
diff --git a/TurTLE/_code.py b/TurTLE/_code.py
index 8e9a724ea471d1a8ca1fd72a75633a1c0c8560e9..7c9f4ce169dfb0bf827e50fd5ca8d007a89bdddc 100644
--- a/TurTLE/_code.py
+++ b/TurTLE/_code.py
@@ -179,24 +179,24 @@ class _code(_base):
         if os.path.exists(os.path.join(self.work_dir, 'stop_' + self.simname)):
             warnings.warn("Found stop_<simname> file, I will remove it.")
             os.remove(os.path.join(self.work_dir, 'stop_' + self.simname))
-        if self.check_current_vorticity_exists:
-            with h5py.File(os.path.join(self.work_dir, self.simname + '.h5'), 'r') as data_file:
-                iter0 = data_file['iteration'][...]
-                checkpoint0 = data_file['checkpoint'][()]
-                cp_fname = os.path.join(self.work_dir, self.simname + '_checkpoint_{0}.h5'.format(checkpoint0))
+        with h5py.File(os.path.join(self.work_dir, self.simname + '.h5'), 'r') as data_file:
+            iter0 = data_file['iteration'][...]
+            checkpoint0 = data_file['checkpoint'][()]
+            cp_fname = os.path.join(
+                    self.work_dir,
+                    self.simname + '_checkpoint_{0}.h5'.format(checkpoint0))
+            if (self.check_current_vorticity_exists
+                    or self.check_current_velocity_exists):
+                if self.check_current_vorticity_exists:
+                    checkpoint_field = 'vorticity'
+                elif self.check_current_velocity_exists:
+                    checkpoint_field = 'velocity'
+                else:
+                    checkpoint_field = 'none'
                 assert os.path.exists(cp_fname)
                 with h5py.File(cp_fname, 'r') as checkpoint_file:
-                    assert '{0}'.format(iter0) in checkpoint_file['vorticity/complex'].keys()
-        elif self.check_current_velocity_exists:
-            with h5py.File(os.path.join(self.work_dir, self.simname + '.h5'), 'r') as data_file:
-                iter0 = data_file['iteration'][...]
-                checkpoint0 = data_file['checkpoint'][()]
-                cp_fname = os.path.join(self.work_dir, self.simname + '_checkpoint_{0}.h5'.format(checkpoint0))
-                assert os.path.exists(cp_fname)
-                with h5py.File(cp_fname, 'r') as checkpoint_file:
-                    assert '{0}'.format(iter0) in checkpoint_file['velocity/complex'].keys()
-        else:
-            iter0 = 0
+                    assert '{0}'.format(iter0) in checkpoint_file[
+                            checkpoint_field + '/complex'].keys()
         if not os.path.isdir(self.work_dir):
             os.makedirs(self.work_dir)
         assert (self.compile_code(no_debug = no_debug) == 0)
diff --git a/TurTLE/test/test_collisions.py b/TurTLE/test/test_collisions.py
index 31d97ec432ca1562911805b81cb2d8e13f1a3ad0..31eacbb847292aacf033069256b53ee458a5c65c 100644
--- a/TurTLE/test/test_collisions.py
+++ b/TurTLE/test/test_collisions.py
@@ -26,6 +26,8 @@ class ADNS_Stokes_test(TurTLE.DNS):
             'stokes_number_evolving' : float(0.1),
             'relative_stokes_distr_width' : float(0),
             'nb_growing_particles' : int(0)}
+        self.check_current_vorticity_exists = True
+        self.check_current_velocity_exists = False
         return None
     def custom_parameter_update(self):
         self.extra_parameters['NSVE_Stokes_growing_subset']['tracers0_integration_steps'] = int(0)
@@ -94,7 +96,6 @@ class ADNS_Stokes_test(TurTLE.DNS):
             outfile.write(self.includes + '\n\n')
             outfile.write(self.definitions + '\n\n')
             outfile.write(self.main + '\n')
-        self.check_current_vorticity_exists = True
         return None
     def generate_tracer_state(
             self,
@@ -209,6 +210,16 @@ class ADNS_Stokes_test(TurTLE.DNS):
                 print('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!WARNING')
         return None
 
+    def generate_initial_condition(
+            self,
+            opt = None):
+        f = h5py.File(self.get_checkpoint_0_fname(), 'a')
+        f['vorticity/complex/0'] = self.generate_vector_field(
+                write_to_file = False,
+                rseed = self.parameters['field_random_seed'])
+        f.close()
+        return None
+
     def launch(
             self,
             args = [],
diff --git a/TurTLE/test/test_parameter_copy.py b/TurTLE/test/test_parameter_copy.py
index c5231b3bb53412d35e6e9071582c6b7507352049..b996780bf316b32551c8e25e0b9de4987dcbd682 100644
--- a/TurTLE/test/test_parameter_copy.py
+++ b/TurTLE/test/test_parameter_copy.py
@@ -35,6 +35,23 @@ def main():
             '--field_random_seed', '0',
             ] + random_parameters)
 
+        # generate fields for child simulation
+        cp_fname = os.path.join(
+                c.work_dir,
+                c.simname + '_checkpoint_0.h5')
+        with h5py.File(cp_fname, 'a') as ofile:
+            data = c.generate_vector_field(
+                    write_to_file = False,
+                    spectra_slope = 2.0,
+                    rseed = c.parameters['field_random_seed'],
+                    amplitude = 0.05)
+            if simulation_type == 'NSVE':
+                ofile['vorticity/complex/0'] = data
+            elif simulation_type == 'NSE':
+                ofile['velocity/complex/0'] = data
+            else:
+                assert 0
+
         # generate child simulation
         c = DNS()
         c.launch([
diff --git a/cpp/field.cpp b/cpp/field.cpp
index 5f384f3cb673203e79b0721739f1710357e132e2..f4242295ffbfeb24485cbb2eaee4f02c498137aa 100644
--- a/cpp/field.cpp
+++ b/cpp/field.cpp
@@ -28,9 +28,7 @@
 #include "shared_array.hpp"
 #include "field_binary_IO.hpp"
 
-#include <sys/stat.h>
-#include <algorithm>
-#include <random>
+#include "field_io.cpp"
 
 template <typename rnumber,
           field_backend be,
@@ -151,42 +149,6 @@ field<rnumber, be, fc>::~field() noexcept(false)
     }
 }
 
-template <typename rnumber,
-          field_backend be,
-          field_components fc>
-int field<rnumber, be, fc>::print_plan(const std::string preamble)
-{
-    char *c2r_plan_information = fftw_interface<rnumber>::sprint(this->c2r_plan);
-    char *r2c_plan_information = fftw_interface<rnumber>::sprint(this->r2c_plan);
-    if (this->myrank == 0)
-    {
-        std::cout << preamble <<
-                     std::endl <<
-                     "----c2r plan is:\n" <<
-                     c2r_plan_information <<
-                     std::endl <<
-                     "----r2c plan is:\n" <<
-                     r2c_plan_information <<
-                     std::endl;
-    }
-#ifdef TURTLE_DEBUG_MESSAGES_ON
-    std::string err_message = (
-            std::string("MPI rank ") +
-            preamble +
-            std::to_string(this->myrank) +
-            std::string("\n----c2r plan is:\n") +
-            std::string(c2r_plan_information) +
-            std::string("\n----r2c plan is:\n") +
-            std::string(r2c_plan_information) +
-            std::string("\n"));
-    std::cerr << err_message;
-#endif//TURTLE_DEBUG_MESSAGES_ON
-
-    free(c2r_plan_information);
-    free(r2c_plan_information);
-    return EXIT_SUCCESS;
-}
-
 template <typename rnumber,
           field_backend be,
           field_components fc>
@@ -211,764 +173,6 @@ 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>
-int field<rnumber, be, fc>::io(
-        const std::string fname,
-        const std::string field_name,
-        const int iteration,
-        const bool read)
-{
-    /* file dataset has same dimensions as field */
-    TIMEZONE("field::io");
-    hid_t file_id, dset_id, plist_id;
-    file_id = H5I_BADID;
-    dset_id = H5I_BADID;
-    plist_id = H5I_BADID;
-    std::string representation = std::string(
-            this->real_space_representation ?
-                "real" : "complex");
-    std::string dset_name = (
-            "/" + field_name +
-            "/" + representation +
-            "/" + std::to_string(iteration));
-
-    /* open/create file */
-    start_mpi_profiling_zone(turtle_mpi_pcontrol::HDF5);
-    plist_id = H5Pcreate(H5P_FILE_ACCESS);
-    H5Pset_fapl_mpio(plist_id, this->comm, MPI_INFO_NULL);
-    bool file_exists = false;
-    {
-        struct stat file_buffer;
-        file_exists = (stat(fname.c_str(), &file_buffer) == 0);
-    }
-    if (read)
-    {
-        DEBUG_MSG("field::io trying to read field %s from file %s\n", dset_name.c_str(), fname.c_str());
-        assert(file_exists);
-        file_id = H5Fopen(fname.c_str(), H5F_ACC_RDONLY, plist_id);
-    }
-    else
-    {
-        if (file_exists)
-            file_id = H5Fopen(fname.c_str(), H5F_ACC_RDWR, plist_id);
-        else
-            file_id = H5Fcreate(fname.c_str(), H5F_ACC_EXCL, H5P_DEFAULT, plist_id);
-    }
-    assert(file_id >= 0);
-    H5Pclose(plist_id);
-
-    /* check what kind of representation is being used */
-    if (read)
-    {
-        dset_id = H5Dopen(
-                file_id,
-                dset_name.c_str(),
-                H5P_DEFAULT);
-        assert(dset_id >= 0);
-        hid_t dset_type = H5Dget_type(dset_id);
-        assert(dset_type >= 0);
-        bool io_for_real = (
-                H5Tequal(dset_type, H5T_IEEE_F32BE) ||
-                H5Tequal(dset_type, H5T_IEEE_F32LE) ||
-                H5Tequal(dset_type, H5T_INTEL_F32) ||
-                H5Tequal(dset_type, H5T_NATIVE_FLOAT) ||
-                H5Tequal(dset_type, H5T_IEEE_F64BE) ||
-                H5Tequal(dset_type, H5T_IEEE_F64LE) ||
-                H5Tequal(dset_type, H5T_INTEL_F64) ||
-                H5Tequal(dset_type, H5T_NATIVE_DOUBLE));
-        variable_used_only_in_assert(io_for_real);
-        H5Tclose(dset_type);
-        assert(this->real_space_representation == io_for_real);
-    }
-    finish_mpi_profiling_zone(turtle_mpi_pcontrol::HDF5);
-
-    /* generic space initialization */
-    hid_t fspace, mspace;
-    hsize_t count[ndim(fc)], offset[ndim(fc)], dims[ndim(fc)];
-    hsize_t memoffset[ndim(fc)], memshape[ndim(fc)];
-
-    if (this->real_space_representation)
-    {
-        for (unsigned int i=0; i<ndim(fc); i++)
-        {
-            count[i] = this->rlayout->subsizes[i];
-            offset[i] = this->rlayout->starts[i];
-            dims[i] = this->rlayout->sizes[i];
-            memshape[i] = this->rmemlayout->subsizes[i];
-            memoffset[i] = 0;
-        }
-    }
-    else
-    {
-        for (unsigned int i=0; i<ndim(fc); i++)
-        {
-            count [i] = this->clayout->subsizes[i];
-            offset[i] = this->clayout->starts[i];
-            dims  [i] = this->clayout->sizes[i];
-            memshape [i] = count[i];
-            memoffset[i] = 0;
-        }
-    }
-    start_mpi_profiling_zone(turtle_mpi_pcontrol::HDF5);
-    mspace = H5Screate_simple(ndim(fc), memshape, NULL);
-    H5Sselect_hyperslab(mspace, H5S_SELECT_SET, memoffset, NULL, count, NULL);
-
-    /* open/create data set */
-    if (read)
-        fspace = H5Dget_space(dset_id);
-    else
-    {
-        if (!H5Lexists(file_id, field_name.c_str(), H5P_DEFAULT))
-        {
-            hid_t gid_tmp = H5Gcreate(
-                    file_id, field_name.c_str(),
-                    H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-            H5Gclose(gid_tmp);
-        }
-
-        if (!H5Lexists(file_id, (field_name + "/" + representation).c_str(), H5P_DEFAULT))
-        {
-            hid_t gid_tmp = H5Gcreate(
-                    file_id, ("/" + field_name + "/" + representation).c_str(),
-                    H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-            H5Gclose(gid_tmp);
-        }
-        if (H5Lexists(file_id, dset_name.c_str(), H5P_DEFAULT))
-        {
-            dset_id = H5Dopen(file_id, dset_name.c_str(), H5P_DEFAULT);
-            fspace = H5Dget_space(dset_id);
-        }
-        else
-        {
-            fspace = H5Screate_simple(
-                    ndim(fc),
-                    dims,
-                    NULL);
-            /* chunking needs to go in here */
-            dset_id = H5Dcreate(
-                    file_id,
-                    dset_name.c_str(),
-                    (this->real_space_representation ? this->rnumber_H5T : this->cnumber_H5T),
-                    fspace,
-                    H5P_DEFAULT,
-                    H5P_DEFAULT,
-                    H5P_DEFAULT);
-            assert(dset_id > 0);
-        }
-    }
-    finish_mpi_profiling_zone(turtle_mpi_pcontrol::HDF5);
-    /* both dset_id and fspace should now have sane values */
-
-    /* if reading, first set local data to 0 */
-    if (read)
-        *this = rnumber(0.0);
-
-    start_mpi_profiling_zone(turtle_mpi_pcontrol::HDF5);
-    /* both dset_id and fspace should now have sane values */
-    /* check file space */
-    int ndims_fspace = H5Sget_simple_extent_dims(fspace, dims, NULL);
-    variable_used_only_in_assert(ndims_fspace);
-    assert(((unsigned int)(ndims_fspace)) == ndim(fc));
-    if (this->real_space_representation)
-    {
-        for (unsigned int i=0; i<ndim(fc); i++)
-        {
-            offset[i] = this->rlayout->starts[i];
-            assert(dims[i] == this->rlayout->sizes[i]);
-        }
-        H5Sselect_hyperslab(fspace, H5S_SELECT_SET, offset, NULL, count, NULL);
-        if (read)
-        {
-            H5Dread(dset_id, this->rnumber_H5T, mspace, fspace, H5P_DEFAULT, this->data);
-        }
-        else
-        {
-            assert(this->real_space_representation);
-            H5Dwrite(dset_id, this->rnumber_H5T, mspace, fspace, H5P_DEFAULT, this->data);
-        }
-        H5Sclose(mspace);
-    }
-    else
-    {
-        for (unsigned int i=0; i<ndim(fc); i++)
-        {
-            offset[i] = this->clayout->starts[i];
-            assert(dims[i] == this->clayout->sizes[i]);
-        }
-        H5Sselect_hyperslab(fspace, H5S_SELECT_SET, offset, NULL, count, NULL);
-        if (read)
-        {
-            H5Dread(dset_id, this->cnumber_H5T, mspace, fspace, H5P_DEFAULT, this->data);
-            this->symmetrize();
-        }
-        else
-        {
-            assert(!this->real_space_representation);
-            H5Dwrite(dset_id, this->cnumber_H5T, mspace, fspace, H5P_DEFAULT, this->data);
-        }
-        H5Sclose(mspace);
-    }
-
-    H5Sclose(fspace);
-    /* close data set */
-    H5Dclose(dset_id);
-    /* close file */
-    H5Fclose(file_id);
-    /* ensure all processes are finished writing before exiting the method */
-    MPI_Barrier(this->comm);
-    finish_mpi_profiling_zone(turtle_mpi_pcontrol::HDF5);
-    return EXIT_SUCCESS;
-}
-
-template <typename rnumber,
-          field_backend be,
-          field_components fc>
-int field<rnumber, be, fc>::io_database(
-        const std::string fname,
-        const std::string field_name,
-        const int toffset,
-        const bool read)
-{
-    /* file dataset is has a time dimension as well */
-    TIMEZONE("field::io_database");
-    hid_t file_id, dset_id, plist_id;
-    dset_id = H5I_BADID;
-    std::string representation = std::string(
-            this->real_space_representation ?
-                "real" : "complex");
-    std::string dset_name = (
-            "/" + field_name +
-            "/" + representation);
-
-    /* open/create file */
-    plist_id = H5Pcreate(H5P_FILE_ACCESS);
-    H5Pset_fapl_mpio(plist_id, this->comm, MPI_INFO_NULL);
-    bool file_exists = false;
-    {
-        struct stat file_buffer;
-        file_exists = (stat(fname.c_str(), &file_buffer) == 0);
-    }
-    if (read)
-    {
-        DEBUG_MSG("field::io trying to read field from file %s\n", fname.c_str());
-        assert(file_exists);
-        file_id = H5Fopen(fname.c_str(), H5F_ACC_RDONLY, plist_id);
-    }
-    else
-    {
-        if (file_exists)
-            file_id = H5Fopen(fname.c_str(), H5F_ACC_RDWR, plist_id);
-        else
-            file_id = H5Fcreate(fname.c_str(), H5F_ACC_EXCL, H5P_DEFAULT, plist_id);
-    }
-    H5Pclose(plist_id);
-
-    /* check what kind of representation is being used */
-    if (read)
-    {
-        dset_id = H5Dopen(
-                file_id,
-                dset_name.c_str(),
-                H5P_DEFAULT);
-        hid_t dset_type = H5Dget_type(dset_id);
-        bool io_for_real = (
-                H5Tequal(dset_type, H5T_IEEE_F32BE) ||
-                H5Tequal(dset_type, H5T_IEEE_F32LE) ||
-                H5Tequal(dset_type, H5T_INTEL_F32) ||
-                H5Tequal(dset_type, H5T_NATIVE_FLOAT) ||
-                H5Tequal(dset_type, H5T_IEEE_F64BE) ||
-                H5Tequal(dset_type, H5T_IEEE_F64LE) ||
-                H5Tequal(dset_type, H5T_INTEL_F64) ||
-                H5Tequal(dset_type, H5T_NATIVE_DOUBLE));
-        variable_used_only_in_assert(io_for_real);
-        H5Tclose(dset_type);
-        assert(this->real_space_representation == io_for_real);
-    }
-
-    /* generic space initialization */
-    hid_t fspace, mspace;
-    hsize_t count[ndim(fc)+1], offset[ndim(fc)+1], dims[ndim(fc)+1];
-    hsize_t memoffset[ndim(fc)+1], memshape[ndim(fc)+1];
-
-    int dim_counter_offset = 1;
-    dim_counter_offset = 1;
-    count[0] = 1;
-    memshape[0] = 1;
-    memoffset[0] = 0;
-    if (this->real_space_representation)
-    {
-        for (unsigned int i=0; i<ndim(fc); i++)
-        {
-            count[i+dim_counter_offset] = this->rlayout->subsizes[i];
-            offset[i+dim_counter_offset] = this->rlayout->starts[i];
-            dims[i+dim_counter_offset] = this->rlayout->sizes[i];
-            memshape[i+dim_counter_offset] = this->rmemlayout->subsizes[i];
-            memoffset[i+dim_counter_offset] = 0;
-        }
-        mspace = H5Screate_simple(dim_counter_offset + ndim(fc), memshape, NULL);
-        H5Sselect_hyperslab(mspace, H5S_SELECT_SET, memoffset, NULL, count, NULL);
-    }
-    else
-    {
-        for (unsigned int i=0; i<ndim(fc); i++)
-        {
-            count[i+dim_counter_offset] = this->clayout->subsizes[i];
-            offset[i+dim_counter_offset] = this->clayout->starts[i];
-            dims[i+dim_counter_offset] = this->clayout->sizes[i];
-            memshape[i+dim_counter_offset] = count[i+dim_counter_offset];
-            memoffset[i+dim_counter_offset] = 0;
-        }
-        mspace = H5Screate_simple(dim_counter_offset + ndim(fc), memshape, NULL);
-        H5Sselect_hyperslab(mspace, H5S_SELECT_SET, memoffset, NULL, count, NULL);
-    }
-
-    /* open/create data set */
-    if (read)
-        fspace = H5Dget_space(dset_id);
-    else
-    {
-        if (!H5Lexists(file_id, field_name.c_str(), H5P_DEFAULT))
-            H5Gcreate(
-                    file_id, field_name.c_str(),
-                    H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-        if (H5Lexists(file_id, dset_name.c_str(), H5P_DEFAULT))
-        {
-            dset_id = H5Dopen(file_id, dset_name.c_str(), H5P_DEFAULT);
-            fspace = H5Dget_space(dset_id);
-        }
-        else
-        {
-            fspace = H5Screate_simple(
-                    ndim(fc),
-                    dims,
-                    NULL);
-            /* chunking needs to go in here */
-            dset_id = H5Dcreate(
-                    file_id,
-                    dset_name.c_str(),
-                    (this->real_space_representation ? this->rnumber_H5T : this->cnumber_H5T),
-                    fspace,
-                    H5P_DEFAULT,
-                    H5P_DEFAULT,
-                    H5P_DEFAULT);
-            assert(dset_id > 0);
-        }
-    }
-    /* both dset_id and fspace should now have sane values */
-
-    if (read)
-        *this = rnumber(0.0);
-
-    /* check file space */
-    int ndims_fspace = H5Sget_simple_extent_dims(fspace, dims, NULL);
-    variable_used_only_in_assert(ndims_fspace);
-    assert(ndims_fspace == int(ndim(fc) + 1));
-    offset[0] = toffset;
-    if (this->real_space_representation)
-    {
-        for (unsigned int i=0; i<ndim(fc); i++)
-        {
-            offset[i+dim_counter_offset] = this->rlayout->starts[i];
-            assert(dims[i+dim_counter_offset] == this->rlayout->sizes[i]);
-        }
-        H5Sselect_hyperslab(fspace, H5S_SELECT_SET, offset, NULL, count, NULL);
-        if (read)
-        {
-            H5Dread(dset_id, this->rnumber_H5T, mspace, fspace, H5P_DEFAULT, this->data);
-            this->real_space_representation = true;
-        }
-        else
-        {
-            assert(this->real_space_representation);
-            H5Dwrite(dset_id, this->rnumber_H5T, mspace, fspace, H5P_DEFAULT, this->data);
-        }
-        H5Sclose(mspace);
-    }
-    else
-    {
-        for (unsigned int i=0; i<ndim(fc); i++)
-        {
-            offset[i+dim_counter_offset] = this->clayout->starts[i];
-            assert(dims[i+dim_counter_offset] == this->clayout->sizes[i]);
-        }
-        H5Sselect_hyperslab(fspace, H5S_SELECT_SET, offset, NULL, count, NULL);
-        if (read)
-        {
-            H5Dread(dset_id, this->cnumber_H5T, mspace, fspace, H5P_DEFAULT, this->data);
-            this->real_space_representation = false;
-            this->symmetrize();
-        }
-        else
-        {
-            assert(!this->real_space_representation);
-            H5Dwrite(dset_id, this->cnumber_H5T, mspace, fspace, H5P_DEFAULT, this->data);
-        }
-        H5Sclose(mspace);
-    }
-
-    H5Sclose(fspace);
-    /* close data set */
-    H5Dclose(dset_id);
-    /* close file */
-    H5Fclose(file_id);
-    return EXIT_SUCCESS;
-}
-
-
-template <typename rnumber,
-          field_backend be,
-          field_components fc>
-int field<rnumber, be, fc>::write_0slice(
-        const hid_t group,
-        const std::string field_name,
-        const int iteration)
-{
-    start_mpi_profiling_zone(turtle_mpi_pcontrol::HDF5);
-    // this should in principle work for any fc
-    TIMEZONE("field::write_0slice");
-    assert(this->real_space_representation);
-    if (this->myrank == 0)
-    {
-        hid_t dset, wspace, mspace;
-        int ndims;
-        hsize_t count[5], offset[5], dims[5];
-        offset[0] = iteration;
-        offset[1] = 0;
-        offset[2] = 0;
-        offset[3] = 0;
-        offset[4] = 0;
-        dset = H5Dopen(
-                group,
-                ("0slices/" + field_name + "/real").c_str(),
-                H5P_DEFAULT);
-        wspace = H5Dget_space(dset);
-        ndims = H5Sget_simple_extent_dims(wspace, dims, NULL);
-        // array in memory has 2 extra x points, because FFTW
-        count[0] = 1;
-        count[1] = this->rmemlayout->sizes[1];
-        count[2] = this->rmemlayout->sizes[2];
-        count[3] = 3;
-        count[4] = 3;
-        mspace = H5Screate_simple(ndims, count, NULL);
-        // array in file should not have the extra 2 points
-        count[1] = this->rlayout->sizes[1];
-        count[2] = this->rlayout->sizes[2];
-        // select right slice in file
-        H5Sselect_hyperslab(
-            wspace,
-            H5S_SELECT_SET,
-            offset,
-            NULL,
-            count,
-            NULL);
-        offset[0] = 0;
-        // select proper regions of memory
-        H5Sselect_hyperslab(
-            mspace,
-            H5S_SELECT_SET,
-            offset,
-            NULL,
-            count,
-            NULL);
-        H5Dwrite(
-            dset,
-            this->rnumber_H5T,
-            mspace,
-            wspace,
-            H5P_DEFAULT,
-            this->data);
-        H5Dclose(dset);
-        H5Sclose(mspace);
-        H5Sclose(wspace);
-    }
-    finish_mpi_profiling_zone(turtle_mpi_pcontrol::HDF5);
-    return EXIT_SUCCESS;
-}
-
-template <typename rnumber,
-          field_backend be,
-          field_components fc>
-int field<rnumber, be, fc>::write_filtered(
-        const std::string fname,
-        const std::string field_name,
-        const int iteration,
-        int nx,
-        int ny,
-        int nz)
-{
-    /* file dataset has same dimensions as field */
-    TIMEZONE("field::write_filtered");
-    // only works in Fourier representation
-    assert(!this->real_space_representation);
-    assert(hsize_t(nx) <= this->rlayout->sizes[2]);
-    assert(hsize_t(ny) <= this->rlayout->sizes[1]);
-    assert(hsize_t(nz) <= this->rlayout->sizes[0]);
-    // current algorithm only works for more than one process
-    assert(this->nprocs >= 2);
-    hid_t file_id, dset_id, plist_id;
-    dset_id = H5I_BADID;
-    std::string dset_name = (
-            "/" + field_name +
-            "/complex" +
-            "/" + std::to_string(iteration));
-
-    /* open/create file */
-    plist_id = H5Pcreate(H5P_FILE_ACCESS);
-    H5Pset_fapl_mpio(plist_id, this->comm, MPI_INFO_NULL);
-    bool file_exists = false;
-    struct stat file_buffer;
-    file_exists = (stat(fname.c_str(), &file_buffer) == 0);
-    if (file_exists)
-        file_id = H5Fopen(fname.c_str(), H5F_ACC_RDWR, plist_id);
-    else
-        file_id = H5Fcreate(fname.c_str(), H5F_ACC_EXCL, H5P_DEFAULT, plist_id);
-    assert(file_id >= 0);
-    H5Pclose(plist_id);
-
-    /* generic space initialization */
-    hid_t fspace, mspace;
-    hsize_t count[ndim(fc)], offset[ndim(fc)], dims[ndim(fc)], fdims[ndim(fc)];
-    hsize_t memoffset[ndim(fc)], memshape[ndim(fc)];
-
-    // set up dimensions
-    for (unsigned int i=3; i<ndim(fc); i++)
-    {
-        count [i] = this->clayout->subsizes[i];
-        offset[i] = this->clayout->starts[i];
-        dims  [i] = this->clayout->sizes[i];
-        memshape [i] = count[i];
-        memoffset[i] = 0;
-    }
-    // these are dimensions of dataset, needed
-    // to create dataset
-    dims[0] = ny;
-    dims[1] = nz;
-    dims[2] = nx/2+1;
-
-    /* open/create data set */
-    if (!H5Lexists(file_id, field_name.c_str(), H5P_DEFAULT))
-    {
-        hid_t gid_tmp = H5Gcreate(
-                file_id, field_name.c_str(),
-                H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-        H5Gclose(gid_tmp);
-    }
-    if (!H5Lexists(file_id, (field_name + "/complex").c_str(), H5P_DEFAULT))
-    {
-        hid_t gid_tmp = H5Gcreate(
-                file_id, ("/" + field_name + "/complex").c_str(),
-                H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-        H5Gclose(gid_tmp);
-    }
-    if (H5Lexists(file_id, dset_name.c_str(), H5P_DEFAULT))
-    {
-        dset_id = H5Dopen(file_id, dset_name.c_str(), H5P_DEFAULT);
-        fspace = H5Dget_space(dset_id);
-    }
-    else
-    {
-        if (nz == 1)
-        {
-            hsize_t temp_dims[ndim(fc)-1];
-            temp_dims[0] = dims[0];
-            for (unsigned int i=1; i<ndim(fc)-1; i++)
-                temp_dims[i] = dims[i+1];
-            fspace = H5Screate_simple(
-                    ndim(fc)-1,
-                    temp_dims,
-                    NULL);
-        }
-        else
-            fspace = H5Screate_simple(
-                    ndim(fc),
-                    dims,
-                    NULL);
-        /* chunking needs to go in here */
-        dset_id = H5Dcreate(
-                file_id,
-                dset_name.c_str(),
-                this->cnumber_H5T,
-                fspace,
-                H5P_DEFAULT,
-                H5P_DEFAULT,
-                H5P_DEFAULT);
-        assert(dset_id > 0);
-    }
-    /* check file space */
-    int ndims_fspace = H5Sget_simple_extent_dims(fspace, fdims, NULL);
-    variable_used_only_in_assert(ndims_fspace);
-    if (nz == 1)
-    {
-        assert(((unsigned int)(ndims_fspace)) == ndim(fc)-1);
-        assert(dims[0] == fdims[0]);
-        for (unsigned int i=1; i<ndim(fc)-1; i++)
-            assert(dims[i+1] == fdims[i]);
-    }
-    else
-    {
-        assert(((unsigned int)(ndims_fspace)) == ndim(fc));
-        for (unsigned int i=0; i<ndim(fc); i++)
-        {
-            assert(dims[i] == fdims[i]);
-        }
-    }
-    /* both dset_id and fspace now have sane values */
-
-    /// set up counts and offsets
-    /// x is easy, since only positive modes are present
-    count [2] = nx/2+1;
-    offset[2] = 0;
-    memshape [2] = this->clayout->subsizes[2];
-    memoffset[2] = 0;
-
-    /// three options for y:
-    /// this->starts[0] <= ny/2
-    /// ny / 2 < this->starts[0] +this->clayout->subsizes[0] < this->sizes[0] - ny/2
-    /// this->starts[0] >= this->sizes[0] - ny/2
-    /// we don't care about saving the ny/2 mode, because of symmetry
-    hsize_t y0 = this->clayout->starts[0];
-    hsize_t y1 = this->clayout->starts[0] + this->clayout->subsizes[0];
-    memshape[0] = this->clayout->subsizes[0];
-    if (y1 <= hsize_t(ny/2))
-    {
-        count[0] = this->clayout->subsizes[0];
-        offset[0] = y0;
-        memoffset[0] = 0;
-    }
-    else
-    {
-        if (y0 < hsize_t(ny)/2)
-        {
-            count[0] = ny/2 - y0;
-            offset[0] = y0;
-            memoffset[0] = 0;
-        }
-        else
-        {
-            if (y1 <= hsize_t(this->clayout->sizes[0] - ny/2 + 1))
-            { // y0 < y1 therefore y0 <= this->clayout->sizes[0] - ny/2
-                count[0] = 0;
-                offset[0] = ny/2;
-                memoffset[0] = 0;
-            }
-            else
-            {
-                if (y0 <= hsize_t(this->clayout->sizes[0] - ny/2))
-                {
-                    count[0] = y1 - (this->clayout->sizes[0] - ny/2);
-                    offset[0] = ny - (this->clayout->sizes[0] - y0);
-                    memoffset[0] = this->clayout->subsizes[0] - count[0];
-                }
-                else
-                {
-                    count[0] = this->clayout->subsizes[0];
-                    offset[0] = ny - (this->clayout->sizes[0] - y0);
-                    memoffset[0] = 0;
-                }
-            }
-        }
-    }
-    if (nz>=2)
-    {
-        assert(nz%2==0);
-        /// for z, we need to take into account that there are
-        /// both positive and negative modes
-        for (int cz = 0; cz < 2; cz++)
-        {
-            count [1] = nz/2;
-            offset[1] = cz*nz/2;
-            memshape [1] = this->clayout->sizes[1];
-            memoffset[1] = cz*(this->clayout->sizes[1] - nz/2);
-
-            //now write data
-            mspace = H5Screate_simple(ndim(fc), memshape, NULL);
-            H5Sselect_hyperslab(mspace, H5S_SELECT_SET, memoffset, NULL, count, NULL);
-            H5Sselect_hyperslab(fspace, H5S_SELECT_SET, offset, NULL, count, NULL);
-            H5Dwrite(dset_id, this->cnumber_H5T, mspace, fspace, H5P_DEFAULT, this->data);
-            H5Sclose(mspace);
-        }
-    }
-    else
-    {
-        assert(nz == 1);
-        count [1] = 1;
-        offset[1] = 0;
-        memshape [1] = this->clayout->sizes[1];
-        memoffset[1] = 0;
-
-        //now write data
-        mspace = H5Screate_simple(ndim(fc), memshape, NULL);
-        H5Sselect_hyperslab(mspace, H5S_SELECT_SET, memoffset, NULL, count, NULL);
-        for (unsigned int i=1; i<ndim(fc)-1; i++)
-            count[i] = count[i+1];
-        H5Sselect_hyperslab(fspace, H5S_SELECT_SET, offset, NULL, count, NULL);
-        H5Dwrite(dset_id, this->cnumber_H5T, mspace, fspace, H5P_DEFAULT, this->data);
-        H5Sclose(mspace);
-    }
-
-
-    /* close file data space */
-    H5Sclose(fspace);
-    /* close data set */
-    H5Dclose(dset_id);
-    /* close file */
-    H5Fclose(file_id);
-    return EXIT_SUCCESS;
-}
-
-
 template <typename rnumber,
           field_backend be,
           field_components fc>
diff --git a/cpp/field_io.cpp b/cpp/field_io.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7d3c55ac4ef32dd43a89e14112f8470c912a3468
--- /dev/null
+++ b/cpp/field_io.cpp
@@ -0,0 +1,834 @@
+/**********************************************************************
+*                                                                     *
+*  Copyright 2024 the TurTLE team                                     *
+*                                                                     *
+*  This file is part of TurTLE.                                       *
+*                                                                     *
+*  TurTLE is free software: you can redistribute it and/or modify     *
+*  it under the terms of the GNU General Public License as published  *
+*  by the Free Software Foundation, either version 3 of the License,  *
+*  or (at your option) any later version.                             *
+*                                                                     *
+*  TurTLE is distributed in the hope that it will be useful,          *
+*  but WITHOUT ANY WARRANTY; without even the implied warranty of     *
+*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the      *
+*  GNU General Public License for more details.                       *
+*                                                                     *
+*  You should have received a copy of the GNU General Public License  *
+*  along with TurTLE.  If not, see <http://www.gnu.org/licenses/>     *
+*                                                                     *
+* Contact: Cristian.Lalescu@mpcdf.mpg.de                              *
+*                                                                     *
+**********************************************************************/
+
+
+#include <sys/stat.h>
+#include <algorithm>
+#include <random>
+
+template <typename rnumber,
+          field_backend be,
+          field_components fc>
+int field<rnumber, be, fc>::print_plan(const std::string preamble)
+{
+    char *c2r_plan_information = fftw_interface<rnumber>::sprint(this->c2r_plan);
+    char *r2c_plan_information = fftw_interface<rnumber>::sprint(this->r2c_plan);
+    if (this->myrank == 0)
+    {
+        std::cout << preamble <<
+                     std::endl <<
+                     "----c2r plan is:\n" <<
+                     c2r_plan_information <<
+                     std::endl <<
+                     "----r2c plan is:\n" <<
+                     r2c_plan_information <<
+                     std::endl;
+    }
+#ifdef TURTLE_DEBUG_MESSAGES_ON
+    std::string err_message = (
+            std::string("MPI rank ") +
+            preamble +
+            std::to_string(this->myrank) +
+            std::string("\n----c2r plan is:\n") +
+            std::string(c2r_plan_information) +
+            std::string("\n----r2c plan is:\n") +
+            std::string(r2c_plan_information) +
+            std::string("\n"));
+    std::cerr << err_message;
+#endif//TURTLE_DEBUG_MESSAGES_ON
+
+    free(c2r_plan_information);
+    free(r2c_plan_information);
+    return EXIT_SUCCESS;
+}
+
+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>
+int field<rnumber, be, fc>::io(
+        const std::string fname,
+        const std::string field_name,
+        const int iteration,
+        const bool read)
+{
+    /* file dataset has same dimensions as field */
+    TIMEZONE("field::io");
+    hid_t file_id, dset_id, plist_id;
+    file_id = H5I_BADID;
+    dset_id = H5I_BADID;
+    plist_id = H5I_BADID;
+    std::string representation = std::string(
+            this->real_space_representation ?
+                "real" : "complex");
+    std::string dset_name = (
+            "/" + field_name +
+            "/" + representation +
+            "/" + std::to_string(iteration));
+
+    /* open/create file */
+    start_mpi_profiling_zone(turtle_mpi_pcontrol::HDF5);
+    plist_id = H5Pcreate(H5P_FILE_ACCESS);
+    H5Pset_fapl_mpio(plist_id, this->comm, MPI_INFO_NULL);
+    bool file_exists = false;
+    {
+        struct stat file_buffer;
+        file_exists = (stat(fname.c_str(), &file_buffer) == 0);
+    }
+    if (read)
+    {
+        DEBUG_MSG("field::io trying to read field %s from file %s\n", dset_name.c_str(), fname.c_str());
+        assert(file_exists);
+        file_id = H5Fopen(fname.c_str(), H5F_ACC_RDONLY, plist_id);
+    }
+    else
+    {
+        if (file_exists)
+            file_id = H5Fopen(fname.c_str(), H5F_ACC_RDWR, plist_id);
+        else
+            file_id = H5Fcreate(fname.c_str(), H5F_ACC_EXCL, H5P_DEFAULT, plist_id);
+    }
+    assert(file_id >= 0);
+    if (file_id < 0) {
+        DEBUG_MSG("couldn't open file");
+        throw std::runtime_error(
+                std::string("Couldn't open file for field I/O.\n")
+                + std::string("file name = ") + fname + "\n");
+    }
+    H5Pclose(plist_id);
+
+    /* check what kind of representation is being used */
+    if (read)
+    {
+        dset_id = H5Dopen(
+                file_id,
+                dset_name.c_str(),
+                H5P_DEFAULT);
+        assert(dset_id >= 0);
+        if (dset_id < 0) {
+            DEBUG_MSG("couldn't open dataset");
+            throw std::runtime_error(
+                    std::string("Couldn't open dataset for field I/O.\n")
+                    + std::string("file name = ") + fname + "\n"
+                    + std::string("dataset name = ") + dset_name + "\n");
+        }
+        hid_t dset_type = H5Dget_type(dset_id);
+        assert(dset_type >= 0);
+        bool io_for_real = (
+                H5Tequal(dset_type, H5T_IEEE_F32BE) ||
+                H5Tequal(dset_type, H5T_IEEE_F32LE) ||
+                H5Tequal(dset_type, H5T_INTEL_F32) ||
+                H5Tequal(dset_type, H5T_NATIVE_FLOAT) ||
+                H5Tequal(dset_type, H5T_IEEE_F64BE) ||
+                H5Tequal(dset_type, H5T_IEEE_F64LE) ||
+                H5Tequal(dset_type, H5T_INTEL_F64) ||
+                H5Tequal(dset_type, H5T_NATIVE_DOUBLE));
+        variable_used_only_in_assert(io_for_real);
+        H5Tclose(dset_type);
+        assert(this->real_space_representation == io_for_real);
+    }
+    finish_mpi_profiling_zone(turtle_mpi_pcontrol::HDF5);
+
+    /* generic space initialization */
+    hid_t fspace, mspace;
+    hsize_t count[ndim(fc)], offset[ndim(fc)], dims[ndim(fc)];
+    hsize_t memoffset[ndim(fc)], memshape[ndim(fc)];
+
+    if (this->real_space_representation)
+    {
+        for (unsigned int i=0; i<ndim(fc); i++)
+        {
+            count[i] = this->rlayout->subsizes[i];
+            offset[i] = this->rlayout->starts[i];
+            dims[i] = this->rlayout->sizes[i];
+            memshape[i] = this->rmemlayout->subsizes[i];
+            memoffset[i] = 0;
+        }
+    }
+    else
+    {
+        for (unsigned int i=0; i<ndim(fc); i++)
+        {
+            count [i] = this->clayout->subsizes[i];
+            offset[i] = this->clayout->starts[i];
+            dims  [i] = this->clayout->sizes[i];
+            memshape [i] = count[i];
+            memoffset[i] = 0;
+        }
+    }
+    start_mpi_profiling_zone(turtle_mpi_pcontrol::HDF5);
+    mspace = H5Screate_simple(ndim(fc), memshape, NULL);
+    H5Sselect_hyperslab(mspace, H5S_SELECT_SET, memoffset, NULL, count, NULL);
+
+    /* open/create data set */
+    if (read)
+        fspace = H5Dget_space(dset_id);
+    else
+    {
+        if (!H5Lexists(file_id, field_name.c_str(), H5P_DEFAULT))
+        {
+            hid_t gid_tmp = H5Gcreate(
+                    file_id, field_name.c_str(),
+                    H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+            H5Gclose(gid_tmp);
+        }
+
+        if (!H5Lexists(file_id, (field_name + "/" + representation).c_str(), H5P_DEFAULT))
+        {
+            hid_t gid_tmp = H5Gcreate(
+                    file_id, ("/" + field_name + "/" + representation).c_str(),
+                    H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+            H5Gclose(gid_tmp);
+        }
+        if (H5Lexists(file_id, dset_name.c_str(), H5P_DEFAULT))
+        {
+            dset_id = H5Dopen(file_id, dset_name.c_str(), H5P_DEFAULT);
+            fspace = H5Dget_space(dset_id);
+        }
+        else
+        {
+            fspace = H5Screate_simple(
+                    ndim(fc),
+                    dims,
+                    NULL);
+            /* chunking needs to go in here */
+            dset_id = H5Dcreate(
+                    file_id,
+                    dset_name.c_str(),
+                    (this->real_space_representation ? this->rnumber_H5T : this->cnumber_H5T),
+                    fspace,
+                    H5P_DEFAULT,
+                    H5P_DEFAULT,
+                    H5P_DEFAULT);
+            assert(dset_id > 0);
+        }
+    }
+    finish_mpi_profiling_zone(turtle_mpi_pcontrol::HDF5);
+    /* both dset_id and fspace should now have sane values */
+
+    /* if reading, first set local data to 0 */
+    if (read)
+        *this = rnumber(0.0);
+
+    start_mpi_profiling_zone(turtle_mpi_pcontrol::HDF5);
+    /* both dset_id and fspace should now have sane values */
+    /* check file space */
+    int ndims_fspace = H5Sget_simple_extent_dims(fspace, dims, NULL);
+    variable_used_only_in_assert(ndims_fspace);
+    assert(((unsigned int)(ndims_fspace)) == ndim(fc));
+    if (this->real_space_representation)
+    {
+        for (unsigned int i=0; i<ndim(fc); i++)
+        {
+            offset[i] = this->rlayout->starts[i];
+            assert(dims[i] == this->rlayout->sizes[i]);
+        }
+        H5Sselect_hyperslab(fspace, H5S_SELECT_SET, offset, NULL, count, NULL);
+        if (read)
+        {
+            H5Dread(dset_id, this->rnumber_H5T, mspace, fspace, H5P_DEFAULT, this->data);
+        }
+        else
+        {
+            assert(this->real_space_representation);
+            H5Dwrite(dset_id, this->rnumber_H5T, mspace, fspace, H5P_DEFAULT, this->data);
+        }
+        H5Sclose(mspace);
+    }
+    else
+    {
+        for (unsigned int i=0; i<ndim(fc); i++)
+        {
+            offset[i] = this->clayout->starts[i];
+            assert(dims[i] == this->clayout->sizes[i]);
+        }
+        H5Sselect_hyperslab(fspace, H5S_SELECT_SET, offset, NULL, count, NULL);
+        if (read)
+        {
+            H5Dread(dset_id, this->cnumber_H5T, mspace, fspace, H5P_DEFAULT, this->data);
+            this->symmetrize();
+        }
+        else
+        {
+            assert(!this->real_space_representation);
+            H5Dwrite(dset_id, this->cnumber_H5T, mspace, fspace, H5P_DEFAULT, this->data);
+        }
+        H5Sclose(mspace);
+    }
+
+    H5Sclose(fspace);
+    /* close data set */
+    H5Dclose(dset_id);
+    /* close file */
+    H5Fclose(file_id);
+    /* ensure all processes are finished writing before exiting the method */
+    MPI_Barrier(this->comm);
+    finish_mpi_profiling_zone(turtle_mpi_pcontrol::HDF5);
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber,
+          field_backend be,
+          field_components fc>
+int field<rnumber, be, fc>::io_database(
+        const std::string fname,
+        const std::string field_name,
+        const int toffset,
+        const bool read)
+{
+    /* file dataset is has a time dimension as well */
+    TIMEZONE("field::io_database");
+    hid_t file_id, dset_id, plist_id;
+    dset_id = H5I_BADID;
+    std::string representation = std::string(
+            this->real_space_representation ?
+                "real" : "complex");
+    std::string dset_name = (
+            "/" + field_name +
+            "/" + representation);
+
+    /* open/create file */
+    plist_id = H5Pcreate(H5P_FILE_ACCESS);
+    H5Pset_fapl_mpio(plist_id, this->comm, MPI_INFO_NULL);
+    bool file_exists = false;
+    {
+        struct stat file_buffer;
+        file_exists = (stat(fname.c_str(), &file_buffer) == 0);
+    }
+    if (read)
+    {
+        DEBUG_MSG("field::io trying to read field from file %s\n", fname.c_str());
+        assert(file_exists);
+        file_id = H5Fopen(fname.c_str(), H5F_ACC_RDONLY, plist_id);
+    }
+    else
+    {
+        if (file_exists)
+            file_id = H5Fopen(fname.c_str(), H5F_ACC_RDWR, plist_id);
+        else
+            file_id = H5Fcreate(fname.c_str(), H5F_ACC_EXCL, H5P_DEFAULT, plist_id);
+    }
+    H5Pclose(plist_id);
+
+    /* check what kind of representation is being used */
+    if (read)
+    {
+        dset_id = H5Dopen(
+                file_id,
+                dset_name.c_str(),
+                H5P_DEFAULT);
+        hid_t dset_type = H5Dget_type(dset_id);
+        bool io_for_real = (
+                H5Tequal(dset_type, H5T_IEEE_F32BE) ||
+                H5Tequal(dset_type, H5T_IEEE_F32LE) ||
+                H5Tequal(dset_type, H5T_INTEL_F32) ||
+                H5Tequal(dset_type, H5T_NATIVE_FLOAT) ||
+                H5Tequal(dset_type, H5T_IEEE_F64BE) ||
+                H5Tequal(dset_type, H5T_IEEE_F64LE) ||
+                H5Tequal(dset_type, H5T_INTEL_F64) ||
+                H5Tequal(dset_type, H5T_NATIVE_DOUBLE));
+        variable_used_only_in_assert(io_for_real);
+        H5Tclose(dset_type);
+        assert(this->real_space_representation == io_for_real);
+    }
+
+    /* generic space initialization */
+    hid_t fspace, mspace;
+    hsize_t count[ndim(fc)+1], offset[ndim(fc)+1], dims[ndim(fc)+1];
+    hsize_t memoffset[ndim(fc)+1], memshape[ndim(fc)+1];
+
+    int dim_counter_offset = 1;
+    dim_counter_offset = 1;
+    count[0] = 1;
+    memshape[0] = 1;
+    memoffset[0] = 0;
+    if (this->real_space_representation)
+    {
+        for (unsigned int i=0; i<ndim(fc); i++)
+        {
+            count[i+dim_counter_offset] = this->rlayout->subsizes[i];
+            offset[i+dim_counter_offset] = this->rlayout->starts[i];
+            dims[i+dim_counter_offset] = this->rlayout->sizes[i];
+            memshape[i+dim_counter_offset] = this->rmemlayout->subsizes[i];
+            memoffset[i+dim_counter_offset] = 0;
+        }
+        mspace = H5Screate_simple(dim_counter_offset + ndim(fc), memshape, NULL);
+        H5Sselect_hyperslab(mspace, H5S_SELECT_SET, memoffset, NULL, count, NULL);
+    }
+    else
+    {
+        for (unsigned int i=0; i<ndim(fc); i++)
+        {
+            count[i+dim_counter_offset] = this->clayout->subsizes[i];
+            offset[i+dim_counter_offset] = this->clayout->starts[i];
+            dims[i+dim_counter_offset] = this->clayout->sizes[i];
+            memshape[i+dim_counter_offset] = count[i+dim_counter_offset];
+            memoffset[i+dim_counter_offset] = 0;
+        }
+        mspace = H5Screate_simple(dim_counter_offset + ndim(fc), memshape, NULL);
+        H5Sselect_hyperslab(mspace, H5S_SELECT_SET, memoffset, NULL, count, NULL);
+    }
+
+    /* open/create data set */
+    if (read)
+        fspace = H5Dget_space(dset_id);
+    else
+    {
+        if (!H5Lexists(file_id, field_name.c_str(), H5P_DEFAULT))
+            H5Gcreate(
+                    file_id, field_name.c_str(),
+                    H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+        if (H5Lexists(file_id, dset_name.c_str(), H5P_DEFAULT))
+        {
+            dset_id = H5Dopen(file_id, dset_name.c_str(), H5P_DEFAULT);
+            fspace = H5Dget_space(dset_id);
+        }
+        else
+        {
+            fspace = H5Screate_simple(
+                    ndim(fc),
+                    dims,
+                    NULL);
+            /* chunking needs to go in here */
+            dset_id = H5Dcreate(
+                    file_id,
+                    dset_name.c_str(),
+                    (this->real_space_representation ? this->rnumber_H5T : this->cnumber_H5T),
+                    fspace,
+                    H5P_DEFAULT,
+                    H5P_DEFAULT,
+                    H5P_DEFAULT);
+            assert(dset_id > 0);
+        }
+    }
+    /* both dset_id and fspace should now have sane values */
+
+    if (read)
+        *this = rnumber(0.0);
+
+    /* check file space */
+    int ndims_fspace = H5Sget_simple_extent_dims(fspace, dims, NULL);
+    variable_used_only_in_assert(ndims_fspace);
+    assert(ndims_fspace == int(ndim(fc) + 1));
+    offset[0] = toffset;
+    if (this->real_space_representation)
+    {
+        for (unsigned int i=0; i<ndim(fc); i++)
+        {
+            offset[i+dim_counter_offset] = this->rlayout->starts[i];
+            assert(dims[i+dim_counter_offset] == this->rlayout->sizes[i]);
+        }
+        H5Sselect_hyperslab(fspace, H5S_SELECT_SET, offset, NULL, count, NULL);
+        if (read)
+        {
+            H5Dread(dset_id, this->rnumber_H5T, mspace, fspace, H5P_DEFAULT, this->data);
+            this->real_space_representation = true;
+        }
+        else
+        {
+            assert(this->real_space_representation);
+            H5Dwrite(dset_id, this->rnumber_H5T, mspace, fspace, H5P_DEFAULT, this->data);
+        }
+        H5Sclose(mspace);
+    }
+    else
+    {
+        for (unsigned int i=0; i<ndim(fc); i++)
+        {
+            offset[i+dim_counter_offset] = this->clayout->starts[i];
+            assert(dims[i+dim_counter_offset] == this->clayout->sizes[i]);
+        }
+        H5Sselect_hyperslab(fspace, H5S_SELECT_SET, offset, NULL, count, NULL);
+        if (read)
+        {
+            H5Dread(dset_id, this->cnumber_H5T, mspace, fspace, H5P_DEFAULT, this->data);
+            this->real_space_representation = false;
+            this->symmetrize();
+        }
+        else
+        {
+            assert(!this->real_space_representation);
+            H5Dwrite(dset_id, this->cnumber_H5T, mspace, fspace, H5P_DEFAULT, this->data);
+        }
+        H5Sclose(mspace);
+    }
+
+    H5Sclose(fspace);
+    /* close data set */
+    H5Dclose(dset_id);
+    /* close file */
+    H5Fclose(file_id);
+    return EXIT_SUCCESS;
+}
+
+
+template <typename rnumber,
+          field_backend be,
+          field_components fc>
+int field<rnumber, be, fc>::write_0slice(
+        const hid_t group,
+        const std::string field_name,
+        const int iteration)
+{
+    start_mpi_profiling_zone(turtle_mpi_pcontrol::HDF5);
+    // this should in principle work for any fc
+    TIMEZONE("field::write_0slice");
+    assert(this->real_space_representation);
+    if (this->myrank == 0)
+    {
+        hid_t dset, wspace, mspace;
+        int ndims;
+        hsize_t count[5], offset[5], dims[5];
+        offset[0] = iteration;
+        offset[1] = 0;
+        offset[2] = 0;
+        offset[3] = 0;
+        offset[4] = 0;
+        dset = H5Dopen(
+                group,
+                ("0slices/" + field_name + "/real").c_str(),
+                H5P_DEFAULT);
+        wspace = H5Dget_space(dset);
+        ndims = H5Sget_simple_extent_dims(wspace, dims, NULL);
+        // array in memory has 2 extra x points, because FFTW
+        count[0] = 1;
+        count[1] = this->rmemlayout->sizes[1];
+        count[2] = this->rmemlayout->sizes[2];
+        count[3] = 3;
+        count[4] = 3;
+        mspace = H5Screate_simple(ndims, count, NULL);
+        // array in file should not have the extra 2 points
+        count[1] = this->rlayout->sizes[1];
+        count[2] = this->rlayout->sizes[2];
+        // select right slice in file
+        H5Sselect_hyperslab(
+            wspace,
+            H5S_SELECT_SET,
+            offset,
+            NULL,
+            count,
+            NULL);
+        offset[0] = 0;
+        // select proper regions of memory
+        H5Sselect_hyperslab(
+            mspace,
+            H5S_SELECT_SET,
+            offset,
+            NULL,
+            count,
+            NULL);
+        H5Dwrite(
+            dset,
+            this->rnumber_H5T,
+            mspace,
+            wspace,
+            H5P_DEFAULT,
+            this->data);
+        H5Dclose(dset);
+        H5Sclose(mspace);
+        H5Sclose(wspace);
+    }
+    finish_mpi_profiling_zone(turtle_mpi_pcontrol::HDF5);
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber,
+          field_backend be,
+          field_components fc>
+int field<rnumber, be, fc>::write_filtered(
+        const std::string fname,
+        const std::string field_name,
+        const int iteration,
+        int nx,
+        int ny,
+        int nz)
+{
+    /* file dataset has same dimensions as field */
+    TIMEZONE("field::write_filtered");
+    // only works in Fourier representation
+    assert(!this->real_space_representation);
+    assert(hsize_t(nx) <= this->rlayout->sizes[2]);
+    assert(hsize_t(ny) <= this->rlayout->sizes[1]);
+    assert(hsize_t(nz) <= this->rlayout->sizes[0]);
+    // current algorithm only works for more than one process
+    assert(this->nprocs >= 2);
+    hid_t file_id, dset_id, plist_id;
+    dset_id = H5I_BADID;
+    std::string dset_name = (
+            "/" + field_name +
+            "/complex" +
+            "/" + std::to_string(iteration));
+
+    /* open/create file */
+    plist_id = H5Pcreate(H5P_FILE_ACCESS);
+    H5Pset_fapl_mpio(plist_id, this->comm, MPI_INFO_NULL);
+    bool file_exists = false;
+    struct stat file_buffer;
+    file_exists = (stat(fname.c_str(), &file_buffer) == 0);
+    if (file_exists)
+        file_id = H5Fopen(fname.c_str(), H5F_ACC_RDWR, plist_id);
+    else
+        file_id = H5Fcreate(fname.c_str(), H5F_ACC_EXCL, H5P_DEFAULT, plist_id);
+    assert(file_id >= 0);
+    H5Pclose(plist_id);
+
+    /* generic space initialization */
+    hid_t fspace, mspace;
+    hsize_t count[ndim(fc)], offset[ndim(fc)], dims[ndim(fc)], fdims[ndim(fc)];
+    hsize_t memoffset[ndim(fc)], memshape[ndim(fc)];
+
+    // set up dimensions
+    for (unsigned int i=3; i<ndim(fc); i++)
+    {
+        count [i] = this->clayout->subsizes[i];
+        offset[i] = this->clayout->starts[i];
+        dims  [i] = this->clayout->sizes[i];
+        memshape [i] = count[i];
+        memoffset[i] = 0;
+    }
+    // these are dimensions of dataset, needed
+    // to create dataset
+    dims[0] = ny;
+    dims[1] = nz;
+    dims[2] = nx/2+1;
+
+    /* open/create data set */
+    if (!H5Lexists(file_id, field_name.c_str(), H5P_DEFAULT))
+    {
+        hid_t gid_tmp = H5Gcreate(
+                file_id, field_name.c_str(),
+                H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+        H5Gclose(gid_tmp);
+    }
+    if (!H5Lexists(file_id, (field_name + "/complex").c_str(), H5P_DEFAULT))
+    {
+        hid_t gid_tmp = H5Gcreate(
+                file_id, ("/" + field_name + "/complex").c_str(),
+                H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+        H5Gclose(gid_tmp);
+    }
+    if (H5Lexists(file_id, dset_name.c_str(), H5P_DEFAULT))
+    {
+        dset_id = H5Dopen(file_id, dset_name.c_str(), H5P_DEFAULT);
+        fspace = H5Dget_space(dset_id);
+    }
+    else
+    {
+        if (nz == 1)
+        {
+            hsize_t temp_dims[ndim(fc)-1];
+            temp_dims[0] = dims[0];
+            for (unsigned int i=1; i<ndim(fc)-1; i++)
+                temp_dims[i] = dims[i+1];
+            fspace = H5Screate_simple(
+                    ndim(fc)-1,
+                    temp_dims,
+                    NULL);
+        }
+        else
+            fspace = H5Screate_simple(
+                    ndim(fc),
+                    dims,
+                    NULL);
+        /* chunking needs to go in here */
+        dset_id = H5Dcreate(
+                file_id,
+                dset_name.c_str(),
+                this->cnumber_H5T,
+                fspace,
+                H5P_DEFAULT,
+                H5P_DEFAULT,
+                H5P_DEFAULT);
+        assert(dset_id > 0);
+    }
+    /* check file space */
+    int ndims_fspace = H5Sget_simple_extent_dims(fspace, fdims, NULL);
+    variable_used_only_in_assert(ndims_fspace);
+    if (nz == 1)
+    {
+        assert(((unsigned int)(ndims_fspace)) == ndim(fc)-1);
+        assert(dims[0] == fdims[0]);
+        for (unsigned int i=1; i<ndim(fc)-1; i++)
+            assert(dims[i+1] == fdims[i]);
+    }
+    else
+    {
+        assert(((unsigned int)(ndims_fspace)) == ndim(fc));
+        for (unsigned int i=0; i<ndim(fc); i++)
+        {
+            assert(dims[i] == fdims[i]);
+        }
+    }
+    /* both dset_id and fspace now have sane values */
+
+    /// set up counts and offsets
+    /// x is easy, since only positive modes are present
+    count [2] = nx/2+1;
+    offset[2] = 0;
+    memshape [2] = this->clayout->subsizes[2];
+    memoffset[2] = 0;
+
+    /// three options for y:
+    /// this->starts[0] <= ny/2
+    /// ny / 2 < this->starts[0] +this->clayout->subsizes[0] < this->sizes[0] - ny/2
+    /// this->starts[0] >= this->sizes[0] - ny/2
+    /// we don't care about saving the ny/2 mode, because of symmetry
+    hsize_t y0 = this->clayout->starts[0];
+    hsize_t y1 = this->clayout->starts[0] + this->clayout->subsizes[0];
+    memshape[0] = this->clayout->subsizes[0];
+    if (y1 <= hsize_t(ny/2))
+    {
+        count[0] = this->clayout->subsizes[0];
+        offset[0] = y0;
+        memoffset[0] = 0;
+    }
+    else
+    {
+        if (y0 < hsize_t(ny)/2)
+        {
+            count[0] = ny/2 - y0;
+            offset[0] = y0;
+            memoffset[0] = 0;
+        }
+        else
+        {
+            if (y1 <= hsize_t(this->clayout->sizes[0] - ny/2 + 1))
+            { // y0 < y1 therefore y0 <= this->clayout->sizes[0] - ny/2
+                count[0] = 0;
+                offset[0] = ny/2;
+                memoffset[0] = 0;
+            }
+            else
+            {
+                if (y0 <= hsize_t(this->clayout->sizes[0] - ny/2))
+                {
+                    count[0] = y1 - (this->clayout->sizes[0] - ny/2);
+                    offset[0] = ny - (this->clayout->sizes[0] - y0);
+                    memoffset[0] = this->clayout->subsizes[0] - count[0];
+                }
+                else
+                {
+                    count[0] = this->clayout->subsizes[0];
+                    offset[0] = ny - (this->clayout->sizes[0] - y0);
+                    memoffset[0] = 0;
+                }
+            }
+        }
+    }
+    if (nz>=2)
+    {
+        assert(nz%2==0);
+        /// for z, we need to take into account that there are
+        /// both positive and negative modes
+        for (int cz = 0; cz < 2; cz++)
+        {
+            count [1] = nz/2;
+            offset[1] = cz*nz/2;
+            memshape [1] = this->clayout->sizes[1];
+            memoffset[1] = cz*(this->clayout->sizes[1] - nz/2);
+
+            //now write data
+            mspace = H5Screate_simple(ndim(fc), memshape, NULL);
+            H5Sselect_hyperslab(mspace, H5S_SELECT_SET, memoffset, NULL, count, NULL);
+            H5Sselect_hyperslab(fspace, H5S_SELECT_SET, offset, NULL, count, NULL);
+            H5Dwrite(dset_id, this->cnumber_H5T, mspace, fspace, H5P_DEFAULT, this->data);
+            H5Sclose(mspace);
+        }
+    }
+    else
+    {
+        assert(nz == 1);
+        count [1] = 1;
+        offset[1] = 0;
+        memshape [1] = this->clayout->sizes[1];
+        memoffset[1] = 0;
+
+        //now write data
+        mspace = H5Screate_simple(ndim(fc), memshape, NULL);
+        H5Sselect_hyperslab(mspace, H5S_SELECT_SET, memoffset, NULL, count, NULL);
+        for (unsigned int i=1; i<ndim(fc)-1; i++)
+            count[i] = count[i+1];
+        H5Sselect_hyperslab(fspace, H5S_SELECT_SET, offset, NULL, count, NULL);
+        H5Dwrite(dset_id, this->cnumber_H5T, mspace, fspace, H5P_DEFAULT, this->data);
+        H5Sclose(mspace);
+    }
+
+
+    /* close file data space */
+    H5Sclose(fspace);
+    /* close data set */
+    H5Dclose(dset_id);
+    /* close file */
+    H5Fclose(file_id);
+    return EXIT_SUCCESS;
+}
+
diff --git a/cpp/full_code/NSE.cpp b/cpp/full_code/NSE.cpp
index 950883094c9af0c9301a645de241b6828a8c34b0..361de8c14f7d85b76d3a331546209c54b5fe78d4 100644
--- a/cpp/full_code/NSE.cpp
+++ b/cpp/full_code/NSE.cpp
@@ -893,6 +893,7 @@ int NSE<rnumber>::get_cpressure(field<rnumber, FFTW, ONE> *destination)
                   / k2);
             });
     destination->symmetrize();
+    this->kk->template dealias<rnumber, ONE>(destination->get_cdata());
     delete temp_scalar;
     delete rvelocity;
     return EXIT_SUCCESS;
diff --git a/cpp/full_code/NSVE.cpp b/cpp/full_code/NSVE.cpp
index 8986842be09db5d8c8f53bd4d912a45395149598..2aeb03f7a8c1903326067f2d3e9ea1106607a091 100644
--- a/cpp/full_code/NSVE.cpp
+++ b/cpp/full_code/NSVE.cpp
@@ -89,29 +89,37 @@ int NSVE<rnumber>::initialize(void)
     this->fs->iteration = this->iteration;
     this->fs->checkpoint = this->checkpoint;
 
+    this->fs->cvorticity->real_space_representation = false;
     if (this->iteration == 0) {
         this->fs->kk->store(stat_file);
     }
-
-    this->fs->cvorticity->real_space_representation = false;
-    if ((this->iteration == 0)
-     && (this->field_random_seed != 0)) {
-        // generate initial condition
-        make_gaussian_random_field(
-            this->fs->kk,
-            this->fs->cvelocity,
-            this->field_random_seed,
-            this->fs->injection_rate,
-            1.0,                        // Lint
-            1.5 / this->fs->kk->kM,     // etaK
-            6.78,
-            0.40,
-            3./2.);
-        compute_curl(this->fs->kk, this->fs->cvelocity, this->fs->cvorticity);
-        this->fs->kk->template project_divfree<rnumber>(
-                this->fs->cvorticity->get_cdata());
-        this->fs->cvorticity->symmetrize();
-        this->fs->io_checkpoint(false);
+    if (this->iteration == 0) {
+        if (!hdf5_tools::field_exists(
+                    this->fs->get_current_fname(),
+                    "vorticity",
+                    "complex",
+                    0)) {
+            // generate initial condition
+            make_gaussian_random_field(
+                this->fs->kk,
+                this->fs->cvelocity,
+                this->field_random_seed,
+                this->fs->injection_rate,
+                1.0,                        // Lint
+                1.5 / this->fs->kk->kM,     // etaK
+                6.78,
+                0.40,
+                3./2.);
+            this->fs->kk->template project_divfree<rnumber>(
+                    this->fs->cvelocity->get_cdata());
+            compute_curl(this->fs->kk, this->fs->cvelocity, this->fs->cvorticity);
+            this->fs->kk->template low_pass<rnumber, THREE>(
+                    this->fs->cvorticity->get_cdata(), this->fs->kk->kM);
+            this->fs->cvorticity->symmetrize();
+            this->fs->io_checkpoint(false);
+        } else {
+            this->fs->io_checkpoint();
+        }
     } else {
         this->fs->io_checkpoint();
     }
diff --git a/cpp/full_code/static_field.cpp b/cpp/full_code/static_field.cpp
index 5d42ccbe80c1ac6b669a4043a60ff01e66c064e6..b4153eb002d8285193323f268ede2e05d05f8769 100644
--- a/cpp/full_code/static_field.cpp
+++ b/cpp/full_code/static_field.cpp
@@ -59,19 +59,50 @@ int static_field<rnumber>::initialize(void)
             this->comm,
             fftw_planner_string_to_flag[this->fftw_plan_rigor]);
     this->velocity->real_space_representation = false;
+
+    this->kk = new kspace<FFTW, SMOOTH>(
+            this->vorticity->clayout, this->dkx, this->dky, this->dkz);
+
     //read static vorticity field from iteration 0
     std::string checkpoint_fname = (
             std::string(this->simname) +
             std::string("_checkpoint_") +
             std::to_string(0) +
             std::string(".h5"));
-    this->vorticity->io(
-	    checkpoint_fname,
-	    "vorticity",
-	    0,
-	    true);
-    this->kk = new kspace<FFTW, SMOOTH>(
-            this->vorticity->clayout, this->dkx, this->dky, this->dkz);
+    if (this->iteration == 0) {
+        if (!hdf5_tools::field_exists(
+                    checkpoint_fname,
+                    "vorticity",
+                    "complex",
+                    0)) {
+            // generate initial condition
+            make_gaussian_random_field(
+                this->kk,
+                this->vorticity,
+                this->field_random_seed,
+                0.4,                        // injection rate
+                1.0,                        // Lint
+                1.5 / this->kk->kM,     // etaK
+                6.78,
+                0.40,
+                3./2.);
+            this->vorticity->symmetrize();
+            this->kk->template project_divfree<rnumber>(
+                    this->vorticity->get_cdata());
+        } else {
+            this->vorticity->io(
+	            checkpoint_fname,
+	            "vorticity",
+	            0,
+	            true);
+        }
+    } else {
+        this->vorticity->io(
+	        checkpoint_fname,
+	        "vorticity",
+	        0,
+	        true);
+    }
 
     // compute the velocity field and store
     invert_curl(
@@ -223,6 +254,7 @@ int static_field<rnumber>::read_parameters(void)
     this->tracers0_integration_steps = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_integration_steps");
     this->tracers0_neighbours = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_neighbours");
     this->tracers0_smoothness = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_smoothness");
+    this->field_random_seed = hdf5_tools::read_value<int>(parameter_file, "parameters/field_random_seed");
     H5Fclose(parameter_file);
     MPI_Barrier(this->comm);
     return EXIT_SUCCESS;
diff --git a/cpp/full_code/static_field.hpp b/cpp/full_code/static_field.hpp
index 1cf40e0c15c98a7458c0b68cebaf1f030807daca..7e7eedfd07168921ad4ba819fad30a7538dca48c 100644
--- a/cpp/full_code/static_field.hpp
+++ b/cpp/full_code/static_field.hpp
@@ -49,6 +49,7 @@ class static_field: public direct_numerical_simulation
         int tracers0_integration_steps;
         int tracers0_neighbours;
         int tracers0_smoothness;
+        int field_random_seed;
 
         /* other stuff */
         kspace<FFTW, SMOOTH> *kk;
diff --git a/cpp/hdf5_tools.cpp b/cpp/hdf5_tools.cpp
index 9653fe9209169e985ef85418ede48ebaeec2a61d..3f4ae9c984098b1c351fc865325d8ea10c36a0d8 100644
--- a/cpp/hdf5_tools.cpp
+++ b/cpp/hdf5_tools.cpp
@@ -1,6 +1,8 @@
 #include "hdf5_tools.hpp"
 #include "scope_timer.hpp"
 
+#include <sys/stat.h>
+
 template <> hid_t hdf5_tools::hdf5_type_id<char>()
 {
     return H5T_NATIVE_CHAR;
@@ -602,6 +604,42 @@ int hdf5_tools::gather_and_write_with_single_rank(
     return EXIT_SUCCESS;
 }
 
+bool hdf5_tools::field_exists(
+            const std::string file_name,
+            const std::string field_name,
+            const std::string representation,
+            const int iteration)
+{
+    /* file dataset has same dimensions as field */
+    TIMEZONE("hdf5_tools::field_exists");
+    hid_t file_id = H5I_BADID;
+    bool all_good = true;
+
+    std::string dset_name = (
+            "/" + field_name +
+            "/" + representation +
+            "/" + std::to_string(iteration));
+
+    /* open/create file */
+    bool file_exists = false;
+    {
+        struct stat file_buffer;
+        file_exists = (stat(file_name.c_str(), &file_buffer) == 0);
+    }
+    all_good = all_good && file_exists;
+    if (not file_exists)
+        return false;
+    file_id = H5Fopen(file_name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+    all_good = all_good && (file_id >= 0);
+    if (file_id >= 0) {
+        all_good = all_good && (H5Lexists(file_id, dset_name.c_str(), H5P_DEFAULT));
+    }
+
+    if (file_id >= 0)
+        H5Fclose(file_id);
+    return all_good;
+}
+
 template
 std::vector<int> hdf5_tools::read_vector<int>(
         const hid_t,
diff --git a/cpp/hdf5_tools.hpp b/cpp/hdf5_tools.hpp
index d9e1bce1cd4c4d6d6c51bb44d6072617ab75e5f7..6be9150fd7c991d4412e4518156d6d29c59596ad 100644
--- a/cpp/hdf5_tools.hpp
+++ b/cpp/hdf5_tools.hpp
@@ -41,6 +41,12 @@ namespace hdf5_tools
             hid_t dset,
             int tincrement);
 
+    bool field_exists(
+            const std::string file_name,
+            const std::string field_name,
+            const std::string representation,
+            const int iteration);
+
     herr_t grow_dataset_visitor(
         hid_t o_id,
         const char *name,
diff --git a/cpp/vorticity_equation.cpp b/cpp/vorticity_equation.cpp
index 496b4434746489f9e492428d82dd9dc67da7df29..043245f8da152dac9df7b11d3e94f8ef2d8748ec 100644
--- a/cpp/vorticity_equation.cpp
+++ b/cpp/vorticity_equation.cpp
@@ -691,7 +691,7 @@ void vorticity_equation<rnumber, be>::compute_pressure(field<rnumber, be, ONE> *
             this->v[1]->rval(rindex,cc) = this->u->rval(rindex,cc)*this->u->rval(rindex,cc);
         });
     this->v[1]->dft();
-    this->kk->template dealias<rnumber, THREE>(this->v[1]->get_cdata());
+    //this->kk->template dealias<rnumber, THREE>(this->v[1]->get_cdata());
     this->kk->CLOOP(
                 [&](const ptrdiff_t cindex,
                     const hsize_t xindex,
@@ -719,7 +719,7 @@ void vorticity_equation<rnumber, be>::compute_pressure(field<rnumber, be, ONE> *
             this->v[1]->rval(rindex,cc) = this->u->rval(rindex,cc)*this->u->rval(rindex,(cc+1)%3);
         });
     this->v[1]->dft();
-    this->kk->template dealias<rnumber, THREE>(this->v[1]->get_cdata());
+    //this->kk->template dealias<rnumber, THREE>(this->v[1]->get_cdata());
     this->kk->CLOOP(
                 [&](const ptrdiff_t cindex,
                     const hsize_t xindex,
diff --git a/examples/Lagrangian_acceleration/acc_check.py b/examples/Lagrangian_acceleration/acc_check.py
index 9c2eac495e418077ddd05cd0a912d4a18033c0d7..0051e0f78bc5cfbc47d24a3a62679687df5e94d1 100644
--- a/examples/Lagrangian_acceleration/acc_check.py
+++ b/examples/Lagrangian_acceleration/acc_check.py
@@ -83,6 +83,44 @@ class acc_test(TurTLE.PP):
         pars = {}
         return pars
 
+def plot_scalars(
+        s0, s1,
+        fname = 'tst.pdf'):
+    import matplotlib.pyplot as plt
+    f = plt.figure(figsize = (9, 9))
+    a = f.add_subplot(331)
+    a.set_axis_off()
+    a.imshow(s0[0])
+    a = f.add_subplot(332)
+    a.set_axis_off()
+    a.imshow(s1[0])
+    a = f.add_subplot(333)
+    a.set_axis_off()
+    a.imshow(s1[0] - s0[0])
+    a = f.add_subplot(334)
+    a.set_axis_off()
+    a.imshow(s0[:, 0])
+    a = f.add_subplot(335)
+    a.set_axis_off()
+    a.imshow(s1[:, 0])
+    a = f.add_subplot(336)
+    a.set_axis_off()
+    a.imshow(s1[:, 0] - s0[:, 0])
+    a = f.add_subplot(337)
+    a.set_axis_off()
+    a.imshow(s0[:, :, 0])
+    a = f.add_subplot(338)
+    a.set_axis_off()
+    a.imshow(s1[:, :, 0])
+    a = f.add_subplot(339)
+    a.set_axis_off()
+    a.imshow(s1[:, :, 0] - s0[:, :, 0])
+    print("error is ", np.max(np.abs(s1 - s0)))
+    f.tight_layout()
+    f.savefig(fname)
+    plt.close(f)
+    return None
+
 def main():
     if not os.path.exists('acceleration_check_data.h5'):
         cc = TurTLE.DNS()
@@ -94,6 +132,17 @@ def main():
              '--ntpp', '1',
              '--simname', 'acceleration_check_data']
             + sys.argv[1:])
+
+    # corresponding C++ code needs to be turned on for this to work
+    ######
+    #dfile = h5py.File('field_test.h5', 'r')
+    #pnse = dfile['nse/real/0'][()]
+    #pnsve = dfile['nsve/real/0'][()]
+    #plot_scalars(pnsve, pnse, 'tst_real.pdf')
+    #pnse = dfile['nse/complex/0'][()]
+    #pnsve = dfile['nsve/complex/0'][()]
+    #plot_scalars(np.abs(pnsve), np.abs(pnse), 'tst_complex.pdf')
+    #dfile.close()
     return None
 
 if __name__ == '__main__':
diff --git a/examples/Lagrangian_acceleration/acceleration_check.cpp b/examples/Lagrangian_acceleration/acceleration_check.cpp
index 5507130c92230559de64d44ace619cde31019ec3..5e9f014e0ff6e08a263135baf44cb6fa058dcde7 100644
--- a/examples/Lagrangian_acceleration/acceleration_check.cpp
+++ b/examples/Lagrangian_acceleration/acceleration_check.cpp
@@ -138,6 +138,7 @@ template <typename rnumber>
 int acceleration_check<rnumber>::work_on_current_iteration(void)
 {
     TIMEZONE("acceleration_check::work_on_current_iteration");
+    DEBUG_MSG("work on current iteration, iteration is %d\n", this->iteration);
 
     // temporary fields
     field<rnumber, FFTW, ONE>* pressure_NSE = new field<rnumber, FFTW, ONE>(
@@ -190,19 +191,37 @@ int acceleration_check<rnumber>::work_on_current_iteration(void)
     nsve->fs->compute_pressure(pressure_NSVE);
     pressure_NSE->ift();
     pressure_NSVE->ift();
+    //pressure_NSE->io("field_test.h5",
+    //                  "nse",
+    //                  this->iteration,
+    //                  false);
+    //pressure_NSVE->io("field_test.h5",
+    //                  "nsve",
+    //                  this->iteration,
+    //                  false);
     get_error(
         pressure_NSE,
         pressure_NSVE,
         pressure_max_diff,
         pressure_L2norm);
     const bool pressure_ok = (pressure_max_diff / pressure_L2norm) < 1e-3;
+    //pressure_NSE->dft();
+    //pressure_NSE->io("field_test.h5",
+    //                  "nse",
+    //                  this->iteration,
+    //                  false);
+    //pressure_NSVE->dft();
+    //pressure_NSVE->io("field_test.h5",
+    //                  "nsve",
+    //                  this->iteration,
+    //                  false);
     /***************************/
 
     /***************************/
     // fixed energy injection rate
     // NSE
     nse->forcing_type = "fixed_energy_injection_rate";
-    DEBUG_MSG("calling acceleration for NSE");
+    DEBUG_MSG("calling acceleration for NSE\n");
     nse->get_cacceleration(acc_NSE);
     // NSVE
     nsve->fs->forcing_type = "fixed_energy_injection_rate";
@@ -246,11 +265,12 @@ int acceleration_check<rnumber>::work_on_current_iteration(void)
     if (acceleration_ok_feir && acceleration_ok_linear && pressure_ok)
         return EXIT_SUCCESS;
     else {
+        DEBUG_MSG("pressure_ok = %d, max_error = %g\n",
+                pressure_ok, pressure_max_diff);
         DEBUG_MSG("acceleration_ok_feir = %d, max error = %g\n",
                 acceleration_ok_feir, err_max_feir);
         DEBUG_MSG("acceleration_ok_linear = %d, max error = %g\n",
                 acceleration_ok_linear, err_max_linear);
-        DEBUG_MSG("pressure_ok = %d\n", pressure_ok);
         return EXIT_FAILURE;
     }
 }