diff --git a/contrib/postprocess_fluid_stats/postprocess_fluid_stats.f90 b/contrib/postprocess_fluid_stats/postprocess_fluid_stats.f90
index 9899f08e380..6b5dfb63816 100644
--- a/contrib/postprocess_fluid_stats/postprocess_fluid_stats.f90
+++ b/contrib/postprocess_fluid_stats/postprocess_fluid_stats.f90
@@ -2,14 +2,12 @@
 !! Martin Karp 27/01-23
 program postprocess_fluid_stats
   use neko
-  use mean_flow
   implicit none
 
-  character(len=NEKO_FNAME_LEN) :: inputchar, mesh_fname, stats_fname, mean_fname
-  type(file_t) :: mean_file, stats_file, output_file, mesh_file
+  character(len=NEKO_FNAME_LEN) :: inputchar, mesh_fname, stats_fname
+  type(file_t) :: stats_file, output_file, mesh_file
   real(kind=rp) :: start_time
-  type(fld_file_data_t) :: stats_data, mean_data
-  type(mean_flow_t) :: avg_flow
+  type(fld_file_data_t) :: stats_data
   type(fluid_stats_t) :: fld_stats
   type(coef_t) :: coef
   type(dofmap_t) :: dof
@@ -25,10 +23,12 @@ program postprocess_fluid_stats
 
   if ((argc .lt. 3) .or. (argc .gt. 3)) then
      if (pe_rank .eq. 0) then
-        write(*,*) 'Usage: ./postprocess_fluid_stats mesh.nmsh mean_field.fld stats.fld'
-        write(*,*) 'Example command: ./postprocess_fluid_stats mesh.nmsh mean_fieldblabla.fld statsblabla.fld'
-        write(*,*) 'Computes the statstics from the fld files described in mean_fielblabla.nek5000 statsblabla.nek5000'
-        write(*,*) 'Currently we output two new fld files reynolds and mean_vei_grad'
+        write(*,*) 'Usage: ./postprocess_fluid_stats mesh.nmsh stats.fld'
+        write(*,*) 'Example command: ./postprocess_fluid_stats mesh.nmsh statsblabla.fld'
+        write(*,*) 'Computes the statstics from the (3D) fld files described in statsblabla.nek5000'
+        write(*,*) 'Recommended to use PyNekTools for more advanced postprocessing.'
+        write(*,*) 'If this does not work, switch to that.'
+        write(*,*) 'Currently we output two new fld files reynolds and mean_vel_grad'
         write(*,*) 'In Reynolds the fields are ordered as:'
         write(*,*) 'x-velocity=<u`u`>'
         write(*,*) 'y-velocity=<v`v`>'
@@ -57,50 +57,45 @@ program postprocess_fluid_stats
   read(inputchar, *) mesh_fname
   mesh_file = file_t(trim(mesh_fname))
   call get_command_argument(2, inputchar)
-  read(inputchar, *) mean_fname
-  mean_file = file_t(trim(mean_fname))
-  call get_command_argument(3, inputchar)
   read(inputchar, *) stats_fname
   stats_file = file_t(trim(stats_fname))
 
   call mesh_file%read(msh)
 
-  call mean_data%init(msh%nelv,msh%offset_el)
   call stats_data%init(msh%nelv,msh%offset_el)
-  call mean_file%read(mean_data)
   call stats_file%read(stats_data)
 
   do i = 1,msh%nelv
-     lx = mean_data%lx
-     msh%elements(i)%e%pts(1)%p%x(1) = mean_data%x%x(linear_index(1,1,1,i,lx,lx,lx))
-     msh%elements(i)%e%pts(2)%p%x(1) = mean_data%x%x(linear_index(lx,1,1,i,lx,lx,lx))
-     msh%elements(i)%e%pts(3)%p%x(1) = mean_data%x%x(linear_index(1,lx,1,i,lx,lx,lx))
-     msh%elements(i)%e%pts(4)%p%x(1) = mean_data%x%x(linear_index(lx,lx,1,i,lx,lx,lx))
-     msh%elements(i)%e%pts(5)%p%x(1) = mean_data%x%x(linear_index(1,1,lx,i,lx,lx,lx))
-     msh%elements(i)%e%pts(6)%p%x(1) = mean_data%x%x(linear_index(lx,1,lx,i,lx,lx,lx))
-     msh%elements(i)%e%pts(7)%p%x(1) = mean_data%x%x(linear_index(1,lx,lx,i,lx,lx,lx))
-     msh%elements(i)%e%pts(8)%p%x(1) = mean_data%x%x(linear_index(lx,lx,lx,i,lx,lx,lx))
-
-     msh%elements(i)%e%pts(1)%p%x(2) = mean_data%y%x(linear_index(1,1,1,i,lx,lx,lx))
-     msh%elements(i)%e%pts(2)%p%x(2) = mean_data%y%x(linear_index(lx,1,1,i,lx,lx,lx))
-     msh%elements(i)%e%pts(3)%p%x(2) = mean_data%y%x(linear_index(1,lx,1,i,lx,lx,lx))
-     msh%elements(i)%e%pts(4)%p%x(2) = mean_data%y%x(linear_index(lx,lx,1,i,lx,lx,lx))
-     msh%elements(i)%e%pts(5)%p%x(2) = mean_data%y%x(linear_index(1,1,lx,i,lx,lx,lx))
-     msh%elements(i)%e%pts(6)%p%x(2) = mean_data%y%x(linear_index(lx,1,lx,i,lx,lx,lx))
-     msh%elements(i)%e%pts(7)%p%x(2) = mean_data%y%x(linear_index(1,lx,lx,i,lx,lx,lx))
-     msh%elements(i)%e%pts(8)%p%x(2) = mean_data%y%x(linear_index(lx,lx,lx,i,lx,lx,lx))
-
-     msh%elements(i)%e%pts(1)%p%x(3) = mean_data%z%x(linear_index(1,1,1,i,lx,lx,lx))
-     msh%elements(i)%e%pts(2)%p%x(3) = mean_data%z%x(linear_index(lx,1,1,i,lx,lx,lx))
-     msh%elements(i)%e%pts(3)%p%x(3) = mean_data%z%x(linear_index(1,lx,1,i,lx,lx,lx))
-     msh%elements(i)%e%pts(4)%p%x(3) = mean_data%z%x(linear_index(lx,lx,1,i,lx,lx,lx))
-     msh%elements(i)%e%pts(5)%p%x(3) = mean_data%z%x(linear_index(1,1,lx,i,lx,lx,lx))
-     msh%elements(i)%e%pts(6)%p%x(3) = mean_data%z%x(linear_index(lx,1,lx,i,lx,lx,lx))
-     msh%elements(i)%e%pts(7)%p%x(3) = mean_data%z%x(linear_index(1,lx,lx,i,lx,lx,lx))
-     msh%elements(i)%e%pts(8)%p%x(3) = mean_data%z%x(linear_index(lx,lx,lx,i,lx,lx,lx))
+     lx = stats_data%lx
+     msh%elements(i)%e%pts(1)%p%x(1) = stats_data%x%x(linear_index(1,1,1,i,lx,lx,lx))
+     msh%elements(i)%e%pts(2)%p%x(1) = stats_data%x%x(linear_index(lx,1,1,i,lx,lx,lx))
+     msh%elements(i)%e%pts(3)%p%x(1) = stats_data%x%x(linear_index(1,lx,1,i,lx,lx,lx))
+     msh%elements(i)%e%pts(4)%p%x(1) = stats_data%x%x(linear_index(lx,lx,1,i,lx,lx,lx))
+     msh%elements(i)%e%pts(5)%p%x(1) = stats_data%x%x(linear_index(1,1,lx,i,lx,lx,lx))
+     msh%elements(i)%e%pts(6)%p%x(1) = stats_data%x%x(linear_index(lx,1,lx,i,lx,lx,lx))
+     msh%elements(i)%e%pts(7)%p%x(1) = stats_data%x%x(linear_index(1,lx,lx,i,lx,lx,lx))
+     msh%elements(i)%e%pts(8)%p%x(1) = stats_data%x%x(linear_index(lx,lx,lx,i,lx,lx,lx))
+
+     msh%elements(i)%e%pts(1)%p%x(2) = stats_data%y%x(linear_index(1,1,1,i,lx,lx,lx))
+     msh%elements(i)%e%pts(2)%p%x(2) = stats_data%y%x(linear_index(lx,1,1,i,lx,lx,lx))
+     msh%elements(i)%e%pts(3)%p%x(2) = stats_data%y%x(linear_index(1,lx,1,i,lx,lx,lx))
+     msh%elements(i)%e%pts(4)%p%x(2) = stats_data%y%x(linear_index(lx,lx,1,i,lx,lx,lx))
+     msh%elements(i)%e%pts(5)%p%x(2) = stats_data%y%x(linear_index(1,1,lx,i,lx,lx,lx))
+     msh%elements(i)%e%pts(6)%p%x(2) = stats_data%y%x(linear_index(lx,1,lx,i,lx,lx,lx))
+     msh%elements(i)%e%pts(7)%p%x(2) = stats_data%y%x(linear_index(1,lx,lx,i,lx,lx,lx))
+     msh%elements(i)%e%pts(8)%p%x(2) = stats_data%y%x(linear_index(lx,lx,lx,i,lx,lx,lx))
+
+     msh%elements(i)%e%pts(1)%p%x(3) = stats_data%z%x(linear_index(1,1,1,i,lx,lx,lx))
+     msh%elements(i)%e%pts(2)%p%x(3) = stats_data%z%x(linear_index(lx,1,1,i,lx,lx,lx))
+     msh%elements(i)%e%pts(3)%p%x(3) = stats_data%z%x(linear_index(1,lx,1,i,lx,lx,lx))
+     msh%elements(i)%e%pts(4)%p%x(3) = stats_data%z%x(linear_index(lx,lx,1,i,lx,lx,lx))
+     msh%elements(i)%e%pts(5)%p%x(3) = stats_data%z%x(linear_index(1,1,lx,i,lx,lx,lx))
+     msh%elements(i)%e%pts(6)%p%x(3) = stats_data%z%x(linear_index(lx,1,lx,i,lx,lx,lx))
+     msh%elements(i)%e%pts(7)%p%x(3) = stats_data%z%x(linear_index(1,lx,lx,i,lx,lx,lx))
+     msh%elements(i)%e%pts(8)%p%x(3) = stats_data%z%x(linear_index(lx,lx,lx,i,lx,lx,lx))
   end do
 
-  call Xh%init(GLL, mean_data%lx, mean_data%ly, mean_data%lz)
+  call Xh%init(GLL, stats_data%lx, stats_data%ly, stats_data%lz)
 
   call dof%init(msh, Xh)
   call gs_h%init(dof)
@@ -116,13 +111,8 @@ program postprocess_fluid_stats
   w => neko_field_registry%get_field('w')
   p => neko_field_registry%get_field('p')
 
-  call avg_flow%init(u, v, w, p)
-  call fld_stats%init(coef,avg_flow%u,avg_flow%v,avg_flow%w,avg_flow%p)
-  n = mean_data%u%n
-  call copy(avg_flow%u%mf%x,mean_data%u%x,n)
-  call copy(avg_flow%v%mf%x,mean_data%v%x,n)
-  call copy(avg_flow%w%mf%x,mean_data%w%x,n)
-  call copy(avg_flow%p%mf%x,mean_data%p%x,n)
+  call fld_stats%init(coef,u,v,w,p)
+  n = stats_data%u%n
 
   call copy(fld_stats%stat_fields%items(1)%ptr%x, stats_data%p%x,n)
   call copy(fld_stats%stat_fields%items(2)%ptr%x, stats_data%u%x,n)
diff --git a/doc/pages/user-guide/simcomps.md b/doc/pages/user-guide/simcomps.md
index 2193fd250af..ee31e81e0b0 100644
--- a/doc/pages/user-guide/simcomps.md
+++ b/doc/pages/user-guide/simcomps.md
@@ -31,6 +31,7 @@ in Neko. The list will be updated as new simcomps are added.
 - Computation of forces and torque on a surface \ref simcomp_force_torque
 - Computation of the weak gradient of a field \ref simcomp_weak_grad
 - User defined components \ref user-file_simcomps
+- Fluid statistics simcomp, "fluid_stats", for more details see the [statistics guide](@ref statistics-guide)
 
 ## Controling execution and file output
 Each simulation component is, by default, executed once per time step to perform
diff --git a/doc/pages/user-guide/statistics-guide.md b/doc/pages/user-guide/statistics-guide.md
index b6ce8f92e89..288144455b0 100644
--- a/doc/pages/user-guide/statistics-guide.md
+++ b/doc/pages/user-guide/statistics-guide.md
@@ -2,8 +2,6 @@
 
 \tableofcontents
 
-Under development, updated incrementally
-
 Statistics in the context of Neko, is the common name for fields that are averaged in time and possible also space.
 
 The statistics module in Neko computes the temporal average of a wide range of fields.
@@ -29,84 +27,91 @@ where \f$ u_0 \f$ is the fields value at \f$ T_0 \f$ and \f$ N \f$ is the number
 In the statistics in Neko, various averages of the the different velocity components, derivatives and pressure are computed. In total, 44 "raw statistics" are computed that are required to compute the Reynolds stress budgets, mean fields, and the different terms in the turbulent kinetic energy equation.
 
 ## Using statistics
-Statistics are enable in the the case file as the following:
+Statistics are enable in the the case file as a simcomp with the added argument `avg_direction`, `set_of_stats`, and `start_time`:
 
 | Name                | Description                                                          | Admissible values | Default value |
 | ------------------- | -------------------------------------------------------------------- | ----------------- | ------------- |
