subroutine pushfki
  use global_parameters
  use particle_array
  use grids_eqdata
  use read_in_parameters
  
  implicit none
  interface
   subroutine pushfk(zpart,zpart0,wppart,wtpart0,wtpart1,&
    jtpart0,jtpart1,wzpart,mesht1d,qpart,apart,mp,kapanp,kapatp,diagpart)
      use precision
      use particle_array, only: mpdiag
      implicit none
  
     !declaration of the dummy arguments
      integer ngyro,mp
      integer,dimension(:,:) :: jtpart0,jtpart1
      real(lk) qpart,apart
      real(lk),dimension(:) :: wzpart
      real(lk),dimension(:) :: mesht1d,kapanp,kapatp
      real(lk),dimension(:,:) :: zpart,zpart0,wppart,wtpart0,wtpart1
      real(lk),dimension(mpdiag) :: diagpart
    end subroutine pushfk
  end interface
  call pushfk(zfki, zfki0, wpfki, wtfki0, wtfki1,&
  jtfki0, jtfki1, wzfki, meshti1d, qfki, afki, mfki, kapanip, kapatip,diagfki)
end subroutine pushfki
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine pushfk(zpart,zpart0,wppart,wtpart0,wtpart1,&
  jtpart0,jtpart1,wzpart,mesht1d,qpart,apart,mp,kapanp,kapatp,diagpart)
  use precision
  use global_parameters, only: istep, irk, torbound, pi2
  use simulation_domain, only: Zsol_domain_width, &
                               psi0_core, psi1_core, psi0_sol, psi1_sol, &
                               sfcore, sfsol
  use spline_function, only: spline2d,sppsiR1d,spline1d
  use grids_fields
  use grids_setup
  use grids_eqdata
  use read_in_parameters, only: nonlinear, tstep, Zsol_domain0, frcfullf
  use particle_array, only: mpdiag
  use spline_function, only: spline2d,sppsiR1d,spline1d
  use cylindrical_function, only: geotheta
  use equilibrium, only: psirz_Eq, r_Eq, z_Eq,&
                         psi_sep, lsr, lsz, spdR, spdZ, zmax, B0R_Eq, B0Z_Eq, &
                         B0mag_Eq
  use mapping_splines, only: lsp_er,dsp_er,phip00_sp
  implicit none

  !declaration of the dummy arguments
  integer ngyro,mp
  integer,dimension(:,:) :: jtpart0,jtpart1
  real(lk) qpart,apart
  real(lk),dimension(:) :: wzpart
  real(lk),dimension(:) :: mesht1d, kapanp,kapatp
  real(lk),dimension(:,:) :: zpart,zpart0,wppart,wtpart0,wtpart1
  real(lk),dimension(mpdiag) :: diagpart


  ! TODO: confirm the variables here
  integer i,j,ij,m,igyro,ierror,ipsi(1),ii,nzsp,nzsp1,nrsp,nrsp1,iout,icount
  real(lk) dtime,wz1,wz0,wp1,wp0,wt11,wt10,wt01,wt00,&
       e1,e2,e3,cmratio,c_inv,wpfk(3,mp),rdot,zdot,wdot,perturb,&
       dpr,dr2,dpz,dz2,acoustic_v(3,mp),&
       br,bz,b,b_inv,bt,dr21,dp1,omegac,&
       spdR_inv,spdZ_inv,rdum_inv,gradphi_dum(3),&
       rtmp,ztmp,rdum0,tdum0,zdum0,rdum1,zdum1,wdrive,t_inv,kappanp,&
       kappatp,t_cy,psidum,vecb(3),vec1(3),vec2(3),vecrho(3),dpsi,vecv(3),&
       angmom,vdr,diagtmp(mpdiag),deltaf
  real(lk) rdum3,energy,&
       tandalpha1,tandalpha2,tandalpha_inv,ur,uz,&
       deltheta1,deltheta2,deltheta3,deltheta4,tdot,&
       ut,ttmp,gqi,gqi_inv,gdtt,gdtt1,gutt,gutt1,gdtt3,gutt3,&
       canmom,psitmp,dr1,gradpsi_r,gradpsi_z,&
       dz1,dz21
  
  cmratio=qpart/apart
  c_inv=1.0_lk/qpart
  perturb=real(nonlinear)
  !if(track_particles==1) perturb=1.0
  spdR_inv=1.0_lk/spdR
  spdZ_inv=1.0_lk/spdZ
  dtime=tstep*0.5_lk*irk
  ngyro=1 !no gyro-averaging for fk particles
