2021218321

Committed 04 Sep 2025 07:44AM UTC coverage: 10.606% (+0.03%) from 10.577%

Build # 2021218321

Build Type

push

gitlab-ci

Committed by David Dickinson

Commit Message

Merged in feature/move_more_initialisation_to_init_levels (pull request #1161)

Run Details

4710 of 44407 relevant lines covered (10.61%)

125698.1 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

1.65

/src/init_g.f90

!> 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, ginitopt_restart_many, restart_file, get_restart_dir, set_restart_dir
  public :: init_init_g, finish_init_g, wnml_init_g, check_init_g, tstart
  public :: new_field_init, initial_condition_is_nonadiabatic_dfn, add_restart_dir_to_file
  public :: init_g_config_type, set_init_g_config, 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_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_eig=35, ginitopt_no_zonal = 36, &
       ginitopt_stationary_mode = 37, ginitopt_cgyro = 38, ginitopt_random_sine = 39, &
       ginitopt_constant = 40

  real :: width0, dphiinit, phiinit, imfac, refac, zf_init, phifrac
  real :: init_spec_pow, k0
  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(3) :: 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.

  interface do_chop_side
     module procedure do_chop_side_slice, do_chop_side_full
  end interface do_chop_side
  
  !> 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(3) :: 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
     !> Initial spectrum power index. with
     !> [[init_g_knobs:ginit_option]] `"noise"`.
     real :: init_spec_pow = 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]     
     !> Used in [[init_g_knobs:ginit_option]] "harris" and "convect".
     real :: k0 = 1.0
     !> 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(3) :: 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
#include "init_g_overrides_and_bound_auto_gen.inc"
     procedure :: set_smart_defaults => set_smart_defaults_local
  end type init_g_config_type
  
  type(init_g_config_type) :: init_g_config
  
contains

  !> FIXME : Add documentation
  subroutine wnml_init_g(unit)
    implicit none
    integer, intent(in) :: unit
    call init_g_config%write(unit)
  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
    use warning_helpers, only: not_exactly_equal
    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 (not_exactly_equal(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_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

  !> Join restart_dir to restart_file, accounting for potential existence
  !> of path in restart_file already.
  subroutine add_restart_dir_to_file(restart_dir, restart_file)
    character(len = *), intent(in out) :: restart_dir, restart_file
    integer :: ind_slash
    ! 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 == 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:))
    end if
  end subroutine add_restart_dir_to_file

  !> FIXME : Add documentation
  subroutine init_init_g(init_g_config_in)
    use gs2_save, only: set_restart_file, proc_to_save_fields
    use mp, only: nproc
    implicit none
    type(init_g_config_type), intent(in), optional :: init_g_config_in    

    if (initialized) return
    initialized = .true.

    call read_parameters(init_g_config_in)
    call add_restart_dir_to_file(restart_dir, restart_file)
    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
    use optionals, only: get_option_with_default
    implicit none
    logical, intent (out) :: restarted
    integer, intent(in), optional :: override_ginitopt_switch
    logical, save :: have_warned_about_hybrid = .false.
    integer :: ginitopt_actual

    ginitopt_actual = get_option_with_default(override_ginitopt_switch, ginitopt_switch)

    restarted = .false.
    select case (ginitopt_actual)
    case (ginitopt_default)
       call ginit_default(use_odd = .false.)
    case (ginitopt_default_odd)
       call ginit_default(use_odd = .true.)
    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(include_apar = .false.)
    case (ginitopt_nl3r)
       call ginit_nl3(include_apar = .true.)
    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)
       ! this line:
       !       tstart = t0
       restarted = .true.
       scale = 1.
    case (ginitopt_nl5)
       call ginit_nl5
       restarted = .true.
       scale = 1.
    case (ginitopt_recon)
       call ginit_recon
       restarted = .true.
       scale = 1.
    case (ginitopt_nl6)
       call ginit_nl6
       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_many)
       call ginit_restart_many 
       call read_t0_from_restart_file (tstart)
       restarted = .true.
       scale = 1.
    case (ginitopt_restart_eig)
       call ginit_restart_eig 
       call read_t0_from_restart_file (tstart)
       restarted = .true.
       scale = 1.
    case (ginitopt_restart_small)
       call ginit_restart_small
       call read_t0_from_restart_file (tstart)
       restarted = .true.
       scale = 1.
    case (ginitopt_zonal_only)
       call ginit_restart_zonal_only
       call read_t0_from_restart_file (tstart)
       restarted = .true.
       scale = 1.
    case (ginitopt_no_zonal)
       call ginit_restart_no_zonal
       call read_t0_from_restart_file (tstart)
       restarted = .true.
       scale = 1.
    case (ginitopt_smallflat)
       call ginit_restart_smallflat
       call read_t0_from_restart_file (tstart)
       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

  !> Set the smart defaults for the init_g_config_type instance
  subroutine set_smart_defaults_local(self)
    use file_utils, only: run_name
    implicit none
    class(init_g_config_type), intent(in out) :: self
    if (self%is_initialised()) return
    if (.not. self%override_restart_file) self%restart_file = trim(run_name)//".nc"
  end subroutine set_smart_defaults_local

  !> FIXME : Add documentation
  subroutine read_parameters(init_g_config_in)
    use file_utils, only: error_unit
    use text_options, only: text_option, get_option_value
    use gs2_save, only: read_many, include_explicit_source_in_restart, proc_to_save_fields
    use mp, only: job
    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_many), &
            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

    if (present(init_g_config_in)) init_g_config = init_g_config_in
    call init_g_config%init(name = 'init_g_knobs', requires_index = .false.)

    ! Copy out internal values into module level parameters
    associate(self => init_g_config)
