! find ion location for interpolation in both gathering (chargei.F90) and scattering (pushi.F90)
subroutine chargefki
    use global_parameters
    use read_in_parameters
    use particle_array
    implicit none

    interface
    subroutine chargefk(qpart,apart,density,flow,zpart,wppart,wtpart0,wtpart1,jtpart0,jtpart1,wzpart,&
                        zonal,zonalc,marker,pmark,mp,pressurepara,pressureperp)
        use global_parameters
        use grids_fields
        use grids_setup
        use read_in_parameters, only: min_r,frcfullf
        use equilibrium, only: psirz_Eq, lsr, lsz, spdR, spdZ, zmax
        use simulation_domain, only: sfcore, sfsol, couplefrc
        use precision
        implicit none
    
        real(lk) :: qpart, apart
        real(lk), dimension(:) :: zonal,zonalc,pmark,wzpart,marker
        integer,dimension(:,:) :: jtpart0,jtpart1
        integer mp
        real(lk),dimension(:,:) :: wppart,wtpart0,wtpart1,zpart
        real(lk),dimension(0:,:) :: density,flow,pressurepara,pressureperp
    
    
        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,&
          nrsp,nzsp
    
        real(lk) :: gyrodum,weight,rdum,rdum_inv,zdum,upara,energypara,energyperp,cmratio,wz1,wz0,wp1,wp0,wt11,wt10,wt01,wt00,&
             sendl(mgrid,4),recvr(mgrid,4),dntmp(0:1,mgrid),djtmp(0:1,mgrid),dpparatmp(0:1,mgrid),dpperptmp(0:1,mgrid),&
             zf0,zc0,temp,adum(mpsi),&
             dpr,dr2,dpz,dz2,gradpsi_r,gradpsi_z,br,bz,b,spdR_inv,spdZ_inv
        real(lk) dbesjn, zeta_gyro_avg, total_density, avg_density
    end subroutine chargefk
    end interface

    call chargefk(qfki,afki,densityfki,flowfki,zfki,wpfki,wtfki0,wtfki1,jtfki0,jtfki1,wzfki,&
                  zonali,zonalci,markerfki,pmarkfki,mfki,pressureipara,pressureiperp)

end subroutine chargefki
subroutine locatefki
    use global_parameters
    use particle_array
    use read_in_parameters
    implicit none
    interface
    subroutine locatefk(zpart,wzpart,wppart,wtpart0,wtpart1,jtpart0,&
    jtpart1,qpart,apart,ngyro,mp,jppart)
    use global_parameters
    use grids_fields
    use grids_setup
    use read_in_parameters
    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 precision
    implicit none
  
    real(lk), dimension(:,:) :: zpart, wppart, wtpart0, wtpart1
    integer, dimension(:,:) :: jtpart0, jtpart1, jppart
    real(lk), dimension(:) :: wzpart
    real(lk) :: qpart, apart
    integer :: mp, ngyro


    integer :: m,nrsp,nzsp,ij,np,i,iz0(1),iz1(1),jj0,jj1,ij0,ij1,is0(1),is1(1)
    real(lk) :: rdum,zdum,zetadum,delz,dpr,dpz,dr2,dz2,psifrc,sign_psi,spdR_inv,spdZ_inv,tdum,&
             sdum,sdum0,sdum1,ns0,ns1,&
             dz0,dz1,ds0,ds1
    real(lk),dimension(:),allocatable :: dels,dels_inv
    end subroutine locatefk
    end interface
    call locatefk(zfki, wzfki, wpfki, wtfki0, wtfki1, jtfki0, jtfki1, qfki, afki, 1, mfki, jpfki)
end subroutine locatefki


subroutine locatefk(zpart,wzpart,wppart,wtpart0,wtpart1,jtpart0,&
    jtpart1,qpart,apart,ngyro,mp,jppart)
    use global_parameters
    use grids_fields
    use grids_setup
    use read_in_parameters
    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 precision
    implicit none
  
    real(lk), dimension(:,:) :: zpart, wppart, wtpart0, wtpart1
    integer, dimension(:,:) :: jtpart0, jtpart1, jppart
    real(lk), dimension(:) :: wzpart
    real(lk) :: qpart, apart
    integer :: mp, ngyro


    integer :: m,nrsp,nzsp,ij,np,i,iz0(1),iz1(1),jj0,jj1,ij0,ij1,is0(1),is1(1)
    real(lk) :: rdum,zdum,zetadum,delz,dpr,dpz,dr2,dz2,psifrc,sign_psi,spdR_inv,spdZ_inv,tdum,&
             sdum,sdum0,sdum1,ns0,ns1,&
             dz0,dz1,ds0,ds1,psidum
    real(lk),dimension(:),allocatable :: dels,dels_inv
    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(:)
    endif