-| `enabled`           | Whether to enable the statistics computation.                        | `true` or `false` | `true`        |
-| `start_time`        | Time at which to start gathering statistics.                         | Positive real     | 0             |
-| `sampling_interval` | Interval, in timesteps, for sampling the flow fields for statistics. | Positive integer  | 10            |
+| `start_time`        | Time at which to start gathering statistics.                        | Positive real     | 0             |
+| `avg_direction`        | Directions to compute spatial average.                         | x,y,z,xy,xz,yz  |  No spatial average           |
+| `set_of_stats`        | What set of stats to compute.                         | basic, full  |  full         |
+| `compute_value` | Interval, in timesteps or simulationtime, depending on compute\_control, for sampling the flow fields for statistics. | Positive real or int  | Not set (but recommended with every 50 timesteps or so  |
 
 In addition to the usual controls for the output, which then outputs the averages computes from the last time the statistics were written to file.
 
-For example, if one wants to sample the fields every 4 time steps and compute the averages in time intervals of 20 and write the output every 20 time units, and start collecting statistics after an initial transient of 50 time units the following would work:
+For example, if one wants to compute only the basic statistics and sample the fields every 4 time steps and compute and output batches every 20 time units and have an initial transient of 60 time units the following would work:
 
 ~~~~~~~~~~~~~~~{.json}
-"statistics": {
-        "enabled": true,
-        "start_time": 50.0,
-        "sampling_interval": 4,
-        "output_control": "simulationtime",
-        "output_value": 20,
-  }
+"simulation_components": 
+  [
+    {
+      "type": "fluid_stats",
+      "compute_control": "tsteps",
+      "compute_value": 4,
+      "output_control": "simulationtime",
+      "output_value": 20,
+      "start_time":60.0,
+      "avg_direction":"xz",
+      "set_of_stats":"basic"
+    }1
+  ]
 ~~~~~~~~~~~~~~~
-When the output is written one obtains two .fld files called mean_field and stats. 
 
-## List of fields in output files 
-In `mean_field` the following averages are stored. The stored in variable column is which field one finds the computed statistic if one opens the file in paraview or visit.
+Preferably set the initial transient to a multiple of output_value as otherwise the first output will be slightly shorter than the rest. The code related to fluid statistics are located in fluid_stats and fluid_stats_simcomp.
 
-| Number | Statistic | Stored in variable |
-| ------ | --------- | ------------------ |
-| 1 | \f$ \langle p \rangle \f$ | Pressure|
-| 2 | \f$ \langle u \rangle \f$ | X-Velocity|
-| 3 | \f$ \langle v \rangle \f$ | Y-Velocity|
-| 4 | \f$ \langle w \rangle \f$ | Z-Velocity|
+The argument "avg_direction" is optional and if ignored we output 3d fields. The statistics are saved in a fld file according to the following in 2D and 3D. Observe that in 2D the mean Z-velocity is stored in a last scalar field. All other fields are kept the same. This is due to a limitation of the fld file format.
+
+For 1D statistics a CSV file is outputted. The first column is the time at which the statistics are collected, the second column the spatial coordinate, and the rest of the data is stored in the order below. In this case all statistics are kept in the same order as in 3D.
 
-In `stats` several other statistics are stored, and while not all might be interesting to your specific use case, with them most different budgets and quantities of interest can be computed. They are stored as the following:
+## List of fields in output files
 
+When only the basic set of stats is enabled, only stats 1-11 are computed. When 2D stats are enabled \f$ \langle w \rangle \f$ is stored in s7 for basic stats and in s40 for the full set of statistics.
 
-| Number | Statistic | Stored in variable |
+| Number | Statistic | Stored in variable (for fld files) |
 | ------ | --------- | ------------------ |
-| 1 | \f$ \langle pp \rangle \f$ | Pressure|
-| 2 | \f$ \langle uu \rangle \f$ | X-Velocity|
-| 3 | \f$ \langle vv \rangle \f$ | Y-Velocity|
-| 4 | \f$ \langle ww \rangle \f$ | Z-Velocity|
-| 5 | \f$ \langle uv \rangle \f$ | Temperature|
-| 6 | \f$ \langle uw \rangle \f$ | Scalar 1 (s1)|
-| 7 | \f$ \langle vw \rangle \f$ | Scalar 2 (s2)|
-| 8 | \f$ \langle uuu \rangle \f$ | s3|
-| 9 | \f$ \langle vvv \rangle \f$ | s4|
-| 10 | \f$ \langle www \rangle \f$ | s5|
-| 11 | \f$ \langle  uuv   \rangle \f$ | s6 |
-| 12 | \f$ \langle  uuw   \rangle \f$ | s7 |
-| 13 | \f$ \langle  uvv   \rangle \f$ | s8 |
-| 14 | \f$ \langle  uvw   \rangle \f$ | s9 |
-| 15 | \f$ \langle  vvw   \rangle \f$ | s10 |
-| 16 | \f$ \langle  uww   \rangle \f$ | s11 |
-| 17 | \f$ \langle  vww   \rangle \f$ | s12 |
-| 18 | \f$ \langle  uuuu  \rangle \f$ | s13 |
-| 19 | \f$ \langle  vvvv  \rangle \f$ | s14 |
-| 20 | \f$ \langle wwww   \rangle \f$ | s15 |
-| 21 | \f$ \langle  ppp   \rangle \f$ | s16 |
-| 22 | \f$ \langle  pppp  \rangle \f$ | s17 |
-| 23 | \f$ \langle  pu    \rangle \f$ | s18 |
-| 24 | \f$ \langle  pv    \rangle \f$ | s19 |
-| 25 | \f$ \langle  pw    \rangle \f$ | s20 |
-| 26 | \f$ \langle  p \frac{\partial u} {\partial x} \rangle \f$ | s21 |
-| 27 | \f$ \langle  p \frac{\partial u} {\partial y}\rangle \f$ | s22 |
-| 28 | \f$ \langle  p \frac{\partial u} {\partial z}\rangle \f$ | s23 |
-| 29 | \f$ \langle  p \frac{\partial v} {\partial x}\rangle \f$ | s24 |
-| 30 | \f$ \langle  p \frac{\partial v} {\partial y}\rangle \f$ | s25 |
-| 31 | \f$ \langle  p \frac{\partial v} {\partial z}\rangle \f$ | s26 |
-| 32 | \f$ \langle  p \frac{\partial w} {\partial x}\rangle \f$ | s27 |
-| 33 | \f$ \langle  p \frac{\partial w} {\partial y}\rangle \f$ | s28 |
-| 34 | \f$ \langle  p \frac{\partial w} {\partial z}\rangle \f$ | s29 |
-| 35 | \f$ \langle  e11   \rangle \f$ | s30 |
-| 36 | \f$ \langle  e22   \rangle \f$ | s31 |
-| 37 | \f$ \langle  e33   \rangle \f$ | s32 |
-| 38 | \f$ \langle  e12   \rangle \f$ | s33 |
-| 39 | \f$ \langle  e13   \rangle \f$ | s34 |
-| 40 | \f$ \langle  e23   \rangle \f$ | s35 |
+| 1 | \f$ \langle p \rangle \f$ | Pressure|
+| 2 | \f$ \langle u \rangle \f$ | X-Velocity|
+| 3 | \f$ \langle v \rangle \f$ | Y-Velocity|
+| 4 | \f$ \langle w \rangle \f$ | Z-Velocity |
+| 5 | \f$ \langle pp \rangle \f$ | Temperature|
+| 6 | \f$ \langle uu \rangle \f$ | Scalar 1 (s1)|
+| 7 | \f$ \langle vv \rangle \f$ | Scalar 2 (s2)|
+| 8 | \f$ \langle ww \rangle \f$ | s3|
+| 9 | \f$ \langle uv \rangle \f$ | s4|
+|10 | \f$ \langle uw \rangle \f$ | s5|
+| 11 | \f$ \langle vw \rangle \f$ | s6|
+| 12| \f$ \langle uuu \rangle \f$ | s7|
+| 13| \f$ \langle vvv \rangle \f$ | s8|
+| 14 | \f$ \langle www \rangle \f$ | s9|
+| 15 | \f$ \langle  uuv   \rangle \f$ | s10 |
+| 16 | \f$ \langle  uuw   \rangle \f$ | s11 |
+| 17 | \f$ \langle  uvv   \rangle \f$ | s12 |
+| 18 | \f$ \langle  uvw   \rangle \f$ | s13 |
+| 19 | \f$ \langle  vvw   \rangle \f$ | s14 |
+| 20 | \f$ \langle  uww   \rangle \f$ | s15 |
+| 21 | \f$ \langle  vww   \rangle \f$ | s16 |
+| 22 | \f$ \langle  uuuu  \rangle \f$ | s17 |
+| 23 | \f$ \langle  vvvv  \rangle \f$ | s18 |
+| 24 | \f$ \langle wwww   \rangle \f$ | s19 |
+| 25 | \f$ \langle  ppp   \rangle \f$ | s20 |
+| 26 | \f$ \langle  pppp  \rangle \f$ | s21 |
+| 27 | \f$ \langle  pu    \rangle \f$ | s22 |
+| 28 | \f$ \langle  pv    \rangle \f$ | s23 |
+| 29 | \f$ \langle  pw    \rangle \f$ | s24 |
+| 30 | \f$ \langle  p \frac{\partial u} {\partial x} \rangle \f$ | s25 |
+| 31 | \f$ \langle  p \frac{\partial u} {\partial y}\rangle \f$ | s26 |
+| 32 | \f$ \langle  p \frac{\partial u} {\partial z}\rangle \f$ | s27 |
+| 33 | \f$ \langle  p \frac{\partial v} {\partial x}\rangle \f$ | s28 |
+| 34 | \f$ \langle  p \frac{\partial v} {\partial y}\rangle \f$ | s29 |
+| 35 | \f$ \langle  p \frac{\partial v} {\partial z}\rangle \f$ | s30 |
+| 36 | \f$ \langle  p \frac{\partial w} {\partial x}\rangle \f$ | s31 |
+| 37 | \f$ \langle  p \frac{\partial w} {\partial y}\rangle \f$ | s32 |
+| 38 | \f$ \langle  p \frac{\partial w} {\partial z}\rangle \f$ | s33 |
+| 39 | \f$ \langle  e11   \rangle \f$ | s34 |
+| 40 | \f$ \langle  e22   \rangle \f$ | s35 |
+| 41 | \f$ \langle  e33   \rangle \f$ | s36 |
+| 42 | \f$ \langle  e12   \rangle \f$ | s37 |
+| 43 | \f$ \langle  e13   \rangle \f$ | s38 |
+| 44 | \f$ \langle  e23   \rangle \f$ | s39 |
 
 where \f$e11,e22...\f$ is computed as:
 $$
@@ -122,50 +127,10 @@ $$
 
 
 # Postprocessing
-Of course, these statistics are only the "raw statistics" in the sense that in general we are not interested in \f$ \langle uu\rangle \f$, but rather say the rms of the velocity fluctuation. FOr this we need to postprocess the statistics. 
-
-There is some rudimentary postprocessing to compute the spatial averages of fld filesa and also to combine the statistics collected from several runs (compute average in time) and also compute both the mean velocity gradient and the Reynolds stresses available among the contrib scripts. By running the contrib scripts without any arguments one gets a hint on their usage, and also the text below gives a guide on how to postprocess the raw statistics. The postprocessing part of Neko is expanding and changing quite a lot at the moment, where we currently envision primarily using python for the postprocessing of the final statistics.
-
-
-
-## Notes on the statistics calculation in Neko
-***Daniele Massaro, Martin Karp (KTH)***
-
-1) Run your simulations and collect mean_field* and stats* files by having the statistics object added to the case file,  and specifying the write interval to something suitable.
-
-2) For each RUN_i, you get a set of mean_field* and stats* files. You can average them for each single RUN_i, or average all of them only once (after re-ordering them properly). If you follow the second approach, go to step 4. 
-Here, for each RUN_i, we compute the averaged means with "average_fields_in_time":
---mean
-`srun --unbuffered /your/location/neko/bin/average_fields_in_time meanXX.fld T0 mean_p.fld`
-where T0 is the initial time. To get some hints on the input for the script one can simply run `./average_fields_in_time` without any arguments. For RUN_1 the time T0 can be taken from the log of the first simulation, or from the header of the first mean_field* file; in this way you discard that file. For RUN_i, with i>1, it can be taken from header of the last file mean_field* of the previous simulation RUN_{i-1}. 
-In the command line, for the name "meanXX.fld", XX indicates the number of the nek5000 file. In mean_fieldXX.nek5000 you set the number of the first mean0* file to read and the number of steps corresponding to the number of files. In this way, the code generates a mean_p0.f00000 and mean_post0.nek5000. It is suggested to rename mean_p0.f00000 as mean_p0.f0000i and move it to a separate folder where you take the average with all the others. 
---stats
-`srun --unbuffered /your/location/neko/bin/average_fields_in_time statXX.fld T0 stat_p.fld`
-T0 is the same as before. In stat0.nek5000 you set the number of the first stat0* file to read and the number of steps corresponds to the number of files. It is suggested to rename stat_p0.f00000 as stat_p0.f0000i and move it to a separate folder where you take the average with all the others. 
-Repeat this for each RUN_i folder. Eventually, given n RUN_i folders, you will get n mean_p* and stat_p* files.
-
-3) Take the average of the averaged runs. Now, the time average over all the n simulations is taken. The procedure is similar, but changing the output name is recommended to  avoid over-writing.
--- mean
-`srun --unbuffered /your/location/neko/bin/average_fields mean_p0.fld T0 mean_post.fld`
-where T0 is the initial time which has been used to compute mean_p* for RUN_1.
--- stats
-`srun --unbuffered /your/location/neko/bin/average_fields stat_p0.fld T0 stat_post.fld`
-where T0 is the initial time which has been used to compute mean_p* for RUN_1.
-
-
-
-
-
-4) Compute Reynolds stress tensors and other statistical moments (see the list).
-`srun --unbuffered /your/location/neko/bin/postprocess_fluid_stats mesh_file.nmsh mean_post0.fld stat_post0.fld`
-
-5) We also provide a tool to average the resulting field in a homogenous direction in `bin/average_field_in_space`. The required arguments are shown if one runs the program without any input. Currently it requires the number of elements in the homogenous direction as an input argument, e.g. 
-`./average_field_in_space mesh.nmsh field.fld x 18 outfield.fld`
-if we want to average a field in the x direction on a mesh with 18 elements in x and output the averaged field in outfield0.nek5000.
-
-
-
+These statistics are only the "raw statistics" in the sense that in general we are not interested in \f$ \langle uu\rangle \f$, but rather say \f$ \langle u'u'\rangle\f$. For this we need to postprocess the statistics. 
 
+There is some rudimentary postprocessing to compute the spatial averages of fld files and also to combine the statistics collected from several runs (compute average in time) and also compute both the mean velocity gradient and the Reynolds stresses available in the contrib scripts. By running the contrib scripts without any arguments one gets a hint on their usage, (e.g. by running `./average_fields_in_time` will give the options). 
 
+To further postprocess the statistics it is suggested to look into PyNekTools which introduces convenient functions for postprocessing, largely based on the PyMech library with the addition of an expanding set of tools for interpolation, computing derivatives and more advanced functionality.  PyNekTools is entirely parallelized in MPI and can also handle large data sets for postprocessing.
 
 
diff --git a/examples/turb_channel/turb_channel.case b/examples/turb_channel/turb_channel.case
index 9cc66018482..b94df7ee963 100644
--- a/examples/turb_channel/turb_channel.case
+++ b/examples/turb_channel/turb_channel.case
@@ -25,14 +25,14 @@
         "velocity_solver": {
             "type": "cg",
             "preconditioner": "jacobi",
-            "projection_space_size": 5,
+            "projection_space_size": 0,
             "absolute_tolerance": 1e-4,
             "max_iterations": 800
         },
         "pressure_solver": {
             "type": "gmres",
             "preconditioner": "hsmg",
-            "projection_space_size": 20,
+            "projection_space_size": 5,
             "absolute_tolerance": 1e-3,
             "max_iterations": 800
         },
@@ -43,26 +43,29 @@
         },
         "boundary_types": ["", "", "w", "w"],
         "output_control": "simulationtime",
-        "output_value": 1.0
-    },
-    "statistics": {
-        "enabled": true,
-        "output_control": "simulationtime",
-        "output_value": 20,
-        "start_time": 70.0,
-        "sampling_interval": 10
+        "output_value": 2.0
     },
  "simulation_components": 
   [
     {
       "type": "force_torque",
       "compute_control": "tsteps",
-      "compute_value": 5,
+      "compute_value": 10,
       "zone_id": 3,
-      "scale": 0.01899,
-      "center": [0,0,0]
+      "zone_name": "wall",
+      "center": [0,0,0],
+      "scale": 0.01899
+    },    
+    { 
+      "type": "fluid_stats",
+      "output_control": "simulationtime",
+      "output_value": 20,
+      "start_time": 60.0,
+      "compute_control": "tsteps",
+      "compute_value": 50,
+      "avg_direction": "xz",
+      "set_of_stats": "full"
     }
   ]   
-
 }
 }
diff --git a/src/.depends b/src/.depends
index 5929f430752..d7c005b40fa 100644
--- a/src/.depends
+++ b/src/.depends
@@ -126,7 +126,7 @@ io/fluid_output.o : io/fluid_output.f90 io/output.o device/device.o config/neko_
 io/fld_file_output.o : io/fld_file_output.f90 io/output.o device/device.o config/neko_config.o field/field_list.o config/num_types.o 
 io/chkp_output.o : io/chkp_output.f90 config/num_types.o io/output.o common/checkpoint.o 
 io/mean_flow_output.o : io/mean_flow_output.f90 io/output.o device/device.o config/num_types.o fluid/mean_flow.o 
-io/fluid_stats_output.o : io/fluid_stats_output.f90 io/output.o device/device.o config/num_types.o config/neko_config.o fluid/fluid_stats.o 
+io/fluid_stats_output.o : io/fluid_stats_output.f90 math/matrix.o io/output.o device/device.o io/fld_file_data.o sem/map_2d.o sem/map_1d.o config/num_types.o config/neko_config.o fluid/fluid_stats.o 
 io/mean_sqr_flow_output.o : io/mean_sqr_flow_output.f90 io/output.o config/num_types.o fluid/mean_sqr_flow.o 
 io/data_streamer.o : io/data_streamer.F90 config/neko_config.o comm/mpi_types.o comm/comm.o device/device.o common/utils.o sem/coef.o field/field.o config/num_types.o 
 common/sampler.o : common/sampler.f90 common/time_based_controller.o config/num_types.o common/profiler.o common/utils.o common/log.o comm/comm.o io/fld_file.o io/output.o 
