diff --git a/runs/cylFlow/compCPUTime.gp b/runs/cylFlow/compCPUTime.gp
index 6ec49f8..2d657fc 100644
--- a/runs/cylFlow/compCPUTime.gp
+++ b/runs/cylFlow/compCPUTime.gp
@@ -1,4 +1,4 @@
-
+## Total time per iteration
 set terminal qt enhanced 1 persist size 1600, 1000 font "Times ,10"
 
 set style line 1  pt 4              lc rgb "#B50427" #Squares   red
@@ -8,14 +8,18 @@ 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
 
-name1   = "Base Case"
-folder1 = "base_case/"
-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/"
+#Name and folder 1 for comparison (include / at the end of the folder)
+name1   = ""
+folder1 = ""
+#Name and folder 2 for comparison (include / at the end of the folder)
+name2   = ""
+folder2 = ""
+#Name and folder 3 for comparison (include / at the end of the folder)
+name3   = ""
+folder3 = ""
+#Name and folder 4 for comparison (include / at the end of the folder)
+name4   = ""
+folder4 = ""
 
 set key box opaque
 set pointsize 1.5
@@ -50,7 +54,7 @@ plot folder1."cpuTime.dat" u 1:(1e3*$6) t name1       ls 1, \
      folder3."cpuTime.dat" u 1:(1e3*$6) t name3       ls 3, \
      folder4."cpuTime.dat" u 1:(1e3*$6) t name4       ls 4
 	 
-set title "Weighting"
+set title "Scattering"
 plot folder1."cpuTime.dat" u 1:(1e3*$7) t name1       ls 1, \
      folder2."cpuTime.dat" u 1:(1e3*$7) t name2       ls 2, \
      folder3."cpuTime.dat" u 1:(1e3*$7) t name3       ls 3, \
@@ -64,7 +68,7 @@ plot folder1."cpuTime.dat" u 1:($2) t name1       ls 1, \
 
 unset multiplot
 	 
-
+## Time per particle
 # set terminal qt enhanced 2 persist size 1600, 1000 font "Times ,10"
 # 
 # set style line 1  pt 4              lc rgb "#B50427" #Squares   red
@@ -107,7 +111,7 @@ unset multiplot
 #      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"
+# set title "Scattering"
 # 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, \
diff --git a/runs/cylFlow/input.json b/runs/cylFlow/input.json
index e50d719..1b56a73 100644
--- a/runs/cylFlow/input.json
+++ b/runs/cylFlow/input.json
@@ -3,7 +3,7 @@
     "path": "./runs/cylFlow/",
 	"trigger": 10,
 	"cpuTime": true,
