diff --git a/runs/cylFlow/compCPUTime.gp b/runs/cylFlow/compCPUTime.gp
index 38da2d6..6ec49f8 100644
--- a/runs/cylFlow/compCPUTime.gp
+++ b/runs/cylFlow/compCPUTime.gp
@@ -10,10 +10,10 @@ set style line 10       lt 1 lw 2.0 lc rgb "black"   #Black solid line
 
 name1   = "Base Case"
 folder1 = "base_case/"
-name2   = "New Version"
-folder2 = "2020-10-17_14.10.29/"
-name3   = "Reset with List"
-folder3 = "2020-10-13_18.29.33/"
+name2   = "OMP LOCK 1"
+folder2 = "OMP_LOCK_1/"
+name3   = "OMP LOCK 2"
+folder3 = "2020-10-17_18.33.45/"
 name4   = "Reset with List Parallel"
 folder4 = "2020-10-11_15.20.09/"
 
@@ -65,52 +65,52 @@ plot folder1."cpuTime.dat" u 1:($2) t name1       ls 1, \
 unset multiplot
 	 
 
-set terminal qt enhanced 2 persist size 1600, 1000 font "Times ,10"
-
-set style line 1  pt 4              lc rgb "#B50427" #Squares   red
-set style line 2  pt 6              lc rgb "#3B4CC1" #Circles   blue
-set style line 3  pt 1              lc rgb "#2CA02C" #Crosses   green
-set style line 4  pt 2              lc rgb "#FE7F0E" #Exes      orange
-set style line 5  pt 8              lc rgb "#D6696B" #Triangles light red
-set style line 10       lt 1 lw 2.0 lc rgb "black"   #Black solid line
-
-set key box opaque
-set pointsize 1.5
-
-set ylabel "Time per particle (micros)"
-set xrange [0:500]
-set xlabel "Iteration"
-
-
-set multiplot layout 2,3
-set title "Total"
-plot folder1."cpuTime.dat" u 1:(1e6*$3/$2) t name1       ls 1, \
-     folder2."cpuTime.dat" u 1:(1e6*$3/$2) t name2       ls 2, \
-     folder3."cpuTime.dat" u 1:(1e6*$3/$2) t name3       ls 3, \
-     folder4."cpuTime.dat" u 1:(1e6*$3/$2) t name4       ls 4
-	 
-set title "Push"
-plot folder1."cpuTime.dat" u 1:(1e6*$4/$2) t name1       ls 1, \
-     folder2."cpuTime.dat" u 1:(1e6*$4/$2) t name2       ls 2, \
-     folder3."cpuTime.dat" u 1:(1e6*$4/$2) t name3       ls 3, \
-     folder4."cpuTime.dat" u 1:(1e6*$4/$2) t name4       ls 4
-	 
-set title "Reset"
-plot folder1."cpuTime.dat" u 1:(1e6*$5/$2) t name1       ls 1, \
-     folder2."cpuTime.dat" u 1:(1e6*$5/$2) t name2       ls 2, \
-     folder3."cpuTime.dat" u 1:(1e6*$5/$2) t name3       ls 3, \
-     folder4."cpuTime.dat" u 1:(1e6*$5/$2) t name4       ls 4
-	 
-set title "Collisions"
-plot folder1."cpuTime.dat" u 1:(1e6*$6/$2) t name1       ls 1, \
-     folder2."cpuTime.dat" u 1:(1e6*$6/$2) t name2       ls 2, \
-     folder3."cpuTime.dat" u 1:(1e6*$6/$2) t name3       ls 3, \
-     folder4."cpuTime.dat" u 1:(1e6*$6/$2) t name4       ls 4
-	 
-set title "Weighting"
-plot folder1."cpuTime.dat" u 1:(1e6*$7/$2) t name1       ls 1, \
-     folder2."cpuTime.dat" u 1:(1e6*$7/$2) t name2       ls 2, \
-     folder3."cpuTime.dat" u 1:(1e6*$7/$2) t name3       ls 3, \
-     folder4."cpuTime.dat" u 1:(1e6*$7/$2) t name4       ls 4
-
-unset multiplot
+# set terminal qt enhanced 2 persist size 1600, 1000 font "Times ,10"
+# 
+# set style line 1  pt 4              lc rgb "#B50427" #Squares   red
+# set style line 2  pt 6              lc rgb "#3B4CC1" #Circles   blue
+# set style line 3  pt 1              lc rgb "#2CA02C" #Crosses   green
+# set style line 4  pt 2              lc rgb "#FE7F0E" #Exes      orange
+# set style line 5  pt 8              lc rgb "#D6696B" #Triangles light red
+# set style line 10       lt 1 lw 2.0 lc rgb "black"   #Black solid line
+# 
+# set key box opaque
+# set pointsize 1.5
+# 
+# set ylabel "Time per particle (micros)"
+# set xrange [0:500]
+# set xlabel "Iteration"
+# 
+# 
+# set multiplot layout 2,3
+# set title "Total"
+# plot folder1."cpuTime.dat" u 1:(1e6*$3/$2) t name1       ls 1, \
+#      folder2."cpuTime.dat" u 1:(1e6*$3/$2) t name2       ls 2, \
+#      folder3."cpuTime.dat" u 1:(1e6*$3/$2) t name3       ls 3, \
+#      folder4."cpuTime.dat" u 1:(1e6*$3/$2) t name4       ls 4
+# 	 
+# set title "Push"
+# plot folder1."cpuTime.dat" u 1:(1e6*$4/$2) t name1       ls 1, \
+#      folder2."cpuTime.dat" u 1:(1e6*$4/$2) t name2       ls 2, \
+#      folder3."cpuTime.dat" u 1:(1e6*$4/$2) t name3       ls 3, \
+#      folder4."cpuTime.dat" u 1:(1e6*$4/$2) t name4       ls 4
+# 	 
+# set title "Reset"
+# plot folder1."cpuTime.dat" u 1:(1e6*$5/$2) t name1       ls 1, \
+#      folder2."cpuTime.dat" u 1:(1e6*$5/$2) t name2       ls 2, \
+#      folder3."cpuTime.dat" u 1:(1e6*$5/$2) t name3       ls 3, \
+#      folder4."cpuTime.dat" u 1:(1e6*$5/$2) t name4       ls 4
+# 	 
+# set title "Collisions"
+# plot folder1."cpuTime.dat" u 1:(1e6*$6/$2) t name1       ls 1, \
+#      folder2."cpuTime.dat" u 1:(1e6*$6/$2) t name2       ls 2, \
+#      folder3."cpuTime.dat" u 1:(1e6*$6/$2) t name3       ls 3, \
+#      folder4."cpuTime.dat" u 1:(1e6*$6/$2) t name4       ls 4
+# 	 
+# set title "Weighting"
+# plot folder1."cpuTime.dat" u 1:(1e6*$7/$2) t name1       ls 1, \
+#      folder2."cpuTime.dat" u 1:(1e6*$7/$2) t name2       ls 2, \
+#      folder3."cpuTime.dat" u 1:(1e6*$7/$2) t name3       ls 3, \
+#      folder4."cpuTime.dat" u 1:(1e6*$7/$2) t name4       ls 4
+# 
+# unset multiplot
diff --git a/src/DSMC_Neutrals.f95 b/src/DSMC_Neutrals.f95
index 310f393..a20d8cf 100644
--- a/src/DSMC_Neutrals.f95
+++ b/src/DSMC_Neutrals.f95
@@ -17,7 +17,6 @@ PROGRAM DSMC_Neutrals
   INTEGER:: t, n, e, i
   CHARACTER(200):: arg1
   CHARACTER(:), ALLOCATABLE:: inputFile
