subroutine charge(species_name)
    use read_in_parameters
    use global_parameters
    use particle_array
    use grids_eqdata
    implicit none
    character(*),intent(in) :: species_name
interface
    subroutine charge_kernel(qpart,apart,density,density0,flow,zpart,wppart,wtpart0,wtpart1,jtpart0,jtpart1,wzpart,&
        zonal,zonalc,marker,pmark,mp,pressurepara,pressureperp,rhogyropart,rhogyropart0, ngyro, fullk)
        use precision
        implicit none
        !dummy parameter
        real(lk) :: qpart, apart
        real(lk), dimension(:) :: zonal,zonalc,pmark,wzpart,marker
        integer,dimension(:,:) :: jtpart0,jtpart1
        integer mp, ngyro
        logical fullk
        real(lk),dimension(:,:) :: wppart,wtpart0,wtpart1,zpart,rhogyropart,rhogyropart0
        real(lk),dimension(0:,:) :: density,density0,flow,pressurepara,pressureperp
    end subroutine charge_kernel
end interface
    if(species_name == 'gki-thermal')then
        call charge_kernel(qion, aion, densityi,densityi0, flowi, zion, wpion, wtion0, wtion1, jtion0, jtion1, wzion,&
            zonali, zonalci, markeri, pmarki, mi, pressureipara, pressureiperp,rhogyroion,rhogyroion0, ngyroi, .false.)
    elseif(species_name == 'electron')then
        call charge_kernel(qelectron, aelectron, densitye,densitye0, flowe, zelectron, wpelectron, wtelectron0, wtelectron1,&
        jtelectron0, jtelectron1, wzelectron, zonale, zonalce, markere, pmarke, me, pressureepara, pressureeperp,&
        rhogyroelectron,rhogyroelectron0,ngyroe, .false.)
    elseif(species_name == 'fki-thermal')then
        call charge_kernel(qfki, afki, densityfki,densityfki0, flowfki, zfki, wpfki, wtfki0,wtfki1, jtfki0, jtfki1,wzfki,&
            zonalfki, zonalcfki, markerfki, pmarkfki, mfki, pressurefkipara, pressurefkiperp, rhogyrofki,rhogyrofki0, 1, .true.)
    endif
end subroutine charge
subroutine locate(species_name)
    use global_parameters
    use read_in_parameters
    use global_parameters
    use particle_array
    implicit none
    character(*),intent(in) :: species_name
interface
    subroutine locate_kernel(zpart,wzpart,wppart,wtpart0,wtpart1,jtpart0,&
        jtpart1,ngyro,mp,jppart,qpart,apart,rhogyropart,fullk)
        use precision
        implicit none
    ! dummy parameters
        real(lk), dimension(:,:) :: zpart, wppart, wtpart0, wtpart1,rhogyropart
        integer, dimension(:,:) :: jtpart0, jtpart1, jppart
        real(lk), dimension(:) :: wzpart
        integer :: mp,ngyro
        logical :: fullk
        real(lk) :: qpart,apart
    end subroutine locate_kernel
end interface
    if(species_name=='gki-thermal')then
        call locate_kernel(zion,wzion,wpion,wtion0,wtion1,jtion0,jtion1,ngyroi,mi,jpion,qion,aion,rhogyroion,.false.)
    elseif(species_name=='electron')then
        call locate_kernel(zelectron,wzelectron,wpelectron,wtelectron0,wtelectron1,jtelectron0,&
        jtelectron1,ngyroe,me,jpelectron,qelectron,aelectron,rhogyroelectron,.false.)
    elseif(species_name=='fki-thermal')then
        call locate_kernel(zfki,wzfki,wpfki,wtfki0,wtfki1,jtfki0,jtfki1,1,mfki,jpfki,qfki,afki,rhogyrofki,.true.)
    endif
end subroutine locate

subroutine charge_kernel(qpart,apart,density,density0,flow,zpart,wppart,wtpart0,wtpart1,jtpart0,jtpart1,wzpart,&
    zonal,zonalc,marker,pmark,mp,pressurepara,pressureperp,rhogyropart,rhogyropart0,ngyro, fullk)
    use global_parameters
    use read_in_parameters, only: min_r, frcfullf, n_modes, sol_para_periodicity, ifluxtube
    use equilibrium, only: spdR, spdZ, lsr, lsz, zmax, psirz_Eq
    use grids_setup, only: mgrid, mpsi, mpsi_sol, mpsi_core, mpara_core, mpara_sol, igrid_core, igrid_sol, &
        mgrid_core, core_para_periodicity, sol_para_period, particle_field_accumulation_sep
    use grids_fields, only: smooth, toroidal_filter
    use simulation_domain, only: sfcore, sfsol, couplefrc
    use precision
    implicit none

    !dummy parameter
    real(lk) :: qpart, apart
    real(lk), dimension(:) :: zonal,zonalc,pmark,wzpart,marker
    integer,dimension(:,:) :: jtpart0,jtpart1
    integer mp, ngyro
    logical fullk
    real(lk),dimension(:,:) :: wppart,wtpart0,wtpart1,zpart,rhogyropart,rhogyropart0
    real(lk),dimension(0:,:) :: density,density0,flow,pressurepara,pressureperp

    !local parameter
    integer :: m,i,im,j,k,ip,jt,ii,j11,j10,j01,j00,ierror,igyro,ij,ipjt,isp,jst,ia,&
        icount,idest,isource,isendtag,irecvtag,istatus(MPI_STATUS_SIZE),ij0,ij1,ijn,&
        nrsp,nzsp
    
    real(lk) :: gyrodum,weight,rdum,rdum_inv,zdum,upara,energypara,energyperp,cmratio,wz1,wz0,wp1,wp0,wt11,wt10,wt01,wt00,&
         sendl(mgrid,5),recvr(mgrid,5),dntmp(0:1,mgrid),dntmp0(0:1,mgrid),djtmp(0:1,mgrid),dpparatmp(0:1,mgrid),dpperptmp(0:1,mgrid),&
         zf0,zc0,temp,adum(mpsi),adum2(0:1,mgrid),&
         dpr,dr2,dpz,dz2,gradpsi_r,gradpsi_z,br,bz,b,spdR_inv,spdZ_inv,wz1_00,wz0_00
    real(lk) total_density,avg_density,igyro_coef


    cmratio=qpart/apart
    spdR_inv=1.0_lk/spdR
    spdZ_inv=1.0_lk/spdZ
