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

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Wrapper subroutine of smooth for cpu and gpu arrays
subroutine smooth_cpu(farray)
  use global_parameters, only: mgrid
  use precision, only: lk
#ifdef _OPENACC
  use openacc
#endif
  implicit none
  logical :: tmp
  real(lk),dimension(0:1,mgrid) :: farray

#ifdef _OPENACC
  if(acc_is_present(farray))then
    tmp = .false.
    !$acc update device(farray)
  else
    tmp=.true.
    !$acc enter data copyin(farray)
  endif
#endif

  call smooth_t(farray, .true.)

#ifdef _OPENACC
  if(tmp)then
    !$acc exit data copyout(farray)
    !$acc exit data delete(farray)
  else
    !$acc update host(farray)
  endif
#endif
end subroutine smooth_cpu

subroutine smooth_gpu(farray)
  call smooth_t(farray, .true.)
end subroutine smooth_gpu
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


subroutine smooth_t(farray,gpu_array)
  use system_env
  use precision, only: lk, machineEpsilon
  use global_parameters
  use field_array
  use interfaces,only: axisExtrapolate
  use particle_array, only: iload
  use boundary_decay, only: field_bound 
  implicit none

  integer i,j,ii,ij,ip,jt,&
       icount,idest,isource,isendtag,&
       irecvtag,istatus(MPI_STATUS_SIZE),parasmooth
  real(lk) sendl(mgrid),sendr(mgrid),recvl(mgrid),recvr(mgrid),phism(mgrid),&
       ptemp,pleft(mthetamax),pright(mthetamax),pright0(mthetamax),pright1(mthetamax),&
       pright2(mthetamax),pright3(mthetamax),farraytmp(mgrid),&
       farray(0:1,mgrid),sigma
  logical :: gpu_array
  !$acc declare create(sendl, sendr, recvl, recvr, phism, pleft, pright, pright0)
  !$acc declare create(pright1,pright2,pright3,farraytmp)
!  if(iload/=9)then
!    ! Remove n=0, m/=0 mode for nonlinear sims with all n's
!    if(nfilter==0 .and. nonlinear>0) call removn0(farray)
!  endif
  parasmooth=1

#ifdef _OPENACC
  !$acc parallel loop gang vector
#else
  !$omp parallel do private(i)
#endif
  do i=1,mgrid
     farraytmp(i)=farray(1,i)
  enddo

  ! Apply Radial Boundary decay on the field array
  if(iload/=9)then
    if(BC_field == 1) call field_bound(farraytmp)
  endif

  do ip=1,ismooth
#ifdef _FRC
#ifdef _OPENACC
       !$acc parallel loop gang vector
#else
       !$omp parallel do private(i)
#endif
     do i=1,mpsi-1
        farraytmp(igrid(i))=farraytmp(igrid(i)+mtheta(i)) ! poloidal BC
     enddo
     ! poloidal smoothing (-0.0625 0.25 0.625 0.25 -0.0625)

#ifdef _OPENACC
     !$acc kernels
#else
     !$omp parallel do private(i,ii,jt,pright,pright0,pright1,pright2,pright3)
#endif
     do i=1,mpsi-1
        ii=igrid(i)
        jt=mtheta(i)
        pright(1:jt)=farraytmp(ii+1:ii+jt)
        pright0(1:jt)=cshift(pright(1:jt),-1)
        pright1(1:jt)=cshift(pright(1:jt),1)
        pright2(1:jt)=cshift(pright0(1:jt),-1)
        pright3(1:jt)=cshift(pright1(1:jt),1)
        farraytmp(ii+1:ii+jt)=0.625*pright(1:jt)+0.25*(pright0(1:jt)+pright1(1:jt))-&
             0.0625*(pright2(1:jt)+pright3(1:jt))
     enddo
     !$acc end kernels
#else
     if(fem>0.and.eta<machineEpsilon)then
       call smooth_perp_fem(farraytmp)
     else
       call smooth_perp(farraytmp)
       do j=1,parasmooth
          call smooth_para(farraytmp)
       enddo
       if(parasmooth>=2)call smooth_perp(farraytmp)
     endif
     call smooth_para(farraytmp)
#endif
  enddo

#ifdef _OPENACC
  !$acc parallel loop gang vector
#else   
  !$omp parallel do private(i)
#endif  
  do i=1,mgrid
     farray(1,i)=farraytmp(i)
  enddo

  if(bcond==1)call axisExtrapolate(farray(1,:))

  call periodicity_gpu(farray)
end subroutine smooth_t

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

subroutine smooth_perp(phitmp)
  use precision, only: lk
  use global_parameters
  use field_array
  use particle_array,only:iload
  implicit none

  integer i,j,ii,ij,ip,jt
  real(lk) phism(mgrid),ptemp,phitmp(mgrid),&
       farray(0:1,mgrid)

  !$acc declare create(phism, phitmp, farray)
