SUBROUTINE stepon
!
!   Advance one time step
!
   USE basic
   USE constants
   USE fields
   USE beam
   USE maxwsrce
   USE celldiag
   USE neutcol
   USE sort
   Use psupply
   use omp_lib
   IMPLICIT NONE

   INTEGER:: 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 BARRIER



! Sort particles for faster rhscon run time
!  DO i=1,nbspecies
!    IF(modulo(step,it2d) .eq. 0) THEN
!      CALL gridsort(partslist(i),1,partslist(i)%Nploc)
!    END IF
!  END DO

! Assemble right hand side of Poisson equation
   CALL rhscon(partslist)
   if (.not. nlfreezephi) THEN
! Solve Poisson equation
      CALL poisson(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

   !$OMP SINGLE
   ! Calculate main physical quantities
   CALL partdiagnostics
   !$OMP END SINGLE NOWAIT

   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

  

   ! update the power supply voltage if necessary
   call psupply_step(the_ps,partslist,cstep)

   !$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 SINGLE
   ! 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
   !$OMP END SINGLE

END SUBROUTINE stepon