init_g.f90 Source File


Contents

Source Code


Source Code

!> This module contains the subroutines which set the initial value of the
!> fields and the distribution function.
module init_g
  use abstract_config, only: abstract_config_type, CONFIG_MAX_NAME_LEN
  use constants, only: run_name_size
  implicit none

  private

  public :: ginit
  public :: init_init_g, finish_init_g, wnml_init_g, check_init_g
  public :: width0
  public :: tstart
  public :: reset_init
  public :: init_vnmult
  public :: new_field_init, initial_condition_is_nonadiabatic_dfn
  public :: ginitopt_restart_many
  public :: ginitopt_restart_memory
  private :: single_initial_kx
  public :: restart_file
  public :: get_restart_dir

  public :: init_g_config_type
  public :: set_init_g_config
  public :: get_init_g_config
  
  ! knobs
  integer :: ginitopt_switch
  integer, parameter :: ginitopt_default = 1,  &
       ginitopt_xi = 3, ginitopt_xi2 = 4, ginitopt_rh = 5, ginitopt_zero = 6, &
       ginitopt_test3 = 7, ginitopt_convect = 8, ginitopt_restart_file = 9, &
       ginitopt_noise = 10, ginitopt_restart_many = 11, ginitopt_continue = 12, &
       ginitopt_nl = 13, ginitopt_kz0 = 14, ginitopt_restart_small = 15, &
       ginitopt_nl2 = 16, ginitopt_nl3 = 17, ginitopt_nl4 = 18, & 
       ginitopt_nl5 = 19, ginitopt_alf = 20, ginitopt_kpar = 21, &
       ginitopt_nl6 = 22, ginitopt_nl7 = 23, ginitopt_gs = 24, ginitopt_recon = 25, &
       ginitopt_nl3r = 26, ginitopt_smallflat = 27, ginitopt_harris = 28, &
       ginitopt_recon3 = 29, ginitopt_ot = 30, &
       ginitopt_zonal_only = 31, ginitopt_single_parallel_mode = 32, &
       ginitopt_all_modes_equal = 33, &
       ginitopt_default_odd=34, ginitopt_restart_memory=35, & !Note ginitopt_restart_memory isn't meant to be a user selected choice.
       ginitopt_restart_eig=36, ginitopt_no_zonal = 37, &
       ginitopt_stationary_mode = 38, ginitopt_cgyro = 39, &
       ginitopt_random_sine = 40, ginitopt_constant = 41

  real :: width0, dphiinit, phiinit, imfac, refac, zf_init, phifrac
  real :: den0, upar0, tpar0, tperp0
  real :: den1, upar1, tpar1, tperp1
  real :: den2, upar2, tpar2, tperp2
  real :: tstart, scale, apar0
  logical :: chop_side, clean_init, left, even, new_field_init
  logical :: initial_condition_is_nonadiabatic_dfn
  character (len = run_name_size) :: restart_file !WARNING: This is also defined in gs2_save. This should be addressed
  character (len=150) :: restart_dir
  integer, dimension(2) :: ikk, itt
  integer, dimension(3) :: ikkk,ittt
  complex, dimension (6) :: phiamp, aparamp
  integer :: restart_eig_id
  integer :: max_mode

  ! These are used for the function ginit_single_parallel_mode, and specify the
  !  kparallel to initialize. In the case of  zero magnetic shear, of course, the box 
  ! is periodic in the parallel direction, and so only harmonics of the box size 
  ! (i.e. ikpar_init) are allowed  EGH</doc>
  integer :: ikpar_init
  real :: kpar_init

  !>  This is used  in linear runs with flow shear  in order to track the
  !! evolution of a single Lagrangian mode.
  integer :: ikx_init

  ! RN> for recon3
  real :: phiinit0 ! amplitude of equilibrium
  real :: phiinit_rand ! amplitude of random perturbation
  real :: a0,b0 ! amplitude of equilibrium Apara or By 
  ! if b0 /= 0, u0 is rescaled to give this By amplitude by overriding a0
  ! if a0 /= 0, u0 is rescaled to give this Apara amplitude
  logical :: null_phi, null_bpar, null_apar ! nullify fields
  ! if use_{Phi,Bpar,Apar} = F, override corresponding flags
  integer :: adj_spec ! adjust input parameter of this spec
  ! if = 0, just warn if given conditions are not satisfied
  character (len=20) :: eq_type = ''
  ! equilibrium type = 'sinusoidal', 'porcelli', 'doubleharris'
  real :: prof_width=-.05
  ! width of porcelli, doubleharris profile
  ! if negative, this gives the ratio to the box size
  integer :: eq_mode_n=0, eq_mode_u=1
  ! mode number in x for sinusoidal equilibrium
  logical :: input_check_recon=.false.
  complex :: ukxy_pt(3)
  ! <RN

  !> Whether to use the constant_random number
  !! generator, which always gives the same
  !! random numbers. Used for testing.
  logical :: constant_random_flag
  
  logical, parameter :: debug = .false.
  logical :: initialized = .false.
  logical :: exist
  
  !> Used to represent the input configuration of init_g
  type, extends(abstract_config_type) :: init_g_config_type
     ! namelist : init_g_knobs
     ! indexed : false     

     !> Amplitude of equilibrium `Apara`: if `a0` is non-zero, \(u_0\) is
     !> rescaled to give this `Apara` amplitude
     real :: a0 = 0.0
     !> Used by [[init_g_knobs:ginit_option]] `"recon3"`: adjust input parameter
     !> of this species; if = 0, just warn if given conditions are not satisfied
     !>
     !> FIXME: What input parameter, what conditions?
     integer :: adj_spec = 0
     !> Used by [[init_g_knobs:ginit_option]] `"nl3r"`: set amplitude of `apar`
     !> if [[init_g_knobs:width0]] is negative
     real :: apar0 = 0.0
     !> Used in initialization for the Orszag-Tang 2D vortex problem
     !> ([[init_g_knobs:ginit_option]] `"ot"`). Controls size of `jpar` term.
     !> FIXME: Reference equations?
     complex, dimension(6) :: aparamp = cmplx(0.0,0.0)
     !> Amplitude of equilibrium `By`: if `b0` is non-zero, \(u_0\) is
     !> rescaled to give this `By` amplitude by overriding [[init_g_knobs:a0]]
     real :: b0 = 0.0
     !> Rarely needed. Forces asymmetry into initial condition. Warning: This
     !> does not behave as one may expect in flux tube simulations (see
     !> clean_init), this can be important as the default is to use chop_side
     logical :: chop_side = .true.
     !> Used with [[init_g_knobs:ginit_option]] `"noise"`. Ensures that in flux
     !> tube simulations the connected boundary points are initialised to the
     !> same value. Also allows for [[init_g_knobs:chop_side]] to behave
     !> correctly in flux tube simulations.
     logical :: clean_init = .true.
     !> Use the constant_random number generator, which always gives the same
     !> random numbers. Useful for testing.
     logical :: constant_random_flag = .false.
     !> Amplitude of zeroth density perturbation for [[init_g_knobs:ginit_option]]
     !> options:
     !>
     !> - `"kpar"` (see [[init_g:ginit_kpar]])
     !> - `"kz0"` (see [[init_g:ginit_kz0]])
     !> - `"nl3"` (see [[init_g:ginit_nl3]])
     !> - `"nl3r"` (see [[init_g:ginit_nl3r]])
     !> - `"nl7"` (see [[init_g:ginit_nl7]])
     !> - `"ot"` (see [[init_g:ginit_ot]])
     !> - `"recon"` (see [[init_g:ginit_recon]])
     real :: den0 = 1.0
     !> Amplitude of first density perturbation for [[init_g_knobs:ginit_option]]
     !> options:
     !>
     !> - `"gs"` (see [[init_g:ginit_gs]])
     !> - `"kpar"` (see [[init_g:ginit_kpar]])
     !> - `"nl3"` (see [[init_g:ginit_nl3]])
     !> - `"nl3r"` (see [[init_g:ginit_nl3r]])
     !> - `"nl7"` (see [[init_g:ginit_nl7]])
     !> - `"recon"` (see [[init_g:ginit_recon]])
     !> - `"recon3"` (see [[init_g:ginit_recon3]])
     real :: den1 = 0.0
     !> Amplitude of second density perturbation for [[init_g_knobs:ginit_option]]
     !> options:
     !>
     !> - `"kpar"` (see [[init_g:ginit_kpar]])
     !> - `"kz0"` (see [[init_g:ginit_kz0]])
     !> - `"nl3"` (see [[init_g:ginit_nl3]])
     !> - `"nl3r"` (see [[init_g:ginit_nl3r]])
     !> - `"nl7"` (see [[init_g:ginit_nl7]])
     !> - `"recon"` (see [[init_g:ginit_recon]])
     !> - `"recon3"` (see [[init_g:ginit_recon3]])
     real :: den2 = 0.0
     !> Amplitude of \(\phi\) for `"nl7"` (see [[init_g:ginit_nl7]])
     real :: dphiinit = 1.0
     !> Mode number in \(x\) for the sinusoidal equilbrium option of
     !> [[init_g_knobs:ginit_option]] `"recon3"` (see [[init_g:ginit_recon3]]),
     !> applies to density and temperatures
     integer :: eq_mode_n = 0
     !> Mode number in \(x\) for the sinusoidal equilbrium option of
     !> [[init_g_knobs:ginit_option]] `"recon3"` (see [[init_g:ginit_recon3]]),
     !> applies to parallel velocity
     integer :: eq_mode_u = 1
     !> Equilbrium choice for [[init_g_knobs:ginit_option]] `"recon3"` (see
     !> [[init_g:ginit_recon3]]). Choices are:
     !>
     !> - `"sinusoidal"`
     !> - `"porcelli"`
     !> - `"doubleharris"`
     character(len = 20) :: eq_type = 'sinusoidal'
     !> Sometimes initial conditions have definite parity; this picks the parity
     !> in those cases. Affects the following choices of
     !> [[init_g_knobs:ginit_option]]:
     !>
     !> - `"nl3"` (see [[init_g:ginit_nl3]])
     !> - `"nl3r"` (see [[init_g:ginit_nl3r]])
     !> - `"nl7"` (see [[init_g:ginit_nl7]])
     !> - `"recon"` (see [[init_g:ginit_recon]])
     !> - `"recon3"` (see [[init_g:ginit_recon3]])
     logical :: even = .true.
     !> Method for initialising `g`. There are many possible options, but the
     !> typical ones are:
     !>
     !>   - `"default"`: Gaussian in \(\theta\). See [[init_g_knobs:phiinit]],
     !>     [[init_g_knobs:width0]], [[init_g_knobs:chop_side]]
     !>   - `"default_odd"`: Guassian in \(\theta\) multiplied by \(\sin(\theta -
     !>     \theta_0)\). See [[init_g_knobs:phiinit]], [[init_g_knobs:width0]],
     !>     [[init_g_knobs:chop_side]]
     !>   - `"noise"`: Noise along field line. The recommended (but not default)
     !>     selection. See [[init_g_knobs:phiinit]], [[init_g_knobs:zf_init]]
     !>   - `"restart", "many"`: Restart, using `g` from restart files.
     !>   - `"single_parallel_mode"`: Initialise only with a single parallel mode
     !>     specified by either [[init_g_knobs:ikpar_init]] for periodic
     !>     boundary conditions or [[init_g_knobs:kpar_init]] for linked
     !>     boundary conditions. Intended for linear calculations.
     !>   - `"eig_restart"`: Uses the restart files written by the
     !>     eigensolver. Also see [[init_g_knobs:restart_eig_id]].
     !>   - `"random_sine"`: This option is notable as it attempts to make sure
     !>     that the initial condition satisfies the parallel boundary condition.
     !>
     !> All the options are detailed further in the [initialisation options
     !> documentation](../user_manual/initialisation.html)
     character(len = 20) :: ginit_option = "default"
     !> \(k_y\) indices of two modes to initialise for certain options
     integer, dimension(2) :: ikk = [1, 2]
     !> \(k_y\) indices of three modes to initialise for certain options
     integer, dimension(3) :: ikkk = [1, 2, 2]
     !> Parallel mode number to initialise for [[init_g_knobs:ginit_option]]
     !> `"single_parallel_mode"` with periodic boundary conditions
     integer :: ikpar_init = 0
     !> If greater than zero, initialise only the given \(k_x\) with
     !> [[init_g_knobs:ginit_option]] `"noise"`. This is useful for linear runs
     !> with flow shear, to track the evolution of a single Lagrangian mode.
     integer :: ikx_init = -1
     !> Amplitude of imaginary part of \(\phi\) for various
     !> [[init_g_knobs:ginit_option]] options
     real :: imfac = 0.0
     !> Do we want to include the explicit source terms and
     !> related timesteps when saving/restoring from the restart file.
     logical :: include_explicit_source_in_restart = .true.
     !> If true then the initial condition is used to set the non-adiabatic
     !> distribution function rather than the `g` evolved by GS2.
     logical :: initial_condition_is_nonadiabatic_dfn = .false.
     !> Print some debug information for [[init_g_knobs:ginit_option]] `"recon3"`
     logical :: input_check_recon = .false.
     !> \(k_x\) indices of two modes to initialise for certain options
     integer, dimension(2) :: itt = [1, 2]     
     !> \(k_x\) indices of three modes to initialise for certain options
     integer, dimension(3) :: ittt = [1, 2, 2]     
     !> Parallel mode number to initialise for [[init_g_knobs:ginit_option]]
     !> `"single_parallel_mode"` with linked boundary conditions
     real :: kpar_init = 0.0
     !> If [[init_g_knobs:chop_side]] is true, chop out left side in theta if
     !> true, right side if false. Applies to: FIXME: list of initial conditions
     logical :: left = .true.
     !> The maximum mode number to initialise for [[init_g_knobs:ginit_option]]
     !> `"random_sine"`. Actual maximum used will be lower of max_mode and `ntheta/8`
     integer :: max_mode = 128
     !> Change fields implicit initialisation somehow
     !>
     !> FIXME: Add documentation
     logical :: new_field_init = .true.
     !> "Nullify fields". Used by [[init_g_knobs:ginit_option]] `"recon3"`
     !>
     !> FIXME: Add documentation
     logical :: null_apar = .false.
     !> "Nullify fields". Used by [[init_g_knobs:ginit_option]] `"recon3"`
     !>
     !> FIXME: Add documentation
     logical :: null_bpar = .false.
     !> "Nullify fields". Used by [[init_g_knobs:ginit_option]] `"recon3"`
     !>
     !> FIXME: Add documentation
     logical :: null_phi = .false.
     !> Used in initialization for the Orszag-Tang 2D vortex problem.
     !>
     !> FIXME: expand
     complex, dimension(6) :: phiamp = cmplx(0.0,0.0)
     !> If doing multiple flux tube calculations, multiply
     !> [[init_g_knobs:phiinit]] by `(job * phifrac + 1)`. Allows for ensemble
     !> averaging of multiple flux tube calculations
     real :: phifrac = 0.1
     !> Average amplitude of initial perturbation of each Fourier mode. Used by
     !> most [[init_g_knobs:ginit_option]] options
     real :: phiinit = 1.0
     !> Amplitude of equilibrium. Used by [[init_g_knobs:ginit_option]] `"recon3"`
     real :: phiinit0 = 0.0
     !> Amplitude of random perturbation. Used by [[init_g_knobs:ginit_option]] `"recon3"`
     real :: phiinit_rand = 0.0
     !> Which processor is responsible for saving/reading potentials to/from restart
     !> files. If <0 then all procs write to restart files. If >= nproc then set to
     !> proc_to_save_fields = mod(proc_to_save_fields, nproc)
     integer :: proc_to_save_fields = 0
     !> Width of "porcelli", "doubleharris" profiles in
     !> [[init_g_knobs:ginit_option]] `"recon3"`. If negative, this gives the
     !> ratio to the box size
     real :: prof_width = -0.1
     !> Restart from many restart files. If false, use a single restart file:
     !> this requires GS2 to have been built with parallel netCDF.
     logical :: read_many = .false.
     !> Amplitude of real part of \(\phi\) for various
     !> [[init_g_knobs:ginit_option]] options
     real :: refac = 1.0
     !> Directory to read/write restart files to/from. Make sure this exists
     !> before running (see [[gs2_diagnostics_knobs:file_safety_check]])
     character(len = 150) :: restart_dir = "./"
     !> Used to select with eigensolver generated restart file to load. Sets
     !> `<id>` in `restart_file//eig_<id>//.<proc>` string used to set
     !> filename. Note that this is zero indexed.
     integer :: restart_eig_id = 0
     !> Basename of restart files
     !>
     !> @note gets a smart default
     character(len = run_name_size) :: restart_file = "restart_file_not_set.nc"
     !> Rescale amplitudes of fields for restart [[init_g_knobs:ginit_option]]
     !> options such as `"restart"`. Note that if `scale` is negative, the
     !> fields are scaled by `-scale / max(abs(phi))`!
     real :: scale = 1.0
     !> Amplitude of zeroth parallel temperature perturbation for
     !> [[init_g_knobs:ginit_option]] options:
     !>
     !> - `"kpar"` (see [[init_g:ginit_kpar]])
     !> - `"kz0"` (see [[init_g:ginit_kz0]])
     !> - `"nl3"` (see [[init_g:ginit_nl3]])
     !> - `"nl3r"` (see [[init_g:ginit_nl3r]])
     !> - `"nl7"` (see [[init_g:ginit_nl7]])
     !> - `"ot"` (see [[init_g:ginit_ot]])
     !> - `"recon"` (see [[init_g:ginit_recon]])
     real :: tpar0 = 0.0
     !> Amplitude of first parallel temperature perturbation for
     !> [[init_g_knobs:ginit_option]] options:
     !>
     !> - `"gs"` (see [[init_g:ginit_gs]])
     !> - `"kpar"` (see [[init_g:ginit_kpar]])
     !> - `"nl3"` (see [[init_g:ginit_nl3]])
     !> - `"nl3r"` (see [[init_g:ginit_nl3r]])
     !> - `"nl7"` (see [[init_g:ginit_nl7]])
     !> - `"recon"` (see [[init_g:ginit_recon]])
     !> - `"recon3"` (see [[init_g:ginit_recon3]])
     real :: tpar1 = 0.0
     !> Amplitude of second parallel temperature perturbation for
     !> [[init_g_knobs:ginit_option]] options:
     !>
     !> - `"kpar"` (see [[init_g:ginit_kpar]])
     !> - `"kz0"` (see [[init_g:ginit_kz0]])
     !> - `"nl3"` (see [[init_g:ginit_nl3]])
     !> - `"nl3r"` (see [[init_g:ginit_nl3r]])
     !> - `"nl7"` (see [[init_g:ginit_nl7]])
     !> - `"recon"` (see [[init_g:ginit_recon]])
     !> - `"recon3"` (see [[init_g:ginit_recon3]])
     real :: tpar2 = 0.0
     !> Amplitude of zeroth perpendicular temperature perturbation for
     !> [[init_g_knobs:ginit_option]] options:
     !>
     !> - `"kpar"` (see [[init_g:ginit_kpar]])
     !> - `"kz0"` (see [[init_g:ginit_kz0]])
     !> - `"nl3"` (see [[init_g:ginit_nl3]])
     !> - `"nl3r"` (see [[init_g:ginit_nl3r]])
     !> - `"nl7"` (see [[init_g:ginit_nl7]])
     !> - `"ot"` (see [[init_g:ginit_ot]])
     !> - `"recon"` (see [[init_g:ginit_recon]])
     real :: tperp0 = 0.0
     !> Amplitude of first perpendicular temperature perturbation for
     !> [[init_g_knobs:ginit_option]] options:
     !>
     !> - `"gs"` (see [[init_g:ginit_gs]])
     !> - `"kpar"` (see [[init_g:ginit_kpar]])
     !> - `"nl3"` (see [[init_g:ginit_nl3]])
     !> - `"nl3r"` (see [[init_g:ginit_nl3r]])
     !> - `"nl7"` (see [[init_g:ginit_nl7]])
     !> - `"recon"` (see [[init_g:ginit_recon]])
     !> - `"recon3"` (see [[init_g:ginit_recon3]])
     real :: tperp1 = 0.0
     !> Amplitude of second perpendicular temperature perturbation for
     !> [[init_g_knobs:ginit_option]] options:
     !>
     !> - `"kpar"` (see [[init_g:ginit_kpar]])
     !> - `"kz0"` (see [[init_g:ginit_kz0]])
     !> - `"nl3"` (see [[init_g:ginit_nl3]])
     !> - `"nl3r"` (see [[init_g:ginit_nl3r]])
     !> - `"nl7"` (see [[init_g:ginit_nl7]])
     !> - `"recon"` (see [[init_g:ginit_recon]])
     !> - `"recon3"` (see [[init_g:ginit_recon3]])
     real :: tperp2 = 0.0
     !> Force `t = tstart` at beginning of run
     real :: tstart = 0.0
     !> Perturbation amplitude of parallel velocity? Used by
     !> [[init_g_knobs:ginit_option]] `"recon3"`
     !>
     !> FIXME: Add documentation
     complex, dimension(3) :: ukxy_pt = cmplx(0.,0.)
     !> Amplitude of zeroth parallel velocity perturbation for
     !> [[init_g_knobs:ginit_option]] options:
     !>
     !> - `"kpar"` (see [[init_g:ginit_kpar]])
     !> - `"kz0"` (see [[init_g:ginit_kz0]])
     !> - `"nl3"` (see [[init_g:ginit_nl3]])
     !> - `"nl3r"` (see [[init_g:ginit_nl3r]])
     !> - `"nl7"` (see [[init_g:ginit_nl7]])
     !> - `"ot"` (see [[init_g:ginit_ot]])
     !> - `"recon"` (see [[init_g:ginit_recon]])
     real :: upar0 = 0.0
     !> Amplitude of first parallel velocity perturbation for
     !> [[init_g_knobs:ginit_option]] options:
     !>
     !> - `"gs"` (see [[init_g:ginit_gs]])
     !> - `"kpar"` (see [[init_g:ginit_kpar]])
     !> - `"nl3"` (see [[init_g:ginit_nl3]])
     !> - `"nl3r"` (see [[init_g:ginit_nl3r]])
     !> - `"nl7"` (see [[init_g:ginit_nl7]])
     !> - `"recon"` (see [[init_g:ginit_recon]])
     !> - `"recon3"` (see [[init_g:ginit_recon3]])
     real :: upar1 = 0.0
     !> Amplitude of second parallel velocity perturbation for
     !> [[init_g_knobs:ginit_option]] options:
     !>
     !> - `"kpar"` (see [[init_g:ginit_kpar]])
     !> - `"kz0"` (see [[init_g:ginit_kz0]])
     !> - `"nl3"` (see [[init_g:ginit_nl3]])
     !> - `"nl3r"` (see [[init_g:ginit_nl3r]])
     !> - `"nl7"` (see [[init_g:ginit_nl7]])
     !> - `"recon"` (see [[init_g:ginit_recon]])
     !> - `"recon3"` (see [[init_g:ginit_recon3]])
     real :: upar2 = 0.0
     !> Width of initial perturbation Gaussian envelope in \(\theta\)
     real :: width0 = -3.5
     !> Amplitude of initial zonal flow perturbations relative to other modes
     real :: zf_init = 1.0
   contains
     procedure, public :: read => read_init_g_config
     procedure, public :: write => write_init_g_config
     procedure, public :: reset => reset_init_g_config
     procedure, public :: broadcast => broadcast_init_g_config
     procedure, public, nopass :: get_default_name => get_default_name_init_g_config
     procedure, public, nopass :: get_default_requires_index => get_default_requires_index_init_g_config
  end type init_g_config_type
  
  type(init_g_config_type) :: init_g_config
  
