Final implementation of ionization process by electron impact.

Possibility to input initial species distributions (density, velocity
    and temperature) via an input file for each species.

New moduleRandom includes function to generate random numbers in
different ways (still uses) the implicit RANDOM_NUMBER().

src/fpakc.f90

@@ -2,13 +2,11 @@
 PROGRAM fpakc
   USE moduleInput
   USE moduleErrors
-  USE OMP_LIB
   USE moduleInject
   USE moduleSolver
-  USE moduleOutput
   USE moduleCompTime
-  USE moduleCollisions
+  USE moduleCaseParam
-  USE moduleMesh
+  USE OMP_LIB
   IMPLICIT NONE
 
   ! t = time step
@@ -28,8 +26,17 @@ PROGRAM fpakc
   !Reads the json configuration file
   CALL readConfig(inputFile)
 
+  !$OMP PARALLEL DEFAULT(SHARED)
+  !$OMP SINGLE
+  CALL verboseError("Initial scatter of particles...")
+  !$OMP END SINGLE
+  CALL doScatter()
+
+  !$OMP SINGLE
   CALL verboseError("Calculating initial EM field...")
+  !$OMP END SINGLE
   CALL doEMField()
+  !$OMP END PARALLEL
 
   tStep = omp_get_wtime() - tStep
 

src/makefile

