diff --git a/src/fpakc.f90 b/src/fpakc.f90
index e90a5e0..ae6e7bb 100644
--- a/src/fpakc.f90
+++ b/src/fpakc.f90
@@ -11,12 +11,12 @@ PROGRAM fpakc
   USE OMP_LIB
   IMPLICIT NONE
 
-  ! t = time step
-  INTEGER:: t
   ! arg1 = Input argument 1 (input file)
   CHARACTER(200):: arg1
   ! inputFile = path+name of input file
   CHARACTER(:), ALLOCATABLE:: inputFile
+  ! generic integer for time step
+  INTEGER:: t
 
   tStep = omp_get_wtime()
   !Gets the input file
@@ -32,10 +32,13 @@ PROGRAM fpakc
   CALL initOutput(inputFile)
 
   !Do '0' iteration
-  t = tInitial
+  timeStep = tInitial
 
   !$OMP PARALLEL DEFAULT(SHARED)
   !$OMP SINGLE
+  ! Initial reset of probes
+  CALL resetProbes()
+
   CALL verboseError("Initial scatter of particles...")
   !$OMP END SINGLE
   CALL doScatter()
@@ -49,19 +52,21 @@ PROGRAM fpakc
   tStep = omp_get_wtime() - tStep
 
   !Output initial state
-  CALL doOutput(t)
+  CALL doOutput()
   CALL verboseError('Starting main loop...')
   !$OMP PARALLEL DEFAULT(SHARED) 
   DO t = tInitial + 1, tFinal
-    !Insert new particles and push them
     !$OMP SINGLE
     tStep = omp_get_wtime()
 
+    ! Update global time step index
+    timeStep = t
+
     !Checks if a species needs to me moved in this iteration
-    CALL solver%updatePushSpecies(t)
+    CALL solver%updatePushSpecies()
 
     !Checks if probes need to be calculated this iteration
-    CALL resetProbes(t)
+    CALL resetProbes()
     tPush = omp_get_wtime()
     !$OMP END SINGLE
 
@@ -79,7 +84,7 @@ PROGRAM fpakc
     !$OMP END SINGLE
 
     IF (doMCCollisions) THEN
-      CALL meshForMCC%doCollisions(t)
+      CALL meshForMCC%doCollisions()
 
     END IF
 
@@ -124,12 +129,12 @@ PROGRAM fpakc
     !$OMP SINGLE
     tEMField = omp_get_wtime() - tEMField
 
-    CALL doAverage(t)
+    CALL doAverage()
 
     tStep = omp_get_wtime() - tStep
 
     !Output data
-    CALL doOutput(t)
+    CALL doOutput()
     !$OMP END SINGLE
 
   END DO
diff --git a/src/modules/common/moduleCaseParam.f90 b/src/modules/common/moduleCaseParam.f90
index 8df3210..551d867 100644
--- a/src/modules/common/moduleCaseParam.f90
+++ b/src/modules/common/moduleCaseParam.f90
@@ -2,8 +2,13 @@
 MODULE moduleCaseParam
   !Final and initial iterations
   INTEGER:: tFinal, tInitial = 0
+  ! Global index of current iteration
+  INTEGER:: timeStep
+  ! Time step for all species
   REAL(8), ALLOCATABLE:: tau(:)
+  ! Minimum time step
   REAL(8):: tauMin
+  ! Time step for Monte-Carlo Collisions
   REAL(8):: tauColl
 
 END MODULE moduleCaseParam
diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90
index 4216c73..6800ac1 100644
--- a/src/modules/init/moduleInput.f90
+++ b/src/modules/init/moduleInput.f90
@@ -264,8 +264,8 @@ MODULE moduleInput
           CALL readEMBoundary(config)
           !Read constant magnetic field
           DO i = 1, 3
