diff --git a/bfps/NSVorticityEquation.py b/bfps/NSVorticityEquation.py
index e3a3c2c376a9b706e19e4ddb54baa182234e3dd2..50e62275111d55d0237cf902a6d38f2a56189573 100644
--- a/bfps/NSVorticityEquation.py
+++ b/bfps/NSVorticityEquation.py
@@ -217,6 +217,42 @@ class NSVorticityEquation(_fluid_particle_base):
             field_H5T = 'H5T_NATIVE_FLOAT'
         elif self.dtype == np.float64:
             field_H5T = 'H5T_NATIVE_DOUBLE'
+        self.variables += 'int checkpoint;\n'
+        self.read_checkpoint = """
+                //begincpp
+                if (myrank == 0)
+                {
+                    hid_t dset = H5Dopen(stat_file, "checkpoint", H5P_DEFAULT);
+                    H5Dread(
+                        dset,
+                        H5T_NATIVE_INT,
+                        H5S_ALL,
+                        H5S_ALL,
+                        H5P_DEFAULT,
+                        &checkpoint);
+                    H5Dclose(dset);
+                }
+                MPI_Bcast(&checkpoint, 1, MPI_INT, 0, MPI_COMM_WORLD);
+                fs->checkpoint = checkpoint;
+                //endcpp
+        """
+        self.store_checkpoint = """
+                //begincpp
+                checkpoint = fs->checkpoint;
+                if (myrank == 0)
+                {
+                    hid_t dset = H5Dopen(stat_file, "checkpoint", H5P_DEFAULT);
+                    H5Dwrite(
+                        dset,
+                        H5T_NATIVE_INT,
+                        H5S_ALL,
+                        H5S_ALL,
+                        H5P_DEFAULT,
+                        &checkpoint);
+                    H5Dclose(dset);
+                }
+                //endcpp
+        """
         self.fluid_start += """
                 //begincpp
                 char fname[512];
@@ -233,6 +269,7 @@ class NSVorticityEquation(_fluid_particle_base):
                         nx, ny, nz,
                         MPI_COMM_WORLD,
                         {1});
+                fs->checkpoints_per_file = checkpoints_per_file;
                 fs->nu = nu;
                 fs->fmode = fmode;
                 fs->famplitude = famplitude;
@@ -240,16 +277,25 @@ class NSVorticityEquation(_fluid_particle_base):
                 fs->fk1 = fk1;
                 strncpy(fs->forcing_type, forcing_type, 128);
                 fs->iteration = iteration;
+                {2}
                 fs->cvorticity->real_space_representation = false;
                 fs->io_checkpoint();
                 //endcpp
-                """.format(self.C_dtype, self.fftw_plan_rigor, field_H5T)
+                """.format(
+                        self.C_dtype,
+                        self.fftw_plan_rigor,
+                        self.read_checkpoint)
         self.fluid_start += self.store_kspace
         self.fluid_loop = 'fs->step(dt);\n'
         self.fluid_loop += ('if (fs->iteration % niter_out == 0)\n{\n' +
-                            self.fluid_output + '\n}\n')
+                            self.fluid_output +
+                            self.store_checkpoint +
+                            '\n}\n')
         self.fluid_end = ('if (fs->iteration % niter_out != 0)\n{\n' +
-                          self.fluid_output + '\n}\n' +
+                          self.fluid_output +
+                          self.store_checkpoint +
+                          'DEBUG_MSG("checkpoint value is %d\\n", checkpoint);\n' +
+                          '\n}\n' +
                           'delete fs;\n' +
                           'delete tmp_vec_field;\n' +
                           'delete tmp_scal_field;\n')
@@ -470,6 +516,7 @@ class NSVorticityEquation(_fluid_particle_base):
                                                  self.parameters['histogram_bins'],
                                                  4),
                                      dtype = np.int64)
+            ofile['checkpoint'] = int(0)
         return None
     def specific_parser_arguments(
             self,
@@ -580,7 +627,6 @@ class NSVorticityEquation(_fluid_particle_base):
                            spectra_slope = 2.0,
                            amplitude = 0.05)
                     f['vorticity/complex/{0}'.format(0)] = data
-                f['fields_stored'] = int(1)
                 f.close()
         self.run(
                 nb_processes = opt.nb_processes,
diff --git a/bfps/_code.py b/bfps/_code.py
index 2f21c16abf9aaff1e6a83e59bf9070a83b61ec0e..3ff3de341b20c578262c1111835ba799359db1f6 100644
--- a/bfps/_code.py
+++ b/bfps/_code.py
@@ -118,8 +118,15 @@ class _code(_base):
                     sprintf(fname, "%s.h5", simname);
                     parameter_file = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
                     Cdset = H5Dopen(parameter_file, "iteration", H5P_DEFAULT);
-                    H5Dread(Cdset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &iteration);
-                    DEBUG_MSG("simname is %s and iteration is %d\\n", simname, iteration);
+                    H5Dread(
+                        Cdset,
+                        H5T_NATIVE_INT,
+                        H5S_ALL,
+                        H5S_ALL,
+                        H5P_DEFAULT,
+                        &iteration);
+                    DEBUG_MSG("simname is %s and iteration is %d\\n",
+                              simname, iteration);
                     H5Dclose(Cdset);
                     H5Fclose(parameter_file);
                     read_parameters();
diff --git a/bfps/cpp/vorticity_equation.cpp b/bfps/cpp/vorticity_equation.cpp
index f40c105c2eedf76be5495357418ac74600a04399..87bcf331ff4d6b6bf215fdc34f713bc9372e87b3 100644
--- a/bfps/cpp/vorticity_equation.cpp
+++ b/bfps/cpp/vorticity_equation.cpp
@@ -46,6 +46,58 @@ void vorticity_equation<rnumber, be>::impose_zero_modes()
     this->v[2]->impose_zero_mode();
 }
 
+template <class rnumber,
+          field_backend be>
+void vorticity_equation<rnumber, be>::update_checkpoint()
+{
+    std::string fname = this->get_current_fname();
+    bool file_exists = false;
+    {
+        struct stat file_buffer;
+        file_exists = (stat(fname.c_str(), &file_buffer) == 0);
+    }
+    if (file_exists)
+    {
+        // check how many fields there are in the checkpoint file
+        // increment checkpoint if needed
+        hsize_t fields_stored;
+        hid_t fid, group_id;
+        fid = H5Fopen(fname.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+        group_id = H5Gopen(fid, "vorticity/complex", H5P_DEFAULT);
+        H5Gget_num_objs(
+                group_id,
+                &fields_stored);
+        H5Gclose(group_id);
+        H5Fclose(fid);
+        if (fields_stored >= this->checkpoints_per_file)
+            this->checkpoint++;
+    }
+    else if (this->cvelocity->myrank == 0)
+    {
+        // create file, create fields_stored dset
+        hid_t fid = H5Fcreate(
+                fname.c_str(),
+                H5F_ACC_EXCL,
+                H5P_DEFAULT,
+                H5P_DEFAULT);
+        hid_t gg = H5Gcreate(
+                fid,
+                "vorticity",
+                H5P_DEFAULT,
+                H5P_DEFAULT,
+                H5P_DEFAULT);
+        hid_t ggg = H5Gcreate(
+                gg,
+                "complex",
+                H5P_DEFAULT,
+                H5P_DEFAULT,
+                H5P_DEFAULT);
+        H5Gclose(ggg);
+        H5Gclose(gg);
+        H5Fclose(fid);
+    }
+}
+
 template <class rnumber,
           field_backend be>
 vorticity_equation<rnumber, be>::vorticity_equation(
diff --git a/bfps/cpp/vorticity_equation.hpp b/bfps/cpp/vorticity_equation.hpp
index aba96d5e2b08149636486e880c7d901eea13bb26..6c2621c98c15179b61f4cf52f6a6a227847574a4 100644
--- a/bfps/cpp/vorticity_equation.hpp
+++ b/bfps/cpp/vorticity_equation.hpp
@@ -103,68 +103,10 @@ class vorticity_equation
                     std::to_string(this->checkpoint) +
                     std::string(".h5"));
         }
-        inline void update_checkpoint()
-        {
-            std::string fname = this->get_current_fname();
-            bool file_exists = false;
-            {
-                struct stat file_buffer;
-                file_exists = (stat(fname.c_str(), &file_buffer) == 0);
-            }
-            if (file_exists)
-            {
-                // check how many fields there are in the checkpoint file
-                // increment checkpoint if needed
-                int fields_stored;
-                hid_t fid, dset_id;
-                fid = H5Fopen(fname.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
-                dset_id = H5Dopen(fid, "fields_stored", H5P_DEFAULT);
-                H5Dread(dset_id,
-                        H5T_NATIVE_INT,
-                        H5S_ALL, H5S_ALL,
-                        H5P_DEFAULT,
-                        &fields_stored);
-                H5Dclose(dset_id);
-                H5Fclose(fid);
-                if (fields_stored >= this->checkpoints_per_file)
-                    this->checkpoint++;
-                else
-                {
-                    // update fields_stored dset
-                    fields_stored++;
-                    if (this->cvelocity->myrank == 0)
-                    {
-                        fid = H5Fopen(fname.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
-                        dset_id = H5Dopen(fid, "fields_stored", H5P_DEFAULT);
-                        H5Dwrite(dset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &fields_stored);
-                        H5Dclose(dset_id);
-                        H5Fclose(fid);
-                    }
-                }
-            }
-            else if (this->cvelocity->myrank == 0)
-            {
-                // create file, create fields_stored dset
-                hid_t fid = H5Fcreate(fname.c_str(), H5F_ACC_EXCL, H5P_DEFAULT, H5P_DEFAULT);
-                hsize_t one[] = {1};
-                hid_t fspace = H5Screate_simple(
-                        1,
-                        one,
-                        NULL);
-                hid_t dset = H5Dcreate(
-                        fid,
-                        "fields_stored",
-                        H5T_NATIVE_INT,
-                        fspace,
-                        H5P_DEFAULT,
-                        H5P_DEFAULT,
-                        H5P_DEFAULT);
-                H5Dclose(dset);
-                H5Fclose(fid);
-            }
-        }
+        void update_checkpoint(void);
         inline void io_checkpoint(bool read = true)
         {
+            assert(!this->cvorticity->real_space_representation);
             if (!read)
                 this->update_checkpoint();
             std::string fname = this->get_current_fname();