subroutine radial_bin(species_name)
  use system_env
  use global_parameters, only: nparami,nparame, nparamf, nparamfe
  use particle_array
  use interfaces, only: binParticles
  implicit none 

  character(*), intent(in) :: species_name

  if (species_name=='thermal-ion') then
     call binParticles(zion, zion0, qion,aion,ngyroi,mi, mimax, nparami)
  else if (species_name=='thermal-electron') then
     call binParticles(zelectron, zelectron0, qelectron, aelectron,ngyroe,me, memax, nparame)
  else if (species_name=='fast-ion') then
     call binParticles(zfast, zfast0, qfast, afast, ngyrof,mf, mfmax, nparamf)
  else if (species_name=='fast-electron') then
     call binParticles(zfaste, zfaste0, qfaste, afaste, ngyrofe, mfe, mfemax, nparamfe)
  else
     write(*,*) 'radial_bin.F90: wrong choice'
     call mpi_barrier(mpi_comm_world, ierror)
  endif
  
end subroutine radial_bin

subroutine binParticles(zpart,zpart0,qpart,apart,ngyro,mp,mpmax,nparam)
  use system_env
  use precision
  use global_parameters, only: psi0,psi1,pi2,fielddir,rho0,rg0,mpsi,nthreads,mype,mgrid
  use field_array,only: deltar,deltat,deltap,deltaz,psimesh,zeta0,mtheta,&
    igrid,pgyro,pgyro2,tgyro,tgyro2,qmesh,qtinv,solvermethod
  use equilibrium,only: lsp,spdpsi_inv,spdpsi,rgpsi
  implicit none

  !declaration of the dummy arguments 
  real(lk), dimension(:,:) :: zpart,zpart0
  real(lk) qpart,apart
  integer ngyro
  integer mp, mpmax, nparam

  !declaration of local varaibles 
  integer i,j,tid,full_bin_count,new_pos,count
  integer m,igyro,ip,jt,ij0,ii,im,j00,j01,isp
  real(lk) pdum0,tdum0,zdum,rho,rg,tflr,pdum,tdum,dpx,delr,delt(0:mpsi),delp(mpsi),delz,rhoc,xdum,ydum,sqrtq0_inv,paxis
  !$acc declare create(delt(0:mpsi),delp(mpsi))
#ifndef _OPENACC
  integer, dimension(:,:), allocatable :: psi_count, psi_offsets
#else
  integer, dimension(:), allocatable :: psi_count, psi_offsets
#endif
  integer ipval(1:mpmax)
  !$acc declare create(psi_count,psi_offsets,ipval)

  delr = 1.0/deltar

#ifndef _OPENACC
  !$omp parallel
  !$omp master 
  nthreads = omp_get_num_threads()
  !$omp end master 
  !$omp end parallel
  
  allocate(psi_count(0:mgrid,0:nthreads-1), psi_offsets(0:mgrid,0:nthreads-1))

  !$omp parallel private(tid,m,pdum,isp,dpx,rg,ip,new_pos)

  tid = omp_get_thread_num()
  psi_count(0:mgrid,tid) = 0
  psi_offsets(0:mgrid,tid) = 0

  !$omp barrier 
#else
  allocate(psi_count(0:mgrid), psi_offsets(0:mgrid))

  !$acc parallel loop gang vector 
  do i=0, mgrid
     psi_count(i) = 0
  enddo
  !$acc end parallel

  do i=0,mgrid
     psi_offsets(i) = 0
  enddo

#endif

  delr=1.0/deltar
  delt=1.0/deltat
  delz=1.0/deltaz
  delp=1.0/deltap
  !$acc update device(delt(0:mpsi),delp(mpsi))
  rhoc=sqrt(apart)/(rho0*qpart) !apart & qpart are needed to calculate ion gyroradius
  if(ngyro==1)rhoc=0.0

  sqrtq0_inv=1.0_lk/sqrt(qmesh(0))
  paxis=0.0
  if(psi0<1e-6_lk)paxis=12.5*rho0*rho0*sqrtq0_inv!Use x-y locate for r < 5 rho0

#ifndef _OPENACC
  !$omp do 
#else
  !$acc parallel loop gang vector 
#endif
  do m=1, mp
     pdum0=zpart(1,m)
     tdum0=zpart(2,m)
     zdum=zpart(3,m)
     rho=zpart(6,m)*rhoc

! Guiding Center location
     isp=max(1,min(lsp-1,ceiling(pdum0*spdpsi_inv)))
     dpx=pdum0-spdpsi*real(isp-1)
