The integer part%sp that referenced the species index has been

substituted for a pointer to the species.

src/modules/mesh/1DCart/moduleMesh1DCart.f90

@@ -395,12 +395,12 @@ MODULE moduleMesh1DCart
       w_p = self%weight(part%xi)
       tensorS = outerProduct(part%v, part%v)
 
-      vertex => self%n1%output(part%sp)
+      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
 
-      vertex => self%n2%output(part%sp)
+      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

src/modules/mesh/1DRad/moduleMesh1DRad.f90

@@ -406,12 +406,12 @@ MODULE moduleMesh1DRad
       w_p = self%weight(part%xi)
       tensorS = outerProduct(part%v, part%v)
 
-      vertex => self%n1%output(part%sp)
+      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
 
-      vertex => self%n2%output(part%sp)
+      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

src/modules/mesh/2DCart/moduleMesh2DCart.f90

@@ -510,22 +510,22 @@ MODULE moduleMesh2DCart
       w_p = self%weight(part%xi)
       tensorS = outerProduct(part%v, part%v)
 
-      vertex => self%n1%output(part%sp)
+      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
 
-      vertex => self%n2%output(part%sp)
+      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
 
-      vertex => self%n3%output(part%sp)
+      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
 
-      vertex => self%n4%output(part%sp)
+      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
@@ -869,17 +869,17 @@ MODULE moduleMesh2DCart
       w_p = self%weight(part%xi)
       tensorS = outerProduct(part%v, part%v)
 
-      vertex => self%n1%output(part%sp)
+      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
 
-      vertex => self%n2%output(part%sp)
+      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
 
-      vertex => self%n3%output(part%sp)
+      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

src/modules/mesh/2DCyl/moduleMesh2DCyl.f90

@@ -531,22 +531,22 @@ MODULE moduleMesh2DCyl
       w_p = self%weight(part%xi)
       tensorS = outerProduct(part%v, part%v)
 
-      vertex => self%n1%output(part%sp)
+      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
 
-      vertex => self%n2%output(part%sp)
+      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
 
-      vertex => self%n3%output(part%sp)
+      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
 
-      vertex => self%n4%output(part%sp)
+      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
@@ -898,17 +898,17 @@ MODULE moduleMesh2DCyl
       w_p = self%weight(part%xi)
       tensorS = outerProduct(part%v, part%v)
 
-      vertex => self%n1%output(part%sp)
+      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
 
-      vertex => self%n2%output(part%sp)
+      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
 
-      vertex => self%n3%output(part%sp)
+      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

src/modules/mesh/3DCart/moduleMesh3DCart.f90

@@ -496,22 +496,22 @@ MODULE moduleMesh3DCart
       w_p = self%weight(part%xi)
       tensorS = outerProduct(part%v, part%v)
 
-      vertex => self%n1%output(part%sp)
+      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
 
-      vertex => self%n2%output(part%sp)
+      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
 
-      vertex => self%n3%output(part%sp)
+      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
 
-      vertex => self%n4%output(part%sp)
+      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

src/modules/mesh/moduleMesh.f90

@@ -375,7 +375,7 @@ MODULE moduleMesh
 
         CLASS IS (meshEdge)
           !Particle encountered a surface, apply boundary
-          CALL nextElement%fBoundary(part%sp)%apply(nextElement,part)
+          CALL nextElement%fBoundary(part%species%n)%apply(nextElement,part)
 
           !If particle is still inside the domain, call findCell
           IF (part%n_in) THEN 
@@ -439,7 +439,7 @@ MODULE moduleMesh
           part_i => partTemp(rnd)%part
           rnd = random(1, nPart)
           part_j => partTemp(rnd)%part
-          ij = interactionIndex(part_i%sp, part_j%sp)
+          ij = interactionIndex(part_i%species%n, part_j%species%n)
           sigmaVrelMaxNew = 0.D0
           DO k = 1, interactionMatrix(ij)%amount
             CALL interactionMatrix(ij)%collisions(k)%obj%collide(self%sigmaVrelMax, sigmaVrelMaxNew, part_i, part_j)