! radial smoothing
#ifdef _OPENACC
    !$acc parallel loop gang vector
#else
    !$omp parallel do private(i)
#endif
     do i=1,mpsi-1
        phitmp(igrid(i))=phitmp(igrid(i)+mtheta(i)) ! poloidal BC
     enddo
  !$acc kernels
     phism=0.0
  !$acc end kernels

#ifdef _OPENACC
    !$acc parallel loop gang
#else
    !$omp parallel do private(i,j,ij)
#endif
     do i=1,mpsi-1
        !$acc loop vector
        do j=1,mtheta(i)
           ij=igrid(i)+j
           phism(ij)=0.25*((1.0-wtp1(1,ij))*phitmp(jtp1(1,ij))+wtp1(1,ij)*phitmp(jtp1(1,ij)+1)+&
                (1.0-wtp1(2,ij))*phitmp(jtp1(2,ij))+wtp1(2,ij)*phitmp(jtp1(2,ij)+1))-&
                0.0625*((1.0-wtp2(1,ij))*phitmp(jtp2(1,ij))+wtp2(1,ij)*phitmp(jtp2(1,ij)+1)+&
                (1.0-wtp2(2,ij))*phitmp(jtp2(2,ij))+wtp2(2,ij)*phitmp(jtp2(2,ij)+1))
        enddo
     enddo

   !todo
   if(.false.)then
!reflective outer boundary
!$omp parallel do private(j,ij)
     do j=1,mtheta(mpsi-1)
         ij=igrid(mpsi-1)+j
       if(iload/=9)then
         phism(ij)=phism(ij)+0.125*((1.0-wtp2(1,ij))*phitmp(jtp2(1,ij))+wtp2(1,ij)*phitmp(jtp2(1,ij)+1))
       else !Non-zero outer boundary
         phism(ij)=phism(ij)+0.125*((1.0-wtp2(1,ij))*phitmp(jtp2(1,ij))+wtp2(1,ij)*phitmp(jtp2(1,ij)+1))-&
                   0.125*((1.0-wtp1(1,ij))*phitmp(jtp1(1,ij))+wtp1(1,ij)*phitmp(jtp1(1,ij)+1))
       endif
     enddo
!reflective inner boundary
!$omp parallel do private(j,ij)
     do j=1,mtheta(1)
         ij=igrid(1)+j
       if(iload/=9)then
         phism(ij)=phism(ij)+0.125*((1.0-wtp2(2,ij))*phitmp(jtp2(2,ij))+wtp2(2,ij)*phitmp(jtp2(2,ij)+1))
       else !Non-zero inner boundary
         phism(ij)=phism(ij)+0.125*((1.0-wtp2(2,ij))*phitmp(jtp2(2,ij))+wtp2(2,ij)*phitmp(jtp2(2,ij)+1))-&
                   0.125*((1.0-wtp1(1,ij))*phitmp(jtp1(2,ij))+wtp1(2,ij)*phitmp(jtp1(2,ij)+1))
       endif
     enddo
   endif

#ifdef _OPENACC
    !$acc parallel loop gang
#else
    !$omp parallel do private(i,j,ij)
#endif
    do i=1,mpsi-1
       !$acc loop vector
       do j=1,mtheta(i)
          ij=igrid(i)+j
          phitmp(ij)=0.625*phitmp(ij)+phism(ij)
       enddo
    enddo

end subroutine smooth_perp
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine smooth_perp_fem(phitmp)
  use precision, only: lk
  use global_parameters
  use field_array
  use particle_array,only:iload
  implicit none

  integer i,j,k,ii,ij,ij_fem,ip,jt
  real(lk) phism(mgrid),ptemp,phitmp(mgrid),&
       smat(mindex_fem,mgrid_fem),rowsum(mgrid_fem),phitmp_fem(mgrid_fem),phism_fem(mgrid_fem)

  !$acc declare create(phism, phitmp, smat, rowsum, phitmp_fem, phism_fem)

  !$acc kernels
   phism=0.0
   !$acc end kernels

#ifdef _OPENACC
   !$acc parallel loop gang
#else
   !$omp parallel do private(i,j,ij,ij_fem)
#endif
   do i=0,mpsi
      !$acc loop vector
      do j=1,mtheta(i)
         ij=igrid(i)+j
         ij_fem=igrid_fem(i)+j
         phitmp_fem(ij_fem)=phitmp(ij)
      enddo
   enddo
#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_fem(i)+j
            rowsum(ij)=real(nindex_fem(ij)-1) 
         enddo
      enddo 
#ifdef _OPENACC
   !$acc parallel loop gang