! electrostatic field
!$omp parallel do private(m,igyro,e1,e2,e3,wz1,wz0,wp1,wp0,wt00,wt10,wt01,wt11,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)then
      wpfk(:,m)=0.0_lk
      cycle
    endif
    if(psidum<psi0_core)then
      wpfk(:,m)=0.0_lk
      cycle
    endif
    if(sfsol==0 .and. psidum>psi1_core)then
      wpfk(:,m)=0.0_lk
      cycle
    endif
    if(sfcore==0 .and. psidum<psi0_sol)then
      wpfk(:,m)=0.0_lk
      cycle
    endif
    e1=0.0_lk
    e2=0.0_lk
    e3=0.0_lk

    wz1=wzpart(m)         !weight for upper toroidal grid
    wz0=1.0_lk-wz1              !weight for lower toroidal grid

    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 poloidal grid on inner flux surface
       wt00=wp0-wt10            !lower poloidal grid on inner flux surface

       wt11=wp1*wtpart1(igyro,m) !upper poloidal grid on outer flux surface
       wt01=wp1-wt11            !lower poloidal grid on outer flux surface  
       ij=jtpart0(igyro,m)       !lower poloidal grid on inner flux surface
       e1=e1+wt00*(wz0*gradphi(1,0,ij)+wz1*gradphi(1,1,ij))
       e2=e2+wt00*(wz0*gradphi(2,0,ij)+wz1*gradphi(2,1,ij))
       e3=e3+wt00*(wz0*gradphi(3,0,ij)+wz1*gradphi(3,1,ij))

       ij=ij+1                  !upper poloidal grid on inner flux surface
       e1=e1+wt10*(wz0*gradphi(1,0,ij)+wz1*gradphi(1,1,ij))
       e2=e2+wt10*(wz0*gradphi(2,0,ij)+wz1*gradphi(2,1,ij))
       e3=e3+wt10*(wz0*gradphi(3,0,ij)+wz1*gradphi(3,1,ij))

       ij=jtpart1(igyro,m)       !lower poloidal grid on outer flux surface
       e1=e1+wt01*(wz0*gradphi(1,0,ij)+wz1*gradphi(1,1,ij))
       e2=e2+wt01*(wz0*gradphi(2,0,ij)+wz1*gradphi(2,1,ij))
       e3=e3+wt01*(wz0*gradphi(3,0,ij)+wz1*gradphi(3,1,ij))

       ij=ij+1                  !upper poloidal grid on outer flux surface
       e1=e1+wt11*(wz0*gradphi(1,0,ij)+wz1*gradphi(1,1,ij))
       e2=e2+wt11*(wz0*gradphi(2,0,ij)+wz1*gradphi(2,1,ij))
       e3=e3+wt11*(wz0*gradphi(3,0,ij)+wz1*gradphi(3,1,ij))
    enddo
    wpfk(1,m)=e1 !0.0
    wpfk(2,m)=e2 !\partial \phi / \partial zeta
    wpfk(3,m)=e3 !\partial \phi / \partial s i.e. (gradphi)_\parallel
 enddo
  if(irk==1)then ! update velocity from t-dt/2 to t+dt/2 and a position guess at t+dt/2
