diff --git a/doc/user-manual/fpakc_UserManual.pdf b/doc/user-manual/fpakc_UserManual.pdf
index 8e8309b..b831bba 100644
Binary files a/doc/user-manual/fpakc_UserManual.pdf and b/doc/user-manual/fpakc_UserManual.pdf differ
diff --git a/doc/user-manual/fpakc_UserManual.tex b/doc/user-manual/fpakc_UserManual.tex
index 3249d02..ae5455f 100644
--- a/doc/user-manual/fpakc_UserManual.tex
+++ b/doc/user-manual/fpakc_UserManual.tex
@@ -10,6 +10,7 @@
 \usepackage[block=ragged,backend=bibtex]{biblatex}
 \usepackage[acronym,toc,automake]{glossaries}
 \usepackage[hidelinks]{hyperref}
+\usepackage[version=4]{mhchem}
 \hypersetup{
   breaklinks = true,                                               % Allows break links in lines
   colorlinks = true,                                               % Colours links instead of ugly boxes
@@ -373,14 +374,22 @@ make
                                \begin{itemize}
                                  \item \textbf{3DCart}: Three-dimensional grid ($x \hyphen y \hyphen z$) in Cartesian coordinates..
                                                         For \Gls{gmsh} mesh format, the coordinates $x$, $y$ and $z$ correspond to $x$, $y$ and $z$ respectively.
-                                 \item \textbf{2DCyl}: Two-dimensional grid ($z \hyphen r$) with symmetry axis at $r = 0$.
-                                                       For \Gls{gmsh} mesh format, the coordinates $x$ and $y$ correspond to  $z$ and $r$ respectively.
+                                 \item \textbf{2DCyl}:  Two-dimensional grid ($z \hyphen r$) with symmetry axis at $r = 0$.
+                                                        For \Gls{gmsh} mesh format, the coordinates $x$ and $y$ correspond to  $z$ and $r$ respectively.
                                  \item \textbf{2DCart}: Two-dimensional grid ($x \hyphen y$) in Cartesian coordinates..
                                                         For \Gls{gmsh} mesh format, the coordinates $x$ and $y$ correspond to  $x$ and $y$ respectively.
-                                 \item \textbf{1DRad}: One-dimensional grid ($r$) in radial coordinates
-                                                       For \Gls{gmsh} mesh format, the coordinates $x$ corresponds to $r$.
+                                 \item \textbf{1DRad}:  One-dimensional grid ($r$) in radial coordinates
+                                                        For \Gls{gmsh} mesh format, the coordinates $x$ corresponds to $r$.
                                  \item \textbf{1DCart}: One-dimensional grid ($x$) in Cartesian coordinates
                                                         For \Gls{gmsh} mesh format, the coordinates $x$ corresponds to $x$.
+                                 \item \textbf{0D}:     Zero dimension ficticius volume.
+                                                        Geometry used mostly to test collisional effects.
+                                                        No boundary or EM field is solved.
+                                                        No injection can be implemented.
+                                                        Initial state must be read from file.
+                                                        No mesh file is required.
+                                                        The optional argument \textbf{geometry.volume} can be used to set a ficticius volume.
+                                                        Otherwise, the volume is set to 1 in non-dimensional units.
                                \end{itemize}
           \item \textbf{meshType}: Character.
                                    Format of mesh file.
@@ -391,6 +400,10 @@ make
           \item \textbf{meshFile}: Character.
                                    Mesh filename.
                                    This file is searched in the path \textbf{output.path} and must contain the file extension.
+          \item \textbf{volume}:   Real
+                                   Units of $\unit{m^-3}$.
+                                   Used to set a ficticius volume for the \textbf{0D} geometry.
+                                   Ignored in the other cases.
         \end{itemize}
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -589,14 +602,16 @@ make
                                  \begin{itemize}
                                    \item \textbf{3DCartNeutral}: Pushes particles in a 3D Cartesian space ($x \hyphen y \hyphen z$) without any external force.
                                    \item \textbf{3DCartCharged}: Pushes particles in a 3D Cartesian space ($x \hyphen y \hyphen z$) including the effect of the electrostatic field.
-                                   \item \textbf{2DCylNeutral}: Pushes particles in a 2D cylindrical space ($z \hyphen r$) without any external force.
-                                   \item \textbf{2DCylCharged}: Pushes particles in a 2D cylindrical space ($z \hyphen r$) including the effect of the electrostatic field.
+                                   \item \textbf{2DCylNeutral}:  Pushes particles in a 2D cylindrical space ($z \hyphen r$) without any external force.
+                                   \item \textbf{2DCylCharged}:  Pushes particles in a 2D cylindrical space ($z \hyphen r$) including the effect of the electrostatic field.
                                    \item \textbf{2DCartNeutral}: Pushes particles in a 2D Cartesian space ($x \hyphen y$) without any external force.
                                    \item \textbf{2DCartCharged}: Pushes particles in a 2D Cartesian space ($x \hyphen y$) including the effect of the electrostatic field.
-                                   \item \textbf{1DRadNeutral}: Pushes particles in a 1D cylindrical space ($r$) without any external force.
-                                   \item \textbf{1DRadCharged}: Pushes particles in a 1D cylindrical space ($r$) accounting the the electrostatic field.
+                                   \item \textbf{1DRadNeutral}:  Pushes particles in a 1D cylindrical space ($r$) without any external force.
+                                   \item \textbf{1DRadCharged}:  Pushes particles in a 1D cylindrical space ($r$) accounting the the electrostatic field.
                                    \item \textbf{1DCartNeutral}: Pushes particles in a 1D Cartesian space ($x$) without any external force.
                                    \item \textbf{1DCartCharged}: Pushes particles in a 1D Cartesian space ($x$) accounting the the electrostatic field.
+                                   \item \textbf{0D}:            Dummy pusher for 0D geometry.
+                                                                 No pushing is actually done.
                                  \end{itemize}
           \item \textbf{WeightingScheme}: Character.
                                           Indicates the variable weighting scheme to be used in the simulation.
@@ -683,11 +698,19 @@ make
         \end{itemize}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   \chapter{Example runs\label{ch:exampleRuns}}
+    This chapter presents a description of the different example files distributed with \acrshort{fpakc}.
+    
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     \section{1D Emissive Cathode (1D\_Cathode)}
       Emission from a 1D cathond in both, cartesian and radial coordinates.
       Both cases insert the same amount of electrons from the minimum coordinate and have the same boundary conditions for particles and the electrostatic field.
       This case is useful to ilustrate hoy \acrshort{fpakc} can deal with different geometries by just modifiying some parameters in the input file.
       The same mesh file (\lstinline|mesh.msh|) is used for both cases.
+      
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    \section{0D \ce{Ar}-\ce{Ar+} Elastic Collision (0D\_Argon)}
+      Test case to check the 0D geometry that includes the elastic collision between \ce{Ar} and \ce{Ar+}.
+      Initial states are readed from the Argon\_Initial.dat and Argon+\_Initial.dat
     
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     \section{ALPHIE Grid system (ALPHIE\_Grid)}
diff --git a/runs/0D_Argon/Argon+_Initial.dat b/runs/0D_Argon/Argon+_Initial.dat
new file mode 100644
index 0000000..75375d1
--- /dev/null
+++ b/runs/0D_Argon/Argon+_Initial.dat
@@ -0,0 +1,2 @@
+ #        t             density                                                  velocity            pressure         temperature
+          0             1.0E+16                                      0.0E+0 0.0E+0 0.0E+0                   0              3000
diff --git a/runs/0D_Argon/Argon_Initial.dat b/runs/0D_Argon/Argon_Initial.dat
new file mode 100644
index 0000000..8aa604f
--- /dev/null
+++ b/runs/0D_Argon/Argon_Initial.dat
@@ -0,0 +1,2 @@
+ #        t             density                                                  velocity            pressure         temperature
+          0             1.0E+16                                      0.0E+0 0.0E+0 0.0E+0                   0               300
diff --git a/runs/0D_Argon/input.json b/runs/0D_Argon/input.json
new file mode 100644
index 0000000..8e62d2a
--- /dev/null
+++ b/runs/0D_Argon/input.json
@@ -0,0 +1,54 @@
+
+{
+  "output": {
+    "path": "./runs/0D_test/",
+    "triggerOutput": 1,
+    "cpuTime": true,
+    "numColl": true,
+    "EMField": true,
+	"folder": "test"
+  },
+  "reference": {
+    "density": 1.0e16,
+    "mass": 6.633e-26,
+    "temperature": 11604.0,
+    "radius": 1.88e-10
+  },
+  "geometry": {
+    "type":     "1DCart",
+    "meshType": "0D",
+    "meshFile": "mesh.msh",
+    "volume": 1e-11
+  },
+  "species": [
+    {"name": "Argon+", "type": "charged", "mass": 6.633e-26, "charge": 1.0, "weight": 1.0e0},
+    {"name": "Argon",  "type": "neutral", "mass": 6.633e-26,                "weight": 1.0e0}
+  ],
+  "case": {
+    "tau":  [1.0e-6, 1.0e-6],
+    "time": 1.0e-3,
+    "pusher": ["0D", "0D"],
+    "initial": [
+      {"speciesName": "Argon+", "initialState": "Argon+_Initial.dat"},
+      {"speciesName": "Argon",  "initialState": "Argon_Initial.dat"}
+    ]
+  },
+  "interactions": {
+    "folderCollisions": "./data/collisions/",
+    "collisions": [
+      {"species_i": "Argon", "species_j": "Argon", 
+       "cTypes": [
+         {"type": "elastic", "crossSection": "EL_Ar-Ar.dat"}
+       ]},
+      {"species_i": "Argon+", "species_j": "Argon", 
+       "cTypes": [
+         {"type": "elastic", "crossSection": "EL_Ar-Ar.dat"}
+       ]}
+    ]
+  },
+  "parallel": {
+    "OpenMP":{
+      "nThreads": 5
+    }
+  }
+}
diff --git a/src/modules/mesh/0D/moduleMesh0D.f90 b/src/modules/mesh/0D/moduleMesh0D.f90
index e945f25..f51efd2 100644
--- a/src/modules/mesh/0D/moduleMesh0D.f90
+++ b/src/modules/mesh/0D/moduleMesh0D.f90
@@ -33,6 +33,7 @@ MODULE moduleMesh0D
     !Init node
     SUBROUTINE initNode0D(self, n, r)
       USE moduleSpecies
+      USE OMP_LIB
       IMPLICIT NONE
 
       CLASS(meshNode0D), INTENT(out):: self
@@ -43,6 +44,8 @@ MODULE moduleMesh0D
 
       ALLOCATE(self%output(1:nSpecies))
 
+      CALL OMP_INIT_LOCK(self%lock)
+
     END SUBROUTINE initNode0D
 
     !Get node coordinates
@@ -113,19 +116,24 @@ MODULE moduleMesh0D
     SUBROUTINE scatter0D(self, part)
       USE moduleMath
       USE moduleSpecies
+      USE OMP_LIB
       IMPLICIT NONE
 
       CLASS(meshVol0D), INTENT(in):: self
       CLASS(particle), INTENT(in):: part
       REAL(8):: tensorS(1:3,1:3)
-      TYPE(outputNode), POINTER:: vertex
+      CLASS(meshNode), POINTER:: node
+      INTEGER:: sp
 
       tensorS = outerProduct(part%v, part%v)
 
-      vertex => self%n1%output(part%species%n)
-      vertex%den          = vertex%den          + part%weight
-      vertex%mom(:)       = vertex%mom(:)       + part%weight*part%v(:)
-      vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*tensorS
+      node => self%n1
+      sp   =  part%species%n
+      CALL OMP_SET_LOCK(node%lock)
+      node%output(sp)%den          = node%output(sp)%den          + part%weight
+      node%output(sp)%mom(:)       = node%output(sp)%mom(:)       + part%weight*part%v(:)
+      node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*tensorS
+      CALL OMP_UNSET_LOCK(node%lock)
 
     END SUBROUTINE scatter0D
 
diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90
index 73a4c1d..e96150c 100644
--- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90
+++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90
@@ -91,6 +91,7 @@ MODULE moduleMesh1DCart
     SUBROUTINE initNode1DCart(self, n, r)
       USE moduleSpecies
       USE moduleRefParam
+      USE OMP_LIB
       IMPLICIT NONE
 
       CLASS(meshNode1DCart), INTENT(out):: self
@@ -105,6 +106,8 @@ MODULE moduleMesh1DCart
       !Allocates output
       ALLOCATE(self%output(1:nSpecies))
 
+      CALL OMP_INIT_LOCK(self%lock)
+
     END SUBROUTINE initNode1DCart
 
     PURE FUNCTION getCoord1DCart(self) RESULT(r)
@@ -378,26 +381,33 @@ MODULE moduleMesh1DCart
     SUBROUTINE scatterSegm(self, part)
       USE moduleMath
       USE moduleSpecies
+      USE OMP_LIB
       IMPLICIT NONE
 
       CLASS(meshVol1DCartSegm), INTENT(in):: self
       CLASS(particle), INTENT(in):: part
-      TYPE(outputNode), POINTER:: vertex
       REAL(8):: w_p(1:2)
       REAL(8):: tensorS(1:3,1:3)
+      CLASS(meshNode), POINTER:: node
+      INTEGER:: sp
 
       w_p = self%weight(part%xi)
       tensorS = outerProduct(part%v, part%v)
 
-      vertex => self%n1%output(part%species%n)
-      vertex%den          = vertex%den          + part%weight*w_p(1)
-      vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(1)*part%v(:)
-      vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS
+      sp   =  part%species%n
+      node => self%n1
+      CALL OMP_SET_LOCK(node%lock)
+      node%output(sp)%den          = node%output(sp)%den          + part%weight*w_p(1)
+      node%output(sp)%mom(:)       = node%output(sp)%mom(:)       + part%weight*w_p(1)*part%v(:)
+      node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(1)*tensorS
+      CALL OMP_UNSET_LOCK(node%lock)
 
-      vertex => self%n2%output(part%species%n)
-      vertex%den          = vertex%den          + part%weight*w_p(2)
-      vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(2)*part%v(:)
-      vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS
+      node => self%n2
+      CALL OMP_SET_LOCK(node%lock)
+      node%output(sp)%den          = node%output(sp)%den          + part%weight*w_p(1)
+      node%output(sp)%mom(:)       = node%output(sp)%mom(:)       + part%weight*w_p(1)*part%v(:)
+      node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(1)*tensorS
+      CALL OMP_UNSET_LOCK(node%lock)
 
     END SUBROUTINE scatterSegm
 
diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90
index d8e5c51..99b3331 100644
--- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90
+++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90
@@ -92,6 +92,7 @@ MODULE moduleMesh1DRad
     SUBROUTINE initNode1DRad(self, n, r)
       USE moduleSpecies
       USE moduleRefParam
+      USE OMP_LIB
       IMPLICIT NONE
 
       CLASS(meshNode1DRad), INTENT(out):: self
@@ -106,6 +107,8 @@ MODULE moduleMesh1DRad
       !Allocates output
       ALLOCATE(self%output(1:nSpecies))
 
+      CALL OMP_INIT_LOCK(self%lock)
+
     END SUBROUTINE initNode1DRad
 
     PURE FUNCTION getCoord1DRad(self) RESULT(r)
@@ -390,26 +393,33 @@ MODULE moduleMesh1DRad
     SUBROUTINE scatterRad(self, part)
       USE moduleMath
       USE moduleSpecies
+      USE OMP_LIB
       IMPLICIT NONE
 
       CLASS(meshVol1DRadSegm), INTENT(in):: self
       CLASS(particle), INTENT(in):: part
-      TYPE(outputNode), POINTER:: vertex
       REAL(8):: w_p(1:2)
       REAL(8):: tensorS(1:3,1:3)
+      CLASS(meshNode), POINTER:: node
+      INTEGER:: sp
 
       w_p = self%weight(part%xi)
       tensorS = outerProduct(part%v, part%v)
 
-      vertex => self%n1%output(part%species%n)
-      vertex%den          = vertex%den          + part%weight*w_p(1)
-      vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(1)*part%v(:)
-      vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS
+      sp   =  part%species%n
+      node => self%n1
+      CALL OMP_SET_LOCK(node%lock)
+      node%output(sp)%den          = node%output(sp)%den          + part%weight*w_p(1)
+      node%output(sp)%mom(:)       = node%output(sp)%mom(:)       + part%weight*w_p(1)*part%v(:)
+      node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(1)*tensorS
+      CALL OMP_UNSET_LOCK(node%lock)
 
-      vertex => self%n2%output(part%species%n)
-      vertex%den          = vertex%den          + part%weight*w_p(2)
-      vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(2)*part%v(:)
-      vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS
+      node => self%n2
+      CALL OMP_SET_LOCK(node%lock)
+      node%output(sp)%den          = node%output(sp)%den          + part%weight*w_p(2)
+      node%output(sp)%mom(:)       = node%output(sp)%mom(:)       + part%weight*w_p(2)*part%v(:)
+      node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(2)*tensorS
+      CALL OMP_UNSET_LOCK(node%lock)
 
     END SUBROUTINE scatterRad
 
diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90
index c830f8d..36a182a 100644
--- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90
+++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90
@@ -130,6 +130,7 @@ MODULE moduleMesh2DCart
     SUBROUTINE initNode2DCart(self, n, r)
       USE moduleSpecies
       USE moduleRefParam
+      USE OMP_LIB
       IMPLICIT NONE
 
       CLASS(meshNode2DCart), INTENT(out):: self
@@ -145,6 +146,8 @@ MODULE moduleMesh2DCart
       !Allocates output:
       ALLOCATE(self%output(1:nSpecies))
 
+      CALL OMP_INIT_LOCK(self%lock)
+
     END SUBROUTINE initNode2DCart
 
     !Get coordinates from node
@@ -494,36 +497,47 @@ MODULE moduleMesh2DCart
     SUBROUTINE scatterQuad(self, part)
       USE moduleMath
       USE moduleSpecies
+      USE OMP_LIB
       IMPLICIT NONE
 
       CLASS(meshVol2DCartQuad), INTENT(in):: self
       CLASS(particle), INTENT(in):: part
-      TYPE(outputNode), POINTER:: vertex
       REAL(8):: w_p(1:4)
       REAL(8):: tensorS(1:3,1:3)
+      CLASS(meshNode), POINTER:: node
+      INTEGER:: sp
 
       w_p = self%weight(part%xi)
       tensorS = outerProduct(part%v, part%v)
 
-      vertex => self%n1%output(part%species%n)
-      vertex%den          = vertex%den          + part%weight*w_p(1)
-      vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(1)*part%v(:)
-      vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS
+      sp   =  part%species%n
+      node => self%n1
+      CALL OMP_SET_LOCK(node%lock)
+      node%output(sp)%den          = node%output(sp)%den          + part%weight*w_p(1)
+      node%output(sp)%mom(:)       = node%output(sp)%mom(:)       + part%weight*w_p(1)*part%v(:)
+      node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(1)*tensorS
+      CALL OMP_UNSET_LOCK(node%lock)
 
-      vertex => self%n2%output(part%species%n)
-      vertex%den          = vertex%den          + part%weight*w_p(2)
-      vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(2)*part%v(:)
-      vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS
+      node => self%n2
+      CALL OMP_SET_LOCK(node%lock)
+      node%output(sp)%den          = node%output(sp)%den          + part%weight*w_p(2)
+      node%output(sp)%mom(:)       = node%output(sp)%mom(:)       + part%weight*w_p(2)*part%v(:)
+      node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(2)*tensorS
+      CALL OMP_UNSET_LOCK(node%lock)
 
-      vertex => self%n3%output(part%species%n)
-      vertex%den          = vertex%den          + part%weight*w_p(3)
-      vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(3)*part%v(:)
-      vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(3)*tensorS
+      node => self%n3
+      CALL OMP_SET_LOCK(node%lock)
+      node%output(sp)%den          = node%output(sp)%den          + part%weight*w_p(3)
+      node%output(sp)%mom(:)       = node%output(sp)%mom(:)       + part%weight*w_p(3)*part%v(:)
+      node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(3)*tensorS
+      CALL OMP_UNSET_LOCK(node%lock)
 
-      vertex => self%n4%output(part%species%n)
-      vertex%den          = vertex%den          + part%weight*w_p(4)
-      vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(4)*part%v(:)
-      vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(4)*tensorS
+      node => self%n4
+      CALL OMP_SET_LOCK(node%lock)
+      node%output(sp)%den          = node%output(sp)%den          + part%weight*w_p(4)
+      node%output(sp)%mom(:)       = node%output(sp)%mom(:)       + part%weight*w_p(4)*part%v(:)
+      node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(4)*tensorS
+      CALL OMP_UNSET_LOCK(node%lock)
 
     END SUBROUTINE scatterQuad
 
@@ -854,31 +868,40 @@ MODULE moduleMesh2DCart
     SUBROUTINE scatterTria(self, part)
       USE moduleMath
       USE moduleSpecies
+      USE OMP_LIB
       IMPLICIT NONE
 
       CLASS(meshVol2DCartTria), INTENT(in):: self
       CLASS(particle), INTENT(in):: part
-      TYPE(outputNode), POINTER:: vertex
       REAL(8):: w_p(1:3)
       REAL(8):: tensorS(1:3,1:3)
+      CLASS(meshNode), POINTER:: node
+      INTEGER:: sp
 
       w_p = self%weight(part%xi)
       tensorS = outerProduct(part%v, part%v)
 
-      vertex => self%n1%output(part%species%n)
-      vertex%den          = vertex%den          + part%weight*w_p(1)
-      vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(1)*part%v(:)
-      vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS
+      sp   =  part%species%n
+      node => self%n1
+      CALL OMP_SET_LOCK(node%lock)
+      node%output(sp)%den          = node%output(sp)%den          + part%weight*w_p(1)
+      node%output(sp)%mom(:)       = node%output(sp)%mom(:)       + part%weight*w_p(1)*part%v(:)
+      node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(1)*tensorS
+      CALL OMP_UNSET_LOCK(node%lock)
 
-      vertex => self%n2%output(part%species%n)
-      vertex%den          = vertex%den          + part%weight*w_p(2)
-      vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(2)*part%v(:)
-      vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS
+      node => self%n2
+      CALL OMP_SET_LOCK(node%lock)
+      node%output(sp)%den          = node%output(sp)%den          + part%weight*w_p(2)
+      node%output(sp)%mom(:)       = node%output(sp)%mom(:)       + part%weight*w_p(2)*part%v(:)
+      node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(2)*tensorS
+      CALL OMP_UNSET_LOCK(node%lock)
 
-      vertex => self%n3%output(part%species%n)
-      vertex%den          = vertex%den          + part%weight*w_p(3)
-      vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(3)*part%v(:)
-      vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(3)*tensorS
+      node => self%n3
+      CALL OMP_SET_LOCK(node%lock)
+      node%output(sp)%den          = node%output(sp)%den          + part%weight*w_p(3)
+      node%output(sp)%mom(:)       = node%output(sp)%mom(:)       + part%weight*w_p(3)*part%v(:)
+      node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(3)*tensorS
+      CALL OMP_UNSET_LOCK(node%lock)
 
     END SUBROUTINE scatterTria
 
diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
index 80c0665..d4e06bd 100644
--- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
+++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
@@ -131,6 +131,7 @@ MODULE moduleMesh2DCyl
     SUBROUTINE initNode2DCyl(self, n, r)
       USE moduleSpecies
       USE moduleRefParam
+      USE OMP_LIB
       IMPLICIT NONE
 
       CLASS(meshNode2DCyl), INTENT(out):: self
@@ -146,6 +147,8 @@ MODULE moduleMesh2DCyl
       !Allocates output:
       ALLOCATE(self%output(1:nSpecies))
 
+      CALL OMP_INIT_LOCK(self%lock)
+
     END SUBROUTINE initNode2DCyl
 
     !Get coordinates from node
@@ -515,36 +518,47 @@ MODULE moduleMesh2DCyl
     SUBROUTINE scatterQuad(self, part)
       USE moduleMath
       USE moduleSpecies
+      USE OMP_LIB
       IMPLICIT NONE
 
       CLASS(meshVol2DCylQuad), INTENT(in):: self
       CLASS(particle), INTENT(in):: part
-      TYPE(outputNode), POINTER:: vertex
       REAL(8):: w_p(1:4)
       REAL(8):: tensorS(1:3,1:3)
+      CLASS(meshNode), POINTER:: node
+      INTEGER:: sp
 
       w_p = self%weight(part%xi)
       tensorS = outerProduct(part%v, part%v)
 
-      vertex => self%n1%output(part%species%n)
-      vertex%den          = vertex%den          + part%weight*w_p(1)
-      vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(1)*part%v(:)
-      vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS
+      sp   =  part%species%n
+      node => self%n1
+      CALL OMP_SET_LOCK(node%lock)
+      node%output(sp)%den          = node%output(sp)%den          + part%weight*w_p(1)
+      node%output(sp)%mom(:)       = node%output(sp)%mom(:)       + part%weight*w_p(1)*part%v(:)
+      node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(1)*tensorS
+      CALL OMP_UNSET_LOCK(node%lock)
 
-      vertex => self%n2%output(part%species%n)
-      vertex%den          = vertex%den          + part%weight*w_p(2)
-      vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(2)*part%v(:)
-      vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS
+      node => self%n2
+      CALL OMP_SET_LOCK(node%lock)
+      node%output(sp)%den          = node%output(sp)%den          + part%weight*w_p(2)
+      node%output(sp)%mom(:)       = node%output(sp)%mom(:)       + part%weight*w_p(2)*part%v(:)
+      node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(2)*tensorS
+      CALL OMP_UNSET_LOCK(node%lock)
 
-      vertex => self%n3%output(part%species%n)
-      vertex%den          = vertex%den          + part%weight*w_p(3)
-      vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(3)*part%v(:)
-      vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(3)*tensorS
+      node => self%n3
+      CALL OMP_SET_LOCK(node%lock)
+      node%output(sp)%den          = node%output(sp)%den          + part%weight*w_p(3)
+      node%output(sp)%mom(:)       = node%output(sp)%mom(:)       + part%weight*w_p(3)*part%v(:)
+      node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(3)*tensorS
+      CALL OMP_UNSET_LOCK(node%lock)
 
-      vertex => self%n4%output(part%species%n)
-      vertex%den          = vertex%den          + part%weight*w_p(4)
-      vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(4)*part%v(:)
-      vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(4)*tensorS
+      node => self%n4
+      CALL OMP_SET_LOCK(node%lock)
+      node%output(sp)%den          = node%output(sp)%den          + part%weight*w_p(4)
+      node%output(sp)%mom(:)       = node%output(sp)%mom(:)       + part%weight*w_p(4)*part%v(:)
+      node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(4)*tensorS
+      CALL OMP_UNSET_LOCK(node%lock)
 
     END SUBROUTINE scatterQuad
 
@@ -884,31 +898,40 @@ MODULE moduleMesh2DCyl
     SUBROUTINE scatterTria(self, part)
       USE moduleMath
       USE moduleSpecies
+      USE OMP_LIB
       IMPLICIT NONE
 
       CLASS(meshVol2DCylTria), INTENT(in):: self
       CLASS(particle), INTENT(in):: part
-      TYPE(outputNode), POINTER:: vertex
       REAL(8):: w_p(1:3)
       REAL(8):: tensorS(1:3,1:3)
+      CLASS(meshNode), POINTER:: node
+      INTEGER:: sp
 
       w_p = self%weight(part%xi)
       tensorS = outerProduct(part%v, part%v)
 
-      vertex => self%n1%output(part%species%n)
-      vertex%den          = vertex%den          + part%weight*w_p(1)
-      vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(1)*part%v(:)
-      vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS
+      sp   =  part%species%n
+      node => self%n1
+      CALL OMP_SET_LOCK(node%lock)
+      node%output(sp)%den          = node%output(sp)%den          + part%weight*w_p(1)
+      node%output(sp)%mom(:)       = node%output(sp)%mom(:)       + part%weight*w_p(1)*part%v(:)
+      node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(1)*tensorS
+      CALL OMP_UNSET_LOCK(node%lock)
 
-      vertex => self%n2%output(part%species%n)
-      vertex%den          = vertex%den          + part%weight*w_p(2)
-      vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(2)*part%v(:)
-      vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS
+      node => self%n2
+      CALL OMP_SET_LOCK(node%lock)
+      node%output(sp)%den          = node%output(sp)%den          + part%weight*w_p(2)
+      node%output(sp)%mom(:)       = node%output(sp)%mom(:)       + part%weight*w_p(2)*part%v(:)
+      node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(2)*tensorS
+      CALL OMP_UNSET_LOCK(node%lock)
 
-      vertex => self%n3%output(part%species%n)
-      vertex%den          = vertex%den          + part%weight*w_p(3)
-      vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(3)*part%v(:)
-      vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(3)*tensorS
+      node => self%n3
+      CALL OMP_SET_LOCK(node%lock)
+      node%output(sp)%den          = node%output(sp)%den          + part%weight*w_p(3)
+      node%output(sp)%mom(:)       = node%output(sp)%mom(:)       + part%weight*w_p(3)*part%v(:)
+      node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(3)*tensorS
+      CALL OMP_UNSET_LOCK(node%lock)
 
     END SUBROUTINE scatterTria
 
diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90
index e65ff25..8c6e350 100644
--- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90
+++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90
@@ -92,6 +92,7 @@ MODULE moduleMesh3DCart
     SUBROUTINE initNode3DCart(self, n, r)
       USE moduleSpecies
       USE moduleRefParam
+      USE OMP_LIB
       IMPLICIT NONE
 
       CLASS(meshNode3DCart), INTENT(out):: self
@@ -108,6 +109,8 @@ MODULE moduleMesh3DCart
       !Allocates output:
       ALLOCATE(self%output(1:nSpecies))
 
+      CALL OMP_INIT_LOCK(self%lock)
+
     END SUBROUTINE initNode3DCart
 
     !Get coordinates from node
@@ -481,36 +484,47 @@ MODULE moduleMesh3DCart
     SUBROUTINE scatterTetra(self, part)
       USE moduleMath
       USE moduleSpecies
+      USE OMP_LIB
       IMPLICIT NONE
 
       CLASS(meshVol3DCartTetra), INTENT(in):: self
       CLASS(particle), INTENT(in):: part
-      TYPE(outputNode), POINTER:: vertex
       REAL(8):: w_p(1:4)
       REAL(8):: tensorS(1:3, 1:3)
+      CLASS(meshNode), POINTER:: node
+      INTEGER:: sp
 
       w_p = self%weight(part%xi)
       tensorS = outerProduct(part%v, part%v)
 
-      vertex => self%n1%output(part%species%n)
-      vertex%den          = vertex%den          + part%weight*w_p(1)
-      vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(1)*part%v(:)
-      vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS
+      sp   =  part%species%n
+      node => self%n1
+      CALL OMP_SET_LOCK(node%lock)
+      node%output(sp)%den          = node%output(sp)%den          + part%weight*w_p(1)
+      node%output(sp)%mom(:)       = node%output(sp)%mom(:)       + part%weight*w_p(1)*part%v(:)
+      node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(1)*tensorS
+      CALL OMP_UNSET_LOCK(node%lock)
 
-      vertex => self%n2%output(part%species%n)
-      vertex%den          = vertex%den          + part%weight*w_p(2)
-      vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(2)*part%v(:)
-      vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS
+      node => self%n2
+      CALL OMP_SET_LOCK(node%lock)
+      node%output(sp)%den          = node%output(sp)%den          + part%weight*w_p(2)
+      node%output(sp)%mom(:)       = node%output(sp)%mom(:)       + part%weight*w_p(2)*part%v(:)
+      node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(2)*tensorS
+      CALL OMP_UNSET_LOCK(node%lock)
 
-      vertex => self%n3%output(part%species%n)
-      vertex%den          = vertex%den          + part%weight*w_p(3)
-      vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(3)*part%v(:)
-      vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(3)*tensorS
+      node => self%n3
+      CALL OMP_SET_LOCK(node%lock)
+      node%output(sp)%den          = node%output(sp)%den          + part%weight*w_p(3)
+      node%output(sp)%mom(:)       = node%output(sp)%mom(:)       + part%weight*w_p(3)*part%v(:)
+      node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(3)*tensorS
+      CALL OMP_UNSET_LOCK(node%lock)
 
-      vertex => self%n4%output(part%species%n)
-      vertex%den          = vertex%den          + part%weight*w_p(4)
-      vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(4)*part%v(:)
-      vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(4)*tensorS
+      node => self%n4
+      CALL OMP_SET_LOCK(node%lock)
+      node%output(sp)%den          = node%output(sp)%den          + part%weight*w_p(4)
+      node%output(sp)%mom(:)       = node%output(sp)%mom(:)       + part%weight*w_p(4)*part%v(:)
+      node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(4)*tensorS
+      CALL OMP_UNSET_LOCK(node%lock)
 
     END SUBROUTINE scatterTetra
 
diff --git a/src/modules/mesh/inout/0D/moduleMeshInput0D.f90 b/src/modules/mesh/inout/0D/moduleMeshInput0D.f90
index 678671a..5ac0682 100644
--- a/src/modules/mesh/inout/0D/moduleMeshInput0D.f90
+++ b/src/modules/mesh/inout/0D/moduleMeshInput0D.f90
@@ -19,6 +19,7 @@ MODULE moduleMeshInput0D
         TYPE IS(meshParticles)
           self%printOutput => printOutput0D
           self%printEM     => printEM0D
+          self%readInitial => readInitial0D
 
         END SELECT
         self%readMesh => read0D
@@ -63,4 +64,37 @@ MODULE moduleMeshInput0D
 
     END SUBROUTINE read0D
 
+    SUBROUTINE readInitial0D(sp, filename, density, velocity, temperature)
+      IMPLICIT NONE
+
+      INTEGER, INTENT(in):: sp
+      CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
+      REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:):: density
+      REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:,:):: velocity
+      REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:):: temperature
+      REAL(8):: dummy
+      INTEGER:: stat
+
+      ALLOCATE(density(1:1))
+      ALLOCATE(velocity(1:1, 1:3))
+      ALLOCATE(temperature(1:1))
+
+      OPEN(10, file = TRIM(filename))
+      
+      !Finds the last line of data
+      stat = 0
+      DO WHILE (stat==0)
+        READ(10, *, iostat = stat)
+
+      END DO
+
+      !Go back two line
+      BACKSPACE(10)
+      BACKSPACE(10)
+
+      !Reads data
+      READ(10, *) dummy, density(1), velocity(1, 1:3), dummy, temperature(1)
+
+    END SUBROUTINE readInitial0D
+
 END MODULE moduleMeshInput0D