#else
   !$omp parallel do private(i,j,ij,k)
#endif

      do i=0,mpsi
         !$acc loop vector
         do j=1,mtheta(i)
            ij=igrid_fem(i)+j
            do k=1,nindex_fem(ij)
              if(ij==indexp_fem(k,ij))then 
                smat(k,ij)=0.6
              else
                smat(k,ij)=0.4/rowsum(ij)
              endif
            enddo
         enddo
      enddo

#ifdef _OPENACC
   !$acc parallel loop gang
#else
   !$omp parallel do private(i,j,ij,k)
#endif
      do i=0,mpsi
         !$acc loop vector
         do j=1,mtheta(i)
            ij=igrid_fem(i)+j
            phism_fem(ij)=0.0_lk
            do k=1,nindex_fem(ij)
               phism_fem(ij)=phism_fem(ij)+smat(k,ij)*phitmp_fem(indexp_fem(k,ij))
            enddo
         enddo
      enddo

#ifdef _OPENACC
   !$acc parallel loop gang
#else
   !$omp parallel do private(i,j,ij,ij_fem)
#endif
      do i=0,mpsi
         !$acc loop vector
         do j=1,mtheta(i)
            ij=igrid(i)+j
            ij_fem=igrid_fem(i)+j
            phism(ij)=phism_fem(ij_fem)
         enddo
      enddo

!reflective outer boundary
#ifdef _OPENACC
   !$acc parallel loop gang vector
#else
   !$omp parallel do private(j,ij)
#endif
     do j=1,mtheta(mpsi-1)
         ij=igrid(mpsi-1)+j
       if(iload/=9)then
         phism(ij)=phism(ij)+0.125*((1.0-wtp2(1,ij))*phitmp(jtp2(1,ij))+wtp2(1,ij)*phitmp(jtp2(1,ij)+1))
       else !Non-zero outer boundary
         phism(ij)=phism(ij)+0.125*((1.0-wtp2(1,ij))*phitmp(jtp2(1,ij))+wtp2(1,ij)*phitmp(jtp2(1,ij)+1))-&
                   0.125*((1.0-wtp1(1,ij))*phitmp(jtp1(1,ij))+wtp1(1,ij)*phitmp(jtp1(1,ij)+1))
       endif
     enddo

#ifdef _OPENACC
   !$acc parallel loop gang
#else
   !$omp parallel do private(i,j,ij)
#endif
     do i=1,mpsi-1
        !$acc loop vector
        do j=1,mtheta(i)
           ij=igrid(i)+j
           phitmp(ij)=0.625*phitmp(ij)+0.375*phism(ij)
        enddo
     enddo

end subroutine smooth_perp_fem

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine smooth_para(phitmp)
  use system_env
  use precision, only: lk, mpi_Rsize
  use global_parameters
  use field_array
#ifdef TIME_MPI
  use mpi_timing
#endif
  implicit none

  integer i,j,ii,ij,ip,jt,&
       icount,idest,isource,isendtag,&
       irecvtag,istatus(MPI_STATUS_SIZE),smoothpoint
  real(lk) sendl(mgrid),sendr(mgrid),recvl(mgrid),recvr(mgrid),phitmp(mgrid),&
       ptemp,phileft(mgrid),phiright(mgrid)
#ifdef _OPENACC
  real(lk) pleft,pright
#else
  real(lk) pleft(mthetamax),pright(mthetamax)
#endif
  !$acc declare create(sendl,sendr,recvl,recvr,phileft,phiright,phitmp)


! parallel smoothing
   smoothpoint=5
   if(smoothpoint==3)then
! send phi to right and receive from left
#ifdef TIME_MPI
       call timer_start()
#endif
     !$acc kernels
     sendr=phitmp
     recvl=0.0
     !$acc end kernels
     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

! send phi to left and receive from right
     !$acc kernels
     sendl=phitmp
     recvr=0.0
     !$acc end kernels
     idest=left_pe
     isource=right_pe
     isendtag=myrank_toroidal
     irecvtag=isource

#ifdef CUDA_MPI
     !$acc host_data use_device(sendl,recvr)
#else
     !$acc update host(sendl)
#endif
     call MPI_SENDRECV(sendl,icount,mpi_Rsize,idest,isendtag,recvr,icount,&
         mpi_Rsize,isource,irecvtag,toroidal_comm,istatus,ierror)
#ifdef CUDA_MPI
     !$acc end host_data
#else
     !$acc update device(recvr)
#endif
#ifdef TIME_MPI
    smooth_sendrecv_time=smooth_sendrecv_time+timer_get_et()
    smooth_sendrecv_count=smooth_sendrecv_count+2
    smooth_sendrecv_size=smooth_sendrecv_size+4.0*icount*lk