!$acc parallel
    do ij=1,mgrid  
    density(:,ij)=0.0_lk
    density0(:,ij)=0.0_lk
    flow(:,ij)=0.0_lk
    pressurepara(:,ij)=0.0_lk
    pressureperp(:,ij)=0.0_lk
    enddo
!$acc end parallel

! scatter ion density for electrostatic simulation
#ifdef _OPENACC
!$acc parallel loop gang vector
#else
!$omp do private(m,rdum,rdum_inv,zdum,nrsp,dpr,dr2,nzsp,dpz,dz2,&
!$omp& igyro_coef,gradpsi_r,gradpsi_z,br,bz,b,&
!$omp& upara,energypara,energyperp,igyro,weight,wz1,wz0,&
!$omp& wz0_00,wz1_00,wp1,wp0,wt10,wt00,wt11,wt01,ij) reduction(+:density,density0,flow,pressurepara,pressureperp)
#endif
    do m=1,mp
        rdum = max(zpart(1,m),min_r)
        rdum_inv=1.0_lk/rdum
        zdum = zpart(3,m)

        nrsp=max(1,min(lsr-1,ceiling(rdum*spdR_inv)))
        dpr=rdum-spdR*real(nrsp-1)
        dr2=dpr*dpr

        nzsp=max(1,min(lsz-1,ceiling((zdum+zmax)*spdZ_inv)))
        dpz=zdum+zmax-spdZ*real(nzsp-1)
        dz2=dpz*dpz

        gradpsi_r= psirz_Eq(2,nrsp,nzsp)+ 2.0*psirz_Eq(3,nrsp,nzsp)*dpr &
             + psirz_Eq(5,nrsp,nzsp)*dpz + 2.0*psirz_Eq(6,nrsp,nzsp)*dpr*dpz &
             + psirz_Eq(8,nrsp,nzsp)*dz2 + 2.0*psirz_Eq(9,nrsp,nzsp)*dpr*dz2

        gradpsi_z=psirz_Eq(4,nrsp,nzsp) + psirz_Eq(5,nrsp,nzsp)*dpr + psirz_Eq(6,nrsp,nzsp)*dr2 &
                 +2.0*psirz_Eq(7,nrsp,nzsp)*dpz + 2.0*psirz_Eq(8,nrsp,nzsp)*dpr*dpz + 2.0*psirz_Eq(9,nrsp,nzsp)*dr2*dpz

        br=-gradpsi_z*rdum_inv
        bz=gradpsi_r*rdum_inv
        b=sqrt(br**2+bz**2)

        upara=zpart(4,m)
        energypara=0.5_lk*apart*upara*upara
        energyperp=zpart(6,m)*zpart(6,m)*b
        weight=zpart(5,m)!*zion0(6,m)

        do igyro=1,ngyro
            wp1=wppart(igyro,m)       !outer flux surface
            wp0=1.0_lk-wp1              !inner flux surface
        
            wt10=wp0*wtpart0(igyro,m) !upper Z grid on inner flux surface
            wt00=wp0-wt10            !lower Z grid on inner flux surface
        
            wt11=wp1*wtpart1(igyro,m) !upper Z grid on outer flux surface
            wt01=wp1-wt11            !lower Z grid on outer flux surface
            
            igyro_coef = rhogyropart(igyro,m)
            wz1=igyro_coef*weight*wzpart(m)         !weight for upper toroidal grid
            wz0=igyro_coef*weight-wz1              !weight for lower toroidal grid
            igyro_coef = rhogyropart0(igyro,m)
            wz1_00=igyro_coef*weight*wzpart(m)         !weight for uppertoroidal grid
            wz0_00=igyro_coef*weight-wz1              !weight for lower toroidal grid