contains

  !> FIXME : Add documentation
  subroutine wnml_init_g(unit)
    use run_parameters, only: k0
    implicit none
    integer, intent(in) :: unit
    if (.not.exist) return
    write (unit, *)
    write (unit, fmt="(' &',a)") "init_g_knobs"
    select case (ginitopt_switch)

    case (ginitopt_default)
       write (unit, fmt="(' ginit_option = ',a)") '"default"'
       write (unit, fmt="(' phiinit = ',e17.10)") phiinit
       write (unit, fmt="(' width0 = ',e17.10)") width0
       write (unit, fmt="(' chop_side = ',L1)") chop_side
       if (chop_side) write (unit, fmt="(' left = ',L1)") left

    case (ginitopt_noise)
       write (unit, fmt="(' ginit_option = ',a)") '"noise"'
       write (unit, fmt="(' phiinit = ',e17.10)") phiinit
       write (unit, fmt="(' zf_init = ',e17.10)") zf_init
       write (unit, fmt="(' chop_side = ',L1)") chop_side
       if (chop_side) write (unit, fmt="(' left = ',L1)") left
       write (unit, fmt="(' clean_init = ',L1)") clean_init

    case (ginitopt_xi)
       write (unit, fmt="(' ginit_option = ',a)") '"xi"'
       write (unit, fmt="(' width0 = ',e17.10)") width0

    case (ginitopt_xi2)
       write (unit, fmt="(' ginit_option = ',a)") '"xi2"'
       write (unit, fmt="(' width0 = ',e17.10)") width0

    case (ginitopt_zero)
       write (unit, fmt="(' ginit_option = ',a)") '"zero"'

    case (ginitopt_test3)
       write (unit, fmt="(' ginit_option = ',a)") '"test3"'

    case (ginitopt_convect)
       write (unit, fmt="(' ginit_option = ',a)") '"convect"'
       write (unit, fmt="(' k0 = ',e17.10)") k0

    case (ginitopt_rh)
       write (unit, fmt="(' ginit_option = ',a)") '"rh"'

    case (ginitopt_restart_many)
       write (unit, fmt="(' ginit_option = ',a)") '"many"'
       write (unit, fmt="(' restart_file = ',a)") '"'//trim(restart_file)//'"'
       write (unit, fmt="(' scale = ',e17.10)") scale

    case (ginitopt_restart_small)
       write (unit, fmt="(' ginit_option = ',a)") '"small"'
       write (unit, fmt="(' phiinit = ',e17.10)") phiinit
       write (unit, fmt="(' zf_init = ',e17.10)") zf_init
       write (unit, fmt="(' chop_side = ',L1)") chop_side
       if (chop_side) write (unit, fmt="(' left = ',L1)") left
       write (unit, fmt="(' restart_file = ',a)") '"'//trim(restart_file)//'"'
       write (unit, fmt="(' scale = ',e17.10)") scale

    case (ginitopt_restart_file)
       write (unit, fmt="(' ginit_option = ',a)") '"file"'
       write (unit, fmt="(' restart_file = ',a)") '"'//trim(restart_file)//'"'
       write (unit, fmt="(' scale = ',e17.10)") scale

    case (ginitopt_continue)
       write (unit, fmt="(' ginit_option = ',a)") '"cont"'

    case (ginitopt_kz0)
       write (unit, fmt="(' ginit_option = ',a)") '"kz0"'
       write (unit, fmt="(' phiinit = ',e17.10)") phiinit
       write (unit, fmt="(' chop_side = ',L1)") chop_side
       if (chop_side) write (unit, fmt="(' left = ',L1)") left

    case (ginitopt_nl)
       write (unit, fmt="(' ginit_option = ',a)") '"nl"'
       write (unit, fmt="(' phiinit = ',e17.10)") phiinit
       write (unit, fmt="(' ikk(1) = ',i3,' itt(1) = ',i3)") ikk(1),itt(1)
       write (unit, fmt="(' ikk(2) = ',i3,' itt(2) = ',i3)") ikk(2), itt(2)
       write (unit, fmt="(' chop_side = ',L1)") chop_side
       if (chop_side) write (unit, fmt="(' left = ',L1)") left

    case (ginitopt_nl2)
       write (unit, fmt="(' ginit_option = ',a)") '"nl2"'
       write (unit, fmt="(' phiinit = ',e17.10)") phiinit
       write (unit, fmt="(' ikk(1) = ',i3,' itt(1) = ',i3)") ikk(1),itt(1)
       write (unit, fmt="(' ikk(2) = ',i3,' itt(2) = ',i3)") ikk(2), itt(2)
       write (unit, fmt="(' chop_side = ',L1)") chop_side
       if (chop_side) write (unit, fmt="(' left = ',L1)") left

    case (ginitopt_nl3)
       write (unit, fmt="(' ginit_option = ',a)") '"nl3"'
       write (unit, fmt="(' phiinit = ',e17.10)") phiinit
       write (unit, fmt="(' width0 = ',e17.10)") width0
       write (unit, fmt="(' refac = ',e17.10)") refac
       write (unit, fmt="(' imfac = ',e17.10)") imfac
       write (unit, fmt="(' ikk(1) = ',i3,' itt(1) = ',i3)") ikk(1),itt(1)
       write (unit, fmt="(' ikk(2) = ',i3,' itt(2) = ',i3)") ikk(2), itt(2)
       write (unit, fmt="(' chop_side = ',L1)") chop_side
       if (chop_side) write (unit, fmt="(' left = ',L1)") left
       write (unit, fmt="(' den0 = ',e17.10)") den0
       write (unit, fmt="(' den1 = ',e17.10)") den1
       write (unit, fmt="(' den2 = ',e17.10)") den2
       write (unit, fmt="(' upar0 = ',e17.10)") upar0
       write (unit, fmt="(' upar1 = ',e17.10)") upar1
       write (unit, fmt="(' upar2 = ',e17.10)") upar2
       write (unit, fmt="(' tpar0 = ',e17.10)") tpar0
       write (unit, fmt="(' tpar1 = ',e17.10)") tpar1
       write (unit, fmt="(' tperp0 = ',e17.10)") tperp0
       write (unit, fmt="(' tperp1 = ',e17.10)") tperp1
       write (unit, fmt="(' tperp2 = ',e17.10)") tperp2

    case (ginitopt_nl4)
       write (unit, fmt="(' ginit_option = ',a)") '"nl4"'
       write (unit, fmt="(' phiinit = ',e17.10)") phiinit
       write (unit, fmt="(' restart_file = ',a)") '"'//trim(restart_file)//'"'
       write (unit, fmt="(' scale = ',e17.10)") scale
       write (unit, fmt="(' ikk(1) = ',i3,' itt(1) = ',i3)") ikk(1),itt(1)
       write (unit, fmt="(' ikk(2) = ',i3,' itt(2) = ',i3)") ikk(2), itt(2)
       write (unit, fmt="(' chop_side = ',L1)") chop_side
       if (chop_side) write (unit, fmt="(' left = ',L1)") left

    case (ginitopt_nl5)
       write (unit, fmt="(' ginit_option = ',a)") '"nl5"'
       write (unit, fmt="(' phiinit = ',e17.10)") phiinit
       write (unit, fmt="(' restart_file = ',a)") '"'//trim(restart_file)//'"'
       write (unit, fmt="(' scale = ',e17.10)") scale
       write (unit, fmt="(' chop_side = ',L1)") chop_side
       if (chop_side) write (unit, fmt="(' left = ',L1)") left

    case (ginitopt_nl6)
       write (unit, fmt="(' ginit_option = ',a)") '"nl6"'
       write (unit, fmt="(' phiinit = ',e17.10)") phiinit
       write (unit, fmt="(' restart_file = ',a)") '"'//trim(restart_file)//'"'
       write (unit, fmt="(' scale = ',e17.10)") scale

    case (ginitopt_alf)
       write (unit, fmt="(' ginit_option = ',a)") '"alf"'
       write (unit, fmt="(' phiinit = ',e17.10)") phiinit

    case (ginitopt_gs)
       write (unit, fmt="(' ginit_option = ',a)") '"gs"'
       write (unit, fmt="(' phiinit = ',e17.10)") phiinit
       write (unit, fmt="(' refac = ',e17.10)") refac
       write (unit, fmt="(' imfac = ',e17.10)") imfac
       write (unit, fmt="(' den1 = ',e17.10)") den1
       write (unit, fmt="(' upar1 = ',e17.10)") upar1
       write (unit, fmt="(' tpar1 = ',e17.10)") tpar1
       write (unit, fmt="(' tperp1 = ',e17.10)") tperp1


    case (ginitopt_kpar)
       write (unit, fmt="(' ginit_option = ',a)") '"kpar"'
       write (unit, fmt="(' phiinit = ',e17.10)") phiinit
       write (unit, fmt="(' width0 = ',e17.10)") width0
       write (unit, fmt="(' refac = ',e17.10)") refac
       write (unit, fmt="(' imfac = ',e17.10)") imfac
       write (unit, fmt="(' den0 = ',e17.10)") den0
       write (unit, fmt="(' den1 = ',e17.10)") den1
       write (unit, fmt="(' den2 = ',e17.10)") den2
       write (unit, fmt="(' upar0 = ',e17.10)") upar0
       write (unit, fmt="(' upar1 = ',e17.10)") upar1
       write (unit, fmt="(' upar2 = ',e17.10)") upar2
       write (unit, fmt="(' tpar0 = ',e17.10)") tpar0
       write (unit, fmt="(' tpar1 = ',e17.10)") tpar1
       write (unit, fmt="(' tperp0 = ',e17.10)") tperp0
       write (unit, fmt="(' tperp1 = ',e17.10)") tperp1
       write (unit, fmt="(' tperp2 = ',e17.10)") tperp2

    case (ginitopt_cgyro)
       write (unit, fmt="(' ginit_option = ',a)") '"cgyro"'
       write (unit, fmt="(' phiinit = ',e17.10)") phiinit
       write (unit, fmt="(' zf_init = ',e17.10)") zf_init
       write (unit, fmt="(' chop_side = ',L1)") chop_side
       if (chop_side) write (unit, fmt="(' left = ',L1)") left
       write (unit, fmt="(' clean_init = ',L1)") clean_init

    case (ginitopt_random_sine)
       write (unit, fmt="(' ginit_option = ',a)") '"random_sine"'
       write (unit, fmt="(' phiinit = ',e17.10)") phiinit
       write (unit, fmt="(' zf_init = ',e17.10)") zf_init
       write (unit, fmt="(' max_mode = ',I0)") max_mode

    case (ginitopt_constant)
       write (unit, fmt="(' ginit_option = ',a)") '"constant"'
       write (unit, fmt="(' phiinit = ',e17.10)") phiinit
       write (unit, fmt="(' zf_init = ',e17.10)") zf_init

    end select
    write (unit, fmt="(' /')")
  end subroutine wnml_init_g

  !> FIXME : Add documentation  
  subroutine check_init_g(report_unit)
    use run_parameters, only : delt_option_switch, delt_option_auto
    use species, only : spec, has_electron_species
    implicit none
    integer, intent(in) :: report_unit
    select case (ginitopt_switch)
    case (ginitopt_default)
       write (report_unit, fmt="('Initial conditions:')")
       write (report_unit, fmt="('  Amplitude:        ',f10.4)") phiinit
       write (report_unit, fmt="('  Width in theta:   ',f10.4)") width0
       if (chop_side) then
          write (report_unit, fmt="('  Parity:   none')") 
       else
          write (report_unit, fmt="('  Parity:   even')") 
       end if

    case (ginitopt_kz0)
       write (report_unit, fmt="('Initial conditions:')")
       write (report_unit, fmt="('  Amplitude:        ',f10.4)") phiinit
       write (report_unit, fmt="('  Constant along field line',f10.4)") width0
       if (chop_side) then
          write (report_unit, *) 
          write (report_unit, fmt="('################# WARNING #######################')")
          write (report_unit, fmt="('  Parity:   none')") 
          write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") 
          write (report_unit, fmt="('Remedy: set chop_side = .false. in init_g_knobs.')") 
          write (report_unit, fmt="('################# WARNING #######################')")
          write (report_unit, *) 
       end if

    case (ginitopt_noise)
       write (report_unit, fmt="('Initial conditions:')")
       write (report_unit, fmt="('  Amplitude:        ',f10.4)") phiinit
       write (report_unit, fmt="('  Noise along field line.')") 
       if (zf_init /= 1.) then
          write (report_unit, fmt="('  Zonal flows adjusted by factor of zf_init = ',f10.4)") zf_init
       end if

    case (ginitopt_kpar)
       write (report_unit, fmt="('Initial conditions:')")
       write (report_unit, fmt="('  Amplitude:             ',f10.4)") phiinit
       write (report_unit, fmt="('  Real part multiplier:  ',f10.4)") refac
       write (report_unit, fmt="('  Imag part multiplier:  ',f10.4)") imfac
       if (width0 > 0.) then
          write (report_unit, fmt="('  Gaussian envelope in theta with width:  ',f10.4)") width0
       end if
       if (chop_side) then
          write (report_unit, *) 
          write (report_unit, fmt="('################# WARNING #######################')")
          write (report_unit, fmt="('  Parity:   none')") 
          write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") 
          write (report_unit, fmt="('Remedy: set chop_side = .false. in init_g_knobs.')") 
          write (report_unit, fmt="('################# WARNING #######################')")
          write (report_unit, *) 
       end if
       if (den0 > epsilon(0.0) .or. den1 > epsilon(0.0) .or. den2 > epsilon(0.0)) then
          write (report_unit, *) 
          write (report_unit, fmt="('Initial density perturbation of the form:')")
          write (report_unit, fmt="('den0   + den1 * cos(theta) + den2 * cos(2.*theta)')")
          write (report_unit, *) 
          write (report_unit, fmt="('with den0 =',f7.4,' den1 = ',f7.4,' den2 = ',f7.4)") den0, den1, den2
       end if
       if (upar0 > epsilon(0.0) .or. upar1 > epsilon(0.0) .or. upar2 > epsilon(0.0)) then
          write (report_unit, *) 
          write (report_unit, fmt="('Initial parallel velocity perturbation of the form:')")
          write (report_unit, fmt="('upar0   + upar1 * sin(theta) + upar2 * sin(2.*theta)')")
          write (report_unit, fmt="('90 degrees out of phase with other perturbations.')")
          write (report_unit, *) 
          write (report_unit, fmt="('with upar0 =',f7.4,' upar1 = ',f7.4,' upar2 = ',f7.4)") upar0, upar1, upar2
       end if
       if (tpar0 > epsilon(0.0) .or. tpar1 > epsilon(0.0) .or. tpar2 > epsilon(0.0)) then
          write (report_unit, *) 
          write (report_unit, fmt="('Initial Tpar perturbation of the form:')")
          write (report_unit, fmt="('tpar0   + tpar1 * cos(theta) + tpar2 * cos(2.*theta)')")
          write (report_unit, *) 
          write (report_unit, fmt="('with tpar0 =',f7.4,' tpar1 = ',f7.4,' tpar2 = ',f7.4)") tpar0, tpar1, tpar2
       end if
       if (tperp0 > epsilon(0.0) .or. tperp1 > epsilon(0.0) .or. tperp2 > epsilon(0.0)) then
          write (report_unit, *) 
          write (report_unit, fmt="('Initial Tperp perturbation of the form:')")
          write (report_unit, fmt="('tperp0   + tperp1 * cos(theta) + tperp2 * cos(2.*theta)')")
          write (report_unit, *) 
          write (report_unit, fmt="('with tperp0 =',f7.4,' tperp1 = ',f7.4,' tperp2 = ',f7.4)") tperp0, tperp1, tperp2
       end if
       if (has_electron_species(spec)) then
          write (report_unit, *) 
          write (report_unit, fmt="('Field line average of g_electron subtracted off.')")
       end if

    case (ginitopt_gs)
       write (report_unit, fmt="('Initial conditions:')")
       write (report_unit, fmt="('  Randomly phased kpar=1 sines and cosines')") 
       write (report_unit, fmt="('  in density, upar, tpar, or tperp.')") 
       write (report_unit, fmt="('  Real part amplitude:  ',f10.4)") refac*phiinit
       write (report_unit, fmt="('  Imag part amplitude:  ',f10.4)") imfac*phiinit
       if (abs( den1)  > epsilon(0.0)) write (report_unit, fmt="('  Density amplitude:  ',f10.4)") den1
       if (abs( upar1) > epsilon(0.0)) write (report_unit, fmt="('  Upar amplitude:  ',f10.4)") upar1
       if (abs( tpar1) > epsilon(0.0)) write (report_unit, fmt="('  Tpar amplitude:  ',f10.4)") tpar1
       if (abs(tperp1) > epsilon(0.0)) write (report_unit, fmt="('  Tperp amplitude:  ',f10.4)") tperp1

    case (ginitopt_nl)
       write (report_unit, fmt="('Initial conditions:')")
       write (report_unit, fmt="('At most two k_perps excited, with amplitude = ',f10.4)") phiinit
       write (report_unit, fmt="(' First k_perp has ik = ',i3,' it = ',i3)") ikk(1), itt(1)
       write (report_unit, fmt="('Second k_perp has ik = ',i3,' it = ',i3)") ikk(2), itt(2)
       if (chop_side) then
          write (report_unit, fmt="('  Parity:   none')") 
       else
          write (report_unit, fmt="('  Parity:   even')") 
       end if
       write (report_unit, fmt="('Reality condition is enforced.')")
       
    case (ginitopt_nl2)
       write (report_unit, fmt="('Initial conditions:')")
       write (report_unit, fmt="('At most two k_perps excited, with amplitude = ',f10.4)") phiinit
       write (report_unit, fmt="(' First k_perp has ik = ',i3,' it = ',i3)") ikk(1), itt(1)
       write (report_unit, fmt="('Second k_perp has ik = ',i3,' it = ',i3)") ikk(2), itt(2)
       if (chop_side) then
          write (report_unit, fmt="('  Parity:   none')") 
       else
          write (report_unit, fmt="('  Parity:   even')") 
       end if
       write (report_unit, fmt="('Reality condition is enforced.')")
       write (report_unit, fmt="('g perturbation proportional to (1+v_parallel)*sin(theta)')")

    case (ginitopt_nl3)
       write (report_unit, fmt="('Initial conditions:')")
       write (report_unit, fmt="('At most two k_perps excited, with amplitude = ',f10.4)") phiinit
       write (report_unit, fmt="(' First k_perp has ik = ',i3,' it = ',i3)") ikk(1), itt(1)
       write (report_unit, fmt="('Second k_perp has ik = ',i3,' it = ',i3)") ikk(2), itt(2)
       write (report_unit, fmt="('  Real part multiplied by:  ',f10.4)") refac
       write (report_unit, fmt="('  Imag part multiplied by:  ',f10.4)") imfac
       if (width0 > 0.) then
          write (report_unit, fmt="('  Gaussian envelope in theta with width:  ',f10.4)") width0
       end if
       if (chop_side) then
          write (report_unit, fmt="('  Parity:   none')") 
       else
          write (report_unit, fmt="('  Parity:   even')") 
       end if
       write (report_unit, fmt="('Reality condition is enforced.')")
       if (den0 > epsilon(0.0) .or. den1 > epsilon(0.0) .or. den2 > epsilon(0.0)) then
          write (report_unit, *) 
          write (report_unit, fmt="('Initial density perturbation of the form:')")
          write (report_unit, fmt="('den0   + den1 * cos(theta) + den2 * cos(2.*theta)')")
          write (report_unit, *) 
          write (report_unit, fmt="('with den0 =',f7.4,' den1 = ',f7.4,' den2 = ',f7.4)") den0, den1, den2
       end if
       if (upar0 > epsilon(0.0) .or. upar1 > epsilon(0.0) .or. upar2 > epsilon(0.0)) then
          write (report_unit, *) 
          write (report_unit, fmt="('Initial parallel velocity perturbation of the form:')")
          write (report_unit, fmt="('upar0   + upar1 * cos(theta) + upar2 * cos(2.*theta)')")
          write (report_unit, fmt="('90 degrees out of phase with other perturbations.')")
          write (report_unit, *) 
          write (report_unit, fmt="('with upar0 =',f7.4,' upar1 = ',f7.4,' upar2 = ',f7.4)") upar0, upar1, upar2
       end if
       if (tpar0 > epsilon(0.0) .or. tpar1 > epsilon(0.0) .or. tpar2 > epsilon(0.0)) then
          write (report_unit, *) 
          write (report_unit, fmt="('Initial Tpar perturbation of the form:')")
          write (report_unit, fmt="('tpar0   + tpar1 * cos(theta) + tpar2 * cos(2.*theta)')")
          write (report_unit, *) 
          write (report_unit, fmt="('with tpar0 =',f7.4,' tpar1 = ',f7.4,' tpar2 = ',f7.4)") tpar0, tpar1, tpar2
       end if
       if (tperp0 > epsilon(0.0) .or. tperp1 > epsilon(0.0) .or. tperp2 > epsilon(0.0)) then
          write (report_unit, *) 
          write (report_unit, fmt="('Initial Tperp perturbation of the form:')")
          write (report_unit, fmt="('tperp0   + tperp1 * cos(theta) + tperp2 * cos(2.*theta)')")
          write (report_unit, *) 
          write (report_unit, fmt="('with tperp0 =',f7.4,' tperp1 = ',f7.4,' tperp2 = ',f7.4)") tperp0, tperp1, tperp2
       end if
       
    case (ginitopt_nl4)
       write (report_unit, fmt="('Initial conditions:')")
       write (report_unit, *) 
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, fmt="('Under development for study of secondary instabilities.')")
       write (report_unit, fmt="('Scale factor:   ',f10.4)") scale
       write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") 
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, *) 

    case (ginitopt_nl5)
       write (report_unit, fmt="('Initial conditions:')")
       write (report_unit, *) 
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, fmt="('Under development for study of secondary instabilities.')")
       write (report_unit, fmt="('Scale factor:   ',f10.4)") scale
       write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") 
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, *) 

    case (ginitopt_nl6)
       write (report_unit, fmt="('Initial conditions:')")
       write (report_unit, *) 
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, fmt="('Change amplitude of a particular mode.')")
       write (report_unit, fmt="('Scale factor:   ',f10.4)") scale
       write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") 
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, *) 

    case (ginitopt_xi)
       write (report_unit, fmt="('Initial conditions:')")
       write (report_unit, *) 
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, fmt="('Perturbation proportional to pitch angle.')")
       write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") 
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, *) 

    case (ginitopt_xi2)
       write (report_unit, fmt="('Initial conditions:')")
       write (report_unit, *) 
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, fmt="('Perturbation proportional to function of pitch angle.')")
       write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") 
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, *) 

    case (ginitopt_rh)
       write (report_unit, fmt="('Initial conditions:')")
       write (report_unit, *) 
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, fmt="('Maxwellian perturbation in ik=1 mode.')")
       write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") 
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, *) 

    case (ginitopt_alf)
       write (report_unit, fmt="('Initial conditions:')")
       write (report_unit, *) 
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, fmt="('Ion dist fn proportional to v_parallel * sin(theta).')")
       write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") 
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, *) 

    case (ginitopt_zero)
       write (report_unit, fmt="('Initial conditions:')")
       write (report_unit, *) 
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, fmt="('Distribution function = 0.')")
       write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") 
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, *) 

    case (ginitopt_test3)
       write (report_unit, fmt="('Initial conditions:')")
       write (report_unit, *) 
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") 
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, *) 

    case (ginitopt_convect)
       write (report_unit, fmt="('Initial conditions:')")
       write (report_unit, *) 
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") 
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, *) 

    case (ginitopt_restart_file)
       write (report_unit, fmt="('Initial conditions:')")
       write (report_unit, *) 
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, fmt="('Restart from a single NetCDF restart file.')") 
       write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") 
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, *) 

    case (ginitopt_restart_many)
       write (report_unit, fmt="('Initial conditions:')")
       write (report_unit, fmt="('Each PE restarts from its own NetCDF restart file.')") 

    case (ginitopt_restart_small)
       write (report_unit, fmt="('Initial conditions:')")
       write (report_unit, fmt="('Each PE restarts from its own NetCDF restart file.')") 
       write (report_unit, fmt="('with amplitudes scaled by factor of scale = ',f10.4)") scale
       write (report_unit, fmt="('Noise added with amplitude = ',f10.4)") phiinit

    case (ginitopt_continue)
       write (report_unit, fmt="('Initial conditions:')")
       write (report_unit, *) 
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") 
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, *) 

    case (ginitopt_cgyro)
       write (report_unit, fmt="('Initial conditions:')")
       write (report_unit, fmt="('  Amplitude:        ',f10.4)") phiinit
       write (report_unit, fmt="('Constant amplitude for ky!=0')")
       write (report_unit, fmt="('Zero amplitude for ky=0')")

    case (ginitopt_random_sine)
       write (report_unit, fmt="('Initial conditions:')")
       write (report_unit, fmt="('  Amplitude:        ',f10.4)") phiinit
       write (report_unit, fmt="('  Zonal amplitude:  ',f10.4)") zf_init
       write (report_unit, fmt="('  Max mode:         ',I0)") max_mode
       write (report_unit, fmt="(' Random sum of sine waves up to max_mode mode number')")

    case (ginitopt_constant)
       write (report_unit, fmt="('Initial conditions:')")
       write (report_unit, fmt="('  Amplitude:        ',f10.4)") phiinit
       write (report_unit, fmt="('  Zonal amplitude:  ',f10.4)") zf_init
       write (report_unit, fmt="(' Constant value')")
    end select

    if (ginitopt_switch == ginitopt_restart_many) then
       if (delt_option_switch == delt_option_auto) then
          write (report_unit, *) 
          write (report_unit, fmt="('This run is a continuation of a previous run.')") 
          write (report_unit, fmt="('The time step at the beginning of this run')") 
          write (report_unit, fmt="('will be taken from the end of the previous run.')") 
       else
          write (report_unit, *) 
          write (report_unit, fmt="('################# WARNING #######################')")
          write (report_unit, fmt="('This run is a continuation of a previous run.')") 
          write (report_unit, fmt="('The time step is being set by hand.')") 
          write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") 
          write (report_unit, fmt="('You probably want to set delt_option to be check_restart in the knobs namelist.')") 
          write (report_unit, fmt="('################# WARNING #######################')")
          write (report_unit, *) 
       end if
    end if

    if (delt_option_switch == delt_option_auto) then
       if (ginitopt_switch /= ginitopt_restart_many) then
          write (report_unit, *) 
          write (report_unit, fmt="('################# WARNING #######################')")
          write (report_unit, fmt="('This is not a normal continuation run.')") 
          write (report_unit, fmt="('You probably want to set delt_option to be default in the knobs namelist.')") 
          write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") 
          write (report_unit, fmt="('################# WARNING #######################')")
          write (report_unit, *) 
       end if
    end if
  end subroutine check_init_g

  !> FIXME : Add documentation  
  subroutine init_init_g(init_g_config_in)
    use gs2_save, only: set_restart_file, read_many, proc_to_save_fields
    use gs2_layouts, only: init_gs2_layouts
    use mp, only: proc0, broadcast, job, nproc
    implicit none
    type(init_g_config_type), intent(in), optional :: init_g_config_in    
    integer :: ind_slash

    if (initialized) return
    initialized = .true.
    call init_gs2_layouts

    call read_parameters(init_g_config_in)

    ! prepend restart_dir to restart_file
    ! append trailing slash if not exists
    if(restart_dir(len_trim(restart_dir):) /= "/") &
         restart_dir=trim(restart_dir)//"/"
