! This file is part of GTC version 4.4 
! GTC version 4.4 is released under the 3-Clause BSD license:

! Copyright (c) 2002,2010,2016, GTC Team (team leader: Zhihong Lin, zhihongl@uci.edu)
! All rights reserved.

! Redistribution and use in source and binary forms, with or without 
! modification, are permitted provided that the following conditions are met:

! 1. Redistributions of source code must retain the above copyright notice, 
!    this list of conditions and the following disclaimer.

! 2. Redistributions in binary form must reproduce the above copyright notice, 
!    this list of conditions and the following disclaimer in the documentation 
!    and/or other materials provided with the distribution.

! 3. Neither the name of the GTC Team nor the names of its contributors may be 
!    used to endorse or promote products derived from this software without 
!    specific prior written permission.
! ==============================================================================

#ifdef _hybrid
#define pushParticle hybridPushParticle
#else
#define pushParticle gkPushParticle
#endif

!> \brief Unified particle push subroutine.
!
!> This subroutine is called through the wrapper subroutine \c push in file \c push.F90.\n
!
!> The particle information array \c zpart is updated based on the particle
!! species, Runge-Kutta step, push type ( \c hybrid or \c gyrokinetic ),
!> subcycling setting, and plasma field evaluated at each particle location.\n
!>\n 
!>
!> Essentially, 5 particle quantities are pushed:
!>  - psi ( \c zpart(1) ): poloidal flux, effective radial coordinate
!>  - theta ( \c zpart(2) ): poloidal angle
!>  - zeta ( \c zpart(3) ): toroidal angle
!>  - rho_para ( \c zpart(4) ): v_para/omega_c, effective parallel gyroradius
!>  - weight ( \c zpart(5) ): distribution function carried by the marker
!>  .
!> \n
!>
!> A push consists of the following parts:
!>  -# set time step based on specific push type
!>    - First ion RK ( \c irk==1 ): half ion step
!>    - Second ion RK ( \c irk==2 ): full ion step
!>    - First electron RK ( \c irksub==1 ): half electron step
!>    - Second electron RK ( \c irksub==2 ): full electron step
!>  -# keep temporary record of initial values if in RK 1
!>  -# evaluate perturbed fields using the chosen gyro-averaging scheme
!>  -# evaluate equilibrium fields at particle GC location
!>  -# calculate the time derivatives of the 5 quantities
!>  -# update the 5 quantities
!>  -# deal with out of boundary particles and restore profiles
!>  -# calculate time diagnostics
!>  .
!> \n
!>
!> @param[in,out] zpart
!> particle information array, stands for zion,zelectron,zfast, or zfaste 
!
!> @param[in,out] zpart0
!> initial values of zpart, set in RK 1, used in RK 2
!
!> @param[in] wppart
!> linear interpolation weight in psi for gyro averaging locations, calculated
!> in subroutine locateParticle
!
!> @param[in] wtpart0
!> linear interpolation weight for upper theta grid on inner flux surface
!
!> @param[in] wtpart1
!> linear interpolation weight for upper theta grid on outer flux surface
!
!> @param[out] data1d 
!> 1D diagnostic data
!
!> @param[in] jtpart0
!> index of the lower theta grid on inner flux surface
!
!> @param[in] jtpart1
!> index of the lower theta grid on outer flux surface
!
!> @param[out] diagpart
!> time diagnostic data
!
!> @param[in] wzpart
!> linear interpolation weight in zeta
!
!> @param[in] rdtem
!> UNKNOWN
!
!> @param[in] kapan
!> density gradient term:\f$ \frac{1}{n}\frac{dn}{d\psi} \f$
!
!> @param[in] kapat
!> temperature gradient term:\f$ \frac{1}{T}\frac{dT}{d\psi} \f$
!
!> @param[in] meshn
!> equilibrium density on simulation grid
!
!> @param[in] mesht
!> equilibrium temperature on simulation grid
!
!> @param[in] pflux
!> UNKNOWN
!
!> @param[in] qpart
!> particle charge
!
!> @param[in] apart
!> particle mass
!
!> @param[in] pload
!> particle load scheme, see iload for details.
!
!> @param[in] ngyro
!> gyro-averaging scheme, normally either 1 or 4
!
!> @param[in] mpmax
!> maximum particle number on single MPI, determines array size
!
!> @param[in] mp
!> actual particle number on current MPI
!
!> @param[in] nparam
!> number of particle information entries, different kinds of particles may
!> have different number of information entries in zpart array.
!
!> @param[in] ncyclep
!> total subcycling number
!
!> @param[in] icycle
!> current subcycle count
!
!> @param[in] irksub
!> current RK step in a subcycle
!
!> @param[in,out] mp1
!> UNKNOWN
!
!> @param[in,out] zpart1
!> UNKNOWN
!
!> @param[in] ihybrid
!> current hybrid loop number
!
!> @param[in] phit
!> time derivative of electic potential used in hybrid model:
!> \f$ \frac{\partial \phi}{\partial t} \f$
!
!> @param[in] dnt
!> time derivative of density used in hybrid model:
!> \f$ \frac{\partial n}{\partial t}\f$
subroutine pushParticle(zpart,zpart0,wppart,wtpart0,wtpart1,&
    data1d,jtpart0,jtpart1,diagpart,wzpart,rdtem,rdden,kapan,kapat,meshn,mesht,&
    pflux,qpart,apart,pload,ngyro,mpmax,mp,nparam,ncyclep,icycle,irksub,mp1,&
    zpart1,ihybrid,phit,dnt)
  use system_env
  use precision, only: lk, mpi_Rsize
  use global_parameters,only: nonlinear,irk,tstep,magnetic,rg0,rg1,ecurrent,&
    paranl,psi1,pi,pi2,psi0,idiag,mpsi,mpsilow,island,gtcout,mype,istep,&
    torbound,rho0,track_particles,nhybrid,numereq,mthetamax,irestore,deltabpara,&
    BC_type, BC_weight, ipoincare, cs_method
  use particle_array,only: mpdiag,mpdata1d,meshne,meshte,kapane,betae,&
    sd_vc,sd_v0,sd_widthInv,sd_l0,dbparat
  use field_array,only: deltar,deltap,deltat,deltaz,gradphi,gradapara,&
    gradphieff,sapara,psimesh,nmodes,nmode,zeta0,phip00,solvermethod,gradpsi,&
    maxbfield,minbfield,thetabmin,thetabmax,thetaupp,thetadwn,bpara,gradbpara,&
    sdelfluidne,sdeltapsi,taparanl_weight,dBdotdE
  use equilibrium,only:&
    lsp,spdpsi_inv,spdpsi,lst,spdtheta_inv,spdtheta,lszeta,spdzeta,spdzeta_inv,&
    nzsp_sec,qpsi,gpsi,cpsi,ropp,erpp,rgpsi,bsp,xsp,spdim,mesher,alphasp,ndim
  use magnetic_island,only: gradalphaIs,alphaIs
  use spline_function, only: sprgpsi,spb
  use boundary_decay, only: weight_bound, decay_coeff 
  use conservative_scheme, only: tapara_na_weight,gradapara_na
  implicit none

  !declaration of the dummy arguments
  integer pload,ngyro,mpmax,mp,nparam
  integer,optional :: ncyclep,icycle,irksub,ihybrid,mp1
  integer,dimension(:,:) :: jtpart0,jtpart1
  real(lk) qpart,apart
  real(lk),dimension(:) :: wzpart
  real(lk),dimension(0:) :: kapan,kapat,meshn,mesht,pflux,rdtem,rdden
  real(lk),dimension(mpdiag) :: diagpart
  real(lk),dimension(:,:) :: zpart,zpart0,wppart,wtpart0,wtpart1
  real(lk),dimension(0:,:) :: data1d
  real(lk),dimension(:,:),optional :: zpart1
  real(lk),dimension(0:,:),optional :: phit,dnt

  !declaration of the local variables
  integer i,m,ii,igyro,ij,irsmooth,icount,iir(mpmax),isp,jst,ksz,b_index,isIhybridOne,& 
    nb,ngyro0,igyro0,&
    iout ! counter for particles going out of radial boundaries per time step
  real(lk) rdum,tdum,dtime,q,q_inv,rg,b,g,gp,ri,rip,dbdp,dbdt,dbdz,dedb,dedb_perp,dedb_para,deni,bcurlb,&
    gqi,deltaf,upara,energy,kappa,dptdp,dptdt,dptdz,vdr,pdot,tdot,zdot,rdot,&
    wdot,wdrift,delr,wz1,wz0,wp1,wp0,wt11,wt10,wt01,wt00,tp_inv,epara,wpara,&
    wdrive,psitmp,thetatmp,zetatmp,cmratio,cinv,dapdp,dapdt,dapdz,dpx,dp2,dtx,&
    dt2,dzx,dx(spdim),lam,dlamdp,dlamdt,dlamdz,plampt,&
    kappan,kappat,b_inv,phipEq,dzonal,&
    angmom,fullf,dp1,upara0,dupara0dp,kappav,dedb0,energy0,majorr,dmajorrdp,&
    gyrodum,gyrodum0,paxis,bdum,vtot,dist_term,eperp,&
    dtem(0:mpsi),dden(0:mpsi),ddum(0:mpsi),vdrtmp(0:mpsi),dmark(0:mpsi),fden(0:mpsi),&
    diagtmp(mpdiag),data1dtmp(0:mpsi,mpdata1d),delp(mpsi),deltab(0:1),&
    pdum,cost,sint,xdot,ydot,xdum,ydum,&
    ! analytic island terms
    wisland,lami,dlamidp,dlamidt,dlamidz,i1,i2,i3,i4,wigc(4,mpmax),&
    ! equilibrium island terms
    ! equilibrium island spline delta in zeta temporarily different to the
    ! geometry spline delta zeta. And the delta array is always 3D,thus 27
    ! elements.
    dzxisland,dxisland(27),lami_eq,dlamidp_eq,dlamidt_eq,dlamidz_eq,bc_decay,vap,&
    ! delta bpara related terms
    domegap,deltaf0,vbpara,ebpara,dbparadp,dbparadt,dbparadz,&
    bpara1,bpara2,bpara3,bpara4,bpara5,wbparagc(5,mpmax),&
    !Non-adiabatic electric field
    c1,c2,c3,c4,wcgc(4,mpmax),epara_na,vap_na,dapara_nadp,dapara_nadt,dapara_nadz,&
    ! perturb field switches
    perturb,&
    ! nonlinear inductive electric field and ponderomotive electric field
    wbgc_nl(2,mpmax),b1_nl,b2_nl,epara_pf,&
#ifdef _hybrid
    dphieffdp,dphieffdt,dphieffdz,kappas,dupara0dt,kappavt,wdrift0,wpara0,&
    qp,gp2,rip2,vet,wdriftind,delz,np_inv,dphiinddp,dphiinddt,dphiinddz,&
    pptpt,dpsidp,dpsidt,dpsidz,e1,e2,e3,e4,b1,b2,b3,b4,b5,b6,b7,b8,&
    b9,b10,b11,b12,delt(0:mpsi),wpgc(4,mpmax),wbgc(12,mpmax)
