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

! Solve fields: phi_eff, phi_ind, fluidue
subroutine field_solver
  use precision
  use system_env
  use global_parameters
  use particle_array
  use field_array
  use equilibrium, only: psiw
  use conservative_scheme, only: phiind
  implicit none

  integer i,j,ij,n,m,win0,win1
  real(lk) gperp2,g2_inv,wa,envelope,tdum,wt,qmin,d2fluidnetmp(0:1,1:mgrid),&
           d4fluidne1(mgrid),d4fluidne2(mgrid),c1,c2
  real(lk) :: buf(mgrid)
  !$acc declare create(d4fluidne1,d4fluidne2,d2fluidnetmp,buf)

  interface
    subroutine laplacian_gpu(scalar,scalar_out,none0bound)
      use precision
      use global_parameters,only:mgrid
      implicit none

      real(lk),dimension(mgrid),intent(in) :: scalar
      real(lk),dimension(mgrid),intent(out) :: scalar_out
      integer, optional, intent (in) :: none0bound
    end subroutine laplacian_gpu
  end interface

  !$acc kernels
  d2fluidnetmp=0.0_lk
  d4fluidne=0.0_lk
  d4fluidne1=0.0_lk
  d4fluidne2=0.0_lk
  !$acc end kernels

! faked k_perp for benchmarking petsc
  gperp2=0.1_lk !k_perp rho0 value
  gperp2=(gperp2/rho0)**2 !k_perp^2 in basic unit
  g2_inv=1.0_lk/gperp2

#ifdef _USING_EXT_SOLVER
  if((istep==1).and.(mype==0).and.(irk==1))write(gtcout,*)'ILAPLACIAN=',ilaplacian
  if(ilaplacian==0) then
  !   call laplacian_petsc(sdelapara(1,:),d2apara) !use Poisson Eq's metric, valid for k_perp rho <<1
      !$acc update host(sapara(1,:))
      call laplacian_petsc(sapara(1,:),d2apara) !d2apara contains non-zonal
      !$acc update device(d2apara)
  elseif(ilaplacian==1 .and. fem==0)then
  !   call laplacian(sdelapara(1,:),d2apara)! finite-difference for Laplacian
      !$acc kernels
        buf=sapara(1,:)
      !$acc end kernels
      call laplacian_gpu(buf,d2apara) ! d2apara contains non-zonal
      !$acc kernels
        buf=phi(1,:)
      !$acc end kernels
      call laplacian_gpu(buf,d2phi)
      call laplacian_gpu(d2phi,d4phi)
  elseif(ilaplacian==1 .and. fem>0)then 
  !   call laplacian_fem(1,sdelapara(1,:),d2apara)!
      !$acc update host(sapara)
      call laplacian_fem(1,sapara(1,:),d2apara)!d2apara contains non-zonal
      call laplacian_fem(1,phi(1,:),d2phi)
      call laplacian_fem(1,d2phi,d4phi)
      !$acc update device(d2apara,d2phi,d4phi)
  endif
#else
  if(mype==0)then
    write(gtcout,*) "Please re-compile GTC with PETSc or other solver to do electromagnetic simulations."
    flush(gtcout)	
  endif
  call MPI_FINALIZE(ierror)
  call exit(1)
#endif

! calculate 4-th order derivative of fluid density in parallel direction for hyperviscosity
! hyperviscosity coefficients (c1 for parallel and c2 for perpendicular)
  c1=-hypr1/(mtoroidal*mtoroidal*mtoroidal*mtoroidal)  !*tstep)
  c2=-hypr2/(mthetamax*mthetamax*mthetamax*mthetamax) !*tstep)
  if (abs(c1)+abs(c2)>1.0e-10) then
   !$acc kernels
   d4fluidne1(:)=sfluidne(1,:)
   d4fluidne2(:)=sfluidne(1,:)
   !$acc end kernels
    do n = 1, 2
      call grad2para(d4fluidne1,d2fluidnetmp)
      !$acc kernels
      d4fluidne1(:)=d2fluidnetmp(1,:)
      !$acc end kernels
      call grad2psi(d4fluidne2,d2fluidnetmp)
      !$acc kernels
      d4fluidne2(:)=d2fluidnetmp(1,:)
      !$acc end kernels
    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=0,mtheta(i)   ! j begins from 0 so that we don't need poloidal BC
        ij=igrid(i)+j
        d4fluidne(ij)=c1*d4fluidne1(ij)+c2*d4fluidne2(ij)*1.0e-12
     enddo
    enddo
  endif

! solver fluid_ue in EM simulation
#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

! Electrostatic potential; , phi=[ B_0(R)/B_0(R_0)]^2 (delta_ne)
! Reduced MHD: k_perp rho_s <<1; GK Poisson Eq. becomes Laplace's Eq.
!        phi(1,ij)=-sfluidne(1,ij)*g2_inv*b*b !use constant k_perp; overide petsc solution

! Electron parallel flow from Ampere's law (Holod09: Eq. 11)
! d2apara is the perpendicular laplacian in GTC unit
        fluidue(1,ij)=2.0_lk*d2apara(ij)*rho0*rho0/(betae*bmesh(ij))+&
             meshni(i)*qion*flowi(1,ij)/bmesh(ij)
        if(izonal>0)fluidue(1,ij)=fluidue(1,ij)+meshni(i)*qion*zonalci(i)/bmesh(ij)
        if(fload>0)then
           fluidue(1,ij)=fluidue(1,ij)+meshnf(i)*qfast*flowf(1,ij)/bmesh(ij)
           if(izonal>0)fluidue(1,ij)=fluidue(1,ij)+meshnf(i)*qfast*zonalcf(i)/bmesh(ij)
        endif

         if(feload>0)then
           fluidue(1,ij)=fluidue(1,ij)+meshnfe(i)*qfaste*flowfe(1,ij)/bmesh(ij)
           if(izonal>0)fluidue(1,ij)=fluidue(1,ij)+meshnfe(i)*qfaste*zonalcfe(i)/bmesh(ij)
        endif

! use constant k_perp in Ampere's law for benchmarking Laplacian operator
!        fluidue(1,ij)=(-gperp2*sapara(1,ij)*rho0*rho0/betae)/bmesh(ij)
     enddo
     fluidue(1,igrid(i))=fluidue(1,igrid(i)+mtheta(i))
  enddo

! smooth potential
  CALL SMOOTH_GPU(fluidue)

! use one or a few toroidal and poloidal modes
  if(nfilter>0)then
     CALL FILTER_GPU(fluidue)
!    if(izonal>0)call filter_keep_zonal_gpu(fluidue,fluidue00) do not keep higher order term: zonal fluidue
  endif

end subroutine field_solver
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine field_gradient
  use global_parameters
  use field_array
  use particle_array
  use magnetic_island
  use equilibrium,only: mesher
  use conservative_scheme,only:sapara_na,gradapara_na,phiind,gradphiind,tapara_na_weight
  use antenna, only: antenna_advance_tstep
  use interfaces
  implicit none

  integer i,j,ij,ibin,ibout,switchepara
  real(lk) rg,q,b_inv,q_inv,ne,ne_inv,te,dnedp,dtedp,deffdp,deffdt,deffdz,ddnedp,ddnedt,ddnedz,dlamdt,dlamdz,&
           dpsidp,dpsidt,dpsidz,ddpdp,ddpdt,ddpdz,phipEq
  real(lk) :: buf1(mgrid),buf2(0:1,mgrid),tmp1_field(0:1,mgrid),tmp2_field(0:1,mgrid),tmp3_field(0:1,mgrid)
  !$acc declare create(buf1,buf2,tmp1_field, tmp2_field, tmp3_field)

  switchepara=1   ! 0: set \delta E_||=0 for singular fluid simulation
  ! prepare the external antenna array used in push-field or poisson
  if(antenna_on>0) then
    call antenna_advance_tstep
  endif