!Determine if restart file contains "/" if so split on this point to give DIR//FILE
    !so restart files are created in DIR//restart_dir//FILE
    ind_slash=index(restart_file,"/",.True.)
    if (ind_slash.EQ.0) then !No slash present
       restart_file=trim(restart_dir)//trim(restart_file)
    else !Slash present
       restart_file=trim(restart_file(1:ind_slash))//trim(restart_dir)//trim(restart_file(ind_slash+1:))
    endif

    ! MAB - allows for ensemble averaging of multiple flux tube calculations
    ! job=0 if not doing multiple flux tube calculations, so phiinit unaffected
    !<DD>Only let proc0 do this and then broadcast as other procs don't know phiinit,phifrac etc.
    if(proc0) phiinit = phiinit * (job*phifrac+1.0)

    call broadcast (ginitopt_switch)
    call broadcast (width0)
    call broadcast (refac)
    call broadcast (imfac)
    call broadcast (den0)
    call broadcast (upar0)
    call broadcast (tpar0)
    call broadcast (tperp0)
    call broadcast (den1)
    call broadcast (upar1)
    call broadcast (tpar1)
    call broadcast (tperp1)
    call broadcast (den2)
    call broadcast (upar2)
    call broadcast (tpar2)
    call broadcast (tperp2)
    call broadcast (phiinit)
    call broadcast (phifrac)
    call broadcast (dphiinit)
    call broadcast (zf_init)
    call broadcast (apar0)
    call broadcast (tstart)
    call broadcast (chop_side)
    call broadcast (even)
    call broadcast (left)
    call broadcast (max_mode)
    call broadcast (clean_init)
    call broadcast (restart_file)
    call broadcast (proc_to_save_fields)
    call broadcast (read_many)
    call broadcast (ikk)
    call broadcast (itt) 
    call broadcast (ikkk)
    call broadcast (ittt) 
    call broadcast (phiamp)
    call broadcast (aparamp)
    call broadcast (scale)
    call broadcast (new_field_init)

    ! RN>
    call broadcast (phiinit0)
    call broadcast (phiinit_rand)
    call broadcast (a0)
    call broadcast (b0)
    call broadcast (null_phi)
    call broadcast (null_bpar)
    call broadcast (null_apar)
    call broadcast (adj_spec)
    call broadcast (eq_type)
    call broadcast (prof_width)
    call broadcast (eq_mode_n)
    call broadcast (eq_mode_u)
    call broadcast (input_check_recon)
    call broadcast (ukxy_pt)
    
    call broadcast (ikpar_init)
    call broadcast (ikx_init)
    call broadcast (kpar_init)
    call broadcast (restart_eig_id)
    call broadcast (constant_random_flag)
    ! <RN
    call set_restart_file (restart_file)
    if (proc_to_save_fields >= 0) proc_to_save_fields = mod(proc_to_save_fields, nproc)
  end subroutine init_init_g

  !> FIXME : Add documentation  
  subroutine ginit (restarted, override_ginitopt_switch)
    use gs2_save, only: read_t0_from_restart_file
    use species, only: has_hybrid_electron_species, spec
    use mp, only: proc0
    implicit none
    logical, intent (out) :: restarted
    integer, intent(in), optional :: override_ginitopt_switch
    logical, save :: have_warned_about_hybrid = .false.
    real :: t0
    integer :: istatus
    integer :: ginitopt_actual

    if (present(override_ginitopt_switch)) then
      ginitopt_actual = override_ginitopt_switch
    else
      ginitopt_actual = ginitopt_switch
    end if

    restarted = .false.
    select case (ginitopt_actual)
    case (ginitopt_default)
       call ginit_default
    case (ginitopt_default_odd)
       call ginit_default_odd
    case (ginitopt_kz0)
       call ginit_kz0
    case (ginitopt_noise)
       call ginit_noise
    case (ginitopt_single_parallel_mode)
       call ginit_single_parallel_mode
    case (ginitopt_all_modes_equal)
       call ginit_all_modes_equal
    case (ginitopt_kpar)
       call ginit_kpar
    case (ginitopt_gs)
       call ginit_gs
    case (ginitopt_nl)
       call ginit_nl
    case (ginitopt_nl2)
       call ginit_nl2
    case (ginitopt_nl3)
       call ginit_nl3
    case (ginitopt_nl3r)
       call ginit_nl3r
    case (ginitopt_harris)
       call ginit_harris
    case (ginitopt_nl4)