#include "init_g_copy_out_auto_gen.inc"
    end associate

    ! MAB - allows for ensemble averaging of multiple flux tube calculations
    ! job=0 if not doing multiple flux tube calculations, so phiinit unaffected
    phiinit = phiinit * (job * phifrac + 1.0)

    call get_option_value &
         (ginit_option, ginitopts, ginitopt_switch, &
         error_unit(), "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_odd)
    use species, only: spec
    use theta_grid, only: ntgrid, theta
    use kt_grids, only: naky, ntheta0, theta0, aky, enforce_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
    use warning_helpers, only: is_zero
    implicit none
    logical, intent(in) :: use_odd
    complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi
    integer :: iglo, ig, ik, it, il, is

    do ik = 1, naky
       do it = 1, ntheta0
          do ig = -ntgrid, ntgrid
             phi(ig, it, ik) = exp(-((theta(ig)-theta0(it, ik))/width0)**2) * cmplx(1.0,1.0)
          end do
       end do
    end do
    if (use_odd) then
       do ik = 1, naky
          do it = 1, ntheta0
             do ig = -ntgrid, ntgrid
                phi(ig, it, ik) = phi(ig, it, ik) * sin((theta(ig)-theta0(it, ik))/width0)
             end do
          end do
       end do
    end if

    call do_chop_side(phi)
    
    if (enforce_reality) then
       phi(:, 1, 1) = 0.0
       if (naky > 1 .and. is_zero(aky(1))) phi(:, :, 1) = 0.0
       call apply_reality(phi)
    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

  !> 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, enforce_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
    use warning_helpers, only: is_zero
    implicit none
    complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi
    integer :: iglo, ik, it, il, is

    phi = cmplx(refac,imfac)
    call do_chop_side(phi)
    
    if (enforce_reality) then
       phi(:,1,1) = 0.0
       if (naky > 1 .and. is_zero(aky(1))) phi(:,:,1) = 0.0
       call apply_reality(phi)
    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

  !> 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: l_links, r_links, naky, ntheta0, aky, enforce_reality, akx, kperp2
    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: pass_right, init_pass_ends, has_linked_boundary
    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
    use constants, only: twopi
    use warning_helpers, only: is_zero
    implicit none
    complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi, phit
    real :: a, b, spec_pow
    integer :: iglo, ig, ik, it, il, is, nn
    logical :: linked_bc

    linked_bc = has_linked_boundary()

    !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
             if (init_spec_pow > 0.) then
                spec_pow = 1./sqrt(1.+kperp2(ig,it,ik)**(0.5*init_spec_pow))
                phi(ig,it,ik) = spec_pow*exp(cmplx(0.,twopi*a))
             else
                phi(ig,it,ik) = cmplx(a,b)
             end if
          end do
          !CMR,28/1/13:
          ! clean_init debrutalises influence of chop_side on phi with linked bc
          if (chop_side) then
             if (clean_init .and. linked_bc) then
                if (left) 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) < 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) == r_links(it, ik)) then
                      phi(:-1,it,ik) = 0.0
                   endif
                else
                   !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) < 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) == r_links(it, ik)) then
                      phi(1:,it,ik) = 0.0
                   endif
                end if
             else
                call do_chop_side(phi(:, it, ik))
             end if
          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 >= 1 .and. is_zero(aky(1))) then
       !Apply scaling factor
       phi(:,:,1) = phi(:,:,1)*zf_init
       
       !Set ky=kx=0.0 mode to zero in amplitude
       if (is_zero(akx(1))) phi(:,1,1) = 0.0
       
       !CMR, 25/1/13: clean_init forces periodic zonal flows
       if (clean_init .and. linked_bc) then
          phi(ntgrid,:,1)=phi(-ntgrid,:,1)
          phit(ntgrid,:,1)=phit(-ntgrid,:,1)
       end if
    end if

    if (enforce_reality) then
       call apply_reality(phi)
       call apply_reality(phit)
    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. linked_bc .and. ik > 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) == 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) == 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. linked_bc .and. sum(r_links+l_links) > 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, shat
    use kt_grids, only: naky, ntheta0, aky, enforce_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
    use warning_helpers, only: is_zero
    implicit none
    complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi, phit
    real :: a, b
    integer :: iglo, 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() .or. is_zero(shat)) 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
             a = cos(kpar_init * theta(ig))
             b = sin(kpar_init * theta(ig))
             phi(ig,it,ik) = cmplx(a,b)
          end do
       end do
    end do

    call do_chop_side(phi)
    if (naky > 1 .and. is_zero(aky(1))) phi(:,:,1) = phi(:,:,1) * zf_init
    if (ikx_init  > 0) call single_initial_kx(phi)

    if (enforce_reality) then
       call apply_reality(phi)
       call apply_reality(phit)
    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, enforce_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
    use warning_helpers, only: is_zero
    implicit none
    complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi, phit
    real :: a, b
    integer :: iglo, 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
             phi(ig,it,ik) = cmplx(a,b)
          end do
       end do
    end do

    call do_chop_side(phi)
    if (ikx_init  > 0) call single_initial_kx(phi)

    if (enforce_reality) then
       call apply_reality(phi)
       call apply_reality(phit)
    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, 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
       call do_chop_side(phi(:, it, ik))
    end do

    call apply_reality(phi)

    ! 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, 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
       call do_chop_side(phi(:, it, ik))
    end do

    call apply_reality(phi)
    
    ! 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
       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(include_apar)
    use theta_grid, only: ntgrid, theta
    use kt_grids, only: naky, ntheta0, theta0, enforce_reality
    use dist_fn_arrays, only: g, gnew
    use array_utils, only: copy, zero_array
    use fields_arrays, only: apar, aparnew
    implicit none
    logical, intent(in) :: include_apar
    complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi
    integer :: ig, ik, it, j
    call zero_array(phi)
    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
          if (include_apar) apar(:, it, ik) = apar0 * phi(:, it, ik)
       end if
       call do_chop_side(phi(:, it, ik))
    end do

    if (enforce_reality .and. include_apar) call apply_reality(apar)
    if (include_apar) call copy(apar, aparnew)

    call set_g_from_moment_facs(g, phiinit, phi, dfac_species_weight = .true., &
         ufac_species_weight = .true., use_even = even, do_reality = enforce_reality)
    call copy(g, gnew)
  end subroutine ginit_nl3

  !> FIXME : Add documentation  
  subroutine ginit_harris
    use species, only: spec, has_electron_species
    use kt_grids, only: naky, ntheta0, enforce_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 array_utils, only: copy, zero_array
    implicit none
    complex, dimension (1,ntheta0,naky) :: phi
    integer :: iglo, 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(1,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(1,j,1) = ff_pr(nx-ntheta0+j)
    end do

    if (enforce_reality) call apply_reality(phi)
    
    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(1,it,ik)/aj0(:,iglo)
          g(:,2,iglo) =  2.* vpa(:,2,iglo)*spec(is)%u0 * phi(1,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 theta_grid, only: ntgrid
    use kt_grids, only: naky, ntheta0, enforce_reality
    use gs2_save, only: gs2_restore
    use dist_fn_arrays, only: g, gnew
    use fields_arrays, only: phi, apar, bpar, phinew, aparnew, bparnew
    use gs2_layouts, only: g_lo, ik_idx
    use run_parameters, only: has_phi, has_apar, has_bpar
    use array_utils, only: copy, zero_array
    implicit none
    complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phiz
    integer :: iglo, ig, ik, it, j
    
    call gs2_restore (g, scale, has_phi, has_apar, has_bpar)

    do iglo = g_lo%llim_proc, g_lo%ulim_proc
       ik = ik_idx(g_lo,iglo)
       if (ik == 1) cycle
       g (:, :, 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)
    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

    call set_g_from_moment_facs(g, phiinit, phiz, dfac_species_weight = .true., &
         ufac_species_weight = .false., use_even = even, do_reality = enforce_reality)
    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 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 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, ig, ik, it, is, il
    
    call gs2_restore (g, scale, has_phi, has_apar, has_bpar)

    do iglo = g_lo%llim_proc, g_lo%ulim_proc
       it = it_idx(g_lo,iglo)
       ik = ik_idx(g_lo,iglo)
       if (it == 1 .and. ik == 2) cycle
       g (:,:,iglo) = 0.
    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
       end do
    end do

    call do_chop_side(phiz)
    call apply_reality(phiz)

    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
       where (forbid(:,il)) g(:,1,iglo) = 0.0
       g(:, 2, iglo) = g(:, 1, iglo)
    end do
    call copy(g, gnew)
  end subroutine ginit_nl4

  !> FIXME : Add documentation. This is near identical to ginit_nl4 aside from initial zeroing
  subroutine ginit_nl5
    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 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, ig, ik, it, is, il
    call gs2_restore (g, scale, has_phi, has_apar, has_bpar)

    do iglo = g_lo%llim_proc, g_lo%ulim_proc
       ik = ik_idx(g_lo, iglo)
       if (ik == 1) cycle
       g (:, :, 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
       end do
    end do

    call do_chop_side(phiz)
    call apply_reality(phiz)

    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
       where (forbid(:,il)) g(:,1,iglo) = 0.0
       g(:, 2, iglo) = g(:, 1, iglo)
    end do
    call copy(g, gnew)
  end subroutine ginit_nl5

  !> FIXME : Add documentation  
  subroutine ginit_nl6
    use gs2_save, only: gs2_restore
    use kt_grids, only: ntheta0
    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 run_parameters, only: has_phi, has_apar, has_bpar
    use array_utils, only: copy, zero_array
    implicit none
    integer :: iglo, ik, it
    
    call gs2_restore (g, scale, has_phi, has_apar, has_bpar)
    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) = 0.
    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 theta_grid, only: ntgrid
    use kt_grids, only: naky, ntheta0, enforce_reality
    use dist_fn_arrays, only: g, gnew
    use array_utils, only: copy
    implicit none
    complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi
    integer :: ig, ik, it, j
    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

    call set_g_from_moment_facs(g, phiinit, phi, dfac_species_weight = .true., &
         ufac_species_weight = .false., use_even = even, do_reality = enforce_reality)
    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, enforce_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
    complex :: 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

    if (enforce_reality) then
       call apply_reality(phi)
       call apply_reality(jpar)
    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
          ! Should this include abs?
          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

    if (enforce_reality) then
       call apply_reality(phi)
       call apply_reality(jpar)
    end if

    ! redefine g
    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 dist_fn_arrays, only: g, gnew
    use array_utils, only: copy
    implicit none
    complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi
    integer :: ig
    
    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

    call do_chop_side(phi)

    call set_g_from_moment_facs(g, phiinit, phi, dfac_species_weight = .false., &
         ufac_species_weight = .false., use_even= .true., do_reality = .false.)

    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
    implicit none
    complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi, odd
    integer :: iglo, ik, it, il, is
    real :: phase
    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
  
    !Should this be triggered for kt_grids::reality=.true. only?
    call apply_reality(phi)
    call apply_reality(odd)

    ! charge dependence keeps initial Phi from being too small
    ! Note not quite suitable for set_g_from_moment_facs as dfac --> den1 etc.
    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, 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, 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, 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, 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
    use warning_helpers, only: is_zero
    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(is_zero(vpa(-ntgrid:ntgrid-1,:,iglo)))) 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 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, enforce_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: gamtot, gamtot1, gamtot2, get_init_field
    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
    use warning_helpers, only: not_exactly_equal, is_not_zero, is_zero
    implicit none
    integer, parameter :: ncnt_mom_coeff=8
    real, dimension(:, :, :, :), allocatable :: mom_coeff
    real, dimension(:, :, :), allocatable :: mom_shift_para, mom_shift_perp
    real :: bsq
    integer :: iglo, ig, ik, it, is, 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), current
    character (len=2) :: cfit
    complex, allocatable :: phi(:,:,:), apar(:,:,:), bpar(:,:,:)
    real :: save_dens0, save_tperp0, save_u0, ratio

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

    complex, allocatable :: by(:,:,:)
    real, allocatable :: byxy(:,:,:)
    integer :: nnx, nny, unit
    real :: xmax, bymax, xa, xl1, xl2, fmin, fmax, a, b

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

    cfit = 'A0'
    check = 0
    current = 0
    nfx = nfxp
    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 (is_not_zero(sum(check(1:2)))) 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 (not_exactly_equal(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 (not_exactly_equal(check(1)/check(2), ratio)) then
             if (proc0) write(6,'(a)') 'WARNING: finite Bpar in initial condition'
          endif
       endif
       if (null_apar) then
          if (is_not_zero(check(3))) then
             if (proc0) write(6,'(a)') 'WARNING: finite Apar in initial condition'
          end if
       end if
    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 (not_exactly_equal(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 (not_exactly_equal(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=='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 (not_exactly_equal(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=='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 (not_exactly_equal(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 (not_exactly_equal(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)=1.
    tparfac(-ntgrid:ntgrid)=1.
    tperpfac(-ntgrid:ntgrid)=1.

    !! equilibrium
    if (is_not_zero(phiinit0)) then
       if (nakx==1 .and. naky==1) then ! grid_option = single case
          nkxy_eq(1,1)=cmplx(.5,0.)
          ukxy_eq(1,1)=cmplx(.5,0.)
       else
          if (debug .and. proc0) write(6,*) 'set equilibrium profile'
          allocate(nxy(nfx,ny),uxy(nfx,ny))
          lx=2.*pi/(akx(2)-akx(1)); ly=2.*pi/(aky(2)-aky(1))
          dx=lx/nfx
          ! if width is negative, it gives the ratio to the box size
          if(prof_width < 0. ) prof_width=-prof_width*lx
          select case (eq_type)
          case ('sinusoidal')
             ! this gives n,Tpara,Tperp \propto cos^2 (2 pi/Lx)
             ! nxy(eq_mode1(1),eq_mode1(2))=cmplx(.25, 0.)
             ! this gives Apara \propto cos(2 pi/Lx), By \propto sin(2 pi x/Lx)
             ! uxy(eq_mode2(1),eq_mode2(2))=cmplx(.5, 0.)
             do it=1,nfx
                xx=dx*(it-1)
                nxy(it,1:ny)=0.
                uxy(it,1:ny)=cos(2.*pi/lx*xx*eq_mode_u)
             end do
          case ('porcelli')
             do it=1,nfx
                xx=dx*(it-1)
                nxy(it,1:ny)=0.
                uxy(it,1:ny)=1./cosh((xx-.5*lx)/prof_width)**2 &
                     & * (tanh(2.*pi/lx*xx)**2+tanh(2.*pi/lx*(xx-lx))**2 &
                     & - tanh(2.*pi)**2) / (2.*tanh(pi)**2-tanh(2.*pi)**2)
             end do
          case ('doubleharris')
             do it=1,nfx
                fmin=tanh(.75*lx/prof_width)-tanh(.25*lx/prof_width)
                fmax=2.*tanh(.25*lx/prof_width)
                xx=dx*(it-1)
                nxy(it,1:ny)=0.
                uxy(it,1:ny)= 2.*( &
                     & tanh((xx-.25*lx)/prof_width) - &
                     & tanh((xx-.75*lx)/prof_width) - fmin)/(fmax-fmin)-1.
             end do
          case default
             if(proc0) write(error_unit(),'(2a)') &
                  & 'ERROR: Invalid equilibrium type', eq_type
             call mp_abort('Invalid equilibrium type')
          end select
          ! subtract (0,0) mode
          ! since it (constant part of potential) does nothing
          do ik=1,ny
             uxy(1:nfx,ik)=uxy(1:nfx,ik) &
                  - sum(uxy(1:nfx,ik))/nfx
          end do
          call inverse2(nxy, nkxy_eq, ny, nfx)
          call inverse2(uxy, ukxy_eq, ny, nfx)
          deallocate(nxy,uxy)
       endif
    endif

    do ig = -ntgrid, ntgrid
       nkxyz(ig,1:nakx,1:naky) = phiinit0*nkxy_eq(1:nakx,1:naky)
       ukxyz(ig,1:nakx,1:naky) = phiinit0*ukxy_eq(1:nakx,1:naky)
    end do
    if(.not.(nakx==1 .and. naky==1)) then ! grid_option /= single case
       nkxyz(-ntgrid:ntgrid,1,1)=cmplx(0.,0.)
       ukxyz(-ntgrid:ntgrid,1,1)=cmplx(0.,0.)
    endif

    !! save equilibrium profile
    nkxyz_eq(-ntgrid:ntgrid,1:nakx,1:naky)= &
         & nkxyz(-ntgrid:ntgrid,1:nakx,1:naky)
    ukxyz_eq(-ntgrid:ntgrid,1:nakx,1:naky)= &
         & ukxyz(-ntgrid:ntgrid,1:nakx,1:naky)

    !! perturbation
    if (is_not_zero(phiinit)) then
       if (debug.and.proc0) write(6,*) 'set perturbation profile'
       do j = 1, 3
          if(ikkk(j)==0) ukxy_pt(j)=.5*ukxy_pt(j) !reality
          ik = ikkk(j)+1
          if(ittt(j) >= 0) then
             it = ittt(j)+1
          else
             it = nakx + ittt(j) + 1
          endif
          do ig = -ntgrid, ntgrid
             ukxyz(ig,it,ik) = ukxyz(ig,it,ik) + phiinit*ukxy_pt(j)
          end do
       end do
    endif

    if (is_not_zero(phiinit_rand)) then
       do it = 1, nakx
          do ik = 1, naky
             a = ranf()-0.5
             b = ranf()-0.5
             ukxyz(:,it,ik) = ukxyz(:,it,ik) + phiinit_rand*cmplx(a,b)
          end do
       end do
    end if

    if (debug.and.proc0) write(6,*) 'set reality condition'
    if (enforce_reality) then
       call apply_reality(nkxyz)
       call apply_reality(ukxyz)
       call apply_reality(nkxyz_eq)
       call apply_reality(ukxyz_eq)
    end if

    !! parallel profile
    if(debug.and.proc0) write(6,*) 'set parallel profile'
    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     = dfac     + den1   * ct + den2   * c2t
    ufac     = ufac     + upar1  * st + upar2  * s2t
    tparfac  = tparfac  + tpar1  * ct + tpar2  * c2t
    tperpfac = tperpfac + tperp1 * ct + tperp2 * c2t

    !! normalization
    if (debug .and. proc0) write(6,*) 'normalization'
    if (is_not_zero(b0)) then
       if (.not.(nakx==1 .and. naky==1)) then ! grid_option /= single case
          if (eq_type == 'porcelli') then
             xmax=.5*prof_width*log(2.+sqrt(3.))+.5*lx
             xmax=get_xmax(xmax)
             xa=(xmax-.5*lx)/prof_width
             xl1=2.*pi/lx*xmax; xl2=xl1-2.*pi
             bymax=(-2./prof_width*sinh(xa)/cosh(xa)**3* &
                  & (tanh(xl1)**2+tanh(xl2)**2-tanh(2.*pi)**2) &
                  & + 1./cosh(xa)**2*4.*pi/lx* &
                  & (sinh(xl1)/cosh(xl1)**3+sinh(xl2)/cosh(xl2)**3) ) &
                  & / (2.*tanh(pi)**2-tanh(2.*pi)**2)
          else if(eq_type == 'doubleharris') then
             bymax=2.*tanh(.5*lx/prof_width)**2/(prof_width*(fmax-fmin))
          else
             bymax=akx(eq_mode_u+1)
          endif
       else
          bymax=sqrt(kperp2(0,1,1))
       endif
       a0 = b0/abs(bymax)
       cfit='B0'
    endif

    allocate(nnrm(nakx,naky,nspec),unrm(nakx,naky,nspec))
    allocate(tparanrm(nakx,naky,nspec),tperpnrm(nakx,naky,nspec))

    if(.not.allocated(mom_coeff)) allocate(mom_coeff(nakx,naky,nspec,ncnt_mom_coeff))
    if(.not.allocated(mom_shift_para)) allocate(mom_shift_para(nakx,naky,nspec))
    if(.not.allocated(mom_shift_perp)) allocate(mom_shift_perp(nakx,naky,nspec))

    do it=1,nakx
       do ik=1,naky
          do is=1,nspec
             bsq=.25*spec(is)%smz**2*kperp2(0,it,ik)
             mom_coeff(it,ik,is,1) = exp(-bsq)
             mom_coeff(it,ik,is,2) = exp(-bsq) *.5
             mom_coeff(it,ik,is,3) = exp(-bsq) *(1.-bsq)
             mom_coeff(it,ik,is,4) = exp(-bsq) *.75
             mom_coeff(it,ik,is,5) = exp(-bsq) *(1.-bsq)*.5
             mom_coeff(it,ik,is,6) = exp(-bsq) *.5
             mom_coeff(it,ik,is,7) = exp(-bsq) *.25
             mom_coeff(it,ik,is,8) = exp(-bsq) *(1.-.5*bsq)
          end do
       end do
    end do
    mom_shift_para=mom_coeff(:,:,:,2)/mom_coeff(:,:,:,1)
    mom_shift_perp=mom_coeff(:,:,:,3)/mom_coeff(:,:,:,1)

    do is=1,nspec
       nnrm(:,:,is)=mom_coeff(:,:,is,1)
       unrm(:,:,is)=2.*mom_coeff(:,:,is,2)
       tparanrm(:,:,is)=2.*(mom_coeff(:,:,is,4)- &
            & 2.*mom_coeff(:,:,is,2)*mom_shift_para(:,:,is))
       tperpnrm(:,:,is)=2.*(mom_coeff(:,:,is,8)- &
            & mom_coeff(:,:,is,6)*mom_shift_perp(:,:,is))

       if (is_not_zero(a0)) then
          ! rescale parallel momentum to obtain a given amplitude of Apar
          current=sum(spec(1:nspec)%z*spec(1:nspec)%dens &
               & *spec(1:nspec)%stm*spec(1:nspec)%u0)
          if (is_zero(current)) then
             if (proc0) write(error_unit(),'(a)') &
                  & 'ERROR in init_g: invalid input a0, u0'
             call mp_abort('invalid input a0, u0')
          end if
          do it=1,nakx
             do ik=1,naky
                if(is_not_zero(cabs(ukxyz(0,it,ik))) .and. is_not_zero(spec(is)%u0)) then
                   fac=2.*beta*current/kperp2(0,it,ik)/a0
                   unrm(it,ik,is)=unrm(it,ik,is)*fac
                   if(proc0) write(6,'(a,a2,a,3(i3,1x),f20.12)') &
                        & 'INFO: u0 is rescaled to fit ',cfit,' input', &
                        & it,ik,is,spec(is)%u0/fac
                end if
             end do
          end do
       end if
    end do

    if (debug .and. proc0) write(6,*) 'calculate g'
    allocate(g_eq(-ntgrid:ntgrid,2,g_lo%llim_proc:g_lo%ulim_alloc))
    do iglo = g_lo%llim_proc, g_lo%ulim_proc
       ik = ik_idx(g_lo,iglo)
       it = it_idx(g_lo,iglo)
       is = is_idx(g_lo,iglo)

       g(:,1,iglo) = ( &
            & ( dfac*spec(is)%dens0 / nnrm(it,ik,is) &
            & + tparfac*spec(is)%tpar0 / tparanrm(it,ik,is) &
            & * (vpa(:,1,iglo)**2-mom_shift_para(it,ik,is)) &
            & + tperpfac*spec(is)%tperp0 / tperpnrm(it,ik,is) &
            & * (vperp2(:,iglo)-mom_shift_perp(it,ik,is)) &
            & ) * nkxyz(:,it,ik) &
            & + ufac*2.*vpa(:,1,iglo)*spec(is)%u0 / unrm(it,ik,is) &
            & * ukxyz(:,it,ik) &
            )
       g(:,2,iglo) = ( &
            & ( dfac*spec(is)%dens0 / nnrm(it,ik,is) &
            & + tparfac*spec(is)%tpar0 / tparanrm(it,ik,is) &
            & * (vpa(:,2,iglo)**2-mom_shift_para(it,ik,is)) &
            & + tperpfac*spec(is)%tperp0 / tperpnrm(it,ik,is) &
            & * (vperp2(:,iglo)-mom_shift_perp(it,ik,is)) &
            & ) * nkxyz(:,it,ik) &
            & + ufac*2.*vpa(:,2,iglo)*spec(is)%u0 / unrm(it,ik,is) &
            & * ukxyz(:,it,ik) &
            )

       g_eq(:,1,iglo) = ( &
            & ( dfac*spec(is)%dens0 / nnrm(it,ik,is) &
            & + tparfac*spec(is)%tpar0 / tparanrm(it,ik,is) &
            & * (vpa(:,1,iglo)**2-mom_shift_para(it,ik,is)) &
            & + tperpfac*spec(is)%tperp0 / tperpnrm(it,ik,is) &
            & * (vperp2(:,iglo)-mom_shift_perp(it,ik,is)) &
            & ) * nkxyz_eq(:,it,ik) &
            & + ufac*2.*vpa(:,1,iglo)*spec(is)%u0 / unrm(it,ik,is) &
            & * ukxyz_eq(:,it,ik) &
            )
       g_eq(:,2,iglo) = ( &
            & ( dfac*spec(is)%dens0 / nnrm(it,ik,is) &
            & + tparfac*spec(is)%tpar0 / tparanrm(it,ik,is) &
            & * (vpa(:,2,iglo)**2-mom_shift_para(it,ik,is)) &
            & + tperpfac*spec(is)%tperp0 / tperpnrm(it,ik,is) &
            & * (vperp2(:,iglo)-mom_shift_perp(it,ik,is)) &
            & ) * nkxyz_eq(:,it,ik) &
            & + ufac*2.*vpa(:,2,iglo)*spec(is)%u0 / unrm(it,ik,is) &
            & * ukxyz_eq(:,it,ik) &
            )
    end do

    deallocate(nnrm,unrm,tparanrm,tperpnrm)
    deallocate(g_eq)

    if(allocated(mom_coeff)) deallocate(mom_coeff)
    if(allocated(mom_shift_para)) deallocate(mom_shift_para)
    if(allocated(mom_shift_perp)) deallocate(mom_shift_perp)

    ! now store the actual g in gnew
    call copy(g, gnew)

    !! check
    if(debug.and.proc0) write(6,*) 'check'
    if(input_check_recon) then

       allocate(phi(-ntgrid:ntgrid,nakx,naky))
       allocate(apar(-ntgrid:ntgrid,nakx,naky))
       allocate(bpar(-ntgrid:ntgrid,nakx,naky))

       call get_init_field(phi,apar,bpar)

       if(nakx==1 .and. naky==1) then ! grid_option = single case
          if(proc0) then
             write(6,'(" phi  = ",2e25.17)') &
                  & real(phi (0,1,1)),aimag(phi (0,1,1))
             write(6,'(" apar = ",2e25.17)') &
                  & real(apar(0,1,1)),aimag(apar(0,1,1))
             write(6,'(" bpar = ",2e25.17)') &
                  & real(bpar(0,1,1)),aimag(bpar(0,1,1))
          endif
       else
          if(nx>nakx) then
             nnx=nx
          else ! nx=nakx=1
             nnx=(3*nakx+1)/2
          endif
          if(ny>naky) then
             nny=ny
          else ! ny=naky=1
             nny=3*naky
          endif
          allocate(phixy(nnx,nny,-ntgrid:ntgrid))
          allocate(aparxy(nnx,nny,-ntgrid:ntgrid))
          allocate(bparxy(nnx,nny,-ntgrid:ntgrid))

          allocate(by(-ntgrid:ntgrid,nakx,naky))
          allocate(byxy(nnx,nny,-ntgrid:ntgrid))

          call transform2(phi,phixy,nny,nnx)
          call transform2(apar,aparxy,nny,nnx)
          call transform2(bpar,bparxy,nny,nnx)
          do it=1,nakx
             by(:,it,:)=-zi*akx(it)*apar(:,it,:)
          end do
          call transform2(by,byxy,nny,nnx)

          if(proc0) then
             write(6,'(" phi : mode=",2i4," amp=",e25.17)') &
                  & eq_mode_n,0, &
                  & real(phi (0,eq_mode_n+1,1))
             write(6,'(" apar: mode=",2i4," amp=",e25.17)') &
                  & eq_mode_u,0, &
                  & real(apar(0,eq_mode_u+1,1))
             write(6,'(" bpar: mode=",2i4," amp=",e25.17)') &
                  & eq_mode_n,0, &
                  & real(bpar(0,eq_mode_n+1,1))

             call open_output_file(unit, '.check_field.dat')
             write(unit,'(a1,1x,5a20)') '#','x','y','phi','apar','bpar','by'
             do it=1,nx+1
                i=mod(it-1,nx)+1
                do ik=1,ny+1
                   j=mod(ik-1,ny)+1
                   write(unit,'(2x,6d20.12)') lx/nx*(it-1),ly/ny*(ik-1), &
                        & phixy(i,j,0),aparxy(i,j,0),bparxy(i,j,0),byxy(i,j,0)
                end do
                write(unit,*)
             end do
             call close_output_file(unit)
          endif
       endif
    endif

    deallocate(phi,apar,bpar)
    if(allocated(phixy)) deallocate(phixy)
    if(allocated(aparxy)) deallocate(aparxy)
    if(allocated(bparxy)) deallocate(bparxy)
    if(allocated(by)) deallocate(by)
    if(allocated(byxy)) deallocate(byxy)

  contains
    !> FIXME : Add documentation    
    function get_xmax(xguess) result(xsol)
      implicit none
      real :: xsol
      real, intent(in) :: xguess
      real, parameter :: tol=1.e-20
      real :: xold
      integer :: i
      integer, parameter :: itmax = 1000
      xold=xguess
      do i=1,itmax
         xsol=xold-f(xold)/fprime(xold)
         if(abs(xsol-xold) > tol) then
            xold=xsol
         else
            return
         endif
      end do
    end function get_xmax

    !> FIXME : Add documentation    
    function f(x)
      implicit none
      real :: f
      real, intent(in) :: x
      real :: xa, xl1, xl2
      xa=(x-.5*lx)/prof_width
      xl1=2.*pi/lx*x
      xl2=xl1-2.*pi

      f= &
           & 2./prof_width**2*(cosh(2.*xa)-2.)/cosh(xa)**4* & ! f''
           & (tanh(xl1)**2+tanh(xl2)**2-tanh(2.*pi)**2) & ! g
           & + 2. * &
           & (-2./prof_width)*sinh(xa)/cosh(xa)**3* & ! f'
           & 4.*pi/lx*(sinh(xl1)/cosh(xl1)**3+sinh(xl2)/cosh(xl2)**3) & ! g'
           & + &
           & 1./cosh(xa)**2* & ! f
           & 8.*(pi/lx)**2*((2.-cosh(2.*xl1))/cosh(xl1)**4 & ! g''
           &               +(2.-cosh(2.*xl2))/cosh(xl2)**4)

    end function f

    !> FIXME : Add documentation
    function fprime(x)
      implicit none
      real :: fprime
      real, intent(in) :: x
      real :: xa, xl1, xl2
      xa=(x-.5*lx)/prof_width
      xl1=2.*pi/lx*x
      xl2=xl1-2.*pi

      fprime = &
           & 8./prof_width**3*(2.*sinh(xa)-sinh(xa)**3)/cosh(xa)**5* & ! f''' 
           & (tanh(xl1)**2+tanh(xl2)**2-tanh(2.*pi)**2) & ! g
           & + 3. * &
           & 2./prof_width**2*(cosh(2.*xa)-2.)/cosh(xa)**4* & ! f''
           & 4.*pi/lx*(sinh(xl1)/cosh(xl1)**3+sinh(xl2)/cosh(xl2)**3) & ! g'
           & + 3. * &
           & (-2./prof_width)*sinh(xa)/cosh(xa)**3* & ! f'
           & 8.*(pi/lx)**2*((2.-cosh(2.*xl1))/cosh(xl1)**4 & ! g''
           &               +(2.-cosh(2.*xl2))/cosh(xl2)**4) &
           & + &
           & 1./cosh(xa)**2* & ! f
           & (-64.)*(pi/lx)**3 * ( & ! g'''
           & (2.*sinh(xl1)-sinh(xl1)**3)/cosh(xl1)**5 + &
           & (2.*sinh(xl2)-sinh(xl2)**3)/cosh(xl2)**5 )
    end function fprime

  end subroutine ginit_recon3

  !> Restore `g` from restart files. Note that the timestep and simulation time
  !> are restored by [[read_t0_from_restart_file]]
  subroutine ginit_restart_many
    use dist_fn_arrays, only: g, gnew
    use fields_arrays, only: phi, apar, bpar, phinew, aparnew, bparnew
    use gs2_save, only: gs2_restore
    use run_parameters, only: has_phi, has_apar, has_bpar
    use array_utils, only: copy, zero_array
    implicit none
    call gs2_restore (g, scale, has_phi, has_apar, has_bpar)
    call copy(g, gnew)
    call copy(phinew, phi) ; call copy(aparnew, apar) ; call copy(bparnew, bpar)
  end subroutine ginit_restart_many

  !> Uses the restart files written by the eigensolver. File to read is set by
  !> [[init_g_knobs:restart_eig_id]]
  subroutine ginit_restart_eig
    use dist_fn_arrays, only: g, gnew
    use fields_arrays, only: phi, apar, bpar, phinew, aparnew, bparnew
    use gs2_save, only: gs2_restore
    use run_parameters, only: has_phi, has_apar, has_bpar
    use array_utils, only: copy, zero_array
    implicit none
    character(len=20) :: restart_unique_string
    write(restart_unique_string,'(".eig_",I0)') restart_eig_id
    call gs2_restore (g, scale, has_phi, has_apar, has_bpar, fileopt=restart_unique_string)
    call copy(g, gnew)
    call copy(phinew, phi) ; call copy(aparnew, apar) ; call copy(bparnew, bpar)
  end subroutine ginit_restart_eig

  !> Restores `g` from restart files _on top of_ a background from [[ginit_noise]].
  !>
  !> Note that this doesnt set `phi`, `apar`, `bpar` -- is this correct?
  subroutine ginit_restart_small
    use dist_fn_arrays, only: g, gnew
    use gs2_save, only: gs2_restore
    use run_parameters, only: has_phi, has_apar, has_bpar
    use array_utils, only: copy, zero_array
    implicit none
    call ginit_noise
    call gs2_restore (g, scale, has_phi, has_apar, has_bpar)
    g = g + gnew
    call copy(g, gnew)
  end subroutine ginit_restart_small

  !> Similar to `"small"`, restores fields from restart files _on top of_ a
  !> random background.
  !>
  !> @Todo this starts with a reduced functionality version of ginit_noise
  !> followed by ginit_restart_many (effectively). This is therefore a reduced
  !> functionality version of ginit_restart_small and we can probably just drop this?
  subroutine ginit_restart_smallflat
    use gs2_save, only: gs2_restore
    use species, only: spec
    use theta_grid, only: ntgrid 
    use kt_grids, only: naky, ntheta0, aky, enforce_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 run_parameters, only: has_phi, has_apar, has_bpar
    use ran, only: ranf
    use array_utils, only: copy, zero_array
    use warning_helpers, only: is_zero
    implicit none
    complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi
    real :: a, b
    integer :: iglo, ik, it, il, is

    do it = 1, ntheta0
       do ik = 1, naky
          a = ranf()-0.5
          b = ranf()-0.5
          phi(:,it,ik) = cmplx(a,b)
       end do
    end do

    if (naky > 1 .and. is_zero(aky(1))) phi(:,:,1) = phi(:,:,1) * zf_init
    if (enforce_reality) call apply_reality(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)
       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)

    call gs2_restore (g, scale, has_phi, has_apar, has_bpar)
    g = g + gnew
    call copy(g, gnew)
  end subroutine ginit_restart_smallflat

  !> Restart but set all non-zonal components \((k_y \ne 0)\) of the potential
  !> and the distribution function to 0. Noise can be added to these other
  !> components by setting [[init_g_knobs:phiinit]] > 0. The size of the zonal
  !> flows can be adjusted using [[init_g_knobs:zf_init]].
  !>
  !> Author EGH
  subroutine ginit_restart_zonal_only
    use gs2_save, only: gs2_restore
    use kt_grids, only: naky, ntheta0, aky
    use dist_fn_arrays, only: g, gnew
    use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx, is_idx
    use fields_arrays, only: phinew
    use run_parameters, only: has_phi, has_apar, has_bpar
    use array_utils, only: copy, zero_array
    use warning_helpers, only: is_zero
    implicit none
    integer :: iglo, ik, it, il, is

    !  Load phi and g from the restart file
    call gs2_restore (g, scale, has_phi, has_apar, has_bpar)
    
    ! Set all non-zonal components of phi to 0
    do it = 1, ntheta0
       do ik = 2, naky ! Starting at 2 is the crucial bit!
          phinew(:,it,ik) = cmplx(0.0,0.0)
       end do
    end do
    
    ! Allow adjustment of the size of the zonal flows via the input parameter zf_init
    if (naky > 1 .and. is_zero(aky(1))) phinew(:,:,1) = phinew(:,:,1) * zf_init
    
    do iglo = g_lo%llim_proc, g_lo%ulim_proc
       ik = ik_idx(g_lo,iglo)
       it = it_idx(g_lo,iglo)
       il = il_idx(g_lo,iglo)
       is = is_idx(g_lo,iglo)
       if (ik > 1) g(:, :, iglo) = 0.0
    end do

    ! If phiinit > 0, add some noise
    if (phiinit >  epsilon(0.0)) then 
      call ginit_noise
      g = g + gnew
    end if
    call copy(g, gnew)
  end subroutine ginit_restart_zonal_only

  !> Restart but remove the zonal flow (ky = 0) component upon restarting.
  !>
  !> Author: FvW (copy of EGH code)
  subroutine ginit_restart_no_zonal
    use gs2_save, only: gs2_restore
    use dist_fn_arrays, only: g, gnew
    use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx, is_idx
    use fields_arrays, only: phinew
    use run_parameters, only: has_phi, has_apar, has_bpar
    use array_utils, only: copy
    implicit none
    integer :: iglo, ik, it, il, is

    !Load phi and g from the restart file
    call gs2_restore (g, scale, has_phi, has_apar, has_bpar)
    
    !Set zonal component of phi to 0
    phinew(:,:,1) = cmplx(0.0,0.0)
    
    !Set zonal components of g to zero using 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)
       if (ik == 1) g(:, :, iglo) = 0.0
    end do
    call copy(g, gnew)
  end subroutine ginit_restart_no_zonal

  !> The stationary mode test case is a method of testing the geometry
  !! implementation by initializing the distribution function in a particular way
  !! For collisionless runs with ky=0, we expect the solution to be
  !! time-independent (see section 3.3.1 of Ball MIT Masters thesis)
  !!
  !! JRB
  subroutine ginit_stationary_mode
    use species, only: spec, nspec
    use theta_grid, only: ntgrid, theta, bmag, rhoc
    use kt_grids, only: akx
    use dist_fn_arrays, only: g, gnew, vpa, aj0
    use gs2_layouts, only: g_lo, it_idx, is_idx
    use geometry, only : surf, geom !< Should be Rplot from theta_grid instead of geom%Rpos?
    use array_utils, only: copy
    use constants, only: zi
    implicit none

    complex, dimension (-ntgrid:ntgrid) :: phi
    integer :: iglo
    integer :: ig, ik, it, is

    real, dimension (-ntgrid:ntgrid) :: numeratorSum
    real :: denominatorSum, qinp, R_geo

    numeratorSum(:) = 0 ; denominatorSum = 0
    qinp = surf%q ; R_geo = surf%r_geo
    it = 1
    ik = 1
    do is = 1, nspec
      do ig = -ntgrid, ntgrid
        numeratorSum(ig) = numeratorSum(ig) + ((spec(is)%z)**2)*spec(is)%dens*exp(-((spec(is)%mass*spec(is)%temp)/(4*(spec(is)%z)**2))*((qinp/rhoc)*akx(it)*geom%Rpos(rhoc,theta(ig)))**2)
      end do
      denominatorSum = denominatorSum + (spec(is)%z)**2*spec(is)%dens/spec(is)%temp
    end do

    phi(:) = numeratorSum(:)/denominatorSum

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

       g(:,1,iglo) = spec(is)%z*exp(-zi*(qinp/rhoc)*akx(it)*vpa(:,1,iglo)*R_geo*(spec(is)%mass/(spec(is)%z*bmag(:)))) &
                   - (spec(is)%z/spec(is)%temp)*aj0(:,iglo)*phi(:)
       g(:,2,iglo) = spec(is)%z*exp(-zi*(qinp/rhoc)*akx(it)*vpa(:,2,iglo)*R_geo*(spec(is)%mass/(spec(is)%z*bmag(:)))) &
                   - (spec(is)%z/spec(is)%temp)*aj0(:,iglo)*phi(:)
    end do
    call copy(g, gnew)
  end subroutine ginit_stationary_mode

  !> Sets g to a constant value everywhere except zonal component which is 0
  !> Default initial conditions for CGYRO
  !>
  !> if phiinit>0:
  !>     g(ky!=0) = constant
  !> else:
  !>     g(ky!=0) = constant/n
  !>
  !> if zfinit!=0:
  !>     g(ky==0) = constant * exp(-|kx|/max(kx))
  !> else:
  !>     g(ky==0) = 0
  !>
  !> Set by having 'ginit_option' = cgyro
  !> Equivalent CGYRO parameters achieved by setting
  !>    GS2   =    CGYRO
  !> phiinit = AMP / sqrt(2)
  !> zfinit   = AMP0 / AMP
  !>
  !> By default AMP0=0 and AMP=0.1 in CGYRO
  !>
  !> Potential issue in AMP=0 in CGYRO due to zfinit defn
  !>
  !> Author: BSP (copy of EGH code)
  subroutine ginit_cgyro
    use mp, only: proc0
    use le_grids, only: forbid
    use kt_grids, only: akx, aky
    use dist_fn_arrays, only: g, gnew
    use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx
    use array_utils, only: copy
    use warning_helpers, only: is_zero
    implicit none
    integer :: iglo
    integer :: ik, it, il, in
    real :: arg, kxrange

    kxrange = maxval(akx) - minval(akx)

    if(debug.and.proc0) write(*,*) 'Entering CGYRO def'
    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)
       if (is_zero(aky(ik))) then
          arg = abs(akx(it))/ kxrange
          g(:,1,iglo) = zf_init * abs(phiinit) * exp(-arg)
       else
          ! Ensure we don't end up with in = 0. This shouldn't
          ! be needed for flux tube (box) simulations as ik == 1
          ! would correspond to the zonal mode, handled in the
          ! earlier branch.
          if (phiinit < 0.0 .and. ik > 1) then
             in = ik - 1
          else
             in = 1
          end if
          g(:,1,iglo) = abs(phiinit) / in**2
       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_cgyro

  !> Construct the initial condition from a sum of sine waves scaled
  !> with random amplitude up to phiinit (ky > 0) or zf_init (ky ==
  !> 0). Use mode numbers up to minimum of max_mode or last resolved
  !> wave (i.e. ntheta/8).
  !>
  !> A sine wave is selected in order to ensure that the initial
  !> distribution satisfies the parallel boundary conditions
  !> (i.e. that either gnew or h = 0).
  subroutine ginit_random_sine
    use theta_grid, only: ntgrid, ntheta, theta
    use kt_grids, only: naky, ntheta0, aky, enforce_reality, akx
    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 constant_random, only: init_constant_random, finish_constant_random
    use mp, only: proc0, broadcast
    use array_utils, only: copy, zero_array
    use warning_helpers, only: is_zero
    use species, only: spec
    implicit none
    complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi
    real, dimension(-ntgrid:ntgrid) :: sin_theta_mode
    real :: a, b
    real, parameter :: sqrt_two = sqrt(2.0)
    integer :: iglo, ik, it, il, imode, the_max_mode

    call zero_array(phi)

    ! Make sure the max mode is at least 1, but no bigger than the smaller
    ! of the input max_mode and ntheta/8
    the_max_mode = max(1, min(max_mode, ntheta/8))

    ! We choose to only initialise on proc0 and then broadcast to make
    ! sure that all procs have the same values. If done on a processor
    ! by processor basis, then it's possible that different processors
    ! would produce different values. This could lead to a non-Maxwellian
    ! perturbation (i.e. something not constant in our velocity space).
    if (proc0) then
       if (constant_random_flag) call init_constant_random
       ! Here we go from 2*min_mode (min_mode == 1) to 2*the_max_mode
       ! but then divide imode by two in order to allow us to capture
       ! n + 1/2 mode numbers and hence provide an even parity
       ! contribution to the initial condition.
       do imode = 2, the_max_mode*2
          sin_theta_mode = sin(imode*(theta-theta(-ntgrid))/2.0)
          do ik = 1, naky
             do it = 1, ntheta0
                call set_a_and_b_uniform_random_weighting(a, b)
                phi(:, it, ik) = phi(:, it, ik) + cmplx(a*sin_theta_mode, b*sin_theta_mode)
             end do
          end do
       end do
       if (constant_random_flag) call finish_constant_random
       if (enforce_reality) call apply_reality(phi)
    end if

    call broadcast(phi)

    ! Now set the distribution function from "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)
       if (is_zero(aky(ik))) then
          g(:, 1, iglo) = zf_init * phiinit * phi(:, it, ik)

          ! Zero the ky=kx=0 mode
          if (is_zero(akx(it))) g(:, 1, iglo) = 0.0
       else
          g(:, 1, iglo) = phiinit * phi(:, it, ik)
       end if
       g(:, 1, iglo) = g(:, 1, iglo) * spec(is_idx(g_lo, iglo))%z
       where (forbid(:, il)) g(:, 1, iglo) = 0.0
       g(:, 2, iglo) = g(:, 1, iglo)
    end do
    call copy(g, gnew)
  contains
    !> Real and imaginary weights are uniformly distributed between
    !> -1/sqrt(2) and 1/sqrt(2) using either the main random number
    !> generator or the constant random one.
    subroutine set_a_and_b_uniform_random_weighting(a, b)
      use ran, only: ranf
      use constant_random, only: constant_ranf=>ranf
      implicit none
      real, intent(out) :: a, b
      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
      a = a * sqrt_two
      b = b * sqrt_two
    end subroutine set_a_and_b_uniform_random_weighting

  end subroutine ginit_random_sine

  !> Initialise the distribution function to a constant but satisfying
  !> forbidden regions. Can scale the zonal component separately as
  !> usual.
  subroutine ginit_constant
    use kt_grids, only: aky, akx
    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
    use warning_helpers, only: is_zero
    use species, only: spec
    implicit none
    integer :: iglo, ik, it, il

    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)
       if (is_zero(aky(ik))) then
          g(:, 1, iglo) = zf_init * phiinit

          ! Zero the ky=kx=0 mode
          if (is_zero(akx(it))) g(:, 1, iglo) = 0.0
       else
          g(:, 1, iglo) = phiinit
       end if
       g(:, 1, iglo) = g(:, 1, iglo) * spec(is_idx(g_lo, iglo))%z
       where (forbid(:, il)) g(:, 1, iglo) = 0.0
       g(:, 2, iglo) = g(:, 1, iglo)
    end do
    call copy(g, gnew)
  end subroutine ginit_constant

  !> Helper function for common approach to setting g from a passed phi array
  subroutine set_g_from_moment_facs(g, phiinit, phi, dfac_species_weight, &
       ufac_species_weight, use_even, do_reality)
    use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx, is_idx
    use theta_grid, only: ntgrid, theta
    use kt_grids, only: naky, ntheta0
    use dist_fn_arrays, only: vpa, vperp2
    use le_grids, only: forbid
    use species, only: spec, nspec
    use constants, only: zi
    complex, dimension(-ntgrid:, :, g_lo%llim_proc:), intent(in out) :: g
    real, intent(in) :: phiinit ! Shadows module level
    complex, dimension (-ntgrid:ntgrid,ntheta0,naky), intent(in out) :: phi
    logical, intent(in) :: dfac_species_weight, ufac_species_weight, use_even, do_reality
    complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: odd
    real, dimension(-ntgrid:ntgrid) :: dfac, ufac, tparfac, tperpfac, ct, st, c2t, s2t
    integer :: iglo, ik, it, il, is
    real, dimension(nspec) :: dfac_weight, ufac_weight
    if (use_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

    dfac_weight = 1.0
    if (dfac_species_weight) dfac_weight = spec%dens0
    ufac_weight = 1.0
    if (ufac_species_weight) ufac_weight = spec%u0

    odd = zi * phi

    ! reality condition for k_theta = 0 component:
    if (do_reality) then
       call apply_reality(phi)
       call apply_reality(odd)
    end if

    !$OMP PARALLEL DO DEFAULT(none) &
    !$OMP PRIVATE(iglo, ik, it, il, is) &
    !$OMP SHARED(g_lo, g, phiinit, dfac, dfac_weight, phi, &
    !$OMP ufac, ufac_weight, vpa, odd, tparfac, tperpfac, vperp2, forbid) &
    !$OMP SCHEDULE(static)
    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*dfac_weight(is)                 * phi(:,it,ik) &
            + 2*ufac*ufac_weight(is)*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*dfac_weight(is)                 * phi(:,it,ik) &
            + 2*ufac*ufac_weight(is)*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
    !$OMP END PARALLEL DO
  end subroutine set_g_from_moment_facs

  !> Apply the reality condition to the zonal mode (kx = conj(-kx))
  subroutine apply_reality(field)
    use kt_grids, only: ntheta0
    complex, dimension(:, :, :), intent(in out) :: field
    integer :: it
    !$OMP PARALLEL DO DEFAULT(none) &
    !$OMP PRIVATE(it) &
    !$OMP SHARED(ntheta0, field) &
    !$OMP SCHEDULE(static)
    do it = 1, ntheta0 / 2
       field(:, it+(ntheta0+1)/2, 1) = conjg(field(:, (ntheta0+1)/2+1-it, 1))
    end do
    !$OMP END PARALLEL DO
  end subroutine apply_reality

  subroutine do_chop_side_slice(field_slice)
    use theta_grid, only: ntgrid
    complex, dimension(-ntgrid:), intent(in out) :: field_slice
    if (chop_side .and. left) field_slice(:-1) = 0.0
    if (chop_side .and. .not. left) field_slice(1:) = 0.0
  end subroutine do_chop_side_slice

  subroutine do_chop_side_full(field)
    use theta_grid, only: ntgrid
    complex, dimension(-ntgrid:,:,:), intent(in out) :: field
    if (chop_side .and. left) field(:-1, :, :) = 0.0
    if (chop_side .and. .not. left) field(1:, :, :) = 0.0
  end subroutine do_chop_side_full

  !> FIXME : Add documentation
  subroutine single_initial_kx(phi)
    use theta_grid, only: ntgrid
    use kt_grids, only: naky, ntheta0
    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 /= ikx_init) then
         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

  !> FIXME : Add documentation
  subroutine flae (g, gavg)
    use species, only: spec, is_electron_species
    use theta_grid, only: ntgrid, field_line_average
    use kt_grids, only: aky
    use gs2_layouts, only: g_lo, is_idx, ik_idx
    use array_utils, only: zero_array
    use warning_helpers, only: is_not_zero
    complex, dimension (-ntgrid:,:,g_lo%llim_proc:), intent (in) :: g
    complex, dimension (-ntgrid:,:,g_lo%llim_proc:), intent (out) :: gavg
    integer :: iglo
    call zero_array(gavg)

    !$OMP PARALLEL DO DEFAULT(NONE)&
    !$OMP PRIVATE(iglo) &
    !$OMP SHARED(g_lo, g, gavg, spec, aky) &
    !$OMP SCHEDULE(static)
    do iglo = g_lo%llim_proc, g_lo%ulim_proc
       if (.not. is_electron_species(spec(is_idx(g_lo, iglo)))) cycle
       if (is_not_zero(aky(ik_idx(g_lo, iglo)))) cycle
       gavg(:, 1, iglo) = field_line_average(g(:, 1, iglo))
       gavg(:, 2, iglo) = field_line_average(g(:, 2, iglo))
    end do
    !$OMP END PARALLEL DO
  end subroutine flae

  !> FIXME : Add documentation    
  subroutine finish_init_g
    implicit none
    initialized = .false.
  end subroutine finish_init_g

  !> Returns the value of [[init_g_knobs:restart_dir]]
  function get_restart_dir()
    character(len=:), allocatable :: get_restart_dir
    get_restart_dir = restart_dir
  end function get_restart_dir

  !> Sets the value of the module level restart_dir
  subroutine set_restart_dir(restart_dir_in)
    character(len=*), intent(in) :: restart_dir_in
    restart_dir = restart_dir_in
  end subroutine set_restart_dir

#include "init_g_auto_gen.inc"
end module init_g

gyrokinetics / gs2 / 2021218321

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous