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..6940f77
--- /dev/null
+++ b/runs/0D_Argon/input.json
@@ -0,0 +1,51 @@
+
+{
+  "output": {
+    "path": "./runs/0D_test/",
+    "triggerOutput": 1,
+    "numColl": true,
+	"folder": "test"
+  },
+  "reference": {
+    "density": 1.0e16,
+    "mass": 6.633e-26,
+    "temperature": 11604.0,
+    "radius": 1.88e-10
+  },
+  "geometry": {
+    "type":     "1DCart",
+    "meshType": "0D",
+    "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": 4
+    }
+  }
+}
diff --git a/src/makefile b/src/makefile
index 75795c8..3a6d68a 100644
--- a/src/makefile
+++ b/src/makefile
@@ -5,11 +5,14 @@ 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)/moduleMeshInputGmsh2.o $(OBJDIR)/moduleMeshOutputGmsh2.o \
+			  	$(OBJDIR)/moduleMeshInput0D.o $(OBJDIR)/moduleMeshOutput0D.o \
 			  	$(OBJDIR)/moduleMesh3DCart.o \
 			  	$(OBJDIR)/moduleMesh2DCyl.o  \
 			  	$(OBJDIR)/moduleMesh2DCart.o \
           $(OBJDIR)/moduleMesh1DRad.o  \
-          $(OBJDIR)/moduleMesh1DCart.o
+          $(OBJDIR)/moduleMesh1DCart.o \
+          $(OBJDIR)/moduleMesh0D.o
+
 
 all: $(OUTPUT)
 