! in an old version, this line was commented out.  Thus, to recover some old
! results, you might need to comment out this line and...
!       t0 = tstart
       call ginit_nl4
       call read_t0_from_restart_file (tstart, istatus)
! this line:
!       tstart = t0
       restarted = .true.
       scale = 1.
    case (ginitopt_nl5)
       t0 = tstart
       call read_t0_from_restart_file (tstart, istatus)
       call ginit_nl5
       tstart = t0
       restarted = .true.
       scale = 1.
    case (ginitopt_recon)
       t0 = tstart
       call read_t0_from_restart_file (tstart, istatus)
       call ginit_recon
       tstart = t0
       restarted = .true.
       scale = 1.
    case (ginitopt_nl6)
       t0 = tstart
       call read_t0_from_restart_file (tstart, istatus)
       call ginit_nl6
       tstart = t0
       restarted = .true.
       scale = 1.
    case (ginitopt_nl7)
       call ginit_nl7
    case (ginitopt_xi)
       call ginit_xi
    case (ginitopt_xi2)
       call ginit_xi2
    case (ginitopt_rh)
       call ginit_rh
    case (ginitopt_alf)
       call ginit_alf
    case (ginitopt_zero)
       call ginit_zero
    case (ginitopt_test3)
       call ginit_test3
    case (ginitopt_convect)
       call ginit_convect
    case (ginitopt_restart_file)
       call ginit_restart_file 
       call read_t0_from_restart_file (tstart, istatus)
       restarted = .true.
       scale = 1.
    case (ginitopt_restart_many)
       call ginit_restart_many 
       call read_t0_from_restart_file (tstart, istatus)
       restarted = .true.
       scale = 1.
    case (ginitopt_restart_eig)
       call ginit_restart_eig 
       call read_t0_from_restart_file (tstart, istatus)
       restarted = .true.
       scale = 1.
    case (ginitopt_restart_small)
       call ginit_restart_small
       call read_t0_from_restart_file (tstart, istatus)
       restarted = .true.
       scale = 1.
    case (ginitopt_restart_memory)
       call ginit_restart_memory
       restarted = .true.       
       scale = 1.
    case (ginitopt_zonal_only)
       call ginit_restart_zonal_only
       call read_t0_from_restart_file (tstart, istatus)
       restarted = .true.
       scale = 1.
    case (ginitopt_no_zonal)
       call ginit_restart_no_zonal
       call read_t0_from_restart_file (tstart, istatus)
       restarted = .true.
       scale = 1.
    case (ginitopt_smallflat)
       call ginit_restart_smallflat
       call read_t0_from_restart_file (tstart, istatus)
       restarted = .true.
       scale = 1.
    case (ginitopt_continue)
       restarted = .true.
       scale = 1.
    case (ginitopt_recon3)
       call ginit_recon3
    case (ginitopt_ot)
       call ginit_ot
    case (ginitopt_stationary_mode)
       call ginit_stationary_mode
    case (ginitopt_cgyro)
       call ginit_cgyro
    case (ginitopt_random_sine)
       call ginit_random_sine
    case (ginitopt_constant)
       call ginit_constant
    end select

    ! If the initial condition is the nonadiabatic dfn then
    ! ensure we've actually set `gnew` correctly.
    ! We avoid doing this if we're restarting as currently we assume
    ! that the restart files always contain our evolved gnew, rather
    ! than the nonadiabatic dfn.
    if (initial_condition_is_nonadiabatic_dfn .and. .not. restarted) then
       call set_gnew_from_nonadiabatic_dfn
    else
       ! If not restarting and have a hybrid electron species then
       ! issue a warning (once only) as recommend initialising the
       ! non-adiabatic dfn directly in this situation as this respects
       ! the hybrid electron h = 0 condition.
       if ((.not. restarted) .and. has_hybrid_electron_species(spec) &
            .and. (.not. have_warned_about_hybrid) .and. proc0) then
          have_warned_about_hybrid = .true.
          write(*,'("Warning : Running with a hybrid electron species but initial_condition_is_nonadiabatic_dfn is false.")')
       end if
    end if
    
  end subroutine ginit

  !> Take the existing gnew value and, assuming it is the nonadiabatic dfn,
  !> calculate the consistent potentials to set the true gnew value.
  subroutine set_gnew_from_nonadiabatic_dfn
    use dist_fn_arrays, only: gnew, g, g_adjust, to_g_gs2
    use dist_fn, only: calculate_potentials_from_nonadiabatic_dfn
    use species, only: spec, has_hybrid_electron_species
    use le_grids, only: is_passing_hybrid_electron
    use gs2_layouts, only: g_lo, ik_idx, il_idx, is_idx
    use fields_arrays, only: phi, apar, bpar, phinew, aparnew, bparnew
    use array_utils, only: copy
    implicit none
    integer :: iglo, ik, il, is

    ! Respect passing hybrid electrons h=0
    if (has_hybrid_electron_species(spec)) then
       do iglo = g_lo%llim_proc, g_lo%ulim_proc
          is = is_idx(g_lo, iglo)
          ik = ik_idx(g_lo, iglo)
          il = il_idx(g_lo, iglo)
          if (is_passing_hybrid_electron(is, ik, il)) gnew(:,:,iglo) = 0.0
       end do
    end if

    ! Calculate the consistent potentials
    call calculate_potentials_from_nonadiabatic_dfn(gnew, phi, apar, bpar)

    ! Adjust the nonadiabatic dfn to give the actual gnew we evolve
    call g_adjust(gnew, phi, bpar, direction = to_g_gs2)

    ! Ensure both timesteps are set correctly
    call copy(gnew, g)
    call copy(phi, phinew)
    call copy(apar, aparnew)
    call copy(bpar, bparnew)
  end subroutine set_gnew_from_nonadiabatic_dfn

  !> FIXME : Add documentation
  subroutine read_parameters(init_g_config_in)
    use file_utils, only: input_unit, error_unit, run_name, input_unit_exist
    use text_options, only: text_option, get_option_value
    use gs2_save, only: read_many, include_explicit_source_in_restart, proc_to_save_fields
    implicit none
    type(init_g_config_type), intent(in), optional :: init_g_config_in

    type (text_option), dimension (*), parameter :: ginitopts = &
         (/ text_option('default', ginitopt_default), &
            text_option('default_odd', ginitopt_default_odd), &
            text_option('noise', ginitopt_noise), &
            text_option('xi', ginitopt_xi), &
            text_option('xi2', ginitopt_xi2), &
            text_option('zero', ginitopt_zero), &
            text_option('test3', ginitopt_test3), &
            text_option('convect', ginitopt_convect), &
            text_option('rh', ginitopt_rh), &
            text_option('restart', ginitopt_restart_many), &
            text_option('many', ginitopt_restart_many), &
            text_option('small', ginitopt_restart_small), &
            text_option('file', ginitopt_restart_file), &
            text_option('cont', ginitopt_continue), &
            text_option('kz0', ginitopt_kz0), &
            text_option('nl', ginitopt_nl), &
            text_option('nl2', ginitopt_nl2), &
            text_option('nl3', ginitopt_nl3), &
            text_option('nl3r', ginitopt_nl3r), &
            text_option('nl4', ginitopt_nl4), &
            text_option('nl5', ginitopt_nl5), &
            text_option('nl6', ginitopt_nl6), &
            text_option('nl7', ginitopt_nl7), &
            text_option('alf', ginitopt_alf), &
            text_option('gs', ginitopt_gs), &
            text_option('kpar', ginitopt_kpar), &
            text_option('smallflat', ginitopt_smallflat), &
            text_option('harris', ginitopt_harris), &
            text_option('recon', ginitopt_recon), &
            text_option('recon3', ginitopt_recon3), &
            text_option('ot', ginitopt_ot), &
            text_option('zonal_only', ginitopt_zonal_only), &
            text_option('no_zonal', ginitopt_no_zonal), &
            text_option('single_parallel_mode', ginitopt_single_parallel_mode), &
            text_option('all_modes_equal', ginitopt_all_modes_equal), &
            text_option('eig_restart', ginitopt_restart_eig), &
            text_option('stationary_mode', ginitopt_stationary_mode), &
            text_option('cgyro', ginitopt_cgyro), &
            text_option('random_sine', ginitopt_random_sine), &
            text_option('constant', ginitopt_constant) &
            /)
    character(20) :: ginit_option

    integer :: ierr

    if (present(init_g_config_in)) init_g_config = init_g_config_in

    ! Smart defaults
    if (.not.init_g_config%is_initialised()) then    
       init_g_config%restart_file = trim(run_name)//".nc"
    end if
    
    call init_g_config%init(name = 'init_g_knobs', requires_index = .false.)

    ! Copy out internal values into module level parameters
    a0 = init_g_config%a0
    adj_spec = init_g_config%adj_spec
    apar0 = init_g_config%apar0
    aparamp = init_g_config%aparamp
    b0 = init_g_config%b0
    chop_side = init_g_config%chop_side
    clean_init = init_g_config%clean_init
    constant_random_flag = init_g_config%constant_random_flag
    den0 = init_g_config%den0
    den1 = init_g_config%den1
    den2 = init_g_config%den2
    dphiinit = init_g_config%dphiinit
    eq_mode_n = init_g_config%eq_mode_n
    eq_mode_u = init_g_config%eq_mode_u
    eq_type = init_g_config%eq_type
    even = init_g_config%even
    ginit_option = init_g_config%ginit_option
    ikk = init_g_config%ikk
    ikkk = init_g_config%ikkk
    ikpar_init = init_g_config%ikpar_init
    ikx_init = init_g_config%ikx_init
    imfac = init_g_config%imfac
    include_explicit_source_in_restart = init_g_config%include_explicit_source_in_restart
    initial_condition_is_nonadiabatic_dfn = init_g_config%initial_condition_is_nonadiabatic_dfn
    input_check_recon = init_g_config%input_check_recon
    itt = init_g_config%itt
    ittt = init_g_config%ittt
    kpar_init = init_g_config%kpar_init
    left = init_g_config%left
    max_mode = init_g_config%max_mode
    new_field_init = init_g_config%new_field_init
    null_apar = init_g_config%null_apar
    null_bpar = init_g_config%null_bpar
    null_phi = init_g_config%null_phi
    phiamp = init_g_config%phiamp
    phifrac = init_g_config%phifrac
    phiinit = init_g_config%phiinit
    phiinit0 = init_g_config%phiinit0
    phiinit_rand = init_g_config%phiinit_rand
    proc_to_save_fields = init_g_config%proc_to_save_fields
    prof_width = init_g_config%prof_width
    read_many = init_g_config%read_many
    refac = init_g_config%refac
    restart_dir = init_g_config%restart_dir
    restart_eig_id = init_g_config%restart_eig_id
    restart_file = init_g_config%restart_file
    scale = init_g_config%scale
    tpar0 = init_g_config%tpar0
    tpar1 = init_g_config%tpar1
    tpar2 = init_g_config%tpar2
    tperp0 = init_g_config%tperp0
    tperp1 = init_g_config%tperp1
    tperp2 = init_g_config%tperp2
    tstart = init_g_config%tstart
    ukxy_pt = init_g_config%ukxy_pt
    upar0 = init_g_config%upar0
    upar1 = init_g_config%upar1
    upar2 = init_g_config%upar2
    width0 = init_g_config%width0
    zf_init = init_g_config%zf_init

    exist = init_g_config%exist
    
    ierr = error_unit()
    call get_option_value &
         (ginit_option, ginitopts, ginitopt_switch, &
         ierr, "ginit_option in ginit_knobs",.true.)
  end subroutine read_parameters

  !> Set initial condition to gaussian in \(\theta\)
  !>
  !> Controlled by:
  !> - Amplitude: [[init_g_knobs:phiinit]]
  !> - Width in \(\theta\): [[init_g_knobs:width0]]
  !> - Parity: even about \(\theta_0\) unless [[init_g_knobs:chop_side]] is true
  subroutine ginit_default
    use species, only: spec
    use theta_grid, only: ntgrid, theta
    use kt_grids, only: naky, ntheta0, theta0, 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
    implicit none
    complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi
    logical :: right
    integer :: iglo
    integer :: ig, ik, it, il, is

    right = .not. left

    do ig = -ntgrid, ntgrid
       phi(ig,:,:) = exp(-((theta(ig)-theta0(:,:))/width0)**2)*cmplx(1.0,1.0)            
    end do
    if (chop_side .and. left) phi(:-1,:,:) = 0.0
    if (chop_side .and. right) phi(1:,:,:) = 0.0
    
    if (reality) then
       phi(:,1,1) = 0.0

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

