solfp_lorentz_standard_layout Subroutine

  subroutine solfp_lorentz_standard_layout (g, gc, gh, diagnostics)
    use theta_grid, only: ntgrid
    use le_grids, only: jend, lambda_map, il_is_wfb, ng2, grid_has_trapped_particles
    use gs2_layouts, only: ig_idx, ik_idx, il_idx, is_idx, it_idx, g_lo, lz_lo
    use redistribute, only: gather, scatter
    use run_parameters, only: ieqzip
    use kt_grids, only: kwork_filter
    use array_utils, only: zero_array
    implicit none
    complex, dimension (-ntgrid:,:,g_lo%llim_proc:), intent (in out) :: g, gc, gh
    integer, optional, intent (in) :: diagnostics
    complex, dimension (:,:), allocatable :: glz, glzc
    complex, dimension (nxi_lim) :: delta
    complex :: fac, gwfb
    integer :: ig, ik, il, is, it, je, nxi_scatt, ilz, cur_jend
    logical :: is_wfb

    allocate (glz(nxi_lim,lz_lo%llim_proc:lz_lo%ulim_alloc))
    call zero_array(glz)

    call gather (lambda_map, g, glz)

    if (heating .and. present(diagnostics)) then
       allocate (glzc(nxi_lim, lz_lo%llim_proc:lz_lo%ulim_alloc))
       call zero_array(glzc)
       !$OMP PARALLEL DO DEFAULT(none) &
       !$OMP PRIVATE(ilz, ig, je, il, fac) &
       !$OMP SHARED(lz_lo, jend, ng2, glz, glzc, d1) &
       !$OMP SCHEDULE(static)
       do ilz = lz_lo%llim_proc, lz_lo%ulim_proc
          ig = ig_idx(lz_lo,ilz)

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

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

          do il = 1, je-1
             fac = glz(il+1,ilz)-glz(il,ilz)
             glzc(il,ilz) = conjg(fac)*fac*d1(il,ilz)  ! d1 accounts for hC(h) entropy
          end do
       end do
       !$OMP END PARALLEL DO

       call scatter (lambda_map, glzc, gc)

       if (hyper_colls) then
          !$OMP PARALLEL DO DEFAULT(none) &
          !$OMP PRIVATE(ilz, ig, je, il, fac) &
          !$OMP SHARED(lz_lo, jend, ng2, glz, glzc, h1) &
          !$OMP SCHEDULE(static)
          do ilz = lz_lo%llim_proc, lz_lo%ulim_proc
             ig = ig_idx(lz_lo,ilz)
             
             je = 2*jend(ig)          
             if (je == 0) then
                je = 2*ng2 
             end if
             
             do il = 1, je-1
                fac = glz(il+1,ilz)-glz(il,ilz)
                glzc(il,ilz) = conjg(fac)*fac*h1(il,ilz)  ! h1 accounts for hH(h) entropy
             end do
          end do
          !$OMP END PARALLEL DO

          call scatter (lambda_map, glzc, gh)
       end if
       if (allocated(glzc)) deallocate (glzc)
    end if

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

       !CMRDDGC, 10/2/2014:
       ! Fixes for wfb treatment below, use same je definition in ALL cases
       !   je  = #physical (unique) xi values + 1
       !       NB +1 above WITH TRAPPED is duplicate xi=vpar=0 point with isign=2
       !          +1 above WITHOUT TRAPPED is entirely unphysical extra grid point
       !  je-1 = #physical xi values removing unphysical/duplicate extra 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 (.not. grid_has_trapped_particles()) nxi_scatt = nxi_scatt - 1

       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.
          glz(cur_jend, ilz) = (glz(cur_jend, ilz) + glz(2 * cur_jend, ilz)) / 2.0
       end if

       ! zero unphysical/duplicate extra point. This shouldn't be needed
       if (grid_has_trapped_particles()) glz(je, ilz) = 0.0

       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 = glz(ng2+1, ilz)
          ! 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
          glz(ng2+1:je-2, ilz) = glz(ng2+2:je-1, ilz)
          ! Zero out the glz value not overwritten in the above. Shouldn't be needed
          glz(je - 1, ilz) = 0.0
          nxi_scatt = nxi_scatt - 1
       end if

       ! right and left sweeps for tridiagonal solve:
       delta(1) = glz(1, ilz)
       do il = 1, nxi_scatt
          delta(il+1) = glz(il+1, ilz) - ql(il+1, ilz) * delta(il)
       end do

       glz(nxi_scatt+1, ilz) = delta(nxi_scatt+1) * betaa(nxi_scatt+1, ilz)
       do il = nxi_scatt, 1, -1
          glz(il, ilz) = (delta(il) - c1(il, ilz) * glz(il+1, ilz)) * betaa(il, ilz)
       end do

       if (is_wfb .and. special_wfb_lorentz) then
          glz(ng2+2:je-1, ilz) = glz(ng2+1:je-2, ilz)
          glz(ng2+1, ilz) = gwfb
       end if

       ! update the duplicate vpar=0 point with the post pitch angle scattering value
       if (cur_jend /= 0) glz(2 * cur_jend, ilz) = glz(cur_jend, ilz)
    end do
    !$OMP END PARALLEL DO

    call scatter (lambda_map, glz, g)
    deallocate (glz)
  end subroutine solfp_lorentz_standard_layout