diff --git a/src/modules/mesh/0D/makefile b/src/modules/mesh/0D/makefile
new file mode 100644
index 0000000..9dcba7e
--- /dev/null
+++ b/src/modules/mesh/0D/makefile
@@ -0,0 +1,5 @@
+all: moduleMesh0D.o
+
+moduleMesh0D.o: moduleMesh0D.f90
+	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
+
diff --git a/src/modules/mesh/0D/moduleMesh0D.f90 b/src/modules/mesh/0D/moduleMesh0D.f90
new file mode 100644
index 0000000..f51efd2
--- /dev/null
+++ b/src/modules/mesh/0D/moduleMesh0D.f90
@@ -0,0 +1,206 @@
+!moduleMesh1D: 0D mesh. No coordinates are relevant. No edges are used
+MODULE moduleMesh0D
+  USE moduleMesh
+  IMPLICIT NONE
+
+  TYPE, PUBLIC, EXTENDS(meshNode):: meshNode0D
+    INTEGER:: n1
+    CONTAINS
+      PROCEDURE, PASS:: init => initNode0D
+      PROCEDURE, PASS:: getCoordinates => getCoord0D
+
+  END TYPE meshNode0D
+
+  TYPE, PUBLIC, EXTENDS(meshVol):: meshVol0D
+    CLASS(meshNode), POINTER:: n1
+    CONTAINS
+      PROCEDURE, PASS:: init => initVol0D
+      PROCEDURE, PASS:: getNodes    => getNodes0D
+      PROCEDURE, PASS:: randPos     => randPos0D
+      PROCEDURE, NOPASS:: fPsi      => fPsi0D
+      PROCEDURE, PASS:: scatter     => scatter0D
+      PROCEDURE, PASS:: gatherEF    => gatherEF0D
+      PROCEDURE, PASS:: elemK       => elemK0D
+      PROCEDURE, PASS:: elemF       => elemF0D
+      PROCEDURE, PASS:: phy2log     => phy2log0D
+      PROCEDURE, NOPASS:: inside    => inside0D
+      PROCEDURE, PASS:: nextElement => nextElement0D
+
+  END TYPE meshVol0D
+
+  CONTAINS
+    !NODE FUNCTIONS
+    !Init node
+    SUBROUTINE initNode0D(self, n, r)
+      USE moduleSpecies
+      USE OMP_LIB
+      IMPLICIT NONE
+
+      CLASS(meshNode0D), INTENT(out):: self
+      INTEGER, INTENT(in):: n
+      REAL(8), INTENT(in):: r(1:3) !Unused variable
+
+      self%n = n
+
+      ALLOCATE(self%output(1:nSpecies))
+
+      CALL OMP_INIT_LOCK(self%lock)
+
+    END SUBROUTINE initNode0D
+
+    !Get node coordinates
+    PURE FUNCTION getCoord0D(self) RESULT(r)
+      IMPLICIT NONE
+
+      CLASS(meshNode0D), INTENT(in):: self
+      REAL(8):: r(1:3)
+
+      r = 0.D0
+
+    END FUNCTION
+
+    !VOLUME FUNCTIONS
+    !Inits dummy 0D volume
+    SUBROUTINE initVol0D(self, n, p, nodes)
+      USE moduleRefParam
+      IMPLICIT NONE
+
+      CLASS(meshVol0D), INTENT(out):: self
+      INTEGER, INTENT(in):: n
+      INTEGER, INTENT(in):: p(:)
+      TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:)
+
+      self%n = n
+
+      self%n1 => nodes(p(1))%obj
+      self%volume = 1.D0
+      self%n1%v = 1.D0
+
+      self%sigmaVrelMax = sigma_ref/L_ref**2
+
+      CALL OMP_INIT_LOCK(self%lock)
+
+    END SUBROUTINE initVol0D
+
+    PURE FUNCTION getNodes0D(self) RESULT(n)
+      IMPLICIT NONE
+
+      CLASS(meshVol0D), INTENT(in):: self
+      INTEGER, ALLOCATABLE:: n(:)
+
+      ALLOCATE(n(1:1))
+
+      n = self%n1%n
+
+    END FUNCTION getNodes0D
+
+    FUNCTION randPos0D(self) RESULT(r)
+      IMPLICIT NONE
+
+      CLASS(meshVol0D), INTENT(in):: self
+      REAL(8):: r(1:3)
+
+      r = 0.D0
+
+    END FUNCTION randPos0D
+
+    PURE FUNCTION fPsi0D(xi) RESULT(fPsi)
+      REAL(8), INTENT(in):: xi(1:3)
+      REAL(8), ALLOCATABLE:: fPsi(:)
+
+      ALLOCATE(fPsi(1:1))
+      fPsi = 1.D0
+
+    END FUNCTION fPsi0D
+
+    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)
+      CLASS(meshNode), POINTER:: node
+      INTEGER:: sp
+
+      tensorS = outerProduct(part%v, part%v)
+
+      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
+
+    PURE FUNCTION gatherEF0D(self, xi) RESULT(EF)
+      IMPLICIT NONE
+
+      CLASS(meshVol0D), INTENT(in):: self
+      REAL(8), INTENT(in):: xi(1:3)
+      REAL(8):: EF(1:3)
+
+      EF = 0.D0
+
+    END FUNCTION gatherEF0D
+
+    PURE FUNCTION elemK0D(self) RESULT(localK)
+      IMPLICIT NONE
+
+      CLASS(meshVol0D), INTENT(in):: self
+      REAL(8), ALLOCATABLE:: localK(:,:)
+
+      ALLOCATE(localK(1:1, 1:1))
+      localK = 0.D0
+
+    END FUNCTION elemK0D
+
+    PURE FUNCTION elemF0D(self, source) RESULT(localF)
+      IMPLICIT NONE
+
+      CLASS(meshVol0D), INTENT(in):: self
+      REAL(8), INTENT(in):: source(1:)
+      REAL(8), ALLOCATABLE:: localF(:)
+
+      ALLOCATE(localF(1:1))
+      localF = 0.D0
+
+    END FUNCTION elemF0D
+
+    PURE FUNCTION phy2log0D(self,r) RESULT(xN) 
+      IMPLICIT NONE
+
+      CLASS(meshVol0D), INTENT(in):: self
+      REAL(8), INTENT(in):: r(1:3)
+      REAL(8):: xN(1:3)
+
+      xN = 0.D0
+
+    END FUNCTION phy2log0D
+
+    PURE FUNCTION inside0D(xi) RESULT(ins)
+      IMPLICIT NONE
+
+      REAL(8), INTENT(in):: xi(1:3)
+      LOGICAL:: ins
+
+      ins = .TRUE.
+
+    END FUNCTION inside0D
+
+    SUBROUTINE nextElement0D(self, xi, nextElement)
+      IMPLICIT NONE
+
+      CLASS(meshVol0D), INTENT(in):: self
+      REAL(8), INTENT(in):: xi(1:3)
+      CLASS(meshElement), POINTER, INTENT(out):: nextElement
+
+      nextElement => NULL()
+
+    END SUBROUTINE nextElement0D
+
+END MODULE moduleMesh0D
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/makefile b/src/modules/mesh/inout/0D/makefile
new file mode 100644
index 0000000..4bf91b7
--- /dev/null
+++ b/src/modules/mesh/inout/0D/makefile
@@ -0,0 +1,7 @@
+all: moduleMeshInput0D.o moduleMeshOutput0D.o
+
+moduleMeshInput0D.o: moduleMeshOutput0D.o moduleMeshInput0D.f90
+	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
+
+%.o: %.f90
+	$(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@
diff --git a/src/modules/mesh/inout/0D/moduleMeshInput0D.f90 b/src/modules/mesh/inout/0D/moduleMeshInput0D.f90
new file mode 100644
index 0000000..5ac0682
--- /dev/null
+++ b/src/modules/mesh/inout/0D/moduleMeshInput0D.f90
@@ -0,0 +1,100 @@
+MODULE moduleMeshInput0D
+  !Creates a 0D mesh. mostly used for testing collisional processes.
+  !This mesh consists on 1 node and 1 volume in which all particles are located.
+  !No boundary conditions and no pusher should be applied to this geometry.
+  !Output should go to a single file (per species) in which each row represents an iteration.
+  !In principle, no EM field is applied.
+
+  CONTAINS
+    !Inits the 0D mesh
+    SUBROUTINE init0D(self)
+      USE moduleMesh
+      USE moduleMeshOutput0D
+      IMPLICIT NONE
+
+      CLASS(meshGeneric), INTENT(inout), TARGET:: self
+
+        IF (ASSOCIATED(meshForMCC, self)) self%printColl => printColl0D
+        SELECT TYPE(self)
+        TYPE IS(meshParticles)
+          self%printOutput => printOutput0D
+          self%printEM     => printEM0D
+          self%readInitial => readInitial0D
+
+        END SELECT
+        self%readMesh => read0D
+
+    END SUBROUTINE init0D
+
+    !Reads a 0D mesh file.
+    !No reading is actually done as the 0D mesh is a fixed one
+    SUBROUTINE read0D(self, filename)
+      USE moduleMesh
+      USE moduleMesh0D
+      IMPLICIT NONE
+
+      CLASS(meshGeneric), INTENT(inout):: self
+      CHARACTER(:), ALLOCATABLE, INTENT(in):: filename !Dummy file, not used
+      REAL(8):: r(1:3)
+
+      !Allocates one node
+      self%numNodes = 1
+      ALLOCATE(self%nodes(1:1))
+      !Allocates one volume
+      self%numVols = 1
+      ALLOCATE(self%vols(1:1))
+      !Allocates matrix K
+      SELECT TYPE(self)
+      TYPE IS(meshParticles)
+        ALLOCATE(self%K(1:1, 1:1))
+        ALLOCATE(self%IPIV(1:1, 1:1))
+        self%K = 0.D0
+        self%IPIV = 0.D0
+
+      END SELECT
+
+      !Creates the node
+      ALLOCATE(meshNode0D:: self%nodes(1)%obj)
+      r = 0.D0
+      CALL self%nodes(1)%obj%init(1, r)
+
+      !Creates the volume
+      ALLOCATE(meshVol0D:: self%vols(1)%obj)
+      CALL self%vols(1)%obj%init(1, (/ 1/), self%nodes)
+
+    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
new file mode 100644
index 0000000..db52de4
--- /dev/null
+++ b/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90
@@ -0,0 +1,76 @@
+MODULE moduleMeshOutput0D
+
+  CONTAINS
+    SUBROUTINE printOutput0D(self, t)
+      USE moduleMesh
+      USE moduleRefParam
+      USE moduleSpecies
+      USE moduleOutput
+      IMPLICIT NONE
+
+      CLASS(meshParticles), INTENT(in):: self
+      INTEGER, INTENT(in):: t
+      INTEGER:: i
+      TYPE(outputFormat):: output
+      CHARACTER(:), ALLOCATABLE:: fileName
+
+      DO i = 1, nSpecies
+        fileName='OUTPUT_' // species(i)%obj%name // '.dat'
+        IF (t == 0) THEN
+          OPEN(20, file = path // folder // '/' // fileName, action = 'write')
+          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)
+
+        END IF
+        
+        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, "(7(ES20.6E3))") REAL(t)*tauMin*ti_ref, output%density, output%velocity, output%pressure, output%temperature
+        CLOSE(20)
+
+      END DO
+
+    END SUBROUTINE printOutput0D
+
+    SUBROUTINE printColl0D(self, t)
+      USE moduleMesh
+      USE moduleRefParam
+      USE moduleCaseParam
+      USE moduleCollisions
+      USE moduleOutput
+      IMPLICIT NONE
+
+      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
+
+    SUBROUTINE printEM0D(self, t)
+      USE moduleMesh
+      USE moduleRefParam
+      USE moduleCaseParam
+      USE moduleOutput
+      IMPLICIT NONE
+
+      CLASS(meshParticles), INTENT(in):: self
+      INTEGER, INTENT(in):: t
+
+    END SUBROUTINE printEM0D
+
+END MODULE moduleMeshOutput0D
diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90
index 35e6fce..d6fe4d7 100644
--- a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90
+++ b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90
@@ -1,4 +1,5 @@
 MODULE moduleMeshInputGmsh2
