subroutine gkpush(zpart,zpart0,wppart,wtpart0,wtpart1,&
    jtpart0,jtpart1,wzpart,mesht1d, jppart,rhogyropart,rhogyropart0,&
    qpart,apart,mp,mpmax,kapanp,kapatp,diagpart, nparam, ngyro, &
    pfluxpart, rdtempart,data1dpart,ncycle, icycle, irksub, zpart1, mp1)
    use precision
    use global_parameters, only: istep, irk, shear_rate, nmode, torbound, pi, pi2, rho0, mype, idiag, gtcout
    use simulation_domain, only: Zsol_domain_width, &
                                 psi0_core, psi1_core, psi0_sol, psi1_sol, &
                                 sfcore, sfsol
    use grids_fields, only: gradphi, gradphi0
    use grids_eqdata, only: phip00
    use grids_setup
    use read_in_parameters, only: nonlinear, irestore, tstep, min_r, frcfullf, Zsol_domain0, Zsol_domain1, sol_para_periodicity, &
                                  Rdomain0, Rdomain1, ifluxtube, uni_profile, izonal, ndiag, irun
    use particle_array, only: mpdiag
    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,rspsol_psi_s,sspsol_psi_z,del_s_sol, del_z_sol, del_psi_sol,&
                           lspsol, lszsol,msp_psi_sol,zspsol_psi_s
    implicit none

!$acc routine(sppsiR1d_gpu) seq
! declaration of the dummy arguments
    integer ngyro,mp, ncycle,nparam,mpmax
    integer,dimension(:,:) :: jtpart0,jtpart1, jppart
    real(lk),dimension(:,:), optional :: zpart1
    real(lk) qpart,apart
    real(lk),dimension(:) :: wzpart
    real(lk),dimension(:) :: mesht1d,kapanp,kapatp
    real(lk),dimension(mpdiag) :: diagpart
    integer, optional :: icycle,irksub,mp1
    real(lk),dimension(:,:) :: zpart,zpart0,wppart,wtpart0,wtpart1,rhogyropart,rhogyropart0,&
                               data1dpart
    real(lk), dimension(:) :: pfluxpart, rdtempart


    integer i,j,m,igyro,ij,nrsp,nzsp,icount
  real(lk) dtime,wz1,wz0,wp1,wp0,wt11,wt10,wt01,wt00,&
       e1,e2,e3,cmratio,c_inv,wpgc(3,mpmax),rdot,zdot,zetadot,vparadot,wdot,energy0,gyrodum,perturb,&
       rdum,rdum_inv,zdum,thdum,vparadum,wdum,mudum,dpr,dr2,dpz,dz2,psidum,psidum0,&
       psifrc,gradpsi_r,gradpsi_z,grad2psi_rr,grad2psi_rz,grad2psi_zr,grad2psi_zz,&
       br,bz,b,b_inv,omegaP,omegaP_inv,&
       grad_norm_br_r,grad_norm_br_z,grad_norm_bz_r,grad_norm_bz_z,&
       grad_value_br_r,grad_value_br_z,grad_value_bz_r,grad_value_bz_z,&
       gradB0R,gradB0Z,&
       spdR_inv,spdZ_inv,psidumtmp,&
       ti_cy,ti_inv,&
       wpara,wdrive,wdrift,&
       kappanp,kappatp,kappanr,kappanz,kappatr,kappatz,kappar,kappaz,vdr,pflux_cy,&
       dedb,&
       phip00_cy,rtmp,ztmp,dpx,dx2,pdum,&
       sdum,aa,bb,cc,c2,c1
  real(lk) dbesjn,kperprho,zeta_gyro_avg,igyro_coef,igyro_coef0
  real(lk) sppsiR1d_gpu
  real(lk) zetadot_Er   !parameters for Er effect
  real(lk) wppart_gc(mp)   !guiding center weight for gyro-averaging
  integer jppart_gc(mp),left,right,np   !guiding center grid index for gyro-averaging
  real(lk) diag_b(mp), diag_psi(mp)   !used for diagnosis
  real(lk) diagtmp(mpdiag),deltaf,energy,angmom    !used for diagnosis
  integer ierror,iout    !used for diagnosis
  logical subcycle

! local parameters for pullback temperature relaxation
    integer :: ii,irsmooth
    real(lk) :: dtem(mpsi),tem_fullf(mpsi),dmark(mpsi),dden(mpsi),dden_df(mpsi),ddum(mpsi), &
                fullf,smtmp1,smtmp2,mesht1d_cy,rdtem_cy,dtest(mpsi),dtest2(mpsi),dtem_gam

#ifndef GPU_UM
  !$acc declare create(dtem, tem_fullf, dmark, dden, dden_df)
  !$acc declare create(wpgc, wppart_gc, jppart_gc)
#endif
  subcycle=.false.
  if(ncycle>0) subcycle=.true.
  wpgc=0.0_lk
  cmratio=qpart/apart
  c_inv=1.0_lk/qpart
  perturb=real(nonlinear)
  spdR_inv=1.0_lk/spdR
  spdZ_inv=1.0_lk/spdZ

  if(.not. subcycle)then
    if(irk==1)then
  ! 1st step of Runge-Kutta method
       dtime=0.5_lk*tstep
  !Notice 1:5 not 1:nparam here (zpart0(6,m) is particle weight, zpart(6,m) is magnetic moment)
#ifdef _OPENACC
  !$acc parallel loop gang vector
#else
  !$omp parallel do private(m)
#endif
       do m=1,mp
          zpart0(1:5,m)=zpart(1:5,m)
          zpart0(7,m)=zpart(7,m)
       enddo
  !$acc end parallel
  ! 2nd step of Runge-Kutta method
    else
       dtime=tstep
    endif
  else
    if(irksub==1)then
      dtime=0.25*tstep/real(ncycle)
      if(icycle==1)then
        if(irk==1)then
          mp1=mp
#ifdef _OPENACC
          !$acc parallel loop gang vector
#else
          !$omp parallel do private(m)
#endif
          do m=1,mp
            zpart1(1:nparam,m)=zpart(1:nparam,m)
            zpart1(nparam+1,m)=zpart0(6,m)
          enddo
          !$acc end parallel
        else
          mp=mp1
#ifdef _OPENACC
          !$acc parallel loop gang vector
#else
          !$omp parallel do private(m)
