! 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.
! ==============================================================================

subroutine pushfield
  use system_env
  use precision, only: lk, machineEpsilon, mpi_Rsize
  use global_parameters
  use field_array
  use particle_array
  use equilibrium
  use conservative_scheme,only:apara_na0,apara_na,tapara_na,sapara_na,gradapara_na,gradphiind
  use antenna, only: grad_phi_ext
  implicit none

  integer i,j,ij,isp,jst,iMHDxi_output
  real(lk) dtime,pdum,tdum,dpx,dp2,dtx,dt2,dbdt,dbdp,dbdz,&
           ri,g,b_inv,q_inv,&
           dptdp,dptdt,dptdz,wpara,wdrive,wstar,wdrift,wdamp,&
           ne,te,dnedp,dtedp,phipEq,wdrift0,web0,ddnedp,ddnedt,ddnedz,&
           ddpdp,ddpdt,ddpdz,dinddp,dinddt,dinddz,dpsidp,dpsidt,dpsidz,dne,&
           daparadp,daparadt,daparadz,dsdp,dphieffdp,dphieffdt,dphieffdz,&
           dkpparadp,dkpparadt,dkpparadz,dkpperpdp,dkpperpdt,dkpperpdz,dkne,dkpe,&
           dbparadp,dbparadt,dbparadz,bpara0,&
           dqdp,dgdp,d2gdp2,dridp,d2ridp2,weqc,rbc(0:mpsi),&
           dudp,dudt,dudz,wstarnl,wperpnl,perturb,&
           gqi_inv,deltapsi00(0:mpsi),apara00_a(0:mpsi),apara00_na(0:mpsi),apara00nl(0:mpsi),dnedt00(0:mpsi),&
           upara0,dupara0dp,majorr,tmp_field(mgrid),tmp_field2(0:1,mgrid),rsfluidne(0:mpsi),rstmp
  !$acc declare create(rbc, apara00_a, apara00_na, apara00nl, dnedt00, deltapsi00,tmp_field,rsfluidne)

  perturb=real(nonlinear)
  rstmp=0.0_lk
  !radial BC
!$acc kernels
  rbc(0:mpsi)=1.0
  rsfluidne(0:mpsi)=0.0_lk
!$acc end kernels

!  do i=0,nbound-1
!     rbc(i)=real(i)/real(nbound)
!  enddo

  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(ij)
#endif
    do ij=1,mgrid
        apara0(1,ij)=apara(1,ij)
        apara_na0(1,ij)=apara_na(1,ij)
        fluidne0(1,ij)=fluidne(1,ij)
        deltapsi0(1,ij)=deltapsi(1,ij)
    enddo

#ifdef _OPENACC
    !$acc parallel loop gang vector
#else
    !$omp parallel do private(ij)
#endif
    do ij=1,mgrid
        aparanl0(1,ij)=aparanl(1,ij) ! Total nonlinear apara
    enddo

! 2nd step of Runge-Kutta method
  else
    dtime=tstep
  endif

! advance del_ne and a_parallel
#ifdef _OPENACC
  !$acc parallel loop gang
#else
  !$omp parallel do private(i,j,ij,pdum,isp,dpx,dp2,g,ri,tdum,jst,dtx,dt2,&
  !$omp& b_inv,q_inv,phipEq,dbdt,dbdp,dbdz,ne,te,dnedp,dtedp,dbparadp,&
  !$omp& wpara,wstar,wdrift,wdrive,wdrift0,web0,dne,wdamp,dbparadt,dbparadz,&
  !$omp& ddnedp,ddnedt,ddnedz,ddpdp,ddpdt,ddpdz,dptdp,dptdt,dptdz,dinddp,dinddt,dinddz,&
  !$omp& daparadp,daparadt,daparadz,dpsidp,dpsidt,dpsidz,dphieffdp,dphieffdt,dphieffdz,&
  !$omp& dkpparadp,dkpparadt,dkpparadz,dkpperpdp,dkpperpdt,dkpperpdz,dkne,dkpe,bpara0,&
  !$omp& dqdp,dgdp,dridp,d2gdp2,d2ridp2,dsdp,weqc,dudp,dudt,dudz,wstarnl,wperpnl,&
  !$omp& upara0,dupara0dp,majorr,&
  !$omp& gqi_inv)