! If no loop-level parallelism, write directly into array "densityi()".
            ij=jtpart0(igyro,m)                         !lower poloidal grid on inner flux surface
            !$acc atomic update
            density(0,ij) = density(0,ij) + wz0*wt00 !lower toroidal grid
            !$acc atomic update
            density(1,ij) = density(1,ij) + wz1*wt00 !upper toroidal grid
            !$acc atomic update
            flow(0,ij) = flow(0,ij) + wz0*wt00*upara !lower toroidal grid
            !$acc atomic update
            flow(1,ij) = flow(1,ij) + wz1*wt00*upara !upper toroidal grid
            !$acc atomic update
            pressurepara(0,ij) = pressurepara(0,ij) + wz0*wt00*energypara/igyro_coef
            !$acc atomic update
            pressurepara(1,ij) = pressurepara(1,ij) + wz1*wt00*energypara/igyro_coef
            !$acc atomic update
            pressureperp(0,ij) = pressureperp(0,ij) + wz0*wt00*energyperp/igyro_coef
            !$acc atomic update
            pressureperp(1,ij) = pressureperp(1,ij) + wz1*wt00*energyperp/igyro_coef
            !$acc atomic update
            density0(0,ij) = density0(0,ij) + wz0_00*wt00
            !$acc atomic update
            density0(1,ij) = density0(1,ij) + wz1_00*wt00

            ij=ij+1                                    !upper poloidal grid on inner flux surface
            !$acc atomic update
            density(0,ij) = density(0,ij) + wz0*wt10
            !$acc atomic update
            density(1,ij) = density(1,ij) + wz1*wt10
            !$acc atomic update
            flow(0,ij) = flow(0,ij) + wz0*wt10*upara
            !$acc atomic update
            flow(1,ij) = flow(1,ij) + wz1*wt10*upara
            !$acc atomic update
            pressurepara(0,ij) = pressurepara(0,ij) + wz0*wt10*energypara/igyro_coef
            !$acc atomic update
            pressurepara(1,ij) = pressurepara(1,ij) + wz1*wt10*energypara/igyro_coef
            !$acc atomic update
            pressureperp(0,ij) = pressureperp(0,ij) + wz0*wt10*energyperp/igyro_coef
            !$acc atomic update
            pressureperp(1,ij) = pressureperp(1,ij) + wz1*wt10*energyperp/igyro_coef
            !$acc atomic update
            density0(0,ij) = density0(0,ij) + wz0_00*wt10
            !$acc atomic update
            density0(1,ij) = density0(1,ij) + wz1_00*wt10

            ij=jtpart1(igyro,m)                         !lower poloidal grid on outer flux surface
            !$acc atomic update
            density(0,ij) = density(0,ij) + wz0*wt01
            !$acc atomic update
            density(1,ij) = density(1,ij) + wz1*wt01
            !$acc atomic update
            flow(0,ij) = flow(0,ij) + wz0*wt01*upara
            !$acc atomic update
            flow(1,ij) = flow(1,ij) + wz1*wt01*upara
            !$acc atomic update
            pressurepara(0,ij) = pressurepara(0,ij) + wz0*wt01*energypara/igyro_coef
            !$acc atomic update
            pressurepara(1,ij) = pressurepara(1,ij) + wz1*wt01*energypara/igyro_coef
            !$acc atomic update
            pressureperp(0,ij) = pressureperp(0,ij) + wz0*wt01*energyperp/igyro_coef
            !$acc atomic update
            pressureperp(1,ij) = pressureperp(1,ij) + wz1*wt01*energyperp/igyro_coef
            !$acc atomic update
            density0(0,ij) = density0(0,ij) + wz0_00*wt01
            !$acc atomic update
            density0(1,ij) = density0(1,ij) + wz1_00*wt01

            ij=ij+1                                    !upper poloidal grid on outer flux surface
            !$acc atomic update
            density(0,ij) = density(0,ij) + wz0*wt11
            !$acc atomic update
            density(1,ij) = density(1,ij) + wz1*wt11
            !$acc atomic update
            flow(0,ij) = flow(0,ij) + wz0*wt11*upara
            !$acc atomic update
            flow(1,ij) = flow(1,ij) + wz1*wt11*upara
            !$acc atomic update
            pressurepara(0,ij) = pressurepara(0,ij) + wz0*wt11*energypara/igyro_coef
            !$acc atomic update
            pressurepara(1,ij) = pressurepara(1,ij) + wz1*wt11*energypara/igyro_coef
            !$acc atomic update
            pressureperp(0,ij) = pressureperp(0,ij) + wz0*wt11*energyperp/igyro_coef
            !$acc atomic update
            pressureperp(1,ij) = pressureperp(1,ij) + wz1*wt11*energyperp/igyro_coef
            !$acc atomic update
            density0(0,ij) = density0(0,ij) + wz0_00*wt11
            !$acc atomic update
            density0(1,ij) = density0(1,ij) + wz1_00*wt11

        enddo
    enddo
    !$acc end parallel
#ifndef GPU_UM
    !$acc update host(density,density0,flow,pressurepara,pressureperp)
#endif

! If we have a particle decomposition on the toroidal domains, do a reduce
! operation to add up all the contributions to charge density on the grid
    if(npartdom>1)then
!$omp parallel do private(ij)
        do ij=1,mgrid
            dntmp(0:1,ij)=density(0:1,ij)
            density(0:1,ij)=0.0_lk
            dntmp0(0:1,ij)=density0(0:1,ij)
            density0(0:1,ij)=0.0_lk
            djtmp(0:1,ij)=flow(0:1,ij)
            flow(0:1,ij)=0.0_lk
            dpparatmp(0:1,ij)=pressurepara(0:1,ij)
            pressurepara(0:1,ij)=0.0_lk
            dpperptmp(0:1,ij)=pressureperp(0:1,ij)
            pressureperp(0:1,ij)=0.0_lk
        enddo
        call MPI_ALLREDUCE(dntmp,density,(mgrid*2),mpi_Rsize,MPI_SUM,partd_comm,ierror)
        call MPI_ALLREDUCE(dntmp0,density0,(mgrid*2),mpi_Rsize,MPI_SUM,partd_comm,ierror)
        call MPI_ALLREDUCE(djtmp,flow,(mgrid*2),mpi_Rsize,MPI_SUM,partd_comm,ierror)
        call MPI_ALLREDUCE(dpparatmp,pressurepara,(mgrid*2),mpi_Rsize,MPI_SUM,partd_comm,ierror)
        call MPI_ALLREDUCE(dpperptmp,pressureperp,(mgrid*2),mpi_Rsize,MPI_SUM,partd_comm,ierror)
    endif

