!> Purpose: calculaing delta B_para from particle information.
!> Input: string indicating whether non-adiabatic electron response is included
!in the calculation.
!"kinetic": with non-adiabatic electron response. "adiabatic": w/o non-adiabatic
!electron response 
!> Output: No output; just updates the array bpara.
!> Options: inversing=0 assumes k_perp*rho_s<<1. inversing=1 uses the Pade-like
!> approximation (Sec IV, G.Dong et al., POP, 2017) 
subroutine bparasolver(scenario)
  use system_env
  use precision
  use global_parameters
  use field_array
  use particle_array
  use petsc_array,only: luserb_bpara,luserx_bpara
  use equilibrium
  use interfaces
#ifdef _PETSc
#include "petscversion.h"
#if PETSC_VERSION_GT(3,8,0)
  use petsc_array,only: petsc_bpara
#endif
#endif
  implicit none

  integer m,i,j,ij,inversing
  real(lk) dtime,pdum,tdum,dpx,dp2,dtx,dt2,dbdt,dbdp,dbdz,&
           ri,g,b_inv,q_inv,&
           dentmp_bpara(mgrid),&
           dptdp,dptdt,dptdz,bpara0,&
           ne,te,dnedp,dtedp,er,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,&
           dkpparadpi,dkpparadti,dkpparadzi,dkpperpdpi,dkpperpdti,dkpperpdzi,gamma00
  character(*),intent(in) :: scenario
#ifdef _PETSc
#if PETSC_VERSION_GT(3,8,0)
  interface
     subroutine lapmat_pssolver(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 lapmat_pssolver
  end interface
#endif
#endif

  !$acc declare create(dentmp_bpara)
  dnedp=0.0_lk
  dtedp=0.0_lk

  inversing=1  !0: assuming k_perp*rho_s<<1 & high beta (Ge Dong 2017)
               !1: using the pade-like approximation & low beta (invalid for fem=1 or fload=1) (Ge Dong 2017)
               !2: assuming k_perp*rho_s<=1 & low beta
  if(fload==1 .or. fem==1)inversing=2
#ifdef _OPENACC
  !$acc parallel loop gang
#else
  !$omp parallel do private(i,j,ij,bpara0)
#endif
  do i=0,mpsi
    !$acc loop vector
    do j=1,mtheta(i)
      ij=igrid(i)+j
      bpara0=bpara(1,ij)
      bpara(1,ij)=-betae*(pressureiperp(1,ij)*meshni(i)+sfluidne(1,ij)*meshte(i)-bpara0*meshne(i)*meshte(i)/rho0**2 &
                  -sdeltapsi(1,ij)*meshte(i)*meshne(i)*kapate(i))/2.0_lk
    enddo
  enddo
  if(scenario=="kinetic")then
#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
          bpara(1,ij)=bpara(1,ij)-betae*meshne(i)*(-densitye(1,ij)*meshte(i)+pressureeperp(1,ij))/2.0_lk
       enddo
    enddo
  endif

  if(inversing==0)then
#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
        bpara(1,ij)=bpara(1,ij)+betae*(densityi(1,ij)*meshni(i)*qion+sfluidne(1,ij)*qelectron)*meshte(1)
      enddo
    enddo
  elseif(inversing==1)then
#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
        dentmp_bpara(ij)=(densityi(1,ij)*meshni(i)*qion+sfluidne(1,ij)*qelectron)
        bpara(1,ij)=bpara(1,ij)+betae*dentmp_bpara(ij)*meshte(1)/2.0
      enddo
      dentmp_bpara(igrid(i)) = dentmp_bpara(igrid(i)+mtheta(i))
    enddo
  else
#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
         bpara(1,ij)=bpara(1,ij)-3.0_lk/4.0_lk*betae*meshni(i)*meshti(i)*aion/(qion*bmesh(ij)**2)*d2phi(ij)&
                     -5.0_lk/8.0_lk*betae*meshni(i)*meshti(i)**2*aion**2/(qion**3*bmesh(ij)**4)*d4phi(ij)
         if(fload==1)then
            bpara(1,ij)=bpara(1,ij)-betae*pressurefperp(1,ij)*meshnf(i)/2.0_lk &
                        -3.0_lk/4.0_lk*betae*meshnf(i)*meshtf(i)*afast/(qfast*bmesh(ij)**2)*d2phi(ij) &
                        -5.0_lk/8.0_lk*betae*meshnf(i)*meshtf(i)**2*afast**2/(qfast**3*bmesh(ij)**4)*d4phi(ij)
         endif
         bpara(1,ij)=bpara(1,ij)/bmesh(ij)
      enddo
      bpara(1,igrid(i)) = bpara(1,igrid(i)+mtheta(i))
    enddo
  endif

  if(inversing==1)then
    !$acc kernels
    if(psi0>0.0_lk)dentmp_bpara(igrid(0):igrid(0)+mtheta(0))=0.0_lk
    dentmp_bpara(igrid(mpsi):igrid(mpsi)+mtheta(mpsi))=0.0_lk
    !$acc end kernels
#ifdef _PETSc
#ifdef _OPENACC
       !$acc parallel loop gang vector
#else
       !$omp parallel do private(i)
#endif
    do i=mgridlow,mgridhigh
      luserb_bpara(i-mgridlow)=dentmp_bpara(i)
    enddo
    !$acc update host(luserb_bpara)
#if PETSC_VERSION_GT(3,8,0)
    call lapmat_pssolver(petsc_bpara,luserx_bpara,luserb_bpara)
#else
    call lapmat_pssolver
#endif
    !$acc update device(luserx_bpara)

#ifdef _OPENACC
    !$acc parallel loop gang vector
#else
    !$omp parallel do private(i)
#endif
     do i=mgridlow,mgridhigh
       bpara(1,i)=bpara(1,i)+betae*meshte(1)*luserx_bpara(i-mgridlow)/2.0
     enddo
#else
  if (mype==0) then
        write(gtcout,*) 'hypre and amgx are not ready for lapmat_pssolver_bpara'
        flush(gtcout)
        call GTC_ABORT
    endif
#endif

#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
        bpara(1,ij)=bpara(1,ij)*bmesh(ij)
      enddo
      bpara(1,igrid(i)) = bpara(1,igrid(i)+mtheta(i))
    enddo
  endif

  !$acc kernels
  bpara(:,igrid(0):igrid(0)+mtheta(0))=0.0
  bpara(:,igrid(mpsi):igrid(mpsi)+mtheta(mpsi))=0.0
  !$acc end kernels
  if(nfilter>0)CALL FILTER_GPU(bpara)
  call smooth_gpu(bpara)
end subroutine bparasolver