#endif
#ifndef _OPENACC
!$omp parallel do private(i,ii,jt,j,ij,ptemp,pleft,pright)
     do i=1,mpsi-1
        ii=igrid(i)
        jt=mtheta(i)
        if(myrank_toroidal==0)then !down-shift for zeta=0
           pleft(1:jt)=cshift(recvl(ii+1:ii+jt),-itran(i))
           pright(1:jt)=recvr(ii+1:ii+jt)
        elseif(myrank_toroidal==mtoroidal-1)then !up-shift for zeta=2*pi
           pright(1:jt)=cshift(recvr(ii+1:ii+jt),itran(i))
           pleft(1:jt)=recvl(ii+1:ii+jt)
        else
           pleft(1:jt)=recvl(ii+1:ii+jt)
           pright(1:jt)=recvr(ii+1:ii+jt)
        endif

        do j=1,mtheta(i)
           ij=igrid(i)+j
           ptemp=phitmp(ij)
           phitmp(ij)=0.5*ptemp+0.25*(pleft(j)+pright(j))
        enddo
    enddo
#else
    !$acc parallel loop gang 
    do i=1,mpsi-1
       ii=igrid(i)
       jt=mtheta(i)
       !$acc loop vector
       do j=1,jt 
          if(myrank_toroidal==0)then !down-shift for zeta=0
             !pleft=cshift(recvl(ii+j),-itran(i))
             pleft=recvl(ii+1+modulo(j-itran(i)-1,jt))
             pright=recvr(ii+j)
          elseif(myrank_toroidal==mtoroidal-1)then !up-shift for zeta=2*pi
             !pright=cshift(recvr(ii+j),itran(i))
             pright=recvr(ii+1+modulo(j+itran(i)-1,jt))
             pleft=recvl(ii+j)
          else
             pleft=recvl(ii+j)
             pright=recvr(ii+j)
          endif

          ij=igrid(i)+j
          phitmp(ij)=0.5*phitmp(ij)+0.25*(pleft+pright)
       enddo
     enddo
#endif
   else ! 5-point parallel smoothing

#ifdef TIME_MPI
       call timer_start()
#endif

     ! send phi to right and receive from left - 1st time
     !$acc kernels
     sendr=phitmp
     recvl=0.0
     !$acc end kernels
     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


#ifndef _OPENACC
!$omp parallel do private(i,j,ii,jt,ij,ptemp,pleft)
     do i=1,mpsi-1
        ii=igrid(i)
        jt=mtheta(i)
        pleft(1:jt)=recvl(ii+1:ii+jt)
        if(myrank_toroidal==0)pleft(1:jt)=cshift(recvl(ii+1:ii+jt),-itran(i)) !down-shift for zeta=0
        do j=1,mtheta(i)
           ij=igrid(i)+j
           ptemp=phitmp(ij)
           phileft(ij)=ptemp-0.25*pleft(j)
        enddo
     enddo
#else
     !$acc parallel loop gang
     do i=1,mpsi-1
        ii=igrid(i)
        jt=mtheta(i)
        !$acc loop vector
        do j=1,jt
           ij=igrid(i)+j
           if(myrank_toroidal==0)then
              pleft=recvl(ii+1+modulo(j-itran(i)-1,jt))
           else
              pleft=recvl(ii+j)
           endif
           phileft(ij)=phitmp(ij)-0.25*pleft
        enddo
     enddo
#endif

     ! send phi to right and receive from left - 2nd time
     !$acc kernels
     sendr=phileft
     recvl=0.0
     !$acc end kernels
     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


     ! send phi to left and receive from right - 1st time
     !$acc kernels
     sendl=phitmp
     recvr=0.0
     !$acc end kernels
     idest=left_pe
     isource=right_pe
     isendtag=myrank_toroidal
     irecvtag=isource

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

#ifndef _OPENACC
!$omp parallel do private(i,j,ii,jt,ij,ptemp,pright)
     do i=1,mpsi-1
        ii=igrid(i)
        jt=mtheta(i)
        pright(1:jt)=recvr(ii+1:ii+jt)
        if(myrank_toroidal==mtoroidal-1)pright(1:jt)=cshift(recvr(ii+1:ii+jt),itran(i)) !up-shift for zeta=2*pi
        do j=1,mtheta(i)
           ij=igrid(i)+j
           ptemp=phitmp(ij)
           phiright(ij)=ptemp-0.25*pright(j)
        enddo
     enddo
#else
     !$acc parallel loop gang
     do i=1,mpsi-1
        ii=igrid(i)
        jt=mtheta(i)
        !$acc loop vector
        do j=1,jt
           if(myrank_toroidal==mtoroidal-1)then
             pright=recvr(ii+1+modulo(j+itran(i)-1,jt))
           else
             pright=recvr(ii+j)
           endif

           ij=igrid(i)+j
           phiright(ij)=phitmp(ij)-0.25*pright
        enddo
     enddo