! prepare the field gradients used in push particle and field

  ! field gradients
  if(magnetic==1)then
     !$acc kernels
     buf1=fluidue(1,:)
     !$acc end kernels
     call gradient_gpu(buf1,gradue)

     !$acc kernels
     buf1=sapara(1,:)
     !$acc end kernels
     call gradient_gpu(buf1,gradapara)

     !$acc kernels
     buf1=sfluidne(1,:)
     !$acc end kernels
     call gradient_gpu(buf1,gradne)

     !$acc kernels
     buf1=sdelfluidne(1,:)
     !$acc end kernels
     call gradient_gpu(buf1,gradne_nz)

     !$acc kernels
     buf1=sdeltapsi(1,:)
     !$acc end kernels
     call gradient_gpu(buf1,gradpsi)

     if (nhybrid>0) then
       !$acc kernels
       buf1=densitye(1,:)
       !$acc end kernels
       call gradient_gpu(buf1,gradkne)
       !$acc kernels
       buf1=pressureepara(1,:)
       !$acc end kernels
       call gradient_gpu(buf1,gradpepara)
       !$acc kernels
       buf1=pressureeperp(1,:)
       !$acc end kernels
       call gradient_gpu(buf1,gradpeperp)
     endif
     if (deltabpara==1)then
        !$acc kernels
        buf1=bpara(1,:)
        !$acc end kernels
        call gradient_gpu(buf1,gradbpara)
        if(nfilter>0)then
           !$acc kernels
           buf2=gradbpara(1,:,:)
           !$acc end kernels
           CALL FILTER_GPU(buf2)
           !$acc kernels
           gradbpara(1,:,:)=buf2
           buf2=gradbpara(2,:,:)
           !$acc end kernels
           CALL FILTER_GPU(buf2)
           !$acc kernels
           gradbpara(2,:,:)=buf2
           buf2=gradbpara(3,:,:)
           !$acc end kernels
           CALL FILTER_GPU(buf2)
           !$acc kernels
           gradbpara(3,:,:)=buf2
           !$acc end kernels
        else
           !$acc kernels
           buf2=gradbpara(1,:,:)
           !$acc end kernels
           call periodicity_gpu(buf2)
           !$acc kernels
           gradbpara(1,:,:)=buf2
           buf2=gradbpara(2,:,:)
           !$acc end kernels
           call periodicity_gpu(buf2)
           !$acc kernels
           gradbpara(2,:,:)=buf2
           buf2=gradbpara(3,:,:)
           !$acc end kernels
           call periodicity_gpu(buf2)
           !$acc kernels
           gradbpara(3,:,:)=buf2
           !$acc end kernels
        endif
     endif

     if(cs_method==1)then
       !$acc kernels
       buf1=sapara_na(1,:)
       !$acc end kernels
       call gradient_gpu(buf1,gradapara_na)
     elseif(cs_method==0)then
       !$acc kernels
       buf1=saparanl(1,:)
       !$acc end kernels
       call gradient_gpu(buf1,gradapara_na)
     endif

     if(switchepara==0)then
!       ideal MHD 
        gradphieff(:,1,ij)=0.0_lk
      else
!       gyrokinetic/tearing
! gradient of phi_eff
#ifdef _OPENACC
        !$acc parallel loop gang 
#else 
        !$omp parallel do private(i,j,ij,b_inv,q_inv,ne,te,dnedp,dtedp,deffdp,deffdt,deffdz,&
        !$omp& ddnedp,ddnedt,ddnedz,dlamdt,dlamdz,dpsidp,dpsidt,dpsidz,ddpdp,ddpdt,ddpdz)
#endif
        do i=0,mpsi
           q_inv=1.0_lk/qmesh(i)
           if(qmesh(i)<1e-6_lk)q_inv=0.0_lk
           !$acc loop vector
           do j=1,mtheta(i)
              ij=igrid(i)+j
              b_inv=1.0_lk/bmesh(ij)

              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

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

              dpsidp=gradpsi(1,1,ij)
              dpsidt=gradpsi(2,1,ij)
              dpsidz=gradpsi(3,1,ij)-gradpsi(2,1,ij)*q_inv

              deffdp=-(ddnedp-dpsidp*dnedp)*qelectron*te/ne+dpsidp*phipEq
              deffdt=-(ddnedt-dpsidt*dnedp)*qelectron*te/ne+dpsidt*phipEq
              deffdz=-(ddnedz-dpsidz*dnedp)*qelectron*te/ne+dpsidz*phipEq

              gradphieff(1,1,ij)=deffdp + gradkne(1,1,ij)*qelectron*te
              gradphieff(2,1,ij)=deffdt + gradkne(2,1,ij)*qelectron*te
              gradphieff(3,1,ij)=deffdz + deffdt*q_inv + (gradkne(3,1,ij)*qelectron*te+&
                                 eta*2.0_lk*d2apara(ij)*rho0*rho0/betae)
              if(deltabpara==1)then  
                gradphieff(1,1,ij)=gradphieff(1,1,ij)+gradbpara(1,1,ij)*te/rho0**2 
                gradphieff(2,1,ij)=gradphieff(2,1,ij)+gradbpara(2,1,ij)*te/rho0**2
                gradphieff(3,1,ij)=gradphieff(3,1,ij)+gradbpara(3,1,ij)*te/rho0**2
              endif
           enddo
        enddo
        !$acc kernels
        tmp1_field=gradphieff(1,:,:)
        tmp2_field=gradphieff(2,:,:)
        tmp3_field=gradphieff(3,:,:)
        !$acc end kernels
        if(nfilter>0)then
           CALL FILTER_GPU(tmp1_field)
           CALL FILTER_GPU(tmp2_field)
           CALL FILTER_GPU(tmp3_field)
        else
           call periodicity_gpu(tmp1_field)
           call periodicity_gpu(tmp2_field)
           call periodicity_gpu(tmp3_field)
        endif
        !$acc kernels
        gradphieff(1,:,:) = tmp1_field
        gradphieff(2,:,:) = tmp2_field
        gradphieff(3,:,:) = tmp3_field
        !$acc end kernels
     endif
  endif

! grad phi_electrostatic
  !$acc kernels
  buf1=phi(1,:)
  !$acc end kernels
  call gradient_gpu(buf1,gradphi)

  if(cs_method == 1)then
#ifdef _OPENACC
     !$acc parallel loop gang
#else
     !$omp parallel do private(i,j,ij,ne,ne_inv,te,dnedp)
