Almost done with EM, not to modify the EM module itself

src/modules/init/moduleInput.f90

@@ -812,15 +812,15 @@ MODULE moduleInput
 
       END DO
 
-      ! Read linking
+      ! Link physical surfaces to boundaries
-      CALL config%info('boundaries.particles.linking', found, n_children = nPhysicalSurfaces)
+      call config%info('boundaries.particles.linking', found, n_children = nPhysicalSurfaces)
       allocate(boundaryParticlesLinking(1:nPhysicalSurfaces))
       allocate(speciesNames(1:nSpecies))
       do b = 1, nPhysicalSurfaces
         write(iString, '(i2)') b
         object = 'boundary.particles.linking(' // trim(iString) // ')'
-        call config%get(object // '.physicalSurfaces', boundaryParticlesLinking(b)%physicalSurface, found) 
+        call config%get(object // '.physicalSurface', boundaryParticlesLinking(b)%physicalSurface, found) 
-        call config%get(object // '.models',           speciesNames,                                found) 
+        call config%get(object // '.model',           speciesNames,                                found) 
         allocate(boundaryParticlesLinking(b)%speciesIndex(1:nSpecies))
         do s = 1, nSpecies
           boundaryParticlesLinking(b)%speciesIndex(s) = boundaryParticlesName_to_Index(speciesNames(s))
@@ -1093,14 +1093,9 @@ MODULE moduleInput
       TYPE(json_file), INTENT(inout):: config
       CHARACTER(:), ALLOCATABLE:: object
       LOGICAL:: found
-      CHARACTER(:), ALLOCATABLE:: typeEM
+      INTEGER:: b
-      REAL(8):: potential
-      INTEGER:: physicalSurface
-      CHARACTER(:), ALLOCATABLE:: temporalProfile, temporalProfilePath
-      INTEGER:: b, s, n, ni
       CHARACTER(2):: bString
-      INTEGER:: info
+      character(len=100), allocatable:: modelName(:)
-      EXTERNAL:: dgetrf
 
       CALL config%info('boundaries.EM.models', found, n_children = nBoundariesEM)
       
@@ -1109,14 +1104,27 @@ MODULE moduleInput
 
       END IF
       
-      DO b = 1, nBoundaryEM
+      do b = 1, nBoundaryEM
-        WRITE(bString, '(I2)') b
+        write(bString, '(I2)') b
         object = 'boundaries.EM.models(' // TRIM(bString) // ')'
 
         call boundariesEM(b)%init(config, object, b)
 
       end do
 
+      ! Link physical surfaces to boundaries
+      call config%info('boundaries.EM.linking', found, n_children = nPhysicalSurfaces)
+      allocate(boundaryEMLinking(1:nPhysicalSurfaces))
+      allocate(speciesNames(1:nSpecies))
+      do b = 1, nPhysicalSurfaces
+        write(iString, '(i2)') b
+        object = 'boundary.EM.linking(' // trim(iString) // ')'
+        call config%get(object // '.physicalSurface', boundaryEMLinking(b)%physicalSurface, found) 
+        call config%get(object // '.model',           modelName,                            found) 
+        boundaryEMLinking(b)%model => boundariesEM(boundaryEMName_to_Index(modelName(s)))%obj
+
+      end do
+
       ! TODO: Move this to the init of species
       ALLOCATE(qSpecies(1:nSpecies))
       DO s = 1, nSpecies
@@ -1131,37 +1139,6 @@ MODULE moduleInput
 
       END DO
 
-      ! TODO: Do this after the mesh has been read
-      ! Modify K matrix due to boundary conditions
-      DO b = 1, nBoundariesEM
-        SELECT TYPE(boundary => boundariesEM(b)%obj)
-        TYPE IS(boundaryEMDirichlet)
-          DO n = 1, boundary%nNodes
-            ni = boundary%nodes(n)%obj%n
-            mesh%K(ni, :)  = 0.D0
-            mesh%K(ni, ni) = 1.D0
-
-          END DO
-
-        TYPE IS(boundaryEMDirichletTime)
-          DO n = 1, boundary%nNodes
-            ni = boundary%nodes(n)%obj%n
-            mesh%K(ni, :)  = 0.D0
-            mesh%K(ni, ni) = 1.D0
-
-          END DO
-
-        END SELECT
-
-      END DO
-
-      !Compute the PLU factorization of K once boundary conditions have been read
-      CALL dgetrf(mesh%numNodes, mesh%numNodes, mesh%K, mesh%numNodes, mesh%IPIV, info)
-      IF (info /= 0) THEN
-        CALL criticalError('Factorization of K matrix failed', 'readBoundaryEM')
-
-      END IF
-
     END SUBROUTINE readBoundaryEM
 
     !Reads the injection of particles from the boundaries

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

@@ -98,7 +98,7 @@ MODULE moduleMesh1DCart
 
     !EDGE FUNCTIONS
     !Init edge element
-    SUBROUTINE initEdge1DCart(self, n, p, bt)
+    SUBROUTINE initEdge1DCart(self, n, p, btPart, btEM)
       USE moduleSpecies, only: nSpecies
       USE moduleErrors
       IMPLICIT NONE
@@ -106,7 +106,8 @@ MODULE moduleMesh1DCart
       CLASS(meshEdge1DCart), INTENT(out):: self
       INTEGER, INTENT(in):: n
       INTEGER, INTENT(in):: p(:)
-      INTEGER, INTENT(in):: bt(1:nSpecies)
+      INTEGER, INTENT(in):: btPart(1:nSpecies)
+      integer, intent(in):: btEM
       REAL(8), DIMENSION(1:3):: r1
       INTEGER:: s
 
@@ -126,11 +127,14 @@ MODULE moduleMesh1DCart
       allocate(self%boundariesParticle(1:nSpecies))
       !Assign functions to boundary
       do s = 1, nSpecies
-        self%boundariesParticle(s)%obj => boundariesParticle(bt(s))%obj
+        self%boundariesParticle(s)%obj => boundariesParticle(btPart(s))%obj
 
       end do
 
-    END SUBROUTINE initEdge1DCart
+      ! Add nodes to EM boundary
+      call boundariesEMLinking(btEM)%model%addNodes(self%getNodes(self%nNodes))
+
+    end subroutine initEdge1DCart
 
     !Get nodes from edge
     PURE FUNCTION getNodes1DCart(self, nNodes) RESULT(n)

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

@@ -98,7 +98,7 @@ MODULE moduleMesh1DRad
 
     !EDGE FUNCTIONS
     !Init edge element
-    SUBROUTINE initEdge1DRad(self, n, p, bt)
+    SUBROUTINE initEdge1DRad(self, n, p, btPart, btEM)
       USE moduleSpecies, only:nSpecies
       USE moduleErrors
       IMPLICIT NONE
@@ -106,7 +106,8 @@ MODULE moduleMesh1DRad
       CLASS(meshEdge1DRad), INTENT(out):: self
       INTEGER, INTENT(in):: n
       INTEGER, INTENT(in):: p(:)
-      INTEGER, INTENT(in):: bt(1:nSpecies)
+      INTEGER, INTENT(in):: btPart(1:nSpecies)
+      integer, intent(in):: btEM
       REAL(8), DIMENSION(1:3):: r1
       INTEGER:: s
 
@@ -126,11 +127,14 @@ MODULE moduleMesh1DRad
       allocate(self%boundariesParticle(1:nSpecies))
       !Assign functions to boundary
       do s = 1, nSpecies
-        self%boundariesParticle(s)%obj => boundariesParticle(bt(s))%obj
+        self%boundariesParticle(s)%obj => boundariesParticle(btPart(s))%obj
 
       end do
 
-    END SUBROUTINE initEdge1DRad
+      ! Add nodes to EM boundary
+      call boundariesEMLinking(btEM)%model%addNodes(self%getNodes(self%nNodes))
+
+    end subroutine initEdge1DRad
 
     !Get nodes from edge
     PURE FUNCTION getNodes1DRad(self, nNodes) RESULT(n)

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

@@ -133,7 +133,7 @@ MODULE moduleMesh2DCart
 
     !EDGE FUNCTIONS
     !Init edge element
-    SUBROUTINE initEdge2DCart(self, n, p, bt)
+    SUBROUTINE initEdge2DCart(self, n, p, btPart, btEm)
       USE moduleSpecies, only:nSpecies
       USE moduleErrors
       IMPLICIT NONE
@@ -141,7 +141,8 @@ MODULE moduleMesh2DCart
       CLASS(meshEdge2DCart), INTENT(out):: self
       INTEGER, INTENT(in):: n
       INTEGER, INTENT(in):: p(:)
-      INTEGER, INTENT(in):: bt(1:nSpecies)
+      INTEGER, INTENT(in):: btPart(1:nSpecies)
+      integer, intent(in):: btEM
       REAL(8), DIMENSION(1:3):: r1, r2
       INTEGER:: s
 
@@ -165,10 +166,13 @@ MODULE moduleMesh2DCart
       allocate(self%boundariesParticle(1:nSpecies))
       !Assign functions to boundary
       do s = 1, nSpecies
-        self%boundariesParticle(s)%obj => boundariesParticle(bt(s))%obj
+        self%boundariesParticle(s)%obj => boundariesParticle(btPart(s))%obj
 
       end do
 
+      ! Add nodes to EM boundary
+      call boundariesEMLinking(btEM)%model%addNodes(self%getNodes(self%nNodes))
+
     END SUBROUTINE initEdge2DCart
 
     !Get nodes from edge

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

@@ -133,7 +133,7 @@ MODULE moduleMesh2DCyl
 
     !EDGE FUNCTIONS
     !Init edge element
-    SUBROUTINE initEdge2DCyl(self, n, p, bt)
+    SUBROUTINE initEdge2DCyl(self, n, p, btPart, btEM)
       USE moduleSpecies
       USE moduleErrors
       USE moduleConstParam, ONLY: PI
@@ -142,7 +142,8 @@ MODULE moduleMesh2DCyl
       CLASS(meshEdge2DCyl), INTENT(out):: self
       INTEGER, INTENT(in):: n
       INTEGER, INTENT(in):: p(:)
-      INTEGER, INTENT(in):: bt(1:nsPecies)
+      INTEGER, INTENT(in):: btPart(1:nsPecies)
+      integer, intent(in):: btEM
       REAL(8), DIMENSION(1:3):: r1, r2
       INTEGER:: s
 
@@ -174,10 +175,13 @@ MODULE moduleMesh2DCyl
       allocate(self%boundariesParticle(1:nSpecies))
       !Assign functions to boundary
       do s = 1, nSpecies
-        self%boundariesParticle(s)%obj => boundariesParticle(bt(s))%obj
+        self%boundariesParticle(s)%obj => boundariesParticle(btPart(s))%obj
 
       end do
 
+      ! Add nodes to EM boundary
+      call boundariesEMLinking(btEM)%model%addNodes(self%getNodes(self%nNodes))
+
     END SUBROUTINE initEdge2DCyl
 
     !Get nodes from edge

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

@@ -104,7 +104,7 @@ MODULE moduleMesh3DCart
 
     !EDGE FUNCTIONS
     !Init surface element
-    SUBROUTINE initEdge3DCartTria(self, n, p, bt)
+    SUBROUTINE initEdge3DCartTria(self, n, p, btPart, btEM)
       USE moduleSpecies, only: nSpecies
       USE moduleErrors
       USE moduleMath
@@ -114,7 +114,8 @@ MODULE moduleMesh3DCart
       CLASS(meshEdge3DCartTria), INTENT(out):: self
       INTEGER, INTENT(in):: n
       INTEGER, INTENT(in):: p(:)
-      INTEGER, INTENT(in):: bt(1:nSpecies)
+      INTEGER, INTENT(in):: btPart(1:nSpecies)
+      integer, intent(in):: btEM
       REAL(8), DIMENSION(1:3):: r1, r2, r3
       REAL(8), DIMENSION(1:3):: vec1, vec2
       INTEGER:: s
@@ -145,12 +146,15 @@ MODULE moduleMesh3DCart
 
       !Boundary particle linking
       allocate(self%boundariesParticle(1:nSpecies))
-      !Assign functions to boundary
+      !Assign functions to particle boundary
       do s = 1, nSpecies
-        self%boundariesParticle(s)%obj => boundariesParticle(bt(s))%obj
+        self%boundariesParticle(s)%obj => boundariesParticle(btPart(s))%obj
 
       end do
 
+      ! Add nodes to EM boundary
+      call boundariesEMLinking(btEM)%model%addNodes(self%getNodes(self%nNodes))
+
     END SUBROUTINE initEdge3DCartTria
 
     !Get nodes from surface

src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90

@@ -36,7 +36,8 @@ MODULE moduleMeshInputGmsh2
       CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
       REAL(8):: r(1:3) !3 generic coordinates
       INTEGER, ALLOCATABLE:: p(:) !Array for nodes
-      INTEGER:: e = 0, n = 0, eTemp = 0, elemType = 0, bt = 0
+      INTEGER:: e = 0, n = 0, eTemp = 0, elemType = 0
+      integer:: btPart = 0, btEM = 0
       INTEGER:: totalNumElem
       INTEGER:: numEdges
       INTEGER:: boundaryType
@@ -199,9 +200,10 @@ MODULE moduleMeshInputGmsh2
           END SELECT
 
           !Associate boundary condition procedure.
-          bt = physicalSurface_to_id(boundariesParticleLinking,boundaryType)
+          btPart = physicalSurface_to_id(boundariesParticleLinking,boundaryType)
+          btEM   = physicalSurface_to_id(boundariesEMLinking,      boundaryType)
 
-          CALL self%edges(e)%obj%init(n, p, boundariesParticleLinking(bt)%speciesIndex)
+          CALL self%edges(e)%obj%init(n, p, boundariesParticleLinking(btPart)%speciesIndex, btEM)
           DEALLOCATE(p)
 
         END DO

src/modules/mesh/inout/text/moduleMeshInputText.f90

@@ -47,7 +47,7 @@ module moduleMeshInputText
       integer:: physicalID
       integer:: n, c
       integer, allocatable:: p(:)
-      integer:: bt
+      integer:: btPart, btEM
 
       fileID = 10
 
@@ -142,8 +142,9 @@ module moduleMeshInputText
             allocate(p(1))
             p(1) = n
             !Associate boundary condition procedure.
-            bt = physicalSurface_to_id(boundariesParticleLinking, physicalID)
+            btPart = physicalSurface_to_id(boundariesParticleLinking, physicalID)
-            call self%edges(physicalID)%obj%init(physicalID, p, boundariesParticleLinking(bt)%speciesIndex)
+            btEM   = physicalSurface_to_id(boundariesEMLinking,       physicalID)
+            call self%edges(physicalID)%obj%init(physicalID, p, boundariesParticleLinking(btPart)%speciesIndex, btEM)
 
             deallocate(p)
 

src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90

@@ -173,7 +173,7 @@ MODULE moduleMeshInputVTU
       REAL(8), ALLOCATABLE, DIMENSION(:):: coordinates
       INTEGER:: n, e, c
       INTEGER, ALLOCATABLE:: p(:)
-      INTEGER:: bt
+      INTEGER:: btPart, btEM
 
       fileID = 10
 
@@ -365,10 +365,12 @@ MODULE moduleMeshInputVTU
           END SELECT
 
           !Associate boundary condition procedure.
-          bt = physicalSurface_to_id(boundariesParticleLinking,entitiesID(n))
+          btPart = physicalSurface_to_id(boundariesParticleLinking, entitiesID(n))
+          btEM   = physicalSurface_to_id(boundariesEMLinking,       entitiesID(n))
 
           !Init edge
-          CALL self%edges(e)%obj%init(n, p, boundariesParticleLinking(bt)%speciesIndex)
+          CALL self%edges(e)%obj%init(n, p, boundariesParticleLinking(btPart)%speciesIndex, btEM)
+
           DEALLOCATE(p)
 
         END DO

src/modules/mesh/moduleMesh.f90

@@ -495,7 +495,7 @@ MODULE moduleMesh
   END TYPE meshParticles
 
   interface
-    pure module subroutine constructGlobalK(self)
+    module subroutine constructGlobalK(self)
       class(meshParticles), intent(inout):: self
 
     end subroutine constructGlobalK
@@ -626,6 +626,12 @@ MODULE moduleMesh
 
     end subroutine initBoundaryParticle
 
+    module function boundaryParticleName_to_Index(boundaryName) result(bp)
+      character(:), allocatable:: boundaryName
+      integer:: bp
+
+    end function boundaryParticleName_to_Index
+
   end interface
 
   abstract interface
@@ -915,12 +921,13 @@ MODULE moduleMesh
 
   type:: boundaryEMLinking
     integer:: physicalSurface
-    class(boundaryEMGeneric), pointer:: boundary => null()
+    class(boundaryEMGeneric), pointer:: model => null()
 
   end type boundaryEMLinking
 
   INTEGER:: nBoundariesEM
   TYPE(boundaryEMCont), ALLOCATABLE, target:: boundariesEM(:)
+  type(boundaryEMLinking), allocatable, dimension(:):: boundariesEMLinking
 
   !Information of charge and reference parameters for rho vector
   REAL(8), ALLOCATABLE:: qSpecies(:)

src/modules/mesh/moduleMesh@boundaryEM.f90

@@ -1,5 +1,29 @@
 submodule(moduleMesh) boundaryEM
   CONTAINS
+    function boundaryEMName_to_Index(boundaryName) result(bp)
+      use moduleErrors
+      implicit none
+
+      character(:), allocatable:: boundaryName
+      integer:: bp
+      integer:: b
+
+      bp = 0
+      do b = 1, nBoundaries
+        if (boundaryName == boundariesEM(b)%obj%name) then
+          bp = boundariesEM(b)%obj%n
+
+        end if
+
+      end do
+
+      if (bp == 0) then
+        call criticalError('Boundary ' // boundaryName // ' not found', 'boundaryEMName_to_Index')
+      
+      end if
+    
+    end function boundaryEMName_to_Index
+
     subroutine addNodes(self, nodes)
       implicit none
 
@@ -159,7 +183,7 @@ submodule(moduleMesh) boundaryEM
 
       CLASS(boundaryEMDirichlet), INTENT(in):: self
       REAL(8), INTENT(inout):: vectorF(:)
-      INTEGER:: n, ni
+      integer:: n
 
       DO n = 1, self%nNodes
         self%nodes(n)%obj%emData%phi  = self%potential
@@ -178,7 +202,7 @@ submodule(moduleMesh) boundaryEM
       CLASS(boundaryEMDirichletTime), INTENT(in):: self
       REAL(8), INTENT(inout):: vectorF(:)
       REAL(8):: timeFactor
-      INTEGER:: n, ni
+      integer:: n
 
       timeFactor = self%temporalProfile%get(DBLE(timeStep)*tauMin)
 
@@ -190,7 +214,7 @@ submodule(moduleMesh) boundaryEM
 
     END SUBROUTINE applyDirichletTime
 
-    function physicalSurface_to_boundaryEMId(boundaryArray,physicalSurface) result(b)
+    function physicalSurface_to_boundaryEMId(boundaryArray, physicalSurface) result(b)
       implicit none
 
       type(boundaryEMLinking):: boundaryArray(:)
@@ -198,9 +222,10 @@ submodule(moduleMesh) boundaryEM
       integer:: b
       integer:: bt
 
+      b = 0
       do bt=1, nPhysicalSurfaces
         if (boundaryArray(bt)%physicalSurface == physicalSurface) then
-          b = boundaryArray(bt)%boundary%n
+          b = boundaryArray(bt)%model%n
           
           exit
 

src/modules/mesh/moduleMesh@boundaryParticle.f90

@@ -1,7 +1,7 @@
 !moduleMeshBoundary: Boundary functions for the mesh edges
 submodule(moduleMesh) boundaryParticle
   contains
-    module function boundaryParticleName_to_Index(boundaryName) result(bp)
+    function boundaryParticleName_to_Index(boundaryName) result(bp)
       use moduleErrors
       implicit none
 
@@ -19,7 +19,7 @@ submodule(moduleMesh) boundaryParticle
       end do
 
       if (bp == 0) then
-        call criticalError('Boundary ' // boundaryName // ' not found', 'boundaryName2Index')
+        call criticalError('Boundary ' // boundaryName // ' not found', 'boundaryParticleName_to_Index')
       
       end if
     
@@ -503,6 +503,7 @@ submodule(moduleMesh) boundaryParticle
       integer:: b
       integer:: bt
 
+      b = 0
       do bt=1, nPhysicalSurfaces
         if (boundaryArray(bt)%physicalSurface == physicalSurface) then
           b = bt

src/modules/mesh/moduleMesh@elements.f90

@@ -1,7 +1,7 @@
 submodule(moduleMesh) elements
   CONTAINS
     !Reset the output of node
-    PURE module SUBROUTINE resetOutput(self)
+    PURE SUBROUTINE resetOutput(self)
       USE moduleSpecies
       USE moduleOutput
       IMPLICIT NONE
@@ -40,40 +40,73 @@ submodule(moduleMesh) elements
     end function meshNodePointer_equal_type_int
 
     !Constructs the global K matrix
-    PURE module SUBROUTINE constructGlobalK(self)
+    SUBROUTINE constructGlobalK(self)
       IMPLICIT NONE
 
       CLASS(meshParticles), INTENT(inout):: self
       INTEGER:: c
-      INTEGER, ALLOCATABLE:: n(:)
+      INTEGER, ALLOCATABLE:: nodes(:)
       REAL(8), ALLOCATABLE:: localK(:,:)
       INTEGER:: i, j
+      integer:: n, b, ni
+      INTEGER:: info
+      EXTERNAL:: dgetrf
 
       DO c = 1, self%numCells
         associate(nNodes => self%cells(c)%obj%nNodes)
-          ALLOCATE(n(1:nNodes))
+          ALLOCATE(nodes(1:nNodes))
           ALLOCATE(localK(1:nNodes, 1:nNodes))
-          n      = self%cells(c)%obj%getNodes(nNodes)
+          nodes  = self%cells(c)%obj%getNodes(nNodes)
           localK = self%cells(c)%obj%elemK(nNodes)
 
           DO i = 1, nNodes
             DO j = 1, nNodes
-              self%K(n(i), n(j)) = self%K(n(i), n(j)) + localK(i, j)
+              self%K(nodes(i), nodes(j)) = self%K(nodes(i), nodes(j)) + localK(i, j)
 
             END DO
 
           END DO
 
-          DEALLOCATE(n, localK)
+          DEALLOCATE(nodes, localK)
 
         end associate
 
       END DO
 
+      ! Modify K matrix due to EM boundary conditions
+      DO b = 1, nBoundariesEM
+        SELECT TYPE(boundary => boundariesEM(b)%obj)
+        TYPE IS(boundaryEMDirichlet)
+          DO n = 1, boundary%nNodes
+            ni = boundary%nodes(n)%obj%n
+            self%K(ni, :)  = 0.D0
+            self%K(ni, ni) = 1.D0
+
+          END DO
+
+        TYPE IS(boundaryEMDirichletTime)
+          DO n = 1, boundary%nNodes
+            ni = boundary%nodes(n)%obj%n
+            self%K(ni, :)  = 0.D0
+            self%K(ni, ni) = 1.D0
+
+          END DO
+
+        END SELECT
+
+      END DO
+
+      !Compute the PLU factorization of K once boundary conditions have been read
+      CALL dgetrf(self%numNodes, self%numNodes, self%K, self%numNodes, self%IPIV, info)
+      IF (info /= 0) THEN
+        CALL criticalError('Factorization of K matrix failed', 'readBoundaryEM')
+
+      END IF
+
     END SUBROUTINE constructGlobalK
 
     ! Gather the value of valNodes at position Xi of an edge
-    pure module function gatherF_edge_scalar(self, Xi, nNodes, valNodes) RESULT(f)
+    pure function gatherF_edge_scalar(self, Xi, nNodes, valNodes) RESULT(f)
       implicit none
 
       class(meshEdge), intent(in):: self
@@ -89,7 +122,7 @@ submodule(moduleMesh) elements
     end function gatherF_edge_scalar
 
     !Gather the value of valNodes (scalar) at position Xi
-    PURE module FUNCTION gatherF_cell_scalar(self, Xi, nNodes, valNodes) RESULT(f)
+    PURE FUNCTION gatherF_cell_scalar(self, Xi, nNodes, valNodes) RESULT(f)
       IMPLICIT NONE
 
       CLASS(meshCell), INTENT(in):: self
@@ -105,7 +138,7 @@ submodule(moduleMesh) elements
     END FUNCTION gatherF_cell_scalar
 
     !Gather the value of valNodes (array) at position Xi
-    PURE module FUNCTION gatherF_cell_array(self, Xi, nNodes, valNodes) RESULT(f)
+    PURE FUNCTION gatherF_cell_array(self, Xi, nNodes, valNodes) RESULT(f)
       IMPLICIT NONE
 
       CLASS(meshCell), INTENT(in):: self
@@ -121,7 +154,7 @@ submodule(moduleMesh) elements
     END FUNCTION gatherF_cell_array
 
     !Gather the spatial derivative of valNodes (scalar) at position Xi
-    PURE module FUNCTION gatherDF_cell_scalar(self, Xi, nNodes, valNodes) RESULT(df)
+    PURE FUNCTION gatherDF_cell_scalar(self, Xi, nNodes, valNodes) RESULT(df)
       IMPLICIT NONE
 
       CLASS(meshCell), INTENT(in):: self
@@ -146,7 +179,7 @@ submodule(moduleMesh) elements
     END FUNCTION gatherDF_cell_scalar
 
     !Scatters particle properties into cell nodes
-    module SUBROUTINE scatter(self, nNodes, part)
+    SUBROUTINE scatter(self, nNodes, part)
       USE moduleMath
       USE moduleSpecies
       USE OMP_LIB
@@ -184,7 +217,7 @@ submodule(moduleMesh) elements
     END SUBROUTINE scatter
 
     !Find next cell for particle
-    RECURSIVE module SUBROUTINE findCell(self, part, oldCell)
+    RECURSIVE SUBROUTINE findCell(self, part, oldCell)
       USE moduleSpecies
       USE moduleErrors
       USE OMP_LIB
@@ -248,7 +281,7 @@ submodule(moduleMesh) elements
     END SUBROUTINE findCell
 
     !If Coll and Particle are the same, simply copy the part%cell into part%cellColl
-    module SUBROUTINE findCellSameMesh(part)
+    SUBROUTINE findCellSameMesh(part)
       USE moduleSpecies
       IMPLICIT NONE
 
@@ -261,7 +294,7 @@ submodule(moduleMesh) elements
     !TODO: try to combine this with the findCell for a regular mesh
     !Find the volume in which particle reside in the mesh for collisions
     !No boundary interaction taken into account
-    module SUBROUTINE findCellCollMesh(part)
+    SUBROUTINE findCellCollMesh(part)
       USE moduleSpecies
       IMPLICIT NONE
 
@@ -311,7 +344,7 @@ submodule(moduleMesh) elements
     !Returns index of volume associated to a position (if any)
     !If no voulme is found, returns 0
     !WARNING: This function is slow and should only be used in initialization phase
-    module FUNCTION findCellBrute(self, r) RESULT(nVol)
+    FUNCTION findCellBrute(self, r) RESULT(nVol)
       USE moduleSpecies
       IMPLICIT NONE
 
@@ -337,7 +370,7 @@ submodule(moduleMesh) elements
     END FUNCTION findCellBrute
 
     !Computes collisions in element
-    module SUBROUTINE doCollisions(self)
+    SUBROUTINE doCollisions(self)
       USE moduleCollisions
       USE moduleSpecies
       USE moduleList
@@ -501,7 +534,7 @@ submodule(moduleMesh) elements
 
     END SUBROUTINE doCollisions
 
-    module SUBROUTINE doCoulomb(self)
+    SUBROUTINE doCoulomb(self)
       USE moduleCoulomb
       USE moduleRandom
       USE moduleOutput