@@ -186,12 +186,13 @@ common/rhs_maker.o : common/rhs_maker.f90 field/field.o field/field_series.o con
 common/rhs_maker_fctry.o : common/rhs_maker_fctry.f90 config/neko_config.o common/bcknd/device/rhs_maker_device.o common/bcknd/sx/rhs_maker_sx.o common/bcknd/cpu/rhs_maker_cpu.o common/rhs_maker.o 
 simulation_components/probes.o : simulation_components/probes.F90 case.o io/csv_file.o io/file.o device/device.o comm/comm.o mesh/point_zone_registry.o mesh/point_zone.o math/tensor.o common/global_interpolation.o common/json_utils.o sem/dofmap.o field/field_registry.o simulation_components/simulation_component.o field/field_list.o common/utils.o common/log.o math/matrix.o config/num_types.o 
 simulation_components/force_torque.o : simulation_components/force_torque.f90 device/device.o math/bcknd/device/device_math.o math/math.o comm/comm.o qoi/drag_torque.o bc/dirichlet.o math/vector.o sem/coef.o simulation_components/field_writer.o common/json_utils.o io/fld_file_output.o case.o math/operators.o field/field.o field/field_registry.o simulation_components/simulation_component.o config/num_types.o 
+simulation_components/fluid_stats_simcomp.o : simulation_components/fluid_stats_simcomp.f90 common/json_utils.o common/log.o comm/comm.o sem/coef.o case.o io/fluid_stats_output.o fluid/fluid_stats.o field/field.o field/field_registry.o simulation_components/simulation_component.o config/num_types.o 
 simulation_components/field_writer.o : simulation_components/field_writer.f90 common/json_utils.o io/fld_file_output.o case.o math/operators.o field/field.o field/field_registry.o simulation_components/simulation_component.o config/num_types.o 
 common/bcknd/cpu/rhs_maker_cpu.o : common/bcknd/cpu/rhs_maker_cpu.f90 field/scratch_registry.o config/num_types.o field/field.o field/field_series.o common/rhs_maker.o 
 common/bcknd/sx/rhs_maker_sx.o : common/bcknd/sx/rhs_maker_sx.f90 field/scratch_registry.o config/num_types.o field/field.o field/field_series.o common/rhs_maker.o 
 common/bcknd/device/rhs_maker_device.o : common/bcknd/device/rhs_maker_device.F90 config/num_types.o field/field.o field/field_series.o common/utils.o device/device.o common/rhs_maker.o 
 config/neko_config.o : config/neko_config.f90 
-case.o : case.f90 mesh/point_zone_registry.o field/scratch_registry.o common/json_utils.o scalar/scalar_pnpn.o common/user_intf.o common/jobctrl.o common/log.o time_schemes/time_scheme_controller.o comm/comm.o mesh/mesh.o common/utils.o io/file.o common/statistics.o scalar/scalar_ic.o fluid/flow_ic.o common/sampler.o comm/redist.o comm/parmetis.o field/mesh_field.o io/fluid_stats_output.o io/mean_flow_output.o io/mean_sqr_flow_output.o io/chkp_output.o io/fluid_output.o fluid/fluid_scheme.o fluid/fluid_pnpn.o config/num_types.o 
+case.o : case.f90 mesh/point_zone_registry.o field/scratch_registry.o common/json_utils.o scalar/scalar_pnpn.o common/user_intf.o common/jobctrl.o common/log.o time_schemes/time_scheme_controller.o comm/comm.o mesh/mesh.o common/utils.o io/file.o common/statistics.o scalar/scalar_ic.o fluid/flow_ic.o common/sampler.o comm/redist.o comm/parmetis.o field/mesh_field.o io/chkp_output.o io/fluid_output.o fluid/fluid_scheme.o fluid/fluid_pnpn.o config/num_types.o 
 common/user_intf.o : common/user_intf.f90 common/log.o common/utils.o common/json_utils.o config/num_types.o bc/field_dirichlet.o bc/usr_scalar.o bc/usr_inflow.o mesh/mesh.o bc/bc.o sem/coef.o scalar/scalar_user_source_term.o fluid/fluid_user_source_term.o field/field_list.o field/field.o 
 fluid/stress_formulation/pnpn_res_stress_fctry.o : fluid/stress_formulation/pnpn_res_stress_fctry.f90 fluid/stress_formulation/bcknd/device/pnpn_res_stress_device.o fluid/stress_formulation/bcknd/cpu/pnpn_res_stress_cpu.o config/neko_config.o fluid/pnpn_res.o 
 fluid/stress_formulation/bcknd/cpu/pnpn_res_stress_cpu.o : fluid/stress_formulation/bcknd/cpu/pnpn_res_stress_cpu.f90 math/math.o sem/space.o config/num_types.o mesh/mesh.o field/scratch_registry.o fluid/pnpn_res.o bc/facet_normal.o sem/coef.o math/ax.o field/field.o math/operators.o gs/gather_scatter.o 
@@ -280,7 +281,7 @@ simulation_components/lambda2.o : simulation_components/lambda2.f90 device/devic
 simulation_components/weak_grad.o : simulation_components/weak_grad.f90 simulation_components/field_writer.o common/json_utils.o io/fld_file_output.o case.o math/operators.o field/field.o field/field_registry.o simulation_components/simulation_component.o config/num_types.o 
 simulation_components/derivative.o : simulation_components/derivative.f90 common/utils.o simulation_components/field_writer.o common/json_utils.o io/fld_file_output.o case.o math/operators.o field/field.o field/field_registry.o simulation_components/simulation_component.o config/num_types.o 
 simulation_components/les_simcomp.o : simulation_components/les_simcomp.f90 simulation_components/field_writer.o common/json_utils.o les/les_model.o case.o simulation_components/simulation_component.o config/num_types.o 
-simulation_components/simulation_component_fctry.o : simulation_components/simulation_component_fctry.f90 simulation_components/derivative.o simulation_components/weak_grad.o simulation_components/field_writer.o common/utils.o simulation_components/les_simcomp.o simulation_components/probes.o simulation_components/lambda2.o simulation_components/force_torque.o simulation_components/vorticity.o simulation_components/simulation_component.o 
+simulation_components/simulation_component_fctry.o : simulation_components/simulation_component_fctry.f90 simulation_components/derivative.o simulation_components/weak_grad.o simulation_components/field_writer.o common/utils.o simulation_components/les_simcomp.o simulation_components/probes.o simulation_components/lambda2.o simulation_components/fluid_stats_simcomp.o simulation_components/force_torque.o simulation_components/vorticity.o simulation_components/simulation_component.o 
 source_terms/source_term.o : source_terms/source_term.f90 field/field_list.o sem/coef.o config/num_types.o 
 source_terms/coriolis_source_term.o : source_terms/coriolis_source_term.f90 source_terms/bcknd/cpu/coriolis_source_term_cpu.o common/utils.o config/neko_config.o sem/coef.o source_terms/source_term.o common/json_utils.o field/field_list.o config/num_types.o 
 source_terms/bcknd/cpu/coriolis_source_term_cpu.o : source_terms/bcknd/cpu/coriolis_source_term_cpu.f90 field/field.o math/math.o field/field_list.o config/num_types.o 
diff --git a/src/Makefile.am b/src/Makefile.am
index a7445ac90af..440752f686b 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -189,6 +189,7 @@ neko_fortran_SOURCES = \
 	common/rhs_maker_fctry.f90\
 	simulation_components/probes.F90\
 	simulation_components/force_torque.f90\
+	simulation_components/fluid_stats_simcomp.f90\
 	simulation_components/field_writer.f90\
 	common/bcknd/cpu/rhs_maker_cpu.f90\
 	common/bcknd/sx/rhs_maker_sx.f90\
diff --git a/src/case.f90 b/src/case.f90
index 0e4e88f0e4e..3342bec7670 100644
--- a/src/case.f90
+++ b/src/case.f90
@@ -37,9 +37,6 @@ module case
   use fluid_scheme, only : fluid_scheme_t, fluid_scheme_factory
   use fluid_output, only : fluid_output_t
   use chkp_output, only : chkp_output_t
-  use mean_sqr_flow_output, only : mean_sqr_flow_output_t
-  use mean_flow_output, only : mean_flow_output_t
-  use fluid_stats_output, only : fluid_stats_output_t
   use mesh_field, only : mesh_fld_t, mesh_field_init, mesh_field_free
   use parmetis, only : parmetis_partmeshkway
   use redist, only : redist_mesh
@@ -71,13 +68,10 @@ module case
      real(kind=rp), dimension(10) :: dtlag
      real(kind=rp) :: dt
      real(kind=rp) :: end_time
+     character(len=:), allocatable :: output_directory
      type(sampler_t) :: s
      type(fluid_output_t) :: f_out
-     type(fluid_stats_output_t) :: f_stats_output
      type(chkp_output_t) :: f_chkp
-     type(mean_flow_output_t) :: f_mf
-     type(mean_sqr_flow_output_t) :: f_msqrf
-     type(stats_t) :: q
      type(user_t) :: usr
      class(fluid_scheme_t), allocatable :: fluid
      type(scalar_pnpn_t), allocatable :: scalar
@@ -136,7 +130,6 @@ end subroutine case_init_from_json
   !> Initialize a case from its (loaded) params object
   subroutine case_init_common(C)
     type(case_t), target, intent(inout) :: C
-    character(len=:), allocatable :: output_directory
     integer :: lx = 0
     logical :: scalar = .false.
     type(file_t) :: msh_file, bdry_file, part_file
@@ -318,14 +311,14 @@ subroutine case_init_common(C)
     ! Get and process output directory
     !
     call json_get_or_default(C%params, 'case.output_directory',&
-                             output_directory, '')
+                             C%output_directory, '')
 
-    output_dir_len = len(trim(output_directory))
+    output_dir_len = len(trim(C%output_directory))
     if (output_dir_len .gt. 0) then
-       if (output_directory(output_dir_len:output_dir_len) .ne. "/") then
-          output_directory = trim(output_directory)//"/"
+       if (C%output_directory(output_dir_len:output_dir_len) .ne. "/") then
+          C%output_directory = trim(C%output_directory)//"/"
           if (pe_rank .eq. 0) then
-             call execute_command_line('mkdir -p '//output_directory)
+             call execute_command_line('mkdir -p '//C%output_directory)
           end if
        end if
     end if
@@ -336,7 +329,7 @@ subroutine case_init_common(C)
     call json_get_or_default(C%params, 'case.output_boundary',&
                              logical_val, .false.)
     if (logical_val) then
-       bdry_file = file_t(trim(output_directory)//'bdry.fld')
+       bdry_file = file_t(trim(C%output_directory)//'bdry.fld')
        call bdry_file%write(C%fluid%bdry)
     end if
 
@@ -348,7 +341,7 @@ subroutine case_init_common(C)
     if (logical_val) then
        call mesh_field_init(msh_part, C%msh, 'MPI_Rank')
        msh_part%data = pe_rank
-       part_file = file_t(trim(output_directory)//'partitions.vtk')
+       part_file = file_t(trim(C%output_directory)//'partitions.vtk')
        call part_file%write(msh_part)
        call mesh_field_free(msh_part)
     end if
@@ -371,10 +364,10 @@ subroutine case_init_common(C)
     call C%s%init(C%end_time)
     if (scalar) then
        C%f_out = fluid_output_t(precision, C%fluid, C%scalar, &
-            path = trim(output_directory))
+            path = trim(C%output_directory))
     else
        C%f_out = fluid_output_t(precision, C%fluid, &
-            path = trim(output_directory))
+            path = trim(C%output_directory))
     end if
 
     call json_get_or_default(C%params, 'case.fluid.output_control',&
@@ -402,7 +395,7 @@ subroutine case_init_common(C)
     if (logical_val) then
        call json_get_or_default(C%params, 'case.checkpoint_format', &
             string_val, "chkp")
-       C%f_chkp = chkp_output_t(C%fluid%chkp, path = output_directory, &
+       C%f_chkp = chkp_output_t(C%fluid%chkp, path = C%output_directory, &
             fmt = trim(string_val))
        call json_get_or_default(C%params, 'case.checkpoint_control', &
             string_val, "simulationtime")
@@ -411,50 +404,6 @@ subroutine case_init_common(C)
        call C%s%add(C%f_chkp, real_val, string_val)
     end if
 
-    !
-    ! Setup statistics
-    !
-
-    ! Always init, so that we can call eval in simulation.f90 with no if.
-    ! Note, don't use json_get_or_default here, because that will break the
-    ! valid_path if statement below (the path will become valid always).
-    call C%params%get('case.statistics.start_time', stats_start_time,&
-                           found)
-    if (.not. found) stats_start_time = 0.0_rp
-
-    call C%params%get('case.statistics.sampling_interval', &
-                           stats_sampling_interval, found)
-    if (.not. found) stats_sampling_interval = 10
-
-    call C%q%init(stats_start_time, stats_sampling_interval)
-
-    found = C%params%valid_path('case.statistics')
-    if (found) then
-       call json_get_or_default(C%params, 'case.statistics.enabled',&
-                                logical_val, .true.)
-       if (logical_val) then
-          call C%q%add(C%fluid%mean%u)
-          call C%q%add(C%fluid%mean%v)
-          call C%q%add(C%fluid%mean%w)
-          call C%q%add(C%fluid%mean%p)
-
-          C%f_mf = mean_flow_output_t(C%fluid%mean, stats_start_time, &
-                                      path = output_directory)
-
-          call json_get(C%params, 'case.statistics.output_control', &
-                        string_val)
-          call json_get(C%params, 'case.statistics.output_value', &
-                        stats_output_val)
-
-          call C%s%add(C%f_mf, stats_output_val, string_val)
-          call C%q%add(C%fluid%stats)
-
-          C%f_stats_output = fluid_stats_output_t(C%fluid%stats, &
-            stats_start_time, path = output_directory)
-          call C%s%add(C%f_stats_output, stats_output_val, string_val)
-       end if
-    end if
-
     !
     ! Setup joblimit
     !
@@ -485,8 +434,6 @@ subroutine case_free(C)
 
     call C%s%free()
 
-    call C%q%free()
-
   end subroutine case_free
 
 end module case
diff --git a/src/field/mean_field.f90 b/src/field/mean_field.f90
index 50240a2a449..3c00e05ab0a 100644
--- a/src/field/mean_field.f90
+++ b/src/field/mean_field.f90
@@ -30,7 +30,7 @@
 ! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 ! POSSIBILITY OF SUCH DAMAGE.
 !
-!> Defines a mean field
+!> Implements mean_field_t.
 !
 module mean_field
   use neko_config, only : NEKO_BCKND_DEVICE
@@ -41,42 +41,53 @@ module mean_field
   implicit none
   private
 
+  !> Computes the temporal mean of a field.
   type, public, extends(stats_quant_t) :: mean_field_t
+     !> Pointer to the averaged field.
      type(field_t), pointer :: f => null()
+     !> Stores the mean field.
      type(field_t) :: mf
+     !> Total time across which the mean has been computed.
      real(kind=rp) :: time
    contains
+     !> Constructor.
      procedure, pass(this) :: init => mean_field_init
+     !> Destructor.
      procedure, pass(this) :: free => mean_field_free
+     !> Updates the mean value with a new sample.
      procedure, pass(this) :: update => mean_field_update
+     !> Resets the mean field.
      procedure, pass(this) :: reset => mean_field_reset
   end type mean_field_t
 
 contains
 
-  !> Initialize a mean field for a field @a f
+  !> Constructor.
+  !! @param f The field that will be averaged.
+  !! @param field_name. Optional name for the mean field. By default the name of
+  !! `f` prepended with `mean_` is used.
   subroutine mean_field_init(this, f, field_name)
     class(mean_field_t), intent(inout) :: this
     type(field_t), intent(inout), target :: f
     character(len=*), optional, intent(in) :: field_name
     character(len=80) :: name
 
-
     call this%free()
 
     this%f => f
     this%time = 0.0_rp
+
     if (present(field_name)) then
        name = field_name
     else
-       write(name, '(A,A)') 'mean_',trim(f%name)
+       write(name, '(A,A)') 'mean_', trim(f%name)
     end if
 
     call this%mf%init(f%dof, name)
 
   end subroutine mean_field_init
 
-  !> Deallocates a mean field
+  !> Destructor.
   subroutine mean_field_free(this)
     class(mean_field_t), intent(inout) :: this
 
@@ -87,7 +98,7 @@ subroutine mean_field_free(this)
 
   end subroutine mean_field_free
 
-  !> Resets a mean field
+  !> Resets a the mean field and the averaging time value to zero.
   subroutine mean_field_reset(this)
     class(mean_field_t), intent(inout) :: this
 
@@ -96,10 +107,11 @@ subroutine mean_field_reset(this)
   end subroutine mean_field_reset
 
 
-  !> Update a mean field
+  !> Update the mean field with a new sample.
+  !! @param k Time since last sample.
   subroutine mean_field_update(this, k)
     class(mean_field_t), intent(inout) :: this
-    real(kind=rp), intent(in) :: k !< Time since last sample
+    real(kind=rp), intent(in) :: k
 
     call field_cmult(this%mf, this%time, size(this%mf%x))
     call field_add2s2(this%mf, this%f, k, size(this%mf%x))
diff --git a/src/fluid/fluid_scheme.f90 b/src/fluid/fluid_scheme.f90
index a7d314434eb..94a1a3a2ee6 100644
--- a/src/fluid/fluid_scheme.f90
+++ b/src/fluid/fluid_scheme.f90
@@ -897,19 +897,6 @@ subroutine fluid_scheme_validate(this)
     !
     call this%chkp%init(this%u, this%v, this%w, this%p)
 
-    !
-    ! Setup mean flow fields if requested
-    !
-    if (this%params%valid_path('case.statistics')) then
-       call json_get_or_default(this%params, 'case.statistics.enabled', &
-                                logical_val, .true.)
-       if (logical_val) then
-          call this%mean%init(this%u, this%v, this%w, this%p)
-          call this%stats%init(this%c_Xh, this%mean%u, &
-               this%mean%v, this%mean%w, this%mean%p)
-       end if
-    end if
-
   end subroutine fluid_scheme_validate
 
   !> Apply all boundary conditions defined for velocity
diff --git a/src/fluid/fluid_stats.f90 b/src/fluid/fluid_stats.f90
index 32668b28400..55f77d614e2 100644
--- a/src/fluid/fluid_stats.f90
+++ b/src/fluid/fluid_stats.f90
@@ -1,4 +1,4 @@
-! Copyright (c) 2022-2023, The Neko Authors
+! Copyright (c) 2022-2024, The Neko Authors
 ! All rights reserved.
 !
 ! Redistribution and use in source and binary forms, with or without
@@ -51,20 +51,23 @@ module fluid_stats
   private
 
   type, public, extends(stats_quant_t) :: fluid_stats_t
-     type(field_t) :: stats_u !< Not reasonable to allocate 20 something fields
-     type(field_t) :: stats_v !< Not reasonable to allocate 20 something fields
-     type(field_t) :: stats_w !< Not reasonable to allocate 20 something fields
-     type(field_t) :: stats_p !< Not reasonable to allocate 20 something fields
-     type(field_t) :: stats_work !< Not reasonable to allocate 20 something fields
+     !> Work fields
+     type(field_t) :: stats_u
+     type(field_t) :: stats_v
+     type(field_t) :: stats_w
+     type(field_t) :: stats_p
+     type(field_t) :: stats_work
+
+     !> Pointers to the instantenious quantities.
      type(field_t), pointer :: u !< u
      type(field_t), pointer :: v !< v
      type(field_t), pointer :: w !< w
      type(field_t), pointer :: p !< p
 
-     type(field_t), pointer :: u_mean !< <u>
-     type(field_t), pointer :: v_mean !< <v>
-     type(field_t), pointer :: w_mean !< <w>
-     type(field_t), pointer :: p_mean !< <p>
+     type(mean_field_t) :: u_mean !< <u>
+     type(mean_field_t) :: v_mean !< <v>
+     type(mean_field_t) :: w_mean !< <w>
+     type(mean_field_t) :: p_mean !< <p>
      !> Velocity squares
      type(mean_field_t) :: uu !< <uu>
      type(mean_field_t) :: vv !< <vv>
@@ -125,276 +128,319 @@ module fluid_stats
      type(field_t) :: dwdy
      type(field_t) :: dwdz
 
+     !> SEM coefficients.
      type(coef_t), pointer :: coef
-     integer :: n_stats = 40
-     !> Used to write stats output
-     !! pressure=pp
-     !! x-vel=uu
-     !! y-vel=vv
-     !! z-vel=ww
-     !! temp=uv
-     !! scalar1=uw
-     !! scalar2=vw
-     !! The rest are stored in the above order in scalar fields
-     type(field_list_t)  :: stat_fields !< Used to write the output
+     !> Number of statistical fields to be computed.
+     integer :: n_stats = 44
+     !> Specifies a subset of the statistics to be collected. All 44 fields by
+     !! default.
+     character(5) :: stat_set
+     !> A list of size n_stats, whith entries pointing to the fields that will
+     !! be output (the field components above.) Used to write the output.
+     type(field_list_t)  :: stat_fields
    contains
+     !> Constructor.
      procedure, pass(this) :: init => fluid_stats_init
+     !> Destructor. 
      procedure, pass(this) :: free => fluid_stats_free
+     !> Update all the mean value fields with a new sample.
      procedure, pass(this) :: update => fluid_stats_update
+     !> Reset all the computed means values and sampling times to zero.
      procedure, pass(this) :: reset => fluid_stats_reset
+     ! Convert computed weak gradients to strong.
      procedure, pass(this) :: make_strong_grad => fluid_stats_make_strong_grad
+     !> Compute certain physical statistical quantities based on existing mean 
+     !! fields.
      procedure, pass(this) :: post_process => fluid_stats_post_process
   end type fluid_stats_t
 
 contains
 
-  !> Initialize the fields associated with fluid_stats
-  subroutine fluid_stats_init(this, coef, u_mf, v_mf, w_mf, p_mf)
-    class(fluid_stats_t), intent(inout), target :: this
+  !> Constructor. Initialize the fields associated with fluid_stats.
+  !! @param coef SEM coefficients. Optional.
+  !! @param u The x component of velocity.
+  !! @param v The y component of velocity.
+  !! @param w The z component of velocity.
+  !! @param p The pressure. 
+  !! @param set Specifies the subset of the statistics to be collected. 
+  !! Optional. Either `basic` or `full`, defaults to `full`.
+  subroutine fluid_stats_init(this, coef, u, v, w, p, set)
+    class(fluid_stats_t), intent(inout), target:: this
     type(coef_t), target, optional :: coef
-    type(mean_field_t), target, intent(inout) :: u_mf, v_mf, w_mf, p_mf
+    type(field_t), target, intent(inout) :: u, v, w, p
+    character(*), intent(in), optional :: set
+
+    call this%free()
     this%coef => coef
 
-    this%u_mean => u_mf%mf
-    this%v_mean => v_mf%mf
-    this%w_mean => w_mf%mf
-    this%p_mean => p_mf%mf
-    this%u => u_mf%f
-    this%v => v_mf%f
-    this%w => w_mf%f
-    this%p => p_mf%f
+    this%u => u
+    this%v => v
+    this%w => w
+    this%p => p
 
+    if (present(set)) then 
+       this%stat_set  = trim(set)
+       if (this%stat_set .eq. 'basic') then
+          this%n_stats = 11
+       end if
+    else
+       this%stat_set = 'full'
+       this%n_stats = 44
+    end if
+     
     call this%stats_work%init(this%u%dof, 'stats')
     call this%stats_u%init(this%u%dof, 'u temp')
     call this%stats_v%init(this%u%dof, 'v temp')
     call this%stats_w%init(this%u%dof, 'w temp')
     call this%stats_p%init(this%u%dof, 'p temp')
-
-    call this%dudx%init(this%u%dof, 'dudx')
-    call this%dudy%init(this%u%dof, 'dudy')
-    call this%dudz%init(this%u%dof, 'dudz')
-    call this%dvdx%init(this%u%dof, 'dvdx')
-    call this%dvdy%init(this%u%dof, 'dvdy')
-    call this%dvdz%init(this%u%dof, 'dvdz')
-    call this%dwdx%init(this%u%dof, 'dwdx')
-    call this%dwdy%init(this%u%dof, 'dwdy')
-    call this%dwdz%init(this%u%dof, 'dwdz')
-
+    call this%u_mean%init(this%u)
+    call this%v_mean%init(this%v)
+    call this%w_mean%init(this%w)
+    call this%p_mean%init(this%p)
     call this%uu%init(this%stats_u, 'uu')
     call this%vv%init(this%stats_v, 'vv')
     call this%ww%init(this%stats_w, 'ww')
     call this%uv%init(this%stats_work, 'uv')
     call this%uw%init(this%stats_work, 'uw')
     call this%vw%init(this%stats_work, 'vw')
-    call this%uuu%init(this%stats_work, 'uuu')!< <uuu>
-    call this%vvv%init(this%stats_work, 'vvv')!< <vvv>
-    call this%www%init(this%stats_work, 'www')!< <www>
-    call this%uuv%init(this%stats_work, 'uuv')!< <uuv>
-    call this%uuw%init(this%stats_work, 'uuw')!< <uuw>
-    call this%uvv%init(this%stats_work, 'uvv')!< <uvv>
-    call this%uvw%init(this%stats_work, 'uvw')!< <uvv>
-    call this%vvw%init(this%stats_work, 'vvw')!< <vvw>
-    call this%uww%init(this%stats_work, 'uww')!< <uww>
-    call this%vww%init(this%stats_work, 'vww')!< <vww>
-    call this%uuuu%init(this%stats_work, 'uuuu') !< <uuuu>
-    call this%vvvv%init(this%stats_work, 'vvvv') !< <vvvv>
-    call this%wwww%init(this%stats_work, 'wwww') !< <wwww>
-    !> Pressure
     call this%pp%init(this%stats_p, 'pp')
-    call this%ppp%init(this%stats_work, 'ppp')
-    call this%pppp%init(this%stats_work, 'pppp')
-    !> Pressure * velocity
-    call this%pu%init(this%stats_work, 'pu')
-    call this%pv%init(this%stats_work, 'pv')
-    call this%pw%init(this%stats_work, 'pw')
-
-    call this%pdudx%init(this%stats_work, 'pdudx')
-    call this%pdudy%init(this%stats_work, 'pdudy')
-    call this%pdudz%init(this%stats_work, 'pdudz')
-    call this%pdvdx%init(this%stats_work, 'pdvdx')
-    call this%pdvdy%init(this%stats_work, 'pdvdy')
-    call this%pdvdz%init(this%stats_work, 'pdvdz')
-    call this%pdwdx%init(this%stats_work, 'pdwdx')
-    call this%pdwdy%init(this%stats_work, 'pdwdy')
-    call this%pdwdz%init(this%stats_work, 'pdwdz')
-
-    call this%e11%init(this%stats_work, 'e11')
-    call this%e22%init(this%stats_work, 'e22')
-    call this%e33%init(this%stats_work, 'e33')
-    call this%e12%init(this%stats_work, 'e12')
-    call this%e13%init(this%stats_work, 'e13')
-    call this%e23%init(this%stats_work, 'e23')
+ 
+    if (this%n_stats .eq. 44) then
+       call this%dudx%init(this%u%dof, 'dudx')
+       call this%dudy%init(this%u%dof, 'dudy')
+       call this%dudz%init(this%u%dof, 'dudz')
+       call this%dvdx%init(this%u%dof, 'dvdx')
+       call this%dvdy%init(this%u%dof, 'dvdy')
+       call this%dvdz%init(this%u%dof, 'dvdz')
+       call this%dwdx%init(this%u%dof, 'dwdx')
+       call this%dwdy%init(this%u%dof, 'dwdy')
+       call this%dwdz%init(this%u%dof, 'dwdz')
+
+       call this%uuu%init(this%stats_work, 'uuu')
+       call this%vvv%init(this%stats_work, 'vvv')
+       call this%www%init(this%stats_work, 'www')
+       call this%uuv%init(this%stats_work, 'uuv')
+       call this%uuw%init(this%stats_work, 'uuw')
+       call this%uvv%init(this%stats_work, 'uvv')
+       call this%uvw%init(this%stats_work, 'uvw')
+       call this%vvw%init(this%stats_work, 'vvw')
+       call this%uww%init(this%stats_work, 'uww')
+       call this%vww%init(this%stats_work, 'vww')
+       call this%uuuu%init(this%stats_work, 'uuuu')
+       call this%vvvv%init(this%stats_work, 'vvvv')
+       call this%wwww%init(this%stats_work, 'wwww')
+       !> Pressure
+       call this%ppp%init(this%stats_work, 'ppp')
+       call this%pppp%init(this%stats_work, 'pppp')
+       !> Pressure * velocity
+       call this%pu%init(this%stats_work, 'pu')
+       call this%pv%init(this%stats_work, 'pv')
+       call this%pw%init(this%stats_work, 'pw')
+
+       call this%pdudx%init(this%stats_work, 'pdudx')
+       call this%pdudy%init(this%stats_work, 'pdudy')
+       call this%pdudz%init(this%stats_work, 'pdudz')
+       call this%pdvdx%init(this%stats_work, 'pdvdx')
+       call this%pdvdy%init(this%stats_work, 'pdvdy')
+       call this%pdvdz%init(this%stats_work, 'pdvdz')
+       call this%pdwdx%init(this%stats_work, 'pdwdx')
+       call this%pdwdy%init(this%stats_work, 'pdwdy')
+       call this%pdwdz%init(this%stats_work, 'pdwdz')
+
+       call this%e11%init(this%stats_work, 'e11')
+       call this%e22%init(this%stats_work, 'e22')
+       call this%e33%init(this%stats_work, 'e33')
+       call this%e12%init(this%stats_work, 'e12')
+       call this%e13%init(this%stats_work, 'e13')
+       call this%e23%init(this%stats_work, 'e23')
+    end if
 
     allocate(this%stat_fields%items(this%n_stats))
 
-    call this%stat_fields%assign_to_field(1, this%pp%mf)
-    call this%stat_fields%assign_to_field(2, this%uu%mf)
-    call this%stat_fields%assign_to_field(3, this%vv%mf)
-    call this%stat_fields%assign_to_field(4, this%ww%mf)
-    call this%stat_fields%assign_to_field(5, this%uv%mf)
-    call this%stat_fields%assign_to_field(6, this%uw%mf)
-    call this%stat_fields%assign_to_field(7, this%vw%mf)
-    call this%stat_fields%assign_to_field(8, this%uuu%mf) !< <uuu>
-    call this%stat_fields%assign_to_field(9, this%vvv%mf) !< <vvv>
-    call this%stat_fields%assign_to_field(10, this%www%mf) !< <www>
-    call this%stat_fields%assign_to_field(11, this%uuv%mf) !< <uuv>
-    call this%stat_fields%assign_to_field(12, this%uuw%mf) !< <uuw>
-    call this%stat_fields%assign_to_field(13, this%uvv%mf) !< <uvv>
-    call this%stat_fields%assign_to_field(14, this%uvw%mf) !< <uvv>
-    call this%stat_fields%assign_to_field(15, this%vvw%mf) !< <vvw>
-    call this%stat_fields%assign_to_field(16, this%uww%mf) !< <uww>
-    call this%stat_fields%assign_to_field(17, this%vww%mf) !< <vww>
-    call this%stat_fields%assign_to_field(18, this%uuuu%mf) !< <uuuu>
-    call this%stat_fields%assign_to_field(19, this%vvvv%mf) !< <vvvv>
-    call this%stat_fields%assign_to_field(20, this%wwww%mf) !< <wwww>
-    call this%stat_fields%assign_to_field(21, this%ppp%mf)
-    call this%stat_fields%assign_to_field(22, this%pppp%mf)
-    call this%stat_fields%assign_to_field(23, this%pu%mf)
-    call this%stat_fields%assign_to_field(24, this%pv%mf)
-    call this%stat_fields%assign_to_field(25, this%pw%mf)
-
-    call this%stat_fields%assign_to_field(26, this%pdudx%mf)
-    call this%stat_fields%assign_to_field(27, this%pdudy%mf)
-    call this%stat_fields%assign_to_field(28, this%pdudz%mf)
-    call this%stat_fields%assign_to_field(29, this%pdvdx%mf)
-    call this%stat_fields%assign_to_field(30, this%pdvdy%mf)
-    call this%stat_fields%assign_to_field(31, this%pdvdz%mf)
-    call this%stat_fields%assign_to_field(32, this%pdwdx%mf)
-    call this%stat_fields%assign_to_field(33, this%pdwdy%mf)
-    call this%stat_fields%assign_to_field(34, this%pdwdz%mf)
-    call this%stat_fields%assign_to_field(35, this%e11%mf)
-    call this%stat_fields%assign_to_field(36, this%e22%mf)
-    call this%stat_fields%assign_to_field(37, this%e33%mf)
-    call this%stat_fields%assign_to_field(38, this%e12%mf)
-    call this%stat_fields%assign_to_field(39, this%e13%mf)
-    call this%stat_fields%assign_to_field(40, this%e23%mf)
-
+    call this%stat_fields%assign_to_field(1, this%p_mean%mf)
+    call this%stat_fields%assign_to_field(2, this%u_mean%mf)
+    call this%stat_fields%assign_to_field(3, this%v_mean%mf)
+    call this%stat_fields%assign_to_field(4, this%w_mean%mf)
+    call this%stat_fields%assign_to_field(5, this%pp%mf)
+    call this%stat_fields%assign_to_field(6, this%uu%mf)
+    call this%stat_fields%assign_to_field(7, this%vv%mf)
+    call this%stat_fields%assign_to_field(8, this%ww%mf)
+    call this%stat_fields%assign_to_field(9, this%uv%mf)
+    call this%stat_fields%assign_to_field(10, this%uw%mf)
+    call this%stat_fields%assign_to_field(11, this%vw%mf)
+
+    if (this%n_stats .eq. 44) then
+       call this%stat_fields%assign_to_field(12, this%uuu%mf)
+       call this%stat_fields%assign_to_field(13, this%vvv%mf)
+       call this%stat_fields%assign_to_field(14, this%www%mf)
+       call this%stat_fields%assign_to_field(15, this%uuv%mf)
+       call this%stat_fields%assign_to_field(16, this%uuw%mf)
+       call this%stat_fields%assign_to_field(17, this%uvv%mf)
+       call this%stat_fields%assign_to_field(18, this%uvw%mf)
+       call this%stat_fields%assign_to_field(19, this%vvw%mf)
+       call this%stat_fields%assign_to_field(20, this%uww%mf)
+       call this%stat_fields%assign_to_field(21, this%vww%mf)
+       call this%stat_fields%assign_to_field(22, this%uuuu%mf)
+       call this%stat_fields%assign_to_field(23, this%vvvv%mf)
+       call this%stat_fields%assign_to_field(24, this%wwww%mf)
+       call this%stat_fields%assign_to_field(25, this%ppp%mf)
+       call this%stat_fields%assign_to_field(26, this%pppp%mf)
+       call this%stat_fields%assign_to_field(27, this%pu%mf)
+       call this%stat_fields%assign_to_field(28, this%pv%mf)
+       call this%stat_fields%assign_to_field(29, this%pw%mf)
+
+       call this%stat_fields%assign_to_field(30, this%pdudx%mf)
+       call this%stat_fields%assign_to_field(31, this%pdudy%mf)
+       call this%stat_fields%assign_to_field(32, this%pdudz%mf)
+       call this%stat_fields%assign_to_field(33, this%pdvdx%mf)
+       call this%stat_fields%assign_to_field(34, this%pdvdy%mf)
+       call this%stat_fields%assign_to_field(35, this%pdvdz%mf)
+       call this%stat_fields%assign_to_field(36, this%pdwdx%mf)
+       call this%stat_fields%assign_to_field(37, this%pdwdy%mf)
+       call this%stat_fields%assign_to_field(38, this%pdwdz%mf)
+       call this%stat_fields%assign_to_field(39, this%e11%mf)
+       call this%stat_fields%assign_to_field(40, this%e22%mf)
+       call this%stat_fields%assign_to_field(41, this%e33%mf)
+       call this%stat_fields%assign_to_field(42, this%e12%mf)
+       call this%stat_fields%assign_to_field(43, this%e13%mf)
+       call this%stat_fields%assign_to_field(44, this%e23%mf)
+    end if
 
   end subroutine fluid_stats_init
 
-  !> Updates all fields
+  !> Updates all fields with a new sample.
+  !! @param k Time elapsed since the last update.
   subroutine fluid_stats_update(this, k)
     class(fluid_stats_t), intent(inout) :: this
     real(kind=rp), intent(in) :: k
     integer :: n
 
-
-
     associate(stats_work => this%stats_work, stats_u => this%stats_u, &
-         stats_v => this%stats_v, stats_w => this%stats_w, &
-         stats_p => this%stats_p)
+              stats_v => this%stats_v, stats_w => this%stats_w, &
+              stats_p => this%stats_p)
       n = stats_work%dof%size()
 
       !> U%f is u and U%mf is <u>
       if (NEKO_BCKND_DEVICE .eq. 1) then
 
-         call device_col3(stats_u%x_d, this%u%x_d, this%u%x_d,n)
-         call device_col3(stats_v%x_d, this%v%x_d, this%v%x_d,n)
-         call device_col3(stats_w%x_d, this%w%x_d, this%w%x_d,n)
-         call device_col3(stats_p%x_d, this%p%x_d, this%p%x_d,n)
+         call this%u_mean%update(k)
+         call this%v_mean%update(k)
+         call this%w_mean%update(k)
+         call this%p_mean%update(k)
+
+         call device_col3(stats_u%x_d, this%u%x_d, this%u%x_d, n)
+         call device_col3(stats_v%x_d, this%v%x_d, this%v%x_d, n)
+         call device_col3(stats_w%x_d, this%w%x_d, this%w%x_d, n)
+         call device_col3(stats_p%x_d, this%p%x_d, this%p%x_d, n)
 
          call this%uu%update(k)
          call this%vv%update(k)
          call this%ww%update(k)
          call this%pp%update(k)
 
-         call device_col3(stats_work%x_d, this%u%x_d, this%v%x_d,n)
+         call device_col3(stats_work%x_d, this%u%x_d, this%v%x_d, n)
          call this%uv%update(k)
-         call device_col3(stats_work%x_d, this%u%x_d, this%w%x_d,n)
+         call device_col3(stats_work%x_d, this%u%x_d, this%w%x_d, n)
          call this%uw%update(k)
-         call device_col3(stats_work%x_d, this%v%x_d, this%w%x_d,n)
+         call device_col3(stats_work%x_d, this%v%x_d, this%w%x_d, n)
          call this%vw%update(k)
-
-         call device_col2(stats_work%x_d, this%u%x_d,n)
+         if (this%n_stats .eq. 11) return
+         call device_col2(stats_work%x_d, this%u%x_d, n)
          call this%uvw%update(k)
-         call device_col3(stats_work%x_d, this%stats_u%x_d, this%u%x_d,n)
+         call device_col3(stats_work%x_d, this%stats_u%x_d, this%u%x_d, n)
          call this%uuu%update(k)
-         call device_col3(stats_work%x_d, this%stats_v%x_d, this%v%x_d,n)
+         call device_col3(stats_work%x_d, this%stats_v%x_d, this%v%x_d, n)
          call this%vvv%update(k)
-         call device_col3(stats_work%x_d, this%stats_w%x_d, this%w%x_d,n)
+         call device_col3(stats_work%x_d, this%stats_w%x_d, this%w%x_d, n)
          call this%www%update(k)
-         call device_col3(stats_work%x_d, this%stats_u%x_d, this%v%x_d,n)
+         call device_col3(stats_work%x_d, this%stats_u%x_d, this%v%x_d, n)
          call this%uuv%update(k)
-         call device_col3(stats_work%x_d, this%stats_u%x_d, this%w%x_d,n)
+         call device_col3(stats_work%x_d, this%stats_u%x_d, this%w%x_d, n)
          call this%uuw%update(k)
-         call device_col3(stats_work%x_d, this%stats_v%x_d, this%u%x_d,n)
+         call device_col3(stats_work%x_d, this%stats_v%x_d, this%u%x_d, n)
          call this%uvv%update(k)
-         call device_col3(stats_work%x_d, this%stats_v%x_d, this%w%x_d,n)
+         call device_col3(stats_work%x_d, this%stats_v%x_d, this%w%x_d, n)
          call this%vvw%update(k)
-         call device_col3(stats_work%x_d, this%stats_w%x_d, this%u%x_d,n)
+         call device_col3(stats_work%x_d, this%stats_w%x_d, this%u%x_d, n)
          call this%uww%update(k)
-         call device_col3(stats_work%x_d, this%stats_w%x_d, this%v%x_d,n)
+         call device_col3(stats_work%x_d, this%stats_w%x_d, this%v%x_d, n)
          call this%vww%update(k)
 
-         call device_col3(stats_work%x_d, this%stats_u%x_d, this%stats_u%x_d,n)
+         call device_col3(stats_work%x_d, this%stats_u%x_d, this%stats_u%x_d, n)
          call this%uuuu%update(k)
-         call device_col3(stats_work%x_d, this%stats_v%x_d, this%stats_v%x_d,n)
+         call device_col3(stats_work%x_d, this%stats_v%x_d, this%stats_v%x_d, n)
          call this%vvvv%update(k)
-         call device_col3(stats_work%x_d, this%stats_w%x_d, this%stats_w%x_d,n)
+         call device_col3(stats_work%x_d, this%stats_w%x_d, this%stats_w%x_d, n)
          call this%wwww%update(k)
 
-         call device_col3(stats_work%x_d, this%stats_p%x_d, this%p%x_d,n)
+         call device_col3(stats_work%x_d, this%stats_p%x_d, this%p%x_d, n)
          call this%ppp%update(k)
-         call device_col3(stats_work%x_d, this%stats_p%x_d, this%stats_p%x_d,n)
+         call device_col3(stats_work%x_d, this%stats_p%x_d, this%stats_p%x_d, n)
          call this%pppp%update(k)
 
-         call device_col3(stats_work%x_d, this%p%x_d, this%u%x_d,n)
+         call device_col3(stats_work%x_d, this%p%x_d, this%u%x_d, n)
          call this%pu%update(k)
-         call device_col3(stats_work%x_d, this%p%x_d, this%v%x_d,n)
+         call device_col3(stats_work%x_d, this%p%x_d, this%v%x_d, n)
          call this%pv%update(k)
-         call device_col3(stats_work%x_d, this%p%x_d, this%w%x_d,n)
+         call device_col3(stats_work%x_d, this%p%x_d, this%w%x_d, n)
          call this%pw%update(k)
 
       else
 
-         call col3(stats_u%x, this%u%x, this%u%x,n)
-         call col3(stats_v%x, this%v%x, this%v%x,n)
-         call col3(stats_w%x, this%w%x, this%w%x,n)
-         call col3(stats_p%x, this%p%x, this%p%x,n)
+         call this%u_mean%update(k)
+         call this%v_mean%update(k)
+         call this%w_mean%update(k)
+         call this%p_mean%update(k)
+         call col3(stats_u%x, this%u%x, this%u%x, n)
+         call col3(stats_v%x, this%v%x, this%v%x, n)
+         call col3(stats_w%x, this%w%x, this%w%x, n)
+         call col3(stats_p%x, this%p%x, this%p%x, n)
 
          call this%uu%update(k)
          call this%vv%update(k)
          call this%ww%update(k)
          call this%pp%update(k)
 
-         call col3(stats_work%x, this%u%x, this%v%x,n)
+         call col3(stats_work%x, this%u%x, this%v%x, n)
          call this%uv%update(k)
-         call col3(stats_work%x, this%u%x, this%w%x,n)
+         call col3(stats_work%x, this%u%x, this%w%x, n)
          call this%uw%update(k)
-         call col3(stats_work%x, this%v%x, this%w%x,n)
+         call col3(stats_work%x, this%v%x, this%w%x, n)
          call this%vw%update(k)
 
-         call col2(stats_work%x, this%u%x,n)
+         if (this%n_stats .eq. 11) return
+
+         call col2(stats_work%x, this%u%x, n)
          call this%uvw%update(k)
-         call col3(stats_work%x, this%stats_u%x, this%u%x,n)
+         call col3(stats_work%x, this%stats_u%x, this%u%x, n)
          call this%uuu%update(k)
-         call col3(stats_work%x, this%stats_v%x, this%v%x,n)
+         call col3(stats_work%x, this%stats_v%x, this%v%x, n)
          call this%vvv%update(k)
-         call col3(stats_work%x, this%stats_w%x, this%w%x,n)
+         call col3(stats_work%x, this%stats_w%x, this%w%x, n)
          call this%www%update(k)
-         call col3(stats_work%x, this%stats_u%x, this%v%x,n)
+         call col3(stats_work%x, this%stats_u%x, this%v%x, n)
          call this%uuv%update(k)
-         call col3(stats_work%x, this%stats_u%x, this%w%x,n)
+         call col3(stats_work%x, this%stats_u%x, this%w%x, n)
          call this%uuw%update(k)
-         call col3(stats_work%x, this%stats_v%x, this%u%x,n)
+         call col3(stats_work%x, this%stats_v%x, this%u%x, n)
          call this%uvv%update(k)
-         call col3(stats_work%x, this%stats_v%x, this%w%x,n)
+         call col3(stats_work%x, this%stats_v%x, this%w%x, n)
          call this%vvw%update(k)
-         call col3(stats_work%x, this%stats_w%x, this%u%x,n)
+         call col3(stats_work%x, this%stats_w%x, this%u%x, n)
          call this%uww%update(k)
-         call col3(stats_work%x, this%stats_w%x, this%v%x,n)
+         call col3(stats_work%x, this%stats_w%x, this%v%x, n)
          call this%vww%update(k)
 
-         call col3(stats_work%x, this%stats_u%x, this%stats_u%x,n)
+         call col3(stats_work%x, this%stats_u%x, this%stats_u%x, n)
          call this%uuuu%update(k)
-         call col3(stats_work%x, this%stats_v%x, this%stats_v%x,n)
+         call col3(stats_work%x, this%stats_v%x, this%stats_v%x, n)
          call this%vvvv%update(k)
-         call col3(stats_work%x, this%stats_w%x, this%stats_w%x,n)
+         call col3(stats_work%x, this%stats_w%x, this%stats_w%x, n)
          call this%wwww%update(k)
 
-         call col3(stats_work%x, this%stats_p%x, this%p%x,n)
+         call col3(stats_work%x, this%stats_p%x, this%p%x, n)
          call this%ppp%update(k)
-         call col3(stats_work%x, this%stats_p%x, this%stats_p%x,n)
+         call col3(stats_work%x, this%stats_p%x, this%stats_p%x, n)
          call this%pppp%update(k)
 
          call col3(stats_work%x, this%p%x, this%u%x,n)
@@ -411,111 +457,118 @@ subroutine fluid_stats_update(this, k)
       call opgrad(this%dwdx%x, this%dwdy%x, this%dwdz%x, this%w%x, this%coef)
 
       if (NEKO_BCKND_DEVICE .eq. 1) then
-         call device_col3(stats_work%x_d, this%dudx%x_d, this%p%x_d,n)
+         call device_col3(stats_work%x_d, this%dudx%x_d, this%p%x_d, n)
          call this%pdudx%update(k)
-         call device_col3(stats_work%x_d, this%dudy%x_d, this%p%x_d,n)
+         call device_col3(stats_work%x_d, this%dudy%x_d, this%p%x_d, n)
          call this%pdudy%update(k)
-         call device_col3(stats_work%x_d, this%dudz%x_d, this%p%x_d,n)
+         call device_col3(stats_work%x_d, this%dudz%x_d, this%p%x_d, n)
          call this%pdudz%update(k)
 
-         call device_col3(stats_work%x_d, this%dvdx%x_d, this%p%x_d,n)
+         call device_col3(stats_work%x_d, this%dvdx%x_d, this%p%x_d, n)
          call this%pdvdx%update(k)
-         call device_col3(stats_work%x_d, this%dvdy%x_d, this%p%x_d,n)
+         call device_col3(stats_work%x_d, this%dvdy%x_d, this%p%x_d, n)
          call this%pdvdy%update(k)
-         call device_col3(stats_work%x_d, this%dvdz%x_d, this%p%x_d,n)
+         call device_col3(stats_work%x_d, this%dvdz%x_d, this%p%x_d, n)
          call this%pdvdz%update(k)
 
-         call device_col3(stats_work%x_d, this%dwdx%x_d, this%p%x_d,n)
+         call device_col3(stats_work%x_d, this%dwdx%x_d, this%p%x_d, n)
          call this%pdwdx%update(k)
-         call device_col3(stats_work%x_d, this%dwdy%x_d, this%p%x_d,n)
+         call device_col3(stats_work%x_d, this%dwdy%x_d, this%p%x_d, n)
          call this%pdwdy%update(k)
-         call device_col3(stats_work%x_d, this%dwdz%x_d, this%p%x_d,n)
+         call device_col3(stats_work%x_d, this%dwdz%x_d, this%p%x_d, n)
          call this%pdwdz%update(k)
 
-         call device_col3(this%stats_work%x_d, this%dudx%x_d, this%dudx%x_d,n)
-         call device_addcol3(this%stats_work%x_d, &
-              this%dudy%x_d, this%dudy%x_d,n)
-         call device_addcol3(this%stats_work%x_d, &
-              this%dudz%x_d, this%dudz%x_d,n)
+         call device_col3(this%stats_work%x_d, this%dudx%x_d, this%dudx%x_d, n)
+         call device_addcol3(this%stats_work%x_d, this%dudy%x_d, &
+              this%dudy%x_d, n)
+         call device_addcol3(this%stats_work%x_d, this%dudz%x_d, &
+              this%dudz%x_d, n)
          call this%e11%update(k)
-         call device_col3(this%stats_work%x_d, &
-              this%dvdx%x_d, this%dvdx%x_d,n)
-         call device_addcol3(this%stats_work%x_d, &
-              this%dvdy%x_d, this%dvdy%x_d,n)
-         call device_addcol3(this%stats_work%x_d, &
-              this%dvdz%x_d, this%dvdz%x_d,n)
+         call device_col3(this%stats_work%x_d, this%dvdx%x_d, this%dvdx%x_d, n)
+         call device_addcol3(this%stats_work%x_d, this%dvdy%x_d, & 
+              this%dvdy%x_d, n)
+         call device_addcol3(this%stats_work%x_d, this%dvdz%x_d, & 
+              this%dvdz%x_d, n)
          call this%e22%update(k)
-         call device_col3(this%stats_work%x_d, this%dwdx%x_d, this%dwdx%x_d,n)
-         call device_addcol3(this%stats_work%x_d, &
-              this%dwdy%x_d, this%dwdy%x_d,n)
-         call device_addcol3(this%stats_work%x_d, &
-              this%dwdz%x_d, this%dwdz%x_d,n)
+         call device_col3(this%stats_work%x_d, this%dwdx%x_d, this%dwdx%x_d, n)
+         call device_addcol3(this%stats_work%x_d, this%dwdy%x_d, &
+              this%dwdy%x_d, n)
+         call device_addcol3(this%stats_work%x_d, this%dwdz%x_d, &
+              this%dwdz%x_d, n)
          call this%e33%update(k)
-         call device_col3(this%stats_work%x_d, this%dudx%x_d, this%dvdx%x_d,n)
-         call device_addcol3(this%stats_work%x_d, &
-              this%dudy%x_d, this%dvdy%x_d,n)
-         call device_addcol3(this%stats_work%x_d, &
-              this%dudz%x_d, this%dvdz%x_d,n)
+         call device_col3(this%stats_work%x_d, this%dudx%x_d, &
+              this%dvdx%x_d, n)
+         call device_addcol3(this%stats_work%x_d, this%dudy%x_d, &
+              this%dvdy%x_d, n)
+         call device_addcol3(this%stats_work%x_d, this%dudz%x_d, &
+              this%dvdz%x_d, n)
          call this%e12%update(k)
-         call device_col3(this%stats_work%x_d, this%dvdx%x_d, this%dwdx%x_d,n)
-         call device_addcol3(this%stats_work%x_d, &
-              this%dvdy%x_d, this%dwdy%x_d,n)
-         call device_addcol3(this%stats_work%x_d, &
-              this%dvdz%x_d, this%dwdz%x_d,n)
+         call device_col3(this%stats_work%x_d, this%dudx%x_d, this%dwdx%x_d, n)
+         call device_addcol3(this%stats_work%x_d, this%dudy%x_d, &
+              this%dwdy%x_d, n)
+         call device_addcol3(this%stats_work%x_d, this%dudz%x_d, &
+              this%dwdz%x_d, n)
+         call this%e13%update(k)
+         call device_col3(this%stats_work%x_d, this%dvdx%x_d, this%dwdx%x_d, n)
+         call device_addcol3(this%stats_work%x_d, this%dvdy%x_d, &
+              this%dwdy%x_d, n)
+         call device_addcol3(this%stats_work%x_d, this%dvdz%x_d, &
+              this%dwdz%x_d, n)
          call this%e23%update(k)
-
-
       else
-         call col3(stats_work%x, this%dudx%x, this%p%x,n)
+         call col3(stats_work%x, this%dudx%x, this%p%x, n)
          call this%pdudx%update(k)
-         call col3(stats_work%x, this%dudy%x, this%p%x,n)
+         call col3(stats_work%x, this%dudy%x, this%p%x, n)
          call this%pdudy%update(k)
-         call col3(stats_work%x, this%dudz%x, this%p%x,n)
+         call col3(stats_work%x, this%dudz%x, this%p%x, n)
          call this%pdudz%update(k)
 
-         call col3(stats_work%x, this%dvdx%x, this%p%x,n)
+         call col3(stats_work%x, this%dvdx%x, this%p%x, n)
          call this%pdvdx%update(k)
-         call col3(stats_work%x, this%dvdy%x, this%p%x,n)
+         call col3(stats_work%x, this%dvdy%x, this%p%x, n)
          call this%pdvdy%update(k)
-         call col3(stats_work%x, this%dvdz%x, this%p%x,n)
+         call col3(stats_work%x, this%dvdz%x, this%p%x, n)
          call this%pdvdz%update(k)
 
-         call col3(stats_work%x, this%dwdx%x, this%p%x,n)
+         call col3(stats_work%x, this%dwdx%x, this%p%x, n)
          call this%pdwdx%update(k)
-         call col3(stats_work%x, this%dwdy%x, this%p%x,n)
+         call col3(stats_work%x, this%dwdy%x, this%p%x, n)
          call this%pdwdy%update(k)
-         call col3(stats_work%x, this%dwdz%x, this%p%x,n)
+         call col3(stats_work%x, this%dwdz%x, this%p%x, n)
          call this%pdwdz%update(k)
 
-         call col3(this%stats_work%x, this%dudx%x, this%dudx%x,n)
-         call addcol3(this%stats_work%x, this%dudy%x, this%dudy%x,n)
-         call addcol3(this%stats_work%x, this%dudz%x, this%dudz%x,n)
+         call col3(this%stats_work%x, this%dudx%x, this%dudx%x, n)
+         call addcol3(this%stats_work%x, this%dudy%x, this%dudy%x, n)
+         call addcol3(this%stats_work%x, this%dudz%x, this%dudz%x, n)
          call this%e11%update(k)
-         call col3(this%stats_work%x, this%dvdx%x, this%dvdx%x,n)
-         call addcol3(this%stats_work%x, this%dvdy%x, this%dvdy%x,n)
-         call addcol3(this%stats_work%x, this%dvdz%x, this%dvdz%x,n)
+         call col3(this%stats_work%x, this%dvdx%x, this%dvdx%x, n)
+         call addcol3(this%stats_work%x, this%dvdy%x, this%dvdy%x, n)
+         call addcol3(this%stats_work%x, this%dvdz%x, this%dvdz%x, n)
          call this%e22%update(k)
-         call col3(this%stats_work%x, this%dwdx%x, this%dwdx%x,n)
-         call addcol3(this%stats_work%x, this%dwdy%x, this%dwdy%x,n)
-         call addcol3(this%stats_work%x, this%dwdz%x, this%dwdz%x,n)
+         call col3(this%stats_work%x, this%dwdx%x, this%dwdx%x, n)
+         call addcol3(this%stats_work%x, this%dwdy%x, this%dwdy%x, n)
+         call addcol3(this%stats_work%x, this%dwdz%x, this%dwdz%x, n)
          call this%e33%update(k)
-         call col3(this%stats_work%x, this%dudx%x, this%dvdx%x,n)
-         call addcol3(this%stats_work%x, this%dudy%x, this%dvdy%x,n)
-         call addcol3(this%stats_work%x, this%dudz%x, this%dvdz%x,n)
+         call col3(this%stats_work%x, this%dudx%x, this%dvdx%x, n)
+         call addcol3(this%stats_work%x, this%dudy%x, this%dvdy%x, n)
+         call addcol3(this%stats_work%x, this%dudz%x, this%dvdz%x, n)
          call this%e12%update(k)
-         call col3(this%stats_work%x, this%dvdx%x, this%dwdx%x,n)
-         call addcol3(this%stats_work%x, this%dvdy%x, this%dwdy%x,n)
-         call addcol3(this%stats_work%x, this%dvdz%x, this%dwdz%x,n)
+         call col3(this%stats_work%x,this%dudx%x, this%dwdx%x,n)
+         call addcol3(this%stats_work%x,this%dudy%x, this%dwdy%x,n)
+         call addcol3(this%stats_work%x,this%dudz%x, this%dwdz%x,n)
+         call this%e13%update(k)
+         call col3(this%stats_work%x, this%dvdx%x, this%dwdx%x, n)
+         call addcol3(this%stats_work%x, this%dvdy%x, this%dwdy%x, n)
+         call addcol3(this%stats_work%x, this%dvdz%x, this%dwdz%x, n)
          call this%e23%update(k)
 
       end if
-
     end associate
 
   end subroutine fluid_stats_update
 
 
-  !> Deallocates a mean flow field
+  !> Destructor.
   subroutine fluid_stats_free(this)
     class(fluid_stats_t), intent(inout) :: this
 
@@ -524,6 +577,11 @@ subroutine fluid_stats_free(this)
     call this%stats_v%free()
     call this%stats_w%free()
 
+    call this%u_mean%free()
+    call this%v_mean%free()
+    call this%w_mean%free()
+    call this%p_mean%free()
+
     call this%uu%free()
     call this%vv%free()
     call this%ww%free()
@@ -544,9 +602,14 @@ subroutine fluid_stats_free(this)
 
   end subroutine fluid_stats_free
 
-  !> Initialize a mean flow field
+  !> Resets all the computed means values and sampling times to zero.
   subroutine fluid_stats_reset(this)
-    class(fluid_stats_t), intent(inout), target :: this
+    class(fluid_stats_t), intent(inout), target:: this
+
+    call this%p_mean%reset()
+    call this%u_mean%reset()
+    call this%v_mean%reset()
+    call this%w_mean%reset()
 
     call this%uu%reset()
     call this%vv%reset()
@@ -554,54 +617,59 @@ subroutine fluid_stats_reset(this)
     call this%uv%reset()
     call this%uw%reset()
     call this%vw%reset()
-    call this%uuu%reset()
-    call this%vvv%reset()
-    call this%www%reset()
-    call this%uuv%reset()
-    call this%uuw%reset()
-    call this%uvv%reset()
-    call this%uvw%reset()
-    call this%vvw%reset()
-    call this%uww%reset()
-    call this%vww%reset()
-    call this%uuuu%reset()
-    call this%vvvv%reset()
-    call this%wwww%reset()
-    !> Pressure
     call this%pp%reset()
-    call this%ppp%reset()
-    call this%pppp%reset()
-    !> Pressure * velocity
-    call this%pu%reset()
-    call this%pv%reset()
-    call this%pw%reset()
-
-    call this%pdudx%reset()
-    call this%pdudy%reset()
-    call this%pdudz%reset()
-    call this%pdvdx%reset()
-    call this%pdvdy%reset()
-    call this%pdvdz%reset()
-    call this%pdwdx%reset()
-    call this%pdwdy%reset()
-    call this%pdwdz%reset()
-
-    call this%e11%reset()
-    call this%e22%reset()
-    call this%e33%reset()
-    call this%e12%reset()
-    call this%e13%reset()
-    call this%e23%reset()
+    if (this%n_stats .eq. 44) then
+       call this%uuu%reset()
+       call this%vvv%reset()
+       call this%www%reset()
+       call this%uuv%reset()
+       call this%uuw%reset()
+       call this%uvv%reset()
+       call this%uvw%reset()
+       call this%vvw%reset()
+       call this%uww%reset()
+       call this%vww%reset()
+       call this%uuuu%reset()
+       call this%vvvv%reset()
+       call this%wwww%reset()
+       call this%ppp%reset()
+       call this%pppp%reset()
+       call this%pu%reset()
+       call this%pv%reset()
+       call this%pw%reset()
+
+       call this%pdudx%reset()
+       call this%pdudy%reset()
+       call this%pdudz%reset()
+       call this%pdvdx%reset()
+       call this%pdvdy%reset()
+       call this%pdvdz%reset()
+       call this%pdwdx%reset()
+       call this%pdwdy%reset()
+       call this%pdwdz%reset()
+
+       call this%e11%reset()
+       call this%e22%reset()
+       call this%e33%reset()
+       call this%e12%reset()
+       call this%e13%reset()
+       call this%e23%reset()
+    end if
 
   end subroutine fluid_stats_reset
 
+  ! Convert computed weak gradients to strong.
   subroutine fluid_stats_make_strong_grad(this)
     class(fluid_stats_t) :: this
     integer :: n
+ 
+    if (this%n_stats .eq. 11) return
+ 
     n = size(this%coef%B)
+ 
     if (NEKO_BCKND_DEVICE .eq. 1) then
-       call device_cfill(this%stats_work%x_d, 1.0_rp,n)
-       call device_invcol2(this%stats_work%x_d, this%coef%B_d,n)
+       call device_cfill(this%stats_work%x_d, 1.0_rp, n)
+       call device_invcol2(this%stats_work%x_d, this%coef%B_d, n)
        call device_col2(this%pdudx%mf%x_d, this%stats_work%x_d, n)
        call device_col2(this%pdudy%mf%x_d, this%stats_work%x_d, n)
        call device_col2(this%pdudz%mf%x_d, this%stats_work%x_d, n)
@@ -612,7 +680,7 @@ subroutine fluid_stats_make_strong_grad(this)
        call device_col2(this%pdwdy%mf%x_d, this%stats_work%x_d, n)
        call device_col2(this%pdwdz%mf%x_d, this%stats_work%x_d, n)
 
-       call device_col2(this%stats_work%x_d, this%stats_work%x_d,n)
+       call device_col2(this%stats_work%x_d, this%stats_work%x_d, n)
        call device_col2(this%e11%mf%x_d, this%stats_work%x_d, n)
        call device_col2(this%e22%mf%x_d, this%stats_work%x_d, n)
        call device_col2(this%e33%mf%x_d, this%stats_work%x_d, n)
@@ -622,7 +690,7 @@ subroutine fluid_stats_make_strong_grad(this)
 
 
     else
-       call invers2(this%stats_work%x, this%coef%B,n)
+       call invers2(this%stats_work%x, this%coef%B, n)
        call col2(this%pdudx%mf%x, this%stats_work%x, n)
        call col2(this%pdudy%mf%x, this%stats_work%x, n)
        call col2(this%pdudz%mf%x, this%stats_work%x, n)
@@ -633,7 +701,7 @@ subroutine fluid_stats_make_strong_grad(this)
        call col2(this%pdwdy%mf%x, this%stats_work%x, n)
        call col2(this%pdwdz%mf%x, this%stats_work%x, n)
 
-       call col2(this%stats_work%x, this%stats_work%x,n)
+       call col2(this%stats_work%x, this%stats_work%x, n)
        call col2(this%e11%mf%x, this%stats_work%x, n)
        call col2(this%e22%mf%x, this%stats_work%x, n)
        call col2(this%e33%mf%x, this%stats_work%x, n)
@@ -645,9 +713,10 @@ subroutine fluid_stats_make_strong_grad(this)
 
   end subroutine fluid_stats_make_strong_grad
 
-  subroutine fluid_stats_post_process(this, mean, reynolds, pressure_flatness, &
-       pressure_skewness, skewness_tensor, &
-       mean_vel_grad, dissipation_tensor)
+  !> Compute certain physical statistical quantities based on existing mean 
+  !! fields.
+  subroutine fluid_stats_post_process(this, mean, reynolds, pressure_flatness,&
+       pressure_skewness, skewness_tensor, mean_vel_grad, dissipation_tensor)
     class(fluid_stats_t) :: this
     type(field_list_t), intent(inout), optional :: mean
     type(field_list_t), intent(inout), optional :: reynolds
@@ -660,69 +729,76 @@ subroutine fluid_stats_post_process(this, mean, reynolds, pressure_flatness, &
 
     if (present(mean)) then
        n = mean%item_size(1)
-       call copy(mean%items(1)%ptr%x, this%u_mean%x, n)
-       call copy(mean%items(2)%ptr%x, this%v_mean%x, n)
-       call copy(mean%items(3)%ptr%x, this%w_mean%x, n)
-       call copy(mean%items(4)%ptr%x, this%p_mean%x, n)
+       call copy(mean%items(1)%ptr%x, this%u_mean%mf%x, n)
+       call copy(mean%items(2)%ptr%x, this%v_mean%mf%x, n)
+       call copy(mean%items(3)%ptr%x, this%w_mean%mf%x, n)
+       call copy(mean%items(4)%ptr%x, this%p_mean%mf%x, n)
     end if
 
     if (present(reynolds)) then
        n = reynolds%item_size(1)
        call copy(reynolds%items(1)%ptr%x, this%pp%mf%x, n)
-       call subcol3(reynolds%items(1)%ptr%x, this%p_mean%x, this%p_mean%x, n)
+       call subcol3(reynolds%items(1)%ptr%x, this%p_mean%mf%x, &
+            this%p_mean%mf%x, n)
 
        call copy(reynolds%items(2)%ptr%x, this%uu%mf%x, n)
-       call subcol3(reynolds%items(2)%ptr%x, this%u_mean%x, this%u_mean%x, n)
+       call subcol3(reynolds%items(2)%ptr%x, this%u_mean%mf%x, &
+            this%u_mean%mf%x, n)
 
        call copy(reynolds%items(3)%ptr%x, this%vv%mf%x, n)
-       call subcol3(reynolds%items(3)%ptr%x, this%v_mean%x, this%v_mean%x,n)
+       call subcol3(reynolds%items(3)%ptr%x, this%v_mean%mf%x, &
+            this%v_mean%mf%x,n)
 
        call copy(reynolds%items(4)%ptr%x, this%ww%mf%x, n)
-       call subcol3(reynolds%items(4)%ptr%x, this%w_mean%x, this%w_mean%x,n)
+       call subcol3(reynolds%items(4)%ptr%x, this%w_mean%mf%x, &
+            this%w_mean%mf%x,n)
 
        call copy(reynolds%items(5)%ptr%x, this%uv%mf%x, n)
-       call subcol3(reynolds%items(5)%ptr%x, this%u_mean%x, this%v_mean%x, n)
+       call subcol3(reynolds%items(5)%ptr%x, this%u_mean%mf%x, &
+            this%v_mean%mf%x, n)
 
        call copy(reynolds%items(6)%ptr%x, this%uw%mf%x, n)
-       call subcol3(reynolds%items(6)%ptr%x, this%u_mean%x, this%w_mean%x, n)
+       call subcol3(reynolds%items(6)%ptr%x, this%u_mean%mf%x, &
+            this%w_mean%mf%x, n)
 
        call copy(reynolds%items(7)%ptr%x, this%vw%mf%x, n)
-       call subcol3(reynolds%items(7)%ptr%x, this%v_mean%x, this%w_mean%x, n)
+       call subcol3(reynolds%items(7)%ptr%x, this%v_mean%mf%x, &
+            this%w_mean%mf%x, n)
     end if
     if (present(pressure_skewness)) then
 
-       call neko_warning('Presssure skewness stat not implemented in &
-            &fluid_stats yet, please help!')
+       call neko_warning('Presssure skewness stat not implemented'// &
+                         ' in fluid_stats, process stats in python instead')
 
     end if
 
     if (present(pressure_flatness)) then
-       call neko_warning('Presssure flatness stat not implemented yet, &
-            &please help!')
+       call neko_warning('Presssure flatness stat not implemented'// &
+                         ' in fluid_stats, process stats in python instead')
 
     end if
 
     if (present(skewness_tensor)) then
-       call neko_warning('Skewness tensor stat not implemented yet, &
-            &please help!')
+       call neko_warning('Skewness tensor stat not implemented'// &
+                         ' in fluid_stats, process stats in python instead')
     end if
 
     if (present(mean_vel_grad)) then
        !Compute gradient of mean flow
        n = mean_vel_grad%item_size(1)
        if (NEKO_BCKND_DEVICE .eq. 1) then
-          call device_memcpy(this%u_mean%x, this%u_mean%x_d, n, &
+          call device_memcpy(this%u_mean%mf%x, this%u_mean%mf%x_d, n, &
                              HOST_TO_DEVICE, sync = .false.)
-          call device_memcpy(this%v_mean%x, this%v_mean%x_d, n, &
+          call device_memcpy(this%v_mean%mf%x, this%v_mean%mf%x_d, n, &
                              HOST_TO_DEVICE, sync = .false.)
-          call device_memcpy(this%w_mean%x, this%w_mean%x_d, n, &
+          call device_memcpy(this%w_mean%mf%x, this%w_mean%mf%x_d, n, &
                              HOST_TO_DEVICE, sync = .false.)
           call opgrad(this%dudx%x, this%dudy%x, this%dudz%x, &
-                      this%u_mean%x, this%coef)
+                      this%u_mean%mf%x, this%coef)
           call opgrad(this%dvdx%x, this%dvdy%x, this%dvdz%x, &
-                      this%v_mean%x, this%coef)
+                      this%v_mean%mf%x, this%coef)
           call opgrad(this%dwdx%x, this%dwdy%x, this%dwdz%x, &
-                      this%w_mean%x, this%coef)
+                      this%w_mean%mf%x, this%coef)
           call device_memcpy(this%dudx%x, this%dudx%x_d, n, &
                              DEVICE_TO_HOST, sync = .false.)
           call device_memcpy(this%dvdx%x, this%dvdx%x_d, n, &
@@ -743,36 +819,37 @@ subroutine fluid_stats_post_process(this, mean, reynolds, pressure_flatness, &
                              DEVICE_TO_HOST, sync = .true.)
        else
           call opgrad(this%dudx%x, this%dudy%x, this%dudz%x, &
-               this%u_mean%x, this%coef)
+                      this%u_mean%mf%x, this%coef)
           call opgrad(this%dvdx%x, this%dvdy%x, this%dvdz%x, &
-               this%v_mean%x, this%coef)
-          call opgrad(this%dwdx%x, this%dwdy%x, this%dwdz%x, &
-               this%w_mean%x, this%coef)
+                      this%v_mean%mf%x, this%coef)
+          call opgrad(this%dwdx%x, this%dwdy%x, this%dwdz%x, & 
+                      this%w_mean%mf%x, this%coef)
        end if
        call invers2(this%stats_work%x, this%coef%B,n)
-       call col3(mean_vel_grad%items(1)%ptr%x, this%dudx%x, &
-            this%stats_work%x, n)
+       call col3(mean_vel_grad%items(1)%ptr%x, this%dudx%x, & 
+                 this%stats_work%x, n)
        call col3(mean_vel_grad%items(2)%ptr%x, this%dudy%x, &
-            this%stats_work%x, n)
+                 this%stats_work%x, n)
        call col3(mean_vel_grad%items(3)%ptr%x, this%dudz%x, &
-            this%stats_work%x, n)
+                 this%stats_work%x, n)
        call col3(mean_vel_grad%items(4)%ptr%x, this%dvdx%x, &
-            this%stats_work%x, n)
+                 this%stats_work%x, n)
        call col3(mean_vel_grad%items(5)%ptr%x, this%dvdy%x, &
-            this%stats_work%x, n)
+                 this%stats_work%x, n)
        call col3(mean_vel_grad%items(6)%ptr%x, this%dvdz%x, &
-            this%stats_work%x, n)
+                 this%stats_work%x, n)
        call col3(mean_vel_grad%items(7)%ptr%x, this%dwdx%x, &
-            this%stats_work%x, n)
+                 this%stats_work%x, n)
        call col3(mean_vel_grad%items(8)%ptr%x, this%dwdy%x, &
-            this%stats_work%x, n)
+                 this%stats_work%x, n)
        call col3(mean_vel_grad%items(9)%ptr%x, this%dwdz%x, &
-            this%stats_work%x, n)
+                 this%stats_work%x, n)
 
     end if
 
     if (present(dissipation_tensor)) then