#endif
     ! send phi to left and receive from right - 2nd time
     !$acc kernels
     sendl=phiright
     recvr=0.0
     !$acc end kernels
     idest=left_pe
     isource=right_pe
     isendtag=myrank_toroidal
     irecvtag=isource
#ifdef CUDA_MPI
     !$acc host_data use_device(sendl,recvr)
#else
     !$acc update host(sendl)
#endif
     call MPI_SENDRECV(sendl,icount,mpi_Rsize,idest,isendtag,recvr,icount,&
         mpi_Rsize,isource,irecvtag,toroidal_comm,istatus,ierror)
#ifdef CUDA_MPI
     !$acc end host_data
#else
     !$acc update device(recvr)
#endif

#ifdef TIME_MPI
    smooth_sendrecv_time=smooth_sendrecv_time+timer_get_et()
    smooth_sendrecv_count=smooth_sendrecv_count+2
    smooth_sendrecv_size=smooth_sendrecv_size+4.0*icount*lk
#endif


#ifndef _OPENACC
!$omp parallel do private(i,ii,jt,j,ij,ptemp,pleft,pright)
     do i=1,mpsi-1
        ii=igrid(i)
        jt=mtheta(i)
        if(myrank_toroidal==0)then !down-shift for zeta=0
           pleft(1:jt)=cshift(recvl(ii+1:ii+jt),-itran(i))
           pright(1:jt)=recvr(ii+1:ii+jt)
        elseif(myrank_toroidal==mtoroidal-1)then !up-shift for zeta=2*pi
           pright(1:jt)=cshift(recvr(ii+1:ii+jt),itran(i))
           pleft(1:jt)=recvl(ii+1:ii+jt)
        else
           pleft(1:jt)=recvl(ii+1:ii+jt)
           pright(1:jt)=recvr(ii+1:ii+jt)
        endif

        do j=1,mtheta(i)
           ij=igrid(i)+j
           ptemp=phitmp(ij)
           phitmp(ij)=0.625*ptemp+0.25*(pleft(j)+pright(j))
        enddo
     enddo
#else
     !$acc parallel loop gang
     do i=1,mpsi-1
        ii=igrid(i)
        jt=mtheta(i)
        !$acc loop vector
        do j=1,jt
           if(myrank_toroidal==0)then !down-shift for zeta=0
              pleft=recvl(ii+1+modulo(j-itran(i)-1,jt))
              pright=recvr(ii+j)
           elseif(myrank_toroidal==mtoroidal-1)then !up-shift for zeta=2*pi
              pright=recvr(ii+1+modulo(j+itran(i)-1,jt))
              pleft=recvl(ii+j)
           else
              pleft=recvl(ii+j)
              pright=recvr(ii+j)
           endif

           ij=igrid(i)+j
           phitmp(ij)=0.625*phitmp(ij)+0.25*(pleft+pright)
        enddo
     enddo
#endif
   endif !3 or 5 points
end subroutine smooth_para

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

! radial interpolation: find poloidal grid index across flux surfaces & along constant theta
subroutine rinterpolation
  use utility, only: check_allocation, mstat
  use global_parameters
  use field_array
  use particle_array,only:iload
  implicit none

  integer i,ip,j,indp,indt,ij,jt,ibin,ibout
  real(lk) tdum,wt

! allocate memory
  allocate(jtp1(2,mgrid),jtp2(2,mgrid),wtp1(2,mgrid),wtp2(2,mgrid), STAT=mstat)
  call check_allocation(mstat, "jtp1 in smooth.F90")

     if(iload==9)then
         ibin=0
         ibout=mpsi
     else
         ibin=1
         ibout=mpsi-1
     endif
!$omp parallel do private(i,ip,j,indp,indt,ij,tdum,jt,wt)
  do i=ibin,ibout
     do ip=1,2
        indp=i+ip
        indt=i-ip

!reflective boundary condition
        if (indp>mpsi) then
            indp=2*mpsi-indp
        endif
        if (indt<0) then
            indt=-indt
        endif

        do j=1,mtheta(i)
           ij=igrid(i)+j
! upward
           if(fielddir==1 .or. fielddir==3)then
             tdum=(real(j)*deltat(i)+(zeta1-pi2)*(qtinv(i)-qtinv(indp)))/deltat(indp)
           else
             tdum=(real(j)*deltat(i)+zeta1*(qtinv(i)-qtinv(indp)))/deltat(indp)
           endif
           jt=floor(tdum)
           wt=tdum-real(jt)
           jt=mod(jt+mtheta(indp),mtheta(indp))
           if(ip==1)then
              wtp1(1,ij)=wt
              jtp1(1,ij)=igrid(indp)+jt
           else
              wtp2(1,ij)=wt
              jtp2(1,ij)=igrid(indp)+jt
           endif
