ingen_mod.fpp Source File


Contents

Source Code


Source Code

!> Reads in namelists, checks for consistency of inputs, writes a report
!> to runname.report, checks for requested variations,
!> generates new input files, and exits.
!>
!> Consistency checks/reports:
module ingen_mod
  use abstract_config, only: abstract_config_type, CONFIG_MAX_NAME_LEN
  use diagnostics_config, only: diagnostics_type
  implicit none

  private

  public :: sweet_spots, run_ingen, init_ingen, finish_ingen, report
  public :: ingen_config_type, set_ingen_config, get_ingen_config
  
  !> FIXME : Add documentation
  type proc_layout_type
    integer :: nproc
    real :: percent_xxf_2_yxf
    logical :: should_use_unbalanced_xxf
    real :: percentage_xxf_unbalanced_amount
    logical :: should_use_unbalanced_yxf
    real :: percentage_yxf_unbalanced_amount
  end type proc_layout_type

  type(proc_layout_type), dimension(:), allocatable :: sweet_spots
  type(diagnostics_type) :: gnostics
  logical, parameter :: debug=.false.
  integer :: interactive_record, interactive_input, ncut, npmax
  logical :: write_nml, scan, stdin, coll_on = .false., initialized = .false.

  !> Used to represent the input configuration of ingen
  type, extends(abstract_config_type) :: ingen_config_type
     ! namelist : ingen_knobs
     ! indexed : false
     !> This sets the minimum number of local elements for which
     !> processor recommendations will be given. In other words, if
     !> the total number of elements is `nmesh =
     !> (2*ntgrid+1)*2*nlambda*negrid*ntheta0*naky*nspec`, then the
     !> largest processor count that ingen will recommend would be
     !> `nmesh/ncut`.
     integer :: ncut = 100000
     !> Sets the maximum processor count considered when calculating
     !> sweetspots.
     integer :: npmax = 100000     
     !> If `true` then alongside standard [[ingen]] report allows the
     !> creation of various types of scans. This will involve user
     !> interaction if `stdin` is `true`.
     logical :: scan = .false.
     !> If `true` (default) then ask for input from stdin when using
     !> ingen to produce parameters for scans (`scan = .true.`). If
     !> `false` then will look in file named `.<run_name>.pythonin`
     !> instead.
     logical :: stdin = .true.
     !> If `true` then write out a subset of the namelists that have
     !> been read in. Original default was `true` but changed to `false`
     !> as utility of output less clear.
     logical :: write_nml = .false.
   contains
     procedure, public :: read => read_ingen_config
     procedure, public :: write => write_ingen_config
     procedure, public :: reset => reset_ingen_config
     procedure, public :: broadcast => broadcast_ingen_config
     procedure, public, nopass :: get_default_name => get_default_name_ingen_config
     procedure, public, nopass :: get_default_requires_index => get_default_requires_index_ingen_config
  end type ingen_config_type
  
  type(ingen_config_type) :: ingen_config