src/modules/mesh/moduleMeshBoundary.f90

@@ -91,7 +91,7 @@ MODULE moduleMeshBoundary
       INTEGER:: i
 
       !Modifies particle velocity according to wall temperature
-      SELECT TYPE(bound => edge%boundary%bTypes(part%sp)%obj)
+      SELECT TYPE(bound => edge%boundary%bTypes(part%species%n)%obj)
       TYPE IS(boundaryWallTemperature)
         DO i = 1, 3
           part%v(i) = part%v(i) + bound%vTh*randomMaxwellian()
@@ -123,9 +123,9 @@ MODULE moduleMeshBoundary
       TYPE(particle), POINTER:: newElectron
       TYPE(particle), POINTER:: newIon
 
-      SELECT TYPE(bound => edge%boundary%bTypes(part%sp)%obj)
+      SELECT TYPE(bound => edge%boundary%bTypes(part%species%n)%obj)
       TYPE IS(boundaryIonization)
-        mRel = (bound%m0*species(part%sp)%obj%m)*(bound%m0+species(part%sp)%obj%m)
+        mRel = (bound%m0*part%species%m)*(bound%m0+part%species%m)
         vRel = SUM(DABS(part%v-bound%v0))
         eRel = mRel*vRel**2*5.D-1
 
@@ -133,7 +133,7 @@ MODULE moduleMeshBoundary
           ionizationRate = part%weight*bound%n0*bound%crossSection%get(eRel)*vRel
 
           !Rounds the number of particles up
-          ionizationPair = NINT(ionizationRate*bound%effectiveTime/species(bound%sp)%obj%weight)
+          ionizationPair = NINT(ionizationRate*bound%effectiveTime/bound%species%weight)
 
           !Create the new pair of particles
           DO p = 1, ionizationPair
@@ -146,8 +146,8 @@ MODULE moduleMeshBoundary
             ALLOCATE(newElectron)
             ALLOCATE(newIon)
 
-            newElectron%sp = part%sp
+            newElectron%species => part%species
-            newIon%sp      = bound%sp
+            newIon%species      => bound%species
 
             newElectron%v = v0 + (1.D0 + bound%deltaV*v0/NORM2(v0))
             newIon%v      = v0
@@ -162,13 +162,13 @@ MODULE moduleMeshBoundary
             newIon%xi      = newElectron%xi
 
             newElectron%qm = part%qm
-            SELECT TYPE(spe => species(bound%sp)%obj)
+            SELECT TYPE(spe => bound%species)
             TYPE IS(speciesCharged)
               newIon%qm = spe%qm
 
             END SELECT
 
-            newElectron%weight = species(bound%sp)%obj%weight
+            newElectron%weight = bound%species%weight
             newIon%weight      = newElectron%weight
 
             newElectron%n_in = .TRUE.

src/modules/moduleBoundary.f90

@@ -1,5 +1,6 @@
 MODULE moduleBoundary
   USE moduleTable
+  USE moduleSpecies
 
   !Generic type for boundaries
   TYPE, PUBLIC:: boundaryGeneric
@@ -36,7 +37,7 @@ MODULE moduleBoundary
   !Ionization boundary
   TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryIonization
     REAL(8):: m0, n0, v0(1:3), vTh !Properties of background neutrals.
-    INTEGER:: sp !Ion species
+    CLASS(speciesGeneric), POINTER:: species !Ion species
     TYPE(table1D):: crossSection
     REAL(8):: effectiveTime
     REAL(8):: eThreshold
@@ -123,7 +124,7 @@ MODULE moduleBoundary
         boundary%n0 = n0 * Vol_ref
         boundary%v0 = v0 / v_ref
         boundary%vTh = DSQRT(kb*T0/m0)/v_ref
