diff --git a/bfps/NSVorticityEquation.py b/bfps/NSVorticityEquation.py
index de0b7012a0461369bffe4e4befe571a76eeb8944..21615f351ad3072bb263ffd8b86409fab283a763 100644
--- a/bfps/NSVorticityEquation.py
+++ b/bfps/NSVorticityEquation.py
@@ -148,7 +148,8 @@ class NSVorticityEquation(_fluid_particle_base):
                     std::string("/tracers0/rhs/")  + std::to_string(fs->iteration), // dataset name for initial input
                     tracers0_neighbours,        // parameter (interpolation no neighbours)
                     tracers0_smoothness,        // parameter
-                    MPI_COMM_WORLD);
+                    MPI_COMM_WORLD,
+                    fs->iteration+1);
             particles_output_hdf5<double,3,3> particles_output_writer_mpi(
                         MPI_COMM_WORLD,
                         "tracers0",
diff --git a/bfps/cpp/particles/particles_system.hpp b/bfps/cpp/particles/particles_system.hpp
index 48b4baba86c84346538b007f14bf8fb316de7cf2..05466fc9313ad5f86c8c8279a407f467ec61d8b1 100644
--- a/bfps/cpp/particles/particles_system.hpp
+++ b/bfps/cpp/particles/particles_system.hpp
@@ -48,7 +48,8 @@ public:
                      const field_rnumber* in_field_data, const std::array<size_t,3>& in_local_field_dims,
                      const std::array<size_t,3>& in_local_field_offset,
                      const std::array<size_t,3>& in_field_memory_dims,
-                     MPI_Comm in_mpi_com)
+                     MPI_Comm in_mpi_com,
+                     const int in_current_iteration = 1)
         : mpi_com(in_mpi_com),
           current_partition_interval({in_local_field_offset[IDX_Z], in_local_field_offset[IDX_Z] + in_local_field_dims[IDX_Z]}),
           partition_interval_size(current_partition_interval.second - current_partition_interval.first),
@@ -58,7 +59,7 @@ public:
                    interpolator, field, in_spatial_box_width, in_spatial_box_offset, in_spatial_partition_width),
           spatial_box_width(in_spatial_box_width), spatial_partition_width(in_spatial_partition_width),
           my_spatial_low_limit(in_my_spatial_low_limit), my_spatial_up_limit(in_my_spatial_up_limit),
-          my_nb_particles(0), step_idx(1){
+          my_nb_particles(0), step_idx(in_current_iteration){
 
         current_my_nb_particles_per_partition.reset(new int[partition_interval_size]);
         current_offset_particles_for_partition.reset(new int[partition_interval_size+1]);
diff --git a/bfps/cpp/particles/particles_system_builder.hpp b/bfps/cpp/particles/particles_system_builder.hpp
index f828343742cec3103942c1c690f95f3942fa3435..ab895d4b2d53ab562dba28e17e404297ab225700 100644
--- a/bfps/cpp/particles/particles_system_builder.hpp
+++ b/bfps/cpp/particles/particles_system_builder.hpp
@@ -118,7 +118,8 @@ struct particles_system_build_container {
              const int nparticles, // to check coherency between parameters and hdf input file
              const std::string& fname_input, // particles input filename
             const std::string& inDatanameState, const std::string& inDatanameRhs, // input dataset names
-             MPI_Comm mpi_comm){
+             MPI_Comm mpi_comm,
+            const int in_current_iteration){
 
         // The size of the field grid (global size) all_size seems
         std::array<size_t,3> field_grid_dim;
@@ -208,7 +209,8 @@ struct particles_system_build_container {
                                                                                                    local_field_dims,
                                                                                                    local_field_offset,
                                                                                                    local_field_mem_size,
-                                                                                                   mpi_comm);
+                                                                                                   mpi_comm,
+                                                                                                   in_current_iteration);
 
         // Load particles from hdf5
         particles_input_hdf5<particles_rnumber, 3,3> generator(mpi_comm, fname_input,
@@ -246,14 +248,15 @@ inline std::unique_ptr<abstract_particles_system<particles_rnumber>> particles_s
         const std::string& inDatanameState, const std::string& inDatanameRhs, // input dataset names
         const int interpolation_size,
         const int spline_mode,
-        MPI_Comm mpi_comm){
+        MPI_Comm mpi_comm,
+        const int in_current_iteration){
     return Template_double_for_if::evaluate<std::unique_ptr<abstract_particles_system<particles_rnumber>>,
                        int, 1, 7, 1, // interpolation_size
                        int, 0, 3, 1, // spline_mode
                        particles_system_build_container<field_rnumber,be,particles_rnumber>>(
                            interpolation_size, // template iterator 1
                            spline_mode, // template iterator 2
-                           fs_field,fs_kk, nsteps, nparticles, fname_input, inDatanameState, inDatanameRhs, mpi_comm);
+                           fs_field,fs_kk, nsteps, nparticles, fname_input, inDatanameState, inDatanameRhs, mpi_comm, in_current_iteration);
 }