! not used:
! reality condition for k_theta = 0 component:
       do it = 1, ntheta0/2
          phi(:,it+(ntheta0+1)/2,1) = conjg(phi(:,(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)
       g(:,1,iglo) = phi(:,it,ik)*spec(is)%z*phiinit
       where (forbid(:,il)) g(:,1,iglo) = 0.0
       g(:,2,iglo) = g(:,1,iglo)
    end do
    call copy(g, gnew)
  end subroutine ginit_default

  !> Set initial condition to gaussian in \(\theta\)
  !>
  !> Controlled by:
  !> - Amplitude: [[init_g_knobs:phiinit]]
  !> - Width in \(\theta\): [[init_g_knobs:width0]]
  !> - Parity: odd about \(\theta_0\) unless [[init_g_knobs:chop_side]] is true
  subroutine ginit_default_odd
    use species, only: spec
    use theta_grid, only: ntgrid, theta
    use kt_grids, only: naky, ntheta0, theta0, 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
    implicit none
    complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi
    logical :: right
    integer :: iglo
    integer :: ig, ik, it, il, is

    right = .not. left

    do ig = -ntgrid, ntgrid
       phi(ig,:,:) = sin((theta(ig)-theta0(:,:))/width0)*exp(-((theta(ig)-theta0(:,:))/width0)**2)*cmplx(1.0,1.0)            
    end do
    if (chop_side .and. left) phi(:-1,:,:) = 0.0
    if (chop_side .and. right) phi(1:,:,:) = 0.0
    
    if (reality) then
       phi(:,1,1) = 0.0

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

! not used:
! reality condition for k_theta = 0 component:
       do it = 1, ntheta0/2
          phi(:,it+(ntheta0+1)/2,1) = conjg(phi(:,(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)
       g(:,1,iglo) = phi(:,it,ik)*spec(is)%z*phiinit
       where (forbid(:,il)) g(:,1,iglo) = 0.0
       g(:,2,iglo) = g(:,1,iglo)
    end do
    call copy(g, gnew)
  end subroutine ginit_default_odd

  !> Initialise with only the kparallel = 0 mode. 
  subroutine ginit_kz0
    use species, only: spec
    use theta_grid, only: ntgrid 
    use kt_grids, only: naky, ntheta0, aky, reality
    use le_grids, only: forbid
    use dist_fn_arrays, only: g, gnew, vpa, vperp2
    use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx, is_idx
    use array_utils, only: copy
    implicit none
    complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi
    logical :: right
    integer :: iglo
    integer :: ik, it, il, is

    right = .not. left

    phi = cmplx(refac,imfac)
    if (chop_side .and. left) phi(:-1,:,:) = 0.0
    if (chop_side .and. right) phi(1:,:,:) = 0.0
    
    if (reality) then
       phi(:,1,1) = 0.0

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

! not used:
! reality condition for k_theta = 0 component:
       do it = 1, ntheta0/2
          phi(:,it+(ntheta0+1)/2,1) = conjg(phi(:,(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)
       g(:,1,iglo) = -phi(:,it,ik)*spec(is)%z*phiinit*(den0 + sqrt(2.)*upar0*vpa(:,1,iglo) + tpar0*(vpa(:,1,iglo)**2-0.5) + tperp0*(vperp2(:,iglo)-1.))
       where (forbid(:,il)) g(:,1,iglo) = 0.0
       g(:,2,iglo) = -phi(:,it,ik)*spec(is)%z*phiinit*(den0 + sqrt(2.)*upar0*vpa(:,2,iglo) + tpar0*(vpa(:,2,iglo)**2-0.5) + tperp0*(vperp2(:,iglo)-1.))
    end do
    call copy(g, gnew)
  end subroutine ginit_kz0

  !> FIXME : Add documentation  
  subroutine single_initial_kx(phi)
    use theta_grid, only: ntgrid 
    use kt_grids, only: naky, ntheta0
!   use kt_grids, only: akx
    use mp, only: mp_abort
    implicit none
    complex, dimension (-ntgrid:ntgrid,ntheta0,naky), intent(inout) :: phi
    integer :: ig, ik, it

    if (ikx_init  < 2 .or. ikx_init > (ntheta0+1)/2) then
      call mp_abort("The subroutine single_initial_kx should only be called when 1 < ikx_init < (ntheta0+1)/2")
    end if

    do it = 1, ntheta0
      if (it .ne. ikx_init) then 
        !write (*,*) "zeroing out kx_index: ", it, "at kx: ", akx(it)
         do ik = 1, naky
            do ig = -ntgrid, ntgrid
               phi(ig,it,ik) = 0.
             end do
         end do
       end if
    end do
  end subroutine single_initial_kx

  !> Initialise the distribution function with random noise. This is the
  !> recommended initialisation option. Each different mode is given a random
  !> amplitude between zero and one.
  subroutine ginit_noise
    use species, only: spec, tracer_species
    use theta_grid, only: ntgrid 
    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, proc_id
    use dist_fn, only: boundary_option_linked, boundary_option_switch
    use dist_fn, only: l_links, r_links
    use dist_fn, only: pass_right, init_pass_ends
    use redistribute, only: fill, delete_redist
    use ran, only: ranf
    use constant_random, only: init_constant_random
    use constant_random, only: finish_constant_random
    use constant_random, only: constant_ranf=>ranf
    use array_utils, only: copy, zero_array
    implicit none
    complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi, phit
    real :: a, b
    integer :: iglo, ig, ik, it, il, is, nn


    !CMR: need to document tracer phit parameter   ;-)
    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

    
    ! keep old (it, ik) loop order to get old results exactly: 

    !Fill phi with random (complex) numbers between -0.5 and 0.5
    if (constant_random_flag) call init_constant_random
    do it = 1, ntheta0
       do ik = 1, naky
          do ig = -ntgrid, ntgrid
            if (constant_random_flag) then 
             a = constant_ranf()-0.5
             b = constant_ranf()-0.5
            else
             a = ranf()-0.5
             b = ranf()-0.5
            end if
             phi(ig,it,ik) = cmplx(a,b)
           end do
!CMR,28/1/13: 
! clean_init debrutalises influence of chop_side on phi with linked bc
          if (chop_side) then
             if (left) then
                if (clean_init.and.(boundary_option_switch .eq. boundary_option_linked)) then
                   !Does this tube have more links to the right than the left?
                   !If so then it's on the left so set to zero
                   if (l_links(it, ik) .lt. r_links(it, ik)) then
                      phi(:,it,ik) = 0.0 
                   !This tube is in the middle so only set the left half to 0
                   else if (l_links(it, ik) .eq. r_links(it, ik)) then
                      phi(:-1,it,ik) = 0.0
                   endif
                else
                   phi(:-1,it,ik) = 0.0
                endif
             else
                if (clean_init.and.(boundary_option_switch .eq. boundary_option_linked)) then
                   !Does this tube have more links to the left than the right?
                   !If so then it's on the right so set to zero
                   if (r_links(it, ik) .lt. l_links(it, ik)) then
                      phi(:,it,ik) = 0.0
                   !This tube is in the middle so only set the right half to 0
                   else if (l_links(it, ik) .eq. r_links(it, ik)) then
                      phi(1:,it,ik) = 0.0
                   endif
                else
                   phi(1:,it,ik) = 0.0
                endif
             endif
          end if
       end do
    end do
    if (constant_random_flag) call finish_constant_random
    
    !Wipe out all but one kx if requested
    if (ikx_init  > 0) call single_initial_kx(phi)
    
    !Sort out the zonal/self-periodic modes
    if (naky .ge. 1 .and. aky(1) == 0.0) then
       !Apply scaling factor
       phi(:,:,1) = phi(:,:,1)*zf_init
       
       !Set ky=kx=0.0 mode to zero in amplitude
       phi(:,1,1) = 0.0
       
       !CMR, 25/1/13: clean_init forces periodic zonal flows
       if (clean_init .and. boundary_option_switch.eq.boundary_option_linked) then
          phi(ntgrid,:,1)=phi(-ntgrid,:,1)
          phit(ntgrid,:,1)=phit(-ntgrid,:,1)
       end if
    end if

    !Apply reality condition (i.e. -kx mode is conjugate of +kx mode)
    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

    !Now set g using data in phi
    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)

       !Handle tracer_species 
       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

       !Make sure there's no g in the forbidden regions
       where (forbid(:,il)) g(:,1,iglo) = 0.0
     
       !Set incoming and outgoing boundary conditions
       if ( clean_init .and. boundary_option_switch .eq. boundary_option_linked .and. ik .gt. 1) then
          !CMR: clean_init sets g=0 for incoming particles, 
          !                         and outgoing particles
          !                      as initial condition sets g(:,1,:) = g(:,2,:) )
          !If no links to the left then we're at the left boundary
          if ( l_links(it, ik) .eq. 0 ) g(-ntgrid,1,iglo)=0.0
          !If no links to the right then we're at the right boundary
          if ( r_links(it, ik) .eq. 0 ) g(ntgrid,1,iglo)=0.0
       endif
    end do

    
    !If clean_init is set and there are any links/connections then make sure repeated points are consistent
    if ( clean_init .and. boundary_option_switch .eq. boundary_option_linked .and. sum(r_links+l_links) .gt. 0 ) then
       !   
       !CMR, 23/1/13:  
       !  clean_init: ensure continuity in || direction with linked bcs
       !  set up redistribute type pass_right 
       !  this passes g(ntgrid,1,iglo) to g(-ntgrid,1,iglo_r) on right neighbour 

       !Initialise communications object !NOTE This should be moved to dist_fn in future
       call init_pass_ends(pass_right,'r',1,'c')

       !Pass right boundaries (ntgrid) of g to linked left boundaries (-ntgrid) of g
       call fill(pass_right,g,g)

       !Deallocate comm object as not used anywhere else !NOTE This will need removing if/when initialisation moved to dist_fn
       call delete_redist(pass_right)

       !Make leftwards and rightwards particles the same
       g(:,2,:)=g(:,1,:)

       !Set gnew to be the "fixed" data
       call copy(g, gnew)
    else 
! finally initialise isign=2
       g(:,2,:) = g(:,1,:)
       call copy(g, gnew)
    endif

  end subroutine ginit_noise

  !> Initialise only with a single parallel mode \(k_\Vert\) specified by either
  !> [[init_g_knobs:ikpar_init]] for periodic boundary conditions or
  !> [[init_g_knobs:kpar_init]] for linked boundary conditions. Only makes sense
  !> for linear calculations.
  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

  !> Initialize with every parallel and perpendicular mode having equal amplitude. 
  !! Only makes sense in a linear calculation. k_parallel is specified with kpar_init 
  !! or with ikpar_init when periodic boundary conditions are used. EGH 
  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

  !> "nl": Two \(k_\perp\) modes set to constant amplitude. Controlled by:
  !>
  !>  - Overall amplitude: [[init_g_knobs:phiinit]]
  !>  - `ky` indices of modes: [[init_g_knobs:ikk]]
  !>  - `kx` indices of modes: [[init_g_knobs:itt]]
  !>  - Parity: even about \(\theta_0\) unless
  !>    [[init_g_knobs:chop_side]] is true, in which case none
  subroutine ginit_nl
    use species, only: spec
    use theta_grid, only: ntgrid 
    use kt_grids, only: naky, ntheta0 
    use le_grids, only: forbid, il_is_wfb
    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
    integer :: iglo
    integer :: ig, ik, it, il, is, j
    
    call zero_array(phi)
    do j = 1, 2
       ik = ikk(j)
       it = itt(j)
       do ig = -ntgrid, ntgrid
!          phi(ig,it,ik) = cmplx(1.0, 0.0)*sin(theta(ig))
          phi(ig,it,ik) = cmplx(1.0, 0.0)!*sin(theta(ig))
       end do
!       do ig = -ntgrid/2, ntgrid/2
!          phi(ig,it,ik) = cmplx(1.0, 1.0)
!       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
    
! reality condition for k_theta = 0 component:
    do it = 1, ntheta0/2
       phi(:,it+(ntheta0+1)/2,1) = conjg(phi(:,(ntheta0+1)/2+1-it,1))
    enddo

! charge dependence keeps initial Phi from being too small
    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) = -phi(:,it,ik)*phiinit*spec(is)%z
       where (forbid(:,il)) g(:,1,iglo) = 0.0
       g(:,2,iglo) = g(:,1,iglo)
       if (il_is_wfb(il)) g(:,:,iglo) = 0.0
    end do
    call copy(g, gnew)
  end subroutine ginit_nl

  !> "nl2": Two \(k_\perp\) modes proportional to \(1 + v_\Vert
  !>     \sin(\theta)\). Controlled by:
  !>
  !>  - Overall amplitude: [[init_g_knobs:phiinit]]
  !>  - `ky` indices of modes: [[init_g_knobs:ikk]]
  !>  - `kx` indices of modes: [[init_g_knobs:itt]]
  !>  - Parity: even about \(\theta_0\) unless
  !>    [[init_g_knobs:chop_side]] is true, in which case none
  subroutine ginit_nl2
    use species, only: spec
    use theta_grid, only: ntgrid, theta
    use kt_grids, only: naky, ntheta0 
    use le_grids, only: forbid, il_is_wfb
    use dist_fn_arrays, only: g, gnew, vpa
    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
    integer :: iglo
    integer :: ig, ik, it, il, is, j
    
    call zero_array(phi)
    do j = 1, 2
       ik = ikk(j)
       it = itt(j)
       do ig = -ntgrid, ntgrid
          phi(ig,it,ik) = cmplx(1.0, 0.0)!*sin(theta(ig))
       end do
!       do ig = -ntgrid/2, ntgrid/2
!          phi(ig,it,ik) = cmplx(1.0, 1.0)
!       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
    
! reality condition for k_theta = 0 component:
    do it = 1, ntheta0/2
       phi(:,it+(ntheta0+1)/2,1) = conjg(phi(:,(ntheta0+1)/2+1-it,1))
    enddo
    
! charge dependence keeps initial Phi from being too small
    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) = -phi(:,it,ik)*phiinit*(1.+vpa(:,1,iglo)*sin(theta))*spec(is)%z
       where (forbid(:,il)) g(:,1,iglo) = 0.0
       g(:,2,iglo) = -phi(:,it,ik)*phiinit*(1.+vpa(:,2,iglo)*sin(theta))*spec(is)%z
!       g(:,1,iglo)*vpa(:,2,iglo)
       if (il_is_wfb(il)) g(:,:,iglo) = 0.0
    end do
    call copy(g, gnew)
  end subroutine ginit_nl2

  !> "nl3": Two \(k_\perp\) modes with a Gaussian envelope, and perturbation in
  !> fields of form \(x_0 + x_1 \cos(\theta) + x_2 \cos(2 \theta)\), except for
  !> parallel velocity which uses is 90 degrees out of phase and uses
  !> \(\sin\). Density, parallel velocity, and perpendicular and parallel
  !> temperatures can all be controlled independently. Controlled by:
  !>
  !>   - Overall amplitude: [[init_g_knobs:phiinit]]
  !>   - Amplitude of real and imaginary parts of \(\phi\):
  !>     [[init_g_knobs:refac]], [[init_g_knobs:imfac]]
  !>   - Gaussian envelope width in \(\theta\): [[init_g_knobs:width0]]
  !>   - `ky` indices of modes: [[init_g_knobs:ikk]]
  !>   - `kx` indices of modes: [[init_g_knobs:itt]]
  !>   - Parity: even about \(\theta_0\) unless
  !>     [[init_g_knobs:chop_side]] is true, in which case none
  !>   - Density perturbation: [[init_g_knobs:den0]],
  !>     [[init_g_knobs:den1]], [[init_g_knobs:den2]]
  !>   - Parallel velocity perturbation: [[init_g_knobs:upar0]],
  !>     [[init_g_knobs:upar1]], [[init_g_knobs:upar2]]
  !>   - Parallel temperature perturbation: [[init_g_knobs:tpar0]],
  !>     [[init_g_knobs:tpar1]], [[init_g_knobs:tpar2]]
  !>   - Perpendicular temperature perturbation: [[init_g_knobs:tperp0]],
  !>     [[init_g_knobs:tperp1]], [[init_g_knobs:tperp2]]
  subroutine ginit_nl3
    use species, only: spec, has_electron_species
    use theta_grid, only: ntgrid, theta
    use kt_grids, only: naky, ntheta0, theta0, reality
    use le_grids, only: forbid
    use dist_fn_arrays, only: g, gnew, vpa, vperp2
    use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx, is_idx
    use constants, only: zi
    use array_utils, only: copy, zero_array
    implicit none
    complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi, odd
    real, dimension (-ntgrid:ntgrid) :: dfac, ufac, tparfac, tperpfac, ct, st, c2t, s2t
    integer :: iglo
    integer :: ig, ik, it, il, is, j
    
    call zero_array(phi) ; call zero_array(odd)
    do j = 1, 2
       ik = ikk(j)
       it = itt(j)
       if (width0 > 0.) then
          do ig = -ntgrid, ntgrid
             phi(ig,it,ik) = exp(-((theta(ig)-theta0(it,ik))/width0)**2)*cmplx(refac, imfac)
          end do
       else
          do ig = -ntgrid, ntgrid
             phi(ig,it,ik) = cmplx(refac, imfac)
          end do
       end if
       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

    odd = zi * phi
    
! reality condition for k_theta = 0 component:
    if (reality) then
       do it = 1, ntheta0/2
          phi(:,it+(ntheta0+1)/2,1) = conjg(phi(:,(ntheta0+1)/2+1-it,1))
          odd(:,it+(ntheta0+1)/2,1) = conjg(odd(:,(ntheta0+1)/2+1-it,1))
       enddo
    end if

    if (even) then
       ct = cos(theta)
       st = sin(theta)

       c2t = cos(2.*theta)
       s2t = sin(2.*theta)
    else
       ct = sin(theta)
       st = cos(theta)

       c2t = sin(2.*theta)
       s2t = cos(2.*theta)
    end if

    dfac     = den0   + den1 * ct + den2 * c2t
    ufac     = upar0  + upar1* st + upar2* s2t
    tparfac  = tpar0  + tpar1* ct + tpar2* c2t
    tperpfac = tperp0 + tperp1*ct + tperp2*c2t


! charge dependence keeps initial Phi from being too small
    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 (is_electron_species(spec(is))) cycle
       g(:,1,iglo) = phiinit* &!spec(is)%z* &
            ( dfac*spec(is)%dens0            * phi(:,it,ik) &
            + 2.*ufac* vpa(:,1,iglo)         * odd(:,it,ik) &
            + tparfac*(vpa(:,1,iglo)**2-0.5) * phi(:,it,ik) &
            +tperpfac*(vperp2(:,iglo)-1.)    * phi(:,it,ik))
       where (forbid(:,il)) g(:,1,iglo) = 0.0

       g(:,2,iglo) = phiinit* &!spec(is)%z* &
            ( dfac*spec(is)%dens0            * phi(:,it,ik) &
            + 2.*ufac* vpa(:,2,iglo)         * odd(:,it,ik) &
            + tparfac*(vpa(:,2,iglo)**2-0.5) * phi(:,it,ik) &
            +tperpfac*(vperp2(:,iglo)-1.)    * phi(:,it,ik))
       where (forbid(:,il)) g(:,2,iglo) = 0.0

!       if (il == ng2+1) g(:,:,iglo) = 0.0
    end do

!    if (has_electron_species(spec)) then
!       call flae (g, gnew)
!       g = g - gnew
!    end if
    call copy(g, gnew)
  end subroutine ginit_nl3

  !> FIXME : Add documentation  
  subroutine ginit_nl3r
    use species, only: spec, has_electron_species
    use theta_grid, only: ntgrid, theta
    use kt_grids, only: naky, ntheta0, theta0, reality
    use le_grids, only: forbid
    use fields_arrays, only: apar
    use fields_arrays, only: aparnew
    use dist_fn_arrays, only: g, gnew, vpa, vperp2
    use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx, is_idx
    use constants, only: zi
    use array_utils, only: copy, zero_array
    implicit none
    complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi, odd
    real, dimension (-ntgrid:ntgrid) :: dfac, ufac, tparfac, tperpfac, ct, st, c2t, s2t
    integer :: iglo
    integer :: ig, ik, it, il, is, j
    
    call zero_array(phi) ; call zero_array(odd)
    do j = 1, 2
       ik = ikk(j)
       it = itt(j)
       if (width0 > 0.) then
          do ig = -ntgrid, ntgrid
             phi(ig,it,ik) = exp(-((theta(ig)-theta0(it,ik))/width0)**2)*cmplx(refac, imfac)
          end do
       else
          do ig = -ntgrid, ntgrid
             phi(ig,it,ik) = cmplx(refac, imfac)
             apar(ig,it,ik) = apar0*cmplx(refac, imfac)
          end do
       end if
       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

    odd = zi * phi
    
! reality condition for k_theta = 0 component:
    if (reality) then
       do it = 1, ntheta0/2
          apar(:,it+(ntheta0+1)/2,1) = conjg(apar(:,(ntheta0+1)/2+1-it,1))
          phi(:,it+(ntheta0+1)/2,1) = conjg(phi(:,(ntheta0+1)/2+1-it,1))
          odd(:,it+(ntheta0+1)/2,1) = conjg(odd(:,(ntheta0+1)/2+1-it,1))
       enddo
    end if

    aparnew = apar

    if (even) then
       ct = cos(theta)
       st = sin(theta)

       c2t = cos(2.*theta)
       s2t = sin(2.*theta)
    else
       ct = sin(theta)
       st = cos(theta)

       c2t = sin(2.*theta)
       s2t = cos(2.*theta)
    end if

    dfac     = den0   + den1 * ct + den2 * c2t
    ufac     = upar0  + upar1* st + upar2* s2t
    tparfac  = tpar0  + tpar1* ct + tpar2* c2t
    tperpfac = tperp0 + tperp1*ct + tperp2*c2t


! charge dependence keeps initial Phi from being too small
    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) = phiinit* &!spec(is)%z* &
            ( dfac*spec(is)%dens0                * phi(:,it,ik) &
            + 2.*ufac* vpa(:,1,iglo)*spec(is)%u0 * odd(:,it,ik) &
            + tparfac*(vpa(:,1,iglo)**2-0.5)     * phi(:,it,ik) &
            +tperpfac*(vperp2(:,iglo)-1.)        * phi(:,it,ik))
       where (forbid(:,il)) g(:,1,iglo) = 0.0

       g(:,2,iglo) = phiinit* &!spec(is)%z* &
            ( dfac*spec(is)%dens0                * phi(:,it,ik) &
            + 2.*ufac* vpa(:,2,iglo)*spec(is)%u0 * odd(:,it,ik) &
            + tparfac*(vpa(:,2,iglo)**2-0.5)     * phi(:,it,ik) &
            +tperpfac*(vperp2(:,iglo)-1.)        * phi(:,it,ik))
       where (forbid(:,il)) g(:,2,iglo) = 0.0

