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

! 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

   ! We sort the particles according to their linear index
!IF(modulo(step,100) .eq. 0) THEN
!   DO i=1,nbspecies
!      ! Boundary conditions for plasma particles outside the plasma region
!      call gridsort(partslist(i), 1, partslist(i)%Nploc)
!   END DO
!end if

   
   ! 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
      !!$OMP MASTER
      CALL poisson(splrz,reducedsol) 
      !!$OMP END MASTER
   end if

   !$OMP PARALLEL DEFAULT(shared), private(i)

   DO i=1,nbspecies
      ! Compute the magnetic field at the particle position
      call comp_mag_p(partslist(i))
   END DO
   if (.not. nlfreezephi) THEN
      call poisson_com(splrz,reducedsol)
      CALL Update_phi(splrz)
   end if


   DO i=1,nbspecies
      ! Compute the electric field at the particle position
      CALL EFieldscompatparts(partslist(i))
   END DO
   !$OMP BARRIER

   DO i=1,nbspecies
      ! 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

   !$OMP BARRIER

   !                    Cell diag quantities
   IF(modulo(step,itcelldiag).eq. 0 .or. nlend) THEN
      CALL celldiag_save(time, fidres)
   END IF

   ! 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)
      Do i=1,nbspecies
         if( partslist(i)%Nploc*2 .lt. size(partslist(i)%pot,1)) then
            call change_parts_allocation(partslist(i),-(size(partslist(i)%pot,1)-partslist(i)%Nploc*2))
         end if
      end do
   END IF

END SUBROUTINE stepon