subroutine load_particle
    use read_in_parameters, only: irun, ifullk
    use global_parameters, only: mype, gtcout
    use precision
    implicit none
    if(mype==0) then
        write(gtcout, *) "============================================="
        write(gtcout, *) "Load Particle Starts:"
        write(gtcout, *) "============================================="
        write(gtcout, *) "" !this is a new line
        flush(gtcout)
    endif
    if(mype==0)write(gtcout,*)"ifullk=",ifullk
    if(ifullk==0)then
        call loadi
    else
        call loadfki
    endif
    call loade
    if(irun>0)call restart_load
end subroutine load_particle
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine loadi
    use precision
    use read_in_parameters, only: track_particles
    use particle_array, only: mimax, zion, zion0, zion1, wzion, wpion, jtion0, jtion1, rhogyroion,&
        wtion0, wtion1, jpion, markeri, pmarki, nonunii, marknumi, micore, misol,pressureipara,&
        pressureiperp, densityi, densityi0, zonali, zonalci, flowi, rhogyroion0, pfluxi, rdtemi,&
        data1di,mpdata1d
    use grids_setup, only: mpsi, mgrid
    use grids_eqdata, only: meshni, meshti1d
    use equilibrium, only: ti_Eq
    use read_in_parameters, only: aion, qion, ngyroi, irun
    use global_parameters, only: mi, nparami,mype,gtcout
    implicit none
    integer ierror, mtest
interface
    subroutine gkload(mp, mpmax, zpart, zpart0, wzpart, wppart, jtpart0, jtpart1,&
    wtpart0, wtpart1, jppart,rhogyropart, t_Eq, apart, qpart, marker,&
    pmark, meshn, mesht1d, nonuni, marknum, mpcore, mpsol, ngyro)
    use precision
! dummy parameters
    real(lk), dimension(:,:,:) :: t_Eq
    real(lk), dimension(:,:) ::zpart, zpart0, wppart, wtpart0,wtpart1,rhogyropart
    integer, dimension(:,:) :: jtpart0, jtpart1, jppart
    real(lk), dimension(:) :: wzpart
    real(lk), dimension(:) :: marker, meshn, mesht1d, nonuni, pmark
    integer, dimension(:) :: marknum
    integer mp, mpmax, mpcore, mpsol, ngyro
    real(lk) apart, qpart
    end subroutine gkload
    subroutine particle_grid_array_alloc(nparam, ngyro, mp, mpmax, marknum, jtpart0, jtpart1, jppart,&
        wzpart, pmark, zonal, zonalc, marker, nonuni, zpart, zpart0, wppart, wtpart0, wtpart1, density,density0,&
        flow, pressurepara, pressureperp,rhogyropart,rhogyropart0,pfluxpart,rdtempart,data1dpart)
        use precision
        implicit none
!dummy parameters
        integer :: nparam, ngyro, mp, mpmax
        integer, dimension(:),allocatable :: marknum
        integer,dimension(:,:),allocatable :: jtpart0,jtpart1,jppart
        real(lk),dimension(:),allocatable :: wzpart,pmark,zonal,zonalc,marker,nonuni
        real(lk),dimension(:,:),allocatable :: zpart,zpart0,wppart,wtpart0,wtpart1,&
            density,density0,flow,pressurepara,pressureperp,rhogyropart,rhogyropart0,&
            data1dpart
        real(lk),dimension(:),allocatable :: pfluxpart,rdtempart
    end subroutine particle_grid_array_alloc