-
+       call neko_warning('Dissipation tensor stat not implemented'// &
+                         ' in fluid_stats, process stats in python instead')
     end if
 
   end subroutine fluid_stats_post_process
diff --git a/src/io/fld_file.f90 b/src/io/fld_file.f90
index 0a204a826e7..b2243ad41fd 100644
--- a/src/io/fld_file.f90
+++ b/src/io/fld_file.f90
@@ -79,6 +79,7 @@ subroutine fld_file_write(this, data, t)
     class(*), target, intent(in) :: data
     real(kind=rp), intent(in), optional :: t
     type(array_ptr_t) :: x, y, z, u, v, w, p, tem
+    real(kind=rp), allocatable, target :: tempo(:)
     type(mesh_t), pointer :: msh
     type(dofmap_t), pointer :: dof
     type(space_t), pointer :: Xh
@@ -126,6 +127,7 @@ subroutine fld_file_write(this, data, t)
        if (data%x%n .gt. 0) x%ptr => data%x%x
        if (data%y%n .gt. 0) y%ptr => data%y%x
        if (data%z%n .gt. 0) z%ptr => data%z%x
+       if (gdim .eq. 2) z%ptr => data%y%x
        if (data%u%n .gt. 0) then
           u%ptr => data%u%x
           write_velocity = .true.
