ginit_all_modes_equal Subroutine

  subroutine ginit_all_modes_equal
    use species, only: spec, tracer_species
    use theta_grid, only: ntgrid, theta, ntheta 
    use kt_grids, only: naky, ntheta0, aky, reality
    use le_grids, only: forbid
    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, ikpar

    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

    do it = 1, ntheta0
       do ik = 1, naky
          do ig = -ntgrid, ntgrid
              ! Set the field to cos(kpar*theta) for all kpar, where we remember that the gridpoints are not necessarily evenly spaced in the parallel direction, so we use theta(ig)</doc>
            a = 0.0 
            b = 0.0
            do ikpar = 0, ntheta - 1 
              a = a + cos(ikpar * theta(ig)) 
              ! we want to include the positive and negative wave numbers in
              ! equal measure, which of course means a real sine wave.
              b = 0.0 !b + cos(ikpar * theta(ig))
            end do

!              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_all_modes_equal