+  !Reads a mesh in the Gmsh v2.0 format
 
   CONTAINS
     !Inits a mesh to use Gmsh2 format
@@ -291,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/inout/gmsh2/moduleMeshOutputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90
index bec5449..43b1ef8 100644
--- a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90
+++ b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90
@@ -96,8 +96,8 @@ MODULE moduleMeshOutputGmsh2
       IMPLICIT NONE
 
       CLASS(meshGeneric), INTENT(in):: self
-      INTEGER:: numEdges
       INTEGER, INTENT(in):: t
+      INTEGER:: numEdges
       INTEGER:: n
       REAL(8):: time
       CHARACTER(:), ALLOCATABLE:: fileName
diff --git a/src/modules/mesh/inout/makefile b/src/modules/mesh/inout/makefile
index 7817608..2f73e73 100644
--- a/src/modules/mesh/inout/makefile
+++ b/src/modules/mesh/inout/makefile
@@ -1,4 +1,7 @@
-all: gmsh2.o
+all: gmsh2.o 0D.o
 
 gmsh2.o:
 	$(MAKE) -C gmsh2 all
+
+0D.o:
+	$(MAKE) -C 0D all
diff --git a/src/modules/mesh/makefile b/src/modules/mesh/makefile
index 0adc8f6..a013878 100644
--- a/src/modules/mesh/makefile
+++ b/src/modules/mesh/makefile
@@ -1,4 +1,4 @@
-all: moduleMesh.o moduleMeshBoundary.o inout.o 3DCart.o 2DCyl.o 2DCart.o 1DRad.o 1DCart.o
+all: moduleMesh.o moduleMeshBoundary.o inout.o 3DCart.o 2DCyl.o 2DCart.o 1DRad.o 1DCart.o 0D.o
 
 3DCart.o: moduleMesh.o
 	$(MAKE) -C 3DCart all