! downward
           if(fielddir==1 .or. fielddir==3)then
             tdum=(real(j)*deltat(i)+(zeta1-pi2)*(qtinv(i)-qtinv(indt)))/deltat(indt)
           else
             tdum=(real(j)*deltat(i)+zeta1*(qtinv(i)-qtinv(indt)))/deltat(indt)
           endif
           jt=floor(tdum)
           wt=tdum-real(jt)
           jt=mod(jt+mtheta(indt),mtheta(indt))
           if(ip==1)then
              wtp1(2,ij)=wt
              jtp1(2,ij)=igrid(indt)+jt
           else
              wtp2(2,ij)=wt
              jtp2(2,ij)=igrid(indt)+jt
           endif
        enddo
     enddo
  enddo
if(iload/=9)then
!!****end points
!! upward
  i=0
!$omp parallel do private(ip,j,indp,ij,tdum,jt,wt)
  do ip=1,2
     indp=min(mpsi,i+ip)
     do j=1,mtheta(i)
           ij=igrid(i)+j
           if(fielddir==1 .or. fielddir==3)then
             tdum=(real(j)*deltat(i)+(zeta1-pi2)*(qtinv(i)-qtinv(indp)))/deltat(indp)
           else
             tdum=(real(j)*deltat(i)+zeta1*(qtinv(i)-qtinv(indp)))/deltat(indp)
           endif
           jt=floor(tdum)
           wt=tdum-real(jt)
           jt=mod(jt+mtheta(indp),mtheta(indp))
           if(ip==1)then
              wtp1(1,ij)=wt
              jtp1(1,ij)=igrid(indp)+jt
           else
              wtp2(1,ij)=wt
              jtp2(1,ij)=igrid(indp)+jt
           endif
     enddo
 enddo
!! downward
  i=mpsi
!$omp parallel do private(ip,j,indt,ij,tdum,jt,wt)
  do ip=1,2
     indt=max(0,i-ip)
     do j=1,mtheta(i)
           ij=igrid(i)+j
           if(fielddir==1 .or. fielddir==3)then
             tdum=(real(j)*deltat(i)+(zeta1-pi2)*(qtinv(i)-qtinv(indt)))/deltat(indt)
           else
             tdum=(real(j)*deltat(i)+zeta1*(qtinv(i)-qtinv(indt)))/deltat(indt)
           endif
           jt=floor(tdum)
           wt=tdum-real(jt)
           jt=mod(jt+mtheta(indt),mtheta(indt))
           if(ip==1)then
              wtp1(2,ij)=wt
              jtp1(2,ij)=igrid(indt)+jt
           else
              wtp2(2,ij)=wt
              jtp2(2,ij)=igrid(indt)+jt
           endif
     enddo
  enddo
endif
!$acc update device(jtp1,jtp2,wtp1,wtp2)
end subroutine rinterpolation
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! d_phieff/d_t for fluid-kinetic hybrid electron model
subroutine dphieff
  use global_parameters
  use particle_array
  use field_array
  implicit none

  integer ij,j,i,jt,iconverged
  real(lk) deltime, dtInv

  deltime=max(1.0e-10,tstep)
  dtInv = 1.0_lk/deltime
  jt=irk+2*(ihybrid-1)

  if (etrap==2) then
    iconverged=1      !0: Use corresponding method to calculate dphi/dt. This method is not numerically
                 !   stable for small time steps when passing electrons are loaded, i.e. etrap = 2.
                 !1: Use converged method. (Lin 2007)
  else
    iconverged=0
  endif

  if(cs_method==1)then
    if(istep==1 .and. mype==0 .and. irk==1 .and. ihybrid==nhybrid)write(gtcout,*)'Using conservative method for kinetic electrons'
  else
    if(istep==1 .and. mype==0 .and. irk==1 .and. ihybrid==nhybrid .and. iconverged==0)write(gtcout,*)'Using corresponding  method for hybrid electrons'
    if(istep==1 .and. mype==0 .and. irk==1 .and. ihybrid==nhybrid .and. iconverged==1)write(gtcout,*)'Using converged method for hybrid electrons'
  endif


  if (magnetic==0) then
     if(iconverged==0)then !corresponding ihybrid method for dphi/dt  
#ifdef _OPENACC
        !$acc parallel loop gang vector
#else 
        !$omp parallel do private(ij)
#endif
        do ij=1,mgrid
           phit(:,ij)=(phi(:,ij)-phisave(:,ij,jt))*dtInv   
           phisave(:,ij,jt)=phi(:,ij) !phisave is stored per ihybrid iteration
        enddo
     elseif(iconverged==1)then !converged nhybrid method for dphi/dt