#else
    e1,e2,e3,&
    b1,b2,b3,b4,b5,wpgc(3,mpmax),wbgc(5,mpmax)
#endif

  !$acc declare create(iir)
  !$acc declare create(delp,wpgc,vdrtmp,wbgc)
  !$acc declare create(dtem,dmark,dden,fden)
  !$acc declare create(wigc,wbparagc)
  !$acc declare create(wcgc)
  !$acc declare create(wbgc_nl)
  logical subcycle, isNhybrid
#ifdef _FRC
  real(8) dbesjn,kperprhoi,zetagyravg
#else
  real(lk) kperprhoi,zetagyravg
#endif
  subcycle=.false.
  if(present(ncyclep))then
    if(ncyclep>0) subcycle=.true.
  endif
!$acc parallel loop gang vector  
  do i=1,mp
    wbparagc(:,i)=0.0_lk
    wcgc(:,i)=0.0_lk
    wbgc_nl(:,i)=0.0_lk
  enddo
!$acc end parallel

  isIhybridOne = 0
  isNhybrid = .false.
  if (present(ihybrid)) then
     if (ihybrid==1) isIhybridOne=1
     if (ihybrid==nhybrid) isNhybrid=.True.
  endif

  !paxis=0.5*(8.0*rho0)**2
  paxis=0.0
  if(psi0<1.0e-8)paxis=psimesh(1)
  delr=1.0/deltar
#ifdef _hybrid
  delt=1.0/deltat !!no 2pi/deltar in version3
  delz=1.0/deltaz
#endif
  cmratio=qpart/apart
  cinv=1.0/qpart
  vdrtmp=0.0
  delp=1.0/deltap

  ! mode parameters
  perturb=real(nonlinear)

  !$acc update device(delp)

! for out-of-boundary particles
  deltab(0)=(maxbfield(0)-minbfield(0))/real(mthetamax/2 -1)
  deltab(1)=(maxbfield(1)-minbfield(1))/real(mthetamax/2 -1)

  if(.not. subcycle)then
    if(irk==1)then
      ! 1st step of Runge-Kutta method
      dtime=0.5*tstep
#ifdef _OPENACC
      !$acc parallel loop gang vector
#else
      !$omp parallel do private(m)
#endif
      do m=1,mp
        zpart0(1:nparam,m)=zpart(1:nparam,m)
      enddo
      !$acc end parallel

    ! 2nd step of Runge-Kutta method
    else
      dtime=tstep
      if(nonlinear==1)vdrtmp=pflux
    endif
  else
    if(irksub==1)then
      ! 1st step of Runge-Kutta method
      dtime=0.25*tstep/real(ncyclep)
      if(icycle==1)then
        if(irk*isIhybridOne==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)
          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)
          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:nparam,m)=zpart(1:nparam,m)
      enddo
      !$acc end parallel

    ! 2nd step of Runge-Kutta method
    else
      dtime=0.5*tstep/real(ncyclep)
      if(nonlinear==1 .and. irk==2 .and. isNhybrid)vdrtmp=pflux
    endif
  endif

  !$acc update device(vdrtmp)

  ! gather e_field using ngyro-point gyro-averaging
  gyrodum=1.0/real(ngyro)
  gyrodum0=1.0/4.0 !For gyro-surface averaging
  ! island modification of lambda and lambda gradient, represented by wigc(1,m) and wigc(2,m) wigc(3,m) wigc(4,m)
  if(island==1)then
#ifdef _OPENACC
    !$acc parallel loop gang vector
#else
    !$omp parallel do private(m,igyro,i1,i2,i3,i4,wz1,wz0,wp1,wp0,wt00,wt10,wt01,wt11,ij)
#endif
    do m=1,mp
      i1=0.0
      i2=0.0
      i3=0.0
      i4=0.0
      wz1=wzpart(m)                !weight for upper toroidal grid
      wz0=1.0-wz1                 !weight for lower toroidal grid

      do igyro=1,ngyro
        wp1=wppart(igyro,m)       !outer flux surface
        wp0=1.0-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
        i1=i1+wt00*(wz0*gradalphaIs(1,0,ij)+wz1*gradalphaIs(1,1,ij))
        i2=i2+wt00*(wz0*gradalphaIs(2,0,ij)+wz1*gradalphaIs(2,1,ij))
        i3=i3+wt00*(wz0*gradalphaIs(3,0,ij)+wz1*gradalphaIs(3,1,ij))
        i4=i4+wt00*(wz0*alphaIs(0,ij)+wz1*alphaIs(1,ij))


        ij=ij+1                  !upper poloidal grid on inner flux surface
        i1=i1+wt10*(wz0*gradalphaIs(1,0,ij)+wz1*gradalphaIs(1,1,ij))
        i2=i2+wt10*(wz0*gradalphaIs(2,0,ij)+wz1*gradalphaIs(2,1,ij))
        i3=i3+wt10*(wz0*gradalphaIs(3,0,ij)+wz1*gradalphaIs(3,1,ij))
        i4=i4+wt10*(wz0*alphaIs(0,ij)+wz1*alphaIs(1,ij))

        ij=jtpart1(igyro,m)       !lower poloidal grid on outer flux surface
        i1=i1+wt01*(wz0*gradalphaIs(1,0,ij)+wz1*gradalphaIs(1,1,ij))
        i2=i2+wt01*(wz0*gradalphaIs(2,0,ij)+wz1*gradalphaIs(2,1,ij))
        i3=i3+wt01*(wz0*gradalphaIs(3,0,ij)+wz1*gradalphaIs(3,1,ij))
        i4=i4+wt01*(wz0*alphaIs(0,ij)+wz1*alphaIs(1,ij))

        ij=ij+1                  !upper poloidal grid on outer flux surface
        i1=i1+wt11*(wz0*gradalphaIs(1,0,ij)+wz1*gradalphaIs(1,1,ij))
        i2=i2+wt11*(wz0*gradalphaIs(2,0,ij)+wz1*gradalphaIs(2,1,ij))
        i3=i3+wt11*(wz0*gradalphaIs(3,0,ij)+wz1*gradalphaIs(3,1,ij))
        i4=i4+wt11*(wz0*alphaIs(0,ij)+wz1*alphaIs(1,ij))
      enddo

      wigc(1,m)=gyrodum*i1
      wigc(2,m)=gyrodum*i2
      wigc(3,m)=gyrodum*i3
      wigc(4,m)=gyrodum*i4
    enddo
    !$acc end parallel
  endif

  ! electrostatic fluctuations
  if(magnetic==0)then
#ifdef _OPENACC
    !$acc parallel loop gang vector
#else
#ifdef _hybrid
    !$omp parallel do private(m,igyro,e1,e2,e3,wz1,wz0,wp0,wp1,wt00,wt10,wt01,wt11,ij,e4)
#else
    !$omp parallel do private(m,igyro,e1,e2,e3,wz1,wz0,wp0,wp1,wt00,wt10,wt01,wt11,ij)
#endif
#endif
    do m=1,mp
      e1=0.0
      e2=0.0
      e3=0.0
#ifdef _hybrid
      e4=0.0
#endif

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

      do igyro=1,ngyro
        wp1=wppart(igyro,m)       !outer flux surface
        wp0=1.0-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))
#ifdef _hybrid
        e4=e4+wt00*(wz0*phit(0,ij)+wz1*phit(1,ij))
#endif

        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))
#ifdef _hybrid
        e4=e4+wt10*(wz0*phit(0,ij)+wz1*phit(1,ij))
#endif

        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))
#ifdef _hybrid
        e4=e4+wt01*(wz0*phit(0,ij)+wz1*phit(1,ij))
#endif

        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))

#ifdef _hybrid
        e4=e4+wt11*(wz0*phit(0,ij)+wz1*phit(1,ij))
#endif
      enddo

      wpgc(1,m)=gyrodum*e1
      wpgc(2,m)=gyrodum*e2
      wpgc(3,m)=gyrodum*e3
#ifdef _hybrid
      wpgc(4,m)=gyrodum*e4
#endif

      wbgc(:,m)=0.0
    enddo
    !$acc end parallel
  else
    ! electromagnetic fields
#ifdef _OPENACC
    !$acc parallel loop gang vector
#else
#ifdef _hybrid
    !$omp parallel do private(m,igyro,e1,e2,e3,b1,b2,b3,b4,b5,wz1,wz0,wp0,wp1,wt00,wt10,wt01,wt11,ij)&
    !$omp& private(e4,b6,b7,b8,b9,b10,b11,b12,b1_nl,b2_nl)
#else
    !$omp parallel do private(m,igyro,e1,e2,e3,b1,b2,b3,b4,b5,b1_nl,b2_nl,wz1,wz0,wp0,wp1,wt00,wt10,wt01,wt11,ij)
#endif
#endif
    do m=1,mp
       !$acc cache(readonly:gradphi,dnt,gradapara,gradpsi,sapara,sdelfluidne,sdeltapsi)
      e1=0.0
      e2=0.0
      e3=0.0
#ifdef _hybrid
      e4=0.0
#endif
      b1=0.0
      b2=0.0
      b3=0.0
      b4=0.0
      b5=0.0
#ifdef _hybrid
      b6=0.0
      b7=0.0
      b8=0.0
      b9=0.0
      b10=0.0
      b11=0.0
      b12=0.0
#endif
      b1_nl=0.0     
      b2_nl=0.0

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

      do igyro=1,ngyro
        wp1=wppart(igyro,m)       !outer flux surface
        wp0=1.0-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))
#ifdef _hybrid
        e4=e4+wt00*(wz0*dnt(0,ij)+wz1*dnt(1,ij))
#endif
        b1=b1+wt00*(wz0*gradapara(1,0,ij)+wz1*gradapara(1,1,ij))
        b2=b2+wt00*(wz0*gradapara(2,0,ij)+wz1*gradapara(2,1,ij))
        b3=b3+wt00*(wz0*gradapara(3,0,ij)+wz1*gradapara(3,1,ij))
        b4=b4+wt00*(wz0*gradphieff(3,0,ij)+wz1*gradphieff(3,1,ij))
        b5=b5+wt00*(wz0*sapara(0,ij)+wz1*sapara(1,ij))
#ifdef _hybrid
        b6=b6+wt00*(wz0*gradphieff(1,0,ij)+wz1*gradphieff(1,1,ij))
        b7=b7+wt00*(wz0*gradphieff(2,0,ij)+wz1*gradphieff(2,1,ij))
        b8=b8+wt00*(wz0*gradpsi(1,0,ij)+wz1*gradpsi(1,1,ij))
        b9=b9+wt00*(wz0*gradpsi(2,0,ij)+wz1*gradpsi(2,1,ij))
        b10=b10+wt00*(wz0*gradpsi(3,0,ij)+wz1*gradpsi(3,1,ij))
        b11=b11+wt00*(wz0*sdelfluidne(0,ij)+wz1*sdelfluidne(1,ij))
        b12=b12+wt00*(wz0*sdeltapsi(0,ij)+wz1*sdeltapsi(1,ij))