@@ -3,7 +3,7 @@ OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleCompTime.o $(OBJDIR)/moduleSolv
 			  	$(OBJDIR)/moduleErrors.o $(OBJDIR)/moduleList.o $(OBJDIR)/moduleOutput.o \
 			  	$(OBJDIR)/moduleBoundary.o $(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.o \
 			  	$(OBJDIR)/moduleCollisions.o $(OBJDIR)/moduleTable.o $(OBJDIR)/moduleParallel.o \
-			  	$(OBJDIR)/moduleEM.o \
+			  	$(OBJDIR)/moduleEM.o $(OBJDIR)/moduleRandom.o \
 			  	$(OBJDIR)/moduleMeshCyl.o	$(OBJDIR)/moduleMeshCylRead.o	$(OBJDIR)/moduleMeshCylBoundary.o \
           $(OBJDIR)/moduleMesh1DCart.o $(OBJDIR)/moduleMesh1DCartRead.o $(OBJDIR)/moduleMesh1DCartBoundary.o \
           $(OBJDIR)/moduleMesh1DRad.o  $(OBJDIR)/moduleMesh1DRadRead.o  $(OBJDIR)/moduleMesh1DRadBoundary.o

src/modules/makefile

@@ -2,20 +2,20 @@
 OBJS = moduleCaseParam.o moduleCompTime.o moduleList.o \
 		   moduleOutput.o moduleInput.o moduleSolver.o \
 			 moduleCollisions.o moduleTable.o moduleParallel.o \
-			 moduleEM.o 
+			 moduleEM.o moduleRandom.o
 
 all: $(OBJS) mesh.o
 
 mesh.o: moduleCollisions.o moduleBoundary.o
 	$(MAKE) -C mesh all
 
-moduleCollisions.o: moduleTable.o moduleSpecies.o moduleRefParam.o moduleConstParam.o moduleCollisions.f90
+moduleCollisions.o: moduleRandom.o moduleTable.o moduleSpecies.o moduleRefParam.o moduleConstParam.o moduleCollisions.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 
 moduleInput.o: moduleParallel.o moduleRefParam.o moduleCaseParam.o moduleSolver.o moduleInject.o moduleBoundary.o moduleErrors.o moduleSpecies.o moduleInput.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 
-moduleInject.o: moduleSpecies.o moduleSolver.o moduleInject.f90
+moduleInject.o: moduleRandom.o moduleSpecies.o moduleSolver.o moduleInject.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 
 moduleList.o: moduleSpecies.o moduleErrors.o moduleList.f90
@@ -24,6 +24,9 @@ moduleList.o: moduleSpecies.o moduleErrors.o moduleList.f90
 moduleOutput.o: moduleSpecies.o moduleRefParam.o moduleOutput.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 
+moduleRandom.o: moduleConstParam.o moduleRandom.f90
+	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
+
 moduleRefParam.o: moduleConstParam.o moduleRefParam.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 

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

@@ -23,7 +23,7 @@ MODULE moduleMesh1DCart
     CONTAINS
       PROCEDURE, PASS:: init     => initEdge1DCart
       PROCEDURE, PASS:: getNodes => getNodes1DCart
-      PROCEDURE, PASS:: randPos  => randPos1DCart
+      PROCEDURE, PASS:: randPos  => randPosEdge
 
   END TYPE meshEdge1DCart
 
@@ -110,6 +110,7 @@ MODULE moduleMesh1DCart
     REAL(8):: arNodes(1:2)
     CONTAINS
       PROCEDURE, PASS::   init        => initVol1DCartSegm
+      PROCEDURE, PASS::   randPos     => randPos1DCartSeg
       PROCEDURE, PASS::   area        => areaSegm
       PROCEDURE, NOPASS:: fPsi        => fPsiSegm
       PROCEDURE, NOPASS:: dPsi        => dPsiSegm
@@ -226,13 +227,13 @@ MODULE moduleMesh1DCart
     END FUNCTION getNodes1DCart
 
     !Calculates a 'random' position in edge
-    FUNCTION randPos1DCart(self) RESULT(r)
+    FUNCTION randPosEdge(self) RESULT(r)
       CLASS(meshEdge1DCart), INTENT(in):: self
       REAL(8):: r(1:3)
 
       r = (/ self%x, 0.D0, 0.D0 /)
 
-    END FUNCTION randPos1DCart
+    END FUNCTION randPosEdge
 
     !VOLUME FUNCTIONS
     !SEGMENT FUNCTIONS
@@ -265,6 +266,24 @@ MODULE moduleMesh1DCart
 
     END SUBROUTINE initVol1DCartSegm
 
+    !Calculates a random position in 1D volume
+    FUNCTION randPos1DCartSeg(self) RESULT(r)
+      USE moduleRandom
+      IMPLICIT NONE
+
+      CLASS(meshVol1DCartSegm), INTENT(in):: self
+      REAL(8):: r(1:3)
+      REAL(8):: xii(1:3)
+      REAL(8), ALLOCATABLE:: fPsi(:)
+
+      xii(1)   = random(-1.D0, 1.D0)
+      xii(2:3) = 0.D0
+
+      fPsi = self%fPsi(xii)
+      r(1) = DOT_PRODUCT(fPsi, self%x)
+
+    END FUNCTION randPos1DCartSeg
+
     !Computes element area
     PURE SUBROUTINE areaSegm(self)
       IMPLICIT NONE
@@ -479,6 +498,7 @@ MODULE moduleMesh1DCart
     END SUBROUTINE nextElementSegm
 
     !COMMON FUNCTIONS FOR 1D VOLUME ELEMENTS
+    !Calculates a random position in 1D volume
     !Computes the element Jacobian determinant
     PURE FUNCTION detJ1DCart(self, xi, dPsi_in) RESULT(dJ)
       IMPLICIT NONE

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

@@ -66,25 +66,16 @@ SUBMODULE (moduleMesh1DCart) moduleMesh1DCartBoundary
     SUBROUTINE wallTemperature(edge, part)
       USE moduleSpecies
       USE moduleBoundary
-      USE moduleConstParam, ONLY: PI
+      USE moduleRandom
       IMPLICIT NONE
 
       CLASS(meshEdge), INTENT(inout):: edge
       CLASS(particle), INTENT(inout):: part
-      REAL(8):: x, y
 
       !Modifies particle velocity according to wall temperature
       SELECT TYPE(bound => edge%boundary%bTypes(part%sp)%obj)
       TYPE IS(boundaryWallTemperature)
-        x = 0.D0
+        part%v(1) = part%v(1) + bound%vTh*randomMaxwellian()
-        DO WHILE (x == 0.D0)
-          CALL RANDOM_NUMBER(x)
-
-        END DO
-
-        CALL RANDOM_NUMBER(y)
-
-        part%v(1) = part%v(1) + bound%vTh*DSQRT(-2.D0*DLOG(x))*DCOS(2.D0*PI*y)
 
       END SELECT
 

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

@@ -75,7 +75,7 @@ MODULE moduleMesh1DRad
       PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi
       PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer
   
-  END TYPE meshVol1Drad
+  END TYPE meshVol1DRad
 
   ABSTRACT INTERFACE
     PURE FUNCTION fPsi_interface(xi) RESULT(fPsi)
@@ -111,6 +111,7 @@ MODULE moduleMesh1DRad
     REAL(8):: arNodes(1:2)
     CONTAINS
       PROCEDURE, PASS::   init        => initVol1DRadSegm
+      PROCEDURE, PASS::   randPos     => randPos1DRadSeg
       PROCEDURE, PASS::   area        => areaRad
       PROCEDURE, NOPASS:: fPsi        => fPsiRad
       PROCEDURE, NOPASS:: dPsi        => dPsiRad
@@ -266,6 +267,24 @@ MODULE moduleMesh1DRad
 
     END SUBROUTINE initVol1DRadSegm
 
+    !Calculates a random position in 1D volume
+    FUNCTION randPos1DRadSeg(self) RESULT(r)
+      USE moduleRandom
+      IMPLICIT NONE
+
+      CLASS(meshVol1DRadSegm), INTENT(in):: self
+      REAL(8):: r(1:3)
+      REAL(8):: xii(1:3)
+      REAL(8), ALLOCATABLE:: fPsi(:)
+
+      xii(1)   = random(-1.D0, 1.D0)
+      xii(2:3) = 0.D0
+
+      fPsi = self%fPsi(xii)
+      r(1) = DOT_PRODUCT(fPsi, self%r)
+
+    END FUNCTION randPos1DRadSeg
+
     !Computes element area
     PURE SUBROUTINE areaRad(self)
       IMPLICIT NONE

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

@@ -68,25 +68,16 @@ SUBMODULE (moduleMesh1DRad) moduleMesh1DRadBoundary
     SUBROUTINE wallTemperature(edge, part)
       USE moduleSpecies
       USE moduleBoundary
-      USE moduleConstParam, ONLY: PI
+      USE moduleRandom
       IMPLICIT NONE
 
       CLASS(meshEdge), INTENT(inout):: edge
       CLASS(particle), INTENT(inout):: part
-      REAL(8):: x, y
 
       !Modifies particle velocity according to wall temperature
       SELECT TYPE(bound => edge%boundary%bTypes(part%sp)%obj)
       TYPE IS(boundaryWallTemperature)
-        x = 0.D0
+        part%v(1) = part%v(1) + bound%vTh*randomMaxwellian()
-        DO WHILE (x == 0.D0)
-          CALL RANDOM_NUMBER(x)
-
-        END DO
-
-        CALL RANDOM_NUMBER(y)
-
-        part%v(1) = part%v(1) + bound%vTh*DSQRT(-2.D0*DLOG(x))*DCOS(2.D0*PI*y)
 
       END SELECT
 

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

@@ -31,7 +31,7 @@ MODULE moduleMeshCyl
     CONTAINS
       PROCEDURE, PASS:: init     => initEdgeCyl
       PROCEDURE, PASS:: getNodes => getNodesCyl
-      PROCEDURE, PASS:: randPos  => randPosCyl
+      PROCEDURE, PASS:: randPos  => randPosEdge
 
   END TYPE meshEdgeCyl
 
@@ -129,6 +129,7 @@ MODULE moduleMeshCyl
 
     CONTAINS
       PROCEDURE, PASS::   init        => initVolQuadCyl
+      PROCEDURE, PASS::   randPos     => randPosVolQuad
       PROCEDURE, PASS::   area        => areaQuad
       PROCEDURE, NOPASS:: fPsi        => fPsiQuad
       PROCEDURE, NOPASS:: dPsi        => dPsiQuad
@@ -161,6 +162,7 @@ MODULE moduleMeshCyl
 
     CONTAINS
       PROCEDURE, PASS::   init        => initVolTriaCyl
+      PROCEDURE, PASS::   randPos     => randPosVolTria
       PROCEDURE, PASS::   area        => areaTria
       PROCEDURE, NOPASS:: fPsi        => fPsiTria
       PROCEDURE, NOPASS:: dPsi        => dPsiTria
@@ -274,6 +276,28 @@ MODULE moduleMeshCyl
 
     END SUBROUTINE initEdgeCyl
 
+    !Random position in quadrilateral volume
+    FUNCTION randPosVolQuad(self) RESULT(r)
+      USE moduleRandom
+      IMPLICIT NONE
+
+      CLASS(meshVolCylQuad), INTENT(in):: self
+      REAL(8):: r(1:3)
+      REAL(8):: xii(1:3)
+      REAL(8), ALLOCATABLE:: fPsi(:)
+
+      xii(1) = random(-1.D0, 1.D0)
+      xii(2) = random(-1.D0, 1.D0)
+      xii(3) = 0.D0
+
+      fPsi = self%fPsi(xii)
+
+      r(1) = DOT_PRODUCT(fPsi, self%z)
+      r(2) = DOT_PRODUCT(fPsi, self%r)
+      r(3) = 0.D0
+
+    END FUNCTION randposVolQuad
+
     !Get nodes from edge
     PURE FUNCTION getNodesCyl(self) RESULT(n)
       IMPLICIT NONE
@@ -287,20 +311,23 @@ MODULE moduleMeshCyl
     END FUNCTION getNodesCyl
 
     !Calculates a random position in edge
-    FUNCTION randPosCyl(self) RESULT(r)
+    FUNCTION randPosEdge(self) RESULT(r)
+      USE moduleRandom
+      IMPLICIT NONE
+
       CLASS(meshEdgeCyl), INTENT(in):: self
       REAL(8):: rnd
       REAL(8):: r(1:3)
       REAL(8):: p1(1:2), p2(1:2)
 
-      CALL RANDOM_NUMBER(rnd)
+      rnd = random()
 
       p1 = (/self%z(1), self%r(1) /)
       p2 = (/self%z(2), self%r(2) /)
       r(1:2) = (1.D0 - rnd)*p1 + rnd*p2
       r(3)   =  0.D0
 
-    END FUNCTION randPosCyl
+    END FUNCTION randPosEdge
 
     !VOLUME FUNCTIONS
     !QUAD FUNCTIONS
@@ -709,6 +736,28 @@ MODULE moduleMeshCyl
 
     END SUBROUTINE initVolTriaCyl
 
+    !Random position in quadrilateral volume
+    FUNCTION randPosVolTria(self) RESULT(r)
+      USE moduleRandom
+      IMPLICIT NONE
+
+      CLASS(meshVolCylTria), INTENT(in):: self
+      REAL(8):: r(1:3)
+      REAL(8):: xii(1:3)
+      REAL(8), ALLOCATABLE:: fPsi(:)
+
+      xii(1) = random( 0.D0, 1.D0)
+      xii(2) = random( 0.D0, 1.D0)
+      xii(3) = 0.D0
+
+      fPsi = self%fPsi(xii)
+
+      r(1) = DOT_PRODUCT(fPsi, self%z)
+      r(2) = DOT_PRODUCT(fPsi, self%r)
+      r(3) = 0.D0
+
+    END FUNCTION randposVolTria
+
     !Calculates area for triangular element
     PURE SUBROUTINE areaTria(self)
       USE moduleConstParam

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

@@ -104,28 +104,19 @@ SUBMODULE (moduleMeshCyl) moduleMeshCylBoundary
     SUBROUTINE wallTemperature(edge, part)
       USE moduleSpecies
       USE moduleBoundary
-      USE moduleConstParam, ONLY: PI
+      USE moduleRandom
       IMPLICIT NONE
 
       CLASS(meshEdge), INTENT(inout):: edge
       CLASS(particle), INTENT(inout):: part
       REAL(8):: edgeNorm, cosT, sinT, rp(1:2), rpp(1:2), vpp(1:2)
       INTEGER:: i
-      REAL(8):: x, y
 
       !Modifies particle velocity according to wall temperature
       SELECT TYPE(bound => edge%boundary%bTypes(part%sp)%obj)
       TYPE IS(boundaryWallTemperature)
         DO i = 1, 3
-          x = 0.D0
+          part%v(i) = part%v(i) + bound%vTh*randomMaxwellian()
-          DO WHILE (x == 0.D0)
-            CALL RANDOM_NUMBER(x)
-
-          END DO
-
-          CALL RANDOM_NUMBER(y)
-
-          part%v(i) = part%v(i) + bound%vTh*DSQRT(-2.D0*DLOG(x))*DCOS(2.D0*PI*y)
 
         END DO
           

src/modules/mesh/moduleMesh.f90

@@ -70,7 +70,7 @@ MODULE moduleMesh
     CONTAINS
       PROCEDURE(initEdge_interface),     DEFERRED, PASS:: init
       PROCEDURE(getNodesEdge_interface), DEFERRED, PASS:: getNodes
-      PROCEDURE(randPos_interface),      DEFERRED, PASS:: randPos
+      PROCEDURE(randPosEdge_interface),      DEFERRED, PASS:: randPos
 
   END TYPE meshEdge
 
@@ -93,12 +93,12 @@ MODULE moduleMesh
 
     END FUNCTION
 
-    FUNCTION randPos_interface(self) RESULT(r)
+    FUNCTION randPosEdge_interface(self) RESULT(r)
       IMPORT:: meshEdge
       CLASS(meshEdge), INTENT(in):: self
       REAL(8):: r(1:3)
 
-    END FUNCTION randPos_interface
+    END FUNCTION randPosEdge_interface
 
   END INTERFACE
 
@@ -138,9 +138,10 @@ MODULE moduleMesh
     REAL(8):: totalWeight = 0.D0
     CONTAINS
       PROCEDURE(initVol_interface),     DEFERRED, PASS::   init
+      PROCEDURE(getNodesVol_interface), DEFERRED, PASS::   getNodes
+      PROCEDURE(randPosVol_interface),  DEFERRED, PASS::   randPos
       PROCEDURE(scatter_interface),     DEFERRED, PASS::   scatter
       PROCEDURE(gatherEF_interface),    DEFERRED, PASS::   gatherEF
-      PROCEDURE(getNodesVol_interface), DEFERRED, PASS::   getNodes
       PROCEDURE(elemF_interface),       DEFERRED, PASS::   elemF
       PROCEDURE,                                  PASS::   findCell
       PROCEDURE(phy2log_interface),     DEFERRED, PASS::   phy2log
@@ -215,6 +216,13 @@ MODULE moduleMesh
 
     END FUNCTION inside_interface
 
+    FUNCTION randPosVol_interface(self) RESULT(r)
+      IMPORT:: meshVol
+      CLASS(meshVol), INTENT(in):: self
+      REAL(8):: r(1:3)
+
+    END FUNCTION randPosVol_interface
+
   END INTERFACE
 
   !Containers for volumes in the mesh
@@ -376,13 +384,13 @@ MODULE moduleMesh
       USE moduleSpecies
       USE moduleList
       use moduleRefParam
+      USE moduleRandom
       IMPLICIT NONE
 
       CLASS(meshVol), INTENT(inout):: self
       INTEGER:: nPart !Number of particles inside the cell
       REAL(8):: pMax !Maximum probability of collision
       INTEGER:: rnd !random index
-      REAL(8):: rndReal
       TYPE(particle), POINTER:: part_i, part_j
       INTEGER:: n !collision
       INTEGER:: ij, k
@@ -403,11 +411,9 @@ MODULE moduleMesh
 
         DO n = 1, self%nColl
           !Select random numbers
-          CALL RANDOM_NUMBER(rndReal)
+          rnd = random(1, nPart)
-          rnd = 1 + FLOOR(nPart*rndReal)
           part_i => partTemp(rnd)%part
-          CALL RANDOM_NUMBER(rndReal)
+          rnd = random(1, nPart)
-          rnd = 1 + FLOOR(nPart*rndReal)
           part_j => partTemp(rnd)%part
           ij = interactionIndex(part_i%sp, part_j%sp)
           sigmaVrelMaxNew = 0.D0

src/modules/moduleCollisions.f90

@@ -148,6 +148,7 @@ MODULE moduleCollisions
                                     part_i, part_j)
       USE moduleSpecies
       USE moduleConstParam