#ifdef _OPENACC
        !$acc parallel loop gang vector
#else
        !$omp parallel do private(ij)
#endif
        do ij=1,mgrid
           phit(:,ij)=(phi(:,ij)-phisave(:,ij,irk))*dtInv
           if(ihybrid==nhybrid) phisave(:,ij,irk)=phi(:,ij) !phisave is stored at ihybrid = nhybrid, needed for etrap=2
        enddo
     endif
  else

     if(cs_method==1)then
   
       call dphieff_cs
  
     else 

       if(iconverged==0)then !corresponding ihybrid method for dne/dt
#ifdef _OPENACC
        !$acc parallel loop gang
#else
        !$omp parallel do private(i,j,ij)
#endif
         do i=0,mpsi
            !$acc loop vector
            do j=0,mtheta(i)
              ij=igrid(i)+j
              if(ihybrid==1)then
                 dnet(:,ij)=(sdelfluidne(:,ij)-dnesave(:,ij,jt))*dtInv
                 dnesave(:,ij,jt)=sdelfluidne(:,ij)
              else
                 dnet(:,ij)=(sdelfluidne(:,ij)-densitye(:,ij)*meshne(i)-dnesave(:,ij,jt))*dtInv
                 dnesave(:,ij,jt)=sdelfluidne(:,ij)-densitye(:,ij)*meshne(i) !dnesave is stored per ihybrid iteration       
              endif
              if(deltabpara==1)then
                dbparat(:,ij)=(bpara(:,ij)-bparasave(:,ij,jt))/meshte(i)*dtInv
                bparasave(:,ij,jt)=bpara(:,ij)
              endif
            enddo
         enddo
       elseif(iconverged==1)then !converged nhybrid method for dne/dt
#ifdef _OPENACC
        !$acc parallel loop gang
#else
        !$omp parallel do private(i,j,ij)
#endif
         do i = 0, mpsi
            !$acc loop vector
            do j = 0, mtheta(i)
              ij=igrid(i)+j 
              dnet(:,ij)=(sdelfluidne(:,ij)-densitye(:,ij)*meshne(i)-dnesave(:,ij,irk))*dtInv
              if(ihybrid==nhybrid)then
                 dnesave(:,ij,irk)=sdelfluidne(:,ij) - densitye(:,ij)*meshne(i)
              endif

              if(deltabpara==1)then
                dbparat(:,ij)=(bpara(:,ij)-bparasave(:,ij,irk))/meshte(i)*dtInv
                if(ihybrid==nhybrid)then
                  bparasave(:,ij,irk)=bpara(:,ij)
                endif  
              endif

            enddo
         enddo
       endif
     endif
  endif

end subroutine dphieff
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine removn0(farray)
  use system_env
  use precision, only: lk, mpi_Rsize
  use global_parameters
  use field_array,only:deltat,qtinv,zeta1,mtheta,igrid,nshift,wtshift
#ifdef TIME_MPI
  use mpi_timing
#endif
  implicit none
  integer i,j,ij,jt
  real(lk) farray(0:1,mgrid),phitmp(mgrid),wtmp
  !$acc declare create(phitmp)

!! map alpha grid to theta grid
#ifndef _OPENACC
  !$omp parallel do private(i,wtmp,j,ij)
  do i=0,mpsi
     phitmp(igrid(i)+1:igrid(i)+mtheta(i))=cshift(farray(1,igrid(i)+1:igrid(i)+mtheta(i)),-nshift(i))
     phitmp(igrid(i))=phitmp(igrid(i)+mtheta(i))
     wtmp=wtshift(i)
! linear interpolation to theta grid
     do j=1,mtheta(i)
        ij=igrid(i)+j
        farray(1,ij)=phitmp(ij-1)*wtmp+phitmp(ij)*(1.0_lk-wtmp)
     enddo
     farray(1,igrid(i))=farray(1,igrid(i)+mtheta(i))
  enddo
#else
  !$acc parallel loop gang
  do i=0,mpsi
     wtmp=wtshift(i)
     jt=mtheta(i)
     !$acc loop vector
     do j=1,jt
        ij=igrid(i)+j
        phitmp(ij)=farray(1,igrid(i)+1+modulo(j-nshift(i)-1,jt))
        farray(1,ij)=phitmp(ij-1)*wtmp+phitmp(ij)*(1.0_lk-wtmp)
     enddo
     farray(1,igrid(i))=farray(1,igrid(i)+mtheta(i))
  enddo
#endif

!! remove n=0 m\=0 mode
#ifdef TIME_MPI
  call timer_start()