!       if (il == ng2+1) g(:,:,iglo) = 0.0
    end do

!    if (has_electron_species(spec)) then
!       call flae (g, gnew)
!       g = g - gnew
!    end if
    call copy(g, gnew)
  end subroutine ginit_nl3r

  !> FIXME : Add documentation  
  subroutine ginit_harris
    use species, only: spec, has_electron_species
    use kt_grids, only: naky, ntheta0, reality, nx
    use le_grids, only: forbid
    use dist_fn_arrays, only: g, gnew, vpa, aj0
    use gs2_layouts, only: g_lo, ik_idx, it_idx, is_idx, il_idx
    use constants, only: zi, pi
    use run_parameters, only: k0
    use array_utils, only: copy, zero_array
    implicit none
    complex, dimension (ntheta0,naky) :: phi
    integer :: iglo
    integer :: ik, it, il, is, j, j1
    real, dimension(nx) :: lx_pr, a_pr
    complex,dimension(nx) :: ff_pr
    real:: L, dx_pr
! 
! Specifying function on x grid, with nx points.  Transforming to kx grid, with 
! nkx < nx points.  Result is function that will be used for initial condition.
! But Fourier transform is the same if you just use nkx points in the first 
! place, right?
!    

! Can specify x0 along with y0; no need to use this construction, although it is okay.
    L=2.*pi*k0
    
! nx is available from kt_grids:
    dx_pr=L/real(nx)
    
    do j = 1, nx
       lx_pr(j)=dx_pr*real(j-1)
    end do
    
    do j = 1,nx
       a_pr(j)=1./cosh((lx_pr(j)-L/4.)/width0)- &
         1./cosh((lx_pr(j)-3.*L/4.)/width0)
    end do
    
    a_pr=a_pr/real(nx)
        
    do j = 1,nx
       ff_pr(j) = 0.
       do j1 = 1,nx
    
          ff_pr(j)=ff_pr(j)+a_pr(j1)*exp(-zi*2.*pi*real(j-1)*real(j1-1)/real(nx))
    
       end do
    end do
    
    call zero_array(phi)

! Dealiasing here:
    do j = 1, ntheta0/2-1
       phi(j,1) = ff_pr(j)
    end do

!
! Note: if the j=1 coefficient is non-zero, it might be a problem, because this 
! coefficient is never used.
!
    
    do j = ntheta0/2, ntheta0
       phi(j,1) = ff_pr(nx-ntheta0+j)
    end do

!
! presumably this is not needed, but leave it in for now:
!
! reality condition for k_theta = 0 component:
    if (reality) then
       do it = 1, ntheta0/2
          phi(it+(ntheta0+1)/2,1) = conjg(phi((ntheta0+1)/2+1-it,1))
       enddo
    end if
    
    call zero_array(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)       
    
! ions, kx/=0:
       if ((is==1) .and.(.not.(it==1))) then
       
          g(:,1,iglo) =  2.* vpa(:,1,iglo)*spec(is)%u0 * phi(it,ik)/aj0(:,iglo)  
          g(:,2,iglo) =  2.* vpa(:,2,iglo)*spec(is)%u0 * phi(it,ik)/aj0(:,iglo)

          where (forbid(:,il)) g(:,1,iglo) = 0.0
          where (forbid(:,il)) g(:,2,iglo) = 0.0
          
       end if
    end do
    call copy(g, gnew)
  end subroutine ginit_harris

  !> FIXME : Add documentation  
  subroutine ginit_recon
    use mp, only: proc0
    use species, only: spec, has_electron_species
    use theta_grid, only: ntgrid, theta
    use kt_grids, only: naky, ntheta0, reality
    use le_grids, only: forbid
    use gs2_save, only: gs2_restore
    use dist_fn_arrays, only: g, gnew, vpa, vperp2
    use fields_arrays, only: phi, apar, bpar
    use fields_arrays, only: phinew, aparnew, bparnew
    use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx, is_idx
    use file_utils, only: error_unit
    use run_parameters, only: has_phi, has_apar, has_bpar
    use constants, only: zi
    use array_utils, only: copy, zero_array
    implicit none
    complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phiz, odd
    real, dimension (-ntgrid:ntgrid) :: dfac, ufac, tparfac, tperpfac, ct, st, c2t, s2t
    integer :: iglo, istatus, ierr
    integer :: ig, ik, it, il, is, j
    
    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
       ik = ik_idx(g_lo,iglo)
       if (ik == 1) cycle

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

    phinew(:,:,2:naky) = 0.
    aparnew(:,:,2:naky) = 0.
    bparnew(:,:,2:naky) = 0.
    phi(:,:,2:naky) = 0.
    apar(:,:,2:naky) = 0.
    bpar(:,:,2:naky) = 0.

    call zero_array(phiz) ; call zero_array(odd)
    do j = 1, 2
       ik = ikk(j)
       it = itt(j)
       do ig = -ntgrid, ntgrid
          phiz(ig,it,ik) = cmplx(refac, imfac)
       end do
    end do

    odd = zi * phiz
    
! reality condition for k_theta = 0 component:
    if (reality) then
       do it = 1, ntheta0/2
          phiz(:,it+(ntheta0+1)/2,1) = conjg(phiz(:,(ntheta0+1)/2+1-it,1))
          odd (:,it+(ntheta0+1)/2,1) = conjg(odd (:,(ntheta0+1)/2+1-it,1))
       enddo
    end if

    if (even) then
       ct = cos(theta)
       st = sin(theta)

       c2t = cos(2.*theta)
       s2t = sin(2.*theta)
    else
       ct = sin(theta)
       st = cos(theta)

       c2t = sin(2.*theta)
       s2t = cos(2.*theta)
    end if

    dfac     = den0   + den1 * ct + den2 * c2t
    ufac     = upar0  + upar1* st + upar2* s2t
    tparfac  = tpar0  + tpar1* ct + tpar2* c2t
    tperpfac = tperp0 + tperp1*ct + tperp2*c2t