#endif
     do i=0,mpsi
        ne=meshne(i)
        ne_inv=1.0_lk/ne
        te=meshte(i)
        dnedp=-ne*kapane(i)
        !$acc loop vector
        do j=1,mtheta(i)
           ij=igrid(i)+j
           gradphiind(1,1,ij)=te*ne_inv*(gradne_nz(1,1,ij)-dnedp*gradpsi(1,1,ij))-gradphi(1,1,ij) !Note that gradphi(1,1,:) is purely non-zonal here
           gradphiind(2,1,ij)=te*ne_inv*(gradne_nz(2,1,ij)-dnedp*gradpsi(2,1,ij))-gradphi(2,1,ij)
           gradphiind(3,1,ij)=te*ne_inv*(gradne_nz(3,1,ij)-dnedp*gradpsi(3,1,ij))-gradphi(3,1,ij)
           gradphieff(1,1,ij)=te*ne_inv*(gradne_nz(1,1,ij)-dnedp*gradpsi(1,1,ij)) !Only used for push hybrid particle in cs_model
           gradphieff(2,1,ij)=te*ne_inv*(gradne_nz(2,1,ij)-dnedp*gradpsi(2,1,ij))
           gradphieff(3,1,ij)=te*ne_inv*(gradne_nz(3,1,ij)-dnedp*gradpsi(3,1,ij))
           if(deltabpara==1)then  
                gradphieff(1,1,ij)=gradphieff(1,1,ij)+gradbpara(1,1,ij) 
                gradphieff(2,1,ij)=gradphieff(2,1,ij)+gradbpara(2,1,ij)
                gradphieff(3,1,ij)=gradphieff(3,1,ij)+gradbpara(3,1,ij) 
                gradphiind(1,1,ij)=gradphiind(1,1,ij)+gradbpara(1,1,ij) 
                gradphiind(2,1,ij)=gradphiind(2,1,ij)+gradbpara(2,1,ij)
                gradphiind(3,1,ij)=gradphiind(3,1,ij)+gradbpara(3,1,ij) 
           endif
        enddo
       gradphiind(1,1,igrid(i))=gradphiind(1,1,igrid(i)+mtheta(i))
       gradphiind(2,1,igrid(i))=gradphiind(2,1,igrid(i)+mtheta(i))
       gradphiind(3,1,igrid(i))=gradphiind(3,1,igrid(i)+mtheta(i))
       gradphieff(1,1,igrid(i))=gradphieff(1,1,igrid(i)+mtheta(i))
       gradphieff(2,1,igrid(i))=gradphieff(2,1,igrid(i)+mtheta(i))
       gradphieff(3,1,igrid(i))=gradphieff(3,1,igrid(i)+mtheta(i))
     enddo
     if(nfilter>0)then
        !$acc kernels
        tmp1_field=gradphiind(1,:,:)
        tmp2_field=gradphiind(2,:,:)
        tmp3_field=gradphiind(3,:,:)
        !$acc end kernels
        CALL FILTER_GPU(tmp1_field)
        CALL FILTER_GPU(tmp2_field)
        CALL FILTER_GPU(tmp3_field)
        !$acc kernels
        gradphiind(1,:,:) = tmp1_field
        gradphiind(2,:,:) = tmp2_field
        gradphiind(3,:,:) = tmp3_field
        !$acc end kernels

        !$acc kernels
        tmp1_field=gradphieff(1,:,:)
        tmp2_field=gradphieff(2,:,:)
        tmp3_field=gradphieff(3,:,:)
        !$acc end kernels
        CALL FILTER_GPU(tmp1_field)
        CALL FILTER_GPU(tmp2_field)
        CALL FILTER_GPU(tmp3_field)
        !$acc kernels
        gradphieff(1,:,:) = tmp1_field
        gradphieff(2,:,:) = tmp2_field
        gradphieff(3,:,:) = tmp3_field
        !$acc end kernels

     else
        !$acc kernels
        tmp1_field=gradphiind(1,:,:)
        tmp2_field=gradphiind(2,:,:)
        tmp3_field=gradphiind(3,:,:)
        !$acc end kernels
        call periodicity_gpu(tmp1_field)
        call periodicity_gpu(tmp2_field)
        call periodicity_gpu(tmp3_field)
        !$acc kernels
        gradphiind(1,:,:) = tmp1_field
        gradphiind(2,:,:) = tmp2_field
        gradphiind(3,:,:) = tmp3_field
        !$acc end kernels

        !$acc kernels
        tmp1_field=gradphieff(1,:,:)
        tmp2_field=gradphieff(2,:,:)
        tmp3_field=gradphieff(3,:,:)
        !$acc end kernels
        call periodicity_gpu(tmp1_field)
        call periodicity_gpu(tmp2_field)
        call periodicity_gpu(tmp3_field)
        !$acc kernels
        gradphieff(1,:,:) = tmp1_field
        gradphieff(2,:,:) = tmp2_field
        gradphieff(3,:,:) = tmp3_field
        !$acc end kernels

     endif
  endif

  ! grad alpha_island
  if(island==1)then
    !$acc kernels
    buf1=alphaIs(1,:)
    !$acc end kernels
    call gradient_gpu(buf1,gradalphaIs)
  endif

! add (0,0) mode to d phi/d psi
  if(izonal==0)then
     !$acc kernels
     phip00=0.0_lk
     phi00=0.0_lk
     !$acc end kernels
  elseif(izonal==1)then
! solve zonal flow use new method 
     if(total_pssolver==1)then
       call gradient_psi(phi00,phip00)
     else
       call zonal  !zonal flow solver for general geometry
     endif
     if(iload==9)then
         ibin=0
         ibout=mpsi
     else
         ibin=1
         ibout=mpsi-1
     endif
#ifdef _OPENACC
     !$acc parallel loop gang
#else
     !$omp parallel do private(i)
#endif
     do i=ibin,ibout
        !$acc loop vector
        do j=0,mtheta(i)
           ij=igrid(i)+j
           gradphi(1,0:1,ij)=gradphi(1,0:1,ij)+phip00(i)
        enddo
     enddo
  else
  ! solve zonal flow
     call zonal_old    !old verson solver
#ifdef _OPENACC
     !$acc parallel loop gang vector
#else
     !$omp parallel do private(i,rg,q)
#endif
     do i=1,mpsi-1
        rg=rg0+deltar*real(i)
        q=qmesh(i)
        gradphi(1,0:1,igrid(i):igrid(i)+mtheta(i))=gradphi(1,0:1,igrid(i):igrid(i)+mtheta(i))+&
           phip00(i)*q/rg
     enddo
  endif
  if(magnetic==1)then
     if(nonlinear == 1)then
        call nonlinear_ohm_hybrid    ! Calculate taparanl for nonlinear simulation of hybrid model and diagnositic in cs_model
     endif
     if(cs_method==1)then
        call nonadiabatic_ohm_cs    ! Calculate tapara_na for cs_method, tapara_na is not zero for linear simulation with kinetic electrons and other non-ideal MHD effects
     endif
! Calculate the non-adiabatic fields in particle's weight equation: dBdotdE,taparanl_weight,tapara_na_weight
     call highorder_field_weight 
  endif
end subroutine field_gradient

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine nonlinear_ohm_hybrid
  use global_parameters
  use field_array
  use particle_array
  use equilibrium,only: lsp,spdpsi_inv,spdpsi,gpsi,cpsi,qpsi
  implicit none

  integer i,j,ij,isp
  real(lk) pdum,q_inv,dpx,dp2,g,ri,b_inv,jbb_inv,&
           daparadp,daparadt,daparadz,dphidp,dphidt,dphidz,dbparadp,dbparadt,dbparadz,&
           dphieffdp,dphieffdt,dphieffdz,wphieffnl,wphinl,wbparanl,&
           tapara00nl(0:mpsi),perturb
  !$acc declare create(tapara00nl)

  perturb=real(nonlinear)

#ifdef _OPENACC
   !$acc parallel loop gang
#else
   !$omp parallel do private(i,j,ij,pdum,q_inv,isp,dpx,dp2,&
   !$omp& g,ri,b_inv,jbb_inv,&
   !$omp& daparadp,daparadt,daparadz,&
   !$omp& dphidp,dphidt,dphidz,dbparadp,dbparadt,dbparadz,&
   !$omp& dphieffdp,dphieffdt,dphieffdz,&
   !$omp& wphieffnl,wphinl,wbparanl)