-  REAL(8), EXTERNAL::omp_get_wtime
 
 
   !Gets the input file
@@ -45,6 +44,7 @@ PROGRAM DSMC_Neutrals
   DO t = 1, tmax
     !Insert new particles
     !$OMP SINGLE
+    n_part_new = 0
     tStep = omp_get_wtime()
     tPush = omp_get_wtime()
     DO i=1, nInject
@@ -62,14 +62,6 @@ PROGRAM DSMC_Neutrals
 
     !$OMP SINGLE
     tPush = omp_get_wtime() - tPush
-    tReset = omp_get_wtime()
-    !$OMP END SINGLE
-
-    !Reset particles
-    CALL resetParticles(part_inj, part_old)
-
-    !$OMP SINGLE
-    tReset = omp_get_wtime() - tReset
     tColl = omp_get_wtime()
     !$OMP END SINGLE
 
@@ -82,6 +74,14 @@ PROGRAM DSMC_Neutrals
 
     !$OMP SINGLE
     tColl = omp_get_wtime() - tColl
+    tReset = omp_get_wtime()
+    !$OMP END SINGLE
+
+    !Reset particles
+    CALL resetParticles(part_inj, part_old)
+
+    !$OMP SINGLE
+    tReset = omp_get_wtime() - tReset
 
     !Weight
     tWeight = omp_get_wtime()
@@ -93,12 +93,14 @@ PROGRAM DSMC_Neutrals
     tWeight = omp_get_wtime() - tWeight
     tStep = omp_get_wtime() - tStep
     !Output data
+    !TODO: Move to subroutine
     counterOutput=counterOutput+1
     IF (counterOutput>=triggerOutput .OR. t == tmax) THEN
       counterOutput=0
       CALL mesh%printOutput(t)
       CALL printTime(t)
       CALL mesh%printColl(t)
