diff --git a/src/modules/mesh/inout/vtk/moduleMeshOutputVTK.f90 b/src/modules/mesh/inout/vtk/moduleMeshOutputVTK.f90
index 203de0c..d962ae9 100644
--- a/src/modules/mesh/inout/vtk/moduleMeshOutputVTK.f90
+++ b/src/modules/mesh/inout/vtk/moduleMeshOutputVTK.f90
@@ -164,23 +164,29 @@ MODULE moduleMeshOutputVTK
 
       END DO
 
+      !Write node data
       WRITE(fileID,"(8X,A)") '<PointData>'
+      !Write density
       WRITE(fileID,"(10X,A)") '<DataArray type="Float64" Name="Density (m^-3)" NumberOfComponents="1">'
       WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))"
       WRITE(fileID,nodeFormat) (output(n)%density, n  = 1, self%numNodes)
       WRITE(fileID,"(10X,A)") '</DataArray>'
+      !Write velocity
       WRITE(fileID,"(10X,A)") '<DataArray type="Float64" Name="Velocity (m s^-1)" NumberOfComponents="3">'
       WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(3(ES20.6E3)))"
-      WRITE(fileID,nodeFormat) (output(n)%velocity(1:3), n  = 1, self%numNodes)
+      WRITE(fileID,nodeFormat) (output(n)%velocity, n  = 1, self%numNodes)
       WRITE(fileID,"(10X,A)") '</DataArray>'
+      !Write pressure
       WRITE(fileID,"(10X,A)") '<DataArray type="Float64" Name="Pressure (Pa)" NumberOfComponents="1">'
       WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))"
       WRITE(fileID,nodeFormat) (output(n)%pressure, n  = 1, self%numNodes)
       WRITE(fileID,"(10X,A)") '</DataArray>'
+      !Write temperature
       WRITE(fileID,"(10X,A)") '<DataArray type="Float64" Name="Temperature (K)" NumberOfComponents="1">'
       WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))"
       WRITE(fileID,nodeFormat) (output(n)%temperature, n  = 1, self%numNodes)
       WRITE(fileID,"(10X,A)") '</DataArray>'
+      !End node data
       WRITE(fileID,"(8X,A)") '</PointData>'
 
     END SUBROUTINE writeSpeciesOutput
@@ -286,6 +292,69 @@ MODULE moduleMeshOutputVTK
 
     END SUBROUTINE writeCollection
     
