diff --git a/src/inital.f90 b/src/inital.f90
index 769a2e7..5b29cc5 100644
--- a/src/inital.f90
+++ b/src/inital.f90
@@ -1,67 +1,67 @@
 SUBROUTINE inital
 !
 USE basic
 USE beam
 USE fields
 USE mpihelper
 USE maxwsrce
 Use geometry
 Use neutcol
 USE splinebound
 !
 !   Set initial conditions
 !
   IMPLICIT NONE
   INTEGER:: i, nbbounds
 !
 !________________________________________________________________________________
   IF(mpirank .eq. 0) WRITE(*,'(a/)') '=== Set initial conditions ==='
 !________________________________________________________________________________
 ! 
   ! Init Electric and Magnetic Fields
   ALLOCATE(partslist(nbspecies))
   CALL fields_init
   call timera(0, "read_geom")
   call read_geom(lu_in, rnorm, splrz, Potinn, Potout)
   call timera(1, "read_geom")
   call mag_init
   
   ! resize nblost array to adapt for correct number of boundaries
   nbbounds=2
   if(the_domain%nbsplines .gt. 0)  nbbounds=the_domain%nbsplines
   Do i=1,nbspecies
     if( allocated(partslist(i)%nblost)) deallocate(partslist(i)%nblost)
     allocate(partslist(i)%nblost(4+nbbounds))
     partslist(i)%nblost=0
   end do
  
   call fields_start
 
   CALL load_parts ! will call the localisation
 
   qnorm=abs(partslist(1)%weight*partslist(1)%q)
   
   IF(mpisize .gt. 1) THEN
     CALL calc_Zbounds(partslist(1), Zbounds, femorder)
   END IF
 
   Do i=1,size(partslist)
     CALL keep_mpi_self_parts(partslist(i), Zbounds)
   END DO
 
  
   CALL fields_comm_init(Zbounds)
   CALL rhscon(partslist)
-  CALL poisson_mpi(splrz)
+  CALL poisson(splrz)
   !$OMP PARALLEL
   CALL Update_phi(splrz)
   !$OMP END PARALLEL
   CALL EFieldscompatparts(partslist(1))
   CALL adapt_vinit(partslist(1))
   
   
   
 
 
 !________________________________________________________________________________
 END SUBROUTINE inital
diff --git a/src/resume.f90 b/src/resume.f90
index c8845ea..47726e2 100644
--- a/src/resume.f90
+++ b/src/resume.f90
@@ -1,86 +1,86 @@
 SUBROUTINE resume
 !
 !   Resume from previous run
 !
   Use beam
   Use basic
   Use fields
   Use sort
   Use maxwsrce
   Use geometry
   Use neutcol
   IMPLICIT NONE
  
 !
 !   Local vars and arrays
  INTEGER:: i, nbbounds
 !________________________________________________________________________________
   WRITE(*,'(a)') '   Resume from previous run'
 !________________________________________________________________________________
 ! 
 !   Open and read initial conditions from restart file
   CALL chkrst(0)
 
 ! If we want to start a new run with a new file
   IF (newres) THEN
     ! we start time and step from 0
     cstep=0
     time=0
   END IF
 
   qnorm=abs(partslist(1)%weight*partslist(1)%q)
   
   CALL fields_init
   call timera(0, "read_geom")
   call read_geom(lu_in, rnorm, splrz, Potinn, Potout)
   call timera(1, "read_geom")
 ! Initialize the magnetic field and the electric field solver
   call mag_init
   CALL fields_start
 
   ! resize nblost array to adapt for correct number of boundaries
   nbbounds=2
   if(the_domain%nbsplines .gt. 0)  nbbounds=the_domain%nbsplines
   Do i=1,nbspecies
     if( allocated(partslist(i)%nblost)) deallocate(partslist(i)%nblost)
     allocate(partslist(i)%nblost(4+nbbounds))
     partslist(i)%nblost=0
   end do
 
   call boundary_loss(partslist(1))
 
   ! Add the requested test particles and read them from input files
   IF (nbaddtestspecies .gt. 0) THEN
     Do i=1,nbaddtestspecies
       CALL load_part_file(partslist(nbspecies+i),addedtestspecfile(i))
       call boundary_loss(partslist(nbspecies+i))
     END DO
     nbspecies=nbspecies+nbaddtestspecies
   END IF
 
   IF(mpisize .gt. 1) THEN
     CALL calc_Zbounds(partslist(1), Zbounds, femorder)
     DO i=1,nbspecies
       CALL keep_mpi_self_parts(partslist(i),Zbounds)
       call collectparts(partslist(i))
     END DO
   END IF
 
   
 
   CALL bound(partslist(1))
   call boundary_loss(partslist(1))
   ! Allocate the variables for MPI communications
   CALL fields_comm_init(Zbounds)
   ! Solve Poisson for the initial conditions
   CALL rhscon(partslist)