#endif
        b1_nl=b1_nl+wt00*(wz0*dBdotdE(0,ij)+wz1*dBdotdE(1,ij))
        b2_nl=b2_nl+wt00*(wz0*taparanl_weight(0,ij)+wz1*taparanl_weight(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))
#ifdef _hybrid
        e4=e4+wt10*(wz0*dnt(0,ij)+wz1*dnt(1,ij))
#endif
        b1=b1+wt10*(wz0*gradapara(1,0,ij)+wz1*gradapara(1,1,ij))
        b2=b2+wt10*(wz0*gradapara(2,0,ij)+wz1*gradapara(2,1,ij))
        b3=b3+wt10*(wz0*gradapara(3,0,ij)+wz1*gradapara(3,1,ij))
        b4=b4+wt10*(wz0*gradphieff(3,0,ij)+wz1*gradphieff(3,1,ij))
        b5=b5+wt10*(wz0*sapara(0,ij)+wz1*sapara(1,ij))
#ifdef _hybrid
        b6=b6+wt10*(wz0*gradphieff(1,0,ij)+wz1*gradphieff(1,1,ij))
        b7=b7+wt10*(wz0*gradphieff(2,0,ij)+wz1*gradphieff(2,1,ij))
        b8=b8+wt10*(wz0*gradpsi(1,0,ij)+wz1*gradpsi(1,1,ij))
        b9=b9+wt10*(wz0*gradpsi(2,0,ij)+wz1*gradpsi(2,1,ij))
        b10=b10+wt10*(wz0*gradpsi(3,0,ij)+wz1*gradpsi(3,1,ij))
        b11=b11+wt10*(wz0*sdelfluidne(0,ij)+wz1*sdelfluidne(1,ij))
        b12=b12+wt10*(wz0*sdeltapsi(0,ij)+wz1*sdeltapsi(1,ij))
#endif
        b1_nl=b1_nl+wt10*(wz0*dBdotdE(0,ij)+wz1*dBdotdE(1,ij))
        b2_nl=b2_nl+wt10*(wz0*taparanl_weight(0,ij)+wz1*taparanl_weight(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))
#ifdef _hybrid
        e4=e4+wt01*(wz0*dnt(0,ij)+wz1*dnt(1,ij))
#endif
        b1=b1+wt01*(wz0*gradapara(1,0,ij)+wz1*gradapara(1,1,ij))
        b2=b2+wt01*(wz0*gradapara(2,0,ij)+wz1*gradapara(2,1,ij))
        b3=b3+wt01*(wz0*gradapara(3,0,ij)+wz1*gradapara(3,1,ij))
        b4=b4+wt01*(wz0*gradphieff(3,0,ij)+wz1*gradphieff(3,1,ij))
        b5=b5+wt01*(wz0*sapara(0,ij)+wz1*sapara(1,ij))
#ifdef _hybrid
        b6=b6+wt01*(wz0*gradphieff(1,0,ij)+wz1*gradphieff(1,1,ij))
        b7=b7+wt01*(wz0*gradphieff(2,0,ij)+wz1*gradphieff(2,1,ij))
        b8=b8+wt01*(wz0*gradpsi(1,0,ij)+wz1*gradpsi(1,1,ij))
        b9=b9+wt01*(wz0*gradpsi(2,0,ij)+wz1*gradpsi(2,1,ij))
        b10=b10+wt01*(wz0*gradpsi(3,0,ij)+wz1*gradpsi(3,1,ij))
        b11=b11+wt01*(wz0*sdelfluidne(0,ij)+wz1*sdelfluidne(1,ij))
        b12=b12+wt01*(wz0*sdeltapsi(0,ij)+wz1*sdeltapsi(1,ij))
#endif
        b1_nl=b1_nl+wt01*(wz0*dBdotdE(0,ij)+wz1*dBdotdE(1,ij))
        b2_nl=b2_nl+wt01*(wz0*taparanl_weight(0,ij)+wz1*taparanl_weight(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))
#ifdef _hybrid
        e4=e4+wt11*(wz0*dnt(0,ij)+wz1*dnt(1,ij))
#endif
        b1=b1+wt11*(wz0*gradapara(1,0,ij)+wz1*gradapara(1,1,ij))
        b2=b2+wt11*(wz0*gradapara(2,0,ij)+wz1*gradapara(2,1,ij))
        b3=b3+wt11*(wz0*gradapara(3,0,ij)+wz1*gradapara(3,1,ij))
        b4=b4+wt11*(wz0*gradphieff(3,0,ij)+wz1*gradphieff(3,1,ij))
        b5=b5+wt11*(wz0*sapara(0,ij)+wz1*sapara(1,ij))
#ifdef _hybrid
        b6=b6+wt11*(wz0*gradphieff(1,0,ij)+wz1*gradphieff(1,1,ij))
        b7=b7+wt11*(wz0*gradphieff(2,0,ij)+wz1*gradphieff(2,1,ij))
        b8=b8+wt11*(wz0*gradpsi(1,0,ij)+wz1*gradpsi(1,1,ij))
        b9=b9+wt11*(wz0*gradpsi(2,0,ij)+wz1*gradpsi(2,1,ij))
        b10=b10+wt11*(wz0*gradpsi(3,0,ij)+wz1*gradpsi(3,1,ij))
        b11=b11+wt11*(wz0*sdelfluidne(0,ij)+wz1*sdelfluidne(1,ij))
        b12=b12+wt11*(wz0*sdeltapsi(0,ij)+wz1*sdeltapsi(1,ij))
#endif
        b1_nl=b1_nl+wt11*(wz0*dBdotdE(0,ij)+wz1*dBdotdE(1,ij))
        b2_nl=b2_nl+wt11*(wz0*taparanl_weight(0,ij)+wz1*taparanl_weight(1,ij))
      enddo

      wpgc(1,m)=gyrodum*e1
      wpgc(2,m)=gyrodum*e2
      wpgc(3,m)=gyrodum*e3
#ifdef _hybrid
      wpgc(4,m)=gyrodum*e4
#endif
      wbgc(1,m)=gyrodum*b1
      wbgc(2,m)=gyrodum*b2
      wbgc(3,m)=gyrodum*b3
      wbgc(4,m)=gyrodum*b4
      wbgc(5,m)=gyrodum*b5
#ifdef _hybrid
      wbgc(6,m)=gyrodum*b6
      wbgc(7,m)=gyrodum*b7
      wbgc(8,m)=gyrodum*b8
      wbgc(9,m)=gyrodum*b9
      wbgc(10,m)=gyrodum*b10
      wbgc(11,m)=gyrodum*b11
      wbgc(12,m)=gyrodum*b12
#endif
      wbgc_nl(1,m)=gyrodum*b1_nl
      wbgc_nl(2,m)=gyrodum*b2_nl
    enddo
    !$acc end parallel
  endif
  
  ! extra flux-surface gyro averaging for bpara formalism
  ! flux surface gyro averaging is always done with 4-point
  igyro0=9
  ngyro0=12
  ! For species without normal gyro averaging, the flux-surface gyro averaging is turned off
  if(ngyro==1)then
    igyro0=1
    ngyro0=1
    gyrodum0=1
  endif

  if(deltabpara==1)then
#ifdef _OPENACC
    !$acc parallel loop gang vector
#else

#ifdef _hybrid
    !$omp parallel do private(m,igyro,bpara1,bpara2,bpara3,bpara4,bpara5,wz1,wz0,wp1,wp0,wt00,wt10,wt01,wt11,ij)
#else
    !$omp parallel do private(m,igyro,bpara1,bpara2,bpara3,wz1,wz0,wp1,wp0,wt00,wt10,wt01,wt11,ij)
#endif

#endif
    do m=1,mp
      bpara1=0.0  
      bpara2=0.0
      bpara3=0.0
#ifdef _hybrid
      bpara4=0.0
      bpara5=0.0
#endif
      wz1=wzpart(m)                !weight for upper toroidal grid
      wz0=1.0-wz1                 !weight for lower toroidal grid

      do igyro=igyro0,ngyro0
        wp1=wppart(igyro,m)       !outer flux surface
        wp0=1.0-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
        bpara1=bpara1+wt00*(wz0*gradbpara(1,0,ij)+wz1*gradbpara(1,1,ij))
        bpara2=bpara2+wt00*(wz0*gradbpara(2,0,ij)+wz1*gradbpara(2,1,ij))
        bpara3=bpara3+wt00*(wz0*gradbpara(3,0,ij)+wz1*gradbpara(3,1,ij))
#ifdef _hybrid
        bpara4=bpara4+wt00*(wz0*bpara(0,ij)+wz1*bpara(1,ij))
        bpara5=bpara5+wt00*(wz0*dbparat(0,ij)+wz1*dbparat(1,ij))
#endif

        ij=ij+1                  !upper poloidal grid on inner flux surface
        bpara1=bpara1+wt10*(wz0*gradbpara(1,0,ij)+wz1*gradbpara(1,1,ij))
        bpara2=bpara2+wt10*(wz0*gradbpara(2,0,ij)+wz1*gradbpara(2,1,ij))
        bpara3=bpara3+wt10*(wz0*gradbpara(3,0,ij)+wz1*gradbpara(3,1,ij))
#ifdef _hybrid
        bpara4=bpara4+wt10*(wz0*bpara(0,ij)+wz1*bpara(1,ij))
        bpara5=bpara5+wt10*(wz0*dbparat(0,ij)+wz1*dbparat(1,ij))
#endif
        ij=jtpart1(igyro,m)       !lower poloidal grid on outer flux surface
        bpara1=bpara1+wt01*(wz0*gradbpara(1,0,ij)+wz1*gradbpara(1,1,ij))
        bpara2=bpara2+wt01*(wz0*gradbpara(2,0,ij)+wz1*gradbpara(2,1,ij))
        bpara3=bpara3+wt01*(wz0*gradbpara(3,0,ij)+wz1*gradbpara(3,1,ij))
#ifdef _hybrid
        bpara4=bpara4+wt01*(wz0*bpara(0,ij)+wz1*bpara(1,ij))
        bpara5=bpara5+wt01*(wz0*dbparat(0,ij)+wz1*dbparat(1,ij))
#endif
        ij=ij+1                  !upper poloidal grid on outer flux surface
        bpara1=bpara1+wt11*(wz0*gradbpara(1,0,ij)+wz1*gradbpara(1,1,ij))
        bpara2=bpara2+wt11*(wz0*gradbpara(2,0,ij)+wz1*gradbpara(2,1,ij))
        bpara3=bpara3+wt11*(wz0*gradbpara(3,0,ij)+wz1*gradbpara(3,1,ij))
#ifdef _hybrid 
        bpara4=bpara4+wt11*(wz0*bpara(0,ij)+wz1*bpara(1,ij))
        bpara5=bpara5+wt11*(wz0*dbparat(0,ij)+wz1*dbparat(1,ij))
#endif
      enddo

      wbparagc(1,m)=gyrodum0*bpara1
      wbparagc(2,m)=gyrodum0*bpara2
      wbparagc(3,m)=gyrodum0*bpara3