@@ -140,7 +142,8 @@ subroutine fld_file_write(this, data, t)
           write_temperature = .true.
           tem%ptr => data%t%x
        end if
-       
+       ! If gdim = 2 and Z-velocity component exists,
+       ! it is stored in last scalar field 
        if (gdim .eq. 2 .and. data%w%n .gt. 0) then
           n_scalar_fields = data%n_scalars + 1
           allocate(scalar_fields(n_scalar_fields))
@@ -156,6 +159,20 @@ subroutine fld_file_write(this, data, t)
           end do
           scalar_fields(n_scalar_fields+1)%ptr => data%w%x
        end if
+       ! This is very stupid...
+       ! Some compilers cannot handle that these pointers dont point to anything
+       ! (although they are not used) this fixes a segfault due to this.
+       if (nelv .eq. 0) then
+          allocate(tempo(1))
+          x%ptr => tempo
+          y%ptr => tempo
+          z%ptr => tempo
+          u%ptr => tempo
+          v%ptr => tempo
+          w%ptr => tempo
+          p%ptr => tempo
+          tem%ptr => tempo
+       end if
 
        allocate(idx(nelv))
        do i = 1, nelv
@@ -262,7 +279,6 @@ subroutine fld_file_write(this, data, t)
     else
        FLD_DATA_SIZE = MPI_REAL_SIZE
     end if