-  CALL poisson_mpi(splrz)
+  CALL poisson(splrz)
   !$OMP PARALLEL
   CALL Update_phi(splrz)
   !$OMP END PARALLEL
   ! Compute the fields at the particles positions
   CALL EFieldscompatparts(partslist(1))
   if(mpirank .eq. 0) WRITE(*,*) "Initial forces computed"
   
 !
 END SUBROUTINE resume
diff --git a/src/stepon.f90 b/src/stepon.f90
index a047dcc..d8dbbcc 100644
--- a/src/stepon.f90
+++ b/src/stepon.f90
@@ -1,123 +1,123 @@
 SUBROUTINE stepon
 !
 !   Advance one time step
 !
    USE basic
    USE constants
    USE fields
    USE maxwsrce
    USE celldiag
    USE neutcol
    USE sort
    Use psupply
    use omp_lib
    USE beam
    IMPLICIT NONE
 
    INTEGER:: i
    
    !$OMP PARALLEL DEFAULT(shared), private(i)
    DO i=1,nbspecies
       ! Boundary conditions for plasma particles outside the plasma region
       CALL bound(partslist(i))
    END DO
    !$OMP BARRIER
    DO i=1,nbspecies
       ! Localisation of particles in cells (calculation of the r and z indices)
       call boundary_loss(partslist(i))
    END DO
    !$OMP BARRIER
 
 
    !                    Cell diag quantities
    IF(modulo(step,itcelldiag).eq. 0 .or. nlend) THEN
       CALL celldiag_save(time, fidres)
    END IF
 ! We compute collisions on the main particles
    IF(modulo(step,itcol).eq. 0) THEN
       CALL neutcol_step(partslist)
    END IF
    !$OMP BARRIER
 
    
 !$OMP SINGLE
 ! The particles are injected by the source
    CALL maxwsrce_inject(time)
 !$OMP END SINGLE
 !$OMP END PARALLEL
 
    ! update the power supply voltage if necessary
    call psupply_step(the_ps,partslist,cstep)
 
 !$OMP PARALLEL
 ! Assemble right hand side of Poisson equation
    CALL rhscon(partslist)
 !$OMP END PARALLEL
 
    if (.not. nlfreezephi) THEN
 ! Solve Poisson equation
-      CALL poisson_mpi(splrz) 
+      CALL poisson(splrz) 
    end if
 
-   !$OMP PARALLEL
+   !$OMP PARALLEL DEFAULT(shared), private(i)
    if (.not. nlfreezephi) THEN
       CALL Update_phi(splrz)
    end if
 
 
    DO i=1,nbspecies
       ! Compute the electric field at the particle position
       CALL EFieldscompatparts(partslist(i))
 
       ! Compute the magnetic field at the particle position
       call comp_mag_p(partslist(i))
    END DO
    !$OMP BARRIER
 
    DO i=1,nbspecies
       ! Solve Newton eq. and advance velocity by delta t
       CALL comp_velocity(partslist(i))
       !$OMP SINGLE
       ! Compute the energy of added particles
       CALL calc_newparts_energy(partslist(i))
       !$OMP END SINGLE NOWAIT
    END DO
 
    !$OMP BARRIER
 
    ! Calculate main physical quantities
    CALL partdiagnostics
 
    IF (modulo(step,it2d).eq. 0 .or. nlend) THEN
       Do i=1,nbspecies
          if(partslist(i)%calc_moments) CALL momentsdiag(partslist(i))
       End do
    END IF
-
+   !$OMP BARRIER
 
   
  
    !$OMP MASTER
    ! Save variables to file
    CALL diagnose(step)
    !$OMP END MASTER
 
    !IF (modulo(step,itparts).eq. 0 .or. modulo(step,ittracer).eq. 0 .or.  modulo(step,itrestart).eq. 0 .or. nlend) THEN
    !$OMP BARRIER
    !END IF
 
    Do i=1,nbspecies
       ! Calculate new positions of particles at time t+delta t
       CALL push(partslist(i))
    END DO
    !$OMP END PARALLEL
 
 
    ! We recalculate the mpi axial boundaries and we adapt them if necessary
    IF(modulo(step,100) .eq. 0) THEN
       CALL calc_Zbounds(partslist(1),Zbounds, femorder)
       CALL fields_comm_init(Zbounds)
       CALL maxwsrce_calcfreq(Zbounds)
    END IF
 
 END SUBROUTINE stepon