!$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& jj0,jj1,ij0,ij1,is0,is1,ds0,ds1,psidum)
    do m=1,mp
        psidum=spline2d(zpart(1,m),zpart(3,m)+zmax,0,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)
        if(psidum>psi1_sol) cycle
        if(psidum<psi0_core) cycle
        if(sfsol==0 .and. psidum>psi1_core) cycle
        if(sfcore==0 .and. psidum<psi0_sol) cycle
        rdum=max(zpart(1,m),min_r)
        zdum=zpart(3,m)
     
        zetadum=zpart(2,m)
  
        wzpart(m)=(zetadum-zeta0)*delz  !weight for upper toroidal grid
  
! psi value at particle position by using spline function from equilibrium  
        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
  
        psifrc=psirz_Eq(1,nrsp,nzsp)     + psirz_Eq(2,nrsp,nzsp)*dpr     + psirz_Eq(3,nrsp,nzsp)*dr2 &
              +psirz_Eq(4,nrsp,nzsp)*dpz + psirz_Eq(5,nrsp,nzsp)*dpr*dpz + psirz_Eq(6,nrsp,nzsp)*dr2*dpz &
              +psirz_Eq(7,nrsp,nzsp)*dz2 + psirz_Eq(8,nrsp,nzsp)*dpr*dz2 + psirz_Eq(9,nrsp,nzsp)*dr2*dz2
  
  
        if(sfsol==1)then
            if((psifrc .ge. psi_sep) .or. sfcore==0)then     
                if(psifrc .ge. psi1_sol)then
                    wppart(1,m)=1.0_lk
                    np=mpsi_sol-1
                    psifrc=psi1_sol
                elseif(psifrc .le. psi0_sol .and. sfcore==0)then
                    wppart(1,m)=0.0_lk
                    np=1
                    psifrc=psi0_sol
                else
                    do i=1,mpsi_sol-1
                        sign_psi=(psifrc-psimesh_sol(i))*(psifrc-psimesh_sol(i+1))
                        if(sign_psi .le. 0.0_lk)then
                            np=i
                            wppart(1,m)=(psifrc-psimesh_sol(np))/(psimesh_sol(np+1)-psimesh_sol(np))
                            exit
                        endif
                    enddo
                endif
  
                jppart(1,m)=np+mpsi_core
       
                sdum0=spline2d(psifrc-msp_psi_sol(1),zdum-Zsol_domain0,0,lspsol,lszsol,del_psi_sol,del_z_sol,sspsol_psi_z,0,0)
                sdum0=max(0.0_lk,min(1.0_lk,sdum0))
                sdum1=sdum0
       
                is0=minloc(abs(sdum0-Ssol(igrid_sol(np)+1:igrid_sol(np)+mpara_sol(np))))
                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
       
                is1=minloc(abs(sdum1-Ssol(igrid_sol(np+1)+1:igrid_sol(np+1)+mpara_sol(np+1))))
                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(1,m)=(sdum0-Ssol(ij0))/(Ssol(ij0+1)-Ssol(ij0))   
                wtpart1(1,m)=(sdum1-Ssol(ij1))/(Ssol(ij1+1)-Ssol(ij1))
                jtpart0(1,m)=mgrid_core+ij0
                jtpart1(1,m)=mgrid_core+ij1
            endif
        endif
  
        if(sfcore==1)then
            if((psifrc .le. psi1_core) .or. (sfsol==0))then
                if((psifrc .ge. psi1_core) .and. sfsol==0)then
                    wppart(1,m)=1.0_lk
                    np=mpsi_core-1
                    psifrc=psi1_core
                elseif(psifrc .le. psi0_core)then
                    wppart(1,m)=0.0_lk
                    np=1
                    psifrc=psi0_core
                else
                    do i=1,mpsi_core-1
                        sign_psi=(psifrc-psimesh_core(i))*(psifrc-psimesh_core(i+1))
                        if(sign_psi .le. 0.0_lk)then
                            np=i
                            wppart(1,m)=(psifrc-psimesh_core(np))/(psimesh_core(np+1)-psimesh_core(np))
                            exit
                        endif
                    enddo
                endif
   
                jppart(1,m)=np
       
                tdum=geotheta(rdum,zdum)
       
                sdum0=spline2d(psifrc-msp_psi_core(1),tdum,0,lspcore,lstcore,del_psi_core,del_theta_core,sspcore_psi_theta,0,0)
                sdum1=sdum0
       
                ns0=sdum0*dels_inv(np)
                ns1=sdum1*dels_inv(np+1)
       
                jtpart0(1,m)=igrid_core(np)+max(1,min(ceiling(ns0),mpara_core(np)-1))
                jtpart1(1,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(1,m)=ns0-int(ns0)
                else
                  wtpart0(1,m)=1.0_lk
                endif
       
                if(ns1 .lt. real(mpara_core(np+1)-1))then
                  wtpart1(1,m)=ns1-int(ns1)
                else
                  wtpart1(1,m)=1.0_lk
                endif
            endif
        endif
    enddo
  
    if(sfcore==1)deallocate(dels,dels_inv)
  
end subroutine locatefk

subroutine chargefk(qpart,apart,density,flow,zpart,wppart,wtpart0,wtpart1,jtpart0,jtpart1,wzpart,&
                    zonal,zonalc,marker,pmark,mp,pressurepara,pressureperp)
    use global_parameters
    use grids_fields
    use grids_setup
    use read_in_parameters, only: min_r,frcfullf
    use equilibrium, only: psirz_Eq, lsr, lsz, spdR, spdZ, zmax
    use simulation_domain, only: sfcore, sfsol, couplefrc,psi0_core,psi1_core,psi0_sol,psi1_sol
    use spline_function, only: spline2d
    use precision
    implicit none

    real(lk) :: qpart, apart
    real(lk), dimension(:) :: zonal,zonalc,pmark,wzpart,marker
    integer,dimension(:,:) :: jtpart0,jtpart1
    integer mp
    real(lk),dimension(:,:) :: wppart,wtpart0,wtpart1,zpart
    real(lk),dimension(0:,:) :: density,flow,pressurepara,pressureperp
    
    
    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,&
      nrsp,nzsp,m_debug
      
    real(lk) :: gyrodum,weight,rdum,rdum_inv,zdum,upara,energypara,energyperp,cmratio,wz1,wz0,wp1,wp0,wt11,wt10,wt01,wt00,&
         sendl(mgrid,4),recvr(mgrid,4),dntmp(0:1,mgrid),djtmp(0:1,mgrid),dpparatmp(0:1,mgrid),dpperptmp(0:1,mgrid),&
         zf0,zc0,temp,adum(mpsi),psidum,&
         dpr,dr2,dpz,dz2,gradpsi_r,gradpsi_z,br,bz,b,spdR_inv,spdZ_inv
    real(lk) dbesjn, zeta_gyro_avg, total_density, avg_density
  
    cmratio=qpart/apart
    spdR_inv=1.0_lk/spdR
    spdZ_inv=1.0_lk/spdZ
  
    density=0.0_lk
    flow=0.0_lk
    pressurepara=0.0_lk
    pressureperp=0.0_lk

#ifdef _OPENMP
! Using OpenMP for gathering operations needs to address issues of memory conflict.
! The following lines are OpenMP directives for loop-level parallelism
! on shared memory machines (see http://www.openmp.org).
! When compiling with the OpenMP option, the preprocessor variable
! _OPENMP is automatically defined and we use it here to add pieces of codes
! needed to avoid contentions between threads. In the following parallel
! loop each thread has its own private array "dntmp()" which is used to
! prevent the thread from writing into the shared array "densityi()"
!$omp parallel private(dntmp,djtmp,dpparatmp,dpperptmp)
    dntmp=0.0_lk   ! Set array elements to zero
    djtmp=0.0_lk   ! Set array elements to zero
    dpparatmp=0.0_lk  ! Set array elements to zero
    dpperptmp=0.0_lk  ! Set array elements to zero
#endif

! scatter ion density for electrostatic simulation
!$omp do private(m,rdum,rdum_inv,zdum,nrsp,dpr,dr2,nzsp,dpz,dz2,&
!$omp& gradpsi_r,gradpsi_z,br,bz,b,&
!$omp& upara,energypara,energyperp,igyro,weight,wz1,wz0,wp1,wp0,wt10,wt00,wt11,wt01,ij,psidum)
    do m=1,mp
        psidum=spline2d(zpart(1,m),zpart(3,m)+zmax,0,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)
        if(psidum>psi1_sol) cycle
        if(psidum<psi0_core) cycle
        if(sfsol==0 .and. psidum>psi1_core) cycle
        if(sfcore==0 .and. psidum<psi0_sol) cycle

        rdum = max(zpart(11,m),min_r)
        rdum_inv=1.0_lk/rdum
        zdum = zpart(13,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)*br + zpart(8,m)*bz
        energypara=apart*upara*upara
        energyperp=0.5_lk*apart*((zpart(6,m)*zpart(11,m))**2+(zpart(8,m)*br/b)**2+&
                   (zpart(4,m)*bz/b)**2-2.0_lk*zpart(4,m)*zpart(8,m)*br*bz/b**2)
        weight=zpart(5,m)!*zpart0(7,m)

        wz1=weight*wzpart(m)         !weight for upper toroidal grid
        wz0=weight-wz1              !weight for lower toroidal grid 
        do igyro=1,1
            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

#ifdef _OPENMP
! Use thread-private temporary array dntmp to store the results
            ij=jtpart0(igyro,m)
            dntmp(0,ij) = dntmp(0,ij) + wz0*wt00
            dntmp(1,ij) = dntmp(1,ij) + wz1*wt00
            djtmp(0,ij) = djtmp(0,ij) + wz0*wt00*upara
            djtmp(1,ij) = djtmp(1,ij) + wz1*wt00*upara
            dpparatmp(0,ij) = dpparatmp(0,ij) + wz0*wt00*energypara
            dpparatmp(1,ij) = dpparatmp(1,ij) + wz1*wt00*energypara
            dpperptmp(0,ij) = dpperptmp(0,ij) + wz0*wt00*energyperp
            dpperptmp(1,ij) = dpperptmp(1,ij) + wz1*wt00*energyperp
 
            ij=ij+1
            dntmp(0,ij) = dntmp(0,ij) + wz0*wt10
            dntmp(1,ij) = dntmp(1,ij) + wz1*wt10
            djtmp(0,ij) = djtmp(0,ij) + wz0*wt10*upara
            djtmp(1,ij) = djtmp(1,ij) + wz1*wt10*upara
            dpparatmp(0,ij) = dpparatmp(0,ij) + wz0*wt10*energypara
            dpparatmp(1,ij) = dpparatmp(1,ij) + wz1*wt10*energypara
            dpperptmp(0,ij) = dpperptmp(0,ij) + wz0*wt10*energyperp
            dpperptmp(1,ij) = dpperptmp(1,ij) + wz1*wt10*energyperp
 
            ij=jtpart1(igyro,m)
            dntmp(0,ij) = dntmp(0,ij) + wz0*wt01
            dntmp(1,ij) = dntmp(1,ij) + wz1*wt01
            djtmp(0,ij) = djtmp(0,ij) + wz0*wt01*upara
            djtmp(1,ij) = djtmp(1,ij) + wz1*wt01*upara
            dpparatmp(0,ij) = dpparatmp(0,ij) + wz0*wt01*energypara
            dpparatmp(1,ij) = dpparatmp(1,ij) + wz1*wt01*energypara
            dpperptmp(0,ij) = dpperptmp(0,ij) + wz0*wt01*energyperp
            dpperptmp(1,ij) = dpperptmp(1,ij) + wz1*wt01*energyperp
 
            ij=ij+1
            dntmp(0,ij) = dntmp(0,ij) + wz0*wt11
            dntmp(1,ij) = dntmp(1,ij) + wz1*wt11
            djtmp(0,ij) = djtmp(0,ij) + wz0*wt11*upara
            djtmp(1,ij) = djtmp(1,ij) + wz1*wt11*upara
            dpparatmp(0,ij) = dpparatmp(0,ij) + wz0*wt11*energypara
            dpparatmp(1,ij) = dpparatmp(1,ij) + wz1*wt11*energypara
            dpperptmp(0,ij) = dpperptmp(0,ij) + wz0*wt11*energyperp
            dpperptmp(1,ij) = dpperptmp(1,ij) + wz1*wt11*energyperp
#else
! If no loop-level parallelism, write directly into array "densityi()".
            ij=jtpart0(igyro,m)                         !lower poloidal grid on inner flux surface
            density(0,ij) = density(0,ij) + wz0*wt00 !lower toroidal grid
            density(1,ij) = density(1,ij) + wz1*wt00 !upper toroidal grid
            flow(0,ij) = flow(0,ij) + wz0*wt00*upara !lower toroidal grid
            flow(1,ij) = flow(1,ij) + wz1*wt00*upara !upper toroidal grid
            pressurepara(0,ij) = pressurepara(0,ij) + wz0*wt00*energypara
            pressurepara(1,ij) = pressurepara(1,ij) + wz1*wt00*energypara
            pressureperp(0,ij) = pressureperp(0,ij) + wz0*wt00*energyperp
            pressureperp(1,ij) = pressureperp(1,ij) + wz1*wt00*energyperp
 
            ij=ij+1                                    !upper poloidal grid on inner flux surface
            density(0,ij) = density(0,ij) + wz0*wt10
            density(1,ij) = density(1,ij) + wz1*wt10
            flow(0,ij) = flow(0,ij) + wz0*wt10*upara
            flow(1,ij) = flow(1,ij) + wz1*wt10*upara
            pressurepara(0,ij) = pressurepara(0,ij) + wz0*wt10*energypara
            pressurepara(1,ij) = pressurepara(1,ij) + wz1*wt10*energypara
            pressureperp(0,ij) = pressureperp(0,ij) + wz0*wt10*energyperp
            pressureperp(1,ij) = pressureperp(1,ij) + wz1*wt10*energyperp
 
            ij=jtpart1(igyro,m)                         !lower poloidal grid on outer flux surface
            density(0,ij) = density(0,ij) + wz0*wt01
            density(1,ij) = density(1,ij) + wz1*wt01
            flow(0,ij) = flow(0,ij) + wz0*wt01*upara
            flow(1,ij) = flow(1,ij) + wz1*wt01*upara
            pressurepara(0,ij) = pressurepara(0,ij) + wz0*wt01*energypara
            pressurepara(1,ij) = pressurepara(1,ij) + wz1*wt01*energypara
            pressureperp(0,ij) = pressureperp(0,ij) + wz0*wt01*energyperp
            pressureperp(1,ij) = pressureperp(1,ij) + wz1*wt01*energyperp
 
            ij=ij+1                                    !upper poloidal grid on outer flux surface
            density(0,ij) = density(0,ij) + wz0*wt11
            density(1,ij) = density(1,ij) + wz1*wt11
            flow(0,ij) = flow(0,ij) + wz0*wt11*upara
            flow(1,ij) = flow(1,ij) + wz1*wt11*upara
            pressurepara(0,ij) = pressurepara(0,ij) + wz0*wt11*energypara
            pressurepara(1,ij) = pressurepara(1,ij) + wz1*wt11*energypara
            pressureperp(0,ij) = pressureperp(0,ij) + wz0*wt11*energyperp
            pressureperp(1,ij) = pressureperp(1,ij) + wz1*wt11*energyperp

#endif
        enddo
    enddo
#ifdef _OPENMP
! For OpenMP, we now write the results accumulated in each thread-private
! array dntmp() back into the shared array densitye(). The loop is enclosed
! in a critical section so that one thread at a time updates densitye(). 
!$omp critical
    do ij=1,mgrid
        density(0:1,ij) = density(0:1,ij) + dntmp(0:1,ij)
        flow(0:1,ij) = flow(0:1,ij) + djtmp(0:1,ij)
        pressurepara(0:1,ij)=pressurepara(0:1,ij)+dpparatmp(0:1,ij)
        pressureperp(0:1,ij)=pressureperp(0:1,ij)+dpperptmp(0:1,ij)
    enddo
!$omp end critical
#endif
!$omp end parallel

! 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
            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(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,:)
    recvr=0.0_lk
    icount=4*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)
    gyrodum=1.0_lk
    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
                zonal(i)=zonal(i)+density(1,ij)
                zonalc(i)=zonalc(i)+flow(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
                zonal(i)=zonal(i)+density(1,ij)
                zonalc(i)=zonalc(i)+flow(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(flow)
        call core_para_periodicity(pressurepara)
        call core_para_periodicity(pressureperp)
    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(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)
                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)
                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(flow)
    CALL TOROIDAL_FILTER(density)
    CALL TOROIDAL_FILTER(flow)
end subroutine chargefk
