diff --git a/src/makefile b/src/makefile
index 2a421f8..32c6569 100644
--- a/src/makefile
+++ b/src/makefile
@@ -5,6 +5,7 @@ OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleMeshBoundary.o $(OBJDIR)/module
 			  	$(OBJDIR)/moduleCollisions.o $(OBJDIR)/moduleTable.o $(OBJDIR)/moduleParallel.o \
 			  	$(OBJDIR)/moduleEM.o $(OBJDIR)/moduleRandom.o $(OBJDIR)/moduleMath.o \
 					$(OBJDIR)/moduleProbe.o $(OBJDIR)/moduleAverage.o \
+			  	$(OBJDIR)/moduleMeshInputVTK.o $(OBJDIR)/moduleMeshOutputVTK.o \
 			  	$(OBJDIR)/moduleMeshInputGmsh2.o $(OBJDIR)/moduleMeshOutputGmsh2.o \
 			  	$(OBJDIR)/moduleMeshInput0D.o $(OBJDIR)/moduleMeshOutput0D.o \
 			  	$(OBJDIR)/moduleMesh3DCart.o \
diff --git a/src/modules/makefile b/src/modules/makefile
index 254d015..743144f 100644
--- a/src/modules/makefile
+++ b/src/modules/makefile
@@ -11,7 +11,7 @@ common.o:
 output.o: moduleSpecies.o common.o
 	$(MAKE) -C output all
 
-mesh.o: moduleCollisions.o moduleBoundary.o output.o
+mesh.o: moduleCollisions.o moduleBoundary.o output.o common.o 
 	$(MAKE) -C mesh all
 
 solver.o: moduleSpecies.o moduleProbe.o common.o output.o mesh.o
diff --git a/src/modules/mesh/inout/makefile b/src/modules/mesh/inout/makefile
index 2f73e73..81e15bc 100644
--- a/src/modules/mesh/inout/makefile
+++ b/src/modules/mesh/inout/makefile
@@ -1,4 +1,7 @@
-all: gmsh2.o 0D.o
+all: vtk.o gmsh2.o 0D.o
+
+vtk.o:
+	$(MAKE) -C vtk all
 
 gmsh2.o:
 	$(MAKE) -C gmsh2 all