#endif
          do m=1,mp
            zpart(1:nparam,m)=zpart1(1:nparam,m)
            zpart0(6,m)=zpart1(nparam+1,m)
          enddo
          !$acc end parallel
        endif
      endif
#ifdef _OPENACC
      !$acc parallel loop gang vector
#else
      !$omp parallel do private(m)
#endif
      do m=1,mp
        zpart0(1:5,m)=zpart(1:5,m)
        zpart0(7,m)=zpart(7,m)
      enddo
      !$acc end parallel
    else
      dtime=0.5*tstep/real(ncycle)
    endif
  endif

! gather e_field using ngyro-point gyro-averaging
  gyrodum=1.0_lk/real(ngyro)

! electrostatic fluctuations
#ifdef _OPENACC
     !$acc parallel loop gang vector
#else
     !$omp parallel do private(m,igyro,e1,e2,e3,wz1,wz0,wp1,wp0,wt00,wt10,wt01,wt11,ij,igyro_coef)
#endif
     do m=1,mp
        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
        !$acc cache(readonly:gradphi)

        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  

           igyro_coef = rhogyropart(igyro,m) ! coefficient needed for each average point
           ij=jtpart0(igyro,m)       !lower poloidal grid on inner flux surface
           e1=e1+wt00*(wz0*gradphi(1,0,ij)+wz1*gradphi(1,1,ij))*igyro_coef
           e2=e2+wt00*(wz0*gradphi(2,0,ij)+wz1*gradphi(2,1,ij))*igyro_coef
           e3=e3+wt00*(wz0*gradphi(3,0,ij)+wz1*gradphi(3,1,ij))*igyro_coef

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

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

           ij=ij+1                  !upper poloidal grid on outer flux surface
           e1=e1+wt11*(wz0*gradphi(1,0,ij)+wz1*gradphi(1,1,ij))*igyro_coef
           e2=e2+wt11*(wz0*gradphi(2,0,ij)+wz1*gradphi(2,1,ij))*igyro_coef
           e3=e3+wt11*(wz0*gradphi(3,0,ij)+wz1*gradphi(3,1,ij))*igyro_coef
        enddo
        wpgc(1,m)=gyrodum*e1
        wpgc(2,m)=gyrodum*e2
        wpgc(3,m)=gyrodum*e3
     enddo
    if(izonal==1)then
#ifdef _OPENACC
      !$acc parallel loop gang vector
#else
      !$omp parallel do private(m,igyro,e1,e2,e3,wz1,wz0,wp1,wp0,wt00,wt10,wt01,wt11,ij,igyro_coef0)
#endif
      do m=1,mp
       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
       !$acc cache(readonly:gradphi0)
         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  
      
            igyro_coef0 = rhogyropart0(igyro,m)
            ij=jtpart0(igyro,m)       !lower poloidal grid on inner flux surface
            e1=e1+wt00*gradphi0(1,0,ij)*igyro_coef0
            e2=e2+wt00*gradphi0(2,0,ij)*igyro_coef0
            e3=e3+wt00*gradphi0(3,0,ij)*igyro_coef0
      
            ij=ij+1                  !upper poloidal grid on inner flux surface
            e1=e1+wt10*gradphi0(1,0,ij)*igyro_coef0
            e2=e2+wt10*gradphi0(2,0,ij)*igyro_coef0
            e3=e3+wt10*gradphi0(3,0,ij)*igyro_coef0
      
            ij=jtpart1(igyro,m)       !lower poloidal grid on outer flux surface
            e1=e1+wt01*gradphi0(1,0,ij)*igyro_coef0
            e2=e2+wt01*gradphi0(2,0,ij)*igyro_coef0
            e3=e3+wt01*gradphi0(3,0,ij)*igyro_coef0
      
            ij=ij+1                  !upper poloidal grid on outer flux surface
            e1=e1+wt11*gradphi0(1,0,ij)*igyro_coef0
            e2=e2+wt11*gradphi0(2,0,ij)*igyro_coef0
            e3=e3+wt11*gradphi0(3,0,ij)*igyro_coef0
         enddo
         wpgc(1,m)=gyrodum*e1+wpgc(1,m)
         wpgc(2,m)=gyrodum*e2+wpgc(2,m)
         wpgc(3,m)=gyrodum*e3+wpgc(3,m)
      enddo
  endif

! update GC position.
#ifdef _OPENACC
!$acc parallel loop gang vector
#else
!$omp parallel do private(m,rdum,rdum_inv,zdum,thdum,vparadum,wdum,mudum,nrsp,dpr,dr2,nzsp,dpz,dz2,&
!$omp& kperprho,zeta_gyro_avg,&
!$omp& psifrc,gradpsi_r,gradpsi_z,grad2psi_rr,grad2psi_rz,grad2psi_zr,grad2psi_zz,br,bz,b,b_inv,&
!$omp& grad_norm_br_r,grad_norm_br_z,grad_norm_bz_r,grad_norm_bz_z,&
!$omp& grad_value_br_r,grad_value_br_z,grad_value_bz_r,grad_value_bz_z,&
!$omp& gradB0R,gradB0Z,energy0,&
!$omp& omegaP,omegaP_inv,rdot,zetadot,zdot,vparadot,wdot,zetadot_Er,&
!$omp& ti_cy,ti_inv,wpara,wdrive,wdrift,&
!$omp& kappanp,kappatp,kappanr,kappanz,kappatr,kappatz,kappar,kappaz,vdr,dedb,phip00_cy,dpx,dx2,left,right,np,i,j)
#endif
  do m=1,mp
!~~~~~~~~~~initial particle quantities~~~~~~~~~~~~~~~~~~~~~~~
    rdum = max(zpart(1,m),min_r)
    rdum_inv=1.0_lk/rdum
    zdum = zpart(3,m)
    thdum = zpart(2,m)
    vparadum = zpart(4,m)
    wdum = zpart(5,m)
    mudum = zpart(6,m)*zpart(6,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(dabs(zdum+zmax)*spdZ_inv)))
    dpz=zdum+zmax-spdZ*real(nzsp-1)
    dz2=dpz*dpz