+      USE moduleRandom
       IMPLICIT NONE
 
       CLASS(collisionBinaryElastic), INTENT(in):: self
@@ -160,27 +161,23 @@ MODULE moduleCollisions
       REAL(8):: vp_i, vp_j, v_i, v_j !post and pre-collision velocities
       REAL(8):: v_ij !sum of velocities modules
       REAL(8):: alpha !random angle of scattering
-      REAL(8):: rnd
 
       !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
       sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel
-      CALL RANDOM_NUMBER(rnd)
+      IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN
-      IF (sigmaVrelMaxNew/sigmaVrelMax > rnd) THEN
         !Applies the collision
         v_i = NORM2(part_i%v)
         v_j = NORM2(part_j%v)
         v_ij = v_i+v_j
         vp_j = v_ij*self%w_i
         vp_i = v_ij*self%w_j
-        CALL RANDOM_NUMBER(rnd)
+        alpha = PI*random()
-        alpha = PI*rnd
         part_i%v(1) = vp_i*DCOS(alpha)
         part_i%v(2) = vp_i*DSIN(alpha)
-        CALL RANDOM_NUMBER(rnd)
+        alpha = PI*random()
-        alpha = PI*rnd
         part_j%v(1) = vp_j*DCOS(alpha)
         part_j%v(2) = vp_j*DSIN(alpha)
 