-        boundary%sp = speciesID
+        boundary%species => species(speciesID)%obj
         boundary%effectiveTime = effTime / ti_ref
         CALL boundary%crossSection%init(crossSection)
         CALL boundary%crossSection%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2)

src/modules/moduleCollisions.f90

@@ -199,8 +199,8 @@ MODULE moduleCollisions
       sigmaVrel = self%crossSec%get(eRel)*vRel
       sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel
       IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN
-        m_i = species(part_i%sp)%obj%m
+        m_i = part_i%species%m
-        m_j = species(part_j%sp)%obj%m
+        m_j = part_j%species%m
         !Applies the collision
         vCM = velocityCM(m_i, part_i%v, m_j, part_j%v)
         vp  = vRel*randomDirectionVHS()
@@ -292,11 +292,11 @@ MODULE moduleCollisions
         sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel
         IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN
           !Find which particle is the ionizing electron
-          IF    (part_i%sp == self%electron%sp) THEN
+          IF    (ASSOCIATED(part_i%species, self%electron)) THEN
             electron => part_i
             neutral  => part_j
 
-          ELSEIF(part_j%sp == self%electron%sp) THEN
+          ELSEIF(ASSOCIATED(part_j%species, self%electron)) THEN
             electron => part_j
             neutral  => part_i
 
@@ -314,7 +314,7 @@ MODULE moduleCollisions
 
           !Creates a new electron from ionization
           ALLOCATE(newElectron)
-          newElectron%sp = electron%sp
+          newElectron%species => electron%species
           newElectron%v       =  vp_n
           newElectron%r       =  neutral%r
           newElectron%xi      = neutral%xi
@@ -324,7 +324,7 @@ MODULE moduleCollisions
           newElectron%qm      = electron%qm
 
           !Ionize neutral particle
-          SELECT TYPE(sp => species(neutral%sp)%obj)
+          SELECT TYPE(sp => neutral%species)
           TYPE IS(speciesNeutral)
             CALL sp%ionize(neutral)
 
@@ -421,11 +421,11 @@ MODULE moduleCollisions
       sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel
       IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN
         !Find which particle is the ionizing electron
-        IF    (part_i%sp == self%electron%sp) THEN
+        IF    (ASSOCIATED(part_i%species, self%electron)) THEN
           electron => part_i
           ion      => part_j
 
-        ELSEIF(part_j%sp == self%electron%sp) THEN
+        ELSEIF(ASSOCIATED(part_j%species, self%electron)) THEN
           electron => part_j
           ion      => part_i
 
@@ -442,7 +442,7 @@ MODULE moduleCollisions
         electron%n_in = .FALSE.
 
         !Neutralize ion particle
-        SELECT TYPE(sp => species(ion%sp)%obj)
+        SELECT TYPE(sp => ion%species)
         TYPE IS(speciesCharged)
           CALL sp%neutralize(ion)
 
@@ -501,7 +501,7 @@ MODULE moduleCollisions
       sigmaVrel = self%crossSec%get(eRel)*vRel
       sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel
       IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN
-        SELECT TYPE(sp => species(part_i%sp)%obj)
+        SELECT TYPE(sp => part_i%species)
         TYPE IS (speciesNeutral)
           !Species i is neutral, ionize particle i
           CALL sp%ionize(part_i)
@@ -512,7 +512,7 @@ MODULE moduleCollisions
 
         END SELECT
 
-        SELECT TYPE(sp => species(part_j%sp)%obj)
+        SELECT TYPE(sp => part_j%species)
         TYPE IS (speciesNeutral)
           !Species j is neutral, ionize particle j
           CALL sp%ionize(part_j)

src/modules/moduleInject.f90

@@ -1,5 +1,6 @@
 !injection of particles
 MODULE moduleInject
+  USE moduleSpecies
 
   !Generic type for velocity distribution function
   TYPE, ABSTRACT:: velDistGeneric