-	"numColl": false
+	"numColl": true
   },
   "geometry": {
     "type":     "2DCyl",
diff --git a/src/DSMC_Neutrals.f95 b/src/DSMC_Neutrals.f95
index a20d8cf..50f4f29 100644
--- a/src/DSMC_Neutrals.f95
+++ b/src/DSMC_Neutrals.f95
@@ -1,24 +1,22 @@
 !  PPartiC neutral DSMC main program.
 PROGRAM DSMC_Neutrals
-  USE moduleMesh
-  USE moduleCompTime
-  USE moduleSolver
-  USE moduleSpecies
-  USE moduleInject
   USE moduleInput
   USE moduleErrors
-  USE moduleList
+  USE OMP_LIB
+  USE moduleInject
+  USE moduleSolver
   USE moduleOutput
-  USE moduleCaseParam
-  USE moduleParallel
+  USE moduleCompTime
+  USE moduleMesh
   IMPLICIT NONE
 
   ! t = time step; n = number of particle; e = number of element; i = number of inject
-  INTEGER:: t, n, e, i
+  INTEGER:: t = 0
   CHARACTER(200):: arg1
   CHARACTER(:), ALLOCATABLE:: inputFile
 
 
+  tStep = omp_get_wtime()
   !Gets the input file
   CALL getarg(1, arg1)
   inputFile = TRIM(arg1)
@@ -28,37 +26,24 @@ PROGRAM DSMC_Neutrals
   !Reads the json configuration file
   CALL readConfig(inputFile)
 
-  t = 0
+  tStep = omp_get_wtime() - tStep
+
   !Output initial state
-  CALL mesh%printOutput(t)
-  CALL printTime(t, .TRUE.)
-  CALL mesh%printColl(t)
-  PRINT *, "t/tmax: ", t, "/", tmax
-  PRINT *, "Particles: ", n_part_old
-
-  !Starts threads for OpenMP parallelization
-  !TODO: make this an input parameter
-  CALL OMP_SET_NUM_THREADS(openMP%nThreads)
-
+  CALL doOutput(t)
+  CALL verboseError('Starting main loop...')
   !$OMP PARALLEL PRIVATE(t) DEFAULT(SHARED) 
   DO t = 1, tmax
-    !Insert new particles
+    !Insert new particles and push them
     !$OMP SINGLE
-    n_part_new = 0
     tStep = omp_get_wtime()
     tPush = omp_get_wtime()
-    DO i=1, nInject
-      CALL inject(i)%addParticles()
-    END DO
+
+    CALL doInjects()
+
     !$OMP END SINGLE NOWAIT
 
-    !$OMP DO
-    DO n=1, n_part_old
-      CALL push(part_old(n))
-      CALL mesh%vols(part_old(n)%e_p)%obj%findCell(part_old(n))
-
-    END DO
-    !$OMP END DO
+    !Push old particles
+    CALL doPushes()
 
     !$OMP SINGLE
     tPush = omp_get_wtime() - tPush
@@ -66,11 +51,7 @@ PROGRAM DSMC_Neutrals
     !$OMP END SINGLE
 
     !Collisions
-    !$OMP DO
-    DO e=1, mesh%numVols
-      CALL mesh%vols(e)%obj%collision()
-    END DO
-    !$OMP END DO
+    CALL doCollisions()
 
     !$OMP SINGLE
     tColl = omp_get_wtime() - tColl
@@ -78,7 +59,7 @@ PROGRAM DSMC_Neutrals
     !$OMP END SINGLE
 
     !Reset particles
-    CALL resetParticles(part_inj, part_old)
+    CALL doReset()
 
     !$OMP SINGLE
     tReset = omp_get_wtime() - tReset
@@ -87,26 +68,13 @@ PROGRAM DSMC_Neutrals
     tWeight = omp_get_wtime()
     !$OMP END SINGLE
 
-    CALL scatterGrid(mesh, part_old)
+    CALL doScatter()
 
     !$OMP SINGLE
     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)"
-      WRITE (*, "(2(F31.3))")  1.D3*tStep,        1.D9*(tPush+tReset+tColl+tWeight)/DBLE(n_part_old)
-
-    END IF
+    CALL doOutput(t)
     !$OMP END SINGLE
 
   END DO
diff --git a/src/modules/moduleErrors.f95 b/src/modules/moduleErrors.f95
index 0a95db8..1d61142 100644
--- a/src/modules/moduleErrors.f95
+++ b/src/modules/moduleErrors.f95
@@ -37,6 +37,7 @@ MODULE moduleErrors
       errorMsg  = msg
 
       WRITE (*, '(A)') errorMsg
+      WRITE (*, *)
 
     END SUBROUTINE verboseError
 
diff --git a/src/modules/moduleInject.f95 b/src/modules/moduleInject.f95
index e0fea6d..fa23b02 100644
--- a/src/modules/moduleInject.f95
+++ b/src/modules/moduleInject.f95
@@ -23,6 +23,18 @@ MODULE moduleInject
   TYPE(injectGeneric), ALLOCATABLE:: inject(:)
 
   CONTAINS
+    !Injection of particles
+    SUBROUTINE doInjects()
+      IMPLICIT NONE
+
+      INTEGER:: i
+
+      DO i=1, nInject
+        CALL inject(i)%addParticles()
+      END DO
+
+    END SUBROUTINE doInjects
+
     SUBROUTINE initInject(self, i, v, n, T, flow, pt, physicalSurface)
       USE moduleMesh
       USE moduleMeshCyl