+    SUBROUTINE writeAverage(self, fileIDMean, fileIDDeviation, speciesIndex)
+      USE moduleMesh
+      USE moduleOutput
+      USE moduleAverage
+      IMPLICIT NONE
+
+      CLASS(meshParticles), INTENT(in):: self
+      INTEGER, INTENT(in):: fileIDMean, fileIDDeviation
+      INTEGER, INTENT(in):: speciesIndex
+      TYPE(outputFormat):: outputMean(1:self%numNodes)
+      TYPE(outputFormat):: outputDeviation(1:self%numNodes)
+      CHARACTER(LEN=25):: nodeFormat
+      INTEGER:: n
+
+      DO n = 1, self%numNodes
+        CALL calculateOutput(averageScheme(n)%mean%output(speciesIndex),      outputMean(n),      &
+                             self%nodes(n)%obj%v, species(speciesIndex)%obj)
+        CALL calculateOutput(averageScheme(n)%deviation%output(speciesIndex), outputDeviation(n), &
+                             self%nodes(n)%obj%v, species(speciesIndex)%obj)
+
+      END DO
+
+      !Write node data
+      WRITE(fileIDMean,     "(8X,A)") '<PointData>'
+      WRITE(fileIDDeviation,"(8X,A)") '<PointData>'
+      !Write density
+      WRITE(fileIDMean,     "(10X,A)") '<DataArray type="Float64" Name="Density, mean (m^-3)"      NumberOfComponents="1">'
+      WRITE(fileIDDeviation,"(10X,A)") '<DataArray type="Float64" Name="Density, deviation (m^-3)" NumberOfComponents="1">'
+      WRITE(nodeFormat,"(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))"
+      WRITE(fileIDMean,     nodeFormat) (outputMean(n)%density,      n  = 1, self%numNodes)
+      WRITE(fileIDDeviation,nodeFormat) (outputDeviation(n)%density, n  = 1, self%numNodes)
+      WRITE(fileIDMean,"(10X,A)")      '</DataArray>'
+      WRITE(fileIDDeviation,"(10X,A)") '</DataArray>'
+      !Write velocity
+      WRITE(fileIDMean,     "(10X,A)") '<DataArray type="Float64" Name="Velocity, mean (m s^-1)"      NumberOfComponents="3">'
+      WRITE(fileIDDeviation,"(10X,A)") '<DataArray type="Float64" Name="Velocity, deviation (m s^-1)" NumberOfComponents="3">'
+      WRITE(nodeFormat,"(A,I10, A)") "(", self%numNodes, "(3(ES20.6E3)))"
+      WRITE(fileIDMean,     nodeFormat) (outputMean(n)%velocity,      n  = 1, self%numNodes)
+      WRITE(fileIDDeviation,nodeFormat) (outputDeviation(n)%velocity, n  = 1, self%numNodes)
+      WRITE(fileIDMean,     "(10X,A)") '</DataArray>'
+      WRITE(fileIDDeviation,"(10X,A)") '</DataArray>'
+      !Write pressure
+      WRITE(fileIDMean,     "(10X,A)") '<DataArray type="Float64" Name="Pressure, mean (Pa)"      NumberOfComponents="1">'
+      WRITE(fileIDDeviation,"(10X,A)") '<DataArray type="Float64" Name="Pressure, deviation (Pa)" NumberOfComponents="1">'
+      WRITE(nodeFormat,"(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))"
+      WRITE(fileIDMean,     nodeFormat) (outputMean(n)%pressure,      n  = 1, self%numNodes)
+      WRITE(fileIDDeviation,nodeFormat) (outputDeviation(n)%pressure, n  = 1, self%numNodes)
+      WRITE(fileIDMean,     "(10X,A)") '</DataArray>'
+      WRITE(fileIDDeviation,"(10X,A)") '</DataArray>'
+      !Write temperature
+      WRITE(fileIDMean,     "(10X,A)") '<DataArray type="Float64" Name="Temperature, mean (K)"      NumberOfComponents="1">'
+      WRITE(fileIDDeviation,"(10X,A)") '<DataArray type="Float64" Name="Temperature, deviation (K)" NumberOfComponents="1">'
+      WRITE(nodeFormat,"(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))"
+      WRITE(fileIDMean,     nodeFormat) (outputMean(n)%temperature,      n  = 1, self%numNodes)
+      WRITE(fileIDDeviation,nodeFormat) (outputDeviation(n)%temperature, n  = 1, self%numNodes)
+      WRITE(fileIDMean,     "(10X,A)") '</DataArray>'
+      WRITE(fileIDDeviation,"(10X,A)") '</DataArray>'
+      !End node data
+      WRITE(fileIDMean,     "(8X,A)") '</PointData>'
+      WRITE(fileIDDeviation,"(8X,A)") '</PointData>'
+
+    END SUBROUTINE writeAverage
+
     SUBROUTINE printOutputVTK(self,t)
       USE moduleMesh
       USE moduleSpecies
@@ -393,4 +462,44 @@ MODULE moduleMeshOutputVTK
 
     END SUBROUTINE printEMVTK
 
+    SUBROUTINE printAverageVTK(self)
+      USE moduleMesh
+      USE moduleSpecies
+      IMPLICIT NONE
+
+      CLASS(meshParticles), INTENT(in):: self
+      INTEGER:: n, i, fileIDMean, fileIDDeviation
+      CHARACTER(:), ALLOCATABLE:: fileNameMean, fileNameDeviation
+      TYPE(outputFormat):: output(1:self%numNodes)
+
+      fileIDMean = 66
+      fileIDDeviation = 67
+
+      DO i = 1, nSpecies
+        fileNameMean = formatFileName('Average_mean', species(i)%obj%name, 'vtu')
+        WRITE(*, "(6X,A15,A)") "Creating file: ", fileNameMean
+        OPEN (fileIDMean, file = path // folder // '/' //  fileNameMean)
+
+        fileNameDeviation = formatFileName('Average_deviation', species(i)%obj%name, 'vtu')
+        WRITE(*, "(6X,A15,A)") "Creating file: ", fileNameDeviation
+        OPEN (fileIDDeviation, file = path // folder // '/' //  fileNameDeviation)
+
+        CALL writeHeader(self%numNodes, self%numCells, fileIDMean)
+        CALL writeHeader(self%numNodes, self%numCells, fileIDDeviation)
+
+        CALL writeMesh(self, fileIDMean)
+        CALL writeMesh(self, fileIDDeviation)
+
+        CALL writeAverage(self, fileIDMean, fileIDDeviation, i)
+
+        CALL writeFooter(fileIDMean)
+        CALL writeFooter(fileIDDeviation)
+
+        CLOSE(fileIDMean)
+        CLOSE(fileIDDeviation)
+
+      END DO
+
+    END SUBROUTINE printAverageVTK
+
 END MODULE moduleMeshOutputVTK
diff --git a/src/modules/solver/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90
index def5ca4..68c1270 100644
--- a/src/modules/solver/moduleSolver.f90
+++ b/src/modules/solver/moduleSolver.f90
@@ -565,6 +565,7 @@ MODULE moduleSolver
       !Output average values
       IF (useAverage .AND. t == tFinal) THEN
         CALL mesh%printAverage()
+        CALL printAverageVTK(mesh) !TEMPORARY TO TEST VTK OUTPUT
 
       END IF