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)
  !$OMP PARALLEL
  CALL poisson(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