ginit_single_parallel_mode Subroutine

  subroutine ginit_single_parallel_mode
    use species, only: spec, tracer_species
    use theta_grid, only: ntgrid, is_effectively_zero_shear, theta
    use kt_grids, only: naky, ntheta0, aky, reality
    use le_grids, only: forbid
    use mp, only: mp_abort
    use dist_fn_arrays, only: g, gnew
    use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx, is_idx
    use array_utils, only: copy, zero_array
    implicit none
    complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi, phit
    real :: a, b
    integer :: iglo
    integer :: ig, ik, it, il, is, nn

    call zero_array(phit)
    do it=2,ntheta0/2+1
       nn = it-1
! extra factor of 4 to reduce the high k wiggles for now
       phit (:, it, 1) = (-1)**nn*exp(-8.*(real(nn)/ntheta0)**2)
    end do

    if (is_effectively_zero_shear()) then
      if (ikpar_init+1 > ntgrid) then
        call mp_abort("Error: this value of k_parallel is too large. Increase ntheta or decrease ikpar_init.")
      end if
      kpar_init = ikpar_init
    end if

    do it = 1, ntheta0
       do ik = 1, naky
          do ig = -ntgrid, ntgrid
              !Set the field to cos(kpar_init*theta), where we remember that the gridpoints are not necessarily evenly spaced in the parallel direction, so we use theta(ig).
               !reality for ky=0 means we must use -kpar for kx < 0
              !if (naky == 1 .and. it > (ntheta0+1)/2) then
                !a = cos(-kpar_init * theta(ig))
                !b = sin(-kpar_init * theta(ig))
              !else
                a = cos(kpar_init * theta(ig))
                b = sin(kpar_init * theta(ig))
              !end if
              

!              a = ranf()-0.5
!              b = ranf()-0.5
!             phi(:,it,ik) = cmplx(a,b)
             phi(ig,it,ik) = cmplx(a,b)
           end do
          if (chop_side) then
             if (left) then
                phi(:-1,it,ik) = 0.0
             else
                phi(1:,it,ik) = 0.0
             end if
          end if
       end do
    end do

    if (naky > 1 .and. aky(1) == 0.0) then
       phi(:,:,1) = phi(:,:,1)*zf_init
    end if


    if (ikx_init  > 0) call single_initial_kx(phi)

    !<doc> reality condition for ky = 0 component: </doc>
    if (reality) then
       do it = 1, ntheta0/2
          phi(:,it+(ntheta0+1)/2,1) = conjg(phi(:,(ntheta0+1)/2+1-it,1))
          phit(:,it+(ntheta0+1)/2,1) = conjg(phit(:,(ntheta0+1)/2+1-it,1))
       enddo
    end if
       
    do iglo = g_lo%llim_proc, g_lo%ulim_proc
       ik = ik_idx(g_lo,iglo)
       it = it_idx(g_lo,iglo)
       il = il_idx(g_lo,iglo)
       is = is_idx(g_lo,iglo)
       if (spec(is)%type == tracer_species) then          
          g(:,1,iglo) =-phit(:,it,ik)*spec(is)%z*phiinit
       else
          g(:,1,iglo) = -phi(:,it,ik)*spec(is)%z*phiinit
       end if
       where (forbid(:,il)) g(:,1,iglo) = 0.0
       g(:,2,iglo) = g(:,1,iglo)
    end do
    call copy(g, gnew)
  end subroutine ginit_single_parallel_mode