#ifdef _hybrid 
      wbparagc(4,m)=gyrodum0*bpara4
      wbparagc(5,m)=gyrodum0*bpara5
#endif      
    enddo
    !$acc end parallel
  endif

  if(cs_method==1)then
#ifdef _OPENACC
    !$acc parallel loop gang vector
#else
    !$omp parallel do private(m,igyro,c1,c2,c3,c4,wz1,wz0,wp1,wp0,wt00,wt10,wt01,wt11,ij)
#endif
    do m=1,mp
      c1=0.0
      c2=0.0
      c3=0.0
      c4=0.0

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

      do igyro=1,ngyro
        wp1=wppart(igyro,m)       !outer flux surface
        wp0=1.0-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
        c1=c1+wt00*(wz0*gradapara_na(1,0,ij)+wz1*gradapara_na(1,1,ij))    
        c2=c2+wt00*(wz0*gradapara_na(2,0,ij)+wz1*gradapara_na(2,1,ij))
        c3=c3+wt00*(wz0*gradapara_na(3,0,ij)+wz1*gradapara_na(3,1,ij))
        c4=c4+wt00*(wz0*tapara_na_weight(0,ij)+wz1*tapara_na_weight(1,ij))

        ij=ij+1                  !upper poloidal grid on inner flux surface
        c1=c1+wt10*(wz0*gradapara_na(1,0,ij)+wz1*gradapara_na(1,1,ij))
        c2=c2+wt10*(wz0*gradapara_na(2,0,ij)+wz1*gradapara_na(2,1,ij))
        c3=c3+wt10*(wz0*gradapara_na(3,0,ij)+wz1*gradapara_na(3,1,ij))
        c4=c4+wt10*(wz0*tapara_na_weight(0,ij)+wz1*tapara_na_weight(1,ij))

        ij=jtpart1(igyro,m)       !lower poloidal grid on outer flux surface
        c1=c1+wt01*(wz0*gradapara_na(1,0,ij)+wz1*gradapara_na(1,1,ij))
        c2=c2+wt01*(wz0*gradapara_na(2,0,ij)+wz1*gradapara_na(2,1,ij))
        c3=c3+wt01*(wz0*gradapara_na(3,0,ij)+wz1*gradapara_na(3,1,ij))
        c4=c4+wt01*(wz0*tapara_na_weight(0,ij)+wz1*tapara_na_weight(1,ij))

        ij=ij+1                  !upper poloidal grid on outer flux surface
        c1=c1+wt11*(wz0*gradapara_na(1,0,ij)+wz1*gradapara_na(1,1,ij))
        c2=c2+wt11*(wz0*gradapara_na(2,0,ij)+wz1*gradapara_na(2,1,ij))
        c3=c3+wt11*(wz0*gradapara_na(3,0,ij)+wz1*gradapara_na(3,1,ij))
        c4=c4+wt11*(wz0*tapara_na_weight(0,ij)+wz1*tapara_na_weight(1,ij))

      enddo

      wcgc(1,m)=gyrodum*c1
      wcgc(2,m)=gyrodum*c2
      wcgc(3,m)=gyrodum*c3
      wcgc(4,m)=gyrodum*c4
    enddo
    !$acc end parallel
  endif

  ! update GC position. Assuming psi0>spdpsi; will be relaxed later
#ifdef _OPENACC
  !$acc parallel loop gang vector private(dx(spdim),dxisland(27))
#else
  !$omp parallel do private(m,psitmp,thetatmp,isp,dpx,dp1,dp2,jst,dtx,dt2,dzx,dx,&
  !$omp& dxisland,rg,q,g,&
  !$omp& gp,ri,rip,b,dbdp,dbdt,dbdz,ii,dedb,dedb_perp,dedb_para,deni,gqi,upara,energy,dptdp,dptdt,bcurlb,&
  !$omp& dptdz,dapdp,dapdt,dapdz,dlamidt,dlamidz,dlamidp_eq,dlamidt_eq,dlamidz_eq,&
  !$omp& epara,vdr,vap,wdrive,wpara,wdrift,wdot,rdot,&
  !$omp& pdot,tdot,zdot,kappa,lam,dlamdp,dlamdt,dlamdz,plampt,kappan,&
  !$omp& kappat,tp_inv,upara0,dupara0dp,kappav,b_inv,majorr,dmajorrdp,dedb0,&
  !$omp& energy0,phipEq,kperprhoi,zetagyravg,vtot,dist_term,eperp,dzonal,&
  !$omp& vbpara,ebpara,dbparadp,dbparadt,dbparadz,deltaf0,&
  !$omp& epara_na,epara_pf,dapara_nadp,dapara_nadt,dapara_nadz,vap_na,&
  !$omp& lami_eq,zetatmp,ksz,lami,dlamidp,wisland,pdum,cost,sint,xdot,ydot,xdum,ydum)&
#ifdef _hybrid
  !$omp& private(kappavt,vet,wdriftind,gp2,rip2,pptpt,wp0,wp1,dupara0dt,&
  !$omp& wdrift0,wpara0,dphieffdp,dphieffdt,dphieffdz,qp,kappas,&
  !$omp& q_inv,np_inv,dphiinddp,dphiinddt,dphiinddz,dpsidp,dpsidt,dpsidz)
#else
  !$omp& private(q_inv)
#endif
#endif
  do m=1,mp
    psitmp=zpart(1,m)
    isp=max(1,min(lsp-1,ceiling(psitmp*spdpsi_inv)))
    ! For 1D profiles, spline is regular near the axis
    dpx=psitmp-spdpsi*real(isp-1)
    dp2=dpx*dpx
    ! poloidal spline is regular
    thetatmp=zpart(2,m)
    jst=max(1,min(lst-1,ceiling(thetatmp*spdtheta_inv)))
    dtx=thetatmp-spdtheta*real(jst-1)
    dt2=dtx*dtx

    dx(1)=1.0
    dx(2)=dpx
    dx(3)=dp2

    ! 1D spline in psi
    q=    qpsi(1,isp)   +qpsi(2,isp)*dx(2)   +qpsi(3,isp)*dx(3)
    g=    gpsi(1,isp)   +gpsi(2,isp)*dx(2)   +gpsi(3,isp)*dx(3)
    gp=                  gpsi(2,isp)       +gpsi(3,isp)*dx(2)*2.0
    ri=   cpsi(1,isp)   +cpsi(2,isp)*dx(2)   +cpsi(3,isp)*dx(3)
    rip=                 cpsi(2,isp)       +cpsi(3,isp)*dx(2)*2.0
    upara0=ropp(1,isp)  +ropp(2,isp)*dx(2)   +ropp(3,isp)*dx(3)
    dupara0dp=           ropp(2,isp)       +ropp(3,isp)*dx(2)*2.0
!    er=   erpp(1,isp)   +erpp(2,isp)*dx(2)   +erpp(3,isp)*dx(3)
    bcurlb=g*rip-ri*gp
#ifdef _hybrid
    ! To make the 2nd derivatives of I and g continuous,
    ! currently we use linear interpolations as a temporary solution.
    ! In the future, higher order spline functions of I and g need to be used
    gp2=2.0*(dpx*spdpsi_inv*(gpsi(3,isp+1)-gpsi(3,isp))+gpsi(3,isp))
    rip2=2.0*(dpx*spdpsi_inv*(cpsi(3,isp+1)-cpsi(3,isp))+cpsi(3,isp))
    qp=                  qpsi(2,isp)       +qpsi(3,isp)*dx(2)*2.0
#endif

    ! Now we deal with the 2D/3D equilibrium quantities: B, X and their derivatives
    ! For coordinates splines, first cell needs to be specially treated.
    ! radaial spline of rg & b avoids sigularity near axis: y=y1+y2*sqrt(x)+y3*x for psi<spdpsi
    if(isp==1)dpx=sqrt(dpx)
    dp2=dpx*dpx

    ! Evaluate rg with special dpsi.
    rg=  rgpsi(1,isp)  +rgpsi(2,isp)*dx(2)  +rgpsi(3,isp)*dx(3)

! 2D and 3D equilibrium spline arrays have different dimensions
! so they need to be treated differently.

! 2D equilibrium coordinate splines
#ifndef _TOROIDAL3D

    dx(1)=1.0_lk
    dx(2)=dpx
    dx(3)=dp2
    dx(4:6)=dx(1:3)*dtx
    dx(7:9)=dx(1:3)*dt2
    b=0.0_lk
    majorr=0.0_lk
    do ii = 1, 9
      b=b +bsp(ii,isp,jst)*dx(ii)
      majorr=majorr +xsp(ii,isp,jst)*dx(ii)
    enddo

    ! zeta-derivative is zero in 2D equilibria
    dbdz=0.0_lk

    ! theta-derivative
    dx(4:6)=dx(1:3)*1.0_lk
    dx(7:9)=dx(1:3)*dtx*2.0_lk
    dbdt=0.0_lk
    do ii = 4, 9
      dbdt=dbdt +bsp(ii,isp,jst)*dx(ii)
    enddo

! psi derivative needs special care when isp==1.
    dx(1)=0.0_lk
    dx(2)=1.0_lk
    dx(3)=dpx*2.0_lk
    if(isp==1) then
      dx(2)=0.5_lk/dpx
      dx(3)=1.0_lk
    endif
    dx(4:6)=dx(1:3)*dtx
    dx(7:9)=dx(1:3)*dt2
    dbdp=0.0_lk
    dmajorrdp=0.0_lk
    do ii = 2, 9 
      dbdp=dbdp +bsp(ii,isp,jst)*dx(ii)
      dmajorrdp=dmajorrdp +xsp(ii,isp,jst)*dx(ii)
    enddo

! 3D equilibrium coordinates splines
#else
    zetatmp = zpart(3,m)
! We first calculate the spline indexes. The psi and theta indexesare the same
! as calculated from other equilibrium splines.
! dzeta needs to be calculated again because of sub-section used in 3D equilibrium
    ksz=max(1,min(nzsp_sec,ceiling((zetatmp-zeta0)*spdzeta_inv)))
    dzx = zetatmp - (zeta0 + real(spdzeta*(ksz-1)))

! Now we evaluate b and majorr using spline coefficients
    dx(1)=1.0_lk
    dx(2)=dpx
    dx(3)=dp2
    dx(4:6)=dx(1:3)*dtx
    dx(7:9)=dx(1:3)*dt2
    dx(10:18)=dx(1:9)*dzx
    dx(19:27)=dx(10:18)*dzx
    b = 0.0_lk
    majorr = 0.0_lk
    do ii=1,27
      b = b + bsp(ii,isp,jst,ksz)*dx(ii)
      majorr = majorr + xsp(ii,isp,jst,ksz)*dx(ii)
    enddo