-
     if (this%dp_precision) then
        allocate(tmp_dp(gdim*n))
     else
@@ -341,15 +357,12 @@ subroutine fld_file_write(this, data, t)
     call MPI_File_write_at_all(fh, byte_offset, idx, nelv, &
          MPI_INTEGER, status, ierr)
     mpi_offset = mpi_offset + int(glb_nelv, i8) * int(MPI_INTEGER_SIZE, i8)
-
     deallocate(idx)
-
     if (write_mesh) then
 
        byte_offset = mpi_offset + int(offset_el, i8) * &
             (int(gdim*lxyz, i8) * &
             int(FLD_DATA_SIZE, i8))
-
        call fld_file_write_vector_field(this, fh, byte_offset, &
             x%ptr, y%ptr, z%ptr, &
             n, gdim, lxyz, nelv)
@@ -357,7 +370,6 @@ subroutine fld_file_write(this, data, t)
             (int(gdim *lxyz, i8) * &
             int(FLD_DATA_SIZE, i8))
     end if
-
     if (write_velocity) then
        byte_offset = mpi_offset + int(offset_el, i8) * &
             (int(gdim * (lxyz), i8) * int(FLD_DATA_SIZE, i8))
@@ -501,7 +513,6 @@ subroutine fld_file_write(this, data, t)
 
     call MPI_File_sync(fh, ierr)
     call MPI_File_close(fh, ierr)
-
     ! Write metadata file
     if (pe_rank .eq. 0) then
        tslash_pos = filename_tslash_pos(this%fname)
diff --git a/src/io/fluid_stats_output.f90 b/src/io/fluid_stats_output.f90
index ac3cd8df382..757f4d68450 100644
--- a/src/io/fluid_stats_output.f90
+++ b/src/io/fluid_stats_output.f90
@@ -1,4 +1,4 @@
-! Copyright (c) 2022, The Neko Authors
+! Copyright (c) 2024, The Neko Authors
 ! All rights reserved.
 !
 ! Redistribution and use in source and binary forms, with or without
@@ -30,68 +30,135 @@
 ! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 ! POSSIBILITY OF SUCH DAMAGE.
 !
-!> Defines an output for a mean flow field
+!> Implements `fluid_stats_ouput_t`.
 module fluid_stats_output
   use fluid_stats, only : fluid_stats_t
   use neko_config, only : NEKO_BCKND_DEVICE
   use num_types, only : rp
+  use map_1d, only : map_1d_t
+  use map_2d, only : map_2d_t
+  use fld_file_data, only : fld_file_data_t
   use device
   use output, only : output_t
+  use matrix, only : matrix_t
   implicit none
   private
 
+  !> Defines an output for the fluid statistics computed using the 
+  !! `fluid_stats_t` object.
   type, public, extends(output_t) :: fluid_stats_output_t
+     !> Pointer to the object computing the statistics.
      type(fluid_stats_t), pointer :: stats
+     !> Space averaging object for 2 homogeneous directions.
+     type(map_1d_t) :: map_1d
+     !> Space averaging object for 1 homogeneous direction.
+     type(map_2d_t) :: map_2d
      real(kind=rp) :: T_begin
+     !> The dimension of the output fields. Either 1, 2, or 3.
+     integer :: output_dim
    contains
+     !> Constructor.
+     procedure, pass(this) :: init => fluid_stats_output_init
+     !> Samples the fields computed by the `stats` component.
      procedure, pass(this) :: sample => fluid_stats_output_sample
   end type fluid_stats_output_t
 
-  interface fluid_stats_output_t
-     module procedure fluid_stats_output_init
-  end interface fluid_stats_output_t
 
 contains
 
-  function fluid_stats_output_init(stats, T_begin, name, path) result(this)
-    type(fluid_stats_t), intent(in), target :: stats
+  !> Constructor.
+  subroutine fluid_stats_output_init(this, stats, T_begin, hom_dir, name, path)
+    type(fluid_stats_t), intent(inout), target :: stats
     real(kind=rp), intent(in) :: T_begin
