diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90
index 19a4369..8ab443f 100644
--- a/src/modules/init/moduleInput.f90
+++ b/src/modules/init/moduleInput.f90
@@ -979,6 +979,11 @@ 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 e257cf2..5fc6f2e 100644
--- a/src/modules/mesh/moduleMesh.f90
+++ b/src/modules/mesh/moduleMesh.f90
@@ -1056,6 +1056,33 @@ 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 26363b5..891e34b 100644
--- a/src/modules/mesh/moduleMesh@boundaryEM.f90
+++ b/src/modules/mesh/moduleMesh@boundaryEM.f90
@@ -30,7 +30,6 @@ submodule(moduleMesh) boundaryEM
 
     ! Apply
     module SUBROUTINE applyDirichlet(self, vectorF)
-      USE moduleMesh
       IMPLICIT NONE
 
       CLASS(boundaryEMDirichlet), INTENT(in):: self
@@ -38,8 +37,11 @@ submodule(moduleMesh) boundaryEM
       integer:: n
 
       DO n = 1, self%nNodes
-        self%nodes(n)%obj%emData%phi  = self%potential
-        vectorF(self%nodes(n)%obj%n)  = self%nodes(n)%obj%emData%phi
+        associate(node => self%nodes(n)%obj)
+          node%emData%phi  = self%potential
+          vectorF(node%n)  = node%emData%phi
+
+        end associate
 
       END DO
 
@@ -123,6 +125,52 @@ 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)
@@ -224,9 +272,6 @@ 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