diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90
index 0a3a557..7832843 100644
--- a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90
+++ b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90
@@ -13,9 +13,10 @@ MODULE moduleMeshInputGmsh2
       IF (ASSOCIATED(meshForMCC, self)) self%printColl => printCollGmsh2
       SELECT TYPE(self)
       TYPE IS(meshParticles)
-        self%printOutput => printOutputGmsh2
-        self%printEM     => printEMGmsh2
-        self%readInitial => readInitialGmsh2
+        self%printOutput  => printOutputGmsh2
+        self%printEM      => printEMGmsh2
+        self%readInitial  => readInitialGmsh2
+        self%printAverage => printAverageGmsh2
 
       END SELECT
       self%readMesh    => readGmsh2
diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90
index 40eda4f..106480a 100644
--- a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90
+++ b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90
@@ -21,7 +21,7 @@ MODULE moduleMeshOutputGmsh2
 
       DO i = 1, nSpecies
         WRITE(tstring, iterationFormat) t
-        fileName='OUTPUT_' // tstring// '_' // species(i)%obj%name // '.msh'
+        fileName= 'OUTPUT_' // tstring// '_' // species(i)%obj%name // '.msh'
         WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
         OPEN (60, file = path // folder // '/' //  fileName)
         WRITE(60, "(A)") '$MeshFormat'
@@ -226,4 +226,84 @@ MODULE moduleMeshOutputGmsh2
 
     END SUBROUTINE printEMGmsh2
 
+    !Prints the average properties of particles into the nodes
+    SUBROUTINE printAverageGmsh2(self)
+      USE moduleMesh
+      USE moduleRefParam
+      USE moduleSpecies
+      USE moduleOutput
+      USE moduleAverage
+      IMPLICIT NONE
+
+      CLASS(meshParticles), INTENT(in):: self
+      INTEGER:: n, i
+      TYPE(outputFormat):: output(1:self%numNodes)
+      CHARACTER(:), ALLOCATABLE:: fileName
+
+      DO i = 1, nSpecies
+        fileName= 'Average_mean_OUTPUT_' // species(i)%obj%name // '.msh'
+        WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
+        OPEN (60, file = path // folder // '/' //  fileName)
+        WRITE(60, "(A)") '$MeshFormat'
+        WRITE(60, "(A)") '2.2 0 8'
+        WRITE(60, "(A)") '$EndMeshFormat'
+        WRITE(60, "(A)") '$NodeData'
+        WRITE(60, "(A)") '1'
+        WRITE(60, "(A)") '"' // species(i)%obj%name // ' density (m^-3)"'
+        WRITE(60, *) 1
+        WRITE(60, *) 0.D0
+        WRITE(60, *) 3
+        WRITE(60, *) 0
+        WRITE(60, *) 1
+        WRITE(60, *) self%numNodes
+        DO n=1, self%numNodes
+          CALL calculateOutput(averageScheme%mean(n)%output(i), output(n), self%nodes(n)%obj%v, species(i)%obj)
+          WRITE(60, "(I6,ES20.6E3)") n, output(n)%density
+        END DO
+        WRITE(60, "(A)") '$EndNodeData'
+        WRITE(60, "(A)") '$NodeData'
+        WRITE(60, "(A)") '1'
+        WRITE(60, "(A)") '"' // species(i)%obj%name // ' velocity (m s^-1)"'
+        WRITE(60, *) 1
+        WRITE(60, *) 0.D0
+        WRITE(60, *) 3
+        WRITE(60, *) 0
+        WRITE(60, *) 3
+        WRITE(60, *) self%numNodes
+        DO n=1, self%numNodes
+          WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%velocity
+        END DO
+        WRITE(60, "(A)") '$EndNodeData'
+        WRITE(60, "(A)") '$NodeData'
+        WRITE(60, "(A)") '1'
+        WRITE(60, "(A)") '"' // species(i)%obj%name // ' pressure (Pa)"'
+        WRITE(60, *) 1
+        WRITE(60, *) 0.D0
+        WRITE(60, *) 3
+        WRITE(60, *) 0
+        WRITE(60, *) 1
+        WRITE(60, *) self%numNodes
+        DO n=1, self%numNodes
+          WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%pressure
+        END DO
+        WRITE(60, "(A)") '$EndNodeData'
+        WRITE(60, "(A)") '$NodeData'
+        WRITE(60, "(A)") '1'
+        WRITE(60, "(A)") '"' // species(i)%obj%name // ' temperature (K)"'
+        WRITE(60, *) 1
+        WRITE(60, *) 0.D0
+        WRITE(60, *) 3
+        WRITE(60, *) 0
+        WRITE(60, *) 1
+        WRITE(60, *) self%numNodes
+        DO n=1, self%numNodes
+          WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%temperature
+        END DO
+        WRITE(60, "(A)") '$EndNodeData'
+        CLOSE (60)
+
+      END DO
+
+    END SUBROUTINE printAverageGmsh2
+
 END MODULE moduleMeshOutputGmsh2
diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90
index ebbbed4..05c4b85 100644
--- a/src/modules/mesh/moduleMesh.f90
+++ b/src/modules/mesh/moduleMesh.f90
@@ -289,7 +289,7 @@ MODULE moduleMesh
     CONTAINS
       PROCEDURE, PASS:: doCollisions
 
-  END TYPE
+  END TYPE meshGeneric
 
   ABSTRACT INTERFACE
     !Reads the mesh from a file
@@ -338,9 +338,10 @@ MODULE moduleMesh
     REAL(8), ALLOCATABLE, DIMENSION(:,:):: K 
     !Permutation matrix for P L U factorization
     INTEGER, ALLOCATABLE, DIMENSION(:,:):: IPIV
-    PROCEDURE(printOutput_interface), POINTER, PASS:: printOutput => NULL()
-    PROCEDURE(printEM_interface),     POINTER, PASS:: printEM     => NULL()
-    PROCEDURE(doCoulomb_interface),   POINTER, PASS:: doCoulomb   => NULL()
+    PROCEDURE(printOutput_interface),  POINTER, PASS:: printOutput  => NULL()
+    PROCEDURE(printEM_interface),      POINTER, PASS:: printEM      => NULL()
+    PROCEDURE(doCoulomb_interface),    POINTER, PASS:: doCoulomb    => NULL()
+    PROCEDURE(printAverage_interface), POINTER, PASS:: printAverage => NULL()
     CONTAINS
       PROCEDURE, PASS:: constructGlobalK
 
@@ -373,6 +374,14 @@ MODULE moduleMesh
 
     END SUBROUTINE printEM_interface
 
+    !Prints average values
+    SUBROUTINE printAverage_interface(self)
+      IMPORT meshParticles
+
+      CLASS(meshParticles), INTENT(in):: self
+
+    END SUBROUTINE printAverage_interface
+
 
   END INTERFACE
 
diff --git a/src/modules/moduleAverage.f90 b/src/modules/moduleAverage.f90
index 9b614c0..b697389 100644
--- a/src/modules/moduleAverage.f90
+++ b/src/modules/moduleAverage.f90
@@ -9,7 +9,7 @@ MODULE moduleAverage
 
   !Generic type for average scheme
   TYPE, PUBLIC:: averageGeneric
-    INTEGER:: tStart !Starting iteartion for average scheme
+    INTEGER:: tStart = 0 !Starting iteartion for average scheme
     TYPE(averageData), ALLOCATABLE:: mean(:)
     TYPE(averageData), ALLOCATABLE:: deviation(:)
     CONTAINS
@@ -20,7 +20,7 @@ MODULE moduleAverage
   TYPE(averageGeneric):: averageScheme
 
   !Logical to determine if average scheme must be used
-  LOGICAL:: useAverage
+  LOGICAL:: useAverage = .FALSE.
 
   CONTAINS
     !Based on Welford's online algorithm
diff --git a/src/modules/moduleSolver.f90 b/src/modules/moduleSolver.f90
index 7182184..4d0ec0d 100644
--- a/src/modules/moduleSolver.f90
+++ b/src/modules/moduleSolver.f90
@@ -818,6 +818,12 @@ MODULE moduleSolver
 
       END IF
 
+      !Output average values
+      IF (useAverage .AND. t == tFinal) THEN
+        CALL mesh%printAverage()
+
+      END IF
+
     END SUBROUTINE doOutput
 
     SUBROUTINE doAverage(t)