@@ -109,7 +121,7 @@ MODULE moduleInject
       INTEGER:: n1 = 0, n2 = 0
       !Edge nodes coordinates
       REAL(8):: p1(1:3) = 0.D0, p2(1:3) = 0.D0
-      INTEGER:: nMin, nMax !Min and Max index in part_inj array
+      INTEGER:: nMin, nMax !Min and Max index in partInj array
       INTEGER:: n
 
       vVec = self%v * self%n
@@ -120,10 +132,10 @@ MODULE moduleInject
       nMin = 1
       nMax = self%nParticles
       !Assign particle type
-      part_inj(nMin:nMax)%pt = self%pt
+      partInj(nMin:nMax)%pt = self%pt
       !Assign weight to particle.
-      part_inj(nMin:nMax)%weight = species(self%pt)%obj%weight
-      part_inj(nMin:nMax)%n_in = .TRUE.
+      partInj(nMin:nMax)%weight = species(self%pt)%obj%weight
+      partInj(nMin:nMax)%n_in = .TRUE.
       DO n = nMin, nMax
         !Select edge randomly from which inject particle
         randomX = RAND()*self%sumWeight
@@ -157,29 +169,29 @@ MODULE moduleInject
 
         !Volume associated to the edge:
         IF (DOT_PRODUCT(self%n, mesh%edges(randomEdge)%obj%normal) >= 0.D0) THEN
-          part_inj(n)%e_p = mesh%edges(randomEdge)%obj%e1%n
+          partInj(n)%e_p = mesh%edges(randomEdge)%obj%e1%n
         
         ELSE
-          part_inj(n)%e_p = mesh%edges(randomEdge)%obj%e2%n
+          partInj(n)%e_p = mesh%edges(randomEdge)%obj%e2%n
 
         END IF
 
