Restart but set all non-zonal components of the potential and the distribution function to 0. Noise can be added to these other components by setting phiinit > 0. The size of the zonal flows can be adjusted using zf_init.
Author EGH
subroutine ginit_restart_zonal_only
use gs2_save, only: gs2_restore
use mp, only: proc0
use file_utils, only: error_unit
use kt_grids, only: naky, ntheta0
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 fields_arrays, only: phinew
use run_parameters, only: has_phi, has_apar, has_bpar
use array_utils, only: copy, zero_array
implicit none
integer :: istatus, ierr
integer :: iglo
integer :: ik, it, il, is
! Load phi and g from the restart file
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
if (proc0) write (*,*) "Initialised g"
! Set all non-zonal components of phi to 0
do it = 1, ntheta0
do ik = 2, naky ! Starting at 2 is the crucial bit!
phinew(:,it,ik) = cmplx(0.0,0.0)
end do
end do
! Allow adjustment of the size of the zonal flows via the input parameter zf_init
phinew(:,:,1) = phinew(:,:,1)*zf_init
!Apply reality condition for k_theta = 0 component!
!if (reality) then
! do it = 1, ntheta0/2
! phinew(:,it+(ntheta0+1)/2,1) = conjg(phinew(:,(ntheta0+1)/2+1-it,1))
! enddo
! end if
!Set non-zonal components of g to zero using phi
if (proc0) write(*,*) "Zeroing g"
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 (ik > 1) then
g(:,1,iglo) = 0.0 !-phi(:,it,ik)*spec(is)%z*phiinit
where (forbid(:,il)) g(:,1,iglo) = 0.0
g(:,2,iglo) = g(:,1,iglo)
end if
end do
! If phiinit > 0, add some noise
if (phiinit > epsilon(0.0)) then
call ginit_noise
g = g + gnew
end if
call copy(g, gnew)
end subroutine ginit_restart_zonal_only