! charge dependence keeps initial Phi from being too small
    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) = phiinit* &!spec(is)%z* &
            ( dfac*spec(is)%dens0            * phiz(:,it,ik) &
            + 2.*ufac* vpa(:,1,iglo)         * odd (:,it,ik) &
            + tparfac*(vpa(:,1,iglo)**2-0.5) * phiz(:,it,ik) &
            +tperpfac*(vperp2(:,iglo)-1.)    * phiz(:,it,ik))
       where (forbid(:,il)) g(:,1,iglo) = 0.0

       g(:,2,iglo) = phiinit* &!spec(is)%z* &
            ( dfac*spec(is)%dens0            * phiz(:,it,ik) &
            + 2.*ufac* vpa(:,2,iglo)         * odd (:,it,ik) &
            + tparfac*(vpa(:,2,iglo)**2-0.5) * phiz(:,it,ik) &
            +tperpfac*(vperp2(:,iglo)-1.)    * phiz(:,it,ik))
       where (forbid(:,il)) g(:,2,iglo) = 0.0

!       if (il == ng2+1) g(:,:,iglo) = 0.0
    end do

!    if (has_electron_species(spec)) then
!       call flae (g, gnew)
!       g = g - gnew
!    end if
    call copy(g, gnew)
  end subroutine ginit_recon

  !> FIXME : Add documentation
  !!
  !! @note nl4 has been monkeyed with over time.  Do not expect to recover old results
  !! that use this startup routine.
  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

  !> FIXME : Add documentation  
  subroutine ginit_nl5
    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
       ik = ik_idx(g_lo,iglo)
       if (ik == 1) cycle

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

    phinew(:,:,2:naky) = 0.
    aparnew(:,:,2:naky) = 0.
    bparnew(:,:,2:naky) = 0.
    phi(:,:,2:naky) = 0.
    apar(:,:,2:naky) = 0.
    bpar(:,:,2:naky) = 0.

    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

    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)*phiinit*spec(is)%z
       g(:,2,iglo) = g(:,2,iglo)-phiz(:,it,ik)*phiinit*spec(is)%z
       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_nl5

  !> FIXME : Add documentation  
  subroutine ginit_nl6
    use mp, only: proc0
    use gs2_save, only: gs2_restore
!    use theta_grid, only: ntgrid
    use kt_grids, only: ntheta0
!   use kt_grids, only: naky
    use dist_fn_arrays, only: g, gnew
    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 array_utils, only: copy, zero_array
    implicit none
    integer :: iglo, istatus
    integer :: ik, it, 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

    call zero_array(gnew)

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

       g (:,:,iglo) = gnew(:,:,iglo)
    end do

    phinew(:,2,1) = 0.   
    aparnew(:,2,1) = 0.
    bparnew(:,2,1) = 0.

    phi(:,2,1) = 0.
    apar(:,2,1) = 0.
    bpar(:,2,1) = 0.

    phinew(:,ntheta0,1) = 0.
    aparnew(:,ntheta0,1) = 0.
    bparnew(:,ntheta0,1) = 0.

    phi(:,ntheta0,1) = 0.
    apar(:,ntheta0,1) = 0.
    bpar(:,ntheta0,1) = 0.

    call copy(g, gnew)
  end subroutine ginit_nl6

  !> FIXME : Add documentation  
  subroutine ginit_nl7
    use species, only: spec, has_electron_species
    use theta_grid, only: ntgrid, theta
    use kt_grids, only: naky, ntheta0, reality
    use le_grids, only: forbid
    use dist_fn_arrays, only: g, gnew, vpa, vperp2
    use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx, is_idx
    use constants, only: zi
    use array_utils, only: copy, zero_array
    implicit none
    complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi, odd
    real, dimension (-ntgrid:ntgrid) :: dfac, ufac, tparfac, tperpfac, ct, st, c2t, s2t
    integer :: iglo
    integer :: ig, ik, it, il, is, j
    
    call zero_array(phi) ; call zero_array(odd)
    do it=1,ntheta0
       do ik=1,naky
          phi(:,it,ik) = cmplx(refac,imfac)*dphiinit
       end do
    end do

    do j = 1, 2
       ik = ikk(j)
       it = itt(j)
       do ig = -ntgrid, ntgrid
          phi(ig,it,ik) = cmplx(refac, imfac)
       end do
    end do

    odd = zi * phi
    
! reality condition for k_theta = 0 component:
    if (reality) then
       do it = 1, ntheta0/2
          phi(:,it+(ntheta0+1)/2,1) = conjg(phi(:,(ntheta0+1)/2+1-it,1))
          odd(:,it+(ntheta0+1)/2,1) = conjg(odd(:,(ntheta0+1)/2+1-it,1))
       enddo
    end if

    if (even) then
       ct = cos(theta)
       st = sin(theta)

       c2t = cos(2.*theta)
       s2t = sin(2.*theta)
    else
       ct = sin(theta)
       st = cos(theta)

       c2t = sin(2.*theta)
       s2t = cos(2.*theta)
    end if

    dfac     = den0   + den1 * ct + den2 * c2t
    ufac     = upar0  + upar1* st + upar2* s2t
    tparfac  = tpar0  + tpar1* ct + tpar2* c2t
    tperpfac = tperp0 + tperp1*ct + tperp2*c2t

    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) = phiinit* &
            ( dfac*spec(is)%dens0            * phi(:,it,ik) &
            + 2.*ufac* vpa(:,1,iglo)         * odd(:,it,ik) &
            + tparfac*(vpa(:,1,iglo)**2-0.5) * phi(:,it,ik) &
            +tperpfac*(vperp2(:,iglo)-1.)    * phi(:,it,ik))
       where (forbid(:,il)) g(:,1,iglo) = 0.0

       g(:,2,iglo) = phiinit* &
            ( dfac*spec(is)%dens0            * phi(:,it,ik) &
            + 2.*ufac* vpa(:,2,iglo)         * odd(:,it,ik) &
            + tparfac*(vpa(:,2,iglo)**2-0.5) * phi(:,it,ik) &
            +tperpfac*(vperp2(:,iglo)-1.)    * phi(:,it,ik))
       where (forbid(:,il)) g(:,2,iglo) = 0.0

    end do
    call copy(g, gnew)
  end subroutine ginit_nl7

  !> Orszag-Tang 2D vortex problem
  subroutine ginit_ot
    use species, only: spec
    use theta_grid, only: ntgrid
    use kt_grids, only: naky, nakx => ntheta0, reality, kperp2
    use dist_fn_arrays, only: g, gnew, vpa
    use gs2_layouts, only: g_lo, ik_idx, it_idx, is_idx
    use fields_arrays, only: phinew, aparnew, bparnew
    use dist_fn, only: get_init_field
    use array_utils, only: copy, zero_array
    implicit none
    integer :: iglo, ik, it, is, i
    real :: fac
    complex, dimension (-ntgrid:ntgrid,nakx,naky) :: phi, jpar ! local !
    real, dimension (-ntgrid:ntgrid) :: dfac, ufac

    !! phi, jpar are local !!
    call zero_array(phi) ; call zero_array(jpar)
!!$    phi(:,1,2) = phiinit * cmplx(2.0, 0.0)  ! 2 cos(y)
!!$    phi(:,2,1) = phiinit * cmplx(1.0, 0.0)  ! 2 cos(x)
!!$    jpar(:,1,2) = apar0 * cmplx(2.0, 0.0) ! 2 cos(y)
!!$    jpar(:,3,1) = apar0 * cmplx(2.0, 0.0) ! 4 cos(2x)
    do i=1, 3
       it = ittt(i)
       ik = ikkk(i)
       phi(:,it,ik) = phiamp(i)
       jpar(:,it,ik) = aparamp(i) * kperp2(:,it,ik)
    end do

! reality condition for ky = 0 component:
    if (reality) then
       do it = 1, nakx/2
          phi(:,it+(nakx+1)/2,1) = conjg(phi(:,(nakx+1)/2+1-it,1))
          jpar(:,it+(nakx+1)/2,1) = conjg(jpar(:,(nakx+1)/2+1-it,1))
       end do
    end if

    dfac     = den0
    ufac     = upar0

    do iglo = g_lo%llim_proc, g_lo%ulim_proc
       it = it_idx(g_lo,iglo)
       ik = ik_idx(g_lo,iglo)
       is = is_idx(g_lo,iglo)

       g(:,1,iglo) = &
            ( dfac*spec(is)%dens0                * phi(:,it,ik) &
            + 2.*ufac* vpa(:,1,iglo)*spec(is)%u0 * jpar(:,it,ik) )

       g(:,2,iglo) = &
            ( dfac*spec(is)%dens0                * phi(:,it,ik) &
            + 2.*ufac* vpa(:,2,iglo)*spec(is)%u0 * jpar(:,it,ik) )

    end do

    call copy(g, gnew)

    ! normalize
    call get_init_field (phinew, aparnew, bparnew)
    do i=1, 3
       it = ittt(i)
       ik = ikkk(i)
       if (abs(phiamp(i)) > epsilon(0.0)) then
          fac = phiamp(i) / phinew(0,it,ik)
          phi(:,it,ik) = phi(:,it,ik) * fac
       end if
       if (abs(aparamp(i)) > epsilon(0.0)) then
          fac = aparamp(i) / aparnew(0,it,ik)
          jpar(:,it,ik) = jpar(:,it,ik) * fac
       end if
    end do

    ! redefine g
    if (reality) then
       do it = 1, nakx/2
          phi(:,it+(nakx+1)/2,1) = conjg(phi(:,(nakx+1)/2+1-it,1))
          jpar(:,it+(nakx+1)/2,1) = conjg(jpar(:,(nakx+1)/2+1-it,1))
       end do
    end if
    do iglo = g_lo%llim_proc, g_lo%ulim_proc
       it = it_idx(g_lo,iglo)
       ik = ik_idx(g_lo,iglo)
       is = is_idx(g_lo,iglo)

       g(:,1,iglo) = &
            ( dfac*spec(is)%dens0                * phi(:,it,ik) &
            + 2.*ufac* vpa(:,1,iglo)*spec(is)%u0 * jpar(:,it,ik) )

       g(:,2,iglo) = &
            ( dfac*spec(is)%dens0                * phi(:,it,ik) &
            + 2.*ufac* vpa(:,2,iglo)*spec(is)%u0 * jpar(:,it,ik) )

    end do
    call copy(g, gnew)
  end subroutine ginit_ot

  !> FIXME : Add documentation
  subroutine ginit_kpar
    use species, only: spec, has_electron_species
    use theta_grid, only: ntgrid, theta
    use kt_grids, only: naky, ntheta0, theta0
    use le_grids, only: forbid
    use dist_fn_arrays, only: g, gnew, vpa, vperp2
    use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx, is_idx
    use constants, only: zi
    use array_utils, only: copy, zero_array
    implicit none
    complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi, odd
    real, dimension (-ntgrid:ntgrid) :: dfac, ufac, tparfac, tperpfac
    integer :: iglo
    integer :: ig, ik, it, il, is
    
    call zero_array(phi) ; call zero_array(odd)
    if (width0 > 0.) then
       do ig = -ntgrid, ntgrid
          phi(ig,:,:) = exp(-((theta(ig)-theta0(:,:))/width0)**2)*cmplx(refac, imfac)
       end do
    else
       do ig = -ntgrid, ntgrid
          phi(ig,:,:) = cmplx(refac, imfac)
       end do
    end if
    if (chop_side) then
       if (left) then
          phi(:-1,:,:) = 0.0
       else
          phi(1:,:,:) = 0.0
       end if
    end if

    odd = zi * phi
        
    dfac     = den0   + den1 * cos(theta) + den2 * cos(2.*theta) 
    ufac     = upar0  + upar1* sin(theta) + upar2* sin(2.*theta) 
    tparfac  = tpar0  + tpar1* cos(theta) + tpar2* cos(2.*theta) 
    tperpfac = tperp0 + tperp1*cos(theta) + tperp2*cos(2.*theta) 

! charge dependence keeps initial Phi from being too small
    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) = phiinit* &!spec(is)%z* &
            ( dfac                           * phi(:,it,ik) &
            + 2.*ufac* vpa(:,1,iglo)         * odd(:,it,ik) &
            + tparfac*(vpa(:,1,iglo)**2-0.5) * phi(:,it,ik) &
            +tperpfac*(vperp2(:,iglo)-1.)    * phi(:,it,ik))
       where (forbid(:,il)) g(:,1,iglo) = 0.0

       g(:,2,iglo) = phiinit* &!spec(is)%z* &
            ( dfac                           * phi(:,it,ik) &
            + 2.*ufac* vpa(:,2,iglo)         * odd(:,it,ik) &
            + tparfac*(vpa(:,2,iglo)**2-0.5) * phi(:,it,ik) &
            +tperpfac*(vperp2(:,iglo)-1.)    * phi(:,it,ik))
       where (forbid(:,il)) g(:,2,iglo) = 0.0

    end do

    if (has_electron_species(spec)) then
       call flae (g, gnew)
       g = g - gnew
    end if
    call copy(g, gnew)
  end subroutine ginit_kpar

  !> FIXME : Add documentation  
  subroutine ginit_gs
    use species, only: spec, has_electron_species
    use theta_grid, only: ntgrid, theta
    use kt_grids, only: naky, ntheta0
    use le_grids, only: forbid
    use dist_fn_arrays, only: g, gnew, vpa, vperp2
    use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx, is_idx
    use constants, only: pi, zi
    use ran, only: ranf
    use array_utils, only: copy, zero_array
    implicit none
    complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi, odd
    integer :: iglo
    integer :: ik, it, il, is
    real :: phase
    
    call zero_array(phi) ; call zero_array(odd)
    do ik=1,naky
       do it=1,ntheta0
          phase = 2.*pi*ranf()
          phi(:,it,ik) = cos(theta+phase)*cmplx(refac,imfac)
          odd(:,it,ik) = sin(theta+phase)*cmplx(refac,imfac) * zi
       end do
    end do
  
!<DD>Should this be triggered for kt_grids::reality=.true. only?
  
! reality condition for k_theta = 0 component:
    do it = 1, ntheta0/2
       phi(:,it+(ntheta0+1)/2,1) = conjg(phi(:,(ntheta0+1)/2+1-it,1))
       odd(:,it+(ntheta0+1)/2,1) = conjg(odd(:,(ntheta0+1)/2+1-it,1))
    enddo