!!!!!!!!!!!!!!!!!!!!!!!!! psi !!!!!!!!!!!!!!!!!!!!!!!!!!
    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

    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

    grad2psi_rr= 2.0_lk*psirz_Eq(3,nrsp,nzsp)+ 2.0_lk*psirz_Eq(6,nrsp,nzsp)*dpz &
                 + 2.0_lk*psirz_Eq(9,nrsp,nzsp)*dz2
    grad2psi_rz= psirz_Eq(5,nrsp,nzsp) + 2.0_lk*psirz_Eq(6,nrsp,nzsp)*dpr &
                 + 2.0_lk*psirz_Eq(8,nrsp,nzsp)*dpz + 4.0_lk*psirz_Eq(9,nrsp,nzsp)*dpr*dpz

    grad2psi_zr= psirz_Eq(5,nrsp,nzsp) + 2.0_lk*psirz_Eq(6,nrsp,nzsp)*dpr &
                 + 2.0_lk*psirz_Eq(8,nrsp,nzsp)*dpz + 4.0_lk*psirz_Eq(9,nrsp,nzsp)*dpr*dpz
    grad2psi_zz= 2.0_lk*psirz_Eq(7,nrsp,nzsp) + 2.0_lk*psirz_Eq(8,nrsp,nzsp)*dpr &
                 + 2.0_lk*psirz_Eq(9,nrsp,nzsp)*dr2
!!!!!!!!!!!!!!!!!!!!!!!!! psi !!!!!!!!!!!!!!!!!!!!!!!!!!


    br=-gradpsi_z*rdum_inv
    bz=gradpsi_r*rdum_inv
    b=sqrt(br**2+bz**2)
    b_inv=1.0_lk/b
!unit magnetic field vector
    br=br*b_inv
    bz=bz*b_inv
!space derivative of Br and Bz
    grad_value_br_r=-grad2psi_zr*rdum_inv + gradpsi_z*rdum_inv*rdum_inv
    grad_value_br_z=-grad2psi_zz*rdum_inv
    grad_value_bz_r=grad2psi_rr*rdum_inv - gradpsi_r*rdum_inv*rdum_inv
    grad_value_bz_z=grad2psi_rz*rdum_inv
!space derivative of B0
    gradB0R=br*grad_value_br_r+bz*grad_value_bz_r
    gradB0Z=br*grad_value_br_z+bz*grad_value_bz_z
!space derivative of unit br and bz
    grad_norm_br_z=grad_value_br_z*b_inv-br*gradB0Z*b_inv
    grad_norm_bz_r=grad_value_bz_r*b_inv-bz*gradB0R*b_inv
    grad_norm_br_r=grad_value_br_r*b_inv-br*gradB0R*b_inv
    grad_norm_bz_z=grad_value_bz_z*b_inv-bz*gradB0Z*b_inv

!cyclotron frequency
    omegaP = cmratio*b
    omegaP_inv = 1.0_lk/omegaP
    energy0 = 0.5_lk*apart*vparadum*vparadum + mudum*b


    !kperprho=(zpart(6,m)*sqrt(2.0_lk*apart/b))*(real(nmode)*rdum_inv)/qpart
    !zeta_gyro_avg=dbesjn(0,kperprho)
    !wpgc(1,m)=wpgc(1,m)*zeta_gyro_avg
    !wpgc(2,m)=wpgc(2,m)*zeta_gyro_avg
    !wpgc(3,m)=wpgc(3,m)*zeta_gyro_avg

