solfp_lorentz_le_layout Subroutine

  subroutine solfp_lorentz_le_layout (gle, diagnostics)
    use le_grids, only: jend, ng2, negrid, integrate_moment, il_is_wfb, grid_has_trapped_particles
    use gs2_layouts, only: ig_idx, ik_idx, is_idx, it_idx
    use run_parameters, only: ieqzip
    use gs2_layouts, only: le_lo
    use kt_grids, only: kwork_filter
    use array_utils, only: zero_array
    implicit none
    complex, dimension (:,:,le_lo%llim_proc:), intent (in out) :: gle
    integer, optional, intent (in) :: diagnostics
    complex, dimension (:), allocatable :: tmp
    complex, dimension (:,:,:), allocatable :: glec
    complex, dimension (nxi_lim) :: delta
    complex :: fac, gwfb
    integer :: ig, ik, il, is, je, it, ie, nxi_scatt, ile, cur_jend
    logical :: is_wfb

    if (heating .and. present(diagnostics)) then
       allocate (tmp(le_lo%llim_proc:le_lo%ulim_alloc)) ; call zero_array(tmp)
       allocate (glec(nxi_lim, negrid+1, le_lo%llim_proc:le_lo%ulim_alloc))
       call zero_array(glec)
       !$OMP PARALLEL DO DEFAULT(none) &
       !$OMP PRIVATE(ile, ig, je, ie, il, fac) &
       !$OMP SHARED(le_lo, jend, ng2, negrid, gle, glec, d1le) &
       !$OMP SCHEDULE(static)
       do ile = le_lo%llim_proc, le_lo%ulim_proc
          ig = ig_idx(le_lo,ile)

          je = 2*jend(ig)
          if (je == 0) then
             je = 2*ng2 
          end if