diff --git a/src/modules/mesh/inout/vtk/makefile b/src/modules/mesh/inout/vtk/makefile
new file mode 100644
index 0000000..18a13a0
--- /dev/null
+++ b/src/modules/mesh/inout/vtk/makefile
@@ -0,0 +1,7 @@
+all: moduleMeshInputVTK.o moduleMeshOutputVTK.o
+
+moduleMeshInputVTK.o: moduleMeshOutputVTK.o moduleMeshInputVTK.f90
+	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
+
+%.o: %.f90
+	$(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@
diff --git a/src/modules/mesh/inout/vtk/moduleMeshInputVTK.f90 b/src/modules/mesh/inout/vtk/moduleMeshInputVTK.f90
new file mode 100644
index 0000000..a044b4f
--- /dev/null
+++ b/src/modules/mesh/inout/vtk/moduleMeshInputVTK.f90
@@ -0,0 +1,3 @@
+MODULE moduleMeshInputVTK
+
+END MODULE moduleMeshInputVTK
diff --git a/src/modules/mesh/inout/vtk/moduleMeshOutputVTK.f90 b/src/modules/mesh/inout/vtk/moduleMeshOutputVTK.f90
new file mode 100644
index 0000000..efaf6e3
--- /dev/null
+++ b/src/modules/mesh/inout/vtk/moduleMeshOutputVTK.f90
@@ -0,0 +1,224 @@
+MODULE moduleMeshOutputVTK
+
+  CONTAINS
+
+    SUBROUTINE writeFileHeader(self, fileID)
+      USE moduleMesh
+      IMPLICIT NONE
+
+      CLASS(meshParticles), INTENT(in):: self
+      INTEGER, INTENT(in):: fileID
+
+      WRITE(fileID,"(A)") '<?xml version="1.0"?>'
+      WRITE(fileID,"(2X, A)") '<VTKFile type="UnstructuredGrid">'
+      WRITE(fileID,"(4X, A,ES20.6E3,A)") '<UnstructuredGrid>'
+      WRITE(fileID,"(6X, A, I10, A, I10, A)") '<Piece NumberOfPoints="', self%numNodes, '" NumberOfCells="', self%numCells, '">'
+
+    END SUBROUTINE writeFileHeader
+
+    SUBROUTINE writeFileFooter(fileID)
+      IMPLICIT NONE
+
+      INTEGER, INTENT(in):: fileID
+
+      WRITE(fileID,"(6X, A)") '</Piece>'
+      WRITE(fileID,"(4X, A)") '</UnstructuredGrid>'
+      WRITE(fileID,"(2X, A)") '</VTKFile>'
+
+    END SUBROUTINE writeFileFooter
+
+    FUNCTION getCellType(cell) RESULT(indexType)
+      USE moduleMesh3DCart
+      USE moduleMesh2DCyl
+      USE moduleMesh2DCart
+      USE moduleMesh1DRad
+      USE moduleMesh1DCart
+      USE moduleMesh0D
+      USE moduleErrors
+      IMPLICIT NONE
+
+      CLASS(meshCell), INTENT(in):: cell
+      INTEGER:: indexType
+      
+      indexType = 0
+
+      SELECT TYPE(cell)
+      TYPE IS(meshCell3DCartTetra)
+        indexType = 10
+
+      TYPE IS(meshCell2DCylQuad)
+        indexType = 9
+
+      TYPE IS(meshCell2DCartQuad)
+        indexType = 9
+
+      TYPE IS(meshCell2DCylTria)
+        indexType = 5
+
+      TYPE IS(meshCell2DCartTria)
+        indexType = 5
+
+      TYPE IS(meshCell1DRadSegm)
+        indexType = 3
+
+      TYPE IS(meshCell1DCartSegm)
+        indexType = 3
+
+      TYPE IS(meshCell0D)
+        indexType = 1
+
+      CLASS DEFAULT
+        CALL criticalError('Cell not valid for VTK output', 'getCellType')
+
+      END SELECT
+
+    END FUNCTION getCellType
+
+    SUBROUTINE writeFileMesh(self, fileID)
+      USE moduleMesh
+      USE moduleRefParam
+      IMPLICIT NONE
+
+      CLASS(meshParticles), INTENT(in):: self
+      INTEGER, INTENT(in):: fileID
+      CHARACTER(LEN=25):: nodeFormat
+      CHARACTER(LEN=25):: cellFormat
+      INTEGER:: e
+      INTEGER:: offset
+
+      !Write nodes coordinates
+      WRITE(fileID, "(8X, A)") '<Points>'
+      WRITE(fileID, "(10X,A)") '<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
+      WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(3(ES20.6E3)))"
+      WRITE(fileID, nodeFormat) (self%nodes(e)%obj%getCoordinates()*L_ref, e = 1, self%numNodes)
+      WRITE(fileID, "(10X, A)") '</DataArray>'
+      WRITE(fileID, "(8X, A)") '</Points>'
+
+      WRITE(fileID, "(8X, A)") '<Cells>'
+      !Write nodes connectivity of each cell
+      WRITE(fileID, "(10X,A)") '<DataArray type="Int32" Name="connectivity"  format="ascii">'
+      DO e = 1, self%numCells
+        WRITE(cellFormat, "(A,I1,A)") "(",self%cells(e)%obj%nNodes,"(I10))"
+        WRITE(fileID, cellFormat, advance="no") self%cells(e)%obj%getNodes(self%cells(e)%obj%nNodes) - 1 !Array starts on 0
+
+      END DO
+      WRITE(fileID, "(10X, A)") '</DataArray>'
+      !Write offset of each cell
+      offset = 0
+      WRITE(fileID, "(10X,A)") '<DataArray type="Int32" Name="offsets"  format="ascii">'
+      DO e = 1, self%numCells
+        WRITE(cellFormat, "(A,I1,A)") "(I10)"
+        offset = offset + self%cells(e)%obj%nNodes
+        WRITE(fileID, cellFormat, advance="no") offset
+
+      END DO
+      WRITE(fileID, "(10X, A)") '</DataArray>'
+      !Write type of each cell
+      WRITE(fileID, "(10X,A)") '<DataArray type="Int32" Name="types"  format="ascii">'
+      DO e = 1, self%numCells
+        WRITE(cellFormat, "(A,I1,A)") "(I10)"
+        WRITE(fileID, cellFormat, advance="no") getCellType(self%cells(e)%obj)
+
+      END DO
+      WRITE(fileID, "(10X, A)") '</DataArray>'
+      WRITE(fileID, "(8X, A)") '</Cells>'
+
+    END SUBROUTINE writeFileMesh
+
+    SUBROUTINE writeNodeData(self, fileID, output)
+      USE moduleMesh
+      USE moduleOutput
+      IMPLICIT NONE
+
+      CLASS(meshParticles), INTENT(in):: self
+      INTEGER, INTENT(in):: fileID
+      TYPE(outputFormat):: output(1:self%numNodes)
+      CHARACTER(LEN=25):: nodeFormat
+      INTEGER:: n
+
+      WRITE(fileID,"(A)") '<PointData>'
+      WRITE(fileID,"(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,"(A)") '</DataArray>'
+      WRITE(fileID,"(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,"(A)") '</DataArray>'
+      WRITE(fileID,"(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,"(A)") '</DataArray>'
+      WRITE(fileID,"(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,"(A)") '</DataArray>'
+      WRITE(fileID,"(A)") '</PointData>'
+
+    END SUBROUTINE writeNodeData
+    
+    SUBROUTINE printOutputVTK(self,t)
+      USE moduleMesh
+      USE moduleRefParam
+      USE moduleSpecies
+      USE moduleOutput
+      USE moduleCaseParam
+      IMPLICIT NONE
+
+      CLASS(meshParticles), INTENT(in):: self
+      INTEGER, INTENT(in):: t
+      INTEGER:: n, i, fileID
+      CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection
+      CHARACTER (LEN=iterationDigits):: tstring
+      TYPE(outputFormat):: output(1:self%numNodes)
+
+      fileID = 60
+
+      DO i = 1, nSpecies
+        WRITE(tstring, iterationFormat) t
+        fileName= 'OUTPUT_' // tstring// '_' // species(i)%obj%name // '.vtu'
+
+        WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
+        OPEN (fileID, file = path // folder // '/' //  fileName)
+
+        fileNameCollection = 'OUTPUT_Collection_' // species(i)%obj%name // '.vtu'
+        IF (t == tInitial) THEN
+          !Create collection file
+          WRITE(*, "(6X,A15,A)") "Creating file: ", fileNameCollection
+          OPEN (fileID + 1, file = path // folder  // '/' // fileNameCollection)
+          WRITE (fileID + 1, "(A)") '<VTKFile type="Collection">'
+          WRITE (fileID + 1, "(2X, A)") '<Collection>'
+          CLOSE(fileID + 1)
+
+        END IF
+
+        OPEN (fileID + 1, file = path // folder  // '/' // fileNameCollection, ACCESS='APPEND')
+        WRITE(fileID + 1, "(4X, A, ES20.6E3, A, A, A)"), '<DataSet timestep="', t*ti_ref,'" file="', fileName,'"/>'
+
+        IF (t == tFinal) THEN
+          WRITE (fileID + 1, "(2X, A)") '</Collection>'
+          WRITE (fileID + 1, "(A)") '</VTKFile>'
+
+        END IF
+        CLOSE(fileID + 1)
+
+        CALL writeFileHeader(self, fileID)
+
+        CALL writeFileMesh(self, fileID)
+
+        DO n = 1, self%numNodes
+          CALL calculateOutput(self%nodes(n)%obj%output(i), output(n), self%nodes(n)%obj%v, species(i)%obj)
+
+        END DO
+
+        CALL writeNodeData(self, fileID, output)
+
+        CALL writeFileFooter(fileID)
+
+        CLOSE(fileID)
+
+      END DO
+
+    END SUBROUTINE printOutputVTK
+
+END MODULE moduleMeshOutputVTK
diff --git a/src/modules/solver/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90
index e557495..813aeb3 100644
--- a/src/modules/solver/moduleSolver.f90
+++ b/src/modules/solver/moduleSolver.f90
@@ -514,6 +514,7 @@ MODULE moduleSolver
       USE moduleSpecies
       USE moduleCompTime
       USE moduleProbe
+      USE moduleMeshOutputVTK !TEMPORARY TO TEST VTK OUTPUT
       IMPLICIT NONE
 
       INTEGER, INTENT(in):: t
@@ -527,6 +528,7 @@ MODULE moduleSolver
 
         CALL outputProbes(t)
         CALL mesh%printOutput(t)
+        CALL printOutputVTK(mesh, t)
         IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%printColl(t)
         CALL mesh%printEM(t)
         WRITE(*, "(5X,A21,I10,A1,I10)") "t/tFinal: ", t, "/", tFinal