! toroidal end cell
    sendl(:,1)=density(0,:)
    sendl(:,2)=flow(0,:)
    sendl(:,3)=pressurepara(0,:)
    sendl(:,4)=pressureperp(0,:)
    sendl(:,5)=density0(0,:)
    recvr=0.0_lk
    icount=5*mgrid
    idest=left_pe
    isource=right_pe
    isendtag=myrank_toroidal
    irecvtag=isource

! send densityi to left and receive from right
    call MPI_SENDRECV(sendl,icount,mpi_Rsize,idest,isendtag,&
         recvr,icount,mpi_Rsize,isource,irecvtag,toroidal_comm,istatus,ierror)
    
    ! B.C. at zeta<2*pi is continuous
    density(1,:)=density(1,:)+recvr(:,1)
    flow(1,:)=flow(1,:)+recvr(:,2)
    pressurepara(1,:)=pressurepara(1,:)+recvr(:,3)
    pressureperp(1,:)=pressureperp(1,:)+recvr(:,4)
    density0(1,:)=density0(1,:)+recvr(:,5)

    gyrodum=1.0_lk/real(ngyro)
    zonal=0.0_lk
    zonalc=0.0_lk

! gyro points normalization
    if(sfcore==1)then
!$omp parallel do private(i,j,ij)
        do i=1,mpsi_core
            do j=1,mpara_core(i)
                ij=igrid_core(i)+j
                density(1,ij)=gyrodum*density(1,ij)
                density0(1,ij)=gyrodum*density0(1,ij)
                zonal(i)=zonal(i)+density0(1,ij)
                flow(1,ij)=gyrodum*flow(1,ij)
                zonalc(i)=zonalc(i)+flow(1,ij)
                pressurepara(1,ij)=gyrodum*pressurepara(1,ij)
                pressureperp(1,ij)=gyrodum*pressureperp(1,ij)
            enddo
        enddo
    endif
    if(sfsol==1)then
!$omp parallel do private(i,j,ij)
        do i=1,mpsi_sol
            do j=1,mpara_sol(i)
                ij=mgrid_core+igrid_sol(i)+j
                density(1,ij)=gyrodum*density(1,ij)
                density0(1,ij)=gyrodum*density0(1,ij)
                zonal(i)=zonal(i)+density0(1,ij)
                flow(1,ij)=gyrodum*flow(1,ij)
                zonalc(i)=zonalc(i)+flow(1,ij)
                pressurepara(1,ij)=gyrodum*pressurepara(1,ij)
                pressureperp(1,ij)=gyrodum*pressureperp(1,ij)
           enddo
        enddo
    endif
! periodic boundary for core region: sum and make equal of particle moments at \theta=0/\theta=2*\pi
    if(sfcore==1)then
        call core_para_periodicity(density)
        call core_para_periodicity(density0)
        call core_para_periodicity(flow)
        call core_para_periodicity(pressurepara)
        call core_para_periodicity(pressureperp)
    endif
! parallel boundary condition for SOL
    if(sfsol==1)then
! periodic
        if(sol_para_periodicity==1)then
            call sol_para_period(density)
            call sol_para_period(density0)
            call sol_para_period(flow)
            call sol_para_period(pressurepara)
            call sol_para_period(pressureperp)
        else
!$omp parallel do private(i,ij1,ijn)
            do i=1,mpsi_sol
                ij1 = mgrid_core+igrid_sol(i)+1
                ijn = mgrid_core+igrid_sol(i)+mpara_sol(i)
                density(0:1,ij1)=0.0_lk
                flow(0:1,ij1)=0.0_lk
                pressurepara(0:1,ij1)=0.0_lk
                pressureperp(0:1,ij1)=0.0_lk
                density(0:1,ijn)=0.0_lk
                flow(0:1,ijn)=0.0_lk
                pressurepara(0:1,ijn)=0.0_lk
                pressureperp(0:1,ijn)=0.0_lk
            enddo
        endif
    endif
! accumulate particle moments on the shared grids at separatrix
    if(couplefrc==1)then
        call particle_field_accumulation_sep(density)
        call particle_field_accumulation_sep(density0)
        call particle_field_accumulation_sep(flow)
        call particle_field_accumulation_sep(pressurepara)
        call particle_field_accumulation_sep(pressureperp)
    endif

!marker normalization
    if(sfcore==1)then
!$omp parallel do private(i,j,ij)
        do i=1,mpsi_core
            do j=1,mpara_core(i)
                ij=igrid_core(i)+j
                density(1,ij)=density(1,ij)*marker(ij)
                density0(1,ij)=density0(1,ij)*marker(ij)
                flow(1,ij)=flow(1,ij)*marker(ij)
                pressurepara(1,ij)=pressurepara(1,ij)*marker(ij)
                pressureperp(1,ij)=pressureperp(1,ij)*marker(ij)
            enddo
        enddo
    endif
    if(sfsol==1)then