! psi derivative needs special care when isp==1.
    dx(1)=0.0_lk
    dx(2)=1.0_lk
    dx(3)=2.0_lk*dpx
    if(isp == 1) then
      dx(2) = 0.5_lk/dpx
      dx(3) = 1.0_lk
    endif
    dx(4:6)=dx(1:3)*dtx
    dx(7:9)=dx(1:3)*dt2
    dx(10:18)=dx(1:9)*dzx
    dx(19:27)=dx(10:18)*dzx
    dbdp=0.0_lk
    dmajorrdp=0.0_lk
    do ii=2,27
      dbdp = dbdp + bsp(ii,isp,jst,ksz)*dx(ii)
      dmajorrdp = dmajorrdp + xsp(ii,isp,jst,ksz)*dx(ii)
    enddo
! theta derivative can be easily evaluated
    dx(1:3)=0.0_lk
    dx(4)=1.0_lk
    dx(5)=dpx
    dx(6)=dp2
    dx(7:9)=dx(4:6)*dtx*2.0_lk
    dx(10:18)=dx(1:9)*dzx
    dx(19:27)=dx(10:18)*dzx
    dbdt = 0.0_lk
    do ii=4,27
      dbdt = dbdt + bsp(ii,isp,jst,ksz)*dx(ii)
    enddo
! zeta derivative can be easily evaluated
    dx(1:9)=0.0_lk
    dx(10)=1.0_lk
    dx(11)=dpx
    dx(12)=dp2
    dx(13:15)=dx(10:12)*dtx
    dx(16:18)=dx(13:15)*dtx
    dx(19:27)=dx(10:18)*dzx*2
    dbdz = 0.0_lk
    do ii=10,27
      dbdz = dbdz + bsp(ii,isp,jst,ksz)*dx(ii)
    enddo
! 2D and 3D equilibrium coordinates spline finishes
#endif

    b_inv=1.0_lk/b

    zetagyravg=1.0_lk

! Equilibrium island quantities. These calculations are hard coded inline to accelerate the execution.
! Created by L. Shi, 03/22/17
! Equilibrium island terms are by default set to 0
    lami_eq = 0.0_lk
    dlamidp_eq = 0.0_lk
    dlamidt_eq = 0.0_lk
    dlamidz_eq = 0.0_lk

! M3DC1 module is only available in 2D compilation
#ifndef _TOROIDAL3D
! We calculate equilibrium island when M3DC1 islands are included
    if(numereq==3 .and. ndim>1) then
! equilibrium island is expressed as a 3D spline
! We first calculate the spline indexes. The psi and theta indexesare the same
! as calculated from other equilibrium splines.
      dxisland(1)=1
      dxisland(2)=dpx
      dxisland(3)=dp2
      dxisland(4:6)=dxisland(1:3)*dtx
      dxisland(7:9)=dxisland(1:3)*dt2
! dzeta needs to be calculated again because of sub-section used in equilibrium
! island spline.
!fujy bugfixed in 4.4 beta version, add definition of zetatmp to M3DC1 RMP simulation
      zetatmp=zpart(3,m)
      ksz=max(1,min(nzsp_sec,ceiling((zetatmp-zeta0)*spdzeta_inv)))
      dzxisland = zetatmp - (zeta0 + real(spdzeta*(ksz-1)))
      dxisland(10:18)=dxisland(1:9)*dzxisland
      dxisland(19:27)=dxisland(10:18)*dzxisland
! Now we evaluate alpha_eq using spline coefficients
      do ii=1,27
        lami_eq = lami_eq+alphasp(ii,isp,jst,ksz)*dxisland(ii)
      enddo
! psi derivative needs special care when isp==1. But we assume here that the
! perturbed magnetic field is zero at the axis, and ignore this part.
!> @todo
!> isp==1 case needs to be taken care correctly when alpha_eq is not sufficiently
!> small near the axis.
      dxisland(1)=0
      dxisland(2)=1
      dxisland(3)=2*dpx
      dxisland(4:6)=dxisland(1:3)*dtx
      dxisland(7:9)=dxisland(1:3)*dt2
      dxisland(10:18)=dxisland(1:9)*dzxisland
      dxisland(19:27)=dxisland(10:18)*dzxisland
      do ii=2,27
        dlamidp_eq = dlamidp_eq+alphasp(ii,isp,jst,ksz)*dxisland(ii)
      enddo
! theta derivative can be easily evaluated
      dxisland(1:3)=0.0_lk
      dxisland(4)=1
      dxisland(5)=dpx
      dxisland(6)=dp2
      dxisland(7:9)=dxisland(4:6)*dtx*2
      dxisland(10:18)=dxisland(1:9)*dzxisland
      dxisland(19:27)=dxisland(10:18)*dzxisland
      do ii=4,27
        dlamidt_eq = dlamidt_eq+alphasp(ii,isp,jst,ksz)*dxisland(ii)
      enddo
! zeta derivative can be easily evaluated
      dxisland(1:9)=0.0_lk
      dxisland(10)=1
      dxisland(11)=dpx
      dxisland(12)=dp2
      dxisland(13:15)=dxisland(10:12)*dtx
      dxisland(16:18)=dxisland(13:15)*dtx
      dxisland(19:27)=dxisland(10:18)*dzxisland*2
      do ii=10,27
        dlamidz_eq = dlamidz_eq+alphasp(ii,isp,jst,ksz)*dxisland(ii)
      enddo
    endif
#endif

#ifdef _FRC
    q_inv=0.0_lk
    !set up toroidal gyroavg by multiplying by besselj(0,k_perp*rho_i)
    !rho_i from sqrt(mu)*sqrt(2)/sqrt(b)/q->zpart(6,m)*sqrt(2)/sqrt(b)
    !k_perp from n/(Rmaj)--> n/majorr
    kperprhoi=(zpart(6,m)*sqrt(2.0_lk*apart/b))*(real(nmode)/majorr)
    zetagyravg=dbesjn(0,kperprhoi)
    wpgc(1,m)=wpgc(1,m)*zetagyravg
    wpgc(2,m)=wpgc(2,m)*zetagyravg
    wpgc(3,m)=wpgc(3,m)*zetagyravg
    !wpgc(4,m)=wpgc(4,m)*zetagyravg
!#elif _CYLINDER
!    q_inv=0.0_lk
#else
    q_inv=1.0_lk/q
#endif
    ! perturbed electric field
    dptdp=wpgc(1,m)
    dptdt=wpgc(2,m)
    dptdz=wpgc(3,m)-wpgc(2,m)*q_inv
#ifdef _hybrid
    pptpt=wpgc(4,m)+wbparagc(5,m)*(1.0_lk-zpart(6,m)*zpart(6,m)*b/meshte(1))*b_inv
    ! effective potential
    dphieffdp=wbgc(6,m)
    dphieffdt=wbgc(7,m)
    dphieffdz=wbgc(4,m)-wbgc(7,m)*q_inv

    ! induced potential
    dphiinddp=dphieffdp-dptdp
    dphiinddt=dphieffdt-dptdt
    dphiinddz=dphieffdz-dptdz

    ! perturbed magnetic flux
    dpsidp=wbgc(8,m)
    dpsidt=wbgc(9,m)
    dpsidz=wbgc(10,m)-wbgc(9,m)*q_inv
#endif

    ! perturbed vector potential
    dapdp=wbgc(1,m)
    dapdt=wbgc(2,m)
    dapdz=wbgc(3,m)-wbgc(2,m)*q_inv

    dapara_nadp=wcgc(1,m)
    dapara_nadt=wcgc(2,m)
    dapara_nadz=wcgc(3,m)-wcgc(2,m)*q_inv

    ! self-consistent lambda term is calculated from delta_Apara
    ! lambda=apara*b_inv
    lam=wbgc(5,m)*b_inv
    dlamdp=(dapdp-lam*dbdp)*b_inv
    dlamdt=(dapdt-lam*dbdt)*b_inv
    dlamdz=(dapdz-lam*dbdz)*b_inv

    ! analytic island perturbed magnetic field
    lami=0.0_lk
    dlamidp=0.0_lk
    dlamidt=0.0_lk
    dlamidz=0.0_lk
    if(island==1)then
      lami=wigc(4,m)
      dlamidp=wigc(1,m)
      dlamidt=wigc(2,m)
    ! Kaisheng & Lei: Add dalpha/dtheta correction to evaluation of
    ! dalpha/dzeta, since gradient evaluates parallel derivative instead of
    ! zeta derivative.
      dlamidz=wigc(3,m)-wigc(2,m)*q_inv
    endif

    ! analytic and equilibrium island fields are summed up
    lami_eq = lami + lami_eq
    dlamidp_eq = dlamidp + dlamidp_eq
    dlamidt_eq = dlamidt + dlamidt_eq
    dlamidz_eq = dlamidz + dlamidz_eq

    ! Now calculate 1D interpolations along psi
    ii=max(1,min(mpsi,ceiling((rg-rg0)*delr)))    !radial grid on inner flux surface
    dp1=(psimesh(ii)-psitmp)*delp(ii) !weight for inner flux surface

    phipEq=-(dp1*mesher(ii-1)+(1.0-dp1)*mesher(ii)) ! mesher = -dphi/dpsi
    kappan=dp1*kapan(ii-1)+(1.0-dp1)*kapan(ii)
    kappat=dp1*kapat(ii-1)+(1.0-dp1)*kapat(ii)
    tp_inv=1.0/(dp1*mesht(ii-1)+(1.0-dp1)*mesht(ii))
    dbparadp=(wbparagc(1,m))*zpart(6,m)*zpart(6,m)*tp_inv/qpart
    dbparadt=(wbparagc(2,m))*zpart(6,m)*zpart(6,m)*tp_inv/qpart
    dbparadz=(wbparagc(3,m)-wbparagc(2,m)*q_inv)*zpart(6,m)*zpart(6,m)*tp_inv/qpart
#ifdef _hybrid
    np_inv=1.0/(dp1*meshn(ii-1)+(1.0-dp1)*meshn(ii))
#endif
    dzonal=dp1*phip00(ii-1)+(1.0-dp1)*phip00(ii)
    deni=1.0/(g*q + ri + (zpart(4,m)+perturb*lam+lami_eq)*bcurlb)
    gqi=1.0/(g*q+ri)

    upara0=upara0*majorr
    ! gradient majorr not included
    dupara0dp=dupara0dp*majorr

#ifdef _hybrid
    ! kappas=(partial S/partial psi/S), S defined in equilibrium current paper
    kappas=(g*rip2-ri*gp2)/bcurlb-(gp*q+g*qp+rip)*gqi

    upara0=upara0-real(ecurrent)*(2.0_lk*rho0*rho0/betae)*b*np_inv*bcurlb*gqi
    dupara0dp=dupara0dp-real(ecurrent)*(2.0_lk*rho0*rho0/betae)*np_inv*bcurlb*gqi*(dbdp+b*(kappan+kappas))
    dupara0dt=-real(ecurrent)*(2.0_lk*rho0*rho0/betae)*dbdt*np_inv*bcurlb*gqi
#endif
    upara=zpart(4,m)*b*cmratio
    dedb=zpart(4,m)*zpart(4,m)*b*cmratio+cinv*zpart(6,m)*zpart(6,m)
    dedb_perp=cinv*zpart(6,m)*zpart(6,m)
    dedb_para=zpart(4,m)*zpart(4,m)*b*cmratio
