ginit_nl4 Subroutine

  subroutine ginit_nl4
    use mp, only: proc0
    use species, only: spec
    use gs2_save, only: gs2_restore
    use theta_grid, only: ntgrid
    use kt_grids, only: naky, ntheta0
    use dist_fn_arrays, only: g, gnew
    use le_grids, only: forbid
    use fields_arrays, only: phi, apar, bpar
    use fields_arrays, only: phinew, aparnew, bparnew
    use gs2_layouts, only: g_lo, ik_idx, it_idx, is_idx, il_idx
    use file_utils, only: error_unit
    use run_parameters, only: has_phi, has_apar, has_bpar
    use ran, only: ranf
    use array_utils, only: copy, zero_array
    implicit none
    complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phiz
    integer :: iglo, istatus
    integer :: ig, ik, it, is, il, ierr
    
    call gs2_restore (g, scale, istatus, has_phi, has_apar, has_bpar)
    if (istatus /= 0) then
       ierr = error_unit()
       if (proc0) write(ierr,*) "Error reading file: ", trim(restart_file)
       call zero_array(g)
    end if

    do iglo = g_lo%llim_proc, g_lo%ulim_proc
       it = it_idx(g_lo,iglo)
       ik = ik_idx(g_lo,iglo)
!       if ((it == 2 .or. it == ntheta0) .and. ik == 1) cycle
!       if (ik == 1) cycle
       if (it == 1 .and. ik == 2) cycle

       g (:,1,iglo) = 0.
       g (:,2,iglo) = 0.
    end do

!    do ik = 1, naky
!       if (ik /= 2) then
!          phinew(:,:,ik) = 0.
!          aparnew(:,:,ik) = 0.
!          bparnew(:,:,ik) = 0.
!          phi(:,:,ik) = 0.
!          apar(:,:,ik) = 0.
!          bpar(:,:,ik) = 0.
!       else
!          phinew(:,2:ntheta0,ik) = 0.
!          aparnew(:,2:ntheta0,ik) = 0.
!          bparnew(:,2:ntheta0,ik) = 0.
!          phi(:,2:ntheta0,ik) = 0.
!          apar(:,2:ntheta0,ik) = 0.
!          bpar(:,2:ntheta0,ik) = 0.
!       end if
!    end do

    do ik = 1, naky
       do it=1,ntheta0
          if (it == 1 .and. ik == 2) cycle
          phinew(:,it,ik) = 0.
          aparnew(:,it,ik) = 0.
          bparnew(:,it,ik) = 0.
          phi(:,it,ik) = 0.
          apar(:,it,ik) = 0.
          bpar(:,it,ik) = 0.          
       end do
    end do
    
    do ik = 1, naky
       do it = 1, ntheta0
          do ig = -ntgrid, ntgrid
             phiz(ig,it,ik) = cmplx(ranf()-0.5,ranf()-0.5)
          end do
          if (chop_side) then
             if (left) then
                phiz(:-1,it,ik) = 0.0
             else
                phiz(1:,it,ik) = 0.0
             end if
          end if
       end do
    end do

! reality condition for k_theta = 0 component:
    do it = 1, ntheta0/2
       phiz(:,it+(ntheta0+1)/2,1) = conjg(phiz(:,(ntheta0+1)/2+1-it,1))
    enddo

    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)
       g(:,1,iglo) = g(:,1,iglo)-phiz(:,it,ik)*spec(is)%z*phiinit
       g(:,2,iglo) = g(:,2,iglo)-phiz(:,it,ik)*spec(is)%z*phiinit
       where (forbid(:,il)) g(:,1,iglo) = 0.0
       where (forbid(:,il)) g(:,2,iglo) = 0.0
    end do
    call copy(g, gnew)
  end subroutine ginit_nl4