!$omp parallel do private(i,j,ij)
        do i=1,mpsi_sol
            do j=1,mpara_sol(i)
                ij=mgrid_core+igrid_sol(i)+j
                density(1,ij)=density(1,ij)*marker(ij)
                density0(1,ij)=density0(1,ij)*marker(ij)
                flow(1,ij)=flow(1,ij)*marker(ij)
                pressurepara(1,ij)=pressurepara(1,ij)*marker(ij)
                pressureperp(1,ij)=pressureperp(1,ij)*marker(ij)
            enddo
        enddo 
    endif

! calculate the zonal modes    
    adum=0.0_lk 
    call MPI_ALLREDUCE(zonal,adum,mpsi,mpi_Rsize,MPI_SUM,toroidal_comm,ierror)
    zonal=adum/pmark
      
    adum=0.0_lk
    call MPI_ALLREDUCE(zonalc,adum,mpsi,mpi_Rsize,MPI_SUM,toroidal_comm,ierror)
    zonalc=adum/pmark
    
    if(sfsol==1)then
!$omp parallel do private(i,j,ij)
        do i=1,mpsi_sol
             do j=1,mpara_sol(i)
                ij=mgrid_core+igrid_sol(i)+j
            density(1,ij)=density(1,ij)-zonal(i)
            flow(1,ij)=flow(1,ij)-zonalc(i)
             enddo
        enddo
    endif
    
    if(frcfullf==1)then
!$omp parallel do private(i)
        do i=1,mgrid
            density(1,i)=density(1,i)-1.0_lk
            flow(1,i)=flow(1,i)-0.0_lk ! U0para=0.0
        enddo
    endif
    
    total_density=0.0_lk
    do i=1,mgrid
        total_density = total_density + density(1,i)
    enddo
    
    avg_density=total_density/mgrid

    CALL SMOOTH(density)
    CALL SMOOTH(density0)
    CALL SMOOTH(flow)
    CALL TOROIDAL_FILTER(density)
    CALL TOROIDAL_FILTER(flow)

    adum2=0.0_lk
    call MPI_ALLREDUCE(density0,adum2,2*mgrid,mpi_Rsize,MPI_SUM,toroidal_comm,ierror)
    density0=adum2/mtoroidal

!Note: density0 is the zonal density
    do i=1,mpsi_core
       do j=1,mpara_core(i)
          ij=igrid_core(i)+j
          density0(1,ij)=zonal(i)
       enddo
    enddo
    do i=1,mpsi_sol
       do j=1,mpara_sol(i)
          ij=mgrid_core+igrid_sol(i)+j
          density0(1,ij)=zonal(mpsi_core+i)
       enddo
    enddo

    if(ifluxtube==0)then
    if(sfsol==1)then
        do j=1,mpara_sol(mpsi_sol)
           density(0:1,mgrid_core+igrid_sol(mpsi_sol)+j)=0.0_lk
        enddo
        if(couplefrc==0)then
            do j=1,mpara_sol(1)
                density(0:1,mgrid_core+igrid_sol(1)+j)=0.0_lk
            enddo
        endif
    endif
    if(sfcore==1)then
        do j=1,mpara_core(1)
            density(0:1,igrid_core(1)+j)=0.0_lk
        enddo
        if(couplefrc==0)then
            do j=1,mpara_core(mpsi_core)
                density(0:1,igrid_core(mpsi_core)+j)=0.0_lk
            enddo
        endif
    endif
    else
    if(sfsol==1)then
    dntmp=0.0_lk
    ! average density in radial direction
    temp=1.0_lk/mpsi_sol
    do i=1,mpsi_sol
       do j=1,mpara_sol(i)
          ij = igrid_sol(i)+j
          dntmp(0:1,igrid_sol(1)+j) = dntmp(0:1,igrid_sol(1)+j) + density(0:1,ij)*temp
       enddo
    enddo
    do i=1,mpsi_sol
       do j=1,mpara_sol(i)
          ij = igrid_sol(i)+j
          density(0:1,ij) = dntmp(0:1,igrid_sol(1)+j)
       enddo
    enddo
    endif
    endif

end subroutine charge_kernel
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! find ion location for interpolation in both gathering (chargei.F90) and scattering (pushi.F90)
subroutine locate_kernel(zpart,wzpart,wppart,wtpart0,wtpart1,jtpart0,&
    jtpart1,ngyro,mp,jppart,qpart,apart,rhogyropart,fullk)
    use global_parameters
    use grids_setup
    use read_in_parameters, only: min_r, Zsol_domain0, n_modes
    use mapping_splines, only: sspsol_psi_z, sspcore_psi_theta, msp_psi_core, msp_psi_sol, &
                               lspsol,lszsol,del_psi_sol,del_z_sol, &
                               lspcore,lstcore,del_psi_core,del_theta_core
    use simulation_domain, only: psi0_sol, psi1_sol, psi0_core, psi1_core, &
                                 sfcore, sfsol, couplefrc
    use spline_function, only: spline1d,spline2d
    use cylindrical_function,only: geotheta
    use equilibrium, only: lsr, lsz, spdR, spdZ, psi_sep, zmax, &
                           psirz_Eq
    use gyro_array, only: pgyro, zetagyro, tgyro
    use precision
    implicit none
! dummy parameters
    real(lk), dimension(:,:) :: zpart, wppart, wtpart0, wtpart1, rhogyropart
    integer, dimension(:,:) :: jtpart0, jtpart1, jppart
    real(lk), dimension(:) :: wzpart
    integer :: mp,ngyro
    logical :: fullk
    real(lk) :: qpart,apart,dbesjn