@@ -51,7 +52,7 @@ MODULE moduleInject
     REAL(8):: T(1:3) !Temperature
     REAL(8):: n(1:3) !Direction of injection
     INTEGER:: nParticles !Number of particles to introduce each time step
-    INTEGER:: sp !Species of injection
+    CLASS(speciesGeneric), POINTER:: species !Species of injection
     INTEGER:: nEdges
     INTEGER, ALLOCATABLE:: edges(:) !Array with edges
     REAL(8), ALLOCATABLE:: cumWeight(:) !Array of cummulative probability
@@ -90,7 +91,7 @@ MODULE moduleInject
       self%vMod       =  v/v_ref
       self%n          =  n
       self%T          =  T/T_ref
-      self%sp  = sp
+      self%species    => species(sp)%obj
       SELECT CASE(units)
       CASE ("sccm")
         !Standard cubic centimeter per minute
@@ -224,14 +225,12 @@ MODULE moduleInject
       !$OMP SINGLE
       nMin = SUM(inject(1:(self%id-1))%nParticles) + 1
       nMax = nMin + self%nParticles - 1
-      !Assign particle type
-      partInj(nMin:nMax)%sp = self%sp
       !Assign weight to particle.
-      partInj(nMin:nMax)%weight = species(self%sp)%obj%weight
+      partInj(nMin:nMax)%weight = self%species%weight
       !Particle is considered to be outside the domain
       partInj(nMin:nMax)%n_in = .FALSE.
       !Assign charge/mass to charged particle.
-      SELECT TYPE(sp => species(self%sp)%obj)
+      SELECT TYPE(sp => self%species)
       TYPE IS(speciesCharged)
         partInj(nMin:nMax)%qm = sp%qm
 
@@ -257,6 +256,9 @@ MODULE moduleInject
 
         END IF
 
+        !Assign particle type
+        partInj(n)%species => self%species
+
         partInj(n)%v = (/ self%v(1)%obj%randomVel(), &
                           self%v(2)%obj%randomVel(), &
                           self%v(3)%obj%randomVel() /)
@@ -264,7 +266,7 @@ MODULE moduleInject
         !Obtain natural coordinates of particle in cell
         partInj(n)%xi = mesh%vols(partInj(n)%vol)%obj%phy2log(partInj(n)%r)
         !Push new particle with the minimum time step
-        CALL solver%pusher(self%sp)%pushParticle(partInj(n), tauMin)
+        CALL solver%pusher(self%species%n)%pushParticle(partInj(n), tauMin)
         !Assign cell to new particle
         CALL solver%updateParticleCell(partInj(n))
 

src/modules/moduleInput.f90

@@ -283,7 +283,7 @@ MODULE moduleInput
             !Allocate new particles
             DO p = 1, nNewPart
               ALLOCATE(partNew)
-              partNew%sp     = sp
+              partNew%species => species(sp)%obj
               partNew%v(1)    =  velocity(1) + vTh*randomMaxwellian()
               partNew%v(2)    =  velocity(2) + vTh*randomMaxwellian()
               partNew%v(3)    =  velocity(3) + vTh*randomMaxwellian()
