init_bessel Subroutine

  subroutine init_bessel
    use dist_fn_arrays, only: aj0, aj1, aj0_gf, aj1_gf
    use kt_grids, only: kperp2
    use species, only: spec, nspec
    use theta_grid, only: ntgrid, bmag
    use le_grids, only: energy,  al,negrid, nlambda, forbid
    use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx, ie_idx, is_idx, gf_lo
    use spfunc, only: j0, j1
    implicit none
    integer :: ik, it, il, ie, is, iglo, igf
    real, dimension(-ntgrid:ntgrid) :: arg

    if (bessinit) return
    bessinit = .true.

    allocate (aj0(-ntgrid:ntgrid,g_lo%llim_proc:g_lo%ulim_alloc))
    allocate (aj1(-ntgrid:ntgrid,g_lo%llim_proc:g_lo%ulim_alloc))

    !$OMP PARALLEL DO DEFAULT(none) &
    !$OMP PRIVATE(iglo, ik, it, il, ie, is, arg) &
    !$OMP SHARED(g_lo, spec, energy, al, bmag, kperp2, aj0, aj1, forbid) &
    !$OMP SCHEDULE(static)
    do iglo = g_lo%llim_proc, g_lo%ulim_proc
       ik = ik_idx(g_lo,iglo)
       it = it_idx(g_lo,iglo)
       il = il_idx(g_lo,iglo)
       ie = ie_idx(g_lo,iglo)
       is = is_idx(g_lo,iglo)
       arg = spec(is)%bess_fac * spec(is)%smz * &
            sqrt(energy(ie,is) * al(il) * kperp2(:,it,ik) / bmag)
       where (forbid(:, il))
          aj0(:, iglo) = 0.0
          aj1(:, iglo) = 0.0
       elsewhere
          aj0(:, iglo) = j0(arg)
          aj1(:, iglo) = j1(arg)
       end where
    end do
    !$OMP END PARALLEL DO

    !AJ For the gf_lo integrate we pre-calculate some of the factors
    !used in the velocity space integration, rather than
    !calculating them on the fly.
    if(gf_lo_integrate) then
       allocate (aj0_gf(-ntgrid:ntgrid,nspec,negrid,nlambda,gf_lo%llim_proc:gf_lo%ulim_alloc))
       allocate (aj1_gf(-ntgrid:ntgrid,nspec,negrid,nlambda,gf_lo%llim_proc:gf_lo%ulim_alloc))

       do igf = gf_lo%llim_proc, gf_lo%ulim_proc
          ik = ik_idx(gf_lo,igf)
          it = it_idx(gf_lo,igf)
          do il = 1,gf_lo%nlambda
             do ie = 1,gf_lo%negrid
                do is = 1,gf_lo%nspec
                   arg = spec(is)%bess_fac * spec(is)%smz * &
                        sqrt(energy(ie,is) * al(il) * kperp2(:,it,ik) / bmag)
                   where (forbid(:, il))
                      aj0_gf(:, is, ie, il, igf) = 0.0
                      aj1_gf(:, is, ie, il, igf) = 0.0
                   else where
                      aj0_gf(:, is, ie, il, igf) = j0(arg)
                      aj1_gf(:, is, ie, il, igf) = j1(arg)
                   end where
                end do
             end do
          end do
       end do
    end if
  end subroutine init_bessel