I think I have it, and I think that the averaring might not be required

src/modules/init/moduleInput.f90

@@ -1351,6 +1351,8 @@ MODULE moduleInput
                 bound%nEdges = size(bound%edges)
                 allocate(bound%electricField(1:bound%nEdges))
                 bound%electricField = 0.0d0
+                allocate(bound%deltaElectricField(1:bound%nEdges))
+                bound%deltaElectricField = 0.0d0
 
               end if
 

src/modules/mesh/moduleMesh.f90

@@ -1118,6 +1118,7 @@ MODULE moduleMesh
     integer:: nEdges
     type(meshEdgePointer), allocatable:: edges(:) ! Edges included in the boundary
     real(8), allocatable:: electricField(:) ! Electric field normal to the edge that must be applied with a Neumann boundary condition
+    real(8), allocatable:: deltaElectricField(:) ! Accumulated change in E
     integer:: counter
     contains
       procedure, pass:: apply => applyFreeCurrent

src/modules/mesh/moduleMesh@boundaryEM.f90

@@ -367,6 +367,7 @@ submodule(moduleMesh) boundaryEM
     ! Update
     subroutine updateFreeCurrent(self)
       use moduleCaseParam, only: tauMin
+      use moduleRefParam, only: Vol_ref, n_ref
       implicit none
 
       class(boundaryEMGeneric), intent(inout):: self
@@ -376,13 +377,19 @@ submodule(moduleMesh) boundaryEM
       real(8), allocatable:: den_nodes(:)
       class(meshNode), pointer:: node
       real(8):: mom_center, den_center, charge_center
-      real(8):: vel_avg_center
+      real(8):: l_c ! Characteristic length
+      integer:: every
+      real(8):: every_inv
+
+      ! TODO: Make this an input parameter
+      every = 10
+      every_inv = 1.0d0 / real(every)
 
       select type(self)
       type is(boundaryEMFreeCurrent)
         self%counter = self%counter + 1
-        do e=1, self%nEdges
 
+        do e=1, self%nEdges
           mom_center    = 0.0d0
           den_center    = 0.0d0
           charge_center = 0.0d0
@@ -392,16 +399,14 @@ submodule(moduleMesh) boundaryEM
             nodes = edge%getNodes(edge%nNodes)
             allocate(mom_nodes(1:edge%nNodes))
             allocate(den_nodes(1:edge%nNodes))
-            mom_nodes = 0.0d0
-            den_nodes = 0.0d0
 
             do s = 1, nSpecies
               do n = 1, self%nNodes
                 node => mesh%nodes(nodes(n))%obj
 
                 ! Minus sign as we look at the values exiting the surface (opposite to the normal) 
-                mom_nodes(n) = - dot_product(node%output(s)%mom, edge%normal)
-                den_nodes(n) =   node%output(s)%den
+                mom_nodes(n) = - dot_product(node%output(s)%mom, edge%normal)/(node%v*Vol_ref*n_ref)
+                den_nodes(n) =   node%output(s)%den/(node%v*Vol_ref*n_ref)
 
               end do
 
@@ -412,20 +417,26 @@ submodule(moduleMesh) boundaryEM
             end do
 
             if (den_center > 1.0e-10) then
-              vel_avg_center = mom_center / den_center
-            else
-              vel_avg_center = 0.0d0
+              self%deltaElectricField(e) = self%deltaElectricField(e) + mom_center * charge_center / den_center * tauMin
 
             end if
 
+            deallocate(nodes, mom_nodes, den_nodes)
+
           end associate
 
-          self%electricField(e) = self%electricField(e) - vel_avg_center*charge_center*tauMin
-
-          deallocate(nodes, mom_nodes, den_nodes)
 
         end do
 
+        if (self%counter == every) then
+
+          self%electricField = self%electricField - self%deltaElectricField*every_inv
+
+          ! reset for the next iteration
+          self%counter = 0
+          self%deltaElectricField = 0.0d0
+
+        end if
 
       end select