! when il=je-1 below, and we have trapped particles, gle is evaluated at gle(2*jend(ig),ie,ile).
! this seems like a bug, since there are only 2*jend(ig)-1 grid points and
! the value gle(2*jend(ig),ie,ile) corresponds to the value of g at xi = 0...this
! doesn't make any sense...MAB

          do ie = 1, negrid
             do il = 1, je-1
                fac = gle(il+1,ie,ile)-gle(il,ie,ile)
                glec(il,ie,ile) = conjg(fac)*fac*d1le(il,ie,ile)  ! d1le accounts for hC(h) entropy
             end do
          end do
       end do
       !$OMP END PARALLEL DO

       call integrate_moment (le_lo, glec, tmp)
       !$OMP PARALLEL DO DEFAULT(none) &
       !$OMP PRIVATE(ile, ig, it, ik, is) &
       !$OMP SHARED(le_lo, c_rate, tmp) &
       !$OMP SCHEDULE(static)
       do ile = le_lo%llim_proc, le_lo%ulim_proc
          ig = ig_idx(le_lo,ile)
          it = it_idx(le_lo,ile)
          ik = ik_idx(le_lo,ile)
          is = is_idx(le_lo,ile)
          c_rate(ig,it,ik,is,1) = tmp(ile)
       end do
       !$OMP END PARALLEL DO

       if (hyper_colls) then
          !$OMP PARALLEL DO DEFAULT(none) &
          !$OMP PRIVATE(ile, ig, je, ie, il, fac) &
          !$OMP SHARED(le_lo, jend, ng2, negrid, gle, glec, h1le) &
          !$OMP SCHEDULE(static)
          do ile = le_lo%llim_proc, le_lo%ulim_proc
             ig = ig_idx(le_lo,ile)

             je = 2*jend(ig)
             if (je == 0) then
                je = 2*ng2
             end if

             do ie = 1, negrid
                do il = 1, je-1
                   fac = gle(il+1,ie,ile)-gle(il,ie,ile)
                   glec(il,ie,ile) = conjg(fac)*fac*h1le(il,ie,ile)  ! h1le accounts for hH(h) entropy
                end do
             end do
          end do
          !$OMP END PARALLEL DO

          call integrate_moment (le_lo, glec, tmp)

          !$OMP PARALLEL DO DEFAULT(none) &
          !$OMP PRIVATE(ile, ig, it, ik, is) &
          !$OMP SHARED(le_lo, c_rate, tmp) &
          !$OMP SCHEDULE(static)
          do ile = le_lo%llim_proc, le_lo%ulim_proc
             ig = ig_idx(le_lo,ile)
             it = it_idx(le_lo,ile)
             ik = ik_idx(le_lo,ile)
             is = is_idx(le_lo,ile)
             c_rate(ig,it,ik,is,2) = tmp(ile)
          end do
          !$OMP END PARALLEL DO
       end if
       deallocate(tmp, glec)
    end if

    ! solve for gle row by row
    !$OMP PARALLEL DO DEFAULT(none) &
    !$OMP PRIVATE(ile, is, it, ik, ig, il, ie, je, nxi_scatt, &
    !$OMP gwfb, delta, is_wfb, cur_jend) &
    !$OMP SHARED(le_lo, kwork_filter, ieqzip, vnew, force_collisions, jend, ng2, special_wfb_lorentz, &
    !$OMP vpar_zero_mean, negrid, gle, qle, betaale, c1le) &
    !$OMP SCHEDULE(static)
    do ile = le_lo%llim_proc, le_lo%ulim_proc
       is = is_idx(le_lo, ile)
       ik = ik_idx(le_lo, ile)
       if ((abs(vnew(ik, 1, is)) < 2.0 * epsilon(0.0)) .and. .not. force_collisions) cycle
       it = it_idx(le_lo, ile)
       if (kwork_filter(it, ik)) cycle
       if (ieqzip(it, ik)) cycle
       ig = ig_idx(le_lo, ile)

       !CMRDDGC, 10/2/1014:
       ! Fixes for wfb treatment below, use same je definition in ALL cases
       !   je  = #physical xi values at location, includes duplicate point at vpar=0
       !  je-1 = #physical xi values removing duplicate vpar=0 point
       cur_jend = jend(ig)
       je = max(2 * cur_jend, 2 * ng2 + 1)
       nxi_scatt = je - 1
       is_wfb = il_is_wfb(cur_jend)
       if ((is_wfb .and. special_wfb_lorentz) .or. .not. grid_has_trapped_particles()) &
            nxi_scatt = nxi_scatt - 1

       do ie = 1, negrid
          if (je == 2 * cur_jend .and. is_wfb .and. vpar_zero_mean) then
             ! MRH if we have trapped particles and hence 2 vpar=0 points on the grid
             ! use the average of both vpar = 0 points in the lorentz diffusion in pitch angle
             ! this turns out to be necessary to suppress numerical instability
             ! when we handle pitch angle scattering physically at theta = +/- pi
             ! by setting special_wfb_lorentz =.false.
             ! Note : This can give large change in g_wfb when it is not treated as a trapped
             ! particle as the input g is unlikely to satisy g_wfb(sigma=1) ~ g_wfb(sigma=2).
             ! Note : We don't apply this treatment to other trapped particles as we assume
             ! g has a unique value at the bounce point for those particles.
             gle(cur_jend, ie, ile) = (gle(cur_jend, ie, ile) + &
                  gle(2 * cur_jend, ie, ile)) / 2.0
          end if

          ! zero redundant duplicate xi, isign=2 for vpar=0! This shouldn't be needed
          if (grid_has_trapped_particles()) gle(je, ie, ile) = 0.0d0

          if (is_wfb .and. special_wfb_lorentz) then
             !CMRDDGC:  special_wfb_lorentz = t  => unphysical handling of wfb at bounce pt:
             !          remove wfb from collisions, reinsert later
             !
             ! first save gwfb for reinsertion later
             ! Note we don't save both sigma values as we force g(v||=0)_+ == g(v||=0)_-
             ! at the end of the loop anyway.
             gwfb = gle(ng2 + 1, ie, ile)
             ! then remove vpa = 0 point, weight 0: (CMR confused by this comment!)
             !The above is referring to the conservative scheme coefficients which involve
             !1/pitch_weight but pitch_weight = 0 for the WFB at the WFB bounce point.
             !Special_wfb_lorentz is a way to avoid this issue by ignoring the WFB pitch angle
             gle(ng2+1:je-2, ie, ile) = gle(ng2+2:je-1, ie, ile)
             ! Zero out the gle value not overwritten in the above. Shouldn't be needed
             gle(je - 1, ie, ile) = 0.0
          end if

          ! right and left sweeps for tridiagonal solve:
          delta(1) = gle(1, ie, ile)
          do il = 1, nxi_scatt
             delta(il+1) = gle(il+1, ie, ile) - qle(il+1, ie, ile) * delta(il)
          end do

          gle(nxi_scatt+1, ie, ile) = delta(nxi_scatt+1) * betaale(nxi_scatt+1, ie, ile)
          do il = nxi_scatt, 1, -1
             gle(il, ie, ile) = (delta(il) - c1le(il, ie, ile) * gle(il+1, ie, ile)) * &
                  betaale(il, ie, ile)
          end do

          if (is_wfb .and. special_wfb_lorentz) then
             gle(ng2+2:je-1, ie, ile) = gle(ng2+1:je-2, ie, ile)
             gle(ng2+1, ie, ile) = gwfb
          end if

          ! next line ensures bounce condition is satisfied after lorentz collisions
          ! this is right thing to do, but it would mask any prior bug in trapping condition!
          if (cur_jend /= 0) gle(2 * cur_jend, ie, ile) = gle(cur_jend, ie, ile)

       end do
    end do
    !$OMP END PARALLEL DO
  end subroutine solfp_lorentz_le_layout