! local parameters
  
    integer :: m,nrsp,nzsp,np,ij,i,j,iz0(1),iz1(1),jj0,jj1,ij0,ij1,is0(1),is1(1),ij0_rg,igyro,left,right,ierror, thread_idx, num_threads
    real(lk) :: rdum,zdum,zetadum,delz,dpr,dpz,dr2,dz2,psifrc,psifrc0,sign_psi,spdR_inv,spdZ_inv,tdum,&
             sdum,sdum0,sdum1,ns0,ns1,&
             dz0,dz1,ds0,ds1,rho_coef,rho_ratio,gc_coef,dpx,dx2,pdum,ds_tmp(mpara_sol(1))
    real(lk),dimension(:),allocatable :: dels,dels_inv
#ifndef GPU_UM
    !$acc declare create(dels, dels_inv,ds_tmp)
#endif


    !rho_coef=sqrt(apart)/(rho0*qpart) !apart & qpart are needed to calculate ion gyroradius
    
    if(fullk .eqv. .false.)then
        rho_coef=1.0_lk*sqrt(apart)/(rho0*qpart) !apart & qpart are needed to calculate ion gyroradius
        if(ngyro == 1)then
            gc_coef=0.0_lk
        else
            gc_coef=1.0_lk
        endif
    else
        gc_coef=0.0_lk
        rho_coef=0.0_lk
    endif

    delz=1.0_lk/deltaz
  
    spdR_inv=1.0_lk/spdR
    spdZ_inv=1.0_lk/spdZ

  
    if(sfcore==1)then
        allocate(dels(mpsi_core),dels_inv(mpsi_core))
        dels(:)=deltas_core1d(:)/totals_core(:)
        dels_inv(:)=1.0_lk/dels(:)
#ifndef GPU_UM
    !$acc update device(dels, dels_inv)
#endif
    endif

! if (fullk == .false.)then
!   gklocate
! else
!   fklocate
! endif
! TODO: implement gklocate based on gyro averaging
#ifdef _OPENACC
!$acc parallel loop gang vector private(is0, is1, ds_tmp)
#else
!$omp parallel do private(m,rdum,zdum,zetadum,nrsp,dpr,dr2,nzsp,dpz,dz2,psifrc,&
!$omp& sign_psi,np,tdum,sdum0,sdum1,ns0,ns1,&
!$omp& igyro,i,psifrc0,ij0_rg,rho_ratio,&
!$omp& jj0,jj1,ij0,ij1,is0,is1,ds0,ds1,pdum,dpx,dx2,j,ds_tmp)
#endif
    do m=1,mp
        rdum=max(zpart(1,m),min_r)
        zdum=zpart(3,m)
        zetadum=zpart(2,m)
        rho_ratio = rho_coef*zpart(6,m)
  
        wzpart(m)=(zetadum-zeta0)*delz  !weight for upper toroidal grid

! psi value at particle position by using spline function from equilibrium  
        psifrc0=zpart(7,m)    
    !!! index of the closest grid near Guiding Center location, i.e. ij0_rg
        ij0_rg = 0
        if(sfsol==1)then
            if((psifrc0 .ge. psi_sep) .or. sfcore==0)then     
                if(psifrc0 .ge. psi1_sol)then
                    np=mpsi_sol-1
                    psifrc=psi1_sol
                elseif(psifrc0 .le. psi0_sol .and. sfcore==0)then
                    np=1
                    psifrc=psi0_sol
                else
                    left = 1
                    right = mpsi_sol
                    do while(right-left>1)
                      i=floor((right-left)/2.0)+left
                      if(psifrc0>=psimesh_sol(i) .and. psifrc0<psimesh_sol(i+1))then
                        left=i
                        exit
                      elseif(psifrc0<psimesh_sol(i))then
                        right=i
                      else
                        left=i+1
                      endif
                    enddo
                    np = left
                    psifrc = psifrc0
                endif

                pdum=max(0.0_lk,psifrc-msp_psi_sol(1))
                zdum=max(0.0_lk,zpart(3,m)-Zsol_domain0)
                i = max(1, min(lspsol-1, ceiling(pdum/del_psi_sol)))
                j = max(1, min(lszsol-1, ceiling(zdum/del_z_sol)))
                dpx = pdum - del_psi_sol*real(i-1, lk)
                dpz = zdum - del_z_sol*real(j-1, lk)
                dx2 = dpx*dpx
                dz2 = dpz*dpz
              
                sdum0 = sspsol_psi_z(1,i,j) + sspsol_psi_z(2,i,j)*dpx + sspsol_psi_z(3,i,j)*dx2 +&
                            (sspsol_psi_z(4,i,j) + sspsol_psi_z(5,i,j)*dpx + sspsol_psi_z(6,i,j)*dx2)*dpz+&
                            (sspsol_psi_z(7,i,j) + sspsol_psi_z(8,i,j)*dpx + sspsol_psi_z(9,i,j)*dx2)*dz2
                sdum0=max(0.0_lk,min(1.0_lk,sdum0))
                ds_tmp = abs(sdum0-Ssol(igrid_sol(np)+1:igrid_sol(np)+mpara_sol(np)))
                is0=minloc(ds_tmp)
                ds0=sdum0-Ssol(igrid_sol(np)+is0(1))
                if(ds0>0.0_lk)then
                    jj0=min(is0(1),mpara_sol(np)-1)
                else
                    jj0=max(is0(1)-1,1)
                endif
                ij0_rg = mgrid_core+igrid_sol(np)+jj0
            endif
        endif