end interface
! initialize particle position and velocity: thermal ions
    allocate(pmarki(mpsi),markeri(mgrid),densityi(0:1,mgrid),densityi0(0:1,mgrid),flowi(0:1,mgrid),&
             pressureipara(0:1,mgrid),pressureiperp(0:1,mgrid),zonali(mpsi),zonalci(mpsi),&
             nonunii(mgrid),marknumi(mgrid),pfluxi(mpsi),rdtemi(mpsi),data1di(mpsi,mpdata1d),&
             STAT=mtest)
    if(mtest /= 0) then
       write(0,*)mype,'*** Cannot allocate ion field array: mtest=',mtest
       call MPI_ABORT(MPI_COMM_WORLD,1,ierror)
    endif

    mimax=mi+1000*ceiling(sqrt(real(mi))) !ions array upper bound
    allocate(zion(1:nparami,mimax),zion0(nparami,mimax),wzion(mimax),wpion(ngyroi,mimax),&
         jtion0(ngyroi,mimax),jtion1(ngyroi,mimax),wtion0(ngyroi,mimax),&
         wtion1(ngyroi,mimax),jpion(ngyroi,mimax),rhogyroion(ngyroi,mimax),rhogyroion0(ngyroi,mimax),&
         STAT=mtest)

    if(mtest /= 0) then
       write(0,*)mype,'*** Cannot allocate ion: mtest=',mtest
       call MPI_ABORT(MPI_COMM_WORLD,1,ierror)
    endif

    pmarki=0.0_lk
    markeri=0.0_lk
    densityi=0.0_lk
    flowi=0.0_lk
    pressureipara=0.0_lk
    pressureiperp=0.0_lk
    zonali=0.0_lk
    zonalci=0.0_lk
    nonunii=0.0_lk
    marknumi=0
    rhogyroion=1.0_lk
    rhogyroion0=1.0_lk
    pfluxi=0.0_lk
    rdtemi=0.0_lk
    data1di=0.0_lk
#ifndef GPU_UM
    !$acc update device(rhogyroion,rhogyroion0)
#endif
    call gkload(mi, mimax, zion, zion0, wzion, wpion, jtion0, jtion1,&
        wtion0, wtion1, jpion,rhogyroion, ti_Eq, aion, qion, markeri, &
        pmarki, meshni, meshti1d, nonunii, marknumi, micore, misol, ngyroi)
    if(track_particles==1)then
        call loadi_tracking
    endif
    allocate(zion1(1:nparami+1,mimax))
    if(irun==0)then
#ifndef GPU_UM
      !$acc update device(zion,zion0)
#endif
    endif

    

end subroutine loadi
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine deprecated_particle_grid_array_alloc(nparam, ngyro, mp, mpmax, marknum, jtpart0, jtpart1, jppart,&
    wzpart, pmark, zonal, zonalc, marker, nonuni, zpart, zpart0, wppart, wtpart0, wtpart1, density,density0,&
    flow, pressurepara, pressureperp, rhogyropart, rhogyropart0, pfluxpart, rdtempart,data1dpart)
    use precision
    ! use particle_array
    ! use read_in_parameters
    use grids_setup, only: mpsi, mgrid
    use global_parameters, only: mype
    use particle_array, only: mpdata1d
    implicit none
    !dummy parameters
    integer :: nparam, ngyro, mp, mpmax
    integer, dimension(:),allocatable :: marknum
    integer,dimension(:,:),allocatable :: jtpart0,jtpart1,jppart
    real(lk),dimension(:),allocatable :: wzpart,pmark,zonal,zonalc,marker,nonuni
    real(lk),dimension(:,:),allocatable :: zpart,zpart0,wppart,wtpart0,wtpart1,&
        density,density0,flow,pressurepara,pressureperp,rhogyropart,rhogyropart0,&
        data1dpart
    real(lk),dimension(:),allocatable :: pfluxpart,rdtempart

    integer mtest,ierror
  ! allocate field memory
    allocate(pmark(mpsi),marker(mgrid),STAT=mtest)
    allocate(density(0:1,mgrid),density0(0:1,mgrid),flow(0:1,mgrid),STAT=mtest)
    allocate(pressurepara(0:1,mgrid),pressureperp(0:1,mgrid),STAT=mtest)
    allocate(zonal(mpsi),zonalc(mpsi),STAT=mtest)
    allocate(nonuni(mgrid),marknum(mgrid),STAT=mtest)
    mpmax=mp+1000*ceiling(sqrt(real(mp))) !ions array upper bound
    !mpmax=3
    !mp=1
    allocate(zpart(1:nparam,mpmax),zpart0(nparam,mpmax),wzpart(mpmax),wppart(ngyro,mpmax),&
         jtpart0(ngyro,mpmax),jtpart1(ngyro,mpmax),wtpart0(ngyro,mpmax),&
         wtpart1(ngyro,mpmax),jppart(ngyro,mpmax),rhogyropart(ngyro,mpmax),rhogyropart0(ngyro,mpmax),&
         pfluxpart(mpsi),rdtempart(mpsi),data1dpart(mpsi,mpdata1d),STAT=mtest)
  
    if(mtest /= 0) then
       write(0,*)mype,'*** Cannot allocate ion: mtest=',mtest
       call MPI_ABORT(MPI_COMM_WORLD,1,ierror)
    endif
    pmark=0.0_lk
    marker=0.0_lk
    density=0.0_lk
    flow=0.0_lk
    pressurepara=0.0_lk
    pressureperp=0.0_lk
    zonal=0.0_lk
    zonalc=0.0_lk
    nonuni=0.0_lk
    marknum=0
    rhogyropart0=1.0_lk
    pfluxpart=0.0_lk
    rdtempart=0.0_lk
    data1dpart=0.0_lk
