diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90
index 8ab443f..19a4369 100644
--- a/src/modules/init/moduleInput.f90
+++ b/src/modules/init/moduleInput.f90
@@ -979,11 +979,6 @@ MODULE moduleInput
 
             call initDirichletTime(self, config, object)
 
-          case ("neumann")
-            allocate(boundaryEMNeumann:: self)
-
-            call initNeumann(self, config, object)
-
           case ("floating")
             allocate(boundaryEMFloating:: self)
 
diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90
index 5fc6f2e..e257cf2 100644
--- a/src/modules/mesh/moduleMesh.f90
+++ b/src/modules/mesh/moduleMesh.f90
@@ -1056,33 +1056,6 @@ MODULE moduleMesh
 
   end interface
 
-  ! Neumann boundary condition (constant electric field normal to edge)
-  type, extends(boundaryEMGeneric):: boundaryEMNeumann
-    real(8):: electricField ! Electric field normal to the edge
-    contains
-      ! boundaryEMGeneric deferred procedures
-      procedure, pass:: apply => applyNeumann
-
-  end type boundaryEMNeumann
-
-  interface
-    module subroutine initNeumann(self, config, object)
-      use json_module
-
-      class(boundaryEMGeneric), allocatable, intent(inout):: self
-      type(json_file), intent(inout):: config
-      character(:), allocatable, intent(in):: object
-
-    end subroutine initNeumann
-
-    module subroutine applyNeumann(self, vectorF)
-      class(boundaryEMNeumann), intent(in):: self
-      real(8), intent(inout):: vectorF(:)
-
-    end subroutine applyNeumann
-
-  end interface
-
   ! Floating: Floating surface, ie, zero net current
   type, extends(boundaryEMGeneric):: boundaryEMFloating
     real(8):: potential
diff --git a/src/modules/mesh/moduleMesh@boundaryEM.f90 b/src/modules/mesh/moduleMesh@boundaryEM.f90
index 891e34b..26363b5 100644
--- a/src/modules/mesh/moduleMesh@boundaryEM.f90
+++ b/src/modules/mesh/moduleMesh@boundaryEM.f90
@@ -30,6 +30,7 @@ submodule(moduleMesh) boundaryEM
 
     ! Apply
     module SUBROUTINE applyDirichlet(self, vectorF)
+      USE moduleMesh
       IMPLICIT NONE
 
       CLASS(boundaryEMDirichlet), INTENT(in):: self
@@ -37,11 +38,8 @@ submodule(moduleMesh) boundaryEM
       integer:: n
 
       DO n = 1, self%nNodes
-        associate(node => self%nodes(n)%obj)
-          node%emData%phi  = self%potential
-          vectorF(node%n)  = node%emData%phi
-
-        end associate
+        self%nodes(n)%obj%emData%phi  = self%potential
+        vectorF(self%nodes(n)%obj%n)  = self%nodes(n)%obj%emData%phi
 
       END DO
 
@@ -125,52 +123,6 @@ submodule(moduleMesh) boundaryEM
 
     end subroutine updateDirichletTime
 
-    ! Neumann constant value
-    ! Init
-    module subroutine initNeumann(self, config, object)
-      use json_module
-      use moduleRefParam, only: EF_ref
-      use moduleErrors, only: criticalError
-      implicit none
-
-      class(boundaryEMGeneric), allocatable, intent(inout):: self
-      type(json_file), intent(inout):: config
-      character(:), allocatable, intent(in):: object
-      real(8):: electricField
-      logical:: found
-
-      select type(self)
-      type is(boundaryEMNeumann)
-        call config%get(object // '.electricField', electricField, found)
-        if (.not. found) then
-          call criticalError('Required parameter "electricField" for Neumann boundary condition not found', 'readBoundaryEM')
-
-        end if
-
-        self%electricField = electricField / EF_ref
-
-      end select
-
-    end subroutine initNeumann
-
-    ! Apply
-    module subroutine applyNeumann(self, vectorF)
-      implicit none
-
-      class(boundaryEMNeumann), intent(in):: self
-      real(8), intent(inout):: vectorF(:)
-      integer:: n
-
-      do n = 1, self%nNodes
-        associate(node => self%nodes(n)%obj)
-          vectorF(node%n) = vectorF(node%n) + self%electricField
-
-        end associate
-
-      end do
-
-    end subroutine applyNeumann
-
     ! Floating surface
     ! Init
     module subroutine initFloating(self, config, object)
@@ -272,6 +224,9 @@ submodule(moduleMesh) boundaryEM
 
             mom_center = edge%gatherF(edge%centerXi(), edge%nNodes, mom_nodes)
 
+            ! Only account for charge exiting the surface
+            mom_center = max(mom_center, 0.0d0)
+            
             edgeDensityCurrent = edgeDensityCurrent + qSpecies(s) * mom_center
 
           end do