!TODO: algorithm for np not valid in core
        if(sfcore==1)then
            if((psifrc0 .le. psi1_core) .or. (sfsol==0))then
                if((psifrc0 .ge. psi1_core) .and. sfsol==0)then
                    np=mpsi_core-1
                    psifrc=psi1_core
                elseif(psifrc0 .le. psi0_core)then
                    np=1
                    psifrc=psi0_core
                else
                    left = 1
                    right = mpsi_core
                    do while(right-left>1)
                      i=floor((right-left)/2.0)+left
                      if(psifrc0>=psimesh_core(i) .and. psifrc0<psimesh_core(i+1))then
                        left=i
                        exit
                      elseif(psifrc0<psimesh_core(i))then
                        right=i
                      else
                        left=i+1
                      endif
                    enddo
                    np = left
                    psifrc = psifrc0
                endif
                tdum=zpart(3,m)/sqrt((rdum-1.0_lk)**2+zpart(3,m)**2+1d-20)
                tdum=dasin(tdum)*sign((rdum-1.0_lk),1.0_lk)-pi*(sign(rdum-1.0_lk,1.0_lk)*0.5_lk+sign(rdum-1.0_lk,1.0_lk)-1.5_lk)
                pdum=max(0.0_lk,psifrc-msp_psi_core(1))
                zdum=max(0.0_lk,tdum)
                i = max(1, min(lspcore-1, ceiling(pdum/del_psi_core)))
                j = max(1, min(lstcore-1, ceiling(zdum/del_theta_core)))
                dpx = pdum - del_psi_core*real(i-1, lk)
                dpz = zdum - del_theta_core*real(j-1, lk)
                dx2 = dpx*dpx
                dz2 = dpz*dpz

                sdum0 = sspcore_psi_theta(1,i,j) + sspcore_psi_theta(2,i,j)*dpx + sspcore_psi_theta(3,i,j)*dx2 +&
                            (sspcore_psi_theta(4,i,j) + sspcore_psi_theta(5,i,j)*dpx + sspcore_psi_theta(6,i,j)*dx2)*dpz+&
                            (sspcore_psi_theta(7,i,j) + sspcore_psi_theta(8,i,j)*dpx + sspcore_psi_theta(9,i,j)*dx2)*dz2

                ! find the closest grid to the guiding center
                ns0=sdum0*dels_inv(np)
                ij0_rg = igrid_core(np)+max(1,min(mpara_core(np)-1,ceiling(ns0-0.5_lk)))
            endif
        endif

    !!! compute weight for each gyro-averaging point
        do igyro=1,ngyro

            ! compute the delta_zeta corresponding to gyroradius 
            !rhogyropart(igyro,m)=zpart(6,m)*sqrt(apart)/(rho0*qpart)*zetagyro(igyro,ij0_rg)/rdum
            rhogyropart(igyro,m)=rho_ratio*zetagyro(igyro,ij0_rg)/rdum

            !rhogyropart(igyro,m)=2.0_lk*dbesjn(0,nmode*rhogyropart(igyro,m))-1.0_lk
            rhogyropart(igyro,m)=bessel_j0(nmode*rhogyropart(igyro,m)*n_modes(1)/0.7071067812_lk)

            ! position of each gyro-averaging point, the gyroradius is approximated at the closest grid computed above.
            psifrc = psifrc0+rho_ratio*gc_coef*pgyro(igyro,ij0_rg)

            if(sfsol==1)then
                if((psifrc .ge. psi_sep) .or. sfcore==0)then     
                    if(psifrc .ge. psi1_sol)then
                        wppart(igyro,m)=1.0_lk
                        np=mpsi_sol-1
                        psifrc=psi1_sol
                    elseif(psifrc .le. psi0_sol .and. sfcore==0)then
                        wppart(igyro,m)=0.0_lk
                        np=1
                        psifrc=psi0_sol
                    else
                        left = 1
                        right = mpsi_sol
                        do while(right-left>1)
                          i=floor((right-left)/2.0)+left
                          if(psifrc>=psimesh_sol(i) .and. psifrc<psimesh_sol(i+1))then
                            left=i
                            exit
                          elseif(psifrc<psimesh_sol(i))then
                            right=i
                          else
                            left=i+1
                          endif
                        enddo
                        np = left
                        wppart(igyro,m) = (psifrc-psimesh_sol(np))/(psimesh_sol(np+1)-psimesh_sol(np))
                    endif

                    jppart(igyro,m)=np+mpsi_core
                    pdum=max(0.0_lk,psifrc-msp_psi_sol(1))
                    zdum=max(0.0_lk,zpart(3,m)-Zsol_domain0)
                    i = max(1, min(lspsol-1, ceiling(pdum/del_psi_sol)))
                    j = max(1, min(lszsol-1, ceiling(zdum/del_z_sol)))
                    dpx = pdum - del_psi_sol*real(i-1, lk)
                    dpz = zdum - del_z_sol*real(j-1, lk)
                    dx2 = dpx*dpx
                    dz2 = dpz*dpz
    
                    sdum0 = sspsol_psi_z(1,i,j) + sspsol_psi_z(2,i,j)*dpx + sspsol_psi_z(3,i,j)*dx2 +&
                                (sspsol_psi_z(4,i,j) + sspsol_psi_z(5,i,j)*dpx + sspsol_psi_z(6,i,j)*dx2)*dpz+&
                                (sspsol_psi_z(7,i,j) + sspsol_psi_z(8,i,j)*dpx + sspsol_psi_z(9,i,j)*dx2)*dz2

                    sdum0=max(0.0_lk,min(1.0_lk,sdum0))
                    sdum1=sdum0
        
                    ds_tmp=abs(sdum0-Ssol(igrid_sol(np)+1:igrid_sol(np)+mpara_sol(np)))
                    is0=minloc(ds_tmp)
                    ds0=sdum0-Ssol(igrid_sol(np)+is0(1))
                    if(ds0>0.0_lk)then
                        jj0=min(is0(1),mpara_sol(np)-1)
                    else
                        jj0=max(is0(1)-1,1)
                    endif

                    ds_tmp=abs(sdum1-Ssol(igrid_sol(np+1)+1:igrid_sol(np+1)+mpara_sol(np+1)))
                    is1=minloc(ds_tmp)
                    ds1=sdum1-Ssol(igrid_sol(np+1)+is1(1))
                    if(ds1>0.0_lk)then
                        jj1=min(is1(1),mpara_sol(np+1)-1)
                    else
                        jj1=max(is1(1)-1,1)
                    endif
                    ij0=igrid_sol(np)+jj0
                    ij1=igrid_sol(np+1)+jj1

                    wtpart0(igyro,m)=(sdum0-Ssol(ij0))/(Ssol(ij0+1)-Ssol(ij0))   
                    wtpart1(igyro,m)=(sdum1-Ssol(ij1))/(Ssol(ij1+1)-Ssol(ij1))
                    jtpart0(igyro,m)=mgrid_core+ij0
                    jtpart1(igyro,m)=mgrid_core+ij1
                endif
            endif