+    character(len=*), intent(in) :: hom_dir
     character(len=*), intent(in), optional :: name
     character(len=*), intent(in), optional :: path
-    type(fluid_stats_output_t) :: this
+    class(fluid_stats_output_t), intent(inout) :: this
     character(len=1024) :: fname
 
-    if (present(name) .and. present(path)) then
-       fname = trim(path) // trim(name) // '.fld'
-    else if (present(name)) then
-       fname = trim(name) // '.fld'
-    else if (present(path)) then
-       fname = trim(path) // 'stats.fld'
+    if (trim(hom_dir) .eq. 'none' .or. &
+        trim(hom_dir) .eq. 'x' .or.&
+        trim(hom_dir) .eq. 'y' .or.&
+        trim(hom_dir) .eq. 'z'&
+       ) then
+       if (present(name) .and. present(path)) then
+          fname = trim(path) // trim(name) // '.fld'
+       else if (present(name)) then
+          fname = trim(name) // '.fld'
+       else if (present(path)) then
+          fname = trim(path) // 'fluid_stats.fld'
+       else
+          fname = 'fluid_stats.fld'
+       end if
+
+       this%output_dim = 3
+
+       if (trim(hom_dir) .eq. 'x' .or.&
+           trim(hom_dir) .eq. 'y' .or.&
+           trim(hom_dir) .eq. 'z' ) then
+          call this%map_2d%init_char(stats%coef, hom_dir, 1e-7_rp)
+          this%output_dim = 2
+       end if
     else
-       fname = 'stats.fld'
+       if (present(name) .and. present(path)) then
+          fname = trim(path) // trim(name) // '.csv'
+       else if (present(name)) then
+          fname = trim(name) // '.csv'
+       else if (present(path)) then
+          fname = trim(path) // 'fluid_stats.csv'
+       else
+          fname = 'fluid_stats.csv'
+       end if
+       call this%map_1d%init_char(stats%coef, hom_dir, 1e-7_rp)
+       this%output_dim = 1
     end if
 
     call this%init_base(fname)
     this%stats => stats
     this%T_begin = T_begin
-  end function fluid_stats_output_init
+  end subroutine fluid_stats_output_init
 
-  !> Sample a mean flow field at time @a t
+  !> Sample fluid_stats at time @a t
   subroutine fluid_stats_output_sample(this, t)
     class(fluid_stats_output_t), intent(inout) :: this
     real(kind=rp), intent(in) :: t
     integer :: i
+    type(matrix_t) :: avg_output_1d
+    type(fld_file_data_t) :: output_2d
+    real(kind=rp) :: u, v, w, p
     associate (out_fields => this%stats%stat_fields%items)
       if (t .ge. this%T_begin) then
          call this%stats%make_strong_grad()
          if ( NEKO_BCKND_DEVICE .eq. 1) then
             do i = 1, size(out_fields)
                call device_memcpy(out_fields(i)%ptr%x, out_fields(i)%ptr%x_d,&
-                  out_fields(i)%ptr%dof%size(), DEVICE_TO_HOST, &
-                  sync=(i .eq. size(out_fields))) ! Sync on last field
+                    out_fields(i)%ptr%dof%size(), DEVICE_TO_HOST, &
+                    sync = (i .eq. size(out_fields))) ! Sync on last field
+            end do
+         end if
+         if (this%output_dim .eq. 1) then
+            call this%map_1d%average_planes(avg_output_1d, &
+                 this%stats%stat_fields)
+            call this%file_%write(avg_output_1d, t)
+         else if (this%output_dim .eq. 2) then
+            call this%map_2d%average(output_2d, this%stats%stat_fields)
+            !Switch around fields to get correct orders
+            !Put average direction mean_vel in scalar45
+            do i = 1, this%map_2d%n_2d
+               u = output_2d%v%x(i)
+               v = output_2d%w%x(i)
+               w = output_2d%p%x(i)
+               p = output_2d%u%x(i)
+               output_2d%p%x(i) = p
+               output_2d%u%x(i) = u
+               output_2d%v%x(i) = v
+               output_2d%w%x(i) = w
             end do
+            
+            call this%file_%write(output_2d, t)
+         else
+            call this%file_%write(this%stats%stat_fields, t)
          end if
-         call this%file_%write(this%stats%stat_fields, t)
          call this%stats%reset()
       end if
     end associate
diff --git a/src/sem/map_1d.f90 b/src/sem/map_1d.f90
index c84708fcd69..b8323404a0d 100644
--- a/src/sem/map_1d.f90
+++ b/src/sem/map_1d.f90
@@ -209,10 +209,6 @@ subroutine map_1d_init(this, coef,  dir, tol)
     end do
     this%n_el_lvls = glimax(this%el_lvl,nelv)
     this%n_gll_lvls = this%n_el_lvls*lx
-    if ( pe_rank .eq. 0) then
-       write(*,*) 'Number of element levels', this%n_el_lvls
-       write(*,*) 'Total number of levels', this%n_gll_lvls
-    end if
     !Numbers the points in each element based on the element level
     !and its orientation
     do e = 1, nelv
diff --git a/src/simulation.f90 b/src/simulation.f90
index b05f579aa11..df7d065067f 100644
--- a/src/simulation.f90
+++ b/src/simulation.f90
@@ -90,7 +90,6 @@ subroutine neko_solve(C)
 
     !> Call stats, samplers and user-init before time loop
     call neko_log%section('Postprocessing')
-    call C%q%eval(t, C%dt, tstep)
     call C%s%sample(t, tstep)
 
     call C%usr%user_init_modules(t, C%fluid%u, C%fluid%v, C%fluid%w,&
@@ -156,7 +155,6 @@ subroutine neko_solve(C)
        ! Execute all simulation components
        call neko_simcomps%compute(t, tstep)
 
-       call C%q%eval(t, C%dt, tstep)
        call C%s%sample(t, tstep)
 
        ! Update material properties
diff --git a/src/simulation_components/fluid_stats_simcomp.f90 b/src/simulation_components/fluid_stats_simcomp.f90
new file mode 100644
index 00000000000..3951862ce90
--- /dev/null
+++ b/src/simulation_components/fluid_stats_simcomp.f90
@@ -0,0 +1,215 @@
+! Copyright (c) 2024, The Neko Authors
+! All rights reserved.
+!
+! Redistribution and use in source and binary forms, with or without
+! modification, are permitted provided that the following conditions
+! are met:
+!
+!   * Redistributions of source code must retain the above copyright
+!     notice, this list of conditions and the following disclaimer.
+!
+!   * Redistributions in binary form must reproduce the above
+!     copyright notice, this list of conditions and the following
+!     disclaimer in the documentation and/or other materials provided
+!     with the distribution.
+!
+!   * Neither the name of the authors nor the names of its
+!     contributors may be used to endorse or promote products derived
+!     from this software without specific prior written permission.
+!
+! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
+! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+! POSSIBILITY OF SUCH DAMAGE.
+!
+!
+!> Implements the `fluid_stats_simcomp_t` type.
+module fluid_stats_simcomp
+  use num_types, only : rp, dp, sp
+  use json_module, only : json_file
+  use simulation_component, only : simulation_component_t
+  use field_registry, only : neko_field_registry
+  use field, only : field_t
+  use fluid_stats, only: fluid_stats_t
+  use fluid_stats_output, only : fluid_stats_output_t
+  use case, only : case_t
+  use coefs, only : coef_t
+  use comm
+  use logger, only : LOG_SIZE, neko_log
+  use json_utils, only : json_get, json_get_or_default
+  implicit none
+  private
+
+  !> A simulation component that computes the velocity and pressure statistics 
+  !! up to 4th order. Can be used to reconstruct the term budget of transport
+  !! equations for, e.g. the Reynolds stresses and the turbulent kinetic energy.
+  !! 
+  !! Similar in functionality to the satistics module in the KTH Framework for
+  !! Nek5000: https://github.com/KTH-Nek5000/KTH_Framework
+  !! See Turbulence Statistics in a Spectral-Element Code: A Toolbox for 
+  !! High-Fidelity Simulations or the origin KTH Nek5000 framework for details. 
+  !!
+  !! For further details see the Neko documentation.
+  type, public, extends(simulation_component_t) :: fluid_stats_simcomp_t
+     !> Backbone object computing the satistics
+     type(fluid_stats_t) :: stats
+     !> Output writer.
+     type(fluid_stats_output_t) :: stats_output
+     !> Time value at which the sampling of statistics is initiated.
+     real(kind=rp) :: start_time
+     real(kind=rp) :: time
+   contains
+     !> Constructor from json, wrapping the actual constructor.
+     procedure, pass(this) :: init => fluid_stats_simcomp_init_from_json
+     !> Actual constructor.
+     procedure, pass(this) :: init_from_attributes => &
+        fluid_stats_simcomp_init_from_attributes
+     !> Destructor.
+     procedure, pass(this) :: free => fluid_stats_simcomp_free
+     !> Does sampling for statistics.
+     procedure, pass(this) :: compute_ => fluid_stats_simcomp_compute
+     !> Write the statistics to disk.
+     procedure, pass(this) :: output_ => fluid_stats_simcomp_compute
+     !> Restart the simcomp.
+     procedure, pass(this) :: restart_ => fluid_stats_simcomp_restart
+  end type fluid_stats_simcomp_t
+
+contains
+
+  !> Constructor from json.
+  !> @param json JSON object with the parameters.
+  !! @param case The case object.
+  subroutine fluid_stats_simcomp_init_from_json(this, json, case)
+    class(fluid_stats_simcomp_t), intent(inout) :: this
+    type(json_file), intent(inout) :: json
+    class(case_t), intent(inout), target :: case
+    character(len=:), allocatable :: filename
+    character(len=:), allocatable :: precision
+    character(len=20), allocatable :: fields(:)
+    character(len=:), allocatable :: hom_dir
+    character(len=:), allocatable :: stat_set
+    real(kind=rp) :: start_time
+    type(field_t), pointer :: u, v, w, p
+    type(coef_t), pointer :: coef
+
+    call this%init_base(json, case)
+    call json_get_or_default(json, 'avg_direction', &
+         hom_dir, 'none')
+    call json_get_or_default(json, 'start_time', &
+         start_time, 0.0_rp)
+    call json_get_or_default(json, 'set_of_stats', &
+         stat_set, 'full')
+
+    u => neko_field_registry%get_field("u")
+    v => neko_field_registry%get_field("v")
+    w => neko_field_registry%get_field("w")
+    p => neko_field_registry%get_field("p")
+    coef => case%fluid%c_Xh
+    call fluid_stats_simcomp_init_from_attributes(this, u, v, w, p, coef, &
+         start_time, hom_dir, stat_set)
+
+  end subroutine fluid_stats_simcomp_init_from_json
+
+  !> Actual constructor.
+  !! @param u x-velocity
+  !! @param v x-velocity
+  !! @param w x-velocity
+  !! @param coef sem coefs
+  !! @param start_time time to start sampling stats
+  !! @param hom_dir directions to average in
+  !! @param stat_set Set of statistics to compute (basic/full)
+  subroutine fluid_stats_simcomp_init_from_attributes(this, u, v, w, p, coef, &
+       start_time, hom_dir, stat_set)
+    class(fluid_stats_simcomp_t), intent(inout) :: this
+    character(len=*), intent(in) :: hom_dir
+    character(len=*), intent(in) :: stat_set
+    real(kind=rp), intent(in) :: start_time
+    type(field_t), intent(inout) :: u, v, w, p !>Should really be intent in I think
+    type(coef_t), intent(in) :: coef
+    character(len=LOG_SIZE) :: log_buf
+    
+    call neko_log%section('Fluid stats')
+    write(log_buf, '(A,E15.7)') 'Start time: ', start_time
+    call neko_log%message(log_buf)
+    write(log_buf, '(A,A)') 'Set of statistics: ', trim(stat_set)
+    call neko_log%message(log_buf)
+    write(log_buf, '(A,A)') 'Averaging in direction: ', trim(hom_dir)
+    call neko_log%message(log_buf)
+
+
+    call this%stats%init(coef, u, v, w, p, stat_set)
+
+    this%start_time = start_time
+    this%time = start_time
+
+    call this%stats_output%init(this%stats, this%start_time, & 
+         hom_dir = hom_dir, path = this%case%output_directory)
+
+    call this%case%s%add(this%stats_output, &
+         this%output_controller%control_value, &
+         this%output_controller%control_mode)
+ 
+    call neko_log%end_section()
+  
+  end subroutine fluid_stats_simcomp_init_from_attributes
+
+  !> Destructor.
+  subroutine fluid_stats_simcomp_free(this)
+    class(fluid_stats_simcomp_t), intent(inout) :: this
+    call this%free_base()
+    call this%stats%free()
+  end subroutine fluid_stats_simcomp_free
+
+  subroutine fluid_stats_simcomp_restart(this, t)
+    class(fluid_stats_simcomp_t), intent(inout) :: this
+    real(kind=rp), intent(in) :: t
+    if (t .gt. this%time) this%time = t
+  end subroutine fluid_stats_simcomp_restart
+
+  !> fluid_stats, called depending on compute_control and compute_value
+  !! @param t The time value.
+  !! @param tstep The current time-step
+  subroutine fluid_stats_simcomp_compute(this, t, tstep)
+    class(fluid_stats_simcomp_t), intent(inout) :: this
+    real(kind=rp), intent(in) :: t
+    integer, intent(in) :: tstep
+    real(kind=rp) :: delta_t
+    real(kind=rp) :: sample_start_time, sample_time
+    character(len=LOG_SIZE) :: log_buf
+    integer :: ierr
+
+    if (t .ge. this%start_time) then
+       delta_t = t - this%time
+
+       call MPI_Barrier(NEKO_COMM, ierr)
+
+       sample_start_time = MPI_WTIME()
+
+       call this%stats%update(delta_t)
+       call MPI_Barrier(NEKO_COMM, ierr)
+       this%time = t
+
+       sample_time = MPI_WTIME() - sample_start_time
+
+       call neko_log%section('Fluid stats')
+       write(log_buf, '(A,E15.7)') 'Sampling at time:', t
+       call neko_log%message(log_buf)
+       write(log_buf, '(A33,E15.7)') 'Simulationtime since last sample:', &
+            delta_t
+       call neko_log%message(log_buf)
+       write(log_buf, '(A,E15.7)') 'Sampling time (s):', sample_time
+       call neko_log%message(log_buf)
+       call neko_log%end_section()
+    end if
+
+  end subroutine fluid_stats_simcomp_compute
+
+end module fluid_stats_simcomp
diff --git a/src/simulation_components/simulation_component_fctry.f90 b/src/simulation_components/simulation_component_fctry.f90
index 2ff2df1885b..5deef0cad4f 100644
--- a/src/simulation_components/simulation_component_fctry.f90
+++ b/src/simulation_components/simulation_component_fctry.f90
@@ -35,6 +35,7 @@
 submodule (simulation_component) simulation_component_fctry
   use vorticity, only : vorticity_t
   use force_torque, only : force_torque_t
+  use fluid_stats_simcomp, only : fluid_stats_simcomp_t
   use lambda2, only : lambda2_t
   use probes, only : probes_t
   use les_simcomp, only : les_simcomp_t
@@ -44,12 +45,13 @@
   use derivative, only : derivative_t
 
   ! List of all possible types created by the factory routine
-  character(len=20) :: SIMCOMPS_KNOWN_TYPES(6) = [character(len=20) :: &
+  character(len=20) :: SIMCOMPS_KNOWN_TYPES(7) = [character(len=20) :: &
      "vorticity", &
      "lambda2", &
      "probes", &
      "les_model", &
      "field_writer", &
+     "fluid_stats", &
      "force_torque"]
 
 contains
@@ -88,6 +90,8 @@ module subroutine simulation_component_factory(object, json, case)
        allocate(derivative_t::object)
     else if (trim(type_name) .eq. "force_torque") then
        allocate(force_torque_t::object)
+    else if (trim(type_name) .eq. "fluid_stats") then
+       allocate(fluid_stats_simcomp_t::object)
     else
        type_string =  concat_string_array(SIMCOMPS_KNOWN_TYPES, &
             NEW_LINE('A') // "-  ",  .true.)