#ifdef _hybrid
    !> @todo 
    !> Check if this term is required for weight equation assuming shifted Maxwellian.
    !> Currently not used in the weight equation.
    dedb0=-cinv*upara*upara0*apart*b_inv  !Not currently implemented in delta h equation
#else
    dedb0=cinv*((upara*upara-upara*upara0)*apart*b_inv+zpart(6,m)*zpart(6,m))
#endif
    energy=0.5*apart*upara*upara+zpart(6,m)*zpart(6,m)*b
    energy0=0.5*apart*(upara-upara0)*(upara-upara0)+zpart(6,m)*zpart(6,m)*b

    ! Additional term for weight equation with slowing distribution
    if(pload==11)then !slowing down
      vtot=min(sqrt(2.0_lk*energy/apart),sd_v0)
      eperp=2*zpart(6,m)*zpart(6,m)/(apart*vtot**2)
      dist_term=-2*qpart*eperp**2*sd_widthInv**2/(zpart(6,m)*zpart(6,m))*(eperp-sd_l0)+3.0*qpart*vtot/(vtot**3+sd_vc**3)/apart
    else ! local maxwellian
        dist_term=qpart*tp_inv
    endif

    !> @note
    !> For slowing down distribution these kappa calculations may need to be modified.
    kappa=kappan+kappat*(energy0*tp_inv-1.5)
    kappav=-(upara-upara0)*dupara0dp*apart*tp_inv
#ifdef _hybrid
    kappavt=-(upara-upara0)*dupara0dt*apart*tp_inv
#endif

    ! dne/n0=e phi/t for electrostatic case
    ! parallel electric field epara, & dlambda/dt=-epara-gradphi(3,m)
    if (magnetic==0) then
#ifdef _FRC
      epara=-(dptdt+q*dptdz)*b*deni
#else
      epara=-wpgc(3,m)*q*b*deni
#endif
    ! static island doesn't have time derivative, so no extra term
    ! due to island here.
    !> @todo
    !> additional terms for non-static island.
      plampt=0.0
    else
      epara=-wbgc(4,m)*b*q*deni
      plampt=-(epara+wpgc(3,m)*q*b*deni)*b_inv
    endif
    ebpara=-(dbparadt+q*dbparadz)*b*deni

    if(cs_method==1)then
      epara_na=-wcgc(4,m)  ! from non-adiabatic ohm's law
    else
      epara_na=-perturb*wbgc_nl(2,m) ! from non-linear ohm's law
    endif
    plampt=plampt-epara_na*b_inv
    vap_na=-(ri*dapara_nadz-g*dapara_nadt)*upara*deni
    epara_pf=wbgc_nl(1,m)  ! Ponderomotive force term: dBdotdE=\delta B \cdot \grad \phi, minus sign is included in dBdotdE(0:1,mgrid)

    !grad B_para drift term
     vbpara=(ri*dbparadz-g*dbparadt)*deni

    ! ExB drift in radial direction for w-dot and flux diagnostics
    vdr=(ri*dptdz-g*dptdt)*deni
    ! vap: the magnetic field line perturbation induced radial drift velocity.
    ! Effect of vap in weight equation and particle trajectory is included
    ! in wisland and lam terms in the calculation of wdot and pdot.
    vap=-b*(ri*(dlamdz+dlamidz_eq)-g*(dlamdt+dlamidt_eq))*upara*deni
#ifdef _hybrid
    deltaf0=wbgc(11,m)+wbgc(12,m)*(kappan-kappa)+(1.0-zpart(6,m)*zpart(6,m)*b*tp_inv)*wbparagc(4,m)*tp_inv*b_inv
    !deltaf0=0.0
    vet=g*dptdp*deni
    ! In wdrive term: vbpara is for finite delta_B_para simulation, vap_na is non-adiabatic magnetic flutter motion in conservative scheme simulation.(vap_na=0.0 for cs_method=0 and vbpara=0.0 for deltabpara=0)
    wdrive=(vdr-vbpara+vap_na)*(kappa+kappav)+vet*kappavt
     !wdrive=vdr*(kappa+kappav)+vet*kappavt
    if (magnetic==0) then
      wpara=pptpt*qpart*tp_inv
      wdrift=(g*dedb*dbdt - ri*dedb*dbdz)*deni*dzonal*qpart*tp_inv+(g*dptdt-ri*dptdz)*deni*(dzonal+phipEq)*qpart*tp_inv
      wdriftind=0.0
    else
    ! The epara_na is non-adiabatic parallel electric field in conservative scheme, and epara_na is zero when CS method is not applied.
      wpara=-pptpt*np_inv+(kappan-kappa)*dphiinddt*q_inv+(epara_na+perturb*epara_pf)*(upara-upara0)*dist_term
      wdriftind=-(g*dbdt*dphiinddp-g*dbdp*dphiinddt+ri*dbdp*dphiinddz-ri*dbdz*dphiinddp)*deni*dedb*tp_inv*qpart-&
      (g*dbdt*dpsidp-ri*dbdz*dpsidp-g*dbdp*dpsidt+ri*dbdp*dpsidz)*deni*dedb*(kappa-qpart*tp_inv*phipEq)+&
      (g*dbdt-ri*dbdz)*dzonal*deni*dedb + apart*upara*upara*tp_inv*deni*(gp*dphiinddt-rip*dphiinddz) &
         +apart*upara*upara*deni*kappa*cinv*(gp*dpsidt-rip*dpsidz)
      wdrift=(g*(dphieffdt*tp_inv*qpart+kappa*dpsidt-qpart*phipEq*tp_inv*dpsidt)-ri*(dphieffdz*tp_inv*qpart+kappa*dpsidz-qpart*phipEq*tp_inv*dpsidz))*deni*(dzonal+phipEq)
    endif

    wdrift0=(g*dbdp*dptdt-g*dbdt*dptdp-ri*dbdp*dptdz+ri*dbdz*dptdp)*gqi*tp_inv*apart*upara*upara0*b_inv &
            - apart*upara*upara0*tp_inv*deni*(gp*dptdt-rip*dptdz)
    wpara0=-upara0*qpart*tp_inv*epara

    wdot=(zpart(7,m)-paranl*(zpart(5,m)+deltaf0))*(wdrive+wpara+wdrift+wdriftind+wdrift0+wpara0)
#else

    ! wdrive: the equilibrium gradients (density, temperature, flow) dot radial drift velocity
    ! The neoclassical drive (gradB and curvature drift induced) is not included, since we normally assume a neoclassical equilibrium, so these terms have been taken into account in the zeroth order equation.
    ! The magnetic flutter term will be included in wisland.
    wdrive=(kappa+kappav)*(vdr+vbpara)
    ! wpara: the parallel electric field accelaration term dot velocity space
    ! gadient. 
    ! wpara: the total parallel electric field interaction term, including the
    ! perturbed magnetic field induced radial interaction term
    ! In wpara: The epara_na is non-adiabatic parallel electric field in conservative scheme, and epara_na is zero when CS method is not applied.
    wpara=(epara+(epara_na+perturb*epara_pf) + b*deni*phipEq*(ri*dlamdz-g*dlamdt))*(upara-upara0)*dist_term
    ! wdrift: radial drift velocity interaction term
    wdrift=(g*(dbdt*dptdp-dbdp*dptdt)+ri*(dbdp*dptdz-dbdz*dptdp))*deni*dedb0*dist_term &
           +apart*(upara*upara-upara*upara0)*tp_inv*deni*(gp*dptdt-rip*dptdz)
    ! wisland contains magnetic island induced field line drift effects.
    ! both convective term and acceleration term are included
    ! Both analytic and equilibrium islands are included.
    wisland=(kappa+kappav)*vap +&
            dist_term*(upara-upara0)*b*deni*(phipEq+dptdp)*(ri*dlamidz_eq-g*dlamidt_eq)
    if(deltabpara==1)then
      !mirror force term--Eq(26),Dong POP 17
      wpara=wpara+ebpara*(upara-upara0)*dist_term
      !v_d dot nabla B_para term--no curvature part
      wdrift=wdrift+(g*(dbdt*dbparadp-dbdp*dbparadt)+ri*(dbdp*dbparadz-dbdz*dbparadp))*deni*dedb0*dist_term &
                   +apart*(upara*upara-upara*upara0)*tp_inv*deni*(gp*dbparadt-rip*dbparadz)
    endif

    wdot=(zpart(7,m)-paranl*zpart(5,m))*(wdrive+wpara+wdrift+wisland)
#endif

    ! Particle orbit calculation
    ! only self-consistent fields are controlled by the non-linear flag
    dptdp=dptdp*perturb + phipEq
    dptdt=dptdt*perturb
    dptdz=dptdz*perturb
    dlamdp=dlamdp*perturb + dlamidp_eq
    dlamdt=dlamdt*perturb + dlamidt_eq
    dlamdz=dlamdz*perturb + dlamidz_eq
    lam=lam*perturb + lami_eq
    plampt=plampt*perturb
    dbparadp=dbparadp*perturb
    dbparadt=dbparadt*perturb
    dbparadz=dbparadz*perturb

    ! particle velocity
    ! rho_para dot partially checked by L.Shi, 03/22/2017
    ! dB/dzeta term due to non-axisymmetric equilibrium not checked.
    ! paranl is added to necessary terms
!==================================================================================================================
 !   rdot = ((gp*(zpart(4,m)+lam)+g*dlamdp-1.0)*(dedb*dbdt+paranl*(dptdt+dbparadt))-&
 !     (q+rip*(zpart(4,m)+lam)+ri*dlamdp)*(dedb*dbdz+paranl*(dptdz+dbparadz))+&
 !     (ri*dlamdz-g*dlamdt)*(dedb*dbdp+paranl*(dptdp-phipEq+dbparadp)+phipEq))*deni-paranl*plampt
!==================================================================================================================
    ! New rdot accurately separates the parallel nonlinearity and the perturbed field in convective term
    rdot = ((gp*(zpart(4,m)+paranl*lam)+paranl*g*dlamdp-1.0)*(dedb_perp*dbdt+paranl*(dptdt+dbparadt))-&
      (q+rip*(zpart(4,m)+paranl*lam)+paranl*ri*dlamdp)*(dedb_perp*dbdz+paranl*(dptdz+dbparadz))+&
      paranl*(ri*dlamdz-g*dlamdt)*(dedb_perp*dbdp+paranl*(dptdp-phipEq+dbparadp)+phipEq))*deni-paranl*plampt+&
      ((gp*(zpart(4,m)+lam)+g*dlamdp-1.0)*dedb_para*dbdt-&
      (q+rip*(zpart(4,m)+lam)+ri*dlamdp)*dedb_para*dbdz+&
      (ri*dlamdz-g*dlamdt)*dedb_para*dbdp+&
      (1.0-paranl)*zpart(4,m)*b_inv*(dptdp*(ri*dbdz-g*dbdt)+(g*dptdt-ri*dptdz)*dbdp))*deni     

    ! psi dot partially checked by L.Shi, 03/22/2017.
    ! dB/dzeta term due to non-axisymmetric equilibrium not checked.
    ! dbpara terms not checked.
    ! vdrtmp is an artificial flux canceling term that suppresses the local
    ! radial flux. vdrtmp is calculated when irestore=1. 
    pdot = (-g*dedb*dbdt +ri*dedb*dbdz-g*(dptdt+dbparadt)+ri*(dptdz+dbparadz)+upara*b*(g*dlamdt-ri*dlamdz))*deni &
           -vdrtmp(ii)
    ! theta dot checked by L.Shi, 03/22/2017
    ! dbpara terms not checked.
    tdot = (upara*b*(1.0-gp*(zpart(4,m)+lam)-g*dlamdp)+g*(dedb*dbdp+dptdp+dbparadp))*deni
    ! zeta dot checked by L.Shi, 03/22/2017
    ! dbpara terms not checked.
    zdot = (upara*b*(q+rip*(zpart(4,m)+lam)+ri*dlamdp)-ri*(dedb*dbdp+dptdp+dbparadp))*deni

