conserve_diffuse_le_layout Subroutine

  subroutine conserve_diffuse_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: wxi, w, negrid
    use run_parameters, only: ieqzip
    use kt_grids, only: kwork_filter
    implicit none
    complex, dimension (:,:,le_lo%llim_proc:), intent (in out) :: gle
    real :: delnu, vnE
    integer :: ig, ik, ie, il, is, it, ile, ixi
    complex :: v0y0, v1y1, v2y2

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

    !$OMP PARALLEL DO DEFAULT(none) &
    !$OMP PRIVATE(ile, it, ik, is, ie, ig, ixi, vnE, v0y0) &
    !$OMP SHARED(le_lo, kwork_filter, negrid, nxi_lim, w, wxi, ieqzip, bz0le, vnmult, vnew_E, aj0le, gle) &
    !$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)
       v0y0 = 0.0
       do ie = 1, negrid
          vnE = vnmult(2) * vnew_E(ik, ie, is) * w(ie, is)
          do ixi = 1, nxi_lim
             v0y0 = v0y0 + vnE * aj0le(ixi, ie, ile) * gle(ixi, ie, ile) * wxi(ixi, ig)
          end do
       end do
       gle(:, :, ile) = gle(:, :, ile) - bz0le(:, :, ile) * v0y0
    end do
    !$OMP END PARALLEL DO

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

    !$OMP PARALLEL DO DEFAULT(none) &
    !$OMP PRIVATE(ile, it, ik, is, ig, delnu, ie, ixi, v1y1) &
    !$OMP SHARED(le_lo, kwork_filter, negrid, vnmult, delvnew, nxi_lim, &
    !$OMP bs0le, vpa_aj0_le, gle, w, wxi, ieqzip) &
    !$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
          delnu = vnmult(2) * delvnew(ik, ie, is) * w(ie, is)
          do ixi = 1, nxi_lim
             v1y1 = v1y1 + vpa_aj0_le(ixi,  ie,  ile) * delnu * gle(ixi, ie, ile) * wxi(ixi, ig)
          end do
       end do
       gle(:, :, ile) = gle(:, :, ile) - bs0le(:, :, ile) * v1y1
    end do
    !$OMP END PARALLEL DO

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

    !$OMP PARALLEL DO DEFAULT(none) &
    !$OMP PRIVATE(ile, it, ik, is, ie, ig, ixi, delnu, v2y2) &
    !$OMP SHARED(le_lo, kwork_filter, negrid, nxi_lim, ieqzip, &
    !$OMP bw0le, delvnew, vperp_aj1le, gle, vnmult, 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)
       v2y2 = 0.0
       do ie = 1, negrid
          delnu = vnmult(2) * delvnew(ik, ie, is) * w(ie, is)
          do ixi = 1, nxi_lim
             v2y2 = v2y2 + delnu * vperp_aj1le(ixi, ie, ile) * gle(ixi, ie, ile) * wxi(ixi, ig)
          end do
       end do
       gle(:, :, ile) = gle(:, :, ile) - bw0le(:, :, ile) * v2y2
    end do
    !$OMP END PARALLEL DO
  end subroutine conserve_diffuse_le_layout