-            WRITE(istring, '(i2)') i
-            CALL config%get(object // '.B(' // istring // ')', B(i), found)
+            WRITE(iString, '(i2)') i
+            CALL config%get(object // '.B(' // iString // ')', B(i), found)
             IF (.NOT. found) THEN
               CALL criticalError('Constant magnetic field not provided in direction ' // iString, 'readSolver')
 
@@ -799,7 +799,7 @@ MODULE moduleInput
 
       TYPE(json_file), INTENT(inout):: config
       INTEGER:: i, s
-      CHARACTER(2):: istring, sString
+      CHARACTER(2):: iString, sString
       CHARACTER(:), ALLOCATABLE:: object, bType
       REAL(8):: Tw, cw !Wall temperature and specific heat
       !Neutral Properties
@@ -815,8 +815,8 @@ MODULE moduleInput
       CALL config%info('boundary', found, n_children = nBoundary)
       ALLOCATE(boundary(1:nBoundary))
       DO i = 1, nBoundary
-        WRITE(istring, '(i2)') i
-        object = 'boundary(' // TRIM(istring) // ')'
+        WRITE(iString, '(i2)') i
+        object = 'boundary(' // TRIM(iString) // ')'
 
         boundary(i)%n = i
         CALL config%get(object // '.name',            boundary(i)%name,            found)
@@ -1100,13 +1100,13 @@ MODULE moduleInput
       TYPE(json_file), INTENT(inout):: config
       CHARACTER(:), ALLOCATABLE:: object
       LOGICAL:: found
-      CHARACTER(2):: istring
+      CHARACTER(2):: iString
       INTEGER:: i
       CHARACTER(:), ALLOCATABLE:: speciesName
       REAL(8), ALLOCATABLE, DIMENSION(:):: r
       REAL(8), ALLOCATABLE, DIMENSION(:):: v1, v2, v3
       INTEGER, ALLOCATABLE, DIMENSION(:):: points
-      REAL(8):: timeStep
+      REAL(8):: everyTimeStep
 
       CALL config%info('output.probes', found, n_children = nProbes)
 
@@ -1114,7 +1114,7 @@ MODULE moduleInput
 
       DO i = 1, nProbes
         WRITE(iString, '(I2)') i
-        object = 'output.probes(' // trim(istring) // ')'
+        object = 'output.probes(' // trim(iString) // ')'
 
         CALL config%get(object // '.species', speciesName, found)
         CALL config%get(object // '.position', r, found)
@@ -1122,16 +1122,14 @@ MODULE moduleInput
         CALL config%get(object // '.velocity_2', v2, found)
         CALL config%get(object // '.velocity_3', v3, found)
         CALL config%get(object // '.points', points, found)
-        CALL config%get(object // '.timeStep', timeStep, found)
+        CALL config%get(object // '.timeStep', everyTimeStep, found)
 
         IF (ANY(points < 2)) CALL criticalError("Number of points in probe " // iString // " incorrect", 'readProbes')
 
-        CALL probe(i)%init(i, speciesName, r, v1, v2, v3, points, timeStep)
+        CALL probe(i)%init(i, speciesName, r, v1, v2, v3, points, everyTimeStep)
 
       END DO
 
-      CALL resetProbes(tInitial)
-
     END SUBROUTINE readProbes
 
     SUBROUTINE readEMBoundary(config)
@@ -1147,40 +1145,56 @@ MODULE moduleInput
       TYPE(json_file), INTENT(inout):: config
       CHARACTER(:), ALLOCATABLE:: object
       LOGICAL:: found
-      CHARACTER(2):: istring
-      INTEGER:: i, e, s
+      CHARACTER(:), ALLOCATABLE:: typeEM
+      REAL(8):: potential
+      INTEGER:: physicalSurface
+      CHARACTER(:), ALLOCATABLE:: timeProfile
+      INTEGER:: b, e, s, n, ni
+      CHARACTER(2):: bString
       INTEGER:: info
       EXTERNAL:: dgetrf
 
       CALL config%info('boundaryEM', found, n_children = nBoundaryEM)
       
-      IF (found) ALLOCATE(boundEM(1:nBoundaryEM))
+      IF (found) ALLOCATE(boundaryEM(1:nBoundaryEM))
       
-      DO i = 1, nBoundaryEM
-        WRITE(istring, '(I2)') i
-        object = 'boundaryEM(' // trim(istring) // ')'
+      DO b = 1, nBoundaryEM
+        WRITE(bString, '(I2)') b
+        object = 'boundaryEM(' // trim(bString) // ')'
 
-        CALL config%get(object // '.type', boundEM(i)%typeEM, found)
+        CALL config%get(object // '.type', typeEM, found)
 
-        SELECT CASE(boundEM(i)%typeEM)
+        SELECT CASE(typeEM)
         CASE ("dirichlet")
-          CALL config%get(object // '.potential', boundEM(i)%potential, found)
-          IF (.NOT. found) &
+          CALL config%get(object // '.potential', potential, found)
+          IF (.NOT. found) THEN
             CALL criticalError('Required parameter "potential" for Dirichlet boundary condition not found', 'readEMBoundary')
-          boundEM(i)%potential = boundEM(i)%potential/Volt_ref
 
-          CALL config%get(object // '.physicalSurface', boundEM(i)%physicalSurface, found)
-          IF (.NOT. found) &
+          END IF
+
+          CALL config%get(object // '.physicalSurface', physicalSurface, found)
+          IF (.NOT. found) THEN
             CALL criticalError('Required parameter "physicalSurface" for Dirichlet boundary condition not found', 'readEMBoundary')
 
+          END IF
+
+          CALL initDirichlet(boundaryEM(b)%obj, physicalSurface, potential)
+
         CASE ("dirichletTime")
-          CALL config%get(object // '.potential', boundEM(i)%potential, found)
-          IF (.NOT. found) &
+          CALL config%get(object // '.potential', potential, found)
+          IF (.NOT. found) THEN
             CALL criticalError('Required parameter "potential" for Dirichlet boundary condition not found', 'readEMBoundary')
-          boundEM(i)%potential = boundEM(i)%potential/Volt_ref
+
+          END IF
+
+          CALL config%get(object // '.physicalSurface', physicalSurface, found)
+          IF (.NOT. found) THEN
+            CALL criticalError('Required parameter "physicalSurface" for Dirichlet boundary condition not found', 'readEMBoundary')
+
+          END IF
 
         CASE DEFAULT
-          CALL criticalError('Boundary type ' // boundEM(i)%typeEM // ' not yet supported', 'readEMBoundary')
+          CALL criticalError('Boundary type ' // typeEM // ' not yet supported', 'readEMBoundary')
 
         END SELECT
 
@@ -1199,18 +1213,28 @@ MODULE moduleInput
 
       END DO
 
-      IF (ALLOCATED(boundEM)) THEN
-        DO e = 1, mesh%numEdges
-          IF (ANY(mesh%edges(e)%obj%physicalSurface == boundEM(:)%physicalSurface)) THEN
-            DO i = 1, nBoundaryEM
-              IF (mesh%edges(e)%obj%physicalSurface == boundEM(i)%physicalSurface) THEN
-                CALL boundEM(i)%apply(mesh%edges(e)%obj)
+      ! Modify K matrix due to boundary conditions
+      DO b = 1, nBoundaryEM
+        SELECT TYPE(boundary => boundaryEM(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 IF
-            END DO
-          END IF
-        END DO
-      END IF
+          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)
@@ -1231,13 +1255,13 @@ MODULE moduleInput
 
       TYPE(json_file), INTENT(inout):: config
       INTEGER:: i
-      CHARACTER(2):: istring
+      CHARACTER(2):: iString
       CHARACTER(:), ALLOCATABLE:: object
       LOGICAL:: found
       CHARACTER(:), ALLOCATABLE:: speciesName
       CHARACTER(:), ALLOCATABLE:: name
       REAL(8):: v
-      REAL(8), ALLOCATABLE:: T(:), normal(:)
+      REAL(8), ALLOCATABLE:: temperature(:), normal(:)
       REAL(8):: flow
       CHARACTER(:), ALLOCATABLE:: units
       INTEGER:: physicalSurface
@@ -1248,8 +1272,8 @@ MODULE moduleInput
       ALLOCATE(inject(1:nInject))
       nPartInj = 0
       DO i = 1, nInject
-        WRITE(istring, '(i2)') i
-        object = 'inject(' // trim(istring) // ')'
+        WRITE(iString, '(i2)') i
+        object = 'inject(' // trim(iString) // ')'
 
         !Find species
         CALL config%get(object // '.species', speciesName, found)
@@ -1257,7 +1281,7 @@ MODULE moduleInput
 
         CALL config%get(object // '.name',            name,            found)
         CALL config%get(object // '.v',               v,               found)
-        CALL config%get(object // '.T',               T,               found)
+        CALL config%get(object // '.T',               temperature,     found)
         CALL config%get(object // '.n',               normal,          found)
         IF (.NOT. found) THEN
           ALLOCATE(normal(1:3))
@@ -1269,7 +1293,7 @@ MODULE moduleInput
         particlesPerEdge = 0
         CALL config%get(object // '.particlesPerEdge', particlesPerEdge, found)
 
-        CALL inject(i)%init(i, v, normal, T, flow, units, sp, physicalSurface, particlesPerEdge)
+        CALL inject(i)%init(i, v, normal, temperature, flow, units, sp, physicalSurface, particlesPerEdge)
 
         CALL readVelDistr(config, inject(i), object)
         
@@ -1329,28 +1353,28 @@ MODULE moduleInput
       TYPE(injectGeneric), INTENT(inout):: inj
       CHARACTER(:), ALLOCATABLE, INTENT(in):: object
       INTEGER:: i
-      CHARACTER(2):: istring
+      CHARACTER(2):: iString
       CHARACTER(:), ALLOCATABLE:: fvType
       LOGICAL:: found
-      REAL(8):: v, T, m
+      REAL(8):: v, temperature, m
 
       !Reads species mass
       m = inj%species%m
       !Reads distribution functions for velocity
       DO i = 1, 3
-        WRITE(istring, '(i2)') i
-        CALL config%get(object // '.velDist('// TRIM(istring) //')', fvType, found)
-        IF (.NOT. found) CALL criticalError("No velocity distribution in direction " // istring // &
+        WRITE(iString, '(i2)') i
+        CALL config%get(object // '.velDist('// TRIM(iString) //')', fvType, found)
+        IF (.NOT. found) CALL criticalError("No velocity distribution in direction " // iString // &
                                             " found for " // object, 'readVelDistr')
 
         SELECT CASE(fvType)
         CASE ("Maxwellian")
-          T = inj%T(i)
-          CALL initVelDistMaxwellian(inj%v(i)%obj, t, m)
+          temperature = inj%temperature(i)
+          CALL initVelDistMaxwellian(inj%v(i)%obj, temperature, m)
 
         CASE ("Half-Maxwellian")
-          T = inj%T(i)
-          CALL initVelDistHalfMaxwellian(inj%v(i)%obj, t, m)
+          temperature = inj%temperature(i)
+          CALL initVelDistHalfMaxwellian(inj%v(i)%obj, temperature, m)
 
         CASE ("Delta")
           v = inj%vMod*inj%n(i)
diff --git a/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 b/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90
index 97ec729..c0dcfbb 100644
--- a/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90
+++ b/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90
@@ -1,22 +1,22 @@
 MODULE moduleMeshOutput0D
 
   CONTAINS
-    SUBROUTINE printOutput0D(self, t)
+    SUBROUTINE printOutput0D(self)
       USE moduleMesh
       USE moduleRefParam
       USE moduleSpecies
       USE moduleOutput
+      USE moduleCaseParam, ONLY: timeStep
       IMPLICIT NONE
 
       CLASS(meshParticles), INTENT(in):: self
-      INTEGER, INTENT(in):: t
       INTEGER:: i
       TYPE(outputFormat):: output
       CHARACTER(:), ALLOCATABLE:: fileName
 
       DO i = 1, nSpecies
         fileName='OUTPUT_' // species(i)%obj%name // '.dat'
-        IF (t == 0) THEN
+        IF (timeStep == 0) THEN
           OPEN(20, file = path // folder // '/' // fileName, action = 'write')
           WRITE(20, "(A1, 14X, A5, A20, 40X, A20, 2(A20))") "#","t (s)","density (m^-3)", "velocity (m/s)", &
                                                                         "pressure (Pa)", "temperature (K)"
@@ -27,14 +27,17 @@ MODULE moduleMeshOutput0D
         
         OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write')
         CALL calculateOutput(self%nodes(1)%obj%output(i), output, self%nodes(1)%obj%v, species(i)%obj)
-        WRITE(20, "(7(ES20.6E3))") REAL(t)*tauMin*ti_ref, output%density, output%velocity, output%pressure, output%temperature
+        WRITE(20, "(7(ES20.6E3))") REAL(timeStep)*tauMin*ti_ref, output%density,    &
+                                                                 output%velocity,   &
+                                                                 output%pressure,   &
+                                                                 output%temperature
         CLOSE(20)
 
       END DO
 
     END SUBROUTINE printOutput0D
 
-    SUBROUTINE printColl0D(self, t)
+    SUBROUTINE printColl0D(self)
       USE moduleMesh
       USE moduleRefParam
       USE moduleCaseParam
@@ -43,12 +46,11 @@ MODULE moduleMeshOutput0D
       IMPLICIT NONE
 
       CLASS(meshGeneric), INTENT(in):: self
-      INTEGER, INTENT(in):: t
       CHARACTER(:), ALLOCATABLE:: fileName
       INTEGER:: k
 
       fileName='OUTPUT_Collisions.dat'
-      IF (t == tInitial) THEN
+      IF (timeStep == tInitial) THEN
         OPEN(20, file = path // folder // '/' // fileName, action = 'write')
         WRITE(20, "(A1, 14X, A5, A20)") "#","t (s)","collisions"
         WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
@@ -57,12 +59,12 @@ MODULE moduleMeshOutput0D
       END IF
       
       OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write')
-      WRITE(20, "(ES20.6E3, 10I20)") REAL(t)*tauMin*ti_ref, (self%cells(1)%obj%tallyColl(k)%tally, k=1,nCollPairs)
+      WRITE(20, "(ES20.6E3, 10I20)") REAL(timeStep)*tauMin*ti_ref, (self%cells(1)%obj%tallyColl(k)%tally, k=1,nCollPairs)
       CLOSE(20)
 
     END SUBROUTINE printColl0D
 
-    SUBROUTINE printEM0D(self, t)
+    SUBROUTINE printEM0D(self)
       USE moduleMesh
       USE moduleRefParam
       USE moduleCaseParam
@@ -70,7 +72,6 @@ MODULE moduleMeshOutput0D
       IMPLICIT NONE
 
       CLASS(meshParticles), INTENT(in):: self
-      INTEGER, INTENT(in):: t
 
     END SUBROUTINE printEM0D
 
diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90
index ccdcf3d..8176bb5 100644
--- a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90
+++ b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90
@@ -80,50 +80,50 @@ MODULE moduleMeshOutputGmsh2
     END SUBROUTINE writeGmsh2FooterElementData
 
     !Prints the scattered properties of particles into the nodes
-    SUBROUTINE printOutputGmsh2(self, t)
+    SUBROUTINE printOutputGmsh2(self)
       USE moduleMesh
       USE moduleRefParam
       USE moduleSpecies
       USE moduleOutput
       USE moduleMeshInoutCommon
+      USE moduleCaseParam, ONLY: timeStep
       IMPLICIT NONE
 
       CLASS(meshParticles), INTENT(in):: self
-      INTEGER, INTENT(in):: t
       INTEGER:: n, i
       TYPE(outputFormat):: output(1:self%numNodes)
       REAL(8):: time
       CHARACTER(:), ALLOCATABLE:: fileName
 
-      time = DBLE(t)*tauMin*ti_ref
+      time = DBLE(timeStep)*tauMin*ti_ref
 
       DO i = 1, nSpecies
-        fileName = formatFileName(prefix, species(i)%obj%name, 'msh', t)
+        fileName = formatFileName(prefix, species(i)%obj%name, 'msh', timeStep)
         WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
         OPEN (60, file = path // folder // '/' //  fileName)
 
         CALL writeGmsh2HeaderMesh(60)
 
-        CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' density (m^-3)', t, time, 1, self%numNodes)
+        CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' density (m^-3)', timeStep, time, 1, self%numNodes)
         DO n=1, self%numNodes
           CALL calculateOutput(self%nodes(n)%obj%output(i), output(n), self%nodes(n)%obj%v, species(i)%obj)
           WRITE(60, "(I6,ES20.6E3)") n, output(n)%density
         END DO
         CALL writeGmsh2FooterNodeData(60)
 
-        CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' velocity (m s^-1)', t, time, 3, self%numNodes)
+        CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' velocity (m s^-1)', timeStep, time, 3, self%numNodes)
         DO n=1, self%numNodes
           WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%velocity
         END DO
         CALL writeGmsh2FooterNodeData(60)
 
-        CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Pressure (Pa)', t, time, 1, self%numNodes)
+        CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Pressure (Pa)', timeStep, time, 1, self%numNodes)
         DO n=1, self%numNodes
           WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%pressure
         END DO
         CALL writeGmsh2FooterNodeData(60)
 
-        CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Temperature (K)', t, time, 1, self%numNodes)
+        CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Temperature (K)', timeStep, time, 1, self%numNodes)
         DO n=1, self%numNodes
           WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%temperature
         END DO
@@ -135,7 +135,7 @@ MODULE moduleMeshOutputGmsh2
     END SUBROUTINE printOutputGmsh2
 
     !Prints the number of collisions into the volumes
-    SUBROUTINE printCollGmsh2(self, t)
+    SUBROUTINE printCollGmsh2(self)
       USE moduleMesh
       USE moduleRefParam
       USE moduleCaseParam
@@ -145,7 +145,6 @@ MODULE moduleMeshOutputGmsh2
       IMPLICIT NONE
 
       CLASS(meshGeneric), INTENT(in):: self
-      INTEGER, INTENT(in):: t
       INTEGER:: numEdges
       INTEGER:: k, c
       INTEGER:: n
@@ -167,9 +166,9 @@ MODULE moduleMeshOutputGmsh2
       END SELECT
 
       IF (collOutput) THEN
-        time = DBLE(t)*tauMin*ti_ref
+        time = DBLE(timeStep)*tauMin*ti_ref
 
-        fileName = formatFileName(prefix, 'Collisions', 'msh', t)
+        fileName = formatFileName(prefix, 'Collisions', 'msh', timeStep)
         WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
         OPEN (60, file = path // folder // '/' //  fileName)
 
@@ -179,7 +178,7 @@ MODULE moduleMeshOutputGmsh2
           DO c = 1, interactionMatrix(k)%amount
             WRITE(cString, "(I2)") c
             title = '"Pair ' // interactionMatrix(k)%sp_i%name // '-' // interactionMatrix(k)%sp_j%name // ' collision ' // cString
-            CALL writeGmsh2HeaderElementData(60, title, t, time, 1, self%numCells)
+            CALL writeGmsh2HeaderElementData(60, title, timeStep, time, 1, self%numCells)
             DO n=1, self%numCells
               WRITE(60, "(I6,I10)") n + numEdges, self%cells(n)%obj%tallyColl(k)%tally(c)
             END DO
@@ -196,7 +195,7 @@ MODULE moduleMeshOutputGmsh2
     END SUBROUTINE printCollGmsh2
 
     !Prints the electrostatic EM properties into the nodes and volumes
-    SUBROUTINE printEMGmsh2(self, t)
+    SUBROUTINE printEMGmsh2(self)
       USE moduleMesh
       USE moduleRefParam
       USE moduleCaseParam
@@ -205,7 +204,6 @@ MODULE moduleMeshOutputGmsh2
       IMPLICIT NONE
 
       CLASS(meshParticles), INTENT(in):: self
-      INTEGER, INTENT(in):: t
       INTEGER:: n, e
       REAL(8):: time
       CHARACTER(:), ALLOCATABLE:: fileName
@@ -214,27 +212,27 @@ MODULE moduleMeshOutputGmsh2
       Xi = (/ 0.D0, 0.D0, 0.D0 /)
 
       IF (emOutput) THEN
-        time = DBLE(t)*tauMin*ti_ref
+        time = DBLE(timeStep)*tauMin*ti_ref
 
-        fileName = formatFileName(prefix, 'EMField', 'msh', t)
+        fileName = formatFileName(prefix, 'EMField', 'msh', timeStep)
         WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
         OPEN (20, file = path // folder // '/' //  fileName)
 
         CALL writeGmsh2HeaderMesh(20)
 
-        CALL writeGmsh2HeaderNodeData(20, 'Potential (V)', t, time, 1, self%numNodes)
+        CALL writeGmsh2HeaderNodeData(20, 'Potential (V)', timeStep, time, 1, self%numNodes)
         DO n=1, self%numNodes
           WRITE(20, *) n, self%nodes(n)%obj%emData%phi*Volt_ref
         END DO
         CALL writeGmsh2FooterNodeData(20)
         
-        CALL writeGmsh2HeaderElementData(20, 'Electric Field (V m^-1)', t, time, 3, self%numCells)
+        CALL writeGmsh2HeaderElementData(20, 'Electric Field (V m^-1)', timeStep, time, 3, self%numCells)
         DO e=1, self%numCells
           WRITE(20, *) e+self%numEdges, self%cells(e)%obj%gatherElectricField(Xi)*EF_ref
         END DO
         CALL writeGmsh2FooterElementData(20)
 
-        CALL writeGmsh2HeaderNodeData(20, 'Magnetic Field (T)', t, time, 3, self%numNodes)
+        CALL writeGmsh2HeaderNodeData(20, 'Magnetic Field (T)', timeStep, time, 3, self%numNodes)
         DO n=1, self%numNodes
           WRITE(20, *) n, self%nodes(n)%obj%emData%B * B_ref
         END DO
diff --git a/src/modules/mesh/inout/moduleMeshInoutCommon.f90 b/src/modules/mesh/inout/moduleMeshInoutCommon.f90
index e4a6c72..7dc2e84 100644
--- a/src/modules/mesh/inout/moduleMeshInoutCommon.f90
+++ b/src/modules/mesh/inout/moduleMeshInoutCommon.f90
@@ -3,17 +3,17 @@ MODULE moduleMeshInoutCommon
   CHARACTER(LEN=4):: prefix = 'Step'
 
   CONTAINS
-    PURE FUNCTION formatFileName(prefix, suffix, extension, t) RESULT(fileName)
+    PURE FUNCTION formatFileName(prefix, suffix, extension, timeStep) RESULT(fileName)
       USE moduleOutput
       IMPLICIT NONE
 
       CHARACTER(*), INTENT(in):: prefix, suffix, extension
-      INTEGER, INTENT(in), OPTIONAL:: t
+      INTEGER, INTENT(in), OPTIONAL:: timeStep
       CHARACTER (LEN=iterationDigits):: tString
       CHARACTER(:), ALLOCATABLE:: fileName
 
-      IF (PRESENT(t)) THEN
-        WRITE(tString, iterationFormat) t
+      IF (PRESENT(timeStep)) THEN
+        WRITE(tString, iterationFormat) timeStep
         fileName = prefix //  '_' // tString // '_' // suffix // '.' // extension
 
       ELSE
diff --git a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90
index 184d645..e07db01 100644
--- a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90
+++ b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90
@@ -167,7 +167,7 @@ MODULE moduleMeshInputVTU
       CLASS(meshGeneric), INTENT(inout):: self
       CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
       REAL(8):: r(1:3) !3 generic coordinates
-      INTEGER:: fileID, error, found
+      INTEGER:: fileID
       CHARACTER(LEN=256):: line
       INTEGER:: numNodes, numElements, numEdges
       INTEGER, ALLOCATABLE, DIMENSION(:):: entitiesID, offsets, connectivity, types
diff --git a/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90 b/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90
index da90b6b..81e4bbf 100644
--- a/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90
+++ b/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90
@@ -215,17 +215,16 @@ MODULE moduleMeshOutputVTU
 
     END SUBROUTINE writeEM
 
-    SUBROUTINE writeCollection(fileID, t, fileNameStep, fileNameCollection)
+    SUBROUTINE writeCollection(fileID, fileNameStep, fileNameCollection)
       USE moduleCaseParam
       USE moduleOutput
       USE moduleRefParam
       IMPLICIT NONE
 
       INTEGER:: fileID
-      INTEGER, INTENT(in):: t
       CHARACTER(*):: fileNameStep, fileNameCollection
 
-      IF (t == tInitial) THEN
+      IF (timeStep == tInitial) THEN
         !Create collection file
         WRITE(*, "(6X,A15,A)") "Creating file: ", fileNameCollection
         OPEN (fileID + 1, file = path // folder  // '/' // fileNameCollection)
@@ -237,10 +236,11 @@ MODULE moduleMeshOutputVTU
 
       !Write iteration file in collection
       OPEN (fileID + 1, file = path // folder  // '/' // fileNameCollection, ACCESS='APPEND')
-      WRITE(fileID + 1, "(4X, A, ES20.6E3, A, A, A)") '<DataSet timestep="', DBLE(t)*tauMin*ti_ref,'" file="', fileNameStep,'"/>'
+      WRITE(fileID + 1, "(4X, A, ES20.6E3, A, A, A)") &
+            '<DataSet timestep="', DBLE(timeStep)*tauMin*ti_ref,'" file="', fileNameStep,'"/>'
 
       !Close collection file
-      IF (t == tFinal) THEN
+      IF (timeStep == tFinal) THEN
         WRITE (fileID + 1, "(2X, A)") '</Collection>'
         WRITE (fileID + 1, "(A)") '</VTKFile>'
 
@@ -307,21 +307,21 @@ MODULE moduleMeshOutputVTU
 
     END SUBROUTINE writeAverage
 
-    SUBROUTINE printOutputVTU(self,t)
+    SUBROUTINE printOutputVTU(self)
       USE moduleMesh
       USE moduleSpecies
       USE moduleMeshInoutCommon
+      USE moduleCaseParam, ONLY: timeStep
       IMPLICIT NONE
 
       CLASS(meshParticles), INTENT(in):: self
-      INTEGER, INTENT(in):: t
       INTEGER:: i, fileID
       CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection
 
       fileID = 60
 
       DO i = 1, nSpecies
-        fileName = formatFileName(prefix, species(i)%obj%name, 'vtu', t)
+        fileName = formatFileName(prefix, species(i)%obj%name, 'vtu', timeStep)
         WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
         OPEN (fileID, file = path // folder // '/' //  fileName)
 
@@ -337,28 +337,27 @@ MODULE moduleMeshOutputVTU
 
         !Write collection file for time plotting
         fileNameCollection = formatFileName('Collection', species(i)%obj%name, 'pvd')
-        CALL writeCollection(fileID, t, fileName, filenameCollection)
+        CALL writeCollection(fileID, fileName, filenameCollection)
 
       END DO
 
     END SUBROUTINE printOutputVTU
 
-    SUBROUTINE printCollVTU(self,t)
+    SUBROUTINE printCollVTU(self)
       USE moduleMesh
       USE moduleOutput
       USE moduleMeshInoutCommon
+      USE moduleCaseParam, ONLY: timeStep
       IMPLICIT NONE
 
       CLASS(meshGeneric), INTENT(in):: self
-      INTEGER, INTENT(in):: t
       INTEGER:: fileID
       CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection
-      CHARACTER (LEN=iterationDigits):: tstring
 
       fileID = 62
 
       IF (collOutput) THEN
-        fileName = formatFileName(prefix, 'Collisions', 'vtu', t)
+        fileName = formatFileName(prefix, 'Collisions', 'vtu', timeStep)
         WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
         OPEN (fileID, file = path // folder // '/' //  fileName)
 
@@ -374,26 +373,26 @@ MODULE moduleMeshOutputVTU
 
         !Write collection file for time plotting
         fileNameCollection = formatFileName('Collection', 'Collisions', 'pvd')
-        CALL writeCollection(fileID, t, fileName, filenameCollection)
+        CALL writeCollection(fileID, fileName, filenameCollection)
 
       END IF
 
     END SUBROUTINE printCollVTU
 
-    SUBROUTINE printEMVTU(self, t)
+    SUBROUTINE printEMVTU(self)
       USE moduleMesh
       USE moduleMeshInoutCommon
+      USE moduleCaseParam, ONLY: timeStep
       IMPLICIT NONE
 
       CLASS(meshParticles), INTENT(in):: self
-      INTEGER, INTENT(in):: t
       INTEGER:: fileID
       CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection
 
       fileID = 64
 
       IF (emOutput) THEN
-        fileName = formatFileName(prefix, 'EMField', 'vtu', t)
+        fileName = formatFileName(prefix, 'EMField', 'vtu', timeStep)
         WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
         OPEN (fileID, file = path // folder // '/' //  fileName)
 
@@ -409,7 +408,7 @@ MODULE moduleMeshOutputVTU
 
         !Write collection file for time plotting
         fileNameCollection = formatFileName('Collection', 'EMField', 'pvd')
-        CALL writeCollection(fileID, t, fileName, filenameCollection)
+        CALL writeCollection(fileID, fileName, filenameCollection)
 
       END IF
 
diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90
index ae53aa0..f933dd1 100644
--- a/src/modules/mesh/moduleMesh.f90
+++ b/src/modules/mesh/moduleMesh.f90
@@ -372,10 +372,9 @@ MODULE moduleMesh
     END SUBROUTINE connectMesh_interface
 
     !Prints number of collisions in each cell
-    SUBROUTINE printColl_interface(self, t)
+    SUBROUTINE printColl_interface(self)
       IMPORT meshGeneric
       CLASS(meshGeneric), INTENT(in):: self
-      INTEGER, INTENT(in):: t
 
     END SUBROUTINE printColl_interface
 
@@ -403,18 +402,16 @@ MODULE moduleMesh
 
   ABSTRACT INTERFACE
     !Prints Species data
-    SUBROUTINE printOutput_interface(self, t)
+    SUBROUTINE printOutput_interface(self)
       IMPORT meshParticles
       CLASS(meshParticles), INTENT(in):: self
-      INTEGER, INTENT(in):: t
 
     END SUBROUTINE printOutput_interface
 
     !Prints EM info
-    SUBROUTINE printEM_interface(self, t)
+    SUBROUTINE printEM_interface(self)
       IMPORT meshParticles
       CLASS(meshParticles), INTENT(in):: self
-      INTEGER, INTENT(in):: t
 
     END SUBROUTINE printEM_interface
 
@@ -789,7 +786,7 @@ MODULE moduleMesh
     END FUNCTION findCellBrute
 
     !Computes collisions in element
-    SUBROUTINE doCollisions(self, t)
+    SUBROUTINE doCollisions(self)
       USE moduleCollisions
       USE moduleSpecies
       USE moduleList
@@ -797,10 +794,10 @@ MODULE moduleMesh
       USE moduleRandom
       USE moduleOutput
       USE moduleMath
+      USE moduleCaseParam, ONLY: timeStep
       IMPLICIT NONE
 
       CLASS(meshGeneric), INTENT(inout), TARGET:: self
-      INTEGER, INTENT(in):: t
       INTEGER:: e
       CLASS(meshCell), POINTER:: cell
       INTEGER:: k, i, j
@@ -816,7 +813,7 @@ MODULE moduleMesh
       REAL(8):: rnd_real !Random number for collision
       INTEGER:: rnd_int  !Random number for collision
 
-      IF (MOD(t, everyColl) == 0) THEN
+      IF (MOD(timeStep, everyColl) == 0) THEN
         !Collisions need to be performed in this iteration
         !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(part_i, part_j, partTemp_i, partTemp_j)
         DO e=1, self%numCells
diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90
index 18c7fbb..5b96be2 100644
--- a/src/modules/moduleInject.f90
+++ b/src/modules/moduleInject.f90
@@ -54,7 +54,7 @@ MODULE moduleInject
     INTEGER:: id
     CHARACTER(:), ALLOCATABLE:: name
     REAL(8):: vMod !Velocity (module)
-    REAL(8):: T(1:3) !Temperature
+    REAL(8):: temperature(1:3) !Temperature
     REAL(8):: n(1:3) !Direction of injection
     LOGICAL:: fixDirection !The injection of particles has a fix direction defined by n
     INTEGER:: nParticles !Number of particles to introduce each time step
@@ -76,7 +76,7 @@ MODULE moduleInject
 
   CONTAINS
     !Initialize an injection of particles
-    SUBROUTINE initInject(self, i, v, n, T, flow, units, sp, physicalSurface, particlesPerEdge)
+    SUBROUTINE initInject(self, i, v, n, temperature, flow, units, sp, physicalSurface, particlesPerEdge)
       USE moduleMesh
       USE moduleRefParam
       USE moduleConstParam
@@ -87,7 +87,7 @@ MODULE moduleInject
 
       CLASS(injectGeneric), INTENT(inout):: self
       INTEGER, INTENT(in):: i
-      REAL(8), INTENT(in):: v, n(1:3), T(1:3)
+      REAL(8), INTENT(in):: v, n(1:3), temperature(1:3)
       INTEGER, INTENT(in):: sp, physicalSurface, particlesPerEdge
       REAL(8):: tauInject
       REAL(8), INTENT(in):: flow
@@ -97,10 +97,10 @@ MODULE moduleInject
       INTEGER:: nVolColl
       REAL(8):: fluxPerStep = 0.D0
 
-      self%id      =  i
-      self%vMod    =  v / v_ref
-      self%n       =  n / NORM2(n)
-      self%T       =  T / T_ref
+      self%id          =  i
+      self%vMod        =  v / v_ref
+      self%n           =  n / NORM2(n)
+      self%temperature =  temperature / T_ref
       !Gets the edge elements from which particles are injected
       DO e = 1, mesh%numEdges
         phSurface(e) = mesh%edges(e)%obj%physicalSurface
@@ -232,23 +232,23 @@ MODULE moduleInject
 
     END SUBROUTINE doInjects
 
-    SUBROUTINE initVelDistMaxwellian(velDist, T, m)
+    SUBROUTINE initVelDistMaxwellian(velDist, temperature, m)
       IMPLICIT NONE
 
       CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
-      REAL(8), INTENT(in):: T, m
+      REAL(8), INTENT(in):: temperature, m
 
-      velDist = velDistMaxwellian(vTh = DSQRT(T/m))
+      velDist = velDistMaxwellian(vTh = DSQRT(temperature/m))
 
     END SUBROUTINE initVelDistMaxwellian
 
-    SUBROUTINE initVelDistHalfMaxwellian(velDist, T, m)
+    SUBROUTINE initVelDistHalfMaxwellian(velDist, temperature, m)
       IMPLICIT NONE
 
       CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
-      REAL(8), INTENT(in):: T, m
+      REAL(8), INTENT(in):: temperature, m
 
-      velDist = velDistHalfMaxwellian(vTh = DSQRT(T/m))
+      velDist = velDistHalfMaxwellian(vTh = DSQRT(temperature/m))
 
     END SUBROUTINE initVelDistHalfMaxwellian
 
diff --git a/src/modules/moduleProbe.f90 b/src/modules/moduleProbe.f90
index c29eeda..754d56a 100644
--- a/src/modules/moduleProbe.f90
+++ b/src/modules/moduleProbe.f90
@@ -27,7 +27,7 @@ MODULE moduleProbe
 
   CONTAINS
     !Functions for probeDistFunc type
-    SUBROUTINE init(self, id, speciesName, r, v1, v2, v3, points, timeStep)
+    SUBROUTINE init(self, id, speciesName, r, v1, v2, v3, points, everyTimeStep)
       USE moduleCaseParam
       USE moduleRefParam
       USE moduleSpecies
@@ -41,7 +41,7 @@ MODULE moduleProbe
       REAL(8), INTENT(in):: r(1:3)
       REAL(8), INTENT(in):: v1(1:2), v2(1:2), v3(1:2)
       INTEGER, INTENT(in):: points(1:3)
-      REAL(8), INTENT(in):: timeStep
+      REAL(8), INTENT(in):: everyTimeStep
       INTEGER:: sp, i
       REAL(8):: dv(1:3)
 
@@ -91,11 +91,11 @@ MODULE moduleProbe
                       1:self%nv(3)))
 
       !Number of iterations between output
-      IF (timeStep == 0.D0) THEN
+      IF (everyTimeStep == 0.D0) THEN
         self%every = 1
 
       ELSE
-        self%every = NINT(timeStep/ tauMin / ti_ref)
+        self%every = NINT(everyTimeStep/ tauMin / ti_ref)
 
       END IF
       
@@ -189,13 +189,13 @@ MODULE moduleProbe
           
     END SUBROUTINE calculate
 
-    SUBROUTINE output(self, t)
+    SUBROUTINE output(self)
       USE moduleOutput
       USE moduleRefParam
+      USE moduleCaseParam, ONLY: timeStep
       IMPLICIT NONE
 
       CLASS(probeDistFunc), INTENT(inout):: self
-      INTEGER, INTENT(in):: t
       CHARACTER (LEN=iterationDigits):: tstring
       CHARACTER (LEN=3):: pstring
       CHARACTER(:), ALLOCATABLE:: filename
@@ -204,14 +204,14 @@ MODULE moduleProbe
       !Divide by the velocity cube volume
       self%f = self%f * self%dvInv
 
-      WRITE(tstring, iterationFormat) t
+      WRITE(tstring, iterationFormat) timeStep
       WRITE(pstring, "(I3.3)") self%id
       fileName='Probe_' // tstring// '_f_' // pstring // '.dat'
       WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
       OPEN (10, file = path // folder // '/' //  fileName)
       WRITE(10, "(A1, 1X, A)")               "# ", self%species%name
       WRITE(10, "(A6, 3(ES15.6E3), A2)") "# r = ", self%r(:)*L_ref, " m"
-      WRITE(10, "(A6, ES15.6E3, A2)")    "# t = ", REAL(t)*tauMin*ti_ref, " s"
+      WRITE(10, "(A6, ES15.6E3, A2)")    "# t = ", REAL(timeStep)*tauMin*ti_ref, " s"
       WRITE(10, "(A1, A19, 3(A20))")     "#", "v1 (m s^-1)", "v2 (m s^-1)", "v3 (m s^-1)", "f"
       DO i = 1, self%nv(1)
         DO j = 1, self%nv(2)
@@ -252,15 +252,15 @@ MODULE moduleProbe
 
     END SUBROUTINE doProbes
 
-    SUBROUTINE outputProbes(t)
+    SUBROUTINE outputProbes()
+      USE moduleCaseParam, ONLY: timeStep
       IMPLICIT NONE
 
-      INTEGER, INTENT(in):: t
       INTEGER:: i
 
       DO i = 1, nProbes
         IF (probe(i)%update) THEN
-          CALL probe(i)%output(t)
+          CALL probe(i)%output()
 
         END IF
 
@@ -268,15 +268,15 @@ MODULE moduleProbe
 
     END SUBROUTINE outputProbes
 
-    SUBROUTINE resetProbes(t)
+    SUBROUTINE resetProbes()
+      USE moduleCaseParam, ONLY: timeStep
       IMPLICIT NONE
 
-      INTEGER, INTENT(in):: t
       INTEGER:: i
 
       DO i = 1, nProbes
         probe(i)%f = 0.D0
-        probe(i)%update = t == tFinal .OR. t == tInitial .OR. MOD(t, probe(i)%every) == 0
+        probe(i)%update = timeStep == tFinal .OR. timeStep == tInitial .OR. MOD(timeStep, probe(i)%every) == 0
 
       END DO
 
diff --git a/src/modules/output/moduleOutput.f90 b/src/modules/output/moduleOutput.f90
index 18fbb7f..e6dc91f 100644
--- a/src/modules/output/moduleOutput.f90
+++ b/src/modules/output/moduleOutput.f90
@@ -160,12 +160,12 @@ MODULE moduleOutput
 
     END SUBROUTINE calculateOutput
 
-    SUBROUTINE printTime(t, first)
+    SUBROUTINE printTime(first)
       USE moduleSpecies
       USE moduleCompTime
+      USE moduleCaseParam, ONLY: timeStep
       IMPLICIT NONE
 
-      INTEGER, INTENT(in):: t
       LOGICAL, INTENT(in), OPTIONAL:: first
       CHARACTER(:), ALLOCATABLE:: fileName
 
@@ -187,7 +187,7 @@ MODULE moduleOutput
 
         OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write')
 
-        WRITE (20, "(I10, I10, 7(ES20.6E3))") t, nPartOld, tStep, tPush, tReset, tColl, tCoul, tWeight, tEMField
+        WRITE (20, "(I10, I10, 7(ES20.6E3))") timeStep, nPartOld, tStep, tPush, tReset, tColl, tCoul, tWeight, tEMField
 
         CLOSE(20)
 
diff --git a/src/modules/solver/electromagnetic/moduleEM.f90 b/src/modules/solver/electromagnetic/moduleEM.f90
index a141eac..30a23ad 100644
--- a/src/modules/solver/electromagnetic/moduleEM.f90
+++ b/src/modules/solver/electromagnetic/moduleEM.f90
@@ -1,6 +1,7 @@
 !Module to solve the electromagnetic field
 MODULE moduleEM
   USE moduleMesh
+  USE moduleTable
   IMPLICIT NONE
 
   ! Array of pointers to nodes.
@@ -11,57 +12,155 @@ MODULE moduleEM
 
   END TYPE meshNodePointer
 
-
-  ! TODO: Make this a derived type.
-  TYPE:: boundaryEM
-    CHARACTER(:), ALLOCATABLE:: typeEM
-    INTEGER:: physicalSurface
-    REAL(8):: potential
+  ! Generic type for electromagnetic boundary conditions
+  TYPE, PUBLIC, ABSTRACT:: boundaryEMGeneric
+    INTEGER:: nNodes
     TYPE(meshNodePointer), ALLOCATABLE:: nodes(:)
 
     CONTAINS
-      PROCEDURE, PASS:: apply
+      PROCEDURE(applyEM_interface), DEFERRED, PASS:: apply
       !PROCEDURE, PASS:: update !only for time dependent boundary conditions or maybe change apply????? That might be better.
 
-  END TYPE boundaryEM
+  END TYPE boundaryEMGeneric
+
+  ABSTRACT INTERFACE
+    ! Apply boundary condition to the load vector for the Poission equation
+    SUBROUTINE applyEM_interface(self, vectorF)
+      IMPORT boundaryEMGeneric
+      CLASS(boundaryEMGeneric), INTENT(in):: self
+      REAL(8), INTENT(inout):: vectorF(:)
+
+    END SUBROUTINE applyEM_interface
+
+  END INTERFACE
+
+  TYPE, EXTENDS(boundaryEMGeneric):: boundaryEMDirichlet
+    REAL(8):: potential
+
+    CONTAINS
+      ! boundaryEMGeneric DEFERRED PROCEDURES
+      PROCEDURE, PASS:: apply => applyDirichlet
+
+  END TYPE boundaryEMDirichlet
+
+  TYPE, EXTENDS(boundaryEMGeneric):: boundaryEMDirichletTime
+    REAL(8):: potential
+    TYPE(table1D):: temporalProfile
+
+    CONTAINS
+      ! boundaryEMGeneric DEFERRED PROCEDURES
+      PROCEDURE, PASS:: apply => applyDirichletTime
+
+  END TYPE boundaryEMDirichletTime
+
+  ! Container for boundary conditions
+  TYPE:: boundaryEMCont
+    CLASS(boundaryEMGeneric), ALLOCATABLE:: obj
+
+  END TYPE boundaryEMCont
 
   INTEGER:: nBoundaryEM
-  TYPE(boundaryEM), ALLOCATABLE:: boundEM(:)
+  TYPE(boundaryEMCont), ALLOCATABLE:: boundaryEM(:)
 
   !Information of charge and reference parameters for rho vector
   REAL(8), ALLOCATABLE:: qSpecies(:)
 
   CONTAINS
-    !Apply boundary conditions to the K matrix for Poisson's equation
-    SUBROUTINE apply(self, edge)
+    ! Initialize Dirichlet boundary condition
+    SUBROUTINE initDirichlet(self, physicalSurface, potential)
+      USE moduleMesh
+      USE moduleRefParam, ONLY: Volt_ref
+      IMPLICIT NONE
+
+      CLASS(boundaryEMGeneric), ALLOCATABLE, INTENT(out):: self
+      INTEGER:: physicalSurface
+      REAL(8), INTENT(in):: potential
+      CLASS(meshEdge), POINTER:: edge
+      INTEGER, ALLOCATABLE:: nodes(:), nodesEdge(:)
+      INTEGER:: nNodes, nodesNew
+      INTEGER:: e, n
+
+      ! Allocate boundary edge
+      ALLOCATE(boundaryEMDirichlet:: self)
+
+      SELECT TYPE(self)
+      TYPE IS(boundaryEMDirichlet)
+        self%potential = potential / Volt_ref
+
+        !TODO: This is going into a function
+        !Temporal array to hold nodes
+        ALLOCATE(nodes(0))
+
+        ! Loop thorugh the edges and identify those that are part of the boundary
+        DO e = 1, mesh%numEdges
+          edge => mesh%edges(e)%obj
+          IF (edge%physicalSurface == physicalSurface) THEN
+            ! Edge is of the right boundary index
+            ! Get nodes in the edge
+            nNodes = edge%nNodes
+            nodesEdge = edge%getNodes(nNodes)
+            ! Collect all nodes that are not already in the temporal array
+            DO n = 1, nNodes
+              IF (ANY(nodes == nodesEdge(n))) THEN
+                ! Node already in array, skip
+                CYCLE
+
+              ELSE
+                ! If not, add element to array of nodes
+                nodes = [nodes, nodesEdge(n)]
+
+              END IF
+
+            END DO
+
+          END IF
+
+        END DO
+
+        ! Point boundary to nodes
+        nNodes = SIZE(nodes)
+        ALLOCATE(self%nodes(nNodes))
+        DO n = 1, nNodes
+          self%nodes(n)%obj => mesh%nodes(nodes(n))%obj
+
+        END DO
+
+      END SELECT
+
+    END SUBROUTINE initDirichlet
+
+    !Apply Dirichlet boundary condition to the poisson equation
+    SUBROUTINE applyDirichlet(self, vectorF)
       USE moduleMesh
       IMPLICIT NONE
 
-      CLASS(boundaryEM), INTENT(in):: self
-      CLASS(meshEdge):: edge
-      INTEGER:: nNodes
-      INTEGER, ALLOCATABLE:: nodes(:)
-      INTEGER:: n
+      CLASS(boundaryEMDirichlet), INTENT(in):: self
+      REAL(8), INTENT(inout):: vectorF(:)
+      INTEGER:: n, ni
 
-      nNodes = 1
-      nNodes = edge%nNodes
-      nodes  = edge%getNodes(nNodes)
-
-      DO n = 1, nNodes
-        SELECT CASE(self%typeEM)
-        CASE ("dirichlet")
-          mesh%K(nodes(n), :)        = 0.D0
-          mesh%K(nodes(n), nodes(n)) = 1.D0
-        
-          ! TODO: Change this to pointer
-          mesh%nodes(nodes(n))%obj%emData%type = self%typeEM
-          mesh%nodes(nodes(n))%obj%emData%phi  = self%potential
-
-        END SELECT
+      DO n = 1, self%nNodes
+        self%nodes(n)%obj%emData%phi = self%potential
 
       END DO
 
-    END SUBROUTINE apply
+    END SUBROUTINE applyDirichlet
+
+    !Apply Dirichlet boundary condition with time temporal profile
+    SUBROUTINE applyDirichletTime(self, vectorF)
+      USE moduleMesh
+      USE moduleCaseParam, ONLY: timeStep, tauMin
+      IMPLICIT NONE
+
+      CLASS(boundaryEMDirichletTime), INTENT(in):: self
+      REAL(8), INTENT(inout):: vectorF(:)
+      INTEGER:: n, ni
+
+      DO n = 1, self%nNodes
+        self%nodes(n)%obj%emData%phi = self%potential
+
+      END DO
+
+    END SUBROUTINE applyDirichletTime
 
     !Assemble the source vector based on the charge density to solve Poisson's equation
     SUBROUTINE assembleSourceVector(vectorF)
@@ -74,13 +173,14 @@ MODULE moduleEM
       INTEGER, ALLOCATABLE:: nodes(:)
       REAL(8), ALLOCATABLE:: rho(:)
       INTEGER:: nNodes
-      INTEGER:: e, i, ni
+      INTEGER:: e, i, ni, b
       CLASS(meshNode), POINTER:: node
 
       !$OMP SINGLE
       vectorF = 0.D0
       !$OMP END SINGLE
 
+      ! Calculate charge density in each node
       !$OMP DO REDUCTION(+:vectorF)
       DO e = 1, mesh%numCells
         nNodes = mesh%cells(e)%obj%nNodes
@@ -112,18 +212,12 @@ MODULE moduleEM
       !$OMP END DO
 
       !Apply boundary conditions
-      !$OMP DO
-      DO i = 1, mesh%numNodes
-        node => mesh%nodes(i)%obj
-
-        SELECT CASE(node%emData%type)
-        CASE ("dirichlet")
-          vectorF(i) = node%emData%phi
-
-        END SELECT
+      !$OMP SINGLE
+      DO b = 1, nBoundaryEM
+        CALL boundaryEM(b)%obj%apply(vectorF)
 
       END DO
-      !$OMP END DO
+      !$OMP END SINGLE
 
     END SUBROUTINE assembleSourceVector
 
diff --git a/src/modules/solver/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90
index 0b2837f..5928508 100644
--- a/src/modules/solver/moduleSolver.f90
+++ b/src/modules/solver/moduleSolver.f90
@@ -491,47 +491,46 @@ MODULE moduleSolver
     END SUBROUTINE updateParticleCell
 
     !Update the information about if a species needs to be moved this iteration
-    SUBROUTINE updatePushSpecies(self, t)
+    SUBROUTINE updatePushSpecies(self)
       USE moduleSpecies
+      USE moduleCaseparam, ONLY: timeStep
       IMPLICIT NONE
 
       CLASS(solverGeneric), INTENT(inout):: self
-      INTEGER, INTENT(in):: t
       INTEGER:: s
 
       DO s=1, nSpecies
-        self%pusher(s)%pushSpecies = MOD(t, self%pusher(s)%every) == 0
+        self%pusher(s)%pushSpecies = MOD(timeStep, self%pusher(s)%every) == 0
 
       END DO
 
     END SUBROUTINE updatePushSpecies
 
     !Output the different data and information
-    SUBROUTINE doOutput(t)
+    SUBROUTINE doOutput()
       USE moduleMesh
       USE moduleOutput
       USE moduleSpecies
       USE moduleCompTime
       USE moduleProbe
+      USE moduleCaseParam, ONLY: timeStep
       IMPLICIT NONE
 
-      INTEGER, INTENT(in):: t
-
-      CALL outputProbes(t)
+      CALL outputProbes()
 
       counterOutput = counterOutput + 1
       IF (counterOutput >= triggerOutput .OR. &
-          t == tFinal .OR. t == tInitial) THEN
+          timeStep == tFinal .OR. timeStep == tInitial) THEN
 
         !Resets output counter
         counterOutput=0
 
-        CALL mesh%printOutput(t)
-        IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%printColl(t)
-        CALL mesh%printEM(t)
-        WRITE(*, "(5X,A21,I10,A1,I10)") "t/tFinal: ", t, "/", tFinal
+        CALL mesh%printOutput()
+        IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%printColl()
+        CALL mesh%printEM()
+        WRITE(*, "(5X,A21,I10,A1,I10)") "t/tFinal: ", timeStep, "/", tFinal
         WRITE(*, "(5X,A21,I10)")        "Particles: ", nPartOld
-        IF (t == 0) THEN
+        IF (timeStep == 0) THEN
           WRITE(*, "(5X,A21,F8.1,A2)")   "     init time: ",    1.D3*tStep, "ms"
 
         ELSE
@@ -549,34 +548,32 @@ MODULE moduleSolver
 
       counterCPUTime = counterCPUTime + 1
       IF (counterCPUTime >= triggerCPUTime .OR. &
-          t == tFinal .OR. t == tInitial) THEN
+          timeStep == tFinal .OR. timeStep == tInitial) THEN
 
         !Reset CPU Time counter
         counterCPUTime = 0
 
-        CALL printTime(t, t == 0)
+        CALL printTime(timeStep == 0)
 
       END IF
 
       !Output average values
-      IF (useAverage .AND. t == tFinal) THEN
+      IF (useAverage .AND. timeStep == tFinal) THEN
         CALL mesh%printAverage()
 
       END IF
 
     END SUBROUTINE doOutput
 
-    SUBROUTINE doAverage(t)
+    SUBROUTINE doAverage()
       USE moduleAverage
       USE moduleMesh
       IMPLICIT NONE
 
-      INTEGER, INTENT(in):: t
       INTEGER:: tAverage, n
 
-
       IF (useAverage) THEN
-        tAverage = t - tAverageStart
+        tAverage = timeStep - tAverageStart
 
         IF (tAverage == 1) THEN
           !First iteration in which average scheme is used