!TODO: algorithm for np not valid in core
            if(sfcore==1)then
                if((psifrc .le. psi1_core) .or. (sfsol==0))then
                    if((psifrc .ge. psi1_core) .and. sfsol==0)then
                        wppart(igyro,m)=1.0_lk
                        np=mpsi_core-1
                        psifrc=psi1_core
                    elseif(psifrc .le. psi0_core)then
                        wppart(igyro,m)=0.0_lk
                        np=1
                        psifrc=psi0_core
                    else
                        left = 1
                        right = mpsi_core
                        do while(right-left>1)
                          i=floor((right-left)/2.0)+left
                          if(psifrc>=psimesh_core(i) .and. psifrc<psimesh_core(i+1))then
                            left=i
                            exit
                          elseif(psifrc<psimesh_core(i))then
                            right=i
                          else
                            left=i+1
                          endif
                        enddo
                        np=i
                        wppart(igyro,m)=(psifrc-psimesh_core(np))/(psimesh_core(np+1)-psimesh_core(np))
                    endif
    
                    jppart(igyro,m)=np
                    tdum=zpart(3,m)/sqrt((rdum-1.0_lk)**2+zpart(3,m)**2+1d-20)
                    tdum=dasin(tdum)*sign((rdum-1.0_lk),1.0_lk)-pi*(sign(rdum-1.0_lk,1.0_lk)*0.5_lk+sign(rdum-1.0_lk,1.0_lk)-1.5_lk)
                    pdum=max(0.0_lk,psifrc-msp_psi_core(1))
                    zdum=max(0.0_lk,tdum)
                    i = max(1, min(lspcore-1, ceiling(pdum/del_psi_core)))
                    j = max(1, min(lstcore-1, ceiling(zdum/del_theta_core)))
                    dpx = pdum - del_psi_core*real(i-1, lk)
                    dpz = zdum - del_theta_core*real(j-1, lk)
                    dx2 = dpx*dpx
                    dz2 = dpz*dpz
                  
                    sdum0 = sspcore_psi_theta(1,i,j) + sspcore_psi_theta(2,i,j)*dpx + sspcore_psi_theta(3,i,j)*dx2 +&
                                (sspcore_psi_theta(4,i,j) + sspcore_psi_theta(5,i,j)*dpx + sspcore_psi_theta(6,i,j)*dx2)*dpz+&
                                (sspcore_psi_theta(7,i,j) + sspcore_psi_theta(8,i,j)*dpx + sspcore_psi_theta(9,i,j)*dx2)*dz2
                    sdum1=sdum0
        
                    ns0=sdum0*dels_inv(np)
                    ns1=sdum1*dels_inv(np+1)
        
                    jtpart0(igyro,m)=igrid_core(np)+max(1,min(ceiling(ns0),mpara_core(np)-1))
                    jtpart1(igyro,m)=igrid_core(np+1)+max(1,min(ceiling(ns1),mpara_core(np+1)-1))
                    if(ns0 .lt. real(mpara_core(np)-1))then
                    wtpart0(igyro,m)=ns0-int(ns0)
                    else
                    wtpart0(igyro,m)=1.0_lk
                    endif
        
                    if(ns1 .lt. real(mpara_core(np+1)-1))then
                    wtpart1(igyro,m)=ns1-int(ns1)
                    else
                    wtpart1(igyro,m)=1.0_lk
                    endif
                endif
            endif

        enddo
    enddo
    !$acc end parallel

    if(sfcore==1)deallocate(dels,dels_inv)

end subroutine locate_kernel