!$omp parallel do private(m,rdum0,tdum0,zdum0,rdum1,rdum3,energy,dpz,dz2,&
!$omp& tandalpha1,tandalpha2,tandalpha_inv,b_inv,bt,bz,br,dpr,dr2,nrsp,&
!$omp& b,ur,uz,deltheta1,deltheta2,deltheta3,deltheta4,rdot,tdot,nzsp,&
!$omp& zdot,rtmp,ut,ttmp,ztmp,gqi,gqi_inv,gdtt,gutt,gdtt1,gutt1,&
!$omp& gdtt3,gutt3,canmom,nrsp1,nzsp1,zdum1,psitmp,gradpsi_r,gradpsi_z,&
!$omp& dr1,dr21,dz1,dz21,rdum_inv,gradphi_dum,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 .or. psidum<psi0_sol)then
         cycle
      endif
      rdum0=zpart0(11,m) !R
      rdum_inv=1.0_lk/rdum0
      tdum0=zpart0(12,m) !toroidal zeta
      zdum0=zpart0(13,m) !Z
      zdum0=modulo((zdum0+zmax),(zmax*2.0_lk))-zmax
      zpart0(5,m)=zpart(5,m)
      zpart0(4,m)=zpart(4,m)
      zpart0(6,m)=zpart(6,m)
      zpart0(8,m)=zpart(8,m)
      ! nearest grid index in r and z direction
      nrsp=max(1,min(lsr-1,int(ceiling(rdum0*spdR_inv))))
      dpr=rdum0-spdR*real(nrsp-1)
      dr2=dpr*dpr
      nzsp=max(1,min(lsz-1,int(ceiling((zdum0+zmax)*spdZ_inv))))
      dpz=zdum0+zmax-spdZ*real(nzsp-1)
      dz2=dpz*dpz

      ! componets of magnetic field at t
      gradpsi_r= psirz_Eq(2,nrsp,nzsp) + 2.0_lk*psirz_Eq(3,nrsp,nzsp)*dpr &
               + psirz_Eq(5,nrsp,nzsp)*dpz + 2.0_lk*psirz_Eq(6,nrsp,nzsp)*dpr*dpz &
               + psirz_Eq(8,nrsp,nzsp)*dz2 + 2.0_lk*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_lk*psirz_Eq(7,nrsp,nzsp)*dpz + 2.0_lk*psirz_Eq(8,nrsp,nzsp)*dpr*dpz + 2.0_lk*psirz_Eq(9,nrsp,nzsp)*dr2*dpz

      br=-gradpsi_z*rdum_inv
      bz=gradpsi_r*rdum_inv
      bt=0.0_lk
      b=sqrt(br**2+bz**2)
      b_inv=1.0_lk/b
      ! metric tensor
      gdtt=rdum0*rdum0
      gutt=1.0_lk/gdtt
      ! Temporaray storage of position      
      rdum1=zpart(11,m)
      zdum1=zpart(13,m)
      zdum1=modulo((zdum1+zmax),(zmax*2.0_lk))-zmax
      gdtt1=rdum1*rdum1
      gutt1=1.0_lk/gdtt1
      ! particle kinetic energy (t-1/2) for diagnostics
      energy=zpart(4,m)*zpart(4,m)+zpart(6,m)*zpart(6,m)*gdtt1+zpart(8,m)*zpart(8,m)

      ! Find nearest grid point in R to the particle position  
      nrsp1=max(1,min(lsr-1,int(ceiling((rdum1)*spdR_inv))))
      dr1=rdum1-spdR*real(nrsp1-1)
      dr21=dr1*dr1
      ! Find nearest grid point in Z to the particle position  
      nzsp1=max(1,min(lsz-1,int(ceiling((zdum1+zmax)*spdZ_inv))))
      dz1=zdum1+zmax-spdZ*real(nzsp1-1)
      dz21=dz1*dz1

      ! poloidal flux function at (t-1/2) for canonical angular momentum
      psitmp=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

      ! particle canonical angular momentum for diagnostics
      canmom=apart*zpart(6,m)*rdum1*rdum1+qpart*psitmp
      ! components of velocity
      deltheta1=cos(zpart(12,m))*cos(zpart0(12,m))+sin(zpart(12,m))*sin(zpart0(12,m))
      deltheta2=-sin(zpart(12,m))*cos(zpart0(12,m))+cos(zpart(12,m))*sin(zpart0(12,m))
      ! uR      
      rdot=(zpart(4,m)*deltheta1+rdum1*zpart(6,m)*deltheta2) 
      ! uzeta     
      tdot=(-zpart(4,m)*deltheta2*rdum0+zpart(6,m)*deltheta1*rdum0*rdum1)/gdtt  
      ! uZ     
      zdot=zpart(8,m)

      ! half acceleration
      rdot=rdot-cmratio*wpfk(3,m)*br*b_inv*perturb*tstep*0.5
      tdot=tdot-cmratio*wpfk(2,m)*perturb*tstep*0.5*rdum0
      zdot=zdot-cmratio*wpfk(3,m)*bz*b_inv*perturb*tstep*0.5
      ! velocity rotation
      tandalpha1=cmratio*tstep*0.5                                              !Q
      tandalpha2=1.0+tandalpha1*tandalpha1*b*b                              !1+T^2 
      tandalpha_inv=2.0/tandalpha2
      tandalpha_inv=tandalpha_inv*tandalpha1                                   ! P 
      gqi=sqrt(gdtt)
      gqi_inv=1.0/gqi

      ur=(1.0-tandalpha1*tandalpha_inv*(bz*bz+bt*bt))*rdot+&               ! R
         (tandalpha1*tandalpha_inv*br*bt*rdum0+tandalpha_inv*gdtt*gqi_inv*bz)*tdot+&
         (tandalpha1*tandalpha_inv*br*bz-tandalpha_inv*bt*rdum0*gqi_inv)*zdot

      ut=(1.0-tandalpha1*tandalpha_inv*(br*br+bz*bz))*tdot+&           ! vzeta
             (tandalpha1*tandalpha_inv*br*bt*rdum0*gutt-tandalpha_inv*gqi_inv*bz)*rdot+&
             (tandalpha1*tandalpha_inv*bz*bt*rdum0*gutt+tandalpha_inv*gqi_inv*br)*zdot     
   
      uz=(1.0-tandalpha1*tandalpha_inv*(br*br+bt*bt))*zdot+&           ! vZ
         (tandalpha1*tandalpha_inv*bz*bt*rdum0-tandalpha_inv*br*gdtt*gqi_inv)*tdot+&     
         (tandalpha1*tandalpha_inv*br*bz+tandalpha_inv*rdum0*bt*gqi_inv)*rdot

      ! half acceleration
      ur=ur-cmratio*wpfk(3,m)*br*b_inv*perturb*tstep*0.5
      ut=ut-cmratio*wpfk(2,m)*perturb*tstep*0.5*rdum0
      uz=uz-cmratio*wpfk(3,m)*bz*b_inv*perturb*tstep*0.5

      !! guess of particle position (t+1/2)
      ! R
      rtmp=zpart0(11,m)+ur*tstep*0.5
      ! toroidal zeta  
      ttmp=zpart0(12,m)+(ut)*tstep*0.5
      ! Z      
      ztmp=zpart0(13,m)+uz*tstep*0.5

      ! Global storage of particle position for next timestep initialization
      zpart(11,m)=rtmp
      zpart(12,m)=ttmp
      zpart(13,m)=ztmp
      zpart(13,m)=modulo((zpart(13,m)-Zsol_domain0), Zsol_domain_width)+Zsol_domain0
      rdum3=rtmp
      ! Metric tensor (t+1/2)
      gdtt3=rdum3*rdum3
      gutt3=1.0/gdtt3

      ! Return to components of v from components of u at (t+1/2)
      deltheta3=cos(zpart0(12,m))*cos(ttmp)+sin(zpart0(12,m))*sin(ttmp)
      deltheta4=-sin(zpart0(12,m))*cos(ttmp)+cos(zpart0(12,m))*sin(ttmp) 

      ! vR      
      zpart(4,m)=(ur*deltheta3+ut*deltheta4*rdum0) 
      ! vzeta    
      zpart(6,m)=(-ur*deltheta4*rdum3+ut*deltheta3*rdum0*rdum3)/gdtt3    
      ! vZ      
      zpart(8,m)=uz
      
      ! when irk==1, use v(t) as the velocity
      ! average velocity (t)
      ! vR
      zpart0(4,m)=(zpart0(4,m)+zpart(4,m))*0.5
      ! vzeta    
      zpart0(6,m)=(zpart0(6,m)+zpart(6,m))*0.5
      ! vZ    
      zpart0(8,m)=(zpart0(8,m)+zpart(8,m))*0.5
      zpart0(1:3,m)=zpart0(11:13,m)
      zpart0(10,m)=zpart(10,m)
      zpart(9,m)=energy
    enddo
  else ! when irk==2, use v(t+dt/2) as the velocity
    !$omp parallel do private(m)
    do m=1,mp
      zpart0(4,m)=zpart(4,m)
      zpart0(6,m)=zpart(6,m)
      zpart0(8,m)=zpart(8,m)
    enddo
  endif
  ! update particle weight
  if(frcfullf==0)then
  !$omp parallel do private(m, energy, wdrive, wdot, t_inv, kappanp,kappatp,rdum_inv,&
  !$omp& t_cy,vecb,vec1,vec2,vecv,vecrho,psidum,ipsi,dpsi,nrsp,dpr,dr2,nzsp,dpz,dz2,&
  !$omp& gradpsi_r, gradpsi_z,br,bz,b,b_inv,ii,dp1,rdum0,zdum0,omegac)
  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 .or. psidum<psi0_sol)then
      zpart(5,m)=0.0_lk
      cycle
    endif
    nrsp=max(1,min(lsr-1,int(ceiling(zpart0(11,m)*spdR_inv))))
    dpr=zpart0(11,m)-spdR*real(nrsp-1)
    dr2=dpr*dpr
    nzsp=max(1,min(lsz-1,int(ceiling((zpart0(13,m)+zmax)*spdZ_inv))))
    dpz=zpart0(13,m)+zmax-spdZ*real(nzsp-1)
    dz2=dpz*dpz
    rdum_inv=1.0_lk/zpart0(11,m)
    ! componets of magnetic field at t
    gradpsi_r= psirz_Eq(2,nrsp,nzsp) + 2.0_lk*psirz_Eq(3,nrsp,nzsp)*dpr &
             + psirz_Eq(5,nrsp,nzsp)*dpz + 2.0_lk*psirz_Eq(6,nrsp,nzsp)*dpr*dpz &
             + psirz_Eq(8,nrsp,nzsp)*dz2 + 2.0_lk*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_lk*psirz_Eq(7,nrsp,nzsp)*dpz + 2.0_lk*psirz_Eq(8,nrsp,nzsp)*dpr*dpz + 2.0_lk*psirz_Eq(9,nrsp,nzsp)*dr2*dpz
    br=-gradpsi_z*rdum_inv
    bz=gradpsi_r*rdum_inv
    bt=0.0_lk
    b=sqrt(br**2+bz**2)
    b_inv=1.0_lk/b
    vecb=(/ br*b_inv, bt*b_inv, bz*b_inv /)
    vecv=(/ zpart0(4,m), zpart0(6,m)*zpart(11,m), zpart0(8,m)/)
    omegac=cmratio*b
    call cross_product3(vecv, vecb, vecrho)
    vecrho=vecrho/omegac
    rdum0=sqrt((zpart0(11,m)+vecrho(1))**2+vecrho(2)**2)
    zdum0=zpart0(13,m)+vecrho(3)
    zdum0=modulo(zdum0+zmax,(zmax*2.0_lk))-zmax
    psidum=spline2d(rdum0,zdum0+zmax,0,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)
    if(sfsol==0 .or. (sfcore==1 .and. psidum<psi1_core))then
      ipsi=minloc(abs(psimesh_core-psidum))
      dpsi=psidum-psimesh_core(ipsi(1))
      if(dpsi>0.0)then
        ii=min(ipsi(1),mpsi_core-1)
        psidum=min(psidum,psimesh_core(mpsi_core))
      else
        ii=max(ipsi(1)-1,1)
        psidum=max(psidum,psimesh_core(1))
      endif
      dp1=(psidum-psimesh_core(ii))/(psimesh_core(ii+1)-psimesh_core(ii)) !weight for outer flux surface
      t_cy=dp1*mesht1d(ii+1)+(1-dp1)*mesht1d(ii)
      kappanp=dp1*kapanp(ii+1)+(1-dp1)*kapanp(ii)
      kappatp=dp1*kapatp(ii+1)+(1-dp1)*kapatp(ii)
    elseif(sfcore==0 .or. (sfsol==1 .and. psidum>psi0_sol))then
      ipsi=minloc(abs(psimesh_sol-psidum))
      dpsi=psidum-psimesh_sol(ipsi(1))
      
      if(dpsi>0.0)then
        ii=min(ipsi(1),mpsi_sol-1)
        psidum=min(psidum,psimesh_sol(mpsi_sol))
      else
        ii=max(ipsi(1)-1,1)
        psidum=max(psidum,psimesh_sol(1))
      endif
      dp1=(psidum-psimesh_sol(ii))/(psimesh_sol(ii+1)-psimesh_sol(ii)) !weight for outer flux surface
      t_cy=dp1*meshti1d(ii+1+mpsi_core)+(1-dp1)*meshti1d(ii+mpsi_core)
      kappanp=dp1*kapanp(ii+1+mpsi_core)+(1-dp1)*kapanp(ii+mpsi_core)
      kappatp=dp1*kapatp(ii+1+mpsi_core)+(1-dp1)*kapatp(ii+mpsi_core)
    endif
    t_inv=1.0_lk/t_cy

    energy=zpart(9,m)
    wdrive=wpfk(2,m)*(kappanp+(energy*t_inv-1.5_lk)*kappatp)
    wdot=zpart0(7,m)*(-qpart*t_inv*(wpfk(3,m)*br*b_inv*zpart0(4,m)+wpfk(2,m)*zpart0(6,m)+wpfk(3,m)*bz*b_inv*zpart0(8,m)+zpart0(4,m)*bz*zpart(11,m)*wpfk(1,m)-zpart0(8,m)*br*zpart(11,m)*wpfk(1,m))+wdrive)
    zpart(5,m)=zpart0(5,m)+wdot*dtime
  enddo
  endif

  ! update particle position when irk==1 or irk==2
  !$omp parallel do private(m)
  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 .or. psidum<psi0_sol)then
      cycle
    endif
    zpart0(11,m)=zpart0(1,m)+dtime*zpart0(4,m)
    zpart0(12,m)=zpart0(2,m)+dtime*zpart0(6,m)
    zpart0(13,m)=zpart0(3,m)+dtime*zpart0(8,m)
    zpart(10,m)=zpart0(10,m)+dtime*zpart0(6,m)
    zpart(10,m)=modulo(zpart(10,m),pi2)
    zpart0(13,m)=modulo((zpart0(13,m)-Zsol_domain0), Zsol_domain_width)+Zsol_domain0
  enddo
  if(irk==2)then
    !$omp parallel do private(m)
    do m=1,mp
      zpart(12,m)=zpart(12,m)-(zpart0(12,m)-modulo(zpart0(12,m),torbound))
    enddo
  endif
  !$omp parallel do private(m)
  do m=1,mp
    zpart0(12,m)=modulo(zpart0(12,m),torbound)
  enddo

  !$omp parallel do private(m, energy, wdot, t_inv,kappanp,kappatp,rdum_inv,&
  !$omp& t_cy,vecb,vec1,vec2,vecv,vecrho,psidum,ipsi,dpsi,nrsp,dpr,dr2,nzsp,dpz,dz2,&
  !$omp& gradpsi_r, gradpsi_z,br,bz,b,b_inv,ii,dp1,rdum0,zdum0,omegac)
  do m=1,mp
    nrsp=max(1,min(lsr-1,int(ceiling(zpart0(11,m)*spdR_inv))))
    dpr=zpart0(11,m)-spdR*real(nrsp-1)
    dr2=dpr*dpr
    nzsp=max(1,min(lsz-1,int(ceiling((zpart0(13,m)+zmax)*spdZ_inv))))
    dpz=zpart0(13,m)+zmax-spdZ*real(nzsp-1)
    dz2=dpz*dpz
    rdum_inv=1.0_lk/zpart0(11,m)
    gradpsi_r= psirz_Eq(2,nrsp,nzsp) + 2.0_lk*psirz_Eq(3,nrsp,nzsp)*dpr &
             + psirz_Eq(5,nrsp,nzsp)*dpz + 2.0_lk*psirz_Eq(6,nrsp,nzsp)*dpr*dpz&
             + psirz_Eq(8,nrsp,nzsp)*dz2 + 2.0_lk*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_lk*psirz_Eq(7,nrsp,nzsp)*dpz +2.0_lk*psirz_Eq(8,nrsp,nzsp)*dpr*dpz + 2.0_lk*psirz_Eq(9,nrsp,nzsp)*dr2*dpz
    br=-gradpsi_z*rdum_inv
    bz=gradpsi_r*rdum_inv
    bt=0.0_lk
    b=sqrt(br**2+bz**2)
    b_inv=1.0_lk/b
    vecb=(/ br*b_inv, bt*b_inv, bz*b_inv /)
    vecv=(/ zpart0(4,m), zpart0(6,m)*zpart(11,m), zpart0(8,m)/)
    omegac=cmratio*b
    call cross_product3(vecv, vecb, vecrho)
    vecrho=vecrho/omegac
    zpart(1,m)=zpart0(11,m)
    zpart(2,m)=zpart0(12,m)
    zpart(3,m)=zpart0(13,m)
    zpart(3,m)=modulo((zpart(3,m)-Zsol_domain0), Zsol_domain_width)+Zsol_domain0
  enddo

  diagpart=0.0_lk
  vdr=0.0
  iout=0
  !$omp parallel do private(m,deltaf,angmom,energy,psidum,vdr,iout)
  do m=1,mp
    psidum=spline2d(zpart(1,m),zpart(3,m)+zmax,0,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)
    energy=0.5_lk*apart*(zpart0(4,m)**2+zpart0(6,m)**2*zpart(11,m)**2+zpart0(8,m)**2)
    deltaf=zpart(5,m)
    angmom=2*zpart0(8,m)*zpart0(11,m)*zpart0(11,m)+spline2d(zpart0(11,m),zpart0(13,m)+zmax,0,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)
    diagpart(1)=diagpart(1)+deltaf
    diagpart(2)=diagpart(2)+deltaf*deltaf
    diagpart(3)=diagpart(3)+angmom
    diagpart(4)=diagpart(4)+angmom*deltaf
    diagpart(5)=diagpart(5)+energy
    diagpart(6)=diagpart(6)+energy*deltaf
    diagpart(7)=diagpart(7)+vdr*deltaf
    diagpart(8)=diagpart(8)+vdr*angmom*deltaf
    diagpart(9)=diagpart(9)+vdr*energy*deltaf
    if(psidum>psi1_sol .or. psidum<psi0_sol) iout=iout+1
  enddo
  diagpart(10)=real(iout)
  diagpart(11)=real(mp)


  icount=mpdiag
  diagtmp=0.0
  call MPI_REDUCE(diagpart,diagtmp,icount,mpi_Rsize,MPI_SUM,0,MPI_COMM_WORLD,ierror)
  diagpart=diagtmp

  
end subroutine pushfk
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine boundary_particle

    use particle_array, only: zfki
    use equilibrium, only: psirz_Eq
    use spline_function, only: spline2d
    use equilibrium, only: lsr, lsz, spdR, spdZ, zmax
    use simulation_domain, only: psi0_core, psi1_core, psi0_sol, psi1_sol, &
                               sfcore, sfsol
    use global_parameters
    implicit none

    integer m, mout,globalmp,globalmout,ierror
    real(lk) psidum

    mout=0
    do m=1,mfki
      psidum=spline2d(zfki(1,m),zfki(3,m)+zmax,0,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)
      if(psidum>psi1_sol .or. psidum<psi0_core)then
         mout=mout+1
      elseif(sfsol==0 .and. psidum>psi1_core)then
         mout=mout+1
      elseif(sfcore==0 .and. psidum<psi0_sol)then
         mout=mout+1
      endif
    enddo
    call MPI_ALLREDUCE(mfki,globalmp,1,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,ierror)
    call MPI_ALLREDUCE(mout,globalmout,1,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,ierror)
    if(mype==0)then
      write(*,*)'total mfki:', globalmp, 'out of boundary', globalmout
    endif
end subroutine boundary_particle