!=============================================================================

    if (ngyro==1)then
        wppart_gc(m)=wppart(1,m)
        jppart_gc(m)=jppart(1,m)
    else
        if(sfsol==1)then
            if((psifrc .ge. psi_sep) .or. sfcore==0)then     
                if(psifrc .ge. psi1_sol)then
                    wppart_gc(m)=1.0_lk
                    np=mpsi_sol-1
                    psifrc=psi1_sol
                elseif(psifrc .le. psi0_sol .and. sfcore==0)then
                    wppart_gc(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_gc(m) = (psifrc-psimesh_sol(left))/(psimesh_sol(left+1)-psimesh_sol(left))
                endif
                jppart_gc(m) = np + mpsi_core
            endif
        endif

!TODO: algorithm for np not valid in core region
        if(sfcore==1)then
            if((psifrc .le. psi1_core) .or. (sfsol==0))then
                if((psifrc .ge. psi1_core) .and. sfsol==0)then
                    wppart_gc(m)=1.0_lk
                    np=mpsi_core-1
                    psifrc=psi1_core
                elseif(psifrc .le. psi0_core)then
                    wppart_gc(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 = left
                    wppart_gc(m) = (psifrc-psimesh_core(left))/(psimesh_core(left+1)-psimesh_core(left))
                endif
                jppart_gc(m)=np
            endif
        endif

    endif
    ! ti_cy=wppart(1,m)*mesht1d(jppart(1,m)+1)+(1.0_lk-wppart(1,m))*mesht1d(jppart(1,m))
    ti_cy=wppart_gc(m)*mesht1d(jppart_gc(m)+1)+(1.0_lk-wppart_gc(m))*mesht1d(jppart_gc(m))
    ti_inv=1.0_lk/ti_cy


    ! kappanp=wppart(1,m)*kapanp(jppart(1,m)+1)+(1.0_lk-wppart(1,m))*kapanp(jppart(1,m))
    ! kappatp=wppart(1,m)*kapatp(jppart(1,m)+1)+(1.0_lk-wppart(1,m))*kapatp(jppart(1,m))
    kappanp=wppart_gc(m)*kapanp(jppart_gc(m)+1)+(1.0_lk-wppart_gc(m))*kapanp(jppart_gc(m))
    kappatp=wppart_gc(m)*kapatp(jppart_gc(m)+1)+(1.0_lk-wppart_gc(m))*kapatp(jppart_gc(m))

    kappanr=kappanp*gradpsi_r
    kappanz=kappanp*gradpsi_z
    kappatr=kappatp*gradpsi_r
    kappatz=kappatp*gradpsi_z


    kappar=kappanr+kappatr*(energy0*ti_inv-1.5_lk)
    kappaz=kappanz+kappatz*(energy0*ti_inv-1.5_lk)

    vdr=-wpgc(2,m)*b_inv*rdum_inv
    dedb=apart*vparadum*vparadum*b_inv+mudum

    wdrive=vdr*(br*kappaz-bz*kappar)
    wpara=-qpart*vparadum*ti_inv*wpgc(3,m)

  !  wdrift=-qpart*ti_inv*c_inv*dedb*(bz*gradB0R-br*gradB0Z)*b_inv*wpgc(2,m)*rdum_inv

    wdrift=-qpart*ti_inv*rdum_inv*b_inv*wpgc(2,m)*(c_inv*dedb*(bz*gradB0R - br*gradB0Z)&
            + vparadum*vparadum*omegap_inv*(grad_value_br_z - grad_value_bz_r))

    if(uni_profile .ne. 1)then
        wdot=(zpart0(6,m)-1.0_lk*perturb*zpart(5,m))*(1.0_lk*wpara + 1.0_lk*wdrive + 1.0_lk*wdrift)
    else
        wdot=(zpart0(6,m)-0.0_lk*zpart(5,m))*(1.0_lk*wpara + 1.0_lk*wdrive + 1.0_lk*wdrift)
    endif

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!! Er effect implementation
!   Above is the induced angular velocity in zeta, w.r.t. c_s
!   In gtc.in, the unit for tstep has the relation:
!       tstep_unit = R_0/c_s = sqrt(apart)/rho0/omega_0
!   where c_s is the sound speed, omega_0 is proton cyclotron freq. 
!   In the code, utime = 1.0/omegacp, thus 
!       tstep_unit_code = tstep_unit*omega_0 = sqrt(apart)/rho0
!   i.e. tstep_unit is rewritten in unit of utime, and the frequency has a
!   reference value of 1.0/tstep_unit.
    phip00_cy=wppart_gc(m)*phip00(jppart_gc(m)+1)+(1.0_lk-wppart_gc(m))*phip00(jppart_gc(m))
    zetadot_Er = phip00_cy*b_inv*(bz*gradpsi_r-br*gradpsi_z)*rdum_inv

! the averaged flux in \psi direction, to be subtracted in nonlinear run.
    pflux_cy=wppart_gc(m)*pfluxpart(jppart_gc(m)+1)+(1.0_lk-wppart_gc(m))*pfluxpart(jppart_gc(m))

    rdot=vparadum*br + 1.0_lk*perturb*(vdr-pflux_cy)*bz
    zetadot=((bz*gradB0R-br*gradB0Z)*mudum/(omegaP*apart)& !gradB drift
            +(grad_norm_br_z-grad_norm_bz_r)*vparadum*vparadum/omegaP)*rdum_inv& !curvature B drift
            + zetadot_Er &   !Er drift
            +1.0_lk*perturb*wpgc(1,m)

    zdot=vparadum*bz - 1.0_lk*perturb*(vdr-pflux_cy)*br
    vparadot=-(br*gradB0R+bz*gradB0Z)*mudum/apart-1.0_lk*perturb*cmratio*wpgc(3,m)   &
             -1.0_lk*perturb*cmratio*vparadum*rdum_inv/omegaP*wpgc(2,m)*(bz*gradB0R-br*gradB0Z+grad_value_br_z-grad_value_bz_r)


    zpart(1,m) = zpart0(1,m)+dtime*rdot
    zpart(3,m) = zpart0(3,m)+dtime*zdot
    zpart(4,m) = zpart0(4,m)+dtime*vparadot
    if(frcfullf==0)zpart(5,m)=zpart0(5,m)+dtime*wdot
   
    zpart(2,m) = zpart0(2,m)+dtime*zetadot
    zpart(2,m) = modulo(zpart(2,m),torbound)
 
    if(zpart(1,m)<0.0_lk)zpart(1,m)=-zpart(1,m)
    wpgc(1,m)=energy0
    wpgc(2,m)=2*zetadot*zpart0(1,m)*zpart0(1,m)
    wpgc(3,m)=vdr
    rdum=max(0.0_lk,zpart(1,m))
    zdum=max(0.0_lk,zpart(3,m)+zmax)
    i = max(1, min(lsr-1, ceiling(rdum/spdR)))
    j = max(1, min(lsz-1, ceiling(zdum/spdZ)))
    dpr = rdum - spdR*real(i-1, lk)
    dpz = zdum - spdZ*real(j-1, lk)
    dr2 = dpr*dpr
    dz2 = dpz*dpz
  
    zpart(7,m) = psirz_Eq(1,i,j) + psirz_Eq(2,i,j)*dpr + psirz_Eq(3,i,j)*dr2 +&
                (psirz_Eq(4,i,j) + psirz_Eq(5,i,j)*dpr + psirz_Eq(6,i,j)*dr2)*dpz+&
                (psirz_Eq(7,i,j) + psirz_Eq(8,i,j)*dpr + psirz_Eq(9,i,j)*dr2)*dz2

  enddo
  !$acc end parallel

if((subcycle .and. irksub==2) .or.((.not. subcycle) .and. irk==2))then

! out of boundary particle
#ifdef _OPENACC
!$acc parallel loop gang vector
#else
!$omp parallel do private(m,psidum0,psidum,psidumtmp,c1,c2,aa,bb,cc,ii,zdum,pdum,&
!$omp& dpx,dx2,dpz,dz2,i,j,sdum,rdum)
#endif
    do m=1,mp
      psidum0 = zpart(7,m)

      if(sfcore==1 .and. sfsol==0)then

!DONT CHANGE ANYTHING(repeat part should be in ifelse control)
         if(psidum0 < psi0_core)then
             zpart(1,m)=zpart0(1,m)
             zpart(3,m)=zpart0(3,m)
             zpart(4,m)=-zpart0(4,m)   ! negative sign for offsetting the numerical error
             zpart(5,m)=0.0_lk
         elseif(psidum0 > psi1_core)then
             zpart(1,m)=zpart0(1,m)
             zpart(3,m)=zpart0(3,m)
             zpart(4,m)=-zpart0(4,m)   ! negative sign for offsetting the numerical error
             zpart(5,m)=0.0_lk
         endif

       elseif(sfcore==0 .and. sfsol==1)then
!DONT CHANGE ANYTHING(repeat part should be in ifelse control)
!         if(psidum0 > psi1_sol)then
!           sdum=spline2d(psidum0-msp_psi_sol(1),zpart(3,m)-Zsol_domain0,0,lspsol,lszsol,del_psi_sol,del_z_sol,sspsol_psi_z,0,0)
!           if(ifluxtube==0)then
!             psidum = 2.0_lk * psi1_sol - psidum0
!           else
!             psidum = psi0_sol + psidum0-psi1_sol
!           endif
!           zpart(1,m)=spline2d(psidum-msp_psi_sol(1),sdum,0,lspsol,lszsol,del_psi_sol,del_s_sol,rspsol_psi_s,0,0)
!           zpart(3,m)=spline2d(psidum-msp_psi_sol(1),sdum,0,lspsol,lszsol,del_psi_sol,del_s_sol,zspsol_psi_s,0,0)
!         elseif(psidum0 < psi0_sol)then
!           sdum=spline2d(psidum0-msp_psi_sol(1),zpart(3,m)-Zsol_domain0,0,lspsol,lszsol,del_psi_sol,del_z_sol,sspsol_psi_z,0,0)
!           if(ifluxtube==0)then
!             psidum = 2.0_lk * psi0_sol - psidum0
!           else
!             psidum = psi1_sol - (psi0_sol-psidum0)
!           endif
!           zpart(1,m)=spline2d(psidum-msp_psi_sol(1),sdum,0,lspsol,lszsol,del_psi_sol,del_s_sol,rspsol_psi_s,0,0)
!           zpart(3,m)=spline2d(psidum-msp_psi_sol(1),sdum,0,lspsol,lszsol,del_psi_sol,del_s_sol,zspsol_psi_s,0,0)
!         endif

         if(psidum0 > psi1_sol)then
           if(ifluxtube==1)then
             psidum = psi0_sol + psidum0 - psi1_sol
           else
             psidum=2.0_lk*psi1_sol-psidum0
           endif
           zdum=zpart0(3,m)+zmax
           j=max(1,min(lsz-1,ceiling(zdum/spdZ)))
           dpz=zdum-spdZ*real(j-1)
           dz2=dpz*dpz


           do i=1,lsr-1
              c1=psirz_Eq(1,i,j)+psirz_Eq(4,i,j)*dpz+psirz_Eq(7,i,j)*dz2
              c2=psirz_Eq(1,i+1,j)+psirz_Eq(4,i+1,j)*dpz+psirz_Eq(7,i+1,j)*dz2
   
              if((psidum-c1)*(psidum-c2) .lt. 0.0_lk)then
               ii=i
               aa=psirz_Eq(3,ii,j)+psirz_Eq(6,ii,j)*dpz+psirz_Eq(9,ii,j)*dz2
               bb=psirz_Eq(2,ii,j)+psirz_Eq(5,ii,j)*dpz+psirz_Eq(8,ii,j)*dz2
               cc=psirz_Eq(1,ii,j)+psirz_Eq(4,ii,j)*dpz+psirz_Eq(7,ii,j)*dz2
               exit
              endif
           enddo
           zpart(1,m)=sppsiR1d_gpu(aa,bb,cc,spdR,zdum,c2) + (ii-1)*spdR
           zpart(3,m)=zpart0(3,m)
           zpart(5,m)=0.0_lk
         elseif(psidum0 < psi0_sol)then
           if(ifluxtube==1)then
             psidum = psi1_sol - (psi0_sol-psidum0)
           else
             psidum=2.0_lk*psi0_sol-psidum0
           endif
           zdum=zpart0(3,m)+zmax
           j=max(1,min(lsz-1,ceiling(zdum/spdZ)))
           dpz=zdum-spdZ*real(j-1)
           dz2=dpz*dpz


           do i=1,lsr-1
              c1=psirz_Eq(1,i,j)+psirz_Eq(4,i,j)*dpz+psirz_Eq(7,i,j)*dz2
              c2=psirz_Eq(1,i+1,j)+psirz_Eq(4,i+1,j)*dpz+psirz_Eq(7,i+1,j)*dz2

              if((psidum-c1)*(psidum-c2) .lt. 0.0_lk)then
               ii=i
               aa=psirz_Eq(3,ii,j)+psirz_Eq(6,ii,j)*dpz+psirz_Eq(9,ii,j)*dz2
               bb=psirz_Eq(2,ii,j)+psirz_Eq(5,ii,j)*dpz+psirz_Eq(8,ii,j)*dz2
               cc=psirz_Eq(1,ii,j)+psirz_Eq(4,ii,j)*dpz+psirz_Eq(7,ii,j)*dz2
               exit
              endif
           enddo
           zpart(1,m)=sppsiR1d_gpu(aa,bb,cc,spdR,zdum,c2) + (ii-1)*spdR

           zpart(3,m)=zpart0(3,m)
           zpart(5,m)=0.0_lk
         endif

       elseif(sfcore==1 .and. sfsol==1)then

         if(psidum0 .ge. psi_sep)then
           if(psidum0 > psi1_sol)then
             psidum=2.0_lk*psi1_sol-psidum0
             zdum=zpart0(3,m)+zmax
             j=max(1,min(lsz-1,ceiling(zdum/spdZ)))
             dpz=zdum-spdZ*real(j-1)
             dz2=dpz*dpz
  
  
             do i=1,lsr-1
                c1=psirz_Eq(1,i,j)+psirz_Eq(4,i,j)*dpz+psirz_Eq(7,i,j)*dz2
                c2=psirz_Eq(1,i+1,j)+psirz_Eq(4,i+1,j)*dpz+psirz_Eq(7,i+1,j)*dz2
  
                if((psidum-c1)*(psidum-c2) .lt. 0.0_lk)then
                 ii=i
                 aa=psirz_Eq(3,ii,j)+psirz_Eq(6,ii,j)*dpz+psirz_Eq(9,ii,j)*dz2
                 bb=psirz_Eq(2,ii,j)+psirz_Eq(5,ii,j)*dpz+psirz_Eq(8,ii,j)*dz2
                 cc=psirz_Eq(1,ii,j)+psirz_Eq(4,ii,j)*dpz+psirz_Eq(7,ii,j)*dz2
                 exit
                endif
             enddo
             zpart(1,m)=sppsiR1d_gpu(aa,bb,cc,spdR,zdum,c2) + (ii-1)*spdR

             zpart(3,m)=zpart0(3,m)
             zpart(5,m)=0.0_lk
           endif
         else
           if(psidum0<psi0_core)then
             zpart(1,m)=zpart0(1,m)
             zpart(3,m)=zpart0(3,m)
             zpart(4,m)=-zpart0(4,m)   ! negative sign for offsetting the numerical error
             zpart(5,m)=0.0_lk
           endif
         endif
       endif

       if(sfsol==1)then   
           if(sol_para_periodicity==1)then   ! periodic particle boundary condition
               if(zpart(3,m) > Zsol_domain1)then
                   psidum=zpart0(7,m)
                   zpart(3,m)=modulo(zpart(3,m)-Zsol_domain_width,-Zsol_domain_width)
                   pdum=max(0.0_lk,zpart0(7,m)-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
                   sdum = 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
                   sdum=max(0.0_lk,sdum)
                   j = max(1, min(lszsol-1, ceiling(sdum/del_s_sol)))
                   dpz = sdum - del_s_sol*real(j-1, lk)
                   dz2 = dpz*dpz
                   zpart(1,m) = rspsol_psi_s(1,i,j) + rspsol_psi_s(2,i,j)*dpx + rspsol_psi_s(3,i,j)*dx2 +&
                               (rspsol_psi_s(4,i,j) + rspsol_psi_s(5,i,j)*dpx + rspsol_psi_s(6,i,j)*dx2)*dpz+&
                               (rspsol_psi_s(7,i,j) + rspsol_psi_s(8,i,j)*dpx + rspsol_psi_s(9,i,j)*dx2)*dz2
               elseif(zpart(3,m) < Zsol_domain0)then
                   psidum=zpart0(7,m)
                   zpart(3,m)=modulo(zpart(3,m)+Zsol_domain_width,Zsol_domain_width)
                   pdum=max(0.0_lk,zpart0(7,m)-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
                   sdum = 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
                   sdum=max(0.0_lk,sdum)
                   j = max(1, min(lszsol-1, ceiling(sdum/del_s_sol)))
                   dpz = sdum - del_s_sol*real(j-1, lk)
                   dz2 = dpz*dpz
                   zpart(1,m) = rspsol_psi_s(1,i,j) + rspsol_psi_s(2,i,j)*dpx + rspsol_psi_s(3,i,j)*dx2 +&
                               (rspsol_psi_s(4,i,j) + rspsol_psi_s(5,i,j)*dpx + rspsol_psi_s(6,i,j)*dx2)*dpz+&
                               (rspsol_psi_s(7,i,j) + rspsol_psi_s(8,i,j)*dpx + rspsol_psi_s(9,i,j)*dx2)*dz2

                endif

           elseif(sol_para_periodicity==0)then    ! reflected particle boundary condition

               if(zpart(3,m) > Zsol_domain1)then
                   zpart(1,m)=zpart0(1,m)
                   zpart(2,m)=zpart0(2,m)
                   zpart(3,m)=zpart0(3,m)
                   zpart(4,m)=-zpart0(4,m)
                   zpart(5,m)=zpart0(5,m)
               elseif(zpart(3,m) < Zsol_domain0)then
                   zpart(1,m)=zpart0(1,m)
                   zpart(2,m)=zpart0(2,m)
                   zpart(3,m)=zpart0(3,m)
                   zpart(4,m)=-zpart0(4,m)
                   zpart(5,m)=zpart0(5,m)
               endif

           endif

       endif
       rdum=max(0.0_lk,zpart(1,m))
       zdum=max(0.0_lk,zpart(3,m)+zmax)
       i = max(1, min(lsr-1, ceiling(rdum/spdR)))
       j = max(1, min(lsz-1, ceiling(zdum/spdZ)))
       dpr = rdum - spdR*real(i-1, lk)
       dpz = zdum - spdZ*real(j-1, lk)
       dr2 = dpr*dpr
       dz2 = dpz*dpz
   
       zpart(7,m) = psirz_Eq(1,i,j) + psirz_Eq(2,i,j)*dpr + psirz_Eq(3,i,j)*dr2 +&
                   (psirz_Eq(4,i,j) + psirz_Eq(5,i,j)*dpr + psirz_Eq(6,i,j)*dr2)*dpz+&
                   (psirz_Eq(7,i,j) + psirz_Eq(8,i,j)*dpr + psirz_Eq(9,i,j)*dr2)*dz2
    enddo
    !$acc end parallel
endif

if (idiag==0 .or. istep+irun+irk==2) then
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! implement temperature restore and output particle diagnosis when irk==1, i.e. idiag=0.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!! pull back the temperature relaxation during the simulation
! initialization of pullback arrays
    !$acc parallel loop
    do i=1,mpsi
        dtem(i) = 0.0_lk
        dmark(i) = 0.0_lk
        dden(i) = 0.0_lk
        dden_df(i) = 0.0_lk
        tem_fullf(i) = 0.0_lk
    enddo
    !$acc end parallel

! gather particle contribution onto 1d psi arrays
#ifdef _OPENACC
    !$acc parallel loop
#else
!TODO: Known issues, if one puts tem_fullf to the 3rd place, it will be all 0. Why?
    !$omp parallel do private(m,ii,fullf,deltaf) reduction(+:tem_fullf,dtem,dmark,dden,dden_df)
#endif
    do m=1,mp
        ii = jppart_gc(m)
        fullf = zpart0(6,m)
        deltaf = zpart(5,m)

        !$acc atomic update
        dtem(ii) = dtem(ii)+wpgc(1,m)*deltaf*(1.0_lk-wppart_gc(m))
        !$acc atomic update
        dmark(ii) = dmark(ii)+wpgc(3,m)*fullf*(1.0_lk-wppart_gc(m))
        !$acc atomic update
        tem_fullf(ii) = tem_fullf(ii)+wpgc(1,m)*fullf*(1.0_lk-wppart_gc(m))
        !$acc atomic update
        dden(ii) = dden(ii)+fullf*(1.0_lk-wppart_gc(m))
        !$acc atomic update
        dden_df(ii) = dden_df(ii)+deltaf*(1.0_lk-wppart_gc(m))

        !$acc atomic update
        dtem(ii+1) = dtem(ii+1)+wpgc(1,m)*deltaf*wppart_gc(m)
        !$acc atomic update
        dmark(ii+1) = dmark(ii+1)+wpgc(3,m)*fullf*wppart_gc(m)
        !$acc atomic update
        tem_fullf(ii+1) = tem_fullf(ii+1)+wpgc(1,m)*fullf*wppart_gc(m)
        !$acc atomic update
        dden(ii+1) = dden(ii+1)+fullf*wppart_gc(m)
        !$acc atomic update
        dden_df(ii+1) = dden_df(ii+1)+deltaf*wppart_gc(m)
    enddo
    !$acc end parallel
#ifndef GPU_UM
    !$acc update host(dden, dmark, tem_fullf, dden_df, dtem)
#endif

! sum over all mpi to get particle contribution from all domains
    icount = mpsi
    ddum = 0.0_lk
    call MPI_REDUCE(dtem,ddum,icount,mpi_Rsize,MPI_SUM,0,MPI_COMM_WORLD,ierror)
    dtem = ddum
    call MPI_REDUCE(dmark,ddum,icount,mpi_Rsize,MPI_SUM,0,MPI_COMM_WORLD,ierror)
    dmark = ddum
    call MPI_REDUCE(dden,ddum,icount,mpi_Rsize,MPI_SUM,0,MPI_COMM_WORLD,ierror)
    dden = ddum
    call MPI_REDUCE(dden_df,ddum,icount,mpi_Rsize,MPI_SUM,0,MPI_COMM_WORLD,ierror)
    dden_df = ddum
    call MPI_REDUCE(tem_fullf,ddum,icount,mpi_Rsize,MPI_SUM,0,MPI_COMM_WORLD,ierror)
    tem_fullf = ddum

    tem_fullf=tem_fullf/max(1.0_lk,dden)/1.5_lk !full temperature, at each 1d \psi location
    dmark=dmark/max(1.0_lk,dden) !radial marker flux, vdr averaged at each 1d \psi location
    dtem=dtem/(mesht1d*max(1.0_lk,dden))/1.5_lk !dT/T averaged at each 1d \psi location
     
!!! output dtem and dn information into data1d array.
    data1dpart=0.0_lk
    data1dpart(:,6)=tem_fullf(:)
    data1dpart(:,7)=dtem(:)*mesht1d(:)
    data1dpart(:,8)=dden(:)
    data1dpart(:,9)=dden_df(:)

    pfluxpart = 0.0_lk
    rdtempart = 0.0_lk
    if(irestore==1) then

        !!! smooth out turbulence scale variance of particle flux
            irsmooth=0!int(sqrt(real(mpsi)))
            do i=1,irsmooth
                smtmp1 = dmark(2)
                smtmp2 = dmark(mpsi-1)
                dmark(2:mpsi-1)=0.5*dmark(2:mpsi-1)+0.25*(dmark(1:mpsi-2)+dmark(3:mpsi))
                dmark(1)=0.5*(dmark(1)+smtmp1)
                dmark(mpsi)=0.5*(dmark(mpsi)+smtmp2)
            enddo
            icount=mpsi
            call MPI_BCAST(dmark,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)   

            smtmp1 = 0.08_lk   !0.1_lk
            pfluxpart = (1.0_lk-smtmp1)*pfluxpart + smtmp1*dmark
!!! comment out the following code if you want to use the flux restoration
            pfluxpart = 0.0_lk

        !!! smooth out turbulence scale variance of temperature flux
            irsmooth=int(mpsi)  !int(mpsi/2)
            do i=1,irsmooth
                smtmp1 = dtem(2)
                smtmp2 = dtem(mpsi-1)
                dtem(2:mpsi-1)=0.5*dtem(2:mpsi-1)+0.25*(dtem(1:mpsi-2)+dtem(3:mpsi))
                dtem(1)=0.5*(dtem(1)+smtmp1)
                dtem(mpsi)=0.5*(dtem(mpsi)+smtmp2)
            enddo
            icount=mpsi
            call MPI_BCAST(dtem,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)   

            ! dtem_gam is currently chosen as the linear growth rate of default 
            ! constructed profile. After 1/(tstep*dtem_gam) steps, we completely 
            ! restore dtem. This assures the physical scale of T restoration 
            ! won't change when changing tstep.
            ! By replacing tstep by 0.22_lk, dtem is fully restored in 1291 steps.
            dtem_gam = 0.00352_lk 
            smtmp1 = min(1.0_lk*(ndiag*0.22_lk*dtem_gam),1.0_lk)   !0.01_lk
            if (smtmp1>0.2_lk)then
                write(gtcout,*)'The time averaging for T restoration is less than 5 tstep, &
                    meaning that you have less than 5 steps for dtem_gam in gkpush.F90, &
                    please use an appropriate time scale for dtem average, & 
                    or a smaller time step, or turn off irestore.'
                call flush(gtcout)
                call GTC_ABORT
            endif
            rdtempart = (1.0_lk-smtmp1)*rdtempart + smtmp1*dtem

        !!! smooth the temp. obtained from fullf before subtracting it from particle velocity.
            irsmooth=int(mpsi)  !int(mpsi/2)
            tem_fullf(:) = tem_fullf(:)-mesht1d(:)
            do i=1,irsmooth
                smtmp1 = tem_fullf(2)
                smtmp2 = tem_fullf(mpsi-1)
                tem_fullf(2:mpsi-1)=0.5*tem_fullf(2:mpsi-1)+0.25*(tem_fullf(1:mpsi-2)+tem_fullf(3:mpsi))
                tem_fullf(1)=0.5*(tem_fullf(1)+smtmp1)
                tem_fullf(mpsi)=0.5*(tem_fullf(mpsi)+smtmp2)
            enddo
            tem_fullf(:) = mesht1d(:)+tem_fullf(:)

            icount=mpsi
            call MPI_BCAST(tem_fullf,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
#ifndef GPU_UM 
            !$acc update device(pfluxpart, rdtempart, tem_fullf)
#endif

        !!! subtract the deltaf and fullf changes due to the temperature relaxation
#ifdef _OPENACC
            !$acc parallel loop gang vector
#else
            !$omp parallel do private(m,ii,mesht1d_cy,rdtem_cy,fullf)
#endif
            do m=1,mp
                ii = jppart_gc(m)
                mesht1d_cy = wppart_gc(m)*mesht1d(jppart_gc(m)+1)+(1.0_lk-wppart_gc(m))*mesht1d(jppart_gc(m))
                rdtem_cy = wppart_gc(m)*rdtempart(jppart_gc(m)+1)+(1.0_lk-wppart_gc(m))*rdtempart(jppart_gc(m))
                
                fullf = zpart0(6,m)
                zpart(5,m) = zpart(5,m)-(wpgc(1,m)/mesht1d_cy-1.5)*rdtem_cy*fullf

        !!! modify the velocity so the temperature dependence is still (v_0^2/T_0),
        !!! by scaling v_s=v*sqrt(T_0/T), so (v^2/T)=(v_s^2/T_0)
                rdtem_cy = wppart_gc(m)*tem_fullf(jppart_gc(m)+1)+(1.0_lk-wppart_gc(m))*tem_fullf(jppart_gc(m))
                rdtem_cy = mesht1d_cy/rdtem_cy
                zpart(4,m) = zpart(4,m)*sqrt(rdtem_cy)
                zpart(6,m) = zpart(6,m)*sqrt(rdtem_cy)
            enddo

        !!! pflux is used in gkpush() function, when pushing particles, the averaged vdr
        !!! is subtracted from the push velocity of the particle.
    endif


!!! compute the particle diagnosis for volume-averaged quantities.
    diagpart=0.0_lk
    vdr=0.0
    iout=0
#ifndef GPU_UM
    !$acc update device(diagpart,data1dpart)
#endif
#ifdef _OPENACC
    !$acc parallel loop gang vector reduction(+:iout)
#else
    !$omp parallel do private(m,deltaf,fullf,angmom,energy,psidum,vdr) reduction(+:diagpart,data1dpart,iout)
#endif
    do m=1,mp
        psidum=zpart(7,m)
        
        mesht1d_cy = wppart_gc(m)*mesht1d(jppart_gc(m)+1)+(1.0_lk-wppart_gc(m))*mesht1d(jppart_gc(m))
        energy=wpgc(1,m)/mesht1d_cy - 1.5_lk
        deltaf=zpart(5,m)
        fullf=zpart0(6,m)
        angmom=psidum
        vdr=wpgc(3,m)
        !$acc atomic update
        diagpart(1)=diagpart(1)+deltaf
        !$acc atomic update
        diagpart(2)=diagpart(2)+deltaf*deltaf
        !$acc atomic update
        diagpart(3)=diagpart(3)+angmom*fullf
        !$acc atomic update
        diagpart(4)=diagpart(4)+angmom*deltaf
        !$acc atomic update
        diagpart(5)=diagpart(5)+energy*fullf
        !$acc atomic update
        diagpart(6)=diagpart(6)+energy*deltaf
        !$acc atomic update
        diagpart(7)=diagpart(7)+vdr*deltaf
        !$acc atomic update
        diagpart(8)=diagpart(8)+vdr*angmom*deltaf
        !$acc atomic update
        diagpart(9)=diagpart(9)+vdr*energy*deltaf
        if(psidum>psi1_sol .or. psidum<psi0_sol) iout=iout+1
  
        ii = jppart_gc(m)
        !$acc atomic update
        data1dpart(ii,1) = data1dpart(ii,1)+vdr*deltaf*(1.0_lk-wppart_gc(m))
        !$acc atomic update
        data1dpart(ii+1,1) = data1dpart(ii+1,1)+vdr*deltaf*wppart_gc(m)
        !$acc atomic update
        data1dpart(ii,2) = data1dpart(ii,2)+vdr*angmom*deltaf*(1.0_lk-wppart_gc(m))
        !$acc atomic update
        data1dpart(ii+1,2) = data1dpart(ii+1,2)+vdr*angmom*deltaf*wppart_gc(m)
        !$acc atomic update
        data1dpart(ii,3) = data1dpart(ii,3)+vdr*energy*deltaf*(1.0_lk-wppart_gc(m))
        !$acc atomic update
        data1dpart(ii+1,3) = data1dpart(ii+1,3)+vdr*energy*deltaf*wppart_gc(m)
    enddo

#ifndef GPU_UM
    !$acc update host(diagpart(1:9), data1dpart(:,1:3))
#endif
    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

! sum over all mpi to get particle contribution from all domains
    icount = mpsi
    call MPI_REDUCE(data1dpart(:,1),ddum,icount,mpi_Rsize,MPI_SUM,0,MPI_COMM_WORLD,ierror)
    data1dpart(:,1) = ddum/max(1.0_lk,dden)
    call MPI_REDUCE(data1dpart(:,2),ddum,icount,mpi_Rsize,MPI_SUM,0,MPI_COMM_WORLD,ierror)
    data1dpart(:,2) = ddum/max(1.0_lk,dden)
    call MPI_REDUCE(data1dpart(:,3),ddum,icount,mpi_Rsize,MPI_SUM,0,MPI_COMM_WORLD,ierror)
    data1dpart(:,3) = ddum/max(1.0_lk,dden)

!    data1dpart(:,4) = tem_fullf(:)
!    data1dpart(:,5) = dtem(:)*mesht1d(:)
!    data1dpart(:,6) = dden(:)
!    data1dpart(:,7) = dden_df(:)

endif
end subroutine gkpush
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#ifdef DOUBLE_PRECISION
real(8) function sppsiR1d_gpu(a,b,c,dx,zi,c2)
#else
real(4) function sppsiR1d_gpu(a,b,c,dx,zi,c2)
#endif
!$acc routine seq
use precision, only: lk
implicit none
real(lk):: zi,a,b,c,dx,c2

real(lk):: phigh, plow,dpx,dx2,zdum,sperr,del_zdum

integer :: s
            phigh=dx
            plow=0.0_lk
            dpx=0.5_lk*(phigh+plow)
            dx2=dpx*dpx
            zdum=a*dx2+b*dpx+c
            sperr=0.000001_lk
            del_zdum=(zi-zdum)/dabs(zi)  ! valid for both negative and positive zi
            s=1
            do while (dabs(del_zdum)>sperr .and. s<1000)
                s=s+1
                if(zi>c2)then      ! case for monotonic decreasing
                    if(del_zdum<0)then
                        plow=dpx
                    else
                        phigh=dpx
                    endif
                elseif(zi<c2)then  ! case for monotonic increasing
                    if(del_zdum>0)then
                        plow=dpx
                    else
                        phigh=dpx
                    endif
                endif

                dpx=0.5_lk*(plow+phigh)
                dx2=dpx*dpx
                zdum=a*dx2+b*dpx+c
                del_zdum=(zi-zdum)/dabs(zi)  ! valid for both negative and positive zi
            enddo

!Recalculate the Z value by using the x and y position
            sppsiR1d_gpu=dpx

end function sppsiR1d_gpu