+      !TODO: Move to subroutine
       PRINT *, "t/tmax: ", t, "/", tmax
       PRINT *, "Particles: ", n_part_old
       WRITE (*, "(2(5X,A26))") "total time (ms)", "avg t/particle (ns)"
diff --git a/src/modules/makefile b/src/modules/makefile
index d8ee4ae..0cf8c24 100644
--- a/src/modules/makefile
+++ b/src/modules/makefile
@@ -8,43 +8,43 @@ OBJS = moduleCaseParam.o moduleCompTime.o moduleList.o\
 all: $(OBJS)
 
 moduleCollisions.o: moduleTable.o moduleSpecies.o moduleRefParam.o moduleConstParam.o moduleMeshCyl.f95
-	$(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@
+	$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
 
 moduleInput.o: moduleParallel.o moduleRefParam.o moduleCaseParam.o moduleSolver.o moduleInject.o moduleBoundary.o moduleMesh.o moduleMeshCylRead.o moduleErrors.o moduleSpecies.o moduleInput.f95
-	$(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@
+	$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
 
 moduleInject.o: moduleSpecies.o moduleSolver.o moduleMesh.o moduleMeshCyl.o moduleInject.f95
-	$(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@
+	$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
 
 moduleList.o: moduleSpecies.o moduleErrors.o moduleList.f95
-	$(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@
+	$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
 
 moduleMesh.o: moduleOutput.o moduleList.o moduleSpecies.o moduleMesh.f95
-	$(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@
+	$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
 
 moduleMeshCyl.o: moduleRefParam.o moduleCollisions.o moduleOutput.o moduleMesh.o moduleMeshCyl.f95
-	$(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@
+	$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
 
 moduleMeshCylBoundary.o: moduleMeshCyl.o moduleMeshCylBoundary.f95
-	$(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@
+	$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
 
 moduleMeshCylRead.o: moduleBoundary.o moduleMeshCyl.o moduleMeshCylBoundary.o moduleMeshCylRead.f95
-	$(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@
+	$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
 
 moduleOutput.o: moduleSpecies.o moduleRefParam.o moduleOutput.f95
-	$(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@
+	$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
 
 moduleRefParam.o: moduleConstParam.o moduleRefParam.f95
-	$(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@
+	$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
 
 moduleSpecies.o: moduleCaseParam.o moduleSpecies.f95
-	$(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@
+	$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
 
 moduleSolver.o: moduleSpecies.o moduleRefParam.o moduleMesh.o moduleOutput.o moduleSolver.f95
-	$(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@
+	$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
 
 moduleTable.o: moduleErrors.o moduleTable.f95
-	$(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@
+	$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
 
 %.o: %.f95
 	$(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@
diff --git a/src/modules/moduleList.f95 b/src/modules/moduleList.f95
index 07a4013..2da0cc4 100644
--- a/src/modules/moduleList.f95
+++ b/src/modules/moduleList.f95
@@ -19,7 +19,6 @@ MODULE moduleList
   END TYPE listNode
 
   CONTAINS
-    !Finalize node
     !Adds element to list
     SUBROUTINE addToList(this,part)
       USE moduleSpecies
diff --git a/src/modules/moduleMesh.f95 b/src/modules/moduleMesh.f95
index 10ea618..e48e3c6 100644
--- a/src/modules/moduleMesh.f95
+++ b/src/modules/moduleMesh.f95
@@ -2,6 +2,7 @@
 MODULE moduleMesh
   USE moduleList
   USE moduleOutput
+  USE OMP_LIB
   IMPLICIT NONE
 
   !Parent of Node element
@@ -90,6 +91,8 @@ MODULE moduleMesh
     REAL(8):: volume = 0.D0
     !List of particles inside the volume
     TYPE(listNode):: listPart_in
+    !Lock indicator for listPart_in
+    INTEGER(KIND=OMP_LOCK_KIND):: lock
     !Number of collisions per volume
     INTEGER:: nColl = 0
     CONTAINS
@@ -131,7 +134,7 @@ MODULE moduleMesh
       IMPORT:: meshVol
       CLASS(meshVol), INTENT(inout):: self
       CLASS(meshVol), OPTIONAL, INTENT(in):: oldCell
-      CLASS(particle), INTENT(inout):: part
+      CLASS(particle), INTENT(inout), TARGET:: part
 
     END SUBROUTINE findCell_interface
 
diff --git a/src/modules/moduleMeshCyl.f95 b/src/modules/moduleMeshCyl.f95
index 61b9fb6..ad9c29e 100644
--- a/src/modules/moduleMeshCyl.f95
+++ b/src/modules/moduleMeshCyl.f95
@@ -175,6 +175,8 @@ MODULE moduleMeshCyl
 
       self%sigmaVrelMax = sigma_ref/L_ref**2
 
+      CALL OMP_INIT_LOCK(self%lock)
+
     END SUBROUTINE initVolQuadCyl
 
     FUNCTION fpsi(xi1,xi2) RESULT(psi) 
@@ -433,7 +435,7 @@ MODULE moduleMeshCyl
 
       CLASS(meshVolCylQuad), INTENT(inout):: self
       CLASS(meshVol), OPTIONAL, INTENT(in):: oldCell
-      CLASS(particle), INTENT(inout):: part
+      CLASS(particle), INTENT(inout), TARGET:: part
       REAL(8):: xLog(1:2)
       REAL(8):: xLogArray(1:4)
       CLASS(*), POINTER:: nextElement
@@ -445,9 +447,15 @@ MODULE moduleMeshCyl
         IF (PRESENT(oldCell)) THEN
           !If oldCell, recalculate particle weight, as particle has entered a new cell.
           part%weight = part%weight*oldCell%volume/self%volume
+
         END IF
         part%e_p  = self%n
         part%xLog = xLog
+        !Assign particle to listPart_in
+        CALL OMP_SET_LOCK(self%lock)
+        CALL self%listPart_in%add(part)
+        CALL OMP_UNSET_LOCK(self%lock)
+
       ELSE
         !If not, searches for a neighbour and repeats the process.
         xLogArray = (/ -xLog(2), xLog(1), xLog(2), -xLog(1) /)
@@ -514,8 +522,6 @@ MODULE moduleMeshCyl
 
 
     SUBROUTINE collision2DCyl(self)
-      USE moduleRefParam
-      USE moduleConstParam
       USE moduleCollisions
       USE moduleSpecies
       USE moduleList
@@ -549,8 +555,11 @@ MODULE moduleMeshCyl
 
         END DO
 
-        !Update maximum cross section times vrelative per each collision
-        IF (sigmaVrelMaxNew > self%sigmaVrelMax) self%sigmaVrelMax = sigmaVrelMaxNew
+        !Update maximum cross section*v_rel per each collision
+        IF (sigmaVrelMaxNew > self%sigmaVrelMax) THEN
+          self%sigmaVrelMax = sigmaVrelMaxNew
+
+        END IF
           
       END DO
 
diff --git a/src/modules/moduleSolver.f95 b/src/modules/moduleSolver.f95
index bf23dc5..8e01835 100644
--- a/src/modules/moduleSolver.f95
+++ b/src/modules/moduleSolver.f95
@@ -1,7 +1,6 @@
 MODULE moduleSolver
 
   CONTAINS
-
     SUBROUTINE scatterGrid(meshIn, partArray)
       USE moduleSpecies
       USE moduleMesh
@@ -53,15 +52,14 @@ MODULE moduleSolver
       part_temp%v(3) = -sin_alpha*v_p_oh_star(2)+cos_alpha*v_p_oh_star(3)
       part_temp%pt  = part%pt
       part_temp%e_p = part%e_p
-      !Assign cell to particle
+      !Copy temporal particle to particle
       part=part_temp
 
     END SUBROUTINE push
 
     SUBROUTINE resetParticles(partInj, partOld)
-      USE moduleSpecies
-      USE moduleList
-      USE moduleMesh
+      USE moduleSpecies, ONLY: particle, &
+                               n_part_old, n_part_new, nPartInj
       IMPLICIT NONE
 
       TYPE(particle), INTENT(in), ALLOCATABLE:: partInj(:)
@@ -88,16 +86,15 @@ MODULE moduleSolver
       !$OMP BARRIER
 
       !$OMP SINGLE
-      CALL MOVE_ALLOC(part_old, partTemp)
+      CALL MOVE_ALLOC(partOld, partTemp)
       n_part_new = n_inj_in + n_old_in
       ALLOCATE(partOld(1:n_part_new))
 
       nn = 0
       DO n = 1, nPartInj
-        IF (part_inj(n)%n_in) THEN
+        IF (partInj(n)%n_in) THEN
           nn = nn + 1
-          part_old(nn) = part_inj(n)
-          CALL mesh%vols(part_old(nn)%e_p)%obj%listPart_in%add(part_old(nn))
+          partOld(nn) = partInj(n)
 
         END IF
 
@@ -107,7 +104,6 @@ MODULE moduleSolver
         IF (partTemp(n)%n_in) THEN
           nn = nn + 1
           partOld(nn) = partTemp(n)
-          CALL mesh%vols(partOld(nn)%e_p)%obj%listPart_in%add(partOld(nn))
 
         END IF