diff --git a/cpp/full_code/kraichnan_field.cpp b/cpp/full_code/kraichnan_field.cpp
index aad4f4020c8b60cf0925abb70829594757af1746..459c9d5f8df2f6de3599f46437bad7c1aae05b21 100644
--- a/cpp/full_code/kraichnan_field.cpp
+++ b/cpp/full_code/kraichnan_field.cpp
@@ -121,7 +121,6 @@ int kraichnan_field<rnumber>::step(void)
         this->spectrum_coefficient * 3./2.); // incompressibility projection factor
     this->velocity->ift();
 
-    // PUT THE INCOMPRESSIBILITY PROJECTION FACTOR sqrt(2/3)
     DEBUG_MSG("L2Norm before: %g\n",
             this->velocity->L2norm(this->kk));
     this->kk->template project_divfree<rnumber>(this->velocity->get_cdata());
@@ -137,6 +136,7 @@ template <typename rnumber>
 int kraichnan_field<rnumber>::write_checkpoint(void)
 {
     TIMEZONE("kraichnan_file::write_checkpoint");
+    this->io_checkpoint(false);
     this->particles_output_writer_mpi->open_file(this->get_current_fname());
     this->particles_output_writer_mpi->template save<3>(
             this->ps->getParticlesState(),
@@ -145,6 +145,7 @@ int kraichnan_field<rnumber>::write_checkpoint(void)
             this->ps->getLocalNbParticles(),
             this->iteration);
     this->particles_output_writer_mpi->close_file();
+    this->write_iteration();
     return EXIT_SUCCESS;
 }
 
@@ -238,6 +239,66 @@ int kraichnan_field<rnumber>::read_parameters(void)
     return EXIT_SUCCESS;
 }
 
+template <class rnumber>
+void kraichnan_field<rnumber>::update_checkpoint()
+{
+    std::string fname = this->get_current_fname();
+    if (this->kk->layout->myrank == 0)
+    {
+        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, "velocity/real", H5P_DEFAULT);
+            H5Gget_num_objs(
+                    group_id,
+                    &fields_stored);
+            bool dset_exists = H5Lexists(
+                    group_id,
+                    std::to_string(this->iteration).c_str(),
+                    H5P_DEFAULT);
+            H5Gclose(group_id);
+            H5Fclose(fid);
+            if ((int(fields_stored) >= this->checkpoints_per_file) &&
+                !dset_exists)
+                this->checkpoint++;
+        }
+        else
+        {
+            // 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,
+                    "velocity",
+                    H5P_DEFAULT,
+                    H5P_DEFAULT,
+                    H5P_DEFAULT);
+            hid_t ggg = H5Gcreate(
+                    gg,
+                    "real",
+                    H5P_DEFAULT,
+                    H5P_DEFAULT,
+                    H5P_DEFAULT);
+            H5Gclose(ggg);
+            H5Gclose(gg);
+            H5Fclose(fid);
+        }
+    }
+    MPI_Bcast(&this->checkpoint, 1, MPI_INT, 0, this->kk->layout->comm);
+}
+
 template class kraichnan_field<float>;
 template class kraichnan_field<double>;
 
diff --git a/cpp/full_code/kraichnan_field.hpp b/cpp/full_code/kraichnan_field.hpp
index 3710a0f7698953e1404c1f9849fd9ee514adb052..992724fb998a69af3ac6eaee6bcdcfa440e13a38 100644
--- a/cpp/full_code/kraichnan_field.hpp
+++ b/cpp/full_code/kraichnan_field.hpp
@@ -86,6 +86,21 @@ class kraichnan_field: public direct_numerical_simulation
         virtual int read_parameters(void);
         int write_checkpoint(void);
         int do_stats(void);
+
+        
+        void update_checkpoint(void);
+        inline void io_checkpoint(bool read = true)
+        {
+            assert(this->velocity->real_space_representation);
+            if (!read)
+                this->update_checkpoint();
+            std::string fname = this->get_current_fname();
+            this->velocity->io(
+                    fname,
+                    "velocity",
+                    this->iteration,
+                    read);
+        }
 };
 
 #endif//KRAICHNAN_FIELD_HPP