! eliminate w's time dependence when doing poincare.
! if we don't do this, the w wouldn't be well controlled 
! and has the possibility of overflow
    if(ipoincare>0) wdot=0

    ! update particle position
    if(zpart0(1,m) < paxis)then
      ! particles close to axis use (x,y) coordinates
      pdum=sqrt(zpart0(1,m))
      cost=pdum*cos(zpart0(2,m))
      sint=pdum*sin(zpart0(2,m))
      pdum=1.0/zpart(1,m)
      xdot   = 0.5*pdot*cost*pdum-sint*tdot
      ydot   = 0.5*pdot*sint*pdum+cost*tdot
      xdum   = cost + dtime*xdot
      ydum   = sint + dtime*ydot
      zpart(1,m) = max(1.0e-8_lk*psi1,xdum*xdum+ydum*ydum)
      zpart(2,m) = sign(1.0_lk,ydum)*acos(max(-1.0_lk,min(1.0_lk,xdum/sqrt(zpart(1,m)))))
    else
      zpart(1,m) = max(1.0e-8_lk*psi1,zpart0(1,m)+dtime*pdot)
      zpart(2,m) = zpart0(2,m)+dtime*tdot
    endif

    zpart(3,m) = zpart0(3,m)+dtime*zdot
    zpart(4,m) = zpart0(4,m)+dtime*rdot
    zpart(5,m) = zpart0(5,m)+dtime*wdot
    zpart(2,m)=modulo(zpart(2,m),pi2)
    zpart(3,m)=modulo(zpart(3,m),torbound)

    ! store GC information for flux measurements
    wpgc(1,m)=vdr
    wpgc(2,m)=energy
    wpgc(3,m)=upara*majorr
    ! store magnetic flutter velocity in island case. L.Shi,04/17/2017
    ! use wbgc(1,m) as temporary storage since wbgc array is not used after
    ! this line
    wbgc(1,m)=vap
  enddo
  !$acc end parallel

  if(pload>99) then
     !$acc parallel loop gang vector
     do m=1,mp
        zpart(5,m)=1.0 ! for full-f simulation
     enddo
     !$acc end parallel
  endif

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

    ! out of boundary particle
    nb=mthetamax/2 !B array index limit, see bminmax in setup.F90
#ifdef _OPENACC
    !$acc parallel loop gang vector private(dx(spdim))
#else
    !$omp parallel do private(m,tdum,zetatmp,isp,jst,ksz,dpx,dp2,dt2,dzx,dx,dtx,bdum,b_index,ii)
#endif
    do m=1,mp
      if(zpart(1,m) > psi1)then
        zpart(1,m)=zpart0(1,m)
        zpart(3,m)=zpart0(3,m)
        zpart(4,m)=zpart0(4,m)
#ifdef _TOROIDAL3D
!!!!! outer boundary condition for stellarator cases
!!!!! updated by HONGYU WANG, JIAN BAO, April, 2018
        zpart(2,m)=pi2-zpart0(2,m)
        zpart(5,m)=0.0_lk
#else
        zpart(5,m)=0.0_lk!zpart0(5,m)
        if(numereq==0)then
          zpart(2,m)=2.0*pi-zpart0(2,m) ! analytic equilibrium assumes up-down symmetry
        else
          tdum=zpart0(2,m)
          isp=max(1,min(lsp-1,ceiling(psi1*spdpsi_inv)))
          jst=max(1,min(lst-1,ceiling(tdum*spdtheta_inv)))
          dpx = psi1-spdpsi*real(isp-1)
          if(isp == 1) dpx = sqrt(dpx)
          dp2 = dpx*dpx
          dtx=tdum-spdtheta*real(jst-1)
          dt2 = dtx*dtx

         ! b-field at leaving point
! 2D and 3D equilibrium spline arrays have different dimensions
! so they need to be treated differently.

! 2D equilibrium coordinate splines

          dx(1)=1.0_lk
          dx(2)=dpx
          dx(3)=dp2
          dx(4:6)=dx(1:3)*dtx
          dx(7:9)=dx(1:3)*dt2
          bdum=0.0_lk
          do ii = 1, 9 
            bdum=bdum +bsp(ii,isp,jst)*dx(ii)
          enddo
          b_index=floor((bdum-minbfield(1))/deltab(1))+1 !index of b-field bin
          ! bdum could be smaller than minbfield or larger than maxbfield due
          ! to low resolution in theta grids
          ! we need to keep ``b_index`` in the correct range
          ! Note that the max index in thetaupp and thetadwn is mthetamax/2
          b_index = min(nb, max(1,b_index))
          tdum=mod(pi2+tdum-thetabmin(1),pi2) ! poloidal angle w.r.t. min b-field angle
          if(tdum<thetabmax(1))then ! return in the  opposite half-plane
             zpart(2,m)=mod(deltat(mpsi)*(real(thetaupp(1,b_index))+0.5),pi2)
          else
             zpart(2,m)=mod(deltat(mpsi)*(real(thetadwn(1,b_index))+0.5),pi2)
          endif
        endif
#endif
      elseif(zpart(1,m) < psi0)then
        zpart(1,m)=zpart0(1,m)
        zpart(3,m)=zpart0(3,m)
        zpart(4,m)=zpart0(4,m)
#ifdef _TOROIDAL3D
!!!!! inner boundary condition for stellarator cases
!!!!! updated by HONGYU WANG, JIAN BAO, April, 2018
        zpart(2,m)=pi2-zpart0(2,m)
        zpart(5,m)=0.0_lk 
#else
        zpart(5,m)=zpart0(5,m)
        if(numereq==0)then
          zpart(2,m)=2.0*pi-zpart0(2,m) ! analytic equilibrium assumes up-down symmetry
        else
          tdum=zpart0(2,m)
          isp=max(1,min(lsp-1,ceiling(psi0*spdpsi_inv)))
          jst=max(1,min(lst-1,ceiling(tdum*spdtheta_inv)))
          dpx = psi0-spdpsi*real(isp-1)
          if(isp == 1) dpx = sqrt(dpx)
          dp2 = dpx*dpx
          dtx=tdum-spdtheta*real(jst-1)
          dt2 = dtx*dtx

        ! b-field at leaving point
! 2D and 3D equilibrium spline arrays have different dimensions
! so they need to be treated differently.

! 2D equilibrium coordinate splines

          dx(1)=1.0_lk
          dx(2)=dpx
          dx(3)=dp2
          dx(4:6)=dx(1:3)*dtx
          dx(7:9)=dx(1:3)*dt2
          bdum=0.0_lk
          do ii = 1, 9 
            bdum=bdum +bsp(ii,isp,jst)*dx(ii)
          enddo
          b_index=floor((bdum-minbfield(0))/deltab(0))+1
          ! bdum could be smaller than minbfield or larger than maxbfield due
          ! to low resolution in theta grids
          ! we need to keep ``b_index`` in the correct range
          ! Note that the max index in thetaupp and thetadwn is mthetamax/2
          b_index = min(nb, max(1, b_index))
          tdum=mod(pi2+tdum-thetabmin(0),pi2)
          if(tdum<thetabmax(0))then ! return in the  opposite half-plane
             zpart(2,m)=mod(deltat(0)*(real(thetaupp(0,b_index))+0.5),pi2)
          else
             zpart(2,m)=mod(deltat(0)*(real(thetadwn(0,b_index))+0.5),pi2)
          endif
        endif
#endif

      endif
    enddo
    !$acc end parallel

    if(irestore==1)then
    ! Restore temperature profile when running a nonlinear calculation
#ifdef _hybrid
      if(nonlinear==1 .and. pload==1 .and. irk==2 .and. isNhybrid .and. icycle==2*ncyclep)then
#else
      if((nonlinear==1 .and. pload==1 .and. irk==2) .and. (.not. subcycle .or. (subcycle .and. &
        icycle==2*ncyclep)))then
#endif

#ifdef _OPENACC
       !$acc parallel loop gang vector
#else
      !$omp parallel do private(m,psitmp,isp,dpx,rg,upara0)
#endif
        do m=1,mp
          psitmp=zpart(1,m)
          isp=max(1,min(lsp-1,ceiling(psitmp*spdpsi_inv)))
          dpx=psitmp-spdpsi*real(isp-1)
          if(isp==1)dpx=sqrt(dpx)
          rg=  rgpsi(1,isp)  +rgpsi(2,isp)*dpx  +rgpsi(3,isp)*dpx*dpx
          iir(m)=max(0,min(mpsi,int((rg-rg0)*delr+0.5)))
        enddo
        !$acc end parallel

        !$acc parallel loop gang vector
        do ii=0,mpsi
          dtem(ii)=0.0_lk
          dden(ii)=0.0_lk
          fden(ii)=0.0_lk
          dmark(ii)=0.0_lk
        enddo
        !$acc end parallel

      !Todo: why cannot use reduction with openacc
      !javierhn: Because in OPENACC 2.5 reduction doesnt work for arrays
      !          Use atomic for now
#ifdef _OPENACC
      !$acc parallel loop gang vector
#else
      !$omp parallel do private(m,ii,fullf,deltaf) reduction(+:dtem,dmark,dden,fden)