@@ -442,7 +442,7 @@ MODULE moduleInput
         !Assign shared parameters for all species
         CALL config%get(object // '.name',   species(i)%obj%name,   found)
         CALL config%get(object // '.weight', species(i)%obj%weight, found)
-        species(i)%obj%sp = i
+        species(i)%obj%n = i
         species(i)%obj%m = mass
 
       END DO
@@ -894,7 +894,7 @@ MODULE moduleInput
       REAL(8):: v, T, m
 
       !Reads species mass
-      m = species(inj%sp)%obj%m
+      m = inj%species%m
       !Reads distribution functions for velocity
       DO i = 1, 3
         WRITE(istring, '(i2)') i

src/modules/moduleSolver.f90

@@ -154,7 +154,7 @@ MODULE moduleSolver
       !$OMP DO
       DO n=1, nPartOld
         !Select species type
-        sp = partOld(n)%sp
+        sp = partOld(n)%species%n
         !Checks if the species sp is update this iteration
         IF (solver%pusher(sp)%pushSpecies) THEN
           !Push particle
@@ -722,7 +722,7 @@ MODULE moduleSolver
 
         part%weight = part%weight * fractionVolume
 
-        fractionWeight = part%weight / species(part%sp)%obj%weight
+        fractionWeight = part%weight / part%species%weight
 
         IF (fractionWeight >= 2.D0) THEN
           nSplit = FLOOR(fractionWeight)

src/modules/moduleSpecies.f90

@@ -7,7 +7,7 @@ MODULE moduleSpecies
   TYPE, ABSTRACT:: speciesGeneric
     CHARACTER(:), ALLOCATABLE:: name
     REAL(8):: m=0.D0, weight=0.D0
-    INTEGER:: sp=0
+    INTEGER:: n=0
   END TYPE speciesGeneric
 
   TYPE, EXTENDS(speciesGeneric):: speciesNeutral
@@ -37,7 +37,7 @@ MODULE moduleSpecies
   TYPE particle
     REAL(8):: r(1:3) !Position
     REAL(8):: v(1:3) !Velocity
-    INTEGER:: sp !Particle species id
+    CLASS(speciesGeneric), POINTER:: species !Pointer to species associated with this particle
     INTEGER:: vol !Index of element in which the particle is located
     REAL(8):: xi(1:3) !Logical coordinates of particle in element e_p.
     LOGICAL:: n_in !Flag that indicates if a particle is in the domain
@@ -48,6 +48,7 @@ MODULE moduleSpecies
 
   !Number of old particles
   INTEGER:: nPartOld
+  !Number of injected particles
   INTEGER:: nPartInj
   !Arrays that contain the particles
   TYPE(particle), ALLOCATABLE, DIMENSION(:), TARGET:: partOld !array of particles from previous iteration
@@ -60,15 +61,18 @@ MODULE moduleSpecies
 
       CHARACTER(:), ALLOCATABLE:: speciesName
       INTEGER:: sp 
-      INTEGER:: n
+      INTEGER:: s
 
       sp = 0
-      DO n = 1, nSpecies
+      DO s = 1, nSpecies
-        IF (speciesName == species(n)%obj%name) THEN
+        IF (speciesName == species(s)%obj%name) THEN
-          sp = species(n)%obj%sp
+          sp = species(s)%obj%n
-          EXIT
+          EXIT !If a species is found, exit the loop
+
         END IF
+
       END DO
+
       !If no species is found, call a critical error
       IF (sp == 0) CALL criticalError('Species ' // speciesName // ' not found.', 'speciesName2Index')
 
@@ -83,7 +87,7 @@ MODULE moduleSpecies
       TYPE(particle), INTENT(inout):: part
 
       IF (ASSOCIATED(self%ion)) THEN
-        part%sp = self%ion%sp
+        part%species => self%ion
 
       ELSE
         CALL criticalError('No ion defined for species' // self%name, 'ionizeNeutral')
@@ -101,7 +105,7 @@ MODULE moduleSpecies
       TYPE(particle), INTENT(inout):: part
 
       IF (ASSOCIATED(self%ion)) THEN
-        part%sp = self%ion%sp
+        part%species => self%ion
 
       ELSE
         CALL criticalError('No ion defined for species' // self%name, 'ionizeCharged')
@@ -119,7 +123,7 @@ MODULE moduleSpecies
       TYPE(particle), INTENT(inout):: part
 
       IF (ASSOCIATED(self%neutral)) THEN
-        part%sp = self%neutral%sp
+        part%species => self%neutral
 
       ELSE
         CALL criticalError('No neutral defined for species' // self%name, 'neutralizeCharged')