conserve_lorentz_le_layout Subroutine

  subroutine conserve_lorentz_le_layout (gle)
    use gs2_layouts, only: ik_idx, it_idx, ie_idx, il_idx, is_idx, ig_idx, le_lo
    use le_grids, only: speed => speed_maxw, w, wxi, negrid
    use run_parameters, only: ieqzip
    use kt_grids, only: kwork_filter
    implicit none
    complex, dimension (:,:,le_lo%llim_proc:), intent (in out) :: gle
    complex :: v0y0, v1y1, v2y2
    real :: nud
    integer :: ig, ik, ie, is, it, ile, ixi

    if (drag) then
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       ! First get v0y0 and y1 = y0 - v0y0 * z0 / (1 + v0z0)

       ! v0 = vpa J0 f0, y0 = gle
       !$OMP PARALLEL DO DEFAULT(none) &
       !$OMP PRIVATE(ile, it, ik, ixi, ie, is, ig, v0y0) &
       !$OMP SHARED(le_lo, kwork_filter, ieqzip, negrid, nxi_lim, &
       !$OMP z0le, w, wxi, vpa_aj0_le, gle) &
       !$OMP SCHEDULE(static)
       do ile = le_lo%llim_proc, le_lo%ulim_proc
          it = it_idx(le_lo,ile)
          ik = ik_idx(le_lo,ile)
          if (kwork_filter(it, ik)) cycle
          if (ieqzip(it, ik)) cycle
          is = is_idx(le_lo, ile)
          ig = ig_idx(le_lo, ile)
          v0y0 = 0.0
          do ie = 1, negrid
             do ixi = 1, nxi_lim
                ! Should we use vpdiff here if conservative?
                v0y0 = v0y0 + vpa_aj0_le(ixi, ie, ile) * gle(ixi, ie, ile) * w(ie, is) * wxi(ixi, ig)
             end do
          end do
          gle(:, :, ile) = gle(:, :, ile) - z0le(:, :, ile) * v0y0
       end do
       !$OMP END PARALLEL DO
    end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Now get v1y1 and y2 = y1 - v1y1 * s1 / (1 + v1s1)

    ! v1 = nud vpa J0 f0, y1 = gle
    if (conservative) then
       !$OMP PARALLEL DO DEFAULT(none) &
       !$OMP PRIVATE(ile, is, it, ik, ig, ixi, ie, nud, v1y1) &
       !$OMP SHARED(le_lo, kwork_filter, ieqzip, negrid, nxi_lim, vpdiffle, &
       !$OMP speed, vnmult, vnew_D, aj0le, gle, s0le, w, wxi) &
       !$OMP SCHEDULE(static)
       do ile = le_lo%llim_proc, le_lo%ulim_proc
          ik = ik_idx(le_lo,ile)
          it = it_idx(le_lo,ile)
          if (kwork_filter(it, ik)) cycle
          if (ieqzip(it, ik)) cycle
          is = is_idx(le_lo,ile)
          ig = ig_idx(le_lo,ile)
          v1y1 = 0.0
          do ie = 1, negrid
             nud = speed(ie) * vnmult(1) * vnew_D(ik, ie, is) * w(ie, is)
             do ixi = 1, nxi_lim
                v1y1 = v1y1 + vpdiffle(ixi, ig) * nud * aj0le(ixi, ie, ile) * &
                     gle(ixi, ie, ile) * wxi(ixi, ig)
             end do
          end do
          gle(:, :, ile) = gle(:, :, ile) - s0le(:, :, ile) * v1y1
       end do
       !$OMP END PARALLEL DO
    else
       !$OMP PARALLEL DO DEFAULT(none) &
       !$OMP PRIVATE(ile, is, it, ik, ig, ixi, ie, nud, v1y1) &
       !$OMP SHARED(le_lo, kwork_filter, ieqzip, negrid, nxi_lim, vpa_aj0_le, vnmult, &
       !$OMP vnew_D, gle, w, wxi, s0le) &
       !$OMP SCHEDULE(static)
       do ile = le_lo%llim_proc, le_lo%ulim_proc
          ik = ik_idx(le_lo, ile)
          it = it_idx(le_lo, ile)
          if (kwork_filter(it, ik)) cycle
          if (ieqzip(it, ik)) cycle
          is = is_idx(le_lo, ile)
          ig = ig_idx(le_lo, ile)
          v1y1 = 0.0
          do ie = 1, negrid
             nud =  vnmult(1) * vnew_D(ik, ie, is) * w(ie, is)
             do ixi = 1, nxi_lim
                v1y1 = v1y1 + nud * vpa_aj0_le(ixi, ie, ile) * &
                     gle(ixi, ie, ile) * wxi(ixi, ig)
             end do
          end do
          gle(:, :, ile) = gle(:, :, ile) - s0le(:, :, ile) * v1y1
       end do
       !$OMP END PARALLEL DO
    end if

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Now get v2y2 and x = y2 - v2y2 * w2 / (1 + v2w2)

    !$OMP PARALLEL DO DEFAULT(none) &
    !$OMP PRIVATE(ile, it, ik, is, ie, ig, ixi, nud, v2y2) &
    !$OMP SHARED(le_lo, kwork_filter, ieqzip, negrid, nxi_lim, &
    !$OMP vnmult, vnew_D, vperp_aj1le, gle, w, wxi, w0le) &
    !$OMP SCHEDULE(static)
    do ile = le_lo%llim_proc, le_lo%ulim_proc
       it = it_idx(le_lo, ile)
       ik = ik_idx(le_lo, ile)
       if (kwork_filter(it, ik)) cycle
       if (ieqzip(it, ik)) cycle
       is = is_idx(le_lo, ile)
       ig = ig_idx(le_lo, ile)
       v2y2 = 0.0
       do ie = 1, negrid
          nud = vnmult(1) * vnew_D(ik, ie, is) * w(ie, is)
          do ixi = 1, nxi_lim
             ! aj1vp2 = 2 * J1(arg)/arg * vperp^2
             v2y2 = v2y2 + nud * vperp_aj1le(ixi, ie, ile) * gle(ixi, ie, ile) * wxi(ixi, ig)
          end do
       end do
       gle(:, :, ile) = gle(:, :, ile) - w0le(:, :, ile) * v2y2
    end do
    !$OMP END PARALLEL DO
  end subroutine conserve_lorentz_le_layout