@@ -15,11 +15,14 @@ all: moduleMesh.o moduleMeshBoundary.o inout.o 3DCart.o 2DCyl.o 2DCart.o 1DRad.o
 1DRad.o: moduleMesh.o
 	$(MAKE) -C 1DRad all
 
+0D.o: moduleMesh.o
+	$(MAKE) -C 0D all
+
 moduleMesh.o: moduleMesh.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 
 moduleMeshBoundary.o: moduleMesh.o moduleMeshBoundary.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 
-inout.o: 3DCart.o 2DCyl.o 2DCart.o 1DRad.o 1DCart.o
+inout.o: 3DCart.o 2DCyl.o 2DCart.o 1DRad.o 1DCart.o 0D.o
 	$(MAKE) -C inout all
diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90
index 1deb764..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
@@ -194,7 +196,6 @@ MODULE moduleMesh
     END SUBROUTINE scatter_interface
 
     PURE FUNCTION gatherEF_interface(self, xi) RESULT(EF)
-      
       IMPORT:: meshVol
       CLASS(meshVol), INTENT(in):: self
       REAL(8), INTENT(in):: xi(1:3)
@@ -672,6 +673,8 @@ MODULE moduleMesh
     END SUBROUTINE doCollisions
 
     SUBROUTINE doCoulomb(self)