diff --git a/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 b/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90
index 6e490b4..db52de4 100644
--- a/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90
+++ b/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90
@@ -18,7 +18,8 @@ MODULE moduleMeshOutput0D
         fileName='OUTPUT_' // species(i)%obj%name // '.dat'
         IF (t == 0) THEN
           OPEN(20, file = path // folder // '/' // fileName, action = 'write')
-          WRITE(20, "(A1, 8X, A1, 6(A20))") "#","t","density", "velocity", "pressure", "temperature"
+          WRITE(20, "(A1, 14X, A5, A20, 40X, A20, 2(A20))") "#","t (s)","density (m^-3)", "velocity (m/s)", &
+                                                                        "pressure (Pa)", "temperature (K)"
           WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
           CLOSE(20)
 
@@ -26,7 +27,7 @@ MODULE moduleMeshOutput0D
         
         OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write')
         CALL calculateOutput(self%nodes(1)%obj%output(i), output, self%nodes(1)%obj%v, species(i)%obj)
-        WRITE(20, "(I10, 6(ES20.6E3))") t, output%density, output%velocity, output%pressure, output%temperature
+        WRITE(20, "(7(ES20.6E3))") REAL(t)*tauMin*ti_ref, output%density, output%velocity, output%pressure, output%temperature
         CLOSE(20)
 
       END DO
@@ -43,6 +44,20 @@ MODULE moduleMeshOutput0D
 
       CLASS(meshGeneric), INTENT(in):: self
       INTEGER, INTENT(in):: t
+      CHARACTER(:), ALLOCATABLE:: fileName
+
+      fileName='OUTPUT_Collisions.dat'
+      IF (t == 0) THEN
+        OPEN(20, file = path // folder // '/' // fileName, action = 'write')
+        WRITE(20, "(A1, 14X, A5, A20)") "#","t (s)","collisions"
+        WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
+        CLOSE(20)
+
+      END IF
+      
+      OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write')
+      WRITE(20, "(ES20.6E3, I20)") REAL(t)*tauMin*ti_ref, mesh%vols(1)%obj%nColl
+      CLOSE(20)
 
     END SUBROUTINE printColl0D
 
diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90
index b446d76..d6fe4d7 100644
--- a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90
+++ b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90
@@ -292,7 +292,6 @@ MODULE moduleMeshInputGmsh2
 
     !Reads the initial information from an output file for an species
     SUBROUTINE readInitialGmsh2(sp, filename, density, velocity, temperature)
-      USE moduleRefParam
       IMPLICIT NONE
 
       INTEGER, INTENT(in):: sp
diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90
index 304dd86..1f6c795 100644
--- a/src/modules/mesh/moduleMesh.f90
+++ b/src/modules/mesh/moduleMesh.f90
@@ -20,6 +20,8 @@ MODULE moduleMesh
     !Output values
     TYPE(outputNode), ALLOCATABLE:: output(:)
     TYPE(emNode):: emData
+    !Lock indicator for scattering
+    INTEGER(KIND=OMP_LOCK_KIND):: lock
     CONTAINS
       PROCEDURE(initNode_interface),    DEFERRED, PASS:: init
       PROCEDURE(getCoord_interface),    DEFERRED, PASS:: getCoordinates
diff --git a/src/modules/moduleCollisions.f90 b/src/modules/moduleCollisions.f90
index 73871ce..a2c48ca 100644
--- a/src/modules/moduleCollisions.f90
+++ b/src/modules/moduleCollisions.f90
@@ -193,7 +193,6 @@ MODULE moduleCollisions
       REAL(8), DIMENSION(1:3):: vCM
       REAL(8):: vp(1:3)
 
-      !eRel (in units of [m][L]^2[t]^-2
       vRel = SUM(DABS(part_i%v-part_j%v)) !TODO make function of norm1
       eRel = self%rMass*vRel**2
       sigmaVrel = self%crossSec%get(eRel)*vRel
diff --git a/src/modules/moduleInput.f90 b/src/modules/moduleInput.f90
index fed9228..27cb16e 100644
--- a/src/modules/moduleInput.f90
+++ b/src/modules/moduleInput.f90
@@ -348,7 +348,6 @@ MODULE moduleInput
               !Assign particle to temporal list of particles
               CALL partInitial%add(partNew)
 
-
             END DO
 
             DEALLOCATE(source)
@@ -768,6 +767,7 @@ MODULE moduleInput
       USE moduleMesh1DCart, ONLY: connectMesh1DCart
       USE moduleErrors
       USE moduleOutput
+      USE moduleRefParam
       USE json_module
       IMPLICIT NONE
 
@@ -776,6 +776,7 @@ MODULE moduleInput
       LOGICAL:: doubleMesh
       CHARACTER(:), ALLOCATABLE:: meshFormat, meshFile
       CHARACTER(:), ALLOCATABLE:: fullPath
+      REAL(8):: volume
 
       !Firstly, indicates if a specific mesh for MC collisions is being use
       doubleMesh = ASSOCIATED(meshForMCC, meshColl)
@@ -792,6 +793,7 @@ MODULE moduleInput
         IF (doubleMesh) CALL initGmsh2(meshColl)
 
       CASE ("0D")
+        CALL config%get('geometry.meshType', meshFormat, found)
         CALL init0D(mesh)
 
       CASE DEFAULT
@@ -812,6 +814,18 @@ MODULE moduleInput
 
       END IF
 
+      !Gets the volume for a 0D mesh
+      !TODO: Try to constrain this to the inout for 0D
+      IF (meshFormat == "0D") THEN
+        CALL config%get('geometry.volume', volume, found)
+        IF (found) THEN
+          mesh%vols(1)%obj%volume = mesh%vols(1)%obj%volume*volume / Vol_ref
+          mesh%nodes(1)%obj%v     = mesh%vols(1)%obj%volume
+
+        END IF
+
+      END IF
+
       !Creates the connectivity between elements
       SELECT CASE(mesh%geometry)
       CASE("3DCart")
diff --git a/src/modules/moduleOutput.f90 b/src/modules/moduleOutput.f90
index 89df330..4a3f4da 100644
--- a/src/modules/moduleOutput.f90
+++ b/src/modules/moduleOutput.f90
@@ -96,7 +96,9 @@ MODULE moduleOutput
         IF (PRESENT(first)) THEN
           IF (first) THEN
             OPEN(20, file = path // folder // '/' // fileName, action = 'write')
-            WRITE(20, "(A1, 8X, A1, 9X, A1, 6(A20))") "#","t","n","total","push","reset","collision","weighting","EMField"
+            WRITE(20, "(A1, 8X, A1, 9X, A1, 7(A20))") "#","t","n","total (s)","push (s)","reset (s)",        &
+                                                                              "collision (s)","coulomb (s)", &
+                                                                              "weighting (s)","EMField (s)"
             WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
             CLOSE(20)