init_ediffuse Subroutine

  subroutine init_ediffuse (vnmult_target)
    use le_grids, only: negrid, forbid, ixi_to_il, speed => speed_maxw
    use gs2_layouts, only: le_lo, e_lo, il_idx
    use gs2_layouts, only: ig_idx, it_idx, ik_idx, is_idx
    use array_utils, only: zero_array
    implicit none
    real, intent (in), optional :: vnmult_target
    integer :: ie, is, ik, il, ig, it, ile, ixi, ielo
    real, dimension (:), allocatable :: aa, bb, cc, xe

    allocate (aa(negrid), bb(negrid), cc(negrid), xe(negrid))

    ! want to use x variables instead of e because we want conservative form
    ! for the x-integration
    xe = speed

    if (present(vnmult_target)) then
       vnmult(2) = max (vnmult_target, 1.0)
    end if

    if(use_le_layout) then

       if (.not.allocated(ec1le)) then
          allocate (ec1le   (nxi_lim, negrid, le_lo%llim_proc:le_lo%ulim_alloc))
          allocate (ebetaale(nxi_lim, negrid, le_lo%llim_proc:le_lo%ulim_alloc))
          allocate (eqle    (nxi_lim, negrid, le_lo%llim_proc:le_lo%ulim_alloc))
          vnmult(2) = max(1.0, vnmult(2))
       end if
       call zero_array(ec1le)
       call zero_array(ebetaale)
       call zero_array(eqle)

       !$OMP PARALLEL DO DEFAULT(none) &
       !$OMP PRIVATE(ile, ik, it, is, ig, ie, ixi, il, aa, bb, cc) &
       !$OMP SHARED(le_lo, nxi_lim, forbid, ec1le, ebetaale, negrid, eqle, ixi_to_il, xe) &
       !$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)
          is = is_idx(le_lo, ile)
          ig = ig_idx(le_lo, ile)
          do ixi = 1, nxi_lim
             il = ixi_to_il(ixi, ig)
             if (forbid(ig, il)) cycle
             call get_ediffuse_matrix (aa, bb, cc, ig, ik, it, il, is, xe)
             ec1le(ixi, :, ile) = cc
             ebetaale(ixi, 1, ile) = 1.0 / bb(1)
             do ie = 1,  negrid - 1
                eqle(ixi, ie+1, ile) = aa(ie+1) * ebetaale(ixi, ie, ile)
                ebetaale(ixi, ie+1, ile) = 1.0 / ( bb(ie+1) - eqle(ixi, ie+1, ile) &
                     * ec1le(ixi, ie, ile) )
             end do
          end do
       end do
       !$OMP END PARALLEL DO
    else

       if (.not.allocated(ec1)) then
          allocate (ec1    (negrid,e_lo%llim_proc:e_lo%ulim_alloc))
          allocate (ebetaa (negrid,e_lo%llim_proc:e_lo%ulim_alloc))
          allocate (eql    (negrid,e_lo%llim_proc:e_lo%ulim_alloc))
          vnmult(2) = max(1.0, vnmult(2))
       endif
       call zero_array(ec1)
       call zero_array(ebetaa)
       call zero_array(eql)

       !$OMP PARALLEL DO DEFAULT(none) &
       !$OMP PRIVATE(ielo, ik, it, is, ig, ie, ixi, il, aa, bb, cc) &
       !$OMP SHARED(e_lo, forbid, xe, ec1, ebetaa, negrid, eql) &
       !$OMP SCHEDULE(static)
       do ielo = e_lo%llim_proc, e_lo%ulim_proc
          il = il_idx(e_lo, ielo)
          ig = ig_idx(e_lo, ielo)
          if (forbid(ig, il)) cycle
          is = is_idx(e_lo, ielo)
          ik = ik_idx(e_lo, ielo)
          it = it_idx(e_lo, ielo)
          call get_ediffuse_matrix (aa, bb, cc, ig, ik, it, il, is, xe)
          ec1(:, ielo) = cc

          ! fill in the arrays for the tridiagonal
          ebetaa(1, ielo) = 1.0 / bb(1)
          do ie = 1, negrid - 1
             eql(ie+1, ielo) = aa(ie+1) * ebetaa(ie, ielo)
             ebetaa(ie+1, ielo) = 1.0 / (bb(ie+1) - eql(ie+1, ielo) * ec1(ie, ielo))
          end do
       end do
       !$OMP END PARALLEL DO
    end if

    deallocate(aa, bb, cc, xe)

  end subroutine init_ediffuse