+      IMPLICIT NONE
+
       CLASS(meshParticles), INTENT(inout):: self
 
     END SUBROUTINE doCoulomb
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 8699f69..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)
@@ -760,6 +759,7 @@ MODULE moduleInput
     SUBROUTINE readGeometry(config)
       USE moduleMesh
       USE moduleMeshInputGmsh2, ONLY: initGmsh2
+      USE moduleMeshInput0D,    ONLY: init0D
       USE moduleMesh3DCart, ONLY: connectMesh3DCart
       USE moduleMesh2DCyl,  ONLY: connectMesh2DCyl
       USE moduleMesh2DCart, ONLY: connectMesh2DCart
@@ -767,6 +767,7 @@ MODULE moduleInput
       USE moduleMesh1DCart, ONLY: connectMesh1DCart
       USE moduleErrors
       USE moduleOutput
+      USE moduleRefParam
       USE json_module
       IMPLICIT NONE
 
@@ -775,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)
@@ -790,6 +792,10 @@ MODULE moduleInput
         CALL initGmsh2(mesh)
         IF (doubleMesh) CALL initGmsh2(meshColl)
 
+      CASE ("0D")
+        CALL config%get('geometry.meshType', meshFormat, found)
+        CALL init0D(mesh)
+
       CASE DEFAULT
         CALL criticalError("Mesh format " // meshFormat // " not recogniced", "readGeometry")
 
@@ -808,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")
@@ -825,8 +843,11 @@ MODULE moduleInput
       CASE("1DCart")
         mesh%connectMesh => connectMesh1DCart
 
+      CASE("0D")
+        mesh%connectMesh => NULL()
+
       END SELECT
-      CALL mesh%connectMesh
+      IF (ASSOCIATED(mesh%connectMesh)) CALL mesh%connectMesh
 
       IF (doubleMesh) THEN
         meshColl%connectMesh => mesh%connectMesh
diff --git a/src/modules/moduleOutput.f90 b/src/modules/moduleOutput.f90
index a9efdc4..4a3f4da 100644
--- a/src/modules/moduleOutput.f90
+++ b/src/modules/moduleOutput.f90
@@ -96,19 +96,18 @@ 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
-
-          ELSE
+            CLOSE(20)
 
           END IF
-            OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write')
-
-        ELSE
-          OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write')
 
         END IF
 
+        OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write')
+
         WRITE (20, "(I10, I10, 7(ES20.6E3))") t, nPartOld, tStep, tPush, tReset, tColl, tCoul, tWeight, tEMField
 
         CLOSE(20)
diff --git a/src/modules/moduleSolver.f90 b/src/modules/moduleSolver.f90
index 089ae71..84fa4e2 100644
--- a/src/modules/moduleSolver.f90
+++ b/src/modules/moduleSolver.f90
@@ -102,6 +102,9 @@ MODULE moduleSolver
       CASE('1DRadCharged')
         self%pushParticle => push1DRadCharged
 
+      CASE('0D')
+        self%pushParticle => push0D
+
       CASE DEFAULT
         CALL criticalError('Pusher ' // pusherType // ' not found','initPusher')
 
@@ -483,6 +486,17 @@ MODULE moduleSolver
 
     END SUBROUTINE push1DRadCharged
 
+    !Dummy pusher for 0D geometry
+    PURE SUBROUTINE push0D(part, tauIn)
+      USE moduleSpecies
+      USE moduleEM
+      IMPLICIT NONE
+
+      TYPE(particle), INTENT(inout):: part
+      REAL(8), INTENT(in):: tauIn
+
+    END SUBROUTINE push0D
+
     SUBROUTINE doReset()
       USE moduleSpecies
       USE moduleMesh