@@ -243,8 +240,7 @@ MODULE moduleCollisions
       USE moduleSpecies
       USE moduleErrors
       USE moduleList
-      USE moduleRefParam !TODO: DELETE
+      USE moduleRandom
-      USE moduleConstParam !TODO: DELETE
       USE OMP_LIB
       IMPLICIT NONE
 
@@ -256,7 +252,6 @@ MODULE moduleCollisions
       TYPE(particle), POINTER:: newElectron
       REAL(8):: vRel, eRel
       REAL(8):: sigmaVrel
-      REAL(8):: rnd
       REAL(8), DIMENSION(1:3):: vp_e, vp_n
 
       !eRel (in units of [m][L]^2[t]^-2
@@ -266,8 +261,7 @@ MODULE moduleCollisions
       IF (eRel > self%eThreshold) THEN
         sigmaVrel = self%crossSec%get(eRel)*vRel
         sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel
-        CALL RANDOM_NUMBER(rnd)
+        IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN
-        IF (sigmaVrelMaxNew/sigmaVrelMax > rnd) THEN
           !Find which particle is the ionizing electron
           IF    (part_i%sp == self%electron%sp) THEN
             electron => part_i
@@ -350,6 +344,7 @@ MODULE moduleCollisions
     SUBROUTINE collideBinaryChargeExchange(self, sigmaVrelMax, sigmaVrelMaxNew, &
                                            part_i, part_j)
       USE moduleSpecies
+      USE moduleRandom
       IMPLICIT NONE
 
       CLASS(collisionBinaryChargeExchange), INTENT(in):: self
@@ -359,15 +354,13 @@ MODULE moduleCollisions
       REAL(8):: sigmaVrel
       REAL(8):: vRel !relative velocity
       REAL(8):: eRel !relative energy
-      REAL(8):: rnd
 
       !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
       sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel
-      CALL RANDOM_NUMBER(rnd)
+      IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN
-      IF (sigmaVrelMaxNew/sigmaVrelMax > rnd) THEN
         SELECT TYPE(sp => species(part_i%sp)%obj)
         TYPE IS (speciesNeutral)
           !Species i is neutral, ionize particle i

src/modules/moduleInject.f90

@@ -174,21 +174,14 @@ MODULE moduleInject
 
     !Random velocity from Maxwellian distribution
     FUNCTION randomVelMaxwellian(self) RESULT (v)
-      USE moduleConstParam, ONLY: PI
+      USE moduleRandom
       IMPLICIT NONE
 
       CLASS(velDistMaxwellian), INTENT(in):: self
       REAL(8):: v
-      REAL(8):: x, y
       v = 0.D0
-      x = 0.D0
 
-      DO WHILE (x == 0.D0)
+      v = self%v + self%vTh*randomMaxwellian()
-        CALL RANDOM_NUMBER(x)
-      END DO
-      CALL RANDOM_NUMBER(y)
-
-      v = self%v + self%vTh*DSQRT(-2.D0*DLOG(x))*DCOS(2.D0*PI*y)
 
     END FUNCTION randomVelMaxwellian
 
@@ -208,6 +201,7 @@ MODULE moduleInject
       USE moduleSpecies
       USE moduleSolver
       USE moduleMesh
+      USE moduleRandom
       IMPLICIT NONE
 
       CLASS(injectGeneric), INTENT(in):: self
@@ -215,7 +209,6 @@ MODULE moduleInject
       INTEGER:: i!, j
       INTEGER, SAVE:: nMin, nMax !Min and Max index in partInj array
       INTEGER:: n
-      REAL(8):: rnd
       CLASS(meshEdge), POINTER:: randomEdge
 
       !Insert particles
@@ -243,8 +236,7 @@ MODULE moduleInject
 
       !$OMP DO
       DO n = nMin, nMax
-        CALL RANDOM_NUMBER(rnd)
+        randomX = random(1, self%nEdges)
-        randomX = INT(DBLE(self%nEdges-1)*rnd) + 1 
 
         randomEdge => mesh%edges(self%edges(randomX))%obj
         !Random position in edge

src/modules/moduleInput.f90

@@ -21,43 +21,68 @@ MODULE moduleInput
       !Loads the config file
       CALL verboseError('Loading input file...')
       CALL config%load(filename = inputFile)
+      CALL checkStatus(config, "load")
 
       !Reads reference parameters
+      CALL verboseError('Reading Reference parameters...')
       CALL readReference(config)
+      CALL checkStatus(config, "readReference")
 
       !Reads output parameters
+      CALL verboseError('Reading Output parameters...')
       CALL readOutput(config)
+      CALL checkStatus(config, "readOutput")
 
       !Read species
       CALL verboseError('Reading species information...')
       CALL readSpecies(config)
+      CALL checkStatus(config, "readSpecies")
 
       !Read interactions between species
       CALL verboseError('Reading interaction between species...')
       CALL readInteractions(config)
+      CALL checkStatus(config, "readInteractions")
 
       !Read boundaries
       CALL verboseError('Reading boundary conditions...')
       CALL readBoundary(config)
+      CALL checkStatus(config, "readBoundary")
 
       !Read Geometry
       CALL verboseError('Reading Geometry...')
       CALL readGeometry(config)
+      CALL checkStatus(config, "readGeometry")
 
       !Reads case parameters
-      CALL verboseError('Reading Case Parameters...')
+      CALL verboseError('Reading Case parameters...')
       CALL readCase(config)
+      CALL checkStatus(config, "readCase")
 
       !Read injection of particles
       CALL verboseError('Reading Interactions between species...')
       CALL readInject(config)
+      CALL checkStatus(config, "readInject")
 
       !Read parallel parameters
       CALL verboseError('Reading Parallel configuration...')
       CALL readParallel(config)
+      CALL checkStatus(config, "readParallel")
 
     END SUBROUTINE readConfig
     
+    !Checks the status of the JSON case file and, if failed, exits the execution.
+    SUBROUTINE checkStatus(config, step)
+      USE moduleErrors
+      USE json_module
+      IMPLICIT NONE
+
+      TYPE(json_file), INTENT(inout):: config
+      CHARACTER(LEN=*), INTENT(in):: step
+
+      IF (config%failed()) CALL criticalError("Error reading the JSON input file", TRIM(step))
+
+    END SUBROUTINE checkStatus
+
     !Reads the reference parameters
     SUBROUTINE readReference(config)
       USE moduleRefParam
@@ -203,8 +228,127 @@ MODULE moduleInput
       WRITE(tString, '(I1)') iterationDigits
       iterationFormat = "(I" // tString // "." // tString // ")"
 
+      !Read initial state for species
+      CALL verboseError('Reading Initial state...')
+      CALL readInitial(config)
+      CALL checkStatus(config, "readInitial")
+
     END SUBROUTINE readCase
 
+    !Reads the initial information for the species
+    SUBROUTINE readInitial(config)
+      USE moduleSpecies
+      USE moduleMesh
+      USE moduleOutput
+      USE moduleRefParam
+      USE moduleRandom
+      USE json_module
+      IMPLICIT NONE
+
+      TYPE(json_file), INTENT(inout):: config
+      LOGICAL:: found
+      CHARACTER(:), ALLOCATABLE:: object
+      INTEGER:: nInitial
+      INTEGER:: i, p, e
+      CHARACTER(LEN=2):: iString
+      CHARACTER(:), ALLOCATABLE:: spName
+      INTEGER:: sp
+      CHARACTER(:), ALLOCATABLE:: spFile
+      INTEGER:: stat
+      CHARACTER(100):: dummy
+      REAL(8):: density, velocity(1:3), temperature
+      INTEGER:: nNewPart = 0.D0
+      TYPE(particle), POINTER:: partNew
+      REAL(8):: vTh
+      TYPE(lNode), POINTER:: partCurr, partNext
+
+      CALL config%info('case.initial', found, n_children = nInitial)
+
+      IF (found) THEN
+        !Reads the information from initial species
+        DO i = 1, nInitial
+          WRITE(iString, '(I2)') i
+          object = 'case.initial(' // iString // ')'
+          CALL config%get(object // '.speciesName', spName, found)
+          sp = speciesName2Index(spName)
+          CALL config%get(object // '.initialState', spFile, found)
+          OPEN (10, FILE = path // spFile, ACTION = 'READ')
+          DO
+            READ(10, '(A)', IOSTAT = stat) dummy
+            !If EoF, exit reading
+            IF (stat /= 0) EXIT
+            !If comment, skip
+            IF (INDEX(dummy,'#') /= 0) CYCLE
+            !Go up one line
+            BACKSPACE(10)
+            !Read information
+            READ(10, *) e, density, velocity, temperature
+            !Scale variables
+            !Particles in cell volume
+            nNewPart = INT(density * (mesh%vols(e)%obj%volume*Vol_ref) / species(sp)%obj%weight)
+            !Non-dimensional velocity
+            velocity = velocity / v_ref
+            !Non-dimensional temperature
+            temperature = temperature / T_ref
+            !Non-dimensional thermal temperature
+            vTh = DSQRT(temperature/species(sp)%obj%m)
+            !Allocate new particles
+            DO p = 1, nNewPart
+              ALLOCATE(partNew)
+              partNew%sp     = sp
+              partNew%v(1)   = velocity(1) + vTh*randomMaxwellian()
+              partNew%v(2)   = velocity(2) + vTh*randomMaxwellian()
+              partNew%v(3)   = velocity(3) + vTh*randomMaxwellian()
+              partNew%vol    = e
+              partNew%r      = mesh%vols(e)%obj%randPos()
+              partNew%xi     = mesh%vols(e)%obj%phy2log(partNew%r)
+              partNew%n_in   = .TRUE.
+              partNew%weight = species(sp)%obj%weight
+              !If charged species, add qm to particle
+              SELECT TYPE(sp => species(sp)%obj)
+              TYPE IS (speciesCharged)
+                partNew%qm = sp%qm
+
+              CLASS DEFAULT
+                partNew%qm = 0.D0
+
+              END SELECT
+
+              !Assign particle to temporal list of particles
+              CALL partInitial%add(partNew)
+
+            END DO
+
+          END DO
+
+        END DO
+
+        !Convert temporal list of particles into initial partOld array
+        !Deallocates the list of initial particles
+        nNewPart = partInitial%amount
+        IF (nNewPart > 0) THEN
+          ALLOCATE(partOld(1:nNewPart))
+          partCurr => partInitial%head
+          DO p = 1, nNewPart
+            partNext => partCurr%next
+            partOld(p) = partCurr%part
+            DEALLOCATE(partCurr)
+            partCurr => partNext
+
+          END DO
+
+          IF (ASSOCIATED(partInitial%head)) NULLIFY(partInitial%head)
+          IF (ASSOCIATED(partInitial%tail)) NULLIFY(partInitial%tail)
+          partInitial%amount = 0
+
+        END IF
+
+        nPartOld = SIZE(partOld)
+
+      END IF
+
+    END SUBROUTINE readInitial
+
     !Reads configuration for the output files
     SUBROUTINE readOutput(config)
       USE moduleErrors

src/modules/moduleList.f90

@@ -24,6 +24,7 @@ MODULE moduleList
   INTEGER(KIND=OMP_LOCK_KIND):: lockWScheme !Lock for the NA list of particles
   TYPE(listNode):: partCollisions !Particles created in collisional process
   INTEGER(KIND=OMP_LOCK_KIND):: lockCollisions !Lock for the NA list of particles
+  TYPE(listNode):: partInitial !Initial distribution of particles
 
   TYPE pointerArray
     TYPE(particle), POINTER:: part

src/modules/moduleRandom.f90

@@ -0,0 +1,69 @@
+MODULE moduleRandom
+  !Interface for random number generator
+  INTERFACE random
+    MODULE PROCEDURE randomReal, randomRealAB, randomIntAB
+
+  END INTERFACE random
+
+  CONTAINS
+    !Returns a Real random number between 0 and 1
+    FUNCTION randomReal() RESULT(rnd)
+      IMPLICIT NONE
+
+      REAL(8):: rnd
+
+      rnd = 0.D0
+      CALL RANDOM_NUMBER(rnd)
+
+    END FUNCTION randomReal
+
+    !Returns a Real random number between a and b
+    FUNCTION randomRealAB(a, b) RESULT(rnd)
+      IMPLICIT NONE
+
+      REAL(8), INTENT(in):: a, b
+      REAL(8):: rnd
+      REAL(8):: rnd01 !random real between 0 and 1
+
+      rnd = 0.D0
+      CALL RANDOM_NUMBER(rnd01)
+
+      rnd = (b - a) * rnd01 + a
+
+    END FUNCTION randomRealAB
+
+    !Returns an Integer random numnber between a and b
+    FUNCTION randomIntAB(a, b) RESULT(rnd)
+      IMPLICIT NONE
+      
+      INTEGER, INTENT(in):: a, b
+      INTEGER:: rnd
+      REAL(8):: rnd01
+
+      rnd = 0.D0
+      CALL RANDOM_NUMBER(rnd01)
+
+      rnd = INT(REAL(b - a) * rnd01) + 1
+
+    END FUNCTION randomIntAB
+
+    !Returns a random number in a Maxwellian distribution of mean 0 and width 1
+    FUNCTION randomMaxwellian() RESULT(rnd)
+      USE moduleConstParam, ONLY: PI
+      IMPLICIT NONE
+
+      REAL(8):: rnd
+      REAL(8):: x, y
+
+      rnd = 0.D0
+      x = 0.D0
+      DO WHILE (x == 0.D0)
+        CALL RANDOM_NUMBER(x)
+      END DO
+      CALL RANDOM_NUMBER(y)
+
+      rnd = DSQRT(-2.D0*DLOG(x))*DCOS(2.D0*PI*y)
+
+    END FUNCTION randomMaxwellian
+
+END MODULE moduleRandom

src/modules/moduleTable.f90

@@ -52,6 +52,7 @@ MODULE moduleTable
       i = 0
       DO 
         READ(id, '(A)', iostat = stat) dummy
+        !TODO: Make this a function
         IF (INDEX(dummy,'#') /= 0) CYCLE
         IF (stat /= 0) EXIT
         !Add data