diff --git a/bfps/DNS.py b/bfps/DNS.py
index 538285e519beae671d6612f0cf3015814a3a2417..90a42c795e2b3ccd6661f3c7ba9c77dad0901664 100644
--- a/bfps/DNS.py
+++ b/bfps/DNS.py
@@ -247,6 +247,7 @@ class DNS(_code):
             ii0 = iter0 // self.parameters['niter_stat']
             ii1 = iter1 // self.parameters['niter_stat']
             self.statistics['kshell'] = data_file['kspace/kshell'].value
+            self.statistics['nshell'] = data_file['kspace/nshell'].value
             for kk in [-1, -2]:
                 if (self.statistics['kshell'][kk] == 0):
                     self.statistics['kshell'][kk] = np.nan
@@ -282,13 +283,13 @@ class DNS(_code):
                     phi_ij[:, :, 1, 1] +
                     phi_ij[:, :, 2, 2])/2
                 pp_file['energy(t)'] = np.sum(energy_tk, axis = 1)
-                pp_file['energy(k)'] = np.mean(energy_tk, axis = 0)
+                pp_file['energy(k)'] = np.mean(energy_tk, axis = 0)*(4*np.pi*self.statistics['kshell']**2) / (self.statistics['dk']*self.statistics['nshell'])
                 enstrophy_tk = (
                     data_file['statistics/spectra/vorticity_vorticity'][ii0:ii1+1, :, 0, 0] +
                     data_file['statistics/spectra/vorticity_vorticity'][ii0:ii1+1, :, 1, 1] +
                     data_file['statistics/spectra/vorticity_vorticity'][ii0:ii1+1, :, 2, 2])/2
                 pp_file['enstrophy(t)'] = np.sum(enstrophy_tk, axis = 1)
-                pp_file['enstrophy(k)'] = np.mean(enstrophy_tk, axis = 0)
+                pp_file['enstrophy(k)'] = np.mean(enstrophy_tk, axis = 0)*(4*np.pi*self.statistics['kshell']**2) / self.statistics['nshell']
                 pp_file['vel_max(t)'] = data_file['statistics/moments/velocity'][ii0:ii1+1, 9, 3]
                 pp_file['renergy(t)'] = data_file['statistics/moments/velocity'][ii0:ii1+1, 2, 3]/2
             for k in ['t',
diff --git a/bfps/NavierStokes.py b/bfps/NavierStokes.py
index d158f5b2f07a1ae154fb75727b4db71e20b0566f..9c12ca201b928566545405c6ecaea12fb9da553d 100644
--- a/bfps/NavierStokes.py
+++ b/bfps/NavierStokes.py
@@ -598,6 +598,7 @@ class NavierStokes(_fluid_particle_base):
             ii0 = iter0 // self.parameters['niter_stat']
             ii1 = iter1 // self.parameters['niter_stat']
             self.statistics['kshell'] = data_file['kspace/kshell'].value
+            self.statistics['nshell'] = data_file['kspace/nshell'].value
             for kk in [-1, -2]:
                 if (self.statistics['kshell'][kk] == 0):
                     self.statistics['kshell'][kk] = np.nan
@@ -633,13 +634,13 @@ class NavierStokes(_fluid_particle_base):
                     phi_ij[:, :, 1, 1] +
                     phi_ij[:, :, 2, 2])/2
                 pp_file['energy(t)'] = np.sum(energy_tk, axis = 1)
-                pp_file['energy(k)'] = np.mean(energy_tk, axis = 0)
+                pp_file['energy(k)'] = np.mean(energy_tk, axis = 0)*(4*np.pi*self.statistics['kshell']**2) / (self.statistics['dk']*self.statistics['nshell'])
                 enstrophy_tk = (
                     data_file['statistics/spectra/vorticity_vorticity'][ii0:ii1+1, :, 0, 0] +
                     data_file['statistics/spectra/vorticity_vorticity'][ii0:ii1+1, :, 1, 1] +
                     data_file['statistics/spectra/vorticity_vorticity'][ii0:ii1+1, :, 2, 2])/2
                 pp_file['enstrophy(t)'] = np.sum(enstrophy_tk, axis = 1)
-                pp_file['enstrophy(k)'] = np.mean(enstrophy_tk, axis = 0)
+                pp_file['enstrophy(k)'] = np.mean(enstrophy_tk, axis = 0)*(4*np.pi*self.statistics['kshell']**2) / (self.statistics['dk']*self.statistics['nshell'])
                 pp_file['vel_max(t)'] = data_file['statistics/moments/velocity'][ii0:ii1+1, 9, 3]
                 pp_file['renergy(t)'] = data_file['statistics/moments/velocity'][ii0:ii1+1, 2, 3]/2
                 if 'trS2_Q_R' in data_file['statistics/moments'].keys():