diff --git a/TurTLE/DNS.py b/TurTLE/DNS.py
index 3aea17babb57b56e02f9bdb5a3635caff925cbf5..ffa839e29f77b4afd3592ae7c00cce2148a6b0fb 100644
--- a/TurTLE/DNS.py
+++ b/TurTLE/DNS.py
@@ -1098,6 +1098,7 @@ class DNS(_code):
                         particle_file.create_group('tracers0/velocity_gradient')
                     if self.dns_type in ['NSVE_Stokes_particles']:
                         particle_file.create_group('tracers0/momentum')
+                        particle_file.create_group('tracers0/vorticity')
                     if self.dns_type in ['NSVEp_extra_sampling']:
                         particle_file.create_group('tracers0/velocity_gradient')
                         particle_file.create_group('tracers0/pressure')
diff --git a/cpp/full_code/NSVE_Stokes_particles.cpp b/cpp/full_code/NSVE_Stokes_particles.cpp
index a3e4dd381c4b2e11e1b5963407f0dd91d2a34793..37e34604fc71c7331c721f8688bf39665766e7c5 100644
--- a/cpp/full_code/NSVE_Stokes_particles.cpp
+++ b/cpp/full_code/NSVE_Stokes_particles.cpp
@@ -183,6 +183,20 @@ int NSVE_Stokes_particles<rnumber>::do_stats()
             this->ps->getLocalNbParticles(),
             this->ps->get_step_idx()-1);
 
+    /// sample vorticity
+    std::fill_n(pdata1.get(), 3*this->ps->getLocalNbParticles(), 0);
+    *this->tmp_vec_field = this->fs->cvorticity->get_cdata();
+    this->tmp_vec_field->ift();
+    this->ps->sample_compute_field(*this->tmp_vec_field, pdata1.get());
+    this->particles_sample_writer_mpi->template save_dataset<3>(
+            "tracers0",
+            "vorticity",
+            pdata0.get(),
+            &pdata1,
+            this->ps->getParticlesIndexes(),
+            this->ps->getLocalNbParticles(),
+            this->ps->get_step_idx()-1);
+
     // deallocate temporary data array
     delete[] pdata0.release();
     delete[] pdata1.release();