#endif
        do m=1,mp
          ii=iir(m)
          fullf=zpart(7,m)
          deltaf=zpart(5,m)
          !$acc atomic update
          dtem(ii)=dtem(ii)+wpgc(2,m)*deltaf
          !$acc atomic update
          dmark(ii)=dmark(ii)+(wpgc(1,m)+wbgc(1,m))*fullf
          !$acc atomic update
          dden(ii)=dden(ii)+deltaf
          !$acc atomic update
          fden(ii)=fden(ii)+fullf
        enddo
        !$acc end parallel

        !$acc update host(dden,dmark,dtem,fden)
        icount=mpsi+1
        ddum=0.0
        call MPI_ALLREDUCE(dtem,ddum,icount,mpi_Rsize,MPI_SUM,MPI_COMM_WORLD,ierror)
        dtem=ddum
        call MPI_ALLREDUCE(dmark,ddum,icount,mpi_Rsize,MPI_SUM,MPI_COMM_WORLD,ierror)
        dmark=ddum
        call MPI_ALLREDUCE(dden,ddum,icount,mpi_Rsize,MPI_SUM,MPI_COMM_WORLD,ierror)
        dden=ddum
        call MPI_ALLREDUCE(fden,ddum,icount,mpi_Rsize,MPI_SUM,MPI_COMM_WORLD,ierror)
        fden=ddum

        irsmooth=int(sqrt(real(mpsi)))
        dmark=dmark/max(1.0_lk,fden) !radial marker flux
        do i=1,irsmooth
          rdum=dmark(1)
          tdum=dmark(mpsi-1)
          dmark(1:mpsi-1)=0.5*dmark(1:mpsi-1)+0.25*(dmark(0:mpsi-2)+dmark(2:mpsi))
          dmark(0)=0.5*(dmark(0)+rdum)
          dmark(mpsi)=0.5*(dmark(mpsi)+tdum)
        enddo
        tdum=0.1
        pflux=(1.0-tdum)*pflux+tdum*dmark

      ! remove small scale temperature perturbation
        irsmooth=mpsi
        dtem=dtem/(mesht*max(1.0_lk,fden))
        do i=1,irsmooth
          rdum=dtem(1)
          tdum=dtem(mpsi-1)
          dtem(1:mpsi-1)=0.5*dtem(1:mpsi-1)+0.25*(dtem(0:mpsi-2)+dtem(2:mpsi))
          dtem(0)=0.5*(dtem(0)+rdum)
          dtem(mpsi)=0.5*(dtem(mpsi)+tdum)
        enddo
        tdum=0.01
        rdtem=(1.0-tdum)*rdtem+tdum*dtem

      ! remove small scale density perturbation
        dden=dden/(max(1.0_lk,fden))
        do i=1,irsmooth
          rdum=dden(1)
          tdum=dden(mpsi-1)
          dden(1:mpsi-1)=0.5*dden(1:mpsi-1)+0.25*(dden(0:mpsi-2)+dden(2:mpsi))
          dden(0)=0.5*(dden(0)+rdum)
          dden(mpsi)=0.5*(dden(mpsi)+tdum)
        enddo
        tdum=0.01
        rdden=(1.0-tdum)*rdden+tdum*dden

      !$acc update device(rdtem,rdden)
#ifdef _OPENACC
      !$acc parallel loop gang vector
#else
      !$omp parallel do private(m,ii)
#endif
        do m=1,mp
          ii=iir(m)
          zpart(5,m)=zpart(5,m)-(wpgc(2,m)/mesht(ii)-1.5)*rdtem(ii)-rdden(ii)
        end do
        !$acc end parallel
      endif
    endif

    ! apply radial decay boundary condition to all the particle weights
    if (BC_weight==1) call weight_bound(zpart, mp)
  endif

#ifdef _hybrid
  if(idiag == 0 .and. isIhybridOne*icycle*irksub == 1)then
#else
  if(idiag==0 .and. (.not. subcycle .or. subcycle .and. icycle*irksub==1))then
#endif
    ! fluxes diagnose at irk=1
#ifdef _OPENACC
    !$acc parallel loop gang vector
#else
    !$omp parallel do private(m)
#endif
    do m=1,mpdiag
      diagpart(m)=0.0_lk
    enddo
    !$acc end parallel

#ifdef _OPENACC
    !$acc parallel loop gang vector
#else
    !$omp parallel do private(m)
#endif
    do m=0,mpsi
      dden(m)=0.0_lk
    enddo
    !$acc end parallel

#ifdef _OPENACC
    !$acc parallel loop gang vector
#else
    !$omp parallel do private(i,m)
#endif
    do i=1,mpdata1d
      do m=0,mpsi
        data1d(m,i)=0.0_lk
      enddo
    enddo
    !$acc end parallel

    iout = 0 !initialize per time step
    !$acc parallel loop gang vector private(dx(spdim))
    do m=1,mp
      psitmp=zpart0(1,m)
      isp=max(1,min(lsp-1,ceiling(psitmp*spdpsi_inv)))
      dpx=psitmp-spdpsi*real(isp-1)
      dp2=dpx*dpx

      upara0=ropp(1,isp)  +ropp(2,isp)*dpx   +ropp(3,isp)*dp2

      if(isp==1)dpx=sqrt(dpx)
      dp2=dpx*dpx
      rg= rgpsi(1,isp)  +rgpsi(2,isp)*dpx  +rgpsi(3,isp)*dp2

 ! poloidal spline is regular
      thetatmp=zpart0(2,m)
      jst=max(1,min(lst-1,ceiling(thetatmp*spdtheta_inv)))
      dtx=thetatmp-spdtheta*real(jst-1)
      dt2=dtx*dtx
      dzx=zpart0(3,m)-zeta0

! 2D and 3D equilibrium spline arrays have different dimensions
! so they need to be treated differently.

! 2D equilibrium coordinate splines
#ifndef _TOROIDAL3D
      dx(1)=1.0_lk
      dx(2)=dpx
      dx(3)=dp2
      dx(4:6)=dx(1:3)*dtx
      dx(7:9)=dx(1:3)*dt2
      majorr=0.0_lk
      do ii = 1, 9
        majorr=majorr +xsp(ii,isp,jst)*dx(ii)
      enddo
#else
      zetatmp = zpart0(3,m)
      ! We first calculate the spline indexes. The psi and theta indexesare the same
      ! as calculated from other equilibrium splines.
      ! dzeta needs to be calculated again because of sub-section used in 3D equilibrium
      ksz=max(1,min(nzsp_sec,ceiling((zetatmp-zeta0)*spdzeta_inv)))
      dzx = zetatmp - (zeta0 + real(spdzeta*(ksz-1)))

      ! Now we evaluate majorr using spline coefficients
      dx(1)=1.0_lk
      dx(2)=dpx
      dx(3)=dp2
      dx(4:6)=dx(1:3)*dtx
      dx(7:9)=dx(1:3)*dt2
      dx(10:18)=dx(1:9)*dzx
      dx(19:27)=dx(10:18)*dzx
      majorr = 0.0_lk
      do ii=1,27
        majorr = majorr + xsp(ii,isp,jst,ksz)*dx(ii)
      enddo
#endif

      upara0=upara0*majorr
      !! r= rpsi(1,isp)  +rpsi(2,isp)*dpx  +rpsi(3,isp)*dp2
      !! q=    qpsi(1,isp)   +qpsi(2,isp)*dpx   +qpsi(3,isp)*dp2

      ii=max(1,min(mpsi,ceiling((rg-rg0)*delr)))    !radial grid on inner flux surface
      dp1=(psimesh(ii)-psitmp)*delp(ii) !weight for inner flux surface

#ifdef _FRC
#ifdef _hybrid
      if(solvermethod==0)then
        if(ii<1)then
          dp1=0.0 !all weight on inner flux surface between [0,1]
        else
          dp1=1.0 !all weight on outer flux surface between [1,2]
        endif
      endif
#endif
#endif
      tp_inv=1.0/(dp1*mesht(ii-1)+(1.0-dp1)*mesht(ii)+0.5*apart*upara0*upara0)

      fullf=zpart(7,m)
      deltaf=zpart0(5,m)
      energy=wpgc(2,m)*tp_inv-1.5
      vdr=wpgc(1,m)!!!*q/r
      angmom=wpgc(3,m)
      vap=wbgc(1,m) ! magnetic flutter velocity

      ! radial profile of particle and energy flux
      !$acc atomic update
      dden(ii-1)=dden(ii-1)+fullf*dp1
      !$acc atomic update
      dden(ii)  =dden(ii)+  fullf*(1.0-dp1)

      ! Use full radial velocity including EXB and magnetic fluttervelocities
      ! to calculate the total particle flux.
      ! L. Shi, 04/17/2017
      !$acc atomic update
      data1d(ii-1,1)=data1d(ii-1,1)+(vdr+vap)*deltaf*dp1
      !$acc atomic update
      data1d(ii,  1)=data1d(ii,  1)+(vdr+vap)*deltaf*(1.0-dp1)

      ! use full radial velocity to calculate the total energy and toroidal
      ! momentum fluxes. L. Shi, 05/02/2017
      !$acc atomic update
      data1d(ii-1,2)=data1d(ii-1,2)+(vdr+vap)*deltaf*energy*dp1
      !$acc atomic update
      data1d(ii,  2)=data1d(ii,  2)+(vdr+vap)*deltaf*energy*(1.0-dp1)
      ! radial profiles of momentum flux
      !$acc atomic update
      data1d(ii-1,3)=data1d(ii-1,3)+(vdr+vap)*deltaf*angmom*dp1
      !$acc atomic update
      data1d(ii,3)=data1d(ii,3)+(vdr+vap)*deltaf*angmom*(1.0-dp1)

      !volume averaged: density,entropy,flow,energy,fluxes of particle,momentum,heat
      !$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
      !$acc atomic update
      diagpart(4)=diagpart(4)+angmom*deltaf
      !$acc atomic update
      diagpart(5)=diagpart(5)+energy
      !$acc atomic update
      diagpart(6)=diagpart(6)+energy*deltaf
      ! use EXB and magnetic flutter velocities to calculate total particle
      ! flux. L. Shi 04/17/2017
      !$acc atomic update
      diagpart(7)=diagpart(7)+(vdr+vap)*deltaf
      ! use total radial velocities to calculate toroidal momentum flux and
      ! energy flux. L. Shi 05/02/2017
      !$acc atomic update
      diagpart(8)=diagpart(8)+(vdr+vap)*angmom*deltaf
      !$acc atomic update
      diagpart(9)=diagpart(9)+(vdr+vap)*energy*deltaf

      ! tally out of bound particles
      if(zpart(1,m) > psi1 .or. zpart(1,m) < psi0) iout = iout + 1
    enddo
    !$acc end parallel
    !$acc update host(dden,data1d,diagpart)
    diagpart(10)=real(iout)
    diagpart(11)=real(mp)

    ! sum over all MPI processes
    icount=mpdiag
    diagtmp=0.0
    call MPI_REDUCE(diagpart,diagtmp,icount,mpi_Rsize,MPI_SUM,0,MPI_COMM_WORLD,ierror)
    diagpart=diagtmp

    icount=(mpsi+1)*mpdata1d
    data1dtmp=0.0
    call MPI_ALLREDUCE(data1d,data1dtmp,icount,mpi_Rsize,MPI_SUM,MPI_COMM_WORLD,ierror)

    icount=mpsi+1
    call MPI_ALLREDUCE(dden,ddum,icount,mpi_Rsize,MPI_SUM,MPI_COMM_WORLD,ierror)

    ! radial profile data normalized by marker #
    !$omp parallel do private(i,m)
    do i=1,mpdata1d
      do m=0,mpsi
        data1d(m,i)=data1dtmp(m,i)/max(1.0,ddum(m))
      enddo
    enddo
  endif
end subroutine pushParticle