#endif
  do i=0,mpsi
     pdum=psimesh(i)
     q_inv=1/qmesh(i)

     isp=max(1,min(lsp-1,ceiling(pdum*spdpsi_inv)))
     dpx=pdum-spdpsi*real(isp-1)
     dp2=dpx*dpx

     g=   gpsi(1,isp)   +gpsi(2,isp)*dpx   +gpsi(3,isp)*dp2
     ri=  cpsi(1,isp)   +cpsi(2,isp)*dpx   +cpsi(3,isp)*dp2

     dqdp=               qpsi(2,isp)       +2.0*qpsi(3,isp)*dpx
     dgdp=               gpsi(2,isp)       +2.0*gpsi(3,isp)*dpx
     dridp=              cpsi(2,isp)       +2.0*cpsi(3,isp)*dpx

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

     gqi_inv=1.0/(g*qmesh(i)+ri)

! 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 introduced
     d2gdp2=2.0*(dpx*spdpsi_inv*(gpsi(3,isp+1)-gpsi(3,isp))+gpsi(3,isp))
     d2ridp2=2.0*(dpx*spdpsi_inv*(cpsi(3,isp+1)-cpsi(3,isp))+cpsi(3,isp))

     dsdp=(g*d2ridp2-ri*d2gdp2)/(g*qmesh(i)+ri) &
        -(g*dridp-ri*dgdp)*(dgdp*qmesh(i)+g*dqdp+dridp)*((gqi_inv)**2)

     ne=meshne(i)
     te=meshte(i)
     dnedp=-ne*kapane(i)
     dtedp=-te*kapate(i)
     phipEq=-mesher(i)! mesher = -dphi/dpsi

     ! Now, for B field spline, when isp==1, dpsi is changed to dpx=sqrt(psi)
     if(isp==1) then
        dpx=sqrt(dpx)
        dp2=dpx*dpx
     endif
     !$acc loop vector
     do j=1,mtheta(i)
        ij=igrid(i)+j
        b_inv=1.0_lk/bmesh(ij)

        if(fielddir==1 .or. fielddir==3)then
          tdum=modulo(deltat(i)*real(j)+(zeta1-pi2)*qtinv(i),pi2)
        else
          tdum=modulo(deltat(i)*real(j)+zeta1*qtinv(i),pi2)
        endif
        jst=max(1,min(lst-1,ceiling(tdum*spdtheta_inv)))
        dtx=tdum-spdtheta*real(jst-1)
        dt2=dtx*dtx
        dbdp = dbmesh(1,ij)
        dbdt = dbmesh(2,ij)
        dbdz = dbmesh(3,ij)
 
        majorr= xsp(1,isp,jst)+xsp(2,isp,jst)*dpx+xsp(3,isp,jst)*dp2+ &
             (xsp(4,isp,jst)+xsp(5,isp,jst)*dpx+xsp(6,isp,jst)*dp2)*dtx+ &
             (xsp(7,isp,jst)+xsp(8,isp,jst)*dpx+xsp(9,isp,jst)*dp2)*dt2

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

        dptdp=gradphi(1,1,ij) ! contain both zonal and non-zonal
        dptdt=gradphi(2,1,ij) ! contain only non-zonal
        dptdz=gradphi(3,1,ij)-gradphi(2,1,ij)*q_inv  ! contain only non-zonal
        dbparadp=0.0
        dbparadt=0.0
        dbparadz=0.0
        bpara0=0.0
        if(deltabpara==1)then
         dbparadp=gradbpara(1,1,ij)
         dbparadt=gradbpara(2,1,ij)
         dbparadz=gradbpara(3,1,ij)-gradbpara(2,1,ij)*q_inv
         bpara0=bpara(1,ij)
        endif
        ddnedp=gradne_nz(1,1,ij)
        ddnedt=gradne_nz(2,1,ij)
        ddnedz=gradne_nz(3,1,ij)-gradne_nz(2,1,ij)*q_inv
        daparadp=gradapara(1,1,ij)
        daparadt=gradapara(2,1,ij)
        daparadz=gradapara(3,1,ij)-gradapara(2,1,ij)*q_inv
        dpsidp=gradpsi(1,1,ij)
        dpsidt=gradpsi(2,1,ij)
        dpsidz=gradpsi(3,1,ij)-gradpsi(2,1,ij)*q_inv

        dphieffdp=gradphieff(1,1,ij)
        dphieffdt=gradphieff(2,1,ij)
        dphieffdz=gradphieff(3,1,ij)-gradphieff(2,1,ij)*q_inv

        if(cs_method==0)then
           dinddp=dphieffdp-dptdp  ! dptdp contains both zonal and non-zonal component (only is used for diagnostic)
           dinddt=dphieffdt-dptdt
           dinddz=dphieffdz-dptdz
        elseif(cs_method==1)then
           dinddp=gradphiind(1,1,ij) ! dinddp only contains non-zonal (only is used for diagnostic)
           dinddt=gradphiind(2,1,ij)
           dinddz=gradphiind(3,1,ij)-gradphiind(2,1,ij)*q_inv
        endif

        dne=sfluidne(1,ij)

        if(eta < machineEpsilon)then
           if(cs_method==0)then
           ! adiabatic pressure (either para or perp because they're equal)
              ddpdp=-qelectron*ne*dphieffdp+(te*dnedp+ne*dtedp+qelectron*ne*phipEq)*dpsidp
              ddpdt=-qelectron*ne*dphieffdt+(te*dnedp+ne*dtedp+qelectron*ne*phipEq)*dpsidt
              ddpdz=-qelectron*ne*dphieffdz+(te*dnedp+ne*dtedp+qelectron*ne*phipEq)*dpsidz
           elseif(cs_method==1)then 
              ddpdp=-qelectron*ne*(dptdp+dinddp)+(te*dnedp+ne*dtedp+qelectron*ne*phipEq)*dpsidp
              ddpdt=-qelectron*ne*(dptdt+dinddt)+(te*dnedp+ne*dtedp+qelectron*ne*phipEq)*dpsidt
              ddpdz=-qelectron*ne*(dptdz+dinddz)+(te*dnedp+ne*dtedp+qelectron*ne*phipEq)*dpsidz
           endif
              
           if(deltabpara==1)then
            ddpdp=ddpdp-2.0_lk*dbparadp*b_inv+dbdp*bpara0*b_inv*b_inv
            ddpdt=ddpdt-2.0_lk*dbparadt*b_inv+dbdt*bpara0*b_inv*b_inv
            ddpdz=ddpdz-2.0_lk*dbparadz*b_inv
           endif
        else
           ! tearing mode
           ddpdp=te*ddnedp+dne*dtedp
           ddpdt=te*ddnedt
           ddpdz=te*ddnedz
        endif

        if(nhybrid>0)then
! parallel non-adiabatic pressure
          dkpparadp=ne*gradpepara(1,1,ij)+pressureepara(1,ij)*dnedp
          dkpparadt=ne*gradpepara(2,1,ij)
          dkpparadz=ne*gradpepara(3,1,ij)-ne*gradpepara(2,1,ij)*q_inv
! perpendicular non-adiabatic pressure
          dkpperpdp=ne*gradpeperp(1,1,ij)+pressureeperp(1,ij)*dnedp
          dkpperpdt=ne*gradpeperp(2,1,ij)
          dkpperpdz=ne*gradpeperp(3,1,ij)-ne*gradpeperp(2,1,ij)*q_inv
          dkne=ne*densitye(1,ij)
          dkpe=ne*pressureeperp(1,ij)
        else
          dkpparadp=0.0
          dkpparadt=0.0
          dkpparadz=0.0
          dkpperpdp=0.0
          dkpperpdt=0.0
          dkpperpdz=0.0
          dkne=0.0
          dkpe=0.0
        endif
        
        dudp=gradue(1,1,ij)
        dudt=gradue(2,1,ij)
        dudz=gradue(3,1,ij)-gradue(2,1,ij)*q_inv

! parallel compression
        wpara=jmesh(ij)*qmesh(i)*gradue(3,1,ij)
! wdrive: perturbed pressure gradient dot grad_B_0, interchange drive
        !wdrive=jmesh(ij)*b_inv*b_inv*b_inv*(ri*dbdp*(ddpdz*2.0_lk+dkpparadz+dkpperpdz) &
        !      -ri*dbdz*(ddpdp*2.0_lk+dkpparadp+dkpperpdp) &
        !      +g*dbdt*(ddpdp*2.0_lk+dkpparadp+dkpperpdp)-g*dbdp*(ddpdt*2.0_lk+dkpparadt+dkpperpdt))
        wdrive=jmesh(ij)*b_inv*b_inv*b_inv &
              *(ri*dbdp*(ddpdz*2.0_lk+real(deltabpara)*dbparadz*b_inv+dkpparadz+dkpperpdz-dbparadz*b_inv+dbdz*bpara0*b_inv*b_inv) &
              -ri*dbdz*(ddpdp*2.0_lk+real(deltabpara)*dbparadp*b_inv+dkpparadp+dkpperpdp-dbparadp*b_inv+dbdp*bpara0*b_inv*b_inv) &
              +g*dbdt*(ddpdp*2.0_lk+real(deltabpara)*dbparadp*b_inv+dkpparadp+dkpperpdp-dbparadp*b_inv+dbdp*bpara0*b_inv*b_inv)&
              -g*dbdp*(ddpdt*2.0_lk+real(deltabpara)*dbparadt*b_inv+dkpparadt+dkpperpdt-dbparadt*b_inv+dbdt*bpara0*b_inv*b_inv))
! wstar: perturbed ExB drift dot grad_n_0, omega_star term
        wstar=(ri*dptdz-g*dptdt)*dnedp*jmesh(ij)*b_inv*b_inv
! wdrift: perturbed ExB drift dot grad_B
        wdrift=2.0_lk*(ne+perturb*dne)*jmesh(ij)*(g*dbdp*dptdt-g*dbdt*dptdp &
                    -ri*dbdp*dptdz+ri*dbdz*dptdp)*b_inv*b_inv*b_inv
     if(deltabpara==1)then 
        wdrift=wdrift-3.0_lk*(ne+perturb*(dne-dkne-bpara0/te+dkpe/te&
       -sdeltapsi(1,ij)*ne*kapate(i)))*jmesh(ij)*(g*dbdp*dbparadt-g*dbdt*dbparadp &
                    -ri*dbdp*dbparadz+ri*dbdz*dbparadp)*b_inv*b_inv*b_inv*b_inv
        wdrift=wdrift+1.0_lk*jmesh(ij)*(g*dbparadt &
                    -ri*dbparadz)*b_inv*b_inv*b_inv*(dnedp*te+dtedp*ne+(ddpdp*1.0+dkpperpdp)*perturb)/(ne*te) &
                    -(ne+perturb*(dne-dkne+dkpe/te-sdeltapsi(1,ij)*ne*kapate(i)))*(dridp*dbparadz-dgdp*dbparadt)*b_inv*gqi_inv
     endif
! wdrift0: equilibrium radial ExB drift dot grad delta_n_e
        wdrift0=-2.0*dne*jmesh(ij)*(g*dbdt-ri*dbdz)*phipEq*b_inv*b_inv*b_inv
! web0: equilibrium radial ExB convection of perturbed density gradient
        web0=jmesh(ij)*(-ri*ddnedz+g*ddnedt)*phipEq*b_inv*b_inv
! weqc: parallel equilibrium current terms (fast ion flow not included yet) 
        weqc=(upara0*b_inv*dnedp-ne*upara0*b_inv**2*dbdp+ne*b_inv*dupara0dp-2.0*real(ecurrent)/betae*rho0*rho0*dsdp) &
             *jmesh(ij)*b_inv*(g*daparadt-ri*daparadz) &
            +ne*upara0*b_inv**2*dbdt*jmesh(ij)*b_inv*g*daparadp & 
            -(dridp*(ddpdz+dkpparadz)-dgdp*(ddpdt+dkpparadt))/(g*qmesh(i)+ri) &
            +perturb*dne*(dridp*dptdz-dgdp*dptdt)/(g*qmesh(i)+ri)
! nonlinear compression along perturbed B-field (\delta B*\nabla(n_0u_\|/B_0))
        wperpnl=-jmesh(ij)*(ri*daparadz*dudp-ri*daparadp*dudz+g*daparadp*dudt-g*daparadt*dudp)*b_inv*perturb
! nonlinear exb convection
        wstarnl=jmesh(ij)*(ri*dptdz*ddnedp-ri*dptdp*ddnedz+g*dptdp*ddnedt-g*dptdt*ddnedp)*b_inv*b_inv*perturb

        if(iload==0)then ! ideal MHD: these two terms cancelled by ions in Poisson Eq.
           wstar=0.0
           wdrift=0.0
        else
           weqc=weqc+ne*(dridp*dptdz-dgdp*dptdt)*gqi_inv
        endif

        wdamp=d4fluidne(ij)
!Dynamic equations are modified by Jian(jianb@uci.edu) for conservative scheme and hybrid model:
!dnedt is the divergence of flux term in continuity equation; fluidne is the same for both conservative scheme and hybrid model.
!In conservative scheme: apara is adiabatic vector potential;apara_na is non-adiabatic vector potential.
!In hybrid model: apara is the linear vector potential in Eq. (12) in Holod et al., Phys. Plasmas, 2009;
!aparanl is the non-linear vector potential advanced by the ponderomotive force
!(initally proposed in Zhixuan's thesis, completed by Jian as term {VI} in Eq. (31) in Bao et al., Phys. Plasmas,2017)

        dnedt(1,ij)=-(wpara+wperpnl+wstarnl+wdrive+web0+wstar+wdrift+wdrift0+weqc-wdamp)
        fluidne(1,ij) = fluidne0(1,ij) + dtime*rbc(i)*dnedt(1,ij)

        apara(1,ij) = apara0(1,ij)+dtime*rbc(i)*(dinddt+qmesh(i)*dinddz)*b_inv*jmesh(ij)+&
                      dtime*jmesh(ij)*qmesh(i)*b_inv*(grad_phi_ext(3,1,ij)*rho0*rho0)     !Normalize grad_phi_ext to GTC units by multiplying by rho0^2

        apara_na(1,ij) = apara_na0(1,ij)+dtime*rbc(i)*tapara_na(1,ij)

        deltapsi(1,ij) = deltapsi0(1,ij)-dtime*rbc(i)*dinddt*q_inv
      
!by Zhixuan: Gauge calculation to get MHD \xi_\perp, not needed for push particles, so no need to do Ronge-Kuta
!only ran when MHD displacement is needed for diagnosis (iMHDxi_output=1)
        iMHDxi_output=0
        if((irk==2) .and. (iMHDxi_output==1))then
            gradgaugef(1,1,ij)=gradgaugef(1,1,ij)+dinddp*dtime
            gradgaugef(2,1,ij)=gradgaugef(2,1,ij)+dinddt*dtime
            gradgaugef(3,1,ij)=gradgaugef(3,1,ij)+dinddz*dtime

!diagnosis only used when the step for snapshot
            if(mod(istep,mstep/msnap)==0 .or. istep*(1-irun)==ndiag)then

                MHDxiperp(1,1,ij)=(ri*gradgaugef(3,1,ij)-g*gradgaugef(2,1,ij))*gqi_inv
                MHDxiperp(2,1,ij)=g*gradgaugef(1,1,ij)*gqi_inv
                MHDxiperp(3,1,ij)=-ri*gradgaugef(1,1,ij)*gqi_inv
                MHDxi_mag(1,ij)=MHDxiperp(1,1,ij)*MHDxiperp(1,1,ij)*gdpp(ij)+&
                                MHDxiperp(2,1,ij)*MHDxiperp(2,1,ij)*gdtt(ij)+&
                                MHDxiperp(3,1,ij)*MHDxiperp(3,1,ij)*gdzz(ij)+&
                                MHDxiperp(1,1,ij)*MHDxiperp(2,1,ij)*gdpt(ij)*2.0
                MHDxi_mag(1,ij)=sqrt(MHDxi_mag(1,ij))
            endif
        endif
       
        if(nonlinear>0)then
           aparanl(1,ij) = aparanl0(1,ij)+dtime*rbc(i)*taparanl(1,ij)
        endif

     enddo
  enddo


  if(irk==2 .and. (iMHDxi_output==1))then
! Zhixuan's diagnostics
    call periodicity_gpu(MHDxi_mag)
    !$acc kernels
    tmp_field2=MHDxiperp(1,:,:)
    !$acc end kernels
    call periodicity_gpu(tmp_field2)
    !$acc kernels
    MHDxiperp(1,:,:)=tmp_field2
    tmp_field2=MHDxiperp(2,:,:) 
    !$acc end kernels
    call periodicity_gpu(tmp_field2)
    !$acc kernels
    MHDxiperp(2,:,:)=tmp_field2
    tmp_field2=MHDxiperp(3,:,:)
    !$acc end kernels
    call periodicity_gpu(tmp_field2)
    !$acc kernels
    MHDxiperp(3,:,:)=tmp_field2
    !$acc end kernels
  endif

!=======   numerical zonal field removement for non-zonal variables =======
  !$acc kernels
  deltapsi00=0.0_lk
  apara00_a=0.0_lk
  !$acc end kernels

! non-tearing parity of perturbed magnetic field in both conservative scheme and hybrid model
  !$acc kernels
  tmp_field = deltapsi(1,:)
  !$acc end kernels
  call zonal_field_gpu(tmp_field,deltapsi00)
! adiabatic A|| in both conservative scheme and hybrid model
  !$acc kernels
  tmp_field=apara(1,:)
  !$acc end kernels
  call zonal_field_gpu(tmp_field,apara00_a)
!Enforce excluding zonal components in non-zonal variable: apara and deltapsi

#ifdef _OPENACC
  !$acc parallel loop gang
#else
  !$omp parallel do private(i,j,ij)
#endif
  do i=0,mpsi
     !$acc loop vector
     do j=1,mtheta(i)
        ij=igrid(i)+j
        apara(1,ij)=apara(1,ij)-apara00_a(i)
        deltapsi(1,ij)=deltapsi(1,ij)-deltapsi00(i)
     enddo
     apara(1,igrid(i))=apara(1,igrid(i)+mtheta(i))
     deltapsi(1,igrid(i))=deltapsi(1,igrid(i)+mtheta(i))
  enddo

! smooth and filter for variables with only non-zonal component:
  !$acc kernels
  sdelapara(:,:)=apara(:,:)
  sdeltapsi(:,:)=deltapsi(:,:)
  !$acc end kernels

  call smooth_gpu(sdelapara)
  call smooth_gpu(sdeltapsi)
  if(nfilter>0)then
    CALL FILTER_GPU(sdelapara)
    CALL FILTER_GPU(sdeltapsi)
  endif

! restore electron density
  if(irestore==1 .and. nonlinear==1 .and. irk==2)then
      !$acc kernels
      rsfluidne(:)=qion*rddeni(:)*meshni(:)+qfast*rddenf(:)*meshnf(:)
      !$acc end kernels
#ifdef _OPENACC
  !$acc parallel loop gang
#else
  !$omp parallel do private(i,j,ij,rstmp)
#endif
      do i=0,mpsi
         rstmp=-rsfluidne(i)/qelectron
         !$acc loop vector
         do j=1,mtheta(i)
            ij=igrid(i)+j
            fluidne(1,ij)=fluidne(1,ij)-rstmp
         enddo
      enddo
  endif

!=======smooth and filter for variables with both zonal and non-zonal components====
  !$acc kernels
  fluidne00=0.0_lk
  apara00_na=0.0_lk
  apara00nl=0.0_lk
  apara00=0.0_lk
  !$acc end kernels

! smooth and filter for variable with both zonal and non-zonal components:
  !$acc kernels
  sfluidne(:,:)=fluidne(:,:)
  sapara_na(:,:)=apara_na(:,:)
  saparanl(:,:)=aparanl(:,:)
  if(cs_method==0)then
     sapara=apara+aparanl
  elseif(cs_method==1)then
     sapara=apara+apara_na
  endif
  !$acc end kernels
  
  call smooth_gpu(sfluidne)
  call smooth_gpu(sapara_na)
  call smooth_gpu(saparanl)
  call smooth_gpu(sapara)

  !$acc kernels
    tmp_field=sfluidne(1,:)
  !$acc end kernels
! sdelfluidne is a non-zonal variable
  call zonal_field_gpu(tmp_field,fluidne00)

  !$acc parallel loop gang 
  do i=0,mpsi
     !$acc loop vector 
     do j=igrid(i),igrid(i)+mtheta(i)
        sdelfluidne(0:1,j)=sfluidne(0:1,j)-fluidne00(i)
     enddo
  enddo

  if(nfilter>0)then
    call filter_gpu(sdelfluidne)
    if(izonal==0)then
       call filter_gpu(sfluidne)
       call filter_gpu(sapara_na)
       call filter_gpu(saparanl)
       call filter_gpu(sapara)
    elseif(izonal>0)then
       call filter_keep_zonal_gpu(sfluidne,fluidne00)
       call filter_keep_zonal_gpu(sapara_na,apara00_na)
       call filter_keep_zonal_gpu(saparanl,apara00nl)
       call filter_keep_zonal_gpu(sapara,apara00)
    endif
  endif   
  
 
!=========================================================================
! The time derivative of fluidne (sdnedt) for pushing kinetic electron weight 
  if(nhybrid>0)then

    !$acc kernels
    sdnedt(:,:)=dnedt(:,:)
    !$acc end kernels
    call smooth_gpu(sdnedt)

    if(nfilter>0)CALL FILTER_GPU(sdnedt)
! Safety Guard: to exclude zonal component in sdnedt for kinetic electron weight equation
    !$acc kernels
       tmp_field=sdnedt(1,:)
    !$acc end kernels
    call zonal_field_gpu(tmp_field,dnedt00)

#ifdef _OPENACC
    !$acc parallel loop gang
#else
    !$omp parallel do private(i,j,ij)
#endif
    do i=0,mpsi
       !$acc loop vector
       do j=1,mtheta(i)
          ij=igrid(i)+j
          sdnedt(1,ij)=sdnedt(1,ij)-dnedt00(i)
       enddo
       sdnedt(1,igrid(i))=sdnedt(1,igrid(i)+mtheta(i))
    enddo
    call periodicity_gpu(sdnedt)
  endif

  !$acc update host(sfluidne,sdeltapsi,fluidne,fluidne0,fluidne00,apara00,sapara_na,&
  !$acc&            saparanl,sapara,sdelfluidne,sdnedt)

end subroutine pushfield
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine continuousbc(farray)
  use global_parameters
  use field_array,only:igrid,mtheta
  implicit none

  integer ::i,j,ij
  real(lk)::farray(0:1,mgrid),ave

  i=1
  ave=0.0_lk
  do j=1,mtheta(i)
    ij=igrid(i)+j
    ave = ave + farray(1,ij)/real(mtheta(i))
  enddo

  i=0
  do j=0,mtheta(i)
    ij=igrid(i)+j
    farray(1,ij)=ave
  enddo

end subroutine continuousbc

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine periodicity_cpu(fgrid)
  use global_parameters, only:mgrid
  use precision, only:lk
#ifdef _OPENACC
  use openacc
#endif
  implicit none
  logical :: tmp
  real(lk) :: fgrid(0:1,mgrid)
#ifdef _OPENACC
  if(acc_is_present(fgrid))then
    tmp = .false.
    !$acc update device(fgrid)
  else
    tmp=.true.
    !$acc enter data copyin(fgrid)
  endif
#endif
  call periodicity(fgrid,.true.)
#ifdef _OPENACC
  if(tmp)then
    !$acc exit data copyout(fgrid)
    !$acc exit data delete(fgrid)
  else
    !$acc update host(fgrid)
  endif
#endif

end subroutine periodicity_cpu
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine periodicity_gpu(fgrid)
  call periodicity(fgrid,.true.)
end subroutine periodicity_gpu
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine periodicity(fgrid,gpu_array)
  use system_env
  use precision, only: lk, mpi_Rsize
  use global_parameters
  use field_array
  implicit none

  integer i,j,ii,jt,icount,idest,isource,isendtag,irecvtag,istatus(MPI_STATUS_SIZE)
  real(lk) sendr(mgrid),recvl(mgrid)
  logical :: gpu_array
  real(lk),dimension(0:1,mgrid) :: fgrid
  !$acc declare create(sendr,recvl)

! toroidal end point, send E to right and receive from left
#ifdef _OPENACC
  !$acc parallel loop gang vector
#else
  !$omp parallel do private(i)
#endif
  do i=1,mgrid
     sendr(i)=fgrid(1,i)
     recvl(i)=0.0
  enddo
  icount=mgrid
  idest=right_pe
  isource=left_pe
  isendtag=myrank_toroidal
  irecvtag=isource

#ifdef CUDA_MPI
    !$acc host_data use_device(sendr,recvl)
#else
    !$acc update host(sendr)
#endif
  call MPI_SENDRECV(sendr,icount,mpi_Rsize,idest,isendtag,&
      recvl,icount,mpi_Rsize,isource,irecvtag,toroidal_comm,istatus,ierror)
#ifdef CUDA_MPI
    !$acc end host_data
#else
    !$acc update device(recvl)
#endif

! unpack end point data for k=0
  if(myrank_toroidal==0)then
#ifndef _OPENACC
!$omp parallel do private(i,ii,jt)
     do i=0,mpsi
        ii=igrid(i)
        jt=mtheta(i)
        fgrid(0,ii+1:ii+jt)=cshift(recvl(ii+1:ii+jt),-itran(i))
     enddo
#else
     !$acc parallel loop gang
     do i=0,mpsi
        ii=igrid(i)
        jt=mtheta(i)
        !$acc loop vector
        do j=1,jt
           fgrid(0,ii+j)=recvl(ii+1+modulo(j-itran(i)-1,jt))
        enddo
     enddo
#endif
  else
#ifdef _OPENACC
     !$acc parallel loop gang vector
#else
     !$omp parallel do private(i)
#endif
     do i=1,mgrid
        fgrid(0,i)=recvl(i)
     enddo
  endif

! poloidal end point
#ifdef _OPENACC
  !$acc parallel loop gang vector
#else
  !$omp parallel do private(i,ii,jt)
#endif
  do i=0,mpsi
     ii=igrid(i)
     jt=mtheta(i)
     fgrid(0:1,ii)=fgrid(0:1,ii+jt)
  enddo

end subroutine periodicity

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
