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

   INTEGER:: i

   DO i=1,nbspecies

      ! Boundary conditions for plasma particles outside the plasma region
      CALL bound(partslist(i))
      ! Localisation of particles in cells (calculation of the r and z indices)
      call boundary_loss(partslist(i))
   END DO


   !                    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
   

! The particles are injected by the source
   CALL maxwsrce_inject(time)

! 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))

      ! Solve Newton eq. and advance velocity by delta t
      CALL comp_velocity(partslist(i))

      ! Compute the energy of added particles
      CALL calc_newparts_energy(partslist(i))
   END DO

   ! 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

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

   ! Save variables to file
   CALL diagnose(step)

   Do i=1,nbspecies
      ! Calculate new positions of particles at time t+delta t
      CALL push(partslist(i))
   END DO



   ! We recalculate the mpi axial boundaries and we adapt them if necessary
   IF(modulo(step,50) .eq. 0) THEN
      CALL calc_Zbounds(partslist(1),Zbounds, femorder)
      CALL fields_comm_init(Zbounds)
      CALL maxwsrce_calcfreq(Zbounds)
   END IF

END SUBROUTINE stepon