get_flux_vs_e Subroutine

  subroutine get_flux_vs_e (g_in, fld, flx)
    use theta_grid, only: ntgrid, grho, theta, jacob
    use kt_grids, only: aky
    use le_grids, only: negrid, wl
    use species, only: nspec
    use mp, only: sum_reduce, proc0
    use gs2_layouts, only: g_lo,ie_idx,il_idx,is_idx,it_idx,ik_idx,isign_idx
    use integration, only: trapezoidal_integration
    use gs2_layouts, only: g_lo
    implicit none
    !> Input weighted distribution
    complex, dimension (-ntgrid:,:,g_lo%llim_proc:), intent (in) :: g_in
    complex, dimension (-ntgrid:,:,:), intent (in) :: fld
    real, dimension (:,:), intent (inout) :: flx
    real, dimension (:,:), allocatable :: total
    real:: wgt,fac
    integer :: ik, it, is, iglo, ie, il

    allocate(total(negrid,nspec))
    total = 0.

    ! This is essentially 2 pi / surfarea
    wgt = 1.0/trapezoidal_integration(theta, grho * jacob)

    do iglo = g_lo%llim_proc, g_lo%ulim_proc
       fac =0.5
       ie = ie_idx(g_lo,iglo)
       il = il_idx(g_lo,iglo)
       is = is_idx(g_lo,iglo)
       it = it_idx(g_lo,iglo)
       ik = ik_idx(g_lo,iglo)
       if (aky(ik) == 0.) fac = 1.0

       total(ie,is) = total(ie,is) + fac*&
            trapezoidal_integration(theta, &
            aimag(g_in(:,1,iglo)*conjg(fld(:,it,ik))) &
            * wl(:,il) * jacob) * aky(ik) * wgt
       total(ie,is) = total(ie,is) + fac*&
            trapezoidal_integration(theta, &
            aimag(g_in(:,2,iglo)*conjg(fld(:,it,ik))) &
            * wl(:,il) * jacob) * aky(ik) * wgt
    end do

    call sum_reduce(total,0)

    if (proc0) flx = 0.5*total

  end subroutine get_flux_vs_e