! radaial spline of rg avoids sigularity near axis: y=y1+y2*sqrt(x)+y3*x for psi<spdpsi
     if(isp==1)dpx=sqrt(dpx)
     rg=  rgpsi(1,isp)  +rgpsi(2,isp)*dpx  +rgpsi(3,isp)*dpx*dpx
     ip=max(0,min(mpsi,floor((rg-rg0)*delr+0.5)))       !radial grid closest to GC
     jt=max(0,min(mtheta(ip),floor(tdum0*delt(ip)+0.5))) !poloidal grid closest to GC
     ij0=igrid(ip)+jt                                   !GC grid index in magnetic coordinates

     igyro=1
     pdum=pdum0+rho*pgyro(igyro,ij0)+rho*rho*pgyro2(igyro,ij0) !particle position
     !if(pdum<0)print *,"pdum,bounded pdum:",pdum,modulo(pdum-psi0,psi1-psi0)+psi0
     !pdum=modulo(pdum-psi0,psi1-psi0)+psi0 !periodic BC
     if(pdum < psi0)pdum=2.0*psi0-pdum !reflective BC
     if(pdum > psi1)pdum=2.0*psi1-pdum
     ! particle position in theta
     tflr=tdum0+rho*tgyro(igyro,ij0)+rho*rho*tgyro2(igyro,ij0)
     ! Near axis gyroaveraging
     if(pdum0<paxis)then
        pdum = sqrt(pdum0)
        xdum = pdum*cos(tdum0)+rho*rho0*sqrtq0_inv*cos(tdum0+pi2*real(igyro-1)/4.0) 
        ydum = pdum*sin(tdum0)+rho*rho0*sqrtq0_inv*sin(tdum0+pi2*real(igyro-1)/4.0) 
        pdum = max(1.0e-8_lk*psi1,xdum*xdum+ydum*ydum)
        tflr = sign(1.0_lk,ydum)*acos(max(-1.0_lk,min(1.0_lk,xdum/&
             sqrt(zpart(1,m)))))
     endif
     isp=max(1,min(lsp-1,ceiling(pdum*spdpsi_inv)))     
     dpx=pdum-spdpsi*real(isp-1)
     rg=rgpsi(1,isp)+rgpsi(2,isp)*dpx+rgpsi(3,isp)*dpx*dpx
     ii=max(0,min(mpsi-1,floor((rg-rg0)*delr)))
     im=ii
     if(fielddir==1 .or. fielddir==3)then
        tdum=modulo(tflr-(zdum-pi2)*qtinv(im),pi2)*delt(im)
     else
        tdum=modulo(tflr-zdum*qtinv(im),pi2)*delt(im)
     endif
     j00=max(0,min(mtheta(im)-1,floor(tdum)))
     ip=igrid(im)+j00
#ifndef _OPENACC
     psi_count(ip,tid) = psi_count(ip,tid)+1
#else
     !$acc atomic update
     psi_count(ip) = psi_count(ip)+1
     !$acc end atomic
#endif
     ipval(m) = ip
  enddo
#ifndef _OPENACC
  !$omp end do 
#else
  !$acc end parallel
#endif

#ifndef _OPENACC
  !$omp barrier 

  if (tid.eq.0) then
     do i=0, mgrid-1
        full_bin_count = 0
        do j=0, nthreads-1
           full_bin_count = full_bin_count + psi_count(i,j)
        enddo
        psi_offsets(i+1,0) = psi_offsets(i,0) + full_bin_count
     enddo 
     
     do i=0, mgrid
        do j=1, nthreads-1
           psi_offsets(i,j) = psi_offsets(i,j-1) + psi_count(i,j-1)
        enddo 
     enddo 
          
     if (psi_offsets(mgrid,nthreads-1)+psi_count(mgrid,nthreads-1).ne.mp) then
        write(*,*) 'mype=', mype, 'psi_offsets(mgrid,nthreads-1)=', psi_offsets(mgrid,nthreads-1), &
             'psi_count(mgrid,nthreads-1)=', psi_count(mgrid,nthreads-1), 'mp=', mp
        write(*,*) 'error in psi_offsets: radial_bin.F90'
        call mpi_barrier(mpi_comm_world, ierror)
     endif
  endif

  !$omp barrier

  psi_count(0:mgrid,tid) = 0

  !$omp barrier 
#else

  !$acc update host(psi_count)
  do i=0, mgrid-1
     psi_offsets(i+1) = psi_offsets(i) + psi_count(i)
  enddo
  if (psi_offsets(mgrid)+psi_count(mgrid).ne.mp) then
     write(*,*) 'mype=', mype, 'psi_offsets(mgrid)=', psi_offsets(mgrid), &
          'psi_count(mgrid)=', psi_count(mgrid), 'mp=', mp
     write(*,*) 'error in psi_offsets: radial_bin.F90'
     stop
  endif
  !$acc update device(psi_offsets)

  !$acc parallel loop gang vector 
  do i=0, mgrid
     psi_count(i)=0
  enddo
  !$acc end parallel

#endif

#ifndef _OPENACC
  !$omp do 
#else
  !$acc parallel loop gang vector 
#endif
  do m=1, mp
     ip = ipval(m)
#ifndef _OPENACC
     psi_count(ip,tid) = psi_count(ip,tid) + 1
     new_pos = psi_offsets(ip,tid) + psi_count(ip,tid) 
#else
     !$acc atomic capture 
     psi_count(ip) = psi_count(ip) + 1
     count = psi_count(ip)
     !$acc end atomic 
     new_pos = psi_offsets(ip) + count
#endif
     zpart0(1:nparam,new_pos) = zpart(1:nparam,m) 
  enddo
#ifndef _OPENACC
  !$omp end do
#else
  !$acc end parallel 
#endif

#ifndef _OPENACC
  !$omp do
#else
  !$acc parallel loop gang vector 
#endif
  do m=1, mp
     zpart(1:nparam,m) = zpart0(1:nparam,m)
  enddo
#ifndef _OPENACC
  !$omp end do
  
  !$omp end parallel

#else
  !$acc end parallel
#endif

  deallocate(psi_count,psi_offsets)
   
end subroutine binParticles
