diff --git a/src/modules/mesh/inout/moduleMeshInoutCommon.f90 b/src/modules/mesh/inout/moduleMeshInoutCommon.f90
index 7dc2e84..8956423 100644
--- a/src/modules/mesh/inout/moduleMeshInoutCommon.f90
+++ b/src/modules/mesh/inout/moduleMeshInoutCommon.f90
@@ -2,26 +2,73 @@ MODULE moduleMeshInoutCommon
 
   CHARACTER(LEN=4):: prefix = 'Step'
 
+  integer, parameter, private:: fileID_boundaryParticle = 30
+  integer, parameter, private:: fileID_boundaryEM = 31
+
   CONTAINS
-    PURE FUNCTION formatFileName(prefix, suffix, extension, timeStep) RESULT(fileName)
+    PURE FUNCTION formatFileName(pref, suff, extension, timeStep) RESULT(fileName)
       USE moduleOutput
       IMPLICIT NONE
 
-      CHARACTER(*), INTENT(in):: prefix, suffix, extension
+      CHARACTER(*), INTENT(in):: pref, suff, extension
       INTEGER, INTENT(in), OPTIONAL:: timeStep
       CHARACTER (LEN=iterationDigits):: tString
       CHARACTER(:), ALLOCATABLE:: fileName
 
       IF (PRESENT(timeStep)) THEN
         WRITE(tString, iterationFormat) timeStep
-        fileName = prefix //  '_' // tString // '_' // suffix // '.' // extension
+        fileName = pref //  '_' // tString // '_' // suff // '.' // extension
 
       ELSE
-        fileName = prefix //  '_' // suffix // '.' // extension
+        fileName = pref //  '_' // suff // '.' // extension
 
       END IF
 
     END FUNCTION formatFileName
 
+    pure function returnFilePath(filename) result(completePath)
+      use moduleOutput, only: path, folder
+      implicit none
+
+      character(:), allocatable, intent(in):: fileName
+      character(:), allocatable:: completePath
+
+      completePath = path // folder // '/' // fileName
+
+    end function returnFilePath
+
+    subroutine informFileCreation(filename)
+      implicit none
+
+      character(:), allocatable, intent(in):: fileName
+
+      write(*, "(6X,A15,A)") "Creating file: ", fileName
+
+    end subroutine informFileCreation
+
+    module subroutine boundariesParticle_print()
+      use moduleCaseparam, only: timeStep
+      use moduleOutput, only: path, folder
+      use moduleMesh, only: nBoundariesParticle, boundariesParticle
+      implicit none
+
+      integer:: b
+      character(:), allocatable:: fileName
+
+      fileName = formatFileName(prefix, 'boundariesParticle', 'csv', timeStep)
+      call informFileCreation(fileName)
+      open(fileID_boundaryParticle, file = path // folder // '/' //  fileName)
+
+      do b = 1, nBoundariesParticle
+        if (associated(boundariesParticle(b)%obj%print)) then
+          call boundariesParticle(b)%obj%print(fileID_boundaryParticle)
+
+        end if
+
+      end do
+
+      close(fileID_boundaryParticle)
+        
+    end subroutine boundariesParticle_print
 
 END MODULE moduleMeshInoutCommon
diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90
index 1cb4dc7..4578652 100644
--- a/src/modules/mesh/moduleMesh.f90
+++ b/src/modules/mesh/moduleMesh.f90
@@ -4,7 +4,6 @@ MODULE moduleMesh
   USE moduleOutput
   USE moduleCollisions
   use moduleSpecies, only: nSpecies
-  IMPLICIT NONE
 
   ! MESH ELEMENTS
   !Generic mesh element
@@ -710,10 +709,11 @@ MODULE moduleMesh
     end subroutine updateParticle_interface
 
     ! Update the values of the particle boundary model
-    subroutine printParticle_interface(self)
+    subroutine printParticle_interface(self, fileID)
       import boundaryParticleGeneric
 
-      class(boundaryParticleGeneric), intent(inout):: self
+      class(boundaryParticleGeneric), intent(in):: self
+      integer, intent(in):: fileID
 
     end subroutine printParticle_interface
 
diff --git a/src/modules/mesh/moduleMesh@boundaryParticle.f90 b/src/modules/mesh/moduleMesh@boundaryParticle.f90
index c30ca1d..7073806 100644
--- a/src/modules/mesh/moduleMesh@boundaryParticle.f90
+++ b/src/modules/mesh/moduleMesh@boundaryParticle.f90
@@ -404,16 +404,21 @@ submodule(moduleMesh) boundaryParticle
 
     end subroutine quasiNeutrality_update
 
-    subroutine quasiNeutrality_print(self)
+    subroutine quasiNeutrality_print(self, fileID)
       implicit none
 
-      class(boundaryParticleGeneric), intent(inout):: self
-
-      print *, 'test'
+      class(boundaryParticleGeneric), intent(in):: self
+      integer, intent(in):: fileID
+      integer:: e
 
+      write(fileID, '(A)') self%name
       select type(self)
       type is(boundaryQuasiNeutrality)
-        print*, self%alpha
+        write(fileID, '(A,",",A)'), '"Edge id"', '"alpha"'
+        do e = 1, size(self%edges)
+          write(fileID, '(I0,",",ES0.6E3)'), self%edges(e)%obj%n, self%alpha(self%edges(e)%obj%n)
+
+        end do
 
       end select
 
@@ -470,6 +475,6 @@ submodule(moduleMesh) boundaryParticle
 
       end do
 
-  end subroutine boundariesParticle_update
+    end subroutine boundariesParticle_update
 
 end submodule boundaryParticle
diff --git a/src/modules/solver/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90
index 4ea051c..0748fb9 100644
--- a/src/modules/solver/moduleSolver.f90
+++ b/src/modules/solver/moduleSolver.f90
@@ -521,10 +521,9 @@ MODULE moduleSolver
       USE moduleCompTime
       USE moduleProbe
       USE moduleCaseParam, ONLY: timeStep
+      use moduleMeshInoutCommon, only: boundariesParticle_print
       IMPLICIT NONE
 
-      integer:: b ! TEMPORARY
-
       CALL outputProbes()
 
       counterOutput = counterOutput + 1
@@ -537,6 +536,10 @@ MODULE moduleSolver
         CALL mesh%printOutput()
         IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%printColl()
         CALL mesh%printEM()
+
+        ! Output of boundaries. TODO: Add an if
+        call boundariesParticle_print()
+
         WRITE(*, "(5X,A21,I10,A1,I10)") "t/tFinal: ", timeStep, "/", tFinal
         WRITE(*, "(5X,A21,I10)")        "Particles: ", nPartOld
         IF (timeStep == 0) THEN
@@ -553,17 +556,6 @@ MODULE moduleSolver
         END IF
         WRITE(*,*)
 
-        ! TEMPORARY
-        ! Output of boundaries. TODO: Move this to mesh module
-        do b = 1, nBoundariesParticle
-          if (associated(boundariesParticle(b)%obj%print)) then
-            print *, b
-            call boundariesParticle(b)%obj%print()
-
-          end if
-
-        end do
-
       END IF
 
       counterCPUTime = counterCPUTime + 1