#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
       !$acc loop vector
       do j=1,mtheta(i)
          ij=igrid(i)+j
          b_inv=1.0_lk/bmesh(ij)
          jbb_inv=jmesh(ij)*b_inv*b_inv

          daparadp=gradapara(1,1,ij)
          daparadt=gradapara(2,1,ij)
          daparadz=gradapara(3,1,ij)-gradapara(2,1,ij)*q_inv

          dphidp=gradphi(1,1,ij)-phip00(i)
          dphidt=gradphi(2,1,ij)
          dphidz=gradphi(3,1,ij)-gradphi(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(deltabpara==1) then
          dbparadp=gradbpara(1,1,ij)
          dbparadt=gradbpara(2,1,ij)
          dbparadz=gradbpara(3,1,ij)-gradbpara(2,1,ij)*q_inv
          else
          dbparadp=0.0_lk
          dbparadt=0.0_lk
          dbparadz=0.0_lk
          endif
          wphieffnl=(daparadp*(ri*dphieffdz - g*dphieffdt) - dphieffdp*(ri*daparadz - g*daparadt))*jbb_inv
          wphinl=-(daparadp*(ri*dphidz - g*dphidt) - dphidp*(ri*daparadz - g*daparadt))*jbb_inv
          wbparanl=(daparadp*(ri*dbparadz - g*dbparadt) - dbparadp*(ri*daparadz - g*daparadt))*jbb_inv

          taparanl(1,ij)=perturb*(wphieffnl+wphinl+wbparanl)
       enddo
       taparanl(1,igrid(i))=taparanl(1,igrid(i)+mtheta(i))
    enddo


    call smooth_gpu(taparanl)
    if(nfilter>0)then
      if(izonal==0)then
         call filter_gpu(taparanl)
      elseif(izonal>0)then
         call filter_keep_zonal_gpu(taparanl,tapara00nl)
      endif
    endif

end subroutine nonlinear_ohm_hybrid
!=======================================================================
subroutine highorder_field_weight
  use global_parameters
  use field_array
  use particle_array
  use equilibrium, only: lsp,spdpsi_inv,spdpsi,gpsi,cpsi,qpsi
  use conservative_scheme, only: tapara_na,tapara_na_weight
  implicit none

  integer i,j,ij,isp
  real(lk) pdum,q_inv,dpx,dp2,g,ri,b_inv,jbb_inv,&
           daparadp,daparadt,daparadz,&
           dphidp,dphidt,dphidz,&
           wphinl,dBdotdE00(0:mpsi),&
           perturb,weight_perturb

!$acc declare create(dBdotdE00)
  perturb=real(nonlinear)

!======================================================================
!        weight_perturb is a control parameter for removing high order 
!        nonlinearity in the weight equation in nonlinear simulation
!        Suggested by Prof. Lin
!======================================================================

! Local control parameter for the nonlinear parallel electric fields used in weight equation for hybrid model and conservative scheme 
  if(nonlinear==0)then
     weight_perturb=1.0_lk   ! keep tapara_na when cs_method=1 to include kinetic electron responses, nonlinear part of tapara_na is removed in nonadiabatic_ohm_na
     !$acc kernels
     taparanl=0.0_lk         ! safety guard for setting taparanl=0.0_lk in hybrid model linear simulation
     !$acc end kernels
  elseif(nonlinear==1)then
     weight_perturb=0.0_lk   ! remove nonlinear parallel electric fields in weight equation
  endif

#ifdef _OPENACC
  !$acc parallel loop gang
#else   
  !$omp parallel do private(i,j,ij,pdum,q_inv,isp,dpx,dp2,&
  !$omp& g,ri,b_inv,jbb_inv,&
  !$omp& daparadp,daparadt,daparadz,&
  !$omp& dphidp,dphidt,dphidz,wphinl)
#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

       !$acc loop vector
       do j=1,mtheta(i)
          ij=igrid(i)+j
          b_inv=1.0_lk/bmesh(ij)
          jbb_inv=jmesh(ij)*b_inv*b_inv

        daparadp=gradapara(1,1,ij)
        daparadt=gradapara(2,1,ij)
        daparadz=gradapara(3,1,ij)-gradapara(2,1,ij)*q_inv

        dphidp=gradphi(1,1,ij)-phip00(i) ! Remove the zonal component in dBdotdE
        dphidt=gradphi(2,1,ij)
        dphidz=gradphi(3,1,ij)-gradphi(2,1,ij)*q_inv

        wphinl=-(daparadp*(ri*dphidz - g*dphidt) - dphidp*(ri*daparadz - g*daparadt))*jbb_inv
        ! dBdotdE is a nonlinear term, by default it is not applied for weight equation in nonlinear simulation
        dBdotdE(1,ij)=perturb*wphinl
     enddo
     dBdotdE(1,igrid(i))=dBdotdE(1,igrid(i)+mtheta(i))
  enddo

!================================================================= 
    call smooth_gpu(dBdotdE)
    if(nfilter>0)then
      if(izonal==0)then
         call filter_gpu(dBdotdE)
      elseif(izonal>0)then
         call filter_keep_zonal_gpu(dBdotdE,dBdotdE00)
      endif
    endif

! By default, dBdotdE=0 in weight equation for nonlinear simulation (weight_perturb=0)
   !$acc kernels
   dBdotdE=weight_perturb*dBdotdE
   !$acc end kernels

! Hybrid model: taparanl is purely nonlinear when nonlinear=1, and taparanl=0.0_lk when nonlinear=0. By default taparanl is not applied for weight equation in nonlinear simulation by ordering estimation
  !$acc kernels
  taparanl_weight=weight_perturb*taparanl
  !$acc end kernels

! CS_method: tapara_na is purely linear when nonlinear=0(fullohm_na=1), and purely nonlinear(fullohm_na=0) when nonlinear=1. By default tapara_na is not applied for weight equation in nonlinear simulation by ordering estimation. 
! NOT TRUE for nonlinear simulation with kinetic electrons and other non-ideal MHD effects when weight_perturb=0.0
   !$acc kernels
   tapara_na_weight=weight_perturb*tapara_na
   !$acc end kernels

end subroutine highorder_field_weight

!=======================================================================
! old version of zonal flow solver
subroutine zonal_old
  use global_parameters
  use particle_array
  use field_array
  implicit none

  integer i
  real(lk) den00(0:mpsi),rg
  !$acc update device(den00)

! zonal density
  !$acc kernels
  phip00=qion*zonali*meshni+qelectron*fluidne00*real(magnetic)

  if(nhybrid>0)phip00=phip00+qelectron*zonale*meshne
  if(fload>0)phip00=phip00+qfast*zonalf*meshnf
  if(feload>0)phip00=phip00+qfaste*zonalfe*meshnfe
  !$acc end kernels

! (-0.0625 0.25 0.625 0.25 -0.0625) radial smoothing of (0,0) mode density phip00
  do i=1,2
     !$acc kernels
     den00(0)=phip00(0)
     den00(mpsi)=phip00(mpsi)
     den00(1)=phip00(3)
     den00(mpsi-1)=phip00(mpsi-3)
     den00(2:mpsi-2)=phip00(0:mpsi-4)+phip00(4:mpsi)
     den00(1:mpsi-1)=0.625_lk*phip00(1:mpsi-1)+0.25_lk*(phip00(0:mpsi-2)+phip00(2:mpsi))-&
          0.0625_lk*den00(1:mpsi-1)
     phip00=den00
     !$acc end kernels
  enddo

! phi00=r*E_r, E_r(a0)=0. Trapezoid rule
  !$acc kernels
  den00=phip00
  phip00=0.0_lk
  !$acc end kernels

  !$acc parallel loop gang vector
  do i=1,mpsi
     rg=rg0+deltar*real(i)
     phip00(i)=phip00(i-1)+0.5_lk*deltar*((rg-deltar)*den00(i-1)+rg*den00(i))
  enddo

! subtract net momentum
!     phip00=phip00-sum(phip00)/real(mpsi+1)

! d phi/dr, in equilibrium unit

  !$acc parallel loop gang vector
  do i=0,mpsi
     rg=rg0+deltar*real(i)
     phip00(i)=-phip00(i)/rg
  enddo

! add FLR contribution using Pade approximation: b*<phi>=(1+b)*<n>

  !$acc kernels
  phi00=den00*meshti/(qion*qion) !!rho0*rho0 ! potential in equilibrium unit
  !$acc end kernels

  !$acc parallel loop gang vector
  do i=1,mpsi-1
     phip00(i)=(phip00(i)/aion+0.5_lk*(phi00(i+1)-phi00(i-1))/deltar)/meshni(i)
  enddo

! (0,0) mode potential store in phi00
  !$acc kernels
  phi00=0.0_lk
  !$acc end kernels

  !$acc parallel loop gang vector
  do i=1,mpsi
     phi00(i)=phi00(i-1)+0.5_lk*deltar*(phip00(i-1)+phip00(i))
  enddo

  !$acc kernels
  if(izonal==0)phip00=0.0_lk
  !$acc end kernels

end subroutine zonal_old

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

!XY new zonal solver

subroutine zonal
  use global_parameters
  use particle_array
  use field_array
  implicit none

  integer i
  real(lk) den00(0:mpsi),h1,h2,delpsi,temp,psidum(0:mpsi)
  !$acc declare create(den00,psidum)


! zonal density
  !$acc kernels
  phip00=qion*zonali*meshni+qelectron*fluidne00*real(magnetic)
  if(nhybrid>0 .and. cs_method==0)phip00=phip00+qelectron*zonale*meshne
  if(fload>0)phip00=phip00+qfast*zonalf*meshnf
  if(feload>0)phip00=phip00+qfaste*zonalfe*meshnfe
  !$acc end kernels

! (-0.0625 0.25 0.625 0.25 -0.0625) radial smoothing of (0,0) mode density phip00
  do i=1,2
     !$acc kernels
     den00(0)=phip00(0)
     den00(mpsi)=phip00(mpsi)
     den00(1)=phip00(3)
     den00(mpsi-1)=phip00(mpsi-3)
     den00(2:mpsi-2)=phip00(0:mpsi-4)+phip00(4:mpsi)
     den00(1:mpsi-1)=0.625_lk*phip00(1:mpsi-1)+0.25_lk*(phip00(0:mpsi-2)+phip00(2:mpsi))-&
          0.0625_lk*den00(1:mpsi-1)
     phip00=den00
     !$acc end kernels
  enddo
  if(iload/=9)then
!   den00=<e dn>, phip00= d_<e dn>/d_psip, nonuniform psip grids
    !$acc kernels
    den00=phip00
    !$acc end kernels
! initial point
    h1=psimesh(1)-psimesh(0)
    h2=psimesh(2)-psimesh(1)
    phip00(0)=(den00(1)-den00(0))*(h1+h2)/(h1*h2)+(den00(0)-den00(2))*h1/(h2*(h1+h2))
    !$acc parallel loop gang vector
    do i=1,mpsi-1
       h1=psimesh(i)-psimesh(i-1)
       h2=psimesh(i+1)-psimesh(i)
       phip00(i)=(den00(i+1)-den00(i))*h1/(h2*(h1+h2))+(den00(i)-den00(i-1))*h2/(h1*(h1+h2))
    enddo
! final point
    h1=psimesh(mpsi-1)-psimesh(mpsi-2)
    h2=psimesh(mpsi)-psimesh(mpsi-1)
    phip00(mpsi)=(den00(mpsi)-den00(mpsi-1))*(h1+h2)/(h1*h2)+(den00(mpsi-2)-den00(mpsi))*h1/(h2*(h1+h2))

! for integration
    !$acc kernels
    phi00=den00*b2m00/(aion*meshni)+gpsi200*dtndpsi*phip00/(qion*qion)
    phi00=phi00*pmarki

! first term in d_<phi>/d_psip
    den00=phip00*meshti/(qion*qion*meshni)

! phi00=r*E_r, E_r(a0)=0. Trapezoid rule
    phip00=0.0_lk
    !$acc end kernels

    !$acc parallel loop gang vector
    do i=1,mpsi
       phip00(i)=phip00(i-1)+(phi00(i)+phi00(i-1))*0.5_lk
    enddo

! d phi/dr, in equilibrium unit
    !$acc parallel loop gang vector
    do i=0,mpsi
       delpsi=0.5_lk*(psimesh(min(mpsi,i+1))-psimesh(max(0,i-1)))
       phip00(i)=den00(i)-phip00(i)*delpsi/(gpsi200(i)*pmarki(i))
    enddo

  else ! Fully kinetic ion
    !$acc kernels
    den00=phip00*jacobianpsi
    phip00=0.0_lk
    !$acc end kernels
    !$acc parallel loop gang vector
    do i=0,mpsi
! normalized by (lambdaD_i)^2/(R_0)^2, lambdaDi, is the Debye-length, and R_0, is the major radius
       temp=(etemp0)/(eden0)*7.43e2*7.43e2/(r0*r0)
       temp=1.0_lk
       den00(i)=-den00(i)/temp

       if(i>0)then
         phip00(i)=phip00(i-1)+(den00(i)+den00(i-1))*(psimesh(i)-psimesh(i-1))*0.5_lk
       endif
       psidum(i)=0.5_lk*(psimesh(min(mpsi,i+1))-psimesh(max(0,i-1)))
    enddo
    !$acc parallel loop gang vector
    do i=0,mpsi
       phip00(i)=phip00(i)/(gpsi200(i)*jacobianpsi(i))
    enddo

! zonal potential are zero at both boundaries
! should be checked (animesh?)
    !$acc kernels
    phip00=phip00-sum(phip00*psidum)/sum(psidum/(gpsi200*jacobianpsi))/(gpsi200*jacobianpsi)
    !$acc end kernels
  endif

! (0,0) mode potential store in phi00
  !$acc kernels
  phi00=0.0_lk
  !$acc end kernels
  !$acc parallel loop gang vector
  do i=1,mpsi
     phi00(i)=phi00(i-1)+(psimesh(i)-psimesh(i-1))*(phip00(i)+phip00(i-1))*0.5_lk
  enddo
  if(iload==9)then
    !$acc parallel loop gang vector
    do i=0,mpsi
       phip00(i)=phip00(i)*rho0*rho0
       phi00(i)=phi00(i)*rho0*rho0
    enddo
  endif

end subroutine zonal


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

subroutine gradient_gpu(scalar,vector)
  use precision
  use global_parameters
  implicit none
  real(lk),dimension(mgrid),intent(in) :: scalar
  real(lk),dimension(3,0:1,mgrid),intent(out) :: vector

  call gradient_t(scalar,vector)
end subroutine gradient_gpu


subroutine gradient_t(scalar,vector)
  use system_env
  use precision, only: lk, mpi_Rsize
  use global_parameters
  use field_array
  use particle_array,only:iload
#ifdef TIME_MPI
  use mpi_timing
#endif
  implicit none

  integer i,ii,ij,j,k,icount,idest,isource,isendtag,irecvtag,jt,&
      istatus(MPI_STATUS_SIZE)
  real(lk) difft(0:mpsi),diffz,&
      sendl(mgrid),recvr(mgrid),sendr(mgrid),recvl(mgrid),&
      q,delq
#ifdef _OPENACC
  real(lk) pleft,pright
#else
  real(lk) pleft(mthetamax),pright(mthetamax)
#endif
  real(lk),dimension(mgrid),intent(in) :: scalar
  real(lk),dimension(3,0:1,mgrid),intent(out) :: vector
  real(lk) h0,h1,wp0,wp1,wt10,wt00,wt11,wt01,sdum0,sdum1,pdum0,pdum1
  integer jm,jp,jt11,jt21,jt31,jt01,jt10,jt20,jt30,jt00,iout,iin,ibin,ibout
  real(lk),dimension(0:1,mgrid) :: tmp_vector1, tmp_vector2, tmp_vector3
  !$acc declare create(tmp_vector1, tmp_vector2,tmp_vector3,sendr,recvl,sendl,recvr, difft)

  ! finite difference for e-field in equilibrium unit
  difft=0.5_lk/deltat
  !$acc update device(difft)
  diffz=0.5_lk/deltaz

#ifdef _OPENACC
  !$acc parallel loop gang vector
#else
  !$omp parallel do private(i,j,k)
#endif
  do i=1,mgrid
     do j=0,1
        do k=1,3
          vector(k,j,i)=0.0_lk
        enddo
     enddo
   enddo
   
!! d_scalar/d_psito higher order
       if(iload==9)then
         ibin=0
         ibout=mpsi
       else
         ibin=1
         ibout=mpsi-1
       endif

#ifdef _OPENACC
   !$acc parallel loop gang 
#else
   !$omp parallel do private(i,j,ij,h0,h1,wp1,wp0,jm,jp,wt10,wt00,wt11,wt01,&
   !$omp& jt11,jt21,jt31,jt01,jt10,jt20,jt30,jt00,sdum1,sdum0,pdum1,pdum0,iin,iout)
#endif
  do i=ibin,ibout

     iout=i+1
     iin=i-1

     if(i==0)iin=1
     if(i==mpsi)iout=mpsi-1

     h0=1.0_lk/(psimesh(i)-psimesh(iin))
     h1=1.0_lk/(psimesh(iout)-psimesh(i))

     if(i==0)h0=-h0
     if(i==mpsi)h1=-h1

     wp1=h1/(h0+h1)
     wp0=1.0_lk-wp1

     jm=igrid(iin)
     jp=igrid(iout)

     !$acc loop vector
     do j=1,mtheta(i)
        ij=igrid(i)+j
        wt10=wtp1(2,ij)   !upper poloidal grid on inner flux surface
        wt00=1.0_lk-wt10       !lower poloidal grid on inner flux surface
        wt11=wtp1(1,ij)   !upper poloidal grid on outer flux surface
        wt01=1.0_lk-wt11       !lower poloidal grid on outer flux surface

        jt11=jtp1(1,ij)   !!lower poloidal grid on outer flux surface
        jt21=jt11+1
        jt31=jt21+1-mtheta(iout)*((jt21-jp)/mtheta(iout))
        jt01=jp+mod(jt11-jp-1+mtheta(iout),mtheta(iout))

        jt10=jtp1(2,ij)   !!lower poloidal grid on the innoer flux surface
        jt20=jt10+1
        jt30=jt20+1-mtheta(iin)*((jt20-jm)/mtheta(iin))
        jt00=jm+mod(jt10-jm-1+mtheta(iin),mtheta(iin))

        sdum1=wt01*scalar(jt11)+wt11*scalar(jt21)
        sdum0=wt00*scalar(jt10)+wt10*scalar(jt20)
        pdum1=sdum1+0.5_lk*wt01*wt11*(sdum1-((1.0_lk+wt11)*scalar(jt31)+(1.0_lk+wt01)*scalar(jt01))/3.0_lk)
        pdum0=sdum0+0.5_lk*wt00*wt10*(sdum0-((1.0_lk+wt10)*scalar(jt30)+(1.0_lk+wt00)*scalar(jt00))/3.0_lk)
        if(i==0)then
          vector(1,1,ij)=h1*(pdum1-scalar(ij))
        elseif(i==mpsi)then
          vector(1,1,ij)=h0*(scalar(ij)-pdum0)
        else
        vector(1,1,ij)=wp1*h1*pdum1-wp0*h0*pdum0+scalar(ij)*(wp0*h0-wp1*h1)
        endif
     enddo
  enddo


! d_scalar/d_theta

#ifdef _OPENACC
  !$acc parallel loop gang
#else
  !$omp parallel do private(i,j,ij,jt)
#endif
  do i=ibin,ibout
     !$acc loop vector
     do j=1,mtheta(i)
        ij=igrid(i)+j
        jt=j+1-mtheta(i)*(j/mtheta(i))
        vector(2,1,ij)=difft(i)*(scalar(igrid(i)+jt)-scalar(igrid(i)+j-1))
     enddo
  enddo


#ifdef TIME_MPI
  call timer_start()
#endif
! send scalar to right and receive from left
  !$acc kernels
  sendr=scalar(:)
  !$acc end kernels
  recvl=0.0_lk
  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 scalar to left and receive from right
  !$acc kernels
  sendl=scalar(:)
  !$acc end kernels
  recvr=0.0_lk
  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
  gradient_sendrecv_time=gradient_sendrecv_time+timer_get_et()
  gradient_sendrecv_count=gradient_sendrecv_count+2
  gradient_sendrecv_size=gradient_sendrecv_size+4.0*icount*lk
#endif

! unpack scalar_boundary and calculate E_zeta at boundaries toroidal
#ifndef _OPENACC
!$omp parallel do private(i,j,ii,jt,ij,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

! d_scalar/d_zeta
     do j=1,mtheta(i)
        ij=igrid(i)+j
        vector(3,1,ij)=(pright(j)-pleft(j))*diffz
     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
        vector(3,1,ij)=(pright-pleft)*diffz
     enddo
  enddo
#endif

! adjust the difference between safety factor q and qtinv for fieldline coordinate
#ifdef _OPENACC
  !$acc parallel loop gang
#else    
  !$omp parallel do private(i,j,q,delq,ij)
#endif   
  do i=1,mpsi-1
     q=qmesh(i)
#ifdef _FRC
        delq=0.0_lk
#else
        delq=(1.0_lk/q-qtinv(i))
#endif
     !$acc loop vector
     do j=1,mtheta(i)
        ij=igrid(i)+j
        vector(3,:,ij)=vector(3,:,ij)+delq*vector(2,:,ij)
     enddo
  enddo

  !$acc kernels
  tmp_vector1 = vector(1,:,:)
  tmp_vector2 = vector(2,:,:)
  tmp_vector3 = vector(3,:,:)
  !$acc end kernels


  call periodicity_gpu(tmp_vector1)
  call periodicity_gpu(tmp_vector2)
  call periodicity_gpu(tmp_vector3)


  !$acc kernels
  vector(1,:,:) = tmp_vector1
  vector(2,:,:) = tmp_vector2
  vector(3,:,:) = tmp_vector3
  !$acc end kernels

end subroutine gradient_t
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine grad2para(scalarin,scalarout)
  use system_env
  use precision, only: lk, mpi_Rsize
  use global_parameters
  use field_array
  implicit none

  integer i,ii,ij,j,k,icount,idest,isource,isendtag,irecvtag,jt,&
      istatus(MPI_STATUS_SIZE)
  real(lk) difft(0:mpsi),diffz,&
      sendl(mgrid),recvr(mgrid),sendr(mgrid),recvl(mgrid),q,delq,d2dt2(mgrid)
#ifdef _OPENACC
  real(lk) pleft,pright
#else
  real(lk)pleft(mthetamax),pright(mthetamax)
#endif
  real(lk),dimension(mgrid),intent(in) :: scalarin
  real(lk),dimension(0:1,mgrid),intent(out) :: scalarout
  !$acc declare create(difft,sendl, recvr, sendr, recvl, d2dt2)

! finite difference for e-field in equilibrium unit
  !$acc kernels
  difft=1.0_lk/deltat
  !$acc end kernels
  diffz=1.0_lk/deltaz
  !$acc kernels
  d2dt2=0.0_lk
  !$acc end kernels

#ifdef _OPENACC
  !$acc parallel loop gang vector
#else
  !$omp parallel do private(i,j,k)
#endif
  do i=1,mgrid
     do j=0,1
          scalarout(j,i)=0.0_lk
     enddo
  enddo

! d2_scalar/d_theta2
#ifdef _OPENACC
  !$acc parallel loop gang
#else
  !$omp parallel do private(i,j,ij,jt)
#endif
  do i=1,mpsi-1
     !$acc loop vector
     do j=1,mtheta(i)
        ij=igrid(i)+j
        jt=j+1-mtheta(i)*(j/mtheta(i))
        d2dt2(ij)=difft(i)*difft(i)*(scalarin(igrid(i)+jt)+scalarin(igrid(i)+j-1)-2.0_lk*scalarin(ij))
     enddo
  enddo

! send scalar to right and receive from left
  !$acc kernels
  sendr=scalarin(:)
  recvl=0.0_lk
  !$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 scalar to left and receive from right
  !$acc kernels
  sendl=scalarin(:)
  recvr=0.0_lk
  !$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

! unpack scalar_boundary and calculate E_zeta at boundaries toroidal
#ifndef _OPENACC
!$omp parallel do private(i,j,ii,jt,ij,pleft,pright,q,delq)
  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

! d2_scalar/d_zeta2
     do j=1,mtheta(i)
        ij=igrid(i)+j
        scalarout(1,ij)=(pright(j)+pleft(j)-2.0_lk*scalarin(ij))*diffz*diffz
     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
        scalarout(1,ij)=(pright+pleft-2.0_lk*scalarin(ij))*diffz*diffz
     enddo
  enddo
#endif

! adjust the difference between safety factor q and qtinv for fieldline coordinate
#ifdef _OPENACC
  !$acc parallel loop gang
#else
  !$omp parallel do private(i,j,q,delq,ij)
#endif
  do i=1,mpsi-1
     q=qmesh(i)
     delq=(1.0_lk/(q*q)-qtinv(i)*qtinv(i))
     !$acc loop vector
     do j=1,mtheta(i)
        ij=igrid(i)+j
        scalarout(1,ij)=scalarout(1,ij)+delq*d2dt2(ij)
     enddo
  enddo
  call periodicity_gpu(scalarout)

end subroutine grad2para
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine grad2psi(scalarin,scalarout)
  use global_parameters
  use field_array
  implicit none

  integer i,ij,j
  real(lk) h0,h1,wp0,wp1,wt10,wt00,wt11,wt01,sdum0,sdum1,pdum0,pdum1
  integer jm,jp,jt11,jt21,jt31,jt01,jt10,jt20,jt30,jt00
  real(lk),dimension(mgrid),intent(in) :: scalarin
  real(lk),dimension(0:1,mgrid),intent(out) :: scalarout
!! d2_scalar/d_psi2 higher order
#ifdef _OPENACC
  !$acc parallel loop gang
#else
  !$omp parallel do private(i,j,ij,h0,h1,wp1,wp0,jm,jp,wt10,wt00,wt11,wt01,&
  !$omp& jt11,jt21,jt31,jt01,jt10,jt20,jt30,jt00,sdum1,sdum0,pdum1,pdum0)
#endif
  do i=1,mpsi-1

     h0=1.0_lk/(psimesh(i)-psimesh(i-1))
     h1=1.0_lk/(psimesh(i+1)-psimesh(i))
     wp1=h1/(h0+h1)
     wp0=1.0_lk-wp1

     jm=igrid(i-1)
     jp=igrid(i+1)

     !$acc loop vector
     do j=1,mtheta(i)
        ij=igrid(i)+j
        wt10=wtp1(2,ij)   !upper poloidal grid on inner flux surface
        wt00=1.0_lk-wt10       !lower poloidal grid on inner flux surface
        wt11=wtp1(1,ij)   !upper poloidal grid on outer flux surface
        wt01=1.0_lk-wt11       !lower poloidal grid on outer flux surface

        jt11=jtp1(1,ij)   !!lower poloidal grid on outer flux surface
        jt21=jt11+1
        jt31=jt21+1-mtheta(i+1)*((jt21-jp)/mtheta(i+1))
        jt01=jp+mod(jt11-jp-1+mtheta(i+1),mtheta(i+1))

        jt10=jtp1(2,ij)   !!lower poloidal grid on the innoer flux surface
        jt20=jt10+1
        jt30=jt20+1-mtheta(i-1)*((jt20-jm)/mtheta(i-1))
        jt00=jm+mod(jt10-jm-1+mtheta(i-1),mtheta(i-1))

        sdum1=wt01*scalarin(jt11)+wt11*scalarin(jt21)
        sdum0=wt00*scalarin(jt10)+wt10*scalarin(jt20)
        pdum1=sdum1+0.5_lk*wt01*wt11*(sdum1-((1.0_lk+wt11)*scalarin(jt31)+(1.0_lk+wt01)*scalarin(jt01))/3.0_lk)
        pdum0=sdum0+0.5_lk*wt00*wt10*(sdum0-((1.0_lk+wt10)*scalarin(jt30)+(1.0_lk+wt00)*scalarin(jt00))/3.0_lk)
        scalarout(1,ij)=2.0_lk*h0*h1*(wp1*pdum1+wp0*pdum0-scalarin(ij))
     enddo
  enddo
  call periodicity_gpu(scalarout)
end subroutine grad2psi
!!!!!!!!!

subroutine laplacian_petsc(scalar,scalar_out)
  use precision
  use global_parameters,only: mgrid,mpsi,rho0,gtcout,mype
  use particle_array,only: aion,qion,meshti,meshni
  use field_array,only: igrid,mtheta,bmesh
  use petsc_array,only: userb,userx
#ifdef _PETSc
#include "petscversion.h"
#if PETSC_VERSION_GT(3,8,0)
  use petsc_array,only: petsc_ps_ilap0
#endif
#endif
  implicit none

  integer :: i,ij,j
  !real(lk),dimension(mpsi+1) :: ddum
  real(lk),dimension(mgrid),intent(in) :: scalar
  real(lk),dimension(mgrid),intent(out) :: scalar_out

#ifdef _PETSc
#if PETSC_VERSION_GT(3,8,0)
  interface 
    subroutine laplacian_core(psdata,luserx_output,luserb_input)
       use precision, only: lk
       use petsc_array
       implicit none
    
       type (petscData),intent(in) :: psdata
       real(lk),dimension(0:),intent(inout) :: luserx_output
       real(lk),dimension(0:),intent(in) :: luserb_input
    end subroutine laplacian_core
  end interface
#endif
#endif

  scalar_out=0.0_lk

!$omp parallel do private(i)
  do i=1,mgrid
     userb(i-1)=scalar(i)
  enddo

#ifdef _PETSc
#if PETSC_VERSION_GT(3,8,0)
  call laplacian_core(petsc_ps_ilap0,userx,userb)
#else
  call laplacian_core
#endif
#endif

#ifdef _AMGX
  if (mype==0) write(gtcout,*)'AmgX solver not implemented for this case'
  stop
#endif
#ifdef _HYPRE
  if (mype==0) write(gtcout,*)'HYPRE solver not implemented for this case'
  stop
#endif

!$omp parallel do private(i,j,ij)
  do i=1,mpsi-1
     do j=0,mtheta(i)
        ij=igrid(i)+j
        scalar_out(ij)=-userx(ij-1)*bmesh(ij)*bmesh(ij)/(aion*rho0*rho0*meshni(i))
     enddo
  enddo

!*Boundary conditions at the inner wall and at the outer wall
!$omp parallel do private(i,j,ij)
  do i=0,mpsi,mpsi
    do j=0,mtheta(i)
      ij=igrid(i)+j
      scalar_out(ij)=0.0_lk
    enddo
  enddo

end subroutine laplacian_petsc
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine laplacian_gpu(scalar, scalar_out, none0bound)
  use precision
  use global_parameters, only: mgrid
  implicit none

  integer,optional, intent(in) :: none0bound
  real(lk),dimension(mgrid),intent(in) :: scalar
  real(lk),dimension(mgrid),intent(out) :: scalar_out
interface
subroutine laplacian_t(scalar,scalar_out,none0bound)
  use system_env
  use global_parameters, only:mgrid
  use precision
  implicit none
  integer, optional, intent (in) :: none0bound
  real(lk),dimension(mgrid),intent(in) :: scalar
  real(lk),dimension(mgrid),intent(out) :: scalar_out
end subroutine laplacian_t
end interface

  if(present(none0bound))then
    call laplacian_t(scalar, scalar_out,none0bound)
  else
    call laplacian_t(scalar, scalar_out)
  endif

end subroutine laplacian_gpu
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine laplacian_t(scalar,scalar_out,none0bound)
  use precision
  use global_parameters,only: mgrid,mpsilow,mpsihigh,mgridlow,mgridhigh,mpsi,mype,istep,&
   pi2,bcond,rg0!,rho0
 !! use particle_array,only: aion,qion,meshti
  use field_array,only: igrid,mtheta,lapmat,indexplap,nindexlap,&
   zeta1,qtinv,deltat,psimesh,deltar,mindexlap
  use interfaces,only: axisExtrapolate
  implicit none

  integer, optional, intent (in) :: none0bound


  integer :: i,ij,j,ii,jj,inone0bound
  real(lk),dimension(mgrid),intent(in) :: scalar
  real(lk),dimension(mgrid),intent(out) :: scalar_out
  real(lk) tdum,rdum,ave,sinave,cosave,ave2,sinave2,cosave2

  if(PRESENT(none0bound))then
      inone0bound=none0bound
  else
      inone0bound=0
  endif

#ifdef _OPENACC
  !$acc parallel loop gang vector present(scalar)
#else   
  !$omp parallel do private(i,j)
#endif  
     do i=mgridlow,mgridhigh
        scalar_out(i)=0.0_lk
        !$acc loop seq
        do j=1,nindexlap(i)
           scalar_out(i)=scalar_out(i)+lapmat(j,i)*scalar(indexplap(j,i))
        enddo
     enddo

!For linear boundary conditions, Fourier decomposes mpsilow+1 flux surface and linear exterpolates
!solution to the magnetic axis.
  if(bcond==1)then
    call axisExtrapolate(scalar_out(:))
!Zero inner radial boundary condition for bcond=0
  else
    !$acc parallel loop gang vector
    do j=0,mtheta(0)
      ij= igrid(0)+j
      scalar_out(ij)=0.0_lk
    enddo
  endif

!Outer radial boundary Condition
  !$acc parallel loop gang vector
   do j=0,mtheta(mpsi)
      ij=igrid(mpsi)+j
      scalar_out(ij)=0.0_lk
   enddo


!correction for the mpihigh flux surface for nonzero outer boundary
   if(inone0bound > 0)then
     !$acc parallel loop gang vector present(scalar)
      do j=0, mtheta(mpsihigh)
         ii = igrid(mpsihigh) + j
         do jj=nindexlap(ii)+1,11
            scalar_out(ii)=scalar_out(ii)+lapmat(jj,ii)*scalar(indexplap(jj,ii))
         enddo
      enddo
   endif

end subroutine laplacian_t
!======================================================================
subroutine filter_keep_zonal_cpu(scalar,scalar00)
  use global_parameters, only: mgrid,mpsi
  use precision, only: lk
#ifdef _OPENACC
  use openacc
#endif
  implicit none
  logical tmp
  real(lk),dimension(0:1,mgrid),intent(inout) :: scalar
  real(lk),dimension(0:mpsi),intent(out) :: scalar00

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

  call filter_keep_zonal_t(scalar,scalar00)

#ifdef _OPENACC
  if(tmp)then
    !$acc exit data copyout(scalar, scalar00)
    !$acc exit data delete(scalar, scalar00)
  else
    !$acc update host(scalar,scalar00)
  endif
#endif
end subroutine filter_keep_zonal_cpu
!======================================================================
subroutine filter_keep_zonal_gpu(scalar,scalar00)
  call filter_keep_zonal_t(scalar,scalar00)
end subroutine filter_keep_zonal_gpu
!======================================================================
subroutine filter_keep_zonal_t(scalar,scalar00)
  use precision
  use system_env
  use global_parameters,only: mpsi,mgrid
  use field_array,only: igrid,mtheta
  implicit none


  integer i,j,ij
  real(lk),dimension(0:1,mgrid),intent(inout) :: scalar
  real(lk),dimension(0:mpsi),intent(out) :: scalar00
  real(lk),dimension(mgrid) :: tmp_scalar

  !$acc declare create(tmp_scalar)
  !$acc kernels
  scalar00=0.0_lk
  !$acc end kernels

  !$acc kernels
  tmp_scalar = scalar(1,:)
  !$acc end kernels
  call zonal_field_gpu(tmp_scalar,scalar00)

#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
        scalar(0,ij)=scalar(0,ij)-scalar00(i)
        scalar(1,ij)=scalar(1,ij)-scalar00(i)
     enddo
     scalar(0,igrid(i))=scalar(0,igrid(i)+mtheta(i))
     scalar(1,igrid(i))=scalar(1,igrid(i)+mtheta(i))
  enddo
  CALL FILTER_GPU(scalar)
  call radial_zonal_smooth_gpu(scalar00)

#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
         scalar(0,ij)=scalar(0,ij)+scalar00(i)
         scalar(1,ij)=scalar(1,ij)+scalar00(i)
      enddo
      scalar(0,igrid(i))=scalar(0,igrid(i)+mtheta(i))
      scalar(1,igrid(i))=scalar(1,igrid(i)+mtheta(i))
   enddo

end subroutine filter_keep_zonal_t

!======================================================================
subroutine zonal_field_gpu(scalar,scalar00)
  call zonal_field_t(scalar,scalar00)
end subroutine zonal_field_gpu
!======================================================================
subroutine zonal_field_cpu(scalar,scalar00)
use global_parameters, only: mgrid,mpsi
  use precision, only: lk
#ifdef _OPENACC
  use openacc
#endif
  implicit none
  logical :: tmp
  real(lk),dimension(mgrid), intent(in) :: scalar
  real(lk),dimension(0:mpsi), intent(out) :: scalar00

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

  call zonal_field_t(scalar, scalar00)

#ifdef _OPENACC
  if(tmp)then
    !$acc exit data copyout(scalar00)
    !$acc exit data delete(scalar,scalar00)
  else
    !$acc update host(scalar00)
  endif
#endif

end subroutine zonal_field_cpu
!======================================================================
subroutine zonal_field_t(scalar,scalar00)
  use precision
  use system_env
  use global_parameters,only: mpsi,mtoroidal,toroidal_comm,mgrid
  use field_array,only: igrid,mtheta,jmesh
  use particle_array,only: jacobianpsi
  implicit none


  integer i,j,ij
  real(lk),dimension(0:mpsi) :: adum
  real(lk),dimension(mgrid),intent(in) :: scalar
  real(lk),dimension(0:mpsi),intent(out) :: scalar00
  !$acc declare create(adum)
  !$acc kernels
  scalar00=0.0_lk
  !$acc end kernels

!calculate zonal <scalar>
#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
        !$acc atomic update
        scalar00(i)=scalar00(i)+scalar(ij)/jmesh(ij)
     enddo
     scalar00(i)=scalar00(i)/jacobianpsi(i)/real(mtheta(i))
  enddo

  adum=0.0_lk
#ifdef CUDA_MPI
  !$acc host_data use_device(scalar00,adum)
#else
  !$acc update host(scalar00)
#endif
  call MPI_ALLREDUCE(scalar00,adum,mpsi+1,mpi_Rsize, MPI_SUM,toroidal_comm,ierror)
#ifdef CUDA_MPI
  !$acc end host_data
#else
  !$acc update device(adum)
#endif

#ifdef _OPENACC
  !$acc kernels
#endif
  scalar00=adum/real(mtoroidal)
#ifdef _OPENACC
  !$acc end kernels
#endif

end subroutine zonal_field_t
!====================================================================
subroutine radial_zonal_smooth_gpu(scalar00)
  call radial_zonal_smooth_t(scalar00)
end subroutine radial_zonal_smooth_gpu

subroutine radial_zonal_smooth_cpu(scalar00)
  use global_parameters, only: mpsi
  use precision, only: lk
#ifdef _OPENACC
  use openacc
#endif
  implicit none
  logical tmp
  real(lk),dimension(0:mpsi),intent(inout) :: scalar00

#ifdef _OPENACC
  if(acc_is_present(scalar00))then
    !$acc update device(scalar00)
    tmp = .false.
  else
    tmp = .true.
    !$acc enter data copyin(scalar00)
  endif
#endif
  call radial_zonal_smooth_t(scalar00)
#ifdef _OPENACC
  if(tmp)then
    !$acc exit data copyout(scalar00)
    !$acc exit data delete(scalar00)
  else
    !$acc update host(scalar00)
  endif
#endif

end subroutine radial_zonal_smooth_cpu

subroutine radial_zonal_smooth_t(scalar00)
  use precision
  use global_parameters,only: mpsi
  implicit none


  integer i
  real(lk),dimension(0:mpsi) :: farray00tmp
  real(lk),dimension(0:mpsi),intent(inout) :: scalar00
  !$acc declare create(farray00tmp, scalar00)

!! (-0.0625 0.25 0.625 0.25 -0.0625) radial smoothing of (0,0) mode
  do i=1,2
     !$acc kernels
     farray00tmp(0)=scalar00(0)
     farray00tmp(mpsi)=scalar00(mpsi)
     farray00tmp(1)=scalar00(3)
     farray00tmp(mpsi-1)=scalar00(mpsi-3)
     farray00tmp(2:mpsi-2)=scalar00(0:mpsi-4)+scalar00(4:mpsi)
     farray00tmp(1:mpsi-1)=0.625_lk*scalar00(1:mpsi-1)+0.25_lk*(scalar00(0:mpsi-2)+scalar00(2:mpsi))-&
          0.0625_lk*farray00tmp(1:mpsi-1)
     scalar00=farray00tmp
     !$acc end kernels
  enddo

end subroutine radial_zonal_smooth_t

!====================================================================
subroutine gradient_psi(scalar00,scalarp00)
  use precision
  use system_env
  use global_parameters,only: mpsi,mgrid
  use field_array,only: igrid,mtheta
  implicit none


  integer i,j,ij
  real(lk),dimension(mgrid) :: scalar00_tmp
  real(lk),dimension(3,0:1,mgrid) :: gradscalar00_tmp
  !$acc declare create(scalar00_tmp,gradscalar00_tmp)
  real(lk),dimension(0:mpsi),intent(in) :: scalar00
  real(lk),dimension(0:mpsi),intent(out) :: scalarp00

  !$acc kernels
  scalarp00=0.0_lk
  scalar00_tmp=0.0_lk
  gradscalar00_tmp=0.0_lk
  !$acc end kernels

#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
        scalar00_tmp(ij)=scalar00(i)
     enddo
  enddo
  call gradient_gpu(scalar00_tmp(:),gradscalar00_tmp)

#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
        scalarp00(i)=gradscalar00_tmp(1,1,ij)
     enddo
  enddo

end subroutine gradient_psi
