get_ediffuse_matrix Subroutine

  subroutine get_ediffuse_matrix (aa, bb, cc, ig, ik, it, il, is, xe)
    use species, only: spec
    use theta_grid, only: bmag
    use le_grids, only: al, negrid, w=>w_maxw
    use kt_grids, only: kperp2
    use gs2_time, only: code_dt, tunits

    implicit none

    integer, intent (in) :: ig, ik, it, il, is
    real, dimension (:), intent (in) :: xe
    real, dimension (:), intent (out) :: aa, bb, cc

    integer :: ie
    real :: vn, slb1, xe0, xe1, xe2, xel, xer, capgr, capgl, ee

    vn = vnmult(2) * spec(is)%vnewk * tunits(ik)

    slb1 = safe_sqrt(1.0 - bmag(ig)*al(il))     ! xi_j

    select case (ediff_switch)
    case (ediff_scheme_default)
       do ie = 2, negrid - 1
          xe0 = xe(ie-1)
          xe1 = xe(ie)
          xe2 = xe(ie+1)

          xel = (xe0 + xe1) * 0.5
          xer = (xe1 + xe2) * 0.5

          capgr = capg(xer)
          capgl = capg(xel)

          ee = 0.125 * (1 - slb1**2) * vnew_s(ik,ie,is) * kperp2(ig,it,ik)*cfac &
               / (bmag(ig) * spec(is)%zstm)**2

          ! coefficients for tridiagonal matrix:
          cc(ie) = -0.25 * vn * code_dt * capgr / (w(ie) * (xe2 - xe1))
          aa(ie) = -0.25 * vn * code_dt * capgl / (w(ie) * (xe1 - xe0))
          bb(ie) = 1.0 - (aa(ie) + cc(ie)) + ee*code_dt
       end do

       ! boundary at v = 0
       xe1 = xe(1)
       xe2 = xe(2)
       xer = (xe1 + xe2) * 0.5

       capgr = capg(xer)

       ee = 0.125 * (1 - slb1**2) * vnew_s(ik,1,is) * kperp2(ig,it,ik)*cfac &
            / (bmag(ig) * spec(is)%zstm)**2

       cc(1) = -0.25 * vn * code_dt * capgr / (w(1) * (xe2 - xe1))
       aa(1) = 0.0
       bb(1) = 1.0 - cc(1) + ee * code_dt

       ! boundary at v = infinity
       xe0 = xe(negrid-1)
       xe1 = xe(negrid)

       xel = (xe1 + xe0) * 0.5

       capgl = capg(xel)

       ee = 0.125 * (1.-slb1**2) * vnew_s(ik,negrid,is) * kperp2(ig,it,ik) * cfac &
            / (bmag(ig) * spec(is)%zstm)**2

       cc(negrid) = 0.0
       aa(negrid) = -0.25 * vn * code_dt * capgl / (w(negrid) * (xe1 - xe0))
       bb(negrid) = 1.0 - aa(negrid) + ee*code_dt

    case (ediff_scheme_old)

       ! non-conservative scheme
       do ie = 2, negrid-1
          xe0 = xe(ie-1)
          xe1 = xe(ie)
          xe2 = xe(ie+1)

          xel = (xe0 + xe1) * 0.5
          xer = (xe1 + xe2) * 0.5

          capgr = capg_old(xe1, xer)
          capgl = capg_old(xe1, xel)

          ee = 0.125 * (1 - slb1**2) * vnew_s(ik,ie,is) * kperp2(ig,it,ik) * cfac &
               / (bmag(ig) * spec(is)%zstm)**2

          ! coefficients for tridiagonal matrix:
          cc(ie) = -vn * code_dt * capgr / ((xer-xel) * (xe2 - xe1))
          aa(ie) = -vn * code_dt * capgl / ((xer-xel) * (xe1 - xe0))
          bb(ie) = 1.0 - (aa(ie) + cc(ie)) + ee * code_dt
       end do

       ! boundary at xe = 0
       xe1 = xe(1)
       xe2 = xe(2)
       xer = (xe1 + xe2) * 0.5

       capgr = capg_old(xe1, xer)

       ee = 0.125 * (1 - slb1**2) * vnew_s(ik,1,is) * kperp2(ig,it,ik)*cfac &
            / (bmag(ig) * spec(is)%zstm)**2

       cc(1) = -vn * code_dt * capgr / (xer * (xe2 - xe1))
       aa(1) = 0.0
       bb(1) = 1.0 - cc(1) + ee * code_dt

       ! boundary at xe = 1
       xe0 = xe(negrid-1)
       xe1 = xe(negrid)
       xel = (xe1 + xe0) * 0.5

       capgl = capg_old(xe1, xel)

       ee = 0.125 * (1 - slb1**2) * vnew_s(ik,negrid,is) * kperp2(ig,it,ik)*cfac &
            / (bmag(ig) * spec(is)%zstm)**2

       cc(negrid) = 0.0
       aa(negrid) = -vn * code_dt * capgl / ((1.0-xel) * (xe1 - xe0))
       bb(negrid) = 1.0 - aa(negrid) + ee * code_dt
    end select

  contains
    elemental real function capg_kernel(xe)
      use constants, only: sqrt_pi
      real, intent(in) :: xe
      ! This is the same as hsg_func(xe) * 2 * xe**2
      capg_kernel =  erf(xe) - 2 * xe * exp(-xe**2) / sqrt_pi
    end function capg_kernel

    elemental real function capg(xe)
      use constants, only: sqrt_pi
      real, intent(in) :: xe
      capg = 2 * exp(-xe**2) * capg_kernel(xe) / (xe * sqrt_pi)
    end function capg

    elemental real function capg_old(xeA, xeB)
      real, intent(in) :: xeA, xeB
      capg_old = (0.5 * exp(xeA**2-xeB**2) / xeA**2) * capg_kernel(xeB) / xeB
    end function capg_old
  end subroutine get_ediffuse_matrix