-        part_inj(n)%v = (/ vBC(vVec(1), vTh(1)), &
+        partInj(n)%v = (/ vBC(vVec(1), vTh(1)), &
                            vBC(vVec(2), vTh(2)), &
                            vBC(vVec(3), vTh(3)) /)
 
         !Random position in edge
         !TODO: Use edge coordinates and transformations for this process
         randomPos = RAND()
-        part_inj(n)%r(1) = p1(1) + randomPos*(p2(1) - p1(1))
-        part_inj(n)%r(2) = p1(2) + randomPos*(p2(2) - p1(2))
-        part_inj(n)%r(3) = p1(3) + randomPos*(p2(3) - p1(3))
+        partInj(n)%r(1) = p1(1) + randomPos*(p2(1) - p1(1))
+        partInj(n)%r(2) = p1(2) + randomPos*(p2(2) - p1(2))
+        partInj(n)%r(3) = p1(3) + randomPos*(p2(3) - p1(3))
 
 
         !Push new particle
-        CALL push(part_inj(n))
+        CALL push(partInj(n))
         !Assign cell to new particle
-        CALL mesh%vols(part_inj(n)%e_p)%obj%findCell(part_inj(n))
+        CALL mesh%vols(partInj(n)%e_p)%obj%findCell(partInj(n))
 
       END DO
 
diff --git a/src/modules/moduleInput.f95 b/src/modules/moduleInput.f95
index fe0d096..dc9cf61 100644
--- a/src/modules/moduleInput.f95
+++ b/src/modules/moduleInput.f95
@@ -19,6 +19,7 @@ MODULE moduleInput
       CALL config%initialize()
 
       !Loads the config file
+      CALL verboseError('Loading input file...')
       CALL config%load(filename = inputFile)
 
       !Reads reference parameters
@@ -40,12 +41,15 @@ MODULE moduleInput
       CALL readBoundary(config)
 
       !Read Geometry
+      CALL verboseError('Reading Geometry...')
       CALL readGeometry(config)
 
       !Read injection of particles
+      CALL verboseError('Reading Interactions between species...')
       CALL readInject(config)
 
       !Read parallel parameters
+      CALL verboseError('Reading Parallel configuration...')
       CALL readParallel(config)
 
     END SUBROUTINE readConfig
@@ -196,7 +200,7 @@ MODULE moduleInput
 
       !Set number of particles to 0 for init state
       !TODO: In a future, this should include the particles from init states
-      n_part_old = 0
+      nPartOld = 0
 
     END SUBROUTINE readSpecies
 
@@ -363,7 +367,7 @@ MODULE moduleInput
 
       !Allocate array for injected particles
       IF (nPartInj > 0) THEN
-        ALLOCATE(part_inj(1:nPartInj))
+        ALLOCATE(partInj(1:nPartInj))
 
       END IF
       
@@ -390,6 +394,8 @@ MODULE moduleInput
 
       END IF
 
+      CALL initParallel()
+
     END SUBROUTINE readParallel
 
 
diff --git a/src/modules/moduleMeshCylRead.f95 b/src/modules/moduleMeshCylRead.f95
index c3e3b69..df6e51b 100644
--- a/src/modules/moduleMeshCylRead.f95
+++ b/src/modules/moduleMeshCylRead.f95
@@ -317,7 +317,7 @@ MODULE moduleMeshCylRead
       DO i = 1, nSpecies
         WRITE(tstring, '(I6.6)') t
         fileName='OUTPUT_' // tstring// '_' // species(i)%obj%name // '.msh'
-        PRINT *, "Creado archivo de datos:", fileName
+        WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
         OPEN (60, file = path // folder // '/' //  fileName)
         WRITE(60, "(A)") '$MeshFormat'
         WRITE(60, "(A)") '2.2 0 8'
@@ -400,7 +400,7 @@ MODULE moduleMeshCylRead
         WRITE(tstring, '(I6.6)') t
 
         fileName='OUTPUT_' // tstring// '_Collisions.msh'
-        PRINT *, "Creado archivo de datos:", fileName
+        WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
         OPEN (60, file = path // folder // '/' //  fileName)
         WRITE(60, "(A)") '$MeshFormat'
         WRITE(60, "(A)") '2.2 0 8'
diff --git a/src/modules/moduleOutput.f95 b/src/modules/moduleOutput.f95
index b9f3ae2..dfe3808 100644
--- a/src/modules/moduleOutput.f95
+++ b/src/modules/moduleOutput.f95
@@ -20,7 +20,6 @@ MODULE moduleOutput
   LOGICAL:: collOutput = .FALSE.
 
   CONTAINS
-
     FUNCTION outerProduct(a,b) RESULT(s)
       IMPLICIT NONE
 
@@ -89,33 +88,34 @@ MODULE moduleOutput
 
       INTEGER, INTENT(in):: t
       LOGICAL, INTENT(in), OPTIONAL:: first
-      CHARACTER(:), ALLOCATABLE:: filename
+      CHARACTER(:), ALLOCATABLE:: fileName
 
-      filename = 'cpuTime.dat'
+      fileName = 'cpuTime.dat'
 
       IF (timeOutput) THEN
         IF (PRESENT(first)) THEN
           IF (first) THEN
-            OPEN(20, file = path // folder // '/' // filename, action = 'write')
+            OPEN(20, file = path // folder // '/' // fileName, action = 'write')
             WRITE(20, "(A1, 8X, A1, 9X, A1, 5(A20))") "#","t","n","total","push","reset","collision","weighting"
+            WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
 
           ELSE
 
           END IF
-            OPEN(20, file = path // folder // '/' // filename, position = 'append', action = 'write')
+            OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write')
 
         ELSE
-          OPEN(20, file = path // folder // '/' // filename, position = 'append', action = 'write')
+          OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write')
 
         END IF
 
-        WRITE (20, "(I10, I10, 5(ES20.6E3))") t, n_part_old, tStep, tPush, tReset, tColl, tWeight
+        WRITE (20, "(I10, I10, 5(ES20.6E3))") t, nPartOld, tStep, tPush, tReset, tColl, tWeight
 
         CLOSE(20)
 
       END IF
 
-    END SUBROUTINE
+    END SUBROUTINE printTime
 
 END MODULE moduleOutput
 
diff --git a/src/modules/moduleParallel.f95 b/src/modules/moduleParallel.f95
index 9c83385..564724c 100644
--- a/src/modules/moduleParallel.f95
+++ b/src/modules/moduleParallel.f95
@@ -1,4 +1,5 @@
 MODULE moduleParallel
+  IMPLICIT NONE
   
   TYPE openMP_type
     INTEGER:: nThreads
@@ -7,4 +8,13 @@ MODULE moduleParallel
 
   TYPE(openMP_type):: openMP
 
+  CONTAINS
+    SUBROUTINE initParallel()
+      IMPLICIT NONE
+
+      !Starts threads for OpenMP parallelization
+      CALL OMP_SET_NUM_THREADS(openMP%nThreads)
+
+    END SUBROUTINE initParallel
+
 END MODULE moduleParallel
diff --git a/src/modules/moduleSolver.f95 b/src/modules/moduleSolver.f95
index 8e01835..59ed245 100644
--- a/src/modules/moduleSolver.f95
+++ b/src/modules/moduleSolver.f95
@@ -1,24 +1,25 @@
 MODULE moduleSolver
 
   CONTAINS
-    SUBROUTINE scatterGrid(meshIn, partArray)
+    SUBROUTINE doPushes()
       USE moduleSpecies
       USE moduleMesh
+      IMPLICIT NONE
 
-      CLASS(meshGeneric), INTENT(in):: meshIn
-      TYPE(particle), INTENT(in):: partArray(:)
       INTEGER:: n
 
-      !Loops over the particles to scatter them
       !$OMP DO
-      DO n=1, n_part_old
-        CALL meshIn%vols(partArray(n)%e_p)%obj%scatter(partArray(n))
+      DO n=1, nPartOld
+        CALL push(partOld(n))
+        CALL mesh%vols(partOld(n)%e_p)%obj%findCell(partOld(n))
 
       END DO
       !$OMP END DO
 
-    END SUBROUTINE scatterGrid
+    END SUBROUTINE doPushes
 
+    !Push one particle. Boris pusher for 2D Cyl
+    !TODO: make general for any pusher
     PURE SUBROUTINE push(part)
       USE moduleSpecies
       USE moduleMesh
@@ -57,28 +58,58 @@ MODULE moduleSolver
 
     END SUBROUTINE push
 
-    SUBROUTINE resetParticles(partInj, partOld)
-      USE moduleSpecies, ONLY: particle, &
-                               n_part_old, n_part_new, nPartInj
+    !Do the collisions in all the cells
+    SUBROUTINE doCollisions()
+      USE moduleMesh
+      IMPLICIT NONE
+
+      INTEGER:: e
+
+      !$OMP DO
+      DO e=1, mesh%numVols
+        CALL mesh%vols(e)%obj%collision()
+      END DO
+      !$OMP END DO
+
+    END SUBROUTINE doCollisions
+    
+    !Scatter particles in the grid
+    SUBROUTINE doScatter
+      USE moduleSpecies
+      USE moduleMesh
+
+      INTEGER:: n
+
+      !Loops over the particles to scatter them
+      !$OMP DO
+      DO n=1, nPartOld
+        CALL mesh%vols(partOld(n)%e_p)%obj%scatter(partOld(n))
+
+      END DO
+      !$OMP END DO
+
+    END SUBROUTINE doScatter
+
+    SUBROUTINE doReset()
+      USE moduleSpecies
       IMPLICIT NONE
 
-      TYPE(particle), INTENT(in), ALLOCATABLE:: partInj(:)
-      TYPE(particle), INTENT(inout), ALLOCATABLE:: partOld(:)
       INTEGER:: nn, n
-      INTEGER, SAVE:: n_inj_in, n_old_in
+      INTEGER:: nPartNew
+      INTEGER, SAVE:: nInjIn, nOldIn
       TYPE(particle), ALLOCATABLE:: partTemp(:)
 
       !$OMP SECTIONS
       !$OMP SECTION
-      n_inj_in = 0
+      nInjIn = 0
       IF (ALLOCATED(partInj)) THEN
-        n_inj_in = COUNT(partInj%n_in)
+        nInjIn = COUNT(partInj%n_in)
 
       END IF
       !$OMP SECTION
-      n_old_in = 0
+      nOldIn = 0
       IF (ALLOCATED(partOld)) THEN
-        n_old_in = COUNT(partOld%n_in)
+        nOldIn = COUNT(partOld%n_in)
 
       END IF
       !$OMP END SECTIONS
@@ -87,8 +118,8 @@ MODULE moduleSolver
 
       !$OMP SINGLE
       CALL MOVE_ALLOC(partOld, partTemp)
-      n_part_new = n_inj_in + n_old_in
-      ALLOCATE(partOld(1:n_part_new))
+      nPartNew = nInjIn + nOldIn
+      ALLOCATE(partOld(1:nPartNew))
 
       nn = 0
       DO n = 1, nPartInj
@@ -100,7 +131,7 @@ MODULE moduleSolver
 
       END DO
 
-      DO n = 1, n_part_old
+      DO n = 1, nPartOld
         IF (partTemp(n)%n_in) THEN
           nn = nn + 1
           partOld(nn) = partTemp(n)
@@ -109,10 +140,50 @@ MODULE moduleSolver
 
       END DO
 
-      n_part_old = n_part_new
+      nPartOld = nPartNew
       !$OMP END SINGLE
 
-    END SUBROUTINE resetParticles
+    END SUBROUTINE doReset
+
+    SUBROUTINE doOutput(t)
+      USE moduleMesh
+      USE moduleOutput
+      USE moduleSpecies
+      USE moduleCompTime
+      IMPLICIT NONE
+
+      INTEGER, INTENT(in):: t
+
+      counterOutput=counterOutput+1
+      IF (counterOutput>=triggerOutput .OR. &
+          t == tmax .OR. t == 0) THEN
+
+        !Resets output counter
+        counterOutput=0
+
+        CALL mesh%printOutput(t)
+        CALL printTime(t, t == 0)
+        CALL mesh%printColl(t)
+        WRITE(*, "(5X,A21,I10,A1,I8)") "t/tmax: ",    t,       "/", tmax
+        WRITE(*, "(5X,A21,I10)")       "Particles: ", nPartOld
+        IF (t == 0) THEN
+          WRITE(*, "(5X,A21,F8.1,A2)")   "     init time: ",    1.D3*tStep, "ms"
+
+        ELSE
+          WRITE(*, "(5X,A21,F8.1,A2)")   "iteration time: ",    1.D3*tStep, "ms"
+
+        END IF
+
+        IF (nPartOld > 0) THEN
+          WRITE(*, "(5X,A21,F8.1,A2)") "avg t/particle: ", 1.D9*(tPush+tReset+tColl+tWeight)/DBLE(nPartOld), "ns"
+
+        END IF
+        WRITE(*,*)
+
+      END IF
+
+    END SUBROUTINE doOutput
+
 
 END MODULE moduleSolver
 
diff --git a/src/modules/moduleSpecies.f95 b/src/modules/moduleSpecies.f95
index 6ca4b87..db2de87 100644
--- a/src/modules/moduleSpecies.f95
+++ b/src/modules/moduleSpecies.f95
@@ -37,12 +37,12 @@ MODULE moduleSpecies
 
   END TYPE particle
 
-  INTEGER:: n_part_old, n_part_new
+  !Number of old particles
+  INTEGER:: nPartOld
   INTEGER:: nPartInj
-  !Arrays that define the particles
-  TYPE(particle), ALLOCATABLE, DIMENSION(:), TARGET:: part_old
-  TYPE(particle), ALLOCATABLE, DIMENSION(:):: part_new
-  TYPE(particle), ALLOCATABLE, DIMENSION(:):: part_inj !array of inject particles
+  !Arrays that contain the particles
+  TYPE(particle), ALLOCATABLE, DIMENSION(:), TARGET:: partOld
+  TYPE(particle), ALLOCATABLE, DIMENSION(:), TARGET:: partInj !array of inject particles
 
   CONTAINS
     FUNCTION speciesName2Index(speciesName) RESULT(pt)