contains

  !> Entry point: initialise needed modules, run checks, and write namelists
  subroutine run_ingen(ingen_config_in, input_file)
    use file_utils, only: init_file_utils
    use standard_header, only: standard_header_type
    use mp, only: init_mp, finish_mp, proc0
    implicit none
    !> Input options. Defaults to reading the 'ingen_knobs' namelist from the input file
    type(ingen_config_type), intent(in), optional :: ingen_config_in
    !> Filename of input file. Default is to read from the command line
    character(len=*), intent(in), optional :: input_file
    type(standard_header_type) :: local_header
    logical :: list
    call init_mp

    local_header = standard_header_type()

    call init_file_utils (list, name="template", input_file=input_file)
    call init_ingen(ingen_config_in)
 
    if (debug) write(6,*) 'ingen: call get_namelists'
    call get_namelists

    if (debug) write(6,*) 'ingen: call report'
    if (proc0) call report(local_header)

    if (debug) write(6,*) 'ingen: call write_namelists=',write_nml
    if (write_nml) call write_namelists(header=local_header)
    
    if (debug) write(6,*) 'ingen: call interactive, scan=',scan
    if (scan) call interactive
    
    call finish_ingen
    call finish_mp
  end subroutine run_ingen

  !> FIXME : Add documentation  
  subroutine init_ingen(ingen_config_in)
    implicit none
    type(ingen_config_type), intent(in), optional :: ingen_config_in
    
    if(initialized) return
    initialized = .true.

    if (present(ingen_config_in)) ingen_config = ingen_config_in

    call ingen_config%init(name = 'ingen_knobs', requires_index = .false.)

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

    allocate(sweet_spots(0))
  end subroutine init_ingen

  !> FIXME : Add documentation  
  subroutine finish_ingen
    if(allocated(sweet_spots)) deallocate(sweet_spots)
    call ingen_config%reset()
    initialized = .false.
  end subroutine finish_ingen

  !> FIXME : Add documentation  
  subroutine interactive(header)
    use species, only: spec, nspec, has_electron_species
    use geometry, only: beta_prime_input, bishop
    use run_parameters, only: beta, fapar, fbpar
    use standard_header, only: standard_header_type
    use file_utils, only: run_name
    !> Header for files with build and run information
    type(standard_header_type), intent(in), optional :: header
    ! Actual value for optional `header` input
    type(standard_header_type) :: local_header
    character(100) :: pythonin
    integer :: sel, nbeta, j, ise, bishop_save, is
    real :: beta_low, beta_high, dbeta, beta_save, alt, aln, fac, beta_prime_save
    real :: fapar_save, fbpar_save, pri, pe, alpi, tpe_save, ptot, alp, dbdr
    real, dimension (:), allocatable :: tp_save, fp_save
    character (500) :: tag1, tag2
    logical, save :: first = .true.
    pythonin = "."//trim(run_name)//".pythonin"
    if (first) then
       if (present(header)) then
         local_header = header
       else
         local_header = standard_header_type()
       end if

       open (newunit=interactive_record, file='.'//trim(run_name)//".record")
       first = .false.

       if (.not. stdin) then
          open (newunit=interactive_input, file=trim(pythonin))
       else
          interactive_input = 5
       end if
    end if

    call tell ('Interactive specification of a parameter scan')

100 continue
    
    call text
    call text ('Choose a parameter that you would like to vary (1-6):')
    call text ('(1) beta            (4) temperature gradient')
    call text ('(2) beta_prime      (5) density gradient')
    call text ('(3) collisionality  (6) Z_effective')
    call text
    call get_choice (6, sel)
    
    select case (sel)
       
    case default
       call tell ('Try again.  Choose an integer between 1 and 6, inclusively.')
       goto 100

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  1.0  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    case (1) 
       call tell ('You have chosen to vary beta.')

101    continue

       call text
       call text ('Choose from the following:')
       call text ('(1) Vary beta self-consistently')
       call text ('(2) Vary beta with all other parameters held fixed (not self-consistently).')
       call text
       call get_choice (2, sel)

       select case (sel)

       case default
          call tell ('Try again.  Choose an integer between 1 and 2, inclusively.')
          goto 101

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  1.1  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

       case (1) 
          call tell ('You have chosen to vary beta self-consistently.')

102       continue
          call text
          call text ('Choose from the following:')
          call text ('(1) Hold beta_prime fixed, vary electron temperature gradient scale length')
          call text ('(2) Hold beta_prime fixed, vary all temperature gradient scale lengths by same factor')
          call text ('(3) Hold beta_prime fixed, vary all density gradient scale lengths by same factor')
          call text

          call get_choice (3, sel)

          select case (sel)
             
          case default
             call tell ('Try again.  Choose an integer between 1 and 2, inclusively.')
             goto 102
             
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  1.1.1  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

          case (1)  
             call tell ('You have chosen to vary beta and electron temperature gradient at fixed beta_prime.')

             call beta_range_low (beta_low, 'le', 0.)
             call beta_range_high (beta_high, 'le', beta_low)
             call num_runs (nbeta)

             call tell ('Preparing a self-consistent beta scan at fixed beta_prime.', &
                  'The electron temperature gradient scale length will be varied', &
                  'to maintain consistency.')

             call run_number (sel, nbeta)

             write (tag1, fmt='(i3," runs prepared with beta_min = ",e16.10,&
                  &" and beta_max = ",e16.10)') nbeta, beta_low, beta_high
             write (tag2, fmt='("Files are numbered from ",i3," to ",i3)') sel, sel+nbeta-1

             call tell (tag1, tag2)

             ptot = 0.
             alp = 0.
             pri = 0.
             pe = 0.
             alpi = 0.
             ise = 0
             do is=1,nspec
                if (spec(is)%type == 2) then
                   pe = spec(is)%dens * spec(is)%temp
                   ise = is
                else
                   pri = pri + spec(is)%dens * spec(is)%temp
                   alpi = alpi + spec(is)%dens * spec(is)%temp *(spec(is)%fprim + spec(is)%tprim)
                endif
                ptot = ptot + spec(is)%dens * spec(is)%temp
                alp = alp + spec(is)%dens * spec(is)%temp *(spec(is)%fprim + spec(is)%tprim)
             end do
             
             if (.not. has_electron_species(spec)) call tell ('You really should use electrons for electromagnetic runs.')

             alp = alp/ptot
             dbdr = - beta*ptot*alp

             dbeta = (beta_high-beta_low)/(nbeta-1)
             do j = sel, sel+nbeta-1
                
                beta_save = beta
                beta = beta_low + (j - sel)*dbeta
                
                tpe_save = spec(ise)%tprim
                spec(ise)%tprim = - (spec(ise)%fprim + alpi/pe + dbdr/beta/pe)

                write (tag1, fmt='("Varying beta and L_Te self-consistently with& 
                     & beta_prime fixed")') 

                write (tag2, fmt='("beta = ",e16.10," and electron tprim = ",e16.10)') beta, spec(ise)%tprim 

                call write_namelists (j, tag1, tag2, local_header)
                spec(ise)%tprim = tpe_save
             end do
             beta = beta_save

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  1.1.2  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

          case (2)

             call tell ('You have chosen to vary beta and all temperature &
                  &gradient scale lengths together at fixed beta_prime.')

             call beta_range_low (beta_low, 'le', 0.)
             call beta_range_high (beta_high, 'le', beta_low)
             call num_runs (nbeta)

             call tell ('Preparing a self-consistent beta scan at fixed beta_prime.', &
                  'All temperature gradient scale lengths will be varied', &
                  'by the same factor to maintain consistency.')

             call run_number (sel, nbeta)

             write (tag1, fmt='(i3," runs prepared with beta_min = ",e16.10,&
                  &" and beta_max = ",e16.10)') nbeta, beta_low, beta_high
             write (tag2, fmt='("Files are numbered from ",i3," to ",i3)') sel, sel+nbeta-1

             call tell (tag1, tag2)

             allocate (tp_save (nspec))

             ptot = 0.
             alt = 0.
             aln = 0.
             do is=1,nspec
                ptot = ptot + spec(is)%dens * spec(is)%temp
                alt = alt + spec(is)%dens * spec(is)%temp *(spec(is)%tprim)
                aln = aln + spec(is)%dens * spec(is)%temp *(spec(is)%fprim)
             end do
             
             if (.not. has_electron_species(spec)) call tell ('You really should use electrons for electromagnetic runs.')

             alp = (alt+aln)/ptot
             dbdr = - beta*ptot*alp

             dbeta = (beta_high-beta_low)/(nbeta-1)
             do j = sel, sel+nbeta-1
                
                beta_save = beta
                beta = beta_low + (j - sel)*dbeta
                
                fac = -(dbdr/beta+aln)/alt
                tp_save = spec%tprim
                spec%tprim = fac*spec%tprim

                write (tag1, fmt='("Varying beta and all L_T values self-consistently")')
                write (tag2, fmt='("beta = ",e16.10," and tprim values scaled by ",e16.10)') beta, fac

                call write_namelists (j, tag1, tag2, local_header)
                spec%tprim = tp_save
             end do
             beta = beta_save

             deallocate (tp_save)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  1.1.3  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

          case (3)

             call tell ('You have chosen to vary beta and all density &
                  &gradient scale lengths together at fixed beta_prime.')

             call beta_range_low (beta_low, 'le', 0.)
             call beta_range_high (beta_high, 'le', beta_low)
             call num_runs (nbeta)

             call tell ('Preparing a self-consistent beta scan at fixed beta_prime.', &
                  'All density gradient scale lengths will be varied', &
                  'by the same factor to maintain consistency.')

             call run_number (sel, nbeta)

             write (tag1, fmt='(i3," runs prepared with beta_min = ",e16.10,&
                  &" and beta_max = ",e16.10)') nbeta, beta_low, beta_high
             write (tag2, fmt='("Files are numbered from ",i3," to ",i3)') sel, sel+nbeta-1

             call tell (tag1, tag2)

             allocate (fp_save (nspec))

             ptot = 0.
             alt = 0.
             aln = 0.
             do is=1,nspec
                ptot = ptot + spec(is)%dens * spec(is)%temp
                alt = alt + spec(is)%dens * spec(is)%temp *(spec(is)%tprim)
                aln = aln + spec(is)%dens * spec(is)%temp *(spec(is)%fprim)
             end do
             
             if (.not. has_electron_species(spec)) call tell ('You really should use electrons for electromagnetic runs.')

             alp = (alt+aln)/ptot
             dbdr = - beta*ptot*alp

             dbeta = (beta_high-beta_low)/(nbeta-1)
             do j = sel, sel+nbeta-1
                
                beta_save = beta
                beta = beta_low + (j - sel)*dbeta
                 
                fac = -(dbdr/beta+alt)/aln
                fp_save = spec%fprim
                spec%fprim = fac*spec%fprim

                write (tag1, fmt='("Varying beta and all L_n values self-consistently")')
                write (tag2, fmt='("beta = ",e16.10," and tprim values scaled by ",e16.10)') beta, fac

                call write_namelists (j, tag1, tag2, local_header)
                spec%fprim = fp_save
             end do
             beta = beta_save

             deallocate (fp_save)
          end select

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  1.2  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

       case (2)
          call tell ('You have selected to vary beta non-self-consistently.')
          
          call beta_range_low (beta_low, 'lt', 0.)
          call beta_range_high (beta_high, 'le', beta_low)
          call num_runs (nbeta)

          call tell ('Preparing a non-self-consistent beta scan.')

          call run_number (sel, nbeta)

          write (tag1, fmt='(i3," runs prepared with beta_min = ",e16.10,&
               &" and beta_max = ",e16.10)') nbeta, beta_low, beta_high
          write (tag2, fmt='("Files are numbered from ",i3," to ",i3)') sel, sel+nbeta-1
          
          call tell (tag1, tag2)

          dbeta = (beta_high-beta_low)/(nbeta-1)
          do j = sel, sel+nbeta-1
             
             beta_save = beta
             fapar_save = fapar 
             fbpar_save = fbpar

             beta = beta_low + (j - sel)*dbeta
             if (beta == 0.) then 
                fapar = 0.
                fbpar = 0.
             else
                if (fapar == 0. .and. fbpar == 0.) then
                   fapar = 1.0 ;  fbpar = 1.0
                end if
             end if
             
             write (tag1, fmt='("Varying beta, all else fixed")')
             write (tag2, fmt='("beta = ",e16.10)') beta

             call write_namelists (j, tag1, tag2, local_header)

             fapar = fapar_save 
             fbpar = fbpar_save
             beta = beta_save
          end do
          
       end select

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  2.0  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    case (2) ! beta_prime
       
       call tell ('You have chosen to vary beta_prime.')

115    continue
       call text
       call text ('Choose from the following:')
       call text ('(1) Vary beta_prime self-consistently')
       call text ('(2) Vary beta_prime with ALL other parameters held fixed (non-self-consistently).')
       call text
       call get_choice (2, sel)

       select case (sel)

       case default
          call tell ('Try again.  Choose an integer between 1 and 2, inclusively.')
          goto 115

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  2.1  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

       case (1)
          call tell ('You have chosen to vary beta_prime self-consistently.')

116       continue
          call text
          call text ('Choose from the following:')
          call text ('(1) Hold gradient scale lengths fixed, vary beta')
          call text

          call get_choice (1, sel)

          select case (sel)

          case default
             call tell ('Try again.  Choose an integer between 1 and 1, inclusively.')
             goto 116

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  2.1.1  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

          case (1)  
             call tell ('You have chosen to vary beta_prime while holding gradient scale lengths fixed.')

             call beta_prime_range_low (beta_low)
             call beta_prime_range_high (beta_high, beta_low)
             call num_runs (nbeta)

             call tell ('Preparing a self-consistent beta_prime scan.', &
                  'Beta will be varied to maintain consistency.')

             call run_number (sel, nbeta)

             write (tag1, fmt='(i3," runs prepared with beta_prime_min = ",e16.10,&
                  &" and beta_prime_max = ",e16.10)') nbeta, beta_low, beta_high
             write (tag2, fmt='("Files are numbered from ",i3," to ",i3)') sel, sel+nbeta-1
             
             call tell (tag1, tag2)

             ptot = 0.
             alp = 0.
             do is=1,nspec
                ptot = ptot + spec(is)%dens * spec(is)%temp
                alp = alp + spec(is)%dens * spec(is)%temp *(spec(is)%fprim + spec(is)%tprim)
             end do
            
             alp = alp/ptot
             dbdr = - beta*ptot*alp

             if (alp == 0.) then
                call tell ('Cannot proceed, because Lp = infinity', &
                     'No input files for scan written')
                return
             end if

             beta_save = beta
             beta_prime_save = dbdr

             fac = -1./(ptot*alp)

             dbeta = (beta_high-beta_low)/(nbeta-1)   ! actually, this is dbeta_prime
             do j = sel, sel+nbeta-1
                
                beta_prime_save = beta_prime_input
                beta_prime_input = beta_low + (j - sel)*dbeta
                
                beta_save = beta
                beta = beta_prime_input*fac

                fapar_save = fapar ; fbpar_save = fbpar
                if (beta == 0.) then
                   fapar = 0.      ; fbpar = 0.
                else
                   if (fapar == 0. .and. fbpar == 0.) then
                      fapar = 1.0 ;  fbpar = 1.0
                   end if
                end if

                select case (bishop)
                case default
                   bishop_save = bishop
                   bishop = 6
                case (4)
                   ! nothing, continue to use bishop = 4
                end select

             
                write (tag1, fmt='("Varying beta_prime and beta self-consistently")')
                write (tag2, fmt='("beta_prime = ",e16.10," and beta = ",e16.10)') beta_prime_input, beta

                call write_namelists (j, tag1, tag2, local_header)

                fapar = fapar_save 
                fbpar = fbpar_save
                beta = beta_save
                beta_prime_input = beta_prime_save

                select case (bishop)
                case default
                   bishop = bishop_save
                case (4)
                   ! nothing, continue to use bishop = 4
                end select

             end do
          end select

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  2.2   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

       case (2)

          call tell ('You have selected to vary beta_prime non-self-consistently.')          
          
          call beta_prime_range_low (beta_low)
          call beta_prime_range_high (beta_high, beta_low)
          call num_runs (nbeta)

          call tell ('Preparing a non-self-consistent beta_prime scan.')

          call run_number (sel, nbeta)
          
          write (tag1, fmt='(i3," runs prepared with beta_prime_min = ",e16.10,& 
               &" and beta_prime_max = ",e16.10)') nbeta, beta_low, beta_high 
          write (tag2, fmt='("Files are numbered from ",i3," to ",i3)') sel, sel+nbeta-1 
          
          call tell (tag1, tag2)

          dbeta = (beta_high-beta_low)/(nbeta-1) 
          do j = sel, sel+nbeta-1
             
             beta_prime_save = beta_prime_input
             beta_prime_input = beta_low + (j - sel)*dbeta
             
             select case (bishop)
             case default
                bishop_save = bishop
                bishop = 6
             case (4)
                ! nothing, continue to use bishop = 4
             end select
             
             write (tag1, fmt='("Varying beta_prime only (non-self-consistently)")')
             write (tag2, fmt='("beta_prime = ",e16.10)') beta_prime_input

             call write_namelists (j, tag1, tag2, local_header)
             
             beta_prime_input = beta_prime_save
             
             select case (bishop)
             case default
                bishop = bishop_save
             case (4)
                ! nothing, continue to use bishop = 4
             end select
             
          end do
          
       end select

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   1.3  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    case (3) ! collisionality
       call tell ('You have chosen to vary collisionality.')
       call text ('Not yet implemented (sorry).')
       call text

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   1.4  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    case (4) ! temperature gradient
       call tell ('You have chosen to vary temperature gradient.')
       call text ('Not yet implemented (sorry).')
       call text

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   1.5  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    case (5) ! density gradient
       call tell ('You have chosen to vary density gradient.')
       call text ('Not yet implemented (sorry).')
       call text

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   1.6  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    case (6) ! Z_effective
       call tell ('You have chosen to vary Z_effective.')
       call text ('Not yet implemented (sorry).')
       call text
       
    end select

  end subroutine interactive

  !> FIXME : Add documentation
  subroutine get_namelists
!CMR, 17/11/2009: use gs2_diagnostics module to pick up public variables
!                 and routines to reduces maintenance of ingen.
!                 Should replicate this for other modules in future.
!CMR, 2/2/2011:   Have extended this to include use of theta_grid module:
!                 Strategy is simply to add two types of routines to modules:
!                      wnml_xxxxx   to write the modules namelists
!                      check_xxxxx  to perform the ingen checking inside the module
!                 More object oriented strategy, easier maintenance of ingen.f90 
!                 which is gradually shrinking.!                  
!
    use antenna, only: init_antenna
    use collisions, only: coll_read_parameters=>read_parameters
    use collisions, only: collision_model_switch, collision_model_none
    use fields, only : fields_read_parameters=>read_parameters
    use hyper, only : hyper_read_parameters=>read_parameters
    use species,only : spec, init_species
    use gs2_layouts, only: init_gs2_layouts
    use gs2_reinit, only: init_reinit
    use dist_fn, only: dist_fn_read_parameters=>read_parameters
    use nonlinear_terms, only: nonlinear_terms_read_parameters=>read_parameters
    use run_parameters, only: init_run_parameters
    use init_g, only: init_init_g
    use run_parameters, only: use_old_diagnostics
    use gs2_diagnostics,only: gs2diag_read_parameters=>read_parameters
    use diagnostics_config, only: init_diagnostics_config
    use theta_grid, only: init_theta_grid
    use kt_grids, only: init_kt_grids
    use le_grids, only: init_le_grids
    use theta_grid, only: init_theta_grid, ntgrid
    use gs2_time, only: init_gs2_time
    implicit none
    logical, parameter :: debug=.false.

    if (debug) write(6,*) 'get_namelists: layouts'
    call init_gs2_layouts

    call init_run_parameters
if (debug) write(6,*) 'get_namelists: collisions'
    call coll_read_parameters

if (debug) write(6,*) 'get_namelists: init_g'
    call init_init_g
    !This needs to come after init_init_g in case we need to read the
    !timestep from restart file (init_init_g sets the restart_dir etc.)
    call init_gs2_time

if (debug) write(6,*) 'get_namelists: fields'
    call fields_read_parameters
    
if (debug) write(6,*) 'get_namelists: gs2_reinit'
    call init_reinit

if (debug) write(6,*) 'get_namelists: hyper'
    call hyper_read_parameters

if (debug) write(6,*) 'get_namelists: kt_grids'
    call init_kt_grids

if (debug) write(6,*) 'get_namelists: le_grids'
    call init_le_grids

if (debug) write(6,*) 'get_namelists: nonlinear terms'
    call nonlinear_terms_read_parameters

if (debug) write(6,*) 'get_namelists: init_species'
    call init_species

    coll_on = any(spec%vnewk > epsilon(0.0)) &
         .and. (collision_model_switch /= collision_model_none)

    call init_antenna

!CMR, 2/2/2011:  reduce much duplication by calling init_theta_grid
if (debug) write(6,*) 'get_namelists: call init_theta_grid, ntgrid=',ntgrid
    call init_theta_grid
if (debug) write(6,*) 'get_namelists: done init_theta_grid, ntgrid=',ntgrid

    call dist_fn_read_parameters
    ! gs2_diagnostics

! CMR 18/11/2009:  reduce duplication by calling gs2diag_read_parameters 
!
   if (use_old_diagnostics) then
      call gs2diag_read_parameters(.false.)
   else
      call init_diagnostics_config(gnostics)
   end if

if (debug) write(6,*) 'get_namelists: returning'

  end subroutine get_namelists

  !> FIXME : Add documentation  
  subroutine write_namelists (jr, tag1, tag2, header)
    use antenna, only: wnml_antenna
    use collisions, only: wnml_collisions
    use dist_fn, only: wnml_dist_fn, wnml_dist_fn_species 
    use fields, only: wnml_fields
    use gs2_reinit, only: wnml_gs2_reinit
    use gs2_layouts, only: wnml_gs2_layouts
    use gs2_diagnostics, only: wnml_gs2_diagnostics
    use run_parameters, only: use_old_diagnostics
    use hyper, only: wnml_hyper
    use init_g, only : wnml_init_g
    use kt_grids, only: wnml_kt
    use le_grids, only: wnml_le_grids
    use nonlinear_terms, only: nonlin, wnml_nonlinear_terms
    use run_parameters, only: wnml_run_parameters
    use species, only: wnml_species
    use theta_grid, only: wnml_theta_grid
    use theta_grid_params, only: wnml_theta_grid_params
    use standard_header, only: standard_header_type
    use file_utils, only: run_name
    integer, intent (in), optional :: jr
    character (*), intent (in), optional :: tag1, tag2
    !> Header for files with build and run information
    type(standard_header_type), intent(in), optional :: header
    ! Actual value for optional `header` input
    type(standard_header_type) :: local_header
    logical :: electrons
    integer :: h, t, u, unit
    character (4) :: suffix
    character(len=:), allocatable :: filename
    
    if (present(header)) then
      local_header = header
    else
      local_header = standard_header_type()
    end if

    if (present(jr)) then

       h = jr / 100
       t = (jr - h * 100) / 10
       u = (jr - h * 100 - t * 10)
       suffix = '_'//achar(48+h)//achar(48+t)//achar(48+u)
       filename = trim(run_name) // suffix // ".in"
    else
       filename = trim(run_name) // ".inp"
    endif

    open (newunit=unit, file=filename)
    write (*, *) "Writing processed input file to '" // filename // "'"

    write (unit, *)

    write (unit=unit, fmt='(a)') local_header%to_string(file_description="Input parameters as set by ingen")

    if (present(tag1)) then
       write (unit, *) '*****************************************************'
       write (unit, *) trim(tag1)
       if (present(tag2)) write (unit, *) trim(tag2)
       write (unit, *) '*****************************************************'
    end if

    call wnml_gs2_layouts(unit)
    call wnml_collisions(unit)
    call wnml_init_g(unit)
    call wnml_dist_fn(unit)
    call wnml_fields(unit)
    call wnml_gs2_diagnostics(unit)
    if (nonlin) call wnml_gs2_reinit(unit)
    call wnml_hyper(unit)
    call wnml_kt(unit)
    call wnml_le_grids(unit)
    call wnml_nonlinear_terms(unit)
    call wnml_run_parameters(unit)
    call wnml_species(unit)
    call wnml_dist_fn_species(unit)

    call wnml_theta_grid_params(unit)
    call wnml_theta_grid(unit)

    call wnml_antenna(unit)

    write(unit, fmt=*)
    close (unit)
    
  end subroutine write_namelists

  !> FIXME : Add documentation   
  subroutine pfactors (n, div)
    integer, intent (in) :: n
    integer, dimension (:), intent (out) :: div
    integer, dimension (50), parameter :: primes = (/ &
         2, 3, 5, 7, 11, &
         13, 17, 19, 23, 29, &
         31, 37, 41, 43, 47, &
         53, 59, 61, 67, 71, &
         73, 79, 83, 89, 97, &
         101, 103, 107, 109, 113, &
         127, 131, 137, 139, 149, &
         151, 157, 163, 167, 173, &
         179, 181, 191, 193, 197, &
         199, 211, 223, 227, 229 /)
    
    integer :: i, ntmp
    
    ntmp = n 
    i=1
    do while (ntmp > 1 .and. i < 51)
       do while (mod(ntmp, primes(i)) == 0)
          if (i < 4) div(i) = div(i) + 1
          if (i > 3) div(4) = primes(i)
          ntmp = ntmp / primes(i)
       end do
       i=i+1
    end do
  end subroutine pfactors

  !> FIXME : Add documentation  
  subroutine report(header)
    use antenna, only: check_antenna
    use collisions, only: check_collisions
    use dist_fn, only: check_dist_fn
    use fields, only: check_fields
    use gs2_diagnostics, only: check_gs2_diagnostics
    use gs2_diagnostics, only: nwrite, dump_fields_periodically, save_for_restart
    use gs2_diagnostics, only: nsave, make_movie, nmovie, exit_when_converged, omegatol
    use gs2_reinit, only: delt_adj
    use gs2_time, only: code_dt, user_dt
    use hyper, only: check_hyper
    use init_g, only: check_init_g
    use kt_grids, only: check_kt_grids, gridopt_switch
    use kt_grids, only: gridopt_box, naky, ntheta0, nx, ny
    use le_grids, only: negrid, nlambda
    use nonlinear_terms, only: nonlin, check_nonlinear_terms
    use run_parameters, only: use_old_diagnostics, check_run_parameters
    use run_parameters, only: beta, tite, nstep, wstar_units
    use species, only: check_species, nspec, has_electron_species
    use theta_grid, only: check_theta_grid, ntgrid
    use collisions, only: collision_model_switch
    use collisions, only: collision_model_full, collision_model_ediffuse
    use collisions, only: collision_model_lorentz, collision_model_lorentz_test
    use collisions, only: use_le_layout
    use standard_header, only: standard_header_type
    use file_utils, only: run_name, open_output_file, close_output_file
    implicit none
    !> Header for files with build and run information
    type(standard_header_type), intent(in), optional :: header
    ! Actual value for optional `header` input
    type(standard_header_type) :: local_header
    real :: alne, dbetadrho_spec, local_omegatol
    logical :: local_exit_when_converged, local_dump_fields_periodically
    logical :: local_save_for_restart, local_make_movie
    integer :: nmesh, local_nwrite, local_nmovie, local_nsave, report_unit, i
    integer, dimension(4) :: pfacs

    if (present(header)) then
      local_header = header
    else
      local_header = standard_header_type()
    end if

    call open_output_file (report_unit, ".report")
    write (*, *) "Writing report to '" // trim(run_name) // ".report'"

    write (report_unit, *)

    write (report_unit, fmt='(a)') local_header%to_string(file_description="Input file report from ingen")
    call write_separator(report_unit)

    if (nonlin) then
       if (gridopt_switch /= gridopt_box) then
          write (report_unit, *) 
          write (report_unit, fmt="('################# WARNING #######################')")
          write (report_unit, fmt="('Nonlinear runs must be carried out in a box.')") 
          write (report_unit, fmt="('Set grid_option to box in the kt_grids_knobs namelist')") 
          write (report_unit, fmt="('or set nonlinear_mode to off in the nonlinear_knobs namelist.')") 
          write (report_unit, fmt="('THIS IS AN ERROR.')") 
          write (report_unit, fmt="('################# WARNING #######################')")
          write (report_unit, *) 
       end if

       write (report_unit, fmt="('nmesh=(2*ntgrid+1)*2*nlambda*negrid*nx*ny*nspec')")
       nmesh = (2*ntgrid+1)*2*nlambda*negrid*nx*ny*nspec
       write (report_unit, fmt="('Number of meshpoints:    ',i16)") nmesh

       !
       ! check that nx, ny have no large prime factors
       !
       pfacs = 0
       call pfactors (ny, pfacs)
       if (pfacs(4) /= 0) then
          write (report_unit, *) 
          write (report_unit, fmt="('################# WARNING #######################')")
          write (report_unit, fmt="('ny is a multiple of ',i4)") pfacs(4)
          write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')")
          write (report_unit, fmt="('ny should have only 2, 3, and/or 5 as prime factors.')")
          write (report_unit, fmt="('################# WARNING #######################')")
          write (report_unit, *) 
       end if
       i = 1
       if (pfacs(1) > 0) i=2**pfacs(1)
       if (pfacs(2) > 0) i=3**pfacs(2)*i
       if (pfacs(3) > 0) i=5**pfacs(3)*i
       if (i /= ny) then
          write (report_unit, *) 
          write (report_unit, fmt="('################# WARNING #######################')")
          write (report_unit, fmt="('ny = ',i3)") ny
          write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')")
          write (report_unit, fmt="('ny should have only 2, 3, and/or 5 as prime factors.')")
          write (report_unit, fmt="('################# WARNING #######################')")
          write (report_unit, *) 
       end if

       pfacs = 0
       call pfactors (nx, pfacs)
       if (pfacs(4) /= 0) then
          write (report_unit, *) 
          write (report_unit, fmt="('################# WARNING #######################')")
          write (report_unit, fmt="('nx is a multiple of ',i3)") pfacs(4)
          write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')")
          write (report_unit, fmt="('nx should have only 2, 3, and/or 5 as prime factors.')")
          write (report_unit, fmt="('################# WARNING #######################')")
          write (report_unit, *) 
       end if
       i = 1
       if (pfacs(1) > 0) i=2**pfacs(1)
       if (pfacs(2) > 0) i=3**pfacs(2)*i
       if (pfacs(3) > 0) i=5**pfacs(3)*i
       if (i /= nx) then
          write (report_unit, *) 
          write (report_unit, fmt="('################# WARNING #######################')")
          write (report_unit, fmt="('nx = ',i3)") nx
          write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')")
          write (report_unit, fmt="('nx should have only 2, 3, and/or 5 as prime factors.')")
          write (report_unit, fmt="('################# WARNING #######################')")
          write (report_unit, *) 
       end if
    else
       write (report_unit, fmt="('nmesh=(2*ntgrid+1)*2*nlambda*negrid*ntheta0*naky*nspec')")
       nmesh = (2*ntgrid+1)*2*nlambda*negrid*ntheta0*naky*nspec
       write (report_unit, fmt="('Number of meshpoints:    ',i14)") nmesh
    end if
    
    write (report_unit, fmt="(T12,' ntgrid=',i12)") ntgrid
    write (report_unit, fmt="(T9,'2*ntgrid+1=',i12)") 2*ntgrid+1
    write (report_unit, fmt="(T12,'nlambda=',i12)") nlambda
    write (report_unit, fmt="(T12,' negrid=',i12)") negrid
    write (report_unit, fmt="(T12,'ntheta0=',i12)") ntheta0
    write (report_unit, fmt="(T12,'     nx=',i12)") nx
    write (report_unit, fmt="(T12,'   naky=',i12)") naky
    write (report_unit, fmt="(T12,'     ny=',i12)") ny
    write (report_unit, fmt="(T12,'  nspec=',i12,/)") nspec

    call nprocs (nmesh, report_unit)
    if (nonlin) then
       write (report_unit, fmt="(/'Nonlinear run => consider #proc sweetspots for xxf+yxf objects!')") 
       call nprocs_xxf(report_unit)
       call nprocs_yxf(report_unit)
    endif

    if(use_le_layout) then
       write (report_unit, fmt="(/'Collisions using le_lo')") 
       call nprocs_le(report_unit)
    else
       select case (collision_model_switch)
       case (collision_model_full)
          write (report_unit, fmt="(/'Collisions using lz_lo')") 
          call nprocs_lz(report_unit)
          write (report_unit, fmt="(/'Collisions using e_lo')") 
          call nprocs_e(report_unit)
       case(collision_model_lorentz,collision_model_lorentz_test)
          write (report_unit, fmt="(/'Collisions using lz_lo')") 
          call nprocs_lz(report_unit)
       case (collision_model_ediffuse)
          write (report_unit, fmt="(/'Collisions using e_lo')") 
          call nprocs_e(report_unit)
       end select
    end if

    call write_separator(report_unit)
    call check_species(report_unit,beta,tite,alne,dbetadrho_spec)
    call write_separator(report_unit)
    call check_theta_grid(report_unit,alne,dbetadrho_spec)
    call write_separator(report_unit)

    if (coll_on) then 
       call check_collisions(report_unit) 
    else
       write (report_unit, fmt="('All collisionality parameters (vnewk) are zero.')")
       write (report_unit, fmt="('No collision operator will be used.')")
    end if

    call write_separator(report_unit)
    call check_fields(report_unit) 
    call write_separator(report_unit)
    call check_init_g(report_unit)
    call write_separator(report_unit)
    call check_nonlinear_terms(report_unit,delt_adj)

    if (.not. nonlin) then
       write (report_unit, *)
       write (report_unit, fmt="('This is a linear calculation.')")
       write (report_unit, *)
       write (report_unit, fmt="('The time step (code_dt) = ',e10.4)") code_dt
       write (report_unit, fmt="('The time step (user_dt) = ',e10.4)") user_dt
       write (report_unit, fmt="('The maximum number of time steps is ',i7)") nstep
       if (use_old_diagnostics) then
          local_exit_when_converged = exit_when_converged
          local_omegatol = omegatol
       else
          local_exit_when_converged = gnostics%exit_when_converged
          local_omegatol = gnostics%omegatol
       end if

       if (local_exit_when_converged) then
          write (report_unit, fmt="('When the frequencies for each k have converged, the run will stop.')")
          write (report_unit, fmt="('The convergence has to be better than one part in ',e10.4)") 1./local_omegatol
       end if

       if (wstar_units) then
          write (report_unit, *)
          write (report_unit, fmt="('The timestep for each ky is scaled by a factor of 1/ky.')")
       end if
    end if

    if (use_old_diagnostics) then
       local_nwrite = nwrite
       local_dump_fields_periodically = dump_fields_periodically
       local_make_movie = make_movie
       local_nmovie = nmovie
       local_save_for_restart = save_for_restart
       local_nsave = nsave
    else
       local_nwrite = gnostics%nwrite
       local_dump_fields_periodically = gnostics%dump_fields_periodically
       local_make_movie = gnostics%make_movie
       local_nmovie = gnostics%nmovie
       local_save_for_restart = gnostics%save_for_restart
       local_nsave = gnostics%nsave
    end if

    write (report_unit, *)
    write (report_unit, fmt="('Data will be written to ',a,' every ',i4,' timesteps.')") trim(run_name)//'.out.nc', local_nwrite
    write (report_unit, *)

    if (local_dump_fields_periodically) then
       write (report_unit, *)
       write (report_unit, fmt="('Data will be written to dump.fields.t=* every ',i4,' timesteps.')") 10*local_nmovie
       write (report_unit, *)
    end if

    if (local_make_movie) then
       write (report_unit, *)
       write (report_unit, fmt="('Movie data will be written to runname.movie.nc every ',i4,' timesteps.')") local_nmovie
       write (report_unit, *)
    end if

    if (local_save_for_restart .and. local_nsave > 0) then
       write (report_unit, *)
       write (report_unit, fmt="('Restart data will be written every ',i4,' timesteps.')") local_nsave
       write (report_unit, *)
    end if

    call write_separator(report_unit)
    call check_kt_grids(report_unit)
    call write_separator(report_unit)
    call check_dist_fn(report_unit)
    call write_separator(report_unit)
    call check_antenna(report_unit)
    call check_hyper(report_unit)
    call write_separator(report_unit)
    call check_run_parameters(report_unit)

    ! diagnostic controls:
    if (use_old_diagnostics) then
       call check_gs2_diagnostics(report_unit)
    end if

    call close_output_file (report_unit)

  end subroutine report

  !> Writes a standard section separator
  subroutine write_separator(report_unit)
    implicit none
    integer, intent(in) :: report_unit
    write (report_unit, *)
    write (report_unit, fmt="('------------------------------------------------------------')")
    write (report_unit, *)
  end subroutine write_separator
  
  !> Writes out sweetspot core counts for a layout
  !>  sym: character string label for each dimension
  !>  sdim:  size of each dimension
  !>  nfac:  #factors of each sdim
  !>  facs(i,j):  ith factor of jth dimension
  subroutine wsweetspots(sym,sdim,nfac,facs,npmax,LUN)
    use optionals, only: get_option_with_default
    implicit none
    character(3), dimension(:), intent(in):: sym
    integer, dimension(:), intent(in):: sdim, nfac
    integer, dimension(:,:), intent(in):: facs
    integer, intent(in):: npmax
    integer, intent(in), optional:: LUN
    integer :: lout, npe, i, j, lcores, nfacs
    lout = get_option_with_default(LUN, 6)
    nfacs=size(nfac)
    lcores=1
    write (lout, fmt="('  npe = ',i8,'  (',a,')')") 1,trim(sym(1))
    do i=1,nfacs
       do j=2,nfac(i)
          npe=facs(j,i)*lcores
          if (npe .gt. npmax) exit
          write (lout, fmt="('  npe = ',i8,'  (',a,')')") npe,trim(sym(i))
       enddo
       lcores=lcores*sdim(i)
    enddo
  end subroutine wsweetspots

  !> FIXME : Add documentation  
  subroutine nprocs_xxf(report_unit)
    use nonlinear_terms, only : nonlin
    use species, only : nspec
    use kt_grids, only: naky
    use le_grids, only: negrid, nlambda
    use theta_grid, only: ntgrid
    use gs2_layouts, only: layout, factors
    implicit none
    integer, intent(in) :: report_unit
    integer :: nefacs, nlfacs, nkyfacs, nsgfacs, nspfacs, ntgfacs
    integer, dimension(:,:), allocatable :: facs
    integer :: maxfacs
    integer, allocatable, dimension(:):: spfacs, efacs, lfacs, sgfacs, tgfacs, kyfacs
    integer, dimension(6) :: nfac, sdim
    character(3), dimension(6):: sym

    if (.not.nonlin) return
    write (report_unit, fmt="('xxf sweetspot #proc up to:',i8)") npmax
    maxfacs=max(nspec,negrid,nlambda,2,2*ntgrid+1,naky)/2+1
    allocate (spfacs(maxfacs),efacs(maxfacs),lfacs(maxfacs),sgfacs(maxfacs),tgfacs(maxfacs),kyfacs(maxfacs),facs(maxfacs,6))
    call factors (nspec, nspfacs, spfacs)
    call factors (negrid, nefacs, efacs)
    call factors (nlambda, nlfacs, lfacs)
    call factors (2, nsgfacs, sgfacs)
    call factors (2*ntgrid+1, ntgfacs, tgfacs)
    call factors (naky, nkyfacs, kyfacs)

    select case (layout)
    case ('lexys','lxyes','lyxes','yxles','xyles')
       sym = (/ "s  ","e  ","l  ","sgn","tg ","y  " /)
       sdim = (/ nspec, negrid, nlambda, 2, 2*ntgrid+1, naky /)
       nfac= (/ nspfacs, nefacs, nlfacs, nsgfacs, ntgfacs, nkyfacs /)
       facs(:,1)=spfacs; facs(:,2)=efacs; facs(:,3)=lfacs
       facs(:,4)=sgfacs; facs(:,5)=tgfacs; facs(:,6)=kyfacs
       call wsweetspots(sym,sdim,nfac,facs,npmax,LUN=report_unit)
    case ('yxels')
       sym = (/ "s  ","l  ","e  ","sgn","tg ","y  " /)
       sdim = (/ nspec, nlambda, negrid, 2, 2*ntgrid+1, naky /)
       nfac= (/ nspfacs, nlfacs, nefacs, nsgfacs, ntgfacs, nkyfacs /)
       facs(:,1)=spfacs; facs(:,2)=lfacs; facs(:,3)=efacs
       facs(:,4)=sgfacs; facs(:,5)=tgfacs; facs(:,6)=kyfacs
       call wsweetspots(sym,sdim,nfac,facs,npmax,LUN=report_unit)
    end select
    deallocate (facs,spfacs,efacs,lfacs,sgfacs,tgfacs,kyfacs)
  end subroutine nprocs_xxf

  !> FIXME : Add documentation    
  subroutine nprocs_yxf(report_unit)
    use nonlinear_terms, only : nonlin
    use species, only : nspec
    use kt_grids, only: naky, nx
    use le_grids, only: negrid, nlambda
    use theta_grid, only: ntgrid
    use gs2_layouts, only: layout, factors
    implicit none
    integer, intent(in) :: report_unit
    integer :: nefacs, nlfacs, nkxfacs, nsgfacs, nspfacs, ntgfacs
    integer, dimension(:,:), allocatable :: facs
    integer :: maxfacs
    integer, allocatable, dimension(:):: spfacs, efacs, lfacs, sgfacs, tgfacs, kxfacs
    integer, dimension(6) :: nfac, sdim
    character(3), dimension(6):: sym

    if (.not.nonlin) return
    write (report_unit, fmt="('yxf sweetspot #proc up to:',i8)") npmax 
    maxfacs=max(nspec,negrid,nlambda,2,2*ntgrid+1,naky,nx)/2+1
    allocate (spfacs(maxfacs),efacs(maxfacs),lfacs(maxfacs),sgfacs(maxfacs),tgfacs(maxfacs),kxfacs(maxfacs),facs(maxfacs,6))
    call factors (nspec, nspfacs, spfacs)
    call factors (negrid, nefacs, efacs)
    call factors (nlambda, nlfacs, lfacs)
    call factors (2, nsgfacs, sgfacs)
    call factors (2*ntgrid+1, ntgfacs, tgfacs)
    call factors (nx, nkxfacs, kxfacs)

    select case (layout)
    case ('lexys','lxyes','lyxes','yxles','xyles')
       sym = (/ "s  ","e  ","l  ","sgn","tg ","x  " /)
       sdim = (/ nspec, negrid, nlambda, 2, 2*ntgrid+1, nx /)
       nfac= (/ nspfacs, nefacs, nlfacs, nsgfacs, ntgfacs, nkxfacs /)
       facs(:,1)=spfacs; facs(:,2)=efacs; facs(:,3)=lfacs
       facs(:,4)=sgfacs; facs(:,5)=tgfacs; facs(:,6)=kxfacs
       call wsweetspots(sym,sdim,nfac,facs,npmax,LUN=report_unit)
    case ('yxels')
       sym = (/ "s  ","l  ","e  ","sgn","tg ","x  " /)
       sdim = (/ nspec, nlambda, negrid, 2, 2*ntgrid+1, nx /)
       nfac= (/ nspfacs, nlfacs, nefacs, nsgfacs, ntgfacs, nkxfacs /)
       facs(:,1)=spfacs; facs(:,2)=lfacs; facs(:,3)=efacs
       facs(:,4)=sgfacs; facs(:,5)=tgfacs; facs(:,6)=kxfacs
       call wsweetspots(sym,sdim,nfac,facs,npmax,LUN=report_unit)
    end select
    deallocate (facs,spfacs,efacs,lfacs,sgfacs,tgfacs,kxfacs)
  end subroutine nprocs_yxf

  !> FIXME : Add documentation  
  subroutine nprocs_e(report_unit)
    use species, only : nspec
    use kt_grids, only: naky, ntheta0
    use le_grids, only: nlambda
    use theta_grid, only: ntgrid
    use gs2_layouts, only: layout, factors
    implicit none
    integer, intent(in) :: report_unit
    integer :: nlfacs, nkxfacs, nkyfacs, nsgfacs, nspfacs, ntgfacs
    integer, dimension(:,:), allocatable :: facs
    integer :: maxfacs
    integer, allocatable, dimension(:):: spfacs, lfacs, sgfacs, tgfacs, kxfacs, kyfacs
    integer, dimension(6) :: nfac, sdim
    character(3), dimension(6):: sym

    write (report_unit, fmt="('#proc sweetspots for e_lo, up to:',i8)") npmax
    maxfacs=max(nspec,nlambda,2,2*ntgrid+1,ntheta0,naky)/2+1
    allocate (spfacs(maxfacs),lfacs(maxfacs),sgfacs(maxfacs),tgfacs(maxfacs),kxfacs(maxfacs),kyfacs(maxfacs),facs(maxfacs,6))
    call factors (nspec, nspfacs, spfacs)
    call factors (nlambda, nlfacs, lfacs)
    call factors (2, nsgfacs, sgfacs)
    call factors (2*ntgrid+1, ntgfacs, tgfacs)
    call factors (ntheta0, nkxfacs, kxfacs)
    call factors (naky, nkyfacs, kyfacs)

    select case (layout)
    case ('lexys','lxyes')
       sym = (/ "s  ","y  ","x  ","l  ","sgn","tg " /)
       sdim = (/ nspec, naky, ntheta0, nlambda, 2, 2*ntgrid+1 /)
       nfac= (/ nspfacs, nkyfacs, nkxfacs, nlfacs, nsgfacs, ntgfacs /)
       facs(:,1)=spfacs; facs(:,2)=kyfacs; facs(:,3)=kxfacs
       facs(:,4)=lfacs; facs(:,5)=sgfacs; facs(:,6)=tgfacs
       call wsweetspots(sym,sdim,nfac,facs,npmax,LUN=report_unit)
    case ('lyxes')
       sym = (/ "s  ","x  ","y  ","l  ","sgn","tg " /)
       sdim = (/ nspec, ntheta0, naky, nlambda, 2, 2*ntgrid+1 /)
       nfac= (/ nspfacs, nkxfacs, nkyfacs, nlfacs, nsgfacs, ntgfacs /)
       facs(:,1)=spfacs; facs(:,2)=kxfacs; facs(:,3)=kyfacs
       facs(:,4)=lfacs; facs(:,5)=sgfacs; facs(:,6)=tgfacs
       call wsweetspots(sym,sdim,nfac,facs,npmax,LUN=report_unit)
    case ('xyles')
       sym = (/ "s  ","l  ","y  ","x  ","sgn","tg " /)
       sdim = (/ nspec, nlambda, naky, ntheta0, 2, 2*ntgrid+1 /)
       nfac= (/ nspfacs, nlfacs, nkyfacs, nkxfacs, nsgfacs, ntgfacs /)
       facs(:,1)=spfacs; facs(:,2)=lfacs; facs(:,3)=kyfacs
       facs(:,4)=kxfacs; facs(:,5)=sgfacs; facs(:,6)=tgfacs
       call wsweetspots(sym,sdim,nfac,facs,npmax,LUN=report_unit)
    case ('yxles','yxels')
       sym = (/ "s  ","l  ","x  ","y  ","sgn","tg " /)
       sdim = (/ nspec, nlambda, ntheta0, naky, 2, 2*ntgrid+1 /)
       nfac= (/ nspfacs, nlfacs, nkxfacs, nkyfacs, nsgfacs, ntgfacs /)
       facs(:,1)=spfacs; facs(:,2)=lfacs; facs(:,3)=kxfacs
       facs(:,4)=kyfacs; facs(:,5)=sgfacs; facs(:,6)=tgfacs
       call wsweetspots(sym,sdim,nfac,facs,npmax,LUN=report_unit)
    end select
    deallocate (facs,spfacs,lfacs,sgfacs,tgfacs,kxfacs,kyfacs)
  end subroutine nprocs_e

  !> FIXME : Add documentation  
  subroutine nprocs_lz(report_unit)
    use species, only : nspec
    use kt_grids, only: naky, ntheta0
    use le_grids, only: negrid
    use theta_grid, only: ntgrid
    use gs2_layouts, only: layout, factors
    implicit none
    integer, intent(in) :: report_unit
    integer :: nefacs, nkxfacs, nkyfacs, nspfacs, ntgfacs
    integer, dimension(:,:), allocatable :: facs
    integer :: maxfacs
    integer, allocatable, dimension(:):: spfacs, efacs, tgfacs, kxfacs, kyfacs
    integer, dimension(5) :: nfac, sdim
    character(3), dimension(5):: sym

    write (report_unit, fmt="('#proc sweetspots for lz_lo, up to:',i8)") npmax

    maxfacs=max(nspec,negrid,2*ntgrid+1,ntheta0,naky)/2+1
    allocate (spfacs(maxfacs),efacs(maxfacs),tgfacs(maxfacs),kxfacs(maxfacs),kyfacs(maxfacs),facs(maxfacs,6))
    call factors (nspec, nspfacs, spfacs)
    call factors (negrid, nefacs, efacs)
    call factors (2*ntgrid+1, ntgfacs, tgfacs)
    call factors (ntheta0, nkxfacs, kxfacs)
    call factors (naky, nkyfacs, kyfacs)

    select case (layout)
    case ('lexys')
       sym = (/ "s  ","y  ","x  ","e  ","tg " /)
       sdim = (/ nspec, naky, ntheta0, negrid, 2*ntgrid+1 /)
       nfac= (/ nspfacs, nkyfacs, nkxfacs, nefacs, ntgfacs /)
       facs(:,1)=spfacs; facs(:,2)=kyfacs; facs(:,3)=kxfacs
       facs(:,4)=efacs; facs(:,5)=tgfacs
       call wsweetspots(sym,sdim,nfac,facs,npmax,LUN=report_unit)
    case ('lyxes','yxles','yxels')
       sym = (/ "s  ","e  ","x  ","y  ","tg " /)
       sdim = (/ nspec, negrid, ntheta0, naky, 2*ntgrid+1 /)
       nfac= (/ nspfacs, nefacs, nkxfacs, nkyfacs, ntgfacs /)
       facs(:,1)=spfacs; facs(:,2)=efacs; facs(:,3)=kxfacs
       facs(:,4)=kyfacs; facs(:,5)=tgfacs
       call wsweetspots(sym,sdim,nfac,facs,npmax,LUN=report_unit)
    case ('xyles','lxyes')
       sym = (/ "s  ","e  ","y  ","x  ","tg " /)
       sdim = (/ nspec, negrid, naky, ntheta0, 2*ntgrid+1 /)
       nfac= (/ nspfacs, nefacs, nkyfacs, nkxfacs, ntgfacs /)
       facs(:,1)=spfacs; facs(:,2)=efacs; facs(:,3)=kyfacs
       facs(:,4)=kxfacs; facs(:,5)=tgfacs
       call wsweetspots(sym,sdim,nfac,facs,npmax,LUN=report_unit)
    end select
    deallocate (spfacs,efacs,tgfacs,kxfacs,kyfacs,facs)
  end subroutine nprocs_lz

  !> FIXME : Add documentation  
  subroutine nprocs_le(report_unit)
    use species, only : nspec
    use kt_grids, only: naky, ntheta0
    use theta_grid, only: ntgrid
    use gs2_layouts, only: layout, factors
    implicit none
    integer, intent(in) :: report_unit
    integer :: nkxfacs, nkyfacs, nspfacs, ntgfacs
    integer, dimension(:,:), allocatable :: facs
    integer :: maxfacs
    integer, allocatable, dimension(:):: spfacs, tgfacs, kxfacs, kyfacs
    integer, dimension(4) :: nfac, sdim
    character(3), dimension(4):: sym

    write (report_unit, fmt="('#proc sweetspots for le_lo, up to:',i8)") npmax

    maxfacs=max(nspec,2*ntgrid+1,ntheta0,naky)/2+1
    allocate (spfacs(maxfacs),tgfacs(maxfacs),kxfacs(maxfacs),kyfacs(maxfacs),facs(maxfacs,4))
    call factors (nspec, nspfacs, spfacs)
    call factors (2*ntgrid+1, ntgfacs, tgfacs)
    call factors (ntheta0, nkxfacs, kxfacs)
    call factors (naky, nkyfacs, kyfacs)

    select case (layout)
    case ('lexys','xyles','lxyes')
       sym = (/ "s  ","y  ","x  ","tg " /)
       sdim = (/ nspec, naky, ntheta0, 2*ntgrid+1 /)
       nfac= (/ nspfacs, nkyfacs, nkxfacs, ntgfacs /)
       facs(:,1)=spfacs; facs(:,2)=kyfacs
       facs(:,3)=kxfacs; facs(:,4)=tgfacs
       call wsweetspots(sym,sdim,nfac,facs,npmax,LUN=report_unit)
    case ('lyxes','yxles','yxels')
       sym = (/ "s  ","x  ","y  ","tg " /)
       sdim = (/ nspec, ntheta0, naky, 2*ntgrid+1 /)
       nfac= (/ nspfacs, nkxfacs, nkyfacs, ntgfacs /)
       facs(:,1)=spfacs; facs(:,2)=kxfacs
       facs(:,3)=kyfacs; facs(:,4)=tgfacs
       call wsweetspots(sym,sdim,nfac,facs,npmax,LUN=report_unit)
    end select
    deallocate (spfacs,tgfacs,kxfacs,kyfacs,facs)
  end subroutine nprocs_le

  !> FIXME : Add documentation
  subroutine nprocs (nmesh, report_unit)
    use nonlinear_terms, only : nonlin
    use species, only : nspec
    use kt_grids, only: naky, ntheta0, nx, ny
    use le_grids, only: negrid, nlambda
    use theta_grid, only: ntgrid
    use gs2_layouts, only: layout, init_x_transform_layouts, init_y_transform_layouts, factors
    use mp, only: nproc, iproc
    implicit none
    integer, intent (in) :: nmesh, report_unit
    integer :: nefacs, nlfacs, nkyfacs, nkxfacs, nspfacs, nkxkyfacs
    integer, dimension(:,:), allocatable :: facs
    integer, dimension(:), allocatable :: mixed_facs, nfac_list, dim_sizes
    integer :: npe, checknpe, i, j, dim_size, nfacs, previous_blocks
    logical :: onlyxoryfac
    character :: layout_char
    write (report_unit, fmt="('Layout = ',a5,/)") layout
    write (report_unit, fmt="('Recommended #proc up to:',i8)") npmax

    ! Initial setup / print headers
    if (nonlin) then
       call init_x_transform_layouts(ntgrid, naky, ntheta0, nlambda, negrid, nspec, nx, nproc, iproc)
       call init_y_transform_layouts(ntgrid, naky, ntheta0, nlambda, negrid, nspec, nx, ny, nproc, iproc)

       ! First write standard sweetspots
       call write_nonlinear_glo_header(report_unit)
    else
       write (report_unit, *)
       write (report_unit, fmt="('Recommended numbers of processors:')")
    end if

    ! Report standard sweet spots
    allocate (facs(max(nspec,naky,ntheta0,negrid,nlambda)/2+1,5))
    allocate (nfac_list(5), dim_sizes(5))
    previous_blocks = 1
    do i = 1, 5
       ! Loop over layout string in reverse
       layout_char = layout(6-i:6-i)
       select case(layout_char)
       case('x')
          dim_size = ntheta0
       case('y')
          dim_size = naky
       case('l')
          dim_size = nlambda
       case('e')
          dim_size = negrid
       case('s')
          dim_size = nspec
       end select
       dim_sizes(i) = dim_size
       call factors (dim_size, nfacs, facs(:,i))
       nfac_list(i) = nfacs
       ! Note special handling for the first call
       if (nonlin) then
          if (i == 1) call write_nonlinear_sweet_spot(layout_char, report_unit, 2, [1, 1], 1, npmax)
          call write_nonlinear_sweet_spot(layout_char,report_unit, nfacs, facs(:, i), &
               previous_blocks, npmax)
       else
          if (i == 1) call write_linear_sweet_spot(layout_char, report_unit, 2, [1, 1], 1, nmesh,ncut)
          call write_linear_sweet_spot(layout_char,report_unit, nfacs, facs(:, i), &
               previous_blocks, nmesh, ncut)
       end if
       previous_blocks = previous_blocks * dim_size
    end do

    ! Now we might want to write sweetspots based on factors of x*y if NL
    if (nonlin) then
       select case (layout)
       case ('yxels', 'yxles', 'xyles')
          write (report_unit, fmt="('------------------------------------------------------------------------------------------------------------------------------------------------')")
          write (report_unit, fmt="(/,'Mixed kx*ky sweetspots (note these are often not recommended:)')")
          call write_nonlinear_glo_header(report_unit)
          allocate (mixed_facs((ntheta0*naky)/2+1))
          call factors (naky*ntheta0, nkxkyfacs, mixed_facs)
          previous_blocks = nlambda * negrid * nspec
          kxky: do i = 2, nkxkyfacs
             npe = mixed_facs(i) * previous_blocks
             if (npe > npmax) exit
             ! Check whether this process count would have been generated by the
             ! plain x or y factorisation or if it is new from the combined x*y
             ! functionality
             onlyxoryfac = .false.
             do j=2, nfac_list(4)
                checknpe = facs(j, 4) * previous_blocks
                if (npe == checknpe) cycle kxky
             end do
             do j = 2, nfac_list(5)
                checknpe = facs(j,5) * dim_sizes(4) * previous_blocks
                if (npe == checknpe) cycle kxky
             end do
             if (.not. onlyxoryfac) call report_idle_procs(npe, 'x*y', report_unit, onlyxoryfac)
          end do kxky
       end select

       write (report_unit, fmt="('------------------------------------------------------------------------------------------------------------------------------------------------')")
       write (report_unit, *)
       write (report_unit, fmt="('(*) denotes process counts that are from factors of the combined kx*ky index rather than the ordered kx or ky indices separately.')")
       write (report_unit, fmt="('To use the unbalanced functionality set unbalanced_xxf = .true. or unbalanced_yxf = .true. in the &layouts_knobs namelist in your GS2 ')")
       write (report_unit, fmt="('input file. You can also set the max_unbalanced_xxf and max_unbalanced_yxf flags in the same namelist in the input file to specify the')")
       write (report_unit, fmt="('maximum amount of computational imbalance allowed. These flags specify the maximum imbalance as 1 with no imbalance as 0, so to allow')")
       write (report_unit, fmt="('50% imbalanced on the xxf decomposition using the following flags in the input file: ')")
       write (report_unit, fmt="('                                                                                     unbalanced_xxf = .true. ')")
       write (report_unit, fmt="('                                                                                     max_unbalanced_xxf = 0.5 ')")
       write (report_unit, fmt="('And likewise for yxf.')")
    end if
    deallocate (facs)
  end subroutine nprocs

  subroutine write_linear_sweet_spot(dim, report_unit, nfacs_dim, facs_dim, &
       previous_dim_block, nmesh, ncut)
    implicit none
    character, intent(in) :: dim
    integer, intent(in) :: report_unit, nfacs_dim, previous_dim_block, nmesh, ncut
    integer, dimension(:), intent(in) :: facs_dim
    integer :: i, npe
    ! We start at 2 to avoid the duplicate factor coming from the previous dimension
    do i = 2, nfacs_dim
       npe = facs_dim(i) * previous_dim_block
       if (nmesh / npe > ncut) write (report_unit, fmt="('  npe = ',i8,'  (',a,')')") npe, dim
    end do
  end subroutine write_linear_sweet_spot

  subroutine write_nonlinear_sweet_spot(dim, report_unit, nfacs_dim, facs_dim, &
       previous_dim_block, npmax)
    implicit none
    character, intent(in) :: dim
    integer, intent(in) :: report_unit, nfacs_dim, previous_dim_block, npmax
    integer, dimension(:), intent(in) :: facs_dim
    integer :: i, npe
    ! We start at 2 to avoid the duplicate factor coming from the previous dimension
    do i = 2, nfacs_dim
       npe = facs_dim(i) * previous_dim_block
       if (npe > npmax) exit
       call report_idle_procs(npe, dim, report_unit)
    end do
  end subroutine write_nonlinear_sweet_spot

  subroutine write_nonlinear_glo_header(report_unit)
    implicit none
    integer, intent(in) :: report_unit
    write (report_unit, *)
    write (report_unit, fmt="('------------------------------------------------------------------------------------------------------------------------------------------------')")
    write (report_unit, fmt="('|   Number of      | Dimension | Percentage data moved | Recommend unbalanced_xxf | xxf unbalanced | Recommend unbalanced_yxf | yxf unbalanced |')")
    write (report_unit, fmt="('|   processes      |   split   | by MPI from xxf->yxf  |    T = true F = false    |   amount (%)   |    T = true F = false    |   amount (%)   |')")
    write (report_unit, fmt="('------------------------------------------------------------------------------------------------------------------------------------------------')")
  end subroutine write_nonlinear_glo_header

  !> This subroutine is used to return the unbalanced suggestions for a given process
  !> count (npe).  The calculation is performed using the calculate_unbalanced_x
  !> and calculate_unbalanced_y subroutines from gs2_layouts.  These return the
  !> unbalanced decomposition size (difference between the suggested small and large
  !> block size for the decomposition) and from this a logical is set to recommend
  !> whether the unbalanced decomposition should be used or not.
  !>
  !> AJ June 2012
  subroutine get_unbalanced_suggestions(npe, percentage_xxf_unbalanced_amount, percentage_yxf_unbalanced_amount, &
       use_unbalanced_xxf, use_unbalanced_yxf)
    use gs2_layouts, only : calculate_unbalanced_x, calculate_unbalanced_y
    implicit none

    integer, intent(in) :: npe
    integer, intent(out) :: percentage_xxf_unbalanced_amount, percentage_yxf_unbalanced_amount
    logical, intent(out) :: use_unbalanced_xxf, use_unbalanced_yxf

    real  :: xxf_unbalanced_amount, yxf_unbalanced_amount

    call calculate_unbalanced_x(npe, 0, xxf_unbalanced_amount)
    if(xxf_unbalanced_amount .gt. 0) then
       use_unbalanced_xxf = .true.
       percentage_xxf_unbalanced_amount = int(xxf_unbalanced_amount * 100.0)
    else
       use_unbalanced_xxf = .false.  
       percentage_xxf_unbalanced_amount = 0
    end if
    call calculate_unbalanced_y(npe, 0, yxf_unbalanced_amount)
    if(yxf_unbalanced_amount .gt. 0) then
       use_unbalanced_yxf = .true.
       percentage_yxf_unbalanced_amount = int(yxf_unbalanced_amount * 100.0)
    else
       use_unbalanced_yxf = .false. 
       percentage_yxf_unbalanced_amount = 0
    end if

  end subroutine get_unbalanced_suggestions

  !> This subroutine is used to return the idle processes from the xxf and yxf layouts.
  !> Idle processes sometimes occur, dependent on the process count used, because
  !> the data domain does not evenly divide by the total number of processes available.
  !> These can cause high communications overheads for the non-linear calculations
  !> for large numbers of processes so it is useful to print this data out in ingen
  !> to let users know which process counts this can happen at for a given input file
  !> We currently use an arbitrary cutoff of 10% difference in idle processes to
  !> suggest that the unbalanced decomposition functionality should be used to
  !> mitigate the impact of this difference.
  !>
  !> AJ June 2012
  subroutine get_idle_processes(npe, idle_percentage, use_unbalanced)
    use gs2_layouts, only : calculate_idle_processes
    implicit none

    integer, intent(in) :: npe
    real, intent(out) :: idle_percentage
    logical, intent(out) :: use_unbalanced

    call calculate_idle_processes(npe, idle_percentage)
    ! This 0.1 represents the arbitrary 10% difference threshold discussed above.
    if(idle_percentage .gt. 0.1) then
       use_unbalanced = .true.
    else
       use_unbalanced = .false.
    end if
    idle_percentage = idle_percentage * 100

  end subroutine get_idle_processes

  !> This subroutine wraps up the output functionality for the code that creates the
  !> list of suggested process counts for the linear computation and checks whether
  !> those process counts work well for the non-linear calculations as well.
  !>
  !> The optional parameter onlyxoryfac (which is a logical) is used to highlight
  !> process counts that have been generated using code which looks for the factors
  !> of kx*ky rather than kx and ky separately.  This is only currently used for
  !> layouts that start with x and y (for instance xyles or yxles).
  !>
  !> AJ June 2012
  subroutine report_idle_procs(npe, distchar, report_unit, onlyxoryfac)
    use optionals, only: get_option_with_default
    implicit none
    integer, intent(in) :: npe, report_unit
    character(*), intent(in) :: distchar
    logical, optional, intent(in) :: onlyxoryfac
    logical :: onlyxoryfac_local
    integer :: percentage_xxf_unbalanced_amount, percentage_yxf_unbalanced_amount
    logical :: use_unbalanced_xxf, use_unbalanced_yxf, use_unbalanced
    real :: idle_percentage
    type(proc_layout_type) :: procl
    character(len = 3) :: xandy_mark
    call get_idle_processes(npe, idle_percentage, use_unbalanced)
    procl%nproc = npe
    procl%percent_xxf_2_yxf = idle_percentage

    onlyxoryfac_local = get_option_with_default(onlyxoryfac, .true.)
    xandy_mark = '   '
    if (.not. onlyxoryfac_local) xandy_mark = '(*)'

    if(use_unbalanced) then
       call get_unbalanced_suggestions(npe, percentage_xxf_unbalanced_amount, percentage_yxf_unbalanced_amount, use_unbalanced_xxf, use_unbalanced_yxf)
       write (report_unit, fmt="('|     ',i8,' ',a3,' |    ',a3,'    |          ',i3,'          |            ',L1,'             |      ',i3,'       |            ',L1,'             |      ',i3,'       |')") & 
            npe, xandy_mark, distchar, INT(idle_percentage), use_unbalanced_xxf, percentage_xxf_unbalanced_amount, use_unbalanced_yxf, percentage_yxf_unbalanced_amount

       procl%should_use_unbalanced_xxf = use_unbalanced_xxf
       procl%percentage_xxf_unbalanced_amount = percentage_xxf_unbalanced_amount
       procl%should_use_unbalanced_yxf = use_unbalanced_yxf
       procl%percentage_yxf_unbalanced_amount = percentage_yxf_unbalanced_amount
    else
       write (report_unit, fmt="('|     ',i8,' ',a3,' |    ',a3,'    |          ',i3,'          |                          |                |                          |                |')") & 
            npe, xandy_mark, distchar, INT(idle_percentage)

       procl%should_use_unbalanced_xxf = .false.
       procl%percentage_xxf_unbalanced_amount = 0.0
       procl%should_use_unbalanced_yxf = .false.
       procl%percentage_yxf_unbalanced_amount = 0.0
    end if

    sweet_spots = [sweet_spots, procl]
  end subroutine report_idle_procs

  !> FIXME : Add documentation
  subroutine tell (a, b, c)
    implicit none
    integer :: j
    character (*), intent (in) :: a
    character (*), intent (in), optional :: b, c
    character (500) :: hack

    j = len_trim(a)

    if (present(b)) then
       j = max(j, len_trim(b))
       if (present(c)) j = max(j, len_trim(c))
    end if

    write (6, *)
    write (6, *)
    hack = repeat('*', j+2)
    write (6, 10) '  '//trim(hack)
    write (6, 10) '   '//trim(a)
    write (interactive_record, 10) trim(a)
    if (present(b)) then
       write (6, 10) '   '//trim(b)
       write (interactive_record, 10) trim(b)
       if (present(c)) then
          write (6, 10) '   '//trim(c)
          write (interactive_record, 10) trim(c)
       end if
    end if
    write (6, 10) '  '//trim(hack)
    write (6, *)
    write (6, *)

10  format (a)

  end subroutine tell

  !> FIXME : Add documentation  
  subroutine text (a)
    use optionals, only: get_option_with_default
    implicit none
    character (*), intent (in), optional :: a
    write(6, '(A)') get_option_with_default(a, '')
  end subroutine text

  !> FIXME : Add documentation  
  subroutine run_number (sel, nbeta)
    implicit none
    integer, intent (out) :: sel
    integer, intent (in) :: nbeta

100 continue
    do 
       call text
       call text ('Lowest run number for this scan (0-999):')
       read (interactive_input, *, err = 100) sel
       if (sel + nbeta - 1 > 999) then
          write (6,*) 'The lowest run number for this scan must be less than ',1000-nbeta
          call try_again
       else if (sel < 0) then
          write (6,*) 'The lowest run number for this scan must be greater than zero.'
          call try_again
       else
          return
       end if
    end do

  end subroutine run_number

  !> FIXME : Add documentation  
  subroutine try_again 
    implicit none
    write (6, *) '*****************'
    write (6, *) 'Please try again.'
    write (6, *) '*****************'

  end subroutine try_again

  !> FIXME : Add documentation  
  subroutine get_choice (in, out)
    implicit none
    integer, intent (in) :: in
    integer, intent (out) :: out

    read (interactive_input, *, err=999) out
    if (out <= 0) out = 0
    if (out > in) out = 0

    return
999 out = 0

  end subroutine get_choice

  !> FIXME : Add documentation  
  subroutine beta_range_low (x, a, x0)
    implicit none
    real, intent (out) :: x
    character (*), intent (in) :: a
    real, intent (in) :: x0

    do
100    continue
       call text
       if (trim(a) == 'le') then
          call text ('Lower limit of beta for this scan (must be > 0):')
       else if(trim(a) == 'lt') then
          call text ('Lower limit of beta for this scan (must be >= 0):')
       end if
       read (interactive_input, *, err=100) x
       if (trim(a) == 'le') then
          if (x <= x0) then
             call try_again
          else
             return
          end if
       else if (trim(a) == 'lt') then
          if (x < x0) then
             call try_again
          else
             return
          end if
       end if
    end do

  end subroutine beta_range_low

  !> FIXME : Add documentation  
  subroutine beta_range_high (x, a, x0)
    implicit none
    real, intent (out) :: x
    character (*), intent (in) :: a
    real, intent (in) :: x0

    do
100    continue
       call text
       call text ('Upper limit of beta for this scan (must be > lower limit):')
       read (interactive_input, *, err=100) x
       if (trim(a) == 'le') then
          if (x <= x0) then
             call try_again
          else
             return
          end if
       end if
    end do

  end subroutine beta_range_high

  !> FIXME : Add documentation  
  subroutine beta_prime_range_low (x)
    implicit none
    real, intent (out) :: x

    do
100    continue
       call tell ('Note: You will be asked to enter:         - d beta/d rho', &
            'Since d beta /d rho <= 0, you should enter a number >= 0')

       call text ('Weakest beta gradient in scan (zero or positive number):')
       read (interactive_input, *, err=100) x
       x = -x
       if (x > 0.) then
          call try_again
       else
          return
       end if
    end do

  end subroutine beta_prime_range_low

  !> FIXME : Add documentation  
  subroutine beta_prime_range_high (x, x0)
    implicit none
    real, intent (out) :: x
    real, intent (in) :: x0

    do
100    continue
       call text
       call text ('Strongest gradient in scan: (positive number)')
       read (interactive_input, *, err=100) x
       x = -x
       if (x >= x0) then
          call try_again
       else
          return
       end if
    end do

  end subroutine beta_prime_range_high

  !> FIXME : Add documentation  
  subroutine num_runs (n)
    implicit none
    integer, intent (out) :: n

    do 
100    continue
       call text
       call text ('How many runs should be done? (n > 1)')
       read (interactive_input, *, err=100) n
       if (n < 2) then
          call try_again
       else
          return
       end if
    end do
  end subroutine num_runs

#include "ingen_auto_gen.inc"  
end module ingen_mod