end subroutine deprecated_particle_grid_array_alloc
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! restart from previous runs
subroutine restart_load
    use precision
    use global_parameters, only: mype, FileExit, irest
    use restart, only: restart_io
    implicit none

    integer iorestart,ierror
    character(len=9)cdum1

    if(mype==0)then
!!get restart directory info
        iorestart=345
        open(iorestart,file="FileExit.out",status="old")
        read(iorestart,"(A9,i1)")cdum1,FileExit
        read(iorestart,"(A9,i5)")cdum1,irest
        close(iorestart)
        irest=irest-1
    endif
    call mpi_bcast(irest,1,mpi_integer,0,mpi_comm_world,ierror)
    call MPI_BARRIER(MPI_COMM_WORLD,ierror)
    call restart_io("read")
    return
end subroutine restart_load
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine loadi_tracking
   use precision
   use particle_array
   use spline_function
   use equilibrium
   use global_parameters
   use read_in_parameters

   integer mtest, ierror
   real(lk) b, vth_norm, vthi, psidum, proton_charge, proton_mass, norm_energy, &
            energy1,energy2,magmom,energy
   
  ! allocate particle memory
    ! mimax=3
    ! mi=1
    ! allocate(zion(1:nparami,mimax),zion0(nparami,mimax),wzion(mimax),wpion(ngyroi,mimax),&
    !      jtion0(ngyroi,mimax),jtion1(ngyroi,mimax),wtion0(ngyroi,mimax),&
    !      wtion1(ngyroi,mimax),jpion(ngyroi,mimax),rhogyroion(ngyroi,mimax),STAT=mtest)
  
     !allocate(pos(3,mimax),STAT=mtest)

    !if(mtest /= 0) then
    !   write(0,*)mype,'*** Cannot allocate ion: mtest=',mtest
    !   call MPI_ABORT(MPI_COMM_WORLD,1,ierror)
    !endif

    zion(1,1)=2.0564_lk 
    zion(2,1)=zeta0+0.5_lk*deltaz
    zion(3,1)=0.0_lk
    
    b = spline2d(zion(1,1),zion(3,1)+zmax,0,lsr,lsz,spdR,spdZ,B0mag_Eq,0,0)
    
    vth_norm = ulength*0.01_lk/utime !m/s
    vthi = 9790.0_lk*sqrt(0.555_lk*etemp0)/sqrt(aion) !m/s
    
    ! uB=mv_{||}^2
    zion(4,1)=vthi/vth_norm
    zion(5,1)=1.0e-8
    zion(6,1)=sqrt(aion*vthi*vthi/(2.0_lk*b)/(vth_norm*vth_norm))
    
    
    psidum = spline2d(zion(1,1),zion(3,1)+zmax,0,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)
    energy=0.5_lk*aion*zion(4,1)*zion(4,1)+zion(6,1)*zion(6,1)*b
    
    proton_charge=1.602e-19 !C
    proton_mass=1.673e-27 !kg
    norm_energy=proton_mass*vth_norm*vth_norm/proton_charge   ! eV
    
    
    if(mype==0)write(gtcout,*)"zion1","R=",zion(1,1),"\zeta=",zion(2,1),"Z=",zion(3,1),"zion4",zion(4,1),"zion6",zion(6,1),"psidum=",psidum,"energy=",energy
    if(mype==0)write(gtcout,*)"vthi=",vthi,"vth_norm=",vth_norm,"b=",b,"norm_energy=",norm_energy 

end subroutine loadi_tracking
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