#endif
#ifdef CUDA_MPI
  !$acc host_data use_device(farray,phitmp)
#else
  !$acc update host(farray(1,:))
#endif
  call MPI_ALLREDUCE(farray(1,:),phitmp,mgrid,mpi_Rsize,MPI_SUM,toroidal_comm,ierror)
#ifdef CUDA_MPI
  !$acc end host_data
#else
  !$acc update device(phitmp)
#endif

#ifdef TIME_MPI
    smooth_allreduce_time=smooth_allreduce_time+timer_get_et()
    smooth_allreduce_count=smooth_allreduce_count+1
    smooth_allreduce_size=smooth_allreduce_size+1.0*mgrid*lk
#endif

  !$acc kernels
  farray(1,:)=farray(1,:)-phitmp/real(mtoroidal)
  !$acc end kernels

!! map theta grid back to alpha grid
#ifndef _OPENACC
  !$omp parallel do private(i,wtmp,j,ij)
  do i=0,mpsi
     wtmp=wtshift(i)
     do j=0,mtheta(i)-1
        ij=igrid(i)+j
        phitmp(ij)=farray(1,ij+1)*wtmp+farray(1,ij)*(1.0_lk-wtmp)
     enddo
     phitmp(igrid(i)+mtheta(i))=phitmp(igrid(i))

     farray(1,igrid(i)+1:igrid(i)+mtheta(i))=cshift(phitmp(igrid(i)+1:igrid(i)+mtheta(i)),nshift(i))
     farray(1,igrid(i))=farray(1,igrid(i)+mtheta(i))
  enddo
#else
  !$acc parallel loop gang
  do i=0,mpsi
     wtmp=wtshift(i)
     !$acc loop vector
     do j=0,mtheta(i)-1
        ij=igrid(i)+j
        phitmp(ij)=farray(1,ij+1)*wtmp+farray(1,ij)*(1.0_lk-wtmp)
     enddo
     phitmp(igrid(i)+mtheta(i))=phitmp(igrid(i))
     !$acc loop vector
     do j=1,mtheta(i)
        ij=igrid(i)+j
        farray(1,ij)=phitmp(igrid(i)+1+modulo(j+nshift(i)+1,jt))
     enddo
     farray(1,igrid(i))=farray(1,igrid(i)+mtheta(i))
  enddo
#endif

end subroutine removn0


subroutine axisExtrapolate(farray)
  use precision
  use global_parameters,only: mgrid,rg0,mpsilow
  use field_array,only: mtheta,igrid,deltar,deltat,zeta1,qtinv
  implicit none

  integer :: i,j,ij
  real(lk),dimension(mgrid) :: farray
  real(lk) ave,cosave,sinave,tdum,rdum

  !Extrapolate to the magnetic axis.
  ave=0.0
  cosave=0.0
  sinave=0.0

#ifdef _OPENACC
  !$acc parallel loop gang
#else
  !$omp parallel do private(i,rdum,j,ij,tdum) reduction(+:ave,cosave,sinave)
#endif
  do i=mpsilow+1, mpsilow+2
    rdum=(rg0+deltar*real(mpsilow+1))/(real(mtheta(i))*(rg0+deltar*real(i)))
    !$acc loop vector reduction(+:ave,cosave,sinave)
    do j=1, mtheta(i)
      ij=igrid(i)+j
      tdum=deltat(i)*real(j)+zeta1*qtinv(i)
      ave=ave+farray(ij)*rdum
      cosave=cosave+2.0*farray(ij)*cos(tdum)*rdum
      sinave=sinave+2.0*farray(ij)*sin(tdum)*rdum
    enddo
  enddo

  sinave=0.5_lk*sinave
  cosave=0.5_lk*cosave
  ave=0.5_lk*ave

!Exterpolation mpsilow flux surface to origin
#ifdef _OPENACC
  !$acc parallel loop gang
#else
  !$omp parallel do private(i,rdum,j,ij,tdum)
#endif
  do i=0, mpsilow
    rdum=(rg0+deltar*real(i))/(rg0+deltar*real(mpsilow+1))
    !$acc loop vector
    do j=1, mtheta(i)
      ij=igrid(i)+j
      tdum=deltat(i)*real(j)+(zeta1)*qtinv(i)
      farray(ij)=ave*rdum !!!
      farray(ij)=farray(ij)+cosave*cos(tdum)*rdum
      farray(ij)=farray(ij)+sinave*sin(tdum)*rdum
    enddo
  enddo

!theta boundary condition
#ifdef _OPENACC
  !$acc parallel loop gang vector
#else
  !$omp parallel do private(i)
#endif
  do i=0, mpsilow+1
      farray(igrid(i))=farray(igrid(i)+mtheta(i))
  enddo
end subroutine axisExtrapolate
