init_le_bessel Subroutine

  subroutine init_le_bessel
    use gs2_layouts, only: g_lo, le_lo, ig_idx
    use dist_fn_arrays, only: aj0, aj1, vpa
    use le_grids, only: negrid, g2le, ixi_to_il, energy => energy_maxw, al
    use theta_grid, only: ntgrid
    use redistribute, only: gather, scatter
    use array_utils, only: zero_array
    implicit none
    complex, dimension (:,:,:), allocatable :: ctmp, z_big
    integer :: ile, ig, ie, ixi, il

    allocate (z_big(-ntgrid:ntgrid, 2, g_lo%llim_proc:g_lo%ulim_alloc))
    allocate (ctmp(nxi_lim, negrid+1, le_lo%llim_proc:le_lo%ulim_alloc))
    ! We need to initialise ctmp as it is used as receiving buffer in
    ! g2le redistribute, which doesn't populate all elements
    call zero_array(ctmp)

    ! next set aj0le & aj1l
    z_big(:,1,:) = cmplx(aj0,aj1)
    z_big(:,2,:) = z_big(:,1,:)

    call gather (g2le, z_big, ctmp)

    if (.not. allocated(aj0le)) then
       allocate (aj0le(nxi_lim, negrid+1, le_lo%llim_proc:le_lo%ulim_alloc))
       allocate (vperp_aj1le(nxi_lim, negrid+1, le_lo%llim_proc:le_lo%ulim_alloc))
    end if

    aj0le = real(ctmp)
    vperp_aj1le = aimag(ctmp) !< Currently just aj1

    !$OMP PARALLEL DO DEFAULT(none) &
    !$OMP PRIVATE(ile, ig, ixi, ie, il) &
    !$OMP SHARED(le_lo, negrid, nxi_lim, ixi_to_il, vperp_aj1le, energy, al) &
    !$OMP SCHEDULE(static)
    do ile = le_lo%llim_proc, le_lo%ulim_proc
       ig = ig_idx(le_lo, ile)
       do ie = 1, negrid
          do ixi = 1, nxi_lim
             il = ixi_to_il(ixi, ig)
             vperp_aj1le(ixi, ie, ile) = vperp_aj1le(ixi, ie, ile) * energy(ie) * al(il)
          end do
       end do
    end do
    !$OMP END PARALLEL DO

    z_big = vpa
    call gather (g2le, z_big, ctmp)
    deallocate(z_big)
    if (.not. allocated(vpa_aj0_le)) then
       allocate (vpa_aj0_le(nxi_lim, negrid+1, le_lo%llim_proc:le_lo%ulim_alloc))
    end if
    vpa_aj0_le = real(ctmp) * aj0le
  end subroutine init_le_bessel