! charge dependence keeps initial Phi from being too small
    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) = phiinit* &!spec(is)%z* &
            ( den1                         * phi(:,it,ik) &
            + 2.*upar1* vpa(:,1,iglo)      * odd(:,it,ik) &
            + tpar1*(vpa(:,1,iglo)**2-0.5) * phi(:,it,ik) &
            + tperp1*(vperp2(:,iglo)-1.)   * phi(:,it,ik))
       where (forbid(:,il)) g(:,1,iglo) = 0.0
       
       g(:,2,iglo) = phiinit* &!spec(is)%z* &
            ( den1                         * phi(:,it,ik) &
            + 2.*upar1* vpa(:,2,iglo)      * odd(:,it,ik) &
            + tpar1*(vpa(:,2,iglo)**2-0.5) * phi(:,it,ik) &
            + tperp1*(vperp2(:,iglo)-1.)   * phi(:,it,ik))
       where (forbid(:,il)) g(:,2,iglo) = 0.0

    end do

    if (has_electron_species(spec)) then
       call flae (g, gnew)
       g = g - gnew
    end if

    call copy(g, gnew)
  end subroutine ginit_gs

  !> FIXME : Add documentation  
  subroutine ginit_xi
    use theta_grid, only: ntgrid, theta, bmag
    use le_grids, only: forbid, al
    use kt_grids, only: theta0
    use dist_fn_arrays, only: g, gnew
    use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx
    use array_utils, only: copy
    implicit none
    integer :: iglo
    integer :: ik, it, il
    real, dimension(-ntgrid:ntgrid) :: xi

    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)
       xi = sqrt(max(1.0-bmag*al(il),0.0))
       g(:,1,iglo) = xi*exp(-((theta-theta0(it,ik))/width0)**2)
       where (forbid(:,il)) g(:,1,iglo) = 0.0
       g(:,2,iglo) = -g(:,1,iglo)
    end do
    call copy(g, gnew)
  end subroutine ginit_xi

  !> FIXME : Add documentation  
  subroutine ginit_xi2
    use theta_grid, only: ntgrid, theta, bmag
    use le_grids, only: forbid, al
    use kt_grids, only: theta0
    use dist_fn_arrays, only: g, gnew
    use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx
    use array_utils, only: copy
    implicit none
    integer :: iglo
    integer :: ik, it, il
    real, dimension(-ntgrid:ntgrid) :: xi

    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)
       xi = sqrt(max(1.0-bmag*al(il),0.0))
       g(:,1,iglo) = (1.0 - 3.0*xi*xi)*exp(-((theta-theta0(it,ik))/width0)**2)
       where (forbid(:,il)) g(:,1,iglo) = 0.0
       g(:,2,iglo) = g(:,1,iglo)
    end do
    call copy(g, gnew)
  end subroutine ginit_xi2

  !> FIXME : Add documentation  
  subroutine ginit_rh
    use le_grids, only: forbid, energy
    use dist_fn_arrays, only: g, gnew
    use gs2_layouts, only: g_lo, il_idx, ie_idx, is_idx
    use array_utils, only: copy    
    implicit none
    integer :: iglo
    integer :: il, ie, is

    do iglo = g_lo%llim_proc, g_lo%ulim_proc
       il = il_idx(g_lo,iglo)
       ie = ie_idx(g_lo,iglo)
       is = is_idx(g_lo,iglo)
       g(:,1,iglo) = exp(-energy(ie,is))
       where (forbid(:,il)) g(:,1,iglo) = 0.0
       g(:,2,iglo) = g(:,1,iglo)
    end do
    call copy(g, gnew)
  end subroutine ginit_rh

  !> FIXME : Add documentation
  subroutine ginit_alf
    use theta_grid, only: theta
    use le_grids, only: forbid
    use dist_fn_arrays, only: g, gnew, vpa
    use gs2_layouts, only: g_lo, il_idx, is_idx
    use species, only: spec, is_electron_species
    use array_utils, only: copy
    implicit none
    integer :: iglo
    integer :: il, is

    do iglo = g_lo%llim_proc, g_lo%ulim_proc
       is = is_idx(g_lo,iglo)
       if (is_electron_species(spec(is))) cycle
       il = il_idx(g_lo,iglo)
       g(:,1,iglo) = sin(theta)*vpa(:,1,iglo)*spec(is)%z
       where (forbid(:,il)) g(:,1,iglo) = 0.0
       g(:,2,iglo) = sin(theta)*vpa(:,2,iglo)*spec(is)%z
       where (forbid(:,il)) g(:,2,iglo) = 0.0
    end do
    g = phiinit * g 
    call copy(g, gnew)
  end subroutine ginit_alf

  !> FIXME : Add documentation
  subroutine ginit_zero
    use dist_fn_arrays, only: g, gnew
    use array_utils, only: zero_array
    implicit none
    call zero_array(g)
    call zero_array(gnew)
  end subroutine ginit_zero

  !> FIXME : Add documentation  
  subroutine ginit_test3
    use dist_fn_arrays, only: g, gnew, vpa
    use theta_grid, only: ntgrid, delthet, bmag
    use kt_grids, only: akx
    use theta_grid_params, only: eps, epsl, pk
    use gs2_layouts, only: g_lo, ik_idx, il_idx
    use mp, only: broadcast
    use constants, only: zi
    use array_utils, only: copy
    implicit none
    integer :: iglo, ik, il
    real :: c1, c2

    call broadcast (epsl)
    call broadcast (eps)
    call broadcast (pk)

    do iglo = g_lo%llim_proc, g_lo%ulim_proc
       ik = ik_idx(g_lo,iglo)
       il = il_idx(g_lo,iglo)
       if (any(vpa(-ntgrid:ntgrid-1,:,iglo) == 0.0)) then
          c1 = 0.0
          c2 = 0.0
       else
          c2 = -akx(ik)*epsl/eps/pk &
               *sum(delthet(-ntgrid:ntgrid-1)/bmag(-ntgrid:ntgrid-1))
          c1 = c2/sum(delthet(-ntgrid:ntgrid-1)/vpa(-ntgrid:ntgrid-1,1,iglo))
          c2 = c2/sum(delthet(-ntgrid:ntgrid-1)/vpa(-ntgrid:ntgrid-1,2,iglo))
       end if
       g(:,1,iglo) = -zi*akx(ik)*epsl/eps/pk*vpa(:,1,iglo)/bmag - zi*c1
       g(:,2,iglo) = -zi*akx(ik)*epsl/eps/pk*vpa(:,2,iglo)/bmag - zi*c2
    end do
    call copy(g, gnew)
  end subroutine ginit_test3

  !> FIXME : Add documentation  
  subroutine ginit_convect
    use dist_fn_arrays, only: g, gnew
    use theta_grid, only: theta
    use kt_grids, only: theta0
    use gs2_layouts, only: g_lo, it_idx, ik_idx
    use run_parameters, only: k0
    use array_utils, only: copy
    implicit none
    integer :: it, ik, iglo

    do iglo = g_lo%llim_proc, g_lo%ulim_proc
       it = it_idx(g_lo,iglo)
       ik = ik_idx(g_lo,iglo)
       g(:,1,iglo) = exp(cmplx(-(theta-theta0(it,ik))**2,k0*theta))
       g(:,2,iglo) = g(:,1,iglo)
    end do
    call copy(g, gnew)
  end subroutine ginit_convect

  !> FIXME : Add documentation  
  subroutine ginit_recon3
    use mp, only: proc0, mp_abort
    use species, only: nspec, spec, has_electron_species
    use theta_grid, only: ntgrid, theta
    use kt_grids, only: naky, nakx => ntheta0, akx, aky, reality
    use kt_grids, only: nx,ny,kperp2
    use dist_fn_arrays, only: g, gnew, vpa, vperp2
    use gs2_layouts, only: g_lo, ik_idx, it_idx, is_idx
    use gs2_transforms, only: inverse2,transform2
    use run_parameters, only: beta
    use le_grids, only: integrate_moment
    use dist_fn, only: getmoms_notgc, get_init_field
    use dist_fn, only: init_mom_coeff
    use dist_fn, only: gamtot, gamtot1, gamtot2
    use fields_arrays, only: phinew, bparnew
    use constants, only: pi, zi
    use file_utils, only: error_unit, open_output_file, close_output_file
    use ran, only: ranf
    use array_utils, only: copy
    implicit none
    real, dimension(:, :, :, :), allocatable :: mom_coeff
    real, dimension(:, :, :), allocatable :: mom_coeff_npara, mom_coeff_nperp
    real, dimension(:, :, :), allocatable :: mom_coeff_tpara, mom_coeff_tperp
    real, dimension(:, :, :), allocatable :: mom_shift_para, mom_shift_perp
    integer :: iglo
    integer :: ig, ik, it, is
    integer :: i,j

    ! nkxyz and ukxyz determine profiles in kx-ky plane
    ! nkxyz is for N, T, and ukxyz is for U
    ! [do not control amplitude by these variables]
    complex :: nkxyz(-ntgrid:ntgrid,nakx,naky)
    complex :: ukxyz(-ntgrid:ntgrid,nakx,naky)
    complex :: nkxyz_eq(-ntgrid:ntgrid,nakx,naky)
    complex :: ukxyz_eq(-ntgrid:ntgrid,nakx,naky)
    complex, allocatable :: g_eq(:,:,:)

    ! equilibrium and perturbation
    complex :: nkxy_eq(nakx,naky), ukxy_eq(nakx,naky)

    ! *fac determine z profile
    ! [do not control amplitude by these variables]
    real :: dfac(-ntgrid:ntgrid), ufac(-ntgrid:ntgrid)
    real :: tparfac(-ntgrid:ntgrid), tperpfac(-ntgrid:ntgrid)
    real :: ct(-ntgrid:ntgrid), st(-ntgrid:ntgrid)
    real :: c2t(-ntgrid:ntgrid), s2t(-ntgrid:ntgrid)

    ! normalization
    real :: fac
    real, allocatable :: nnrm(:,:,:),unrm(:,:,:)
    real, allocatable :: tparanrm(:,:,:),tperpnrm(:,:,:)

    real :: check(3)=0.,current=0.
    character (len=2) :: cfit='A0'
    complex, allocatable :: dens(:,:,:,:),upar(:,:,:,:)
    complex, allocatable :: tpar(:,:,:,:),tper(:,:,:,:)
    complex, allocatable :: phi(:,:,:),apar(:,:,:),bpar(:,:,:)
    complex, allocatable :: jpar(:,:,:)

    real :: save_dens0, save_tperp0, save_u0
    real :: ratio

    ! real space profile to be Fourier transformed
    real :: xx,dx,lx,ly
    integer, parameter :: nfxp=2**10
    integer :: nfx=nfxp
    real, allocatable :: nxy(:,:),uxy(:,:)
    real, allocatable :: phixy(:,:,:),aparxy(:,:,:),bparxy(:,:,:)
    real, allocatable :: jparxy(:,:,:)
    real, allocatable :: densxy(:,:,:,:),uparxy(:,:,:,:)
    real, allocatable :: tparxy(:,:,:,:),tperxy(:,:,:,:)

    complex, allocatable :: by(:,:,:)
    real, allocatable :: byxy(:,:,:)

    integer :: nnx,nny
    integer :: unit
    real :: xmax, bymax, xa, xl1, xl2
    real :: fmin, fmax
    
    real :: a,b

    if(debug.and.proc0) write(6,*) 'Initialization recon3'

    if (nfxp < nx) nfx=nx

!!! adjust input parameters to kill initial field if wanted
    if(debug.and.proc0) write(6,*) 'check parameters'

    do is=1,nspec
       if(adj_spec == is) cycle
       check(1)=check(1)+spec(is)%z*spec(is)%dens*spec(is)%dens0
       check(2)=check(2)+spec(is)%dens*spec(is)%temp* &
            & (spec(is)%dens0+spec(is)%tperp0)
       check(3)=check(3)+spec(is)%z*spec(is)%dens*spec(is)%stm*spec(is)%u0
    end do

    if(adj_spec == 0) then
       ! just warn if fields don't satisfy given conditions
       if(null_phi.and.null_bpar) then
          if(sum(check(1:2)) /= 0.) then
             if(proc0) write(6,'(a)') &
                  'WARNING: finite Phi or Bpar in initial condition'
          endif
       else if(null_bpar.and..not.null_phi) then
          ratio=gamtot1(0,eq_mode_n+1,1)/gamtot(0,eq_mode_n+1,1)
          if(check(1)/check(2) /= ratio) then
             if(proc0) write(6,'(a)') &
                  'WARNING: finite Bpar in initial condition'
          endif
       else if(null_phi.and..not.null_bpar) then
          ratio=-(2.*gamtot2(0,eq_mode_n+1,1)+2./beta) &
               & /gamtot1(0,eq_mode_n+1,1)
          if(check(1)/check(2) /= ratio) then
             if(proc0) write(6,'(a)') &
                  'WARNING: finite Bpar in initial condition'
          endif
       endif
       if(null_apar) then
          if(check(3) /= 0.) then
             if(proc0) write(6,'(a)') &
                  'WARNING: finite Apar in initial condition'
          endif
       endif
    else
       ! adjust input parameter to satisfy given conditions
       if(null_phi.and.null_bpar) then
          save_dens0=spec(adj_spec)%dens0
          save_tperp0=spec(adj_spec)%tperp0
          spec(adj_spec)%dens0=-check(1)/(spec(adj_spec)%z*spec(adj_spec)%dens)
          spec(adj_spec)%tperp0=-spec(adj_spec)%dens0 &
               & -check(2)/(spec(adj_spec)%dens*spec(adj_spec)%temp)

          if(spec(adj_spec)%dens0 /= save_dens0) then
             if(proc0) write(6,'(a,i0,a,f10.2,a)') &
                  & 'WARNING: Initial density of spec=', spec(adj_spec)%type, &
                  & ' is adjusted to ', spec(adj_spec)%dens0, &
                  & ' to kill Phi and Bpar'
          endif
          if(spec(adj_spec)%tperp0 /= save_tperp0) then
             if(proc0) write(6,'(a,i0,a,f10.2,a)') &
                  & 'WARNING: Initial Tperp of spec=', spec(adj_spec)%type, &
                  & ' is adjusted to ', spec(adj_spec)%tperp0, &
                  & ' to kill Phi and Bpar'
          endif
       else if(null_bpar.and..not.null_phi.and.eq_type.eq.'sinusoidal') then
          save_tperp0=spec(adj_spec)%tperp0
          check(1)=check(1)+ &
               & spec(adj_spec)%z*spec(adj_spec)%dens*spec(adj_spec)%dens0
          ratio=gamtot1(0,eq_mode_n+1,1)/gamtot(0,eq_mode_n+1,1)

          spec(adj_spec)%tperp0=(-check(1)*ratio-check(2)) &
               & /(spec(adj_spec)%dens*spec(adj_spec)%temp) &
               & -spec(adj_spec)%dens0
          if(spec(adj_spec)%tperp0 /= save_tperp0) then
             if(proc0) write(6,'(a,i0,a,f10.2,a)') &
                  & 'WARNING: Initial Tperp of spec=', spec(adj_spec)%type, &
                  & ' is adjusted to ', spec(adj_spec)%tperp0, &
                  & ' to kill Bpar'
          endif
       else if(null_phi.and..not.null_bpar.and.eq_type.eq.'sinusoidal') then
          save_tperp0=spec(adj_spec)%tperp0
          check(1)=check(1)+ &
               & spec(adj_spec)%z*spec(adj_spec)%dens*spec(adj_spec)%dens0
          ratio=-(2.*gamtot2(0,eq_mode_n+1,1)+2./beta) &
               & /gamtot1(0,eq_mode_n+1,1)

          spec(adj_spec)%tperp0=(-check(1)*ratio-check(2)) &
               & /(spec(adj_spec)%dens*spec(adj_spec)%temp) &
               & -spec(adj_spec)%dens0

          if(spec(adj_spec)%tperp0 /= save_tperp0) then
             if(proc0) write(6,'(a,i0,a,f10.2,a)') &
                  & 'WARNING: Initial Tperp of spec=', spec(adj_spec)%type, &
                  & ' is adjusted to ', spec(adj_spec)%tperp0, &
                  & ' to kill Phi'
          endif
       endif
    
       if (null_apar) then
          save_u0=spec(adj_spec)%u0
          spec(adj_spec)%u0=-check(3)/ &
               & (spec(adj_spec)%z*spec(adj_spec)%dens*spec(adj_spec)%stm)
          if(spec(adj_spec)%u0 /= save_u0) then
             if(proc0) write(6,'(a,i0,a,f10.2,a)') &
                  & 'WARNING: Initial U of spec=', spec(adj_spec)%type, &
                  & ' is adjusted to ', spec(adj_spec)%u0, &
                  & ' to kill Apar'
          endif
       endif
    endif

!!! initialize
    if(debug.and.proc0) write(6,*) 'initialize variable'
    nkxyz(-ntgrid:ntgrid,1:nakx,1:naky)=cmplx(0.,0.)
    ukxyz(-ntgrid:ntgrid,1:nakx,1:naky)=cmplx(0.,0.)
    nkxyz_eq(-ntgrid:ntgrid,1:nakx,1:naky)=cmplx(0.,0.)
    ukxyz_eq(-ntgrid:ntgrid,1:nakx,1:naky)=cmplx(0.,0.)

    nkxy_eq(1:nakx,1:naky)=cmplx(0.,0.)
    ukxy_eq(1:nakx,1:naky)=cmplx(0.,0.)

    dfac(-ntgrid:ntgrid)=1.
    ufac(-ntgrid:ntgrid)=