species.f90 Source File


Contents

Source Code


Source Code

!> FIXME : Add documentation  
module species
  use abstract_config, only: abstract_config_type, CONFIG_MAX_NAME_LEN
  
  implicit none

  private
  public :: init_species, finish_species, wnml_species, check_species
  public :: write_trinity_parameters
  public :: nspec, specie, spec, tracer_species
  public :: f0_maxwellian, f0_tabulated, f0_sdanalytic
  public :: has_ion_species, has_electron_species, has_nonmaxw_species
  public :: has_hybrid_electron_species, is_electron_species, is_hybrid_electron_species
  public :: calculate_f0_arrays, set_current_f0_species, eval_f0
  public :: nonmaxw_corr, f0_values, dlnf0drho, set_overrides
  public :: species_config_type, species_element_config_type
  public :: set_species_config, get_species_config, get_species_element_config

  !> FIXME : Add documentation 
  type :: specie
     real :: z
     real :: mass
     real :: dens, dens0, u0
     real :: tpar0,tperp0
     real :: temp
     real :: tprim
     real :: fprim
     real :: uprim, uprim2
     real :: vnewk
     real :: nu_h ! nu_h controls a hyperviscous term embedded in the collision operator
     real :: stm, zstm, tz, smz, zt
     real :: bess_fac !Artificial factor multiplying the Bessel function argument
     real :: vcrit
     real :: vcprim
     integer :: type
     integer :: f0type
  end type specie

  integer, parameter :: ion_species = 1
  integer, parameter :: electron_species = 2 ! for collision operator
  integer, parameter :: tracer_species = 3 ! for test particle diffusion studies
  integer, parameter :: hybrid_electron_species = 4 ! Adibatic passing, non-adiabatic trapped

  integer, dimension(*), parameter :: electron_like_species = &
       [electron_species, hybrid_electron_species]

  integer, parameter :: f0_maxwellian = 1, &           !< Regular Maxwellian species
                        f0_tabulated = 2, &            !< f0 determined from table in external input file
                        f0_sdanalytic = 3              !< Analytic Gaffey slowing-down distribution

  ! FIXME
  !> Various quantities needed in the calculation of the
  !> slowing-down distribution in case electrons are adiabatic
  real:: me, ZI_fac
  
  !> Defines a working species so that eval_f0 only needs to take one argument (needed for genquad)
  integer:: is_global = -1

  integer :: nspec
  type (specie), dimension (:), allocatable :: spec
  real, dimension (:,:), allocatable :: nonmaxw_corr, f0_values, dlnf0drho

  logical :: initialized = .false.
  logical :: exist

  !> Used to represent the input configuration of species
  type, extends(abstract_config_type) :: species_config_type
     ! namelist : species_knobs
     ! indexed : false
     !> Specifies the electron mass used in calculating the slowing
     !> down parameters in simulations without kinetic electrons.
     real :: me = -1.0
     !> Number of kinetic species to evolve.
     integer :: nspec = 2
     !> Specify the ion charge factor used in calculating the slowing
     !> down parameters. If not specified that will attempt to
     !> calculate from the kinetic ion species.
     real :: zi_fac = -1.0
   contains
     procedure, public :: read => read_species_config
     procedure, public :: write => write_species_config
     procedure, public :: reset => reset_species_config
     procedure, public :: broadcast => broadcast_species_config
     procedure, public, nopass :: get_default_name => get_default_name_species_config
     procedure, public, nopass :: get_default_requires_index => get_default_requires_index_species_config
  end type species_config_type
  
  type(species_config_type) :: species_config

  !> Used to represent the input configuration of species
  type, extends(abstract_config_type) :: species_element_config_type
     ! namelist : species_parameters
     ! indexed : true
     !> Multiplies the argument of the Bessel function for this
     !> species. Useful for exploring the effect of the Bessel
     !> functions.
     real :: bess_fac = 1.0
     !> Set the normalised density for this species
     real :: dens = 1.0
     !> Parameter used for a few specific initial condition types.
     !>
     !> @todo Consider moving this parameter to [[init_g_knobs]].
     real :: dens0 = 1.0
     !> Select which type of background distribution should be assumed
     !> for this species. Must be one of:
     !>
     !> - 'maxwellian' : Standard maxwellian
     !> - 'tabulated' : A general tabulated form
     !> - 'sdanalytic' : A analytic slowing down form.
     !>
     !> See [G. Wilkie
     !> thesis](https://drum.lib.umd.edu/handle/1903/17302).
     character(len = 20) :: f0type = "maxwellian"
     !> Normalised inverse density gradient: \(-\frac{1}{n}\frac{dn}{d
     !> \rho_N}\) (Note here \(\rho_N\) is the normalised radius
     !> \(\frac{\rho}{L_\textrm{ref}}\).
     real :: fprim = 2.2
     !> Normalised mass of this species
     real :: mass = 1.0
     !> Set the species hyper collision frequency used in the pitch
     !> angle scattering collision operator if
     !> [[collision_knobs::hypermult]] is `true`.
     real :: nu_h = 0.0
     !> Set the normalised species temperature
     real :: temp = 1.0
     !> Parameter used for a few specific initial condition types.
     !>
     !> @todo Consider moving this parameter to [[init_g_knobs]].
     real :: tpar0 = 0.
     !> Parameter used for a few specific initial condition types.
     !>
     !> @todo Consider moving this parameter to [[init_g_knobs]].
     real :: tperp0 = 0.
     !> Normalised inverse temperature gradient: \(-\frac{1}{T}\frac{dT}{d
     !> \rho_N}\) (Note here \(\rho_N\) is the normalised radius
     !> \(\frac{\rho}{L_\textrm{ref}}\).
     real :: tprim = 6.9
     !> Sets the characterisation of the species. Should be one of
     !>
     !> - 'ion' : Thermal ion species
     !> - 'default' :  Same as 'ion'
     !> - 'electron' : Thermal electron species
     !> - 'e' : Same as 'electron'
     !> - 'trace' : Tracer species
     !> - 'hybrid_electron' : Adiabatic passing electrons, full
     !>    trapped electron treatment. *NOTE: This option is currently
     !>    considered experimental and is intended for use in
     !>    electrostatic simulations.*
     !>
     !> This primarily controls the treatment in the collision
     !> operator, the detection of adiabatic species. Tracer species
     !> do not contribute to velocity space integrals.
     character(len = 20) :: type = "default"
     !> Parameter used for a few specific initial condition types.
     !>
     !> @todo Consider moving this parameter to [[init_g_knobs]].
     real :: u0 = 1.0
     !> Controls an energy independent part of a source term.
     !>
     !> @note Documentation needs improving, looks like it's an early
     !> attempt at adding in the parallel flow shear drive term
     real :: uprim = 0.0
     !> Controls an energy dependent part of a source term.
     !>
     !> @note Documentation needs improving, looks like it's an early
     !> attempt at adding in the parallel flow shear drive term
     real :: uprim2 = 0.0
     !> Part of the `sdanalytic` slowing down model. See [G. Wilkie
     !> thesis](https://drum.lib.umd.edu/handle/1903/17302).
     real :: vcprim = -999.9
     !> Part of the `sdanalytic` slowing down model. See [G. Wilkie
     !> thesis](https://drum.lib.umd.edu/handle/1903/17302).
     real :: vcrit = -1.0
     !> Sets the normalised species-species collisionality
     !> frequency. For species s, `vnewk = nu_ss Lref/vref` where Lref
     !> is the reference length, `vref` is the reference thermal speed
     !> (not the thermal speed of species s!), \(\nu_{ss}=\sqrt{2}\pi
     !> n_s Z_s^4 e^4 \log(\lambda)/\sqrt{m_s}T^{3/2}_s\) is a
     !> dimensional collision frequency, e is the proton charge,
     !> \(\log(\lambda)\) is the Coloumb logarithm, and
     !> \(Z_s,T_s,n_s,m_s\) are the (charge, temperature, density,
     !> mass) of species s. (The dimensional collision frequency given
     !> here is in Gaussian units. For SI units, include an extra
     !> factor \((4\pi\epsilon_0)^2\) in the denominator of the
     !> definition of \(\nu_{ss}\).
     real :: vnewk = 0.0
     !> Normalised species charge
     real :: z = 1.0
   contains
     procedure, public :: read => read_species_element_config
     procedure, public :: write => write_species_element_config
     procedure, public :: reset => reset_species_element_config
     procedure, public :: broadcast => broadcast_species_element_config
     procedure, public, nopass :: get_default_name => get_default_name_species_element_config
     procedure, public, nopass :: get_default_requires_index => get_default_requires_index_species_element_config
  end type species_element_config_type
  
  type(species_element_config_type), dimension(:), allocatable :: species_element_config

contains
  !> FIXME : Add documentation   
  subroutine check_species(report_unit,beta,tite,alne,dbetadrho_spec)
    implicit none
    integer, intent(in) :: report_unit
    real, intent(in) :: beta, tite
    real, intent(out) :: alne, dbetadrho_spec
    integer :: is
    real :: aln, alp, charge, ee, ne, ptot, zeff_calc
    write (report_unit, fmt="('Number of species: ',i3)") nspec
    zeff_calc = 0.
    charge = 0.
    aln = 0.
    alne = 0.
    alp = 0.
    ptot = 0.
    do is=1, nspec
       write (report_unit, *) 
       write (report_unit, fmt="('  Species ',i3)") is
       if (spec(is)%type == ion_species) write (report_unit, fmt="('    Type:             Ion')")
       if (spec(is)%type == electron_species) write (report_unit, fmt="('    Type:             Electron')")
       if (spec(is)%type == tracer_species) write (report_unit, fmt="('    Type:             Trace')")
       if (spec(is)%type == hybrid_electron_species) write (report_unit, fmt="('    Type:             Hybrid-Electron')")       
       write (report_unit, fmt="('    Charge:         ',f7.3)") spec(is)%z
       write (report_unit, fmt="('    Mass:             ',es11.4)") spec(is)%mass
       write (report_unit, fmt="('    Density:        ',f7.3)") spec(is)%dens
       write (report_unit, fmt="('    Temperature:    ',f7.3)") spec(is)%temp
       write (report_unit, fmt="('    Collisionality:   ',es11.4)") spec(is)%vnewk
       write (report_unit, fmt="('    Normalized Inverse Gradient Scale Lengths:')")
       write (report_unit, fmt="('      Temperature:  ',f7.3)") spec(is)%tprim
       write (report_unit, fmt="('      Density:      ',f7.3)") spec(is)%fprim
       write (report_unit, fmt="('      Parallel v:   ',f7.3)") spec(is)%uprim
       if (spec(is)%bess_fac.ne.1.0) &
            write (report_unit, fmt="('      Bessel function arg multiplied by:   ',f7.3)") spec(is)%bess_fac
       !        write (report_unit, fmt="('    Ignore this:')")
       !        write (report_unit, fmt="('    D_0: ',es10.4)") spec(is)%dens0
       if (spec(is)%type /= 2) then
          zeff_calc = zeff_calc + spec(is)%dens*spec(is)%z**2
          charge = charge + spec(is)%dens*spec(is)%z
          aln = aln + spec(is)%dens*spec(is)%z*spec(is)%fprim
       else
          alne = alne + spec(is)%dens*spec(is)%z*spec(is)%fprim
          ne = spec(is)%dens
          ee = spec(is)%z
       end if
       alp = alp + spec(is)%dens * spec(is)%temp *(spec(is)%fprim + spec(is)%tprim)
       ptot = ptot + spec(is)%dens * spec(is)%temp
    end do
    
    if (.not. has_electron_species(spec)) then
       ptot = ptot + 1./tite   ! electron contribution to pressure
       alp = alp + aln/tite    ! assuming charge neutrality, electron contribution to alp
    end if
    
    alp = alp / ptot
    
    write (report_unit, *) 
    write (report_unit, fmt="('------------------------------------------------------------')")
    
    write (report_unit, fmt="('Calculated Z_eff: ',f7.3)") zeff_calc
    
    if (has_electron_species(spec)) then
       if (abs(charge+ne*ee) > 1.e-2) then
          if (charge+ne*ee < 0.) then
             write (report_unit, *) 
             write (report_unit, fmt="('################# WARNING #######################')")
             write (report_unit, fmt="('You are neglecting an ion species.')")
             write (report_unit, fmt="('This species has a charge fraction of ',f7.3)") abs(charge+ne*ee)
             write (report_unit, &
                  & fmt="('and a normalized inverse density gradient scale length of ',f7.3)") &
                  (aln+alne)/(charge+ne*ee)
             write (report_unit, fmt="('################# WARNING #######################')")
          else
             write (report_unit, *) 
             write (report_unit, fmt="('################# WARNING #######################')")
             write (report_unit, fmt="('There is an excess ion charge fraction of ',f7.3)") abs(charge+ne*ee)
             write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") 
             write (report_unit, fmt="('################# WARNING #######################')")
          end if
       else
          if (abs(aln+alne) > 1.e-2) then
             write (report_unit, *) 
             write (report_unit, fmt="('################# WARNING #######################')")
             write (report_unit, fmt="('The density gradients are inconsistent'/' a/lni =',e12.4,' but alne =',e12.4)") aln, alne
             write (report_unit, fmt="('################# WARNING #######################')")
          end if
       end if
    else
       if (charge > 1.01) then
          write (report_unit, *) 
          write (report_unit, fmt="('################# WARNING #######################')")
          write (report_unit, fmt="('There is an excess ion charge fraction of ',f7.3)") charge-1.
          write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") 
          write (report_unit, fmt="('################# WARNING #######################')")
       end if
       if (charge < 0.99) then
          write (report_unit, *) 
          write (report_unit, fmt="('################# WARNING #######################')")
          write (report_unit, fmt="('You are neglecting an ion species.')")
          write (report_unit, fmt="('This species has a charge fraction of ',f7.3)") abs(charge-1.)
          write (report_unit, fmt="('################# WARNING #######################')")
       end if
    end if

    if (has_hybrid_electron_species(spec)) then
       write (report_unit, *)
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, fmt="('You are using a hybrid electron species.')")
       write (report_unit, fmt="('This is currently considered experimental.')")
       write (report_unit, fmt="('################# WARNING #######################')")
    end if

    write (report_unit, *)
    write (report_unit, fmt="('------------------------------------------------------------')")
    
    write (report_unit, *) 
    write (report_unit, fmt="('GS2 beta parameter = ',f9.4)") beta
    write (report_unit, fmt="('Total beta = ',f9.4)") beta*ptot
    write (report_unit, *) 
    write (report_unit, fmt="('The total normalized inverse pressure gradient scale length is ',f10.4)") alp
    dbetadrho_spec = -beta*ptot*alp
    write (report_unit, fmt="('corresponding to d beta / d rho = ',f10.4)") dbetadrho_spec
  end subroutine check_species

  !> FIXME : Add documentation   
  subroutine wnml_species(unit)
    implicit none
    integer, intent(in) :: unit
    integer :: i
    character (100) :: line
    if (.not. exist) return
    write (unit, *)
    write (unit, fmt="(' &',a)") "species_knobs"
    write (unit, fmt="(' nspec = ',i2)") nspec
    write (unit, fmt="(' /')")

    do i=1,nspec
       write (unit, *)
       write (line, *) i
       write (unit, fmt="(' &',a)") &
            & trim("species_parameters_"//trim(adjustl(line)))
       write (unit, fmt="(' z = ',e13.6)") spec(i)%z
       write (unit, fmt="(' mass = ',e13.6)") spec(i)%mass
       write (unit, fmt="(' dens = ',e13.6)") spec(i)%dens
       write (unit, fmt="(' temp = ',e13.6)") spec(i)%temp
       write (unit, fmt="(' tprim = ',e13.6)") spec(i)%tprim
       write (unit, fmt="(' fprim = ',e13.6)") spec(i)%fprim
       write (unit, fmt="(' uprim = ',e13.6)") spec(i)%uprim
       write (unit, fmt="(' vcrit = ',e13.6)") spec(i)%vcrit
       write (unit, fmt="(' vcprim = ',e13.6)") spec(i)%vcprim
       if (spec(i)%uprim2 /= 0.) write (unit, fmt="(' uprim2 = ',e13.6)") spec(i)%uprim2
       write (unit, fmt="(' vnewk = ',e13.6)") spec(i)%vnewk
       if (spec(i)%type == ion_species) &
            write (unit, fmt="(a)") ' type = "ion" /'
       if (spec(i)%type == electron_species) &
            write (unit, fmt="(a)") ' type = "electron"  /'
       if (spec(i)%type == hybrid_electron_species) &
            write (unit, fmt="(a)") ' type = "hybrid_electron"  /'
       if (spec(i)%type == tracer_species) &
            write (unit, fmt="(a)") ' type = "trace"  /'
       if (spec(i)%f0type == f0_maxwellian) &
            write (unit, fmt="(a)") ' f0type = "maxwellian" /'
       if (spec(i)%f0type == f0_sdanalytic) &
            write (unit, fmt="(a)") ' f0type = "sdanalytic"  /'
       if (spec(i)%f0type == f0_tabulated) &
            write (unit, fmt="(a)") ' f0type = "tabulated"  /'
       write (unit, fmt="(' dens0 = ',e13.6)") spec(i)%dens0
       if(spec(i)%bess_fac.ne.1.0) &
            write (unit, fmt="(' bess_fac = ',e13.6)") spec(i)%bess_fac
    end do
  end subroutine wnml_species

  !> FIXME : Add documentation   
  subroutine init_species(species_config_in, species_elements_config_in)
    implicit none
    type(species_config_type), intent(in), optional :: species_config_in
    type(species_element_config_type), intent(in), dimension(:), allocatable, optional :: species_elements_config_in

    if (initialized) return
    initialized = .true.

    call read_parameters(species_config_in, species_elements_config_in)

  end subroutine init_species

  !> FIXME : Add documentation   
  subroutine read_parameters(species_config_in, species_elements_config_in)
    use file_utils, only: input_unit, error_unit, get_indexed_namelist_unit, input_unit_exist
    use text_options, only: text_option, get_option_value
    use mp, only: proc0, mp_abort
    implicit none
    type(species_config_type), intent(in), optional :: species_config_in
    type(species_element_config_type), intent(in), dimension(:), allocatable, optional :: species_elements_config_in
    character(20) :: type, f0type
    integer :: is, ierr

    type (text_option), dimension (*), parameter :: typeopts = [ &
         text_option('default', ion_species), &
         text_option('ion', ion_species), &
         text_option('electron', electron_species), &
         text_option('e', electron_species), &
         text_option('trace', tracer_species), &
         text_option('hybrid_electron', hybrid_electron_species) &
         ]

    type (text_option), dimension (3), parameter :: f0_opts = [ &
         text_option('maxwellian', f0_maxwellian), &
         text_option('tabulated', f0_tabulated), &
         text_option('sdanalytic', f0_sdanalytic) &
         ]

    if (present(species_config_in)) species_config = species_config_in

    call species_config%init(name = 'species_knobs', requires_index = .false.)

    ! Copy out internal values into module level parameters
    me = species_config%me
    nspec = species_config%nspec
    zi_fac = species_config%zi_fac

    exist = species_config%exist
    
    if (proc0) then
       if (nspec < 1) then
          ierr = error_unit()
          write (unit=ierr, &
               fmt="('Invalid nspec in species_knobs: ', i5)") nspec
          call mp_abort('Invalid nspec in species_knobs')
       end if
    end if

    allocate (spec(nspec))

    if (present(species_elements_config_in)) species_element_config = species_elements_config_in

    if (.not.allocated(species_element_config)) allocate(species_element_config(nspec))

    if (size(species_element_config) .ne. nspec) then
       if (proc0) print*,"inconsistent number of config elements"
    endif
       
    ierr = error_unit()
    do is = 1, nspec
       call species_element_config(is)%init(name = 'species_parameters', requires_index = .true., index = is)

       ! Copy out internal values into module level parameters
       spec(is)%bess_fac = species_element_config(is)%bess_fac
       spec(is)%dens = species_element_config(is)%dens
       spec(is)%dens0 = species_element_config(is)%dens0
       spec(is)%fprim = species_element_config(is)%fprim
       spec(is)%mass = species_element_config(is)%mass
       spec(is)%nu_h = species_element_config(is)%nu_h
       spec(is)%temp = species_element_config(is)%temp
       spec(is)%tpar0 = species_element_config(is)%tpar0
       spec(is)%tperp0 = species_element_config(is)%tperp0
       spec(is)%tprim = species_element_config(is)%tprim
       spec(is)%u0 = species_element_config(is)%u0
       spec(is)%uprim = species_element_config(is)%uprim
       spec(is)%uprim2 = species_element_config(is)%uprim2
       spec(is)%vcprim = species_element_config(is)%vcprim
       spec(is)%vcrit = species_element_config(is)%vcrit
       spec(is)%vnewk = species_element_config(is)%vnewk
       spec(is)%z = species_element_config(is)%z

       call get_option_value (species_element_config(is)%type, typeopts, spec(is)%type, ierr, "type in species_parameters_x",.true.)
       call get_option_value (species_element_config(is)%f0type, f0_opts, spec(is)%f0type, ierr, "f0type in species_parameters_x",.true.)
    end do

    spec%stm = sqrt(spec%temp / spec%mass)
    spec%zstm = spec%z / sqrt(spec%temp * spec%mass)
    spec%tz = spec%temp / spec%z
    spec%zt = spec%z / spec%temp
    spec%smz = abs(sqrt(spec%temp * spec%mass) / spec%z)
  end subroutine read_parameters

  !> FIXME : Add documentation
  subroutine calculate_slowingdown_parameters(alpha_index,vc,vcprim)
    use constants, only: pi
    use mp, only: proc0
    implicit none
    integer, intent(in):: alpha_index
    real, intent(inout):: vc, vcprim
    integer:: electron_spec, main_ion_species, is,js
    real:: vinj, Te_prim, ZI_fac_local, ne_prim, ni_prim, ne, vte
    
    electron_spec = -1
    main_ion_species = -1

    do is = 1,nspec
       if (is_electron_species(spec(is))) electron_spec = is
       if (main_ion_species < 1 .and. spec(is)%type .eq. ion_species) main_ion_species = is
    end do

    if (electron_spec .GT. 0) then
       me = spec(electron_spec)%mass
       ne = spec(electron_spec)%dens
       ne_prim = spec(electron_spec)%fprim
       Te_prim = spec(electron_spec)%tprim
       vte = sqrt(spec(electron_spec)%temp/spec(electron_spec)%mass)
    else
       if (me .LT. 0.0) then
          if (proc0) write(*,*) "WARNING: vte not specified in species_knobs. Assuming reference species is deuterium, so me = 0.0002723085"
          me = 0.0002723085
       end if             

       ne = sum(spec(:)%z*spec(:)%dens)
       ne_prim = sum(spec(:)%z*spec(:)%dens*spec(:)%fprim)/ne
       Te_prim = spec(main_ion_species)%tprim
       vte = sqrt(spec(main_ion_species)%temp/spec(main_ion_species)%mass)
       if (proc0) then
          write(*,*) "Since electrons are adiabatic, we need to improvise electron parameters to "
          write(*,*) "caluclate alpha species properties. Imposed by global quasineutrality:"
          write(*,*) "  ne = ", ne
          write(*,*) "  vte = ", vte
          write(*,*) "  ne_prim = ", ne_prim
          write(*,*) "  Te_prim = ", Te_prim
       end if
    end if

    ni_prim = 0.0
    ! If ZI_fac is not specified, calculate it from existing ion species:
    if (ZI_fac .LT. 0.0) then
       ZI_fac_local = 0.0
       do js = 1,nspec
          if (spec(js)%type .eq. ion_species) then 
             ZI_fac_local = ZI_fac_local   + spec(alpha_index)%mass*spec(js)%dens*spec(js)%z**2/(spec(js)%mass*ne)
             ni_prim = ni_prim + spec(alpha_index)%mass*spec(js)%dens*spec(js)%z**2*spec(js)%fprim/(spec(js)%mass*ne)
          end if
       end do
    else
       ZI_fac_local = ZI_fac
    end if
    ni_prim = ni_prim / ZI_fac_local

    ! Calculate vc/vinj
    vinj = sqrt(spec(alpha_index)%temp/spec(alpha_index)%mass)
    vc = (0.25*3.0*sqrt(pi)*ZI_fac_local*me/spec(alpha_index)%mass)**(1.0/3.0)*vte/vinj
    vcprim = ni_prim - ne_prim + 1.5*Te_prim

  end subroutine calculate_slowingdown_parameters

  !> FIXME : Add documentation   
  subroutine calculate_f0_arrays(egrid)
     implicit none
     real, dimension (:,:), intent (in) :: egrid
     integer:: negrid,is,ie
     real:: v, f0, df0dE, f0prim
    
     negrid = size(egrid(:,1))

     allocate(nonmaxw_corr(negrid,nspec))
     allocate(f0_values(negrid,nspec))
     allocate(dlnf0drho(negrid,nspec))

     do is = 1,nspec
        is_global = is
        do ie = 1,negrid
           v = sqrt(egrid(ie,is))

           select case (spec(is)%f0type)
           case (f0_maxwellian)
              f0_values(ie,is) = eval_f0(v)
              nonmaxw_corr(ie,is) = 1.0
              dlnf0drho(ie,is) = -spec(is)%fprim - (egrid(ie,is) - 1.5)*spec(is)%tprim
           case (f0_sdanalytic)
              call eval_f0_sdanalytic(is,v,f0,df0dE,f0prim)
              f0_values(ie,is) = f0
              nonmaxw_corr(ie,is) = -spec(is)%temp*df0dE
              dlnf0drho(ie,is) = f0prim
           end select
        end do
        if (spec(is)%f0type .EQ. f0_tabulated) then
           call calculate_f0_arrays_tabulated(is,egrid)
        end if
        
     end do

  end subroutine calculate_f0_arrays

  !> FIXME : Add documentation   
  subroutine set_current_f0_species(is)
     implicit none
     integer, intent(in):: is
     is_global = is
  end subroutine set_current_f0_species

  !> FIXME : Add documentation
  !!
  !! @note Relies upon is_global to determine species. This function only takes one argument so that
  !! it can be passed as an argument to genquad
  function eval_f0(v)
     use constants, only: pi
     use mp, only: mp_abort
     implicit none
     real, intent(in):: v
     real:: eval_f0
     real::f0temp, dummy1,dummy2

     if (is_global .LT. 0)  call mp_abort("is_global not specified when calling eval_f0")

     select case(spec(is_global)%f0type)
     case (f0_maxwellian)
        f0temp = exp(-v**2)/pi**1.5        
     case (f0_sdanalytic)
        call eval_f0_sdanalytic(is_global,v,f0temp,dummy1,dummy2)
     case (f0_tabulated)
        call mp_abort("ERROR: eval_f0 cannot yet handle tabulated option.")
     end select

     eval_f0 = f0temp

  end function eval_f0

  !> FIXME : Add documentation   
  subroutine eval_f0_sdanalytic(is,v,f0,df0dE_out,f0prim)
    use constants, only: pi
    use mp, only: mp_abort
    implicit none
    integer,intent(in):: is
    real,intent(in):: v
    real,intent(inout):: f0,df0dE_out, f0prim
    real:: A, df0dv,vcprim, vcrit

    if (spec(is)%vcrit .LE. 0.0) then
       call calculate_slowingdown_parameters(is,vcrit,vcprim)
       if (abs(spec(is)%vcprim + 999.0) .LT. 1.0e-3)  vcprim = spec(is)%vcprim
    else
       !> If vcrit is specified, there is no way of knowing its radial derivative, so set dvc/drho = 0
       vcrit = spec(is)%vcrit
       if (abs(spec(is)%vcprim + 999.0) .LT. 1.0e-3)   vcprim = 0.0
    end if


    A = (3.0/(4.0*pi))/log( (vcrit**3 + 1.0)/(vcrit**3))

    f0 = A/(vcrit**3 + v**3)
    df0dv = -A*3.0*v**2/(vcrit**3 + v**3)**2
    df0dE_out =  (0.5/v)*df0dv / ( f0 * spec(is)%temp )
    f0prim = -spec(is)%fprim - (-(vcrit**3/(vcrit**3 + v**3)) + (1.0/( log(1.0 + vcrit**(-3)) * (1.0 + vcrit**(3)))))*vcprim

    if (v .GT. 1.0) then
       f0 = 0.0
       f0prim = 0.0
       df0dE_out = 1.0
    end if

!    write(*,*) v, vcrit,vcprim,f0prim

  end subroutine eval_f0_sdanalytic

  !> This subroutine calculates f0 on the grid from an external
  !! input file. The grid on the input file can differ from that 
  !! of gs2. A cubic spline is used to interpolate between the two.
  !! The user can either specify df0/dE or it can be estimated internally.
  !! The radial derivative of f0 as a function of energy should also be given.
  subroutine calculate_f0_arrays_tabulated(is,egrid)
    !
    ! The egrid is already given by get_legendre_grids. 
    !
    ! The input file can take two forms, as controlled by the
    ! control parameter mode, the first integer read from the file
    !  mode = 1: E points in first columns, F0(E) in second. To calculate
    !            generalized temperature, a cubic spline is used to 
    !            interpolate and the slope is taken from that.
    !       = 2: First column is E, second is F0(E), the third is (1/F0)*dF0/dE
    !       = 3: E, F0(E), -(1/F0)*dF0/drho 
    !       = 4: E, F0(E),dF0/dE, -(1/F0)*dF0/drho 
    use mp, only: broadcast
    use constants, only: pi
    use file_utils, only: run_name
    use splines, only: fitp_curvd, fitp_curv1, fitp_curv2
    implicit none
    integer, intent(in) :: is
    real, dimension(:,:),intent(in):: egrid
    integer:: f0in_unit = 21, mode, ie, ierr, numdat, negrid
    logical:: acceptable_spline
    real:: moment0, sigma, sigma_fac, df0dE
    real, dimension(:), allocatable:: f0_values_dat, df0dE_dat, egrid_dat, &
                                      yp, temp, df0drho_dat
    
    ! Open file and read column option
    open(unit=f0in_unit,file=trim(run_name)//'.f0in',status='old',action='read')
    read(f0in_unit,*) mode
    read(f0in_unit,*) numdat

    allocate(f0_values_dat(numdat))
    allocate(df0dE_dat(numdat))
    allocate(df0drho_dat(numdat))
    allocate(egrid_dat(numdat))
    allocate(yp(numdat))
    allocate(temp(numdat))

    negrid = size(egrid(:,1))
    
    sigma = 1.0
    sigma_fac = 2.0
    acceptable_spline = .false.

    if (mode .EQ. 1) then
       ! Read f0 values
       do ie = 1,numdat
          read(f0in_unit,*) egrid_dat(ie), f0_values_dat(ie)
       end do
       close(f0in_unit)

       ! Perform cubic spline to get F0 and its slope at grid points
       ! If F0 ends up being negative or zero anywhere on the grid, double the tension
       ! factor and try again

       do while (.NOT. acceptable_spline)
          ! Generate spline parameters
          call fitp_curv1(numdat,egrid_dat,f0_values_dat,0.0,0.0,3,yp,temp,sigma,ierr)
          if (ierr .NE. 0) then
             write(*,*) "fitp_curv1 returned error code ", ierr
             stop 1
          end if

          do ie = 1,negrid
             ! Interpolate to get f0 at grid points
             f0_values(ie,is) = fitp_curv2(egrid(ie,is),numdat,egrid_dat, &
                             f0_values_dat,yp,sigma)

             ! Calculate d/dE lnF0 to get generalised temperature
             df0dE =  fitp_curvd(egrid(ie,is),numdat,egrid_dat, f0_values_dat,yp,sigma)  / &
                                (f0_values(ie,is) * spec(is)%temp)
             nonmaxw_corr(ie,is) = -spec(is)%temp*df0dE

          end do

          if (minval(f0_values(:,is)) .GT. epsilon(0.0)) then 
             acceptable_spline = .true.
          else
             sigma = sigma * sigma_fac
          end if
       end do

       dlnf0drho(:,is) = - spec(is)%fprim

    else if (mode .EQ. 2) then
       ! Read both f0 and df0/dE
       do ie = 1,numdat
          read(f0in_unit,*) egrid_dat(ie), f0_values_dat(ie), df0dE_dat(ie)
       end do
       close(f0in_unit)

       do while (.NOT. acceptable_spline)
          ! Generate spline parameters for f0
          call fitp_curv1(numdat,egrid_dat,f0_values_dat,0.0,0.0,3,yp,temp,sigma,ierr)

          if (ierr .NE. 0) then
             write(*,*) "fitp_curv1 returned error code ", ierr
             stop 1
          end if

          do ie = 1,negrid

             ! Interpolate to get F0 at grid points
             f0_values(ie,is) = fitp_curv2(egrid(ie,is),numdat,egrid_dat, &
                                f0_values_dat,yp,sigma)
          end do

          ! Generate spline parameters for df0/dE
          call fitp_curv1(numdat,egrid_dat,df0dE_dat,0.0,0.0,3,yp,temp,sigma,ierr)
          if (ierr .NE. 0) then
             write(*,*) "fitp_curv1 returned error code ", ierr
             stop 1
          end if

          do ie = 1,negrid
             ! Interpolate to get f0 at grid points
             df0dE     = fitp_curv2(egrid(ie,is),numdat,egrid_dat, &
                                df0dE_dat,yp,sigma)  / spec(is)%temp
             nonmaxw_corr(ie,is) = -spec(is)%temp*df0dE
          end do
 
          if (minval(f0_values(:,is)) .GT. epsilon(0.0)) then 
             acceptable_spline = .true.
          else
             sigma = sigma * sigma_fac
          end if

       end do

       dlnf0drho(:,is) = - spec(is)%fprim

    else if (mode .EQ. 3) then
       ! Read f0 and df0/drho
       do ie = 1,numdat
          read(f0in_unit,*) egrid_dat(ie), f0_values_dat(ie), df0drho_dat(ie)
       end do
       close(f0in_unit)

       do while (.NOT. acceptable_spline)
          ! Generate spline parameters for f0
          call fitp_curv1(numdat,egrid_dat,f0_values_dat,0.0,0.0,3,yp,temp,sigma,ierr)

          if (ierr .NE. 0) then
             write(*,*) "fitp_curv1 returned error code ", ierr
             stop 1
          end if

          do ie = 1,negrid
             ! Interpolate to get F0 at grid points
             f0_values(ie,is) = fitp_curv2(egrid(ie,is),numdat,egrid_dat, &
                                 f0_values_dat,yp,sigma)
             ! Calculate d/dE lnF0 to get generalised temperature
             df0dE = fitp_curvd(egrid(ie,is),numdat,egrid_dat, &
                             f0_values_dat,yp,sigma) / (f0_values(ie,is) * spec(is)%temp)
             nonmaxw_corr(ie,is) = -spec(is)%temp*df0dE
          end do

          ! Generate spline parameters for df0/drho
          call fitp_curv1(numdat,egrid_dat,df0drho_dat,0.0,0.0,3,yp,temp,sigma,ierr)
          if (ierr .NE. 0) then
             write(*,*) "fitp_curv1 returned error code ", ierr
             stop 1
          end if

          do ie = 1,negrid
             ! Interpolate to get f0prim at grid points
             dlnf0drho(ie,is)    = fitp_curv2(egrid(ie,is),numdat,egrid_dat, &
                                df0drho_dat,yp,sigma)
          end do
 
          if (minval(f0_values(:,is)) .GT. epsilon(0.0)) then 
             acceptable_spline = .true.
          else
             sigma = sigma * sigma_fac
          end if

       end do

    else if (mode .EQ. 4) then
       ! Read f0, df0/dE, df0/drho
       do ie = 1,numdat
          read(f0in_unit,*) egrid_dat(ie), f0_values_dat(ie), df0dE_dat(ie), df0drho_dat(ie)
       end do
       close(f0in_unit)

       do while (.NOT. acceptable_spline)
          ! Generate spline parameters for f0
          call fitp_curv1(numdat,egrid_dat,f0_values_dat,0.0,0.0,3,yp,temp,sigma,ierr)

          if (ierr .NE. 0) then
             write(*,*) "fitp_curv1 returned error code ", ierr
             stop 1
          end if

          do ie = 1,negrid

             ! Interpolate to get F0 at grid points
             f0_values(ie,is) = fitp_curv2(egrid(ie,is),numdat,egrid_dat, &
                                f0_values_dat,yp,sigma)
          end do

          ! Generate spline parameters for df0/dE
          call fitp_curv1(numdat,egrid_dat,df0dE_dat,0.0,0.0,3,yp,temp,sigma,ierr)
          if (ierr .NE. 0) then
             write(*,*) "fitp_curv1 returned error code ", ierr
             stop 1
          end if

          do ie = 1,negrid
             ! Interpolate to get f0 at grid points
             df0dE     = fitp_curv2(egrid(ie,is),numdat,egrid_dat, &
                                df0dE_dat,yp,sigma)  / spec(is)%temp
             nonmaxw_corr(ie,is) = -spec(is)%temp*df0dE
          end do
 
          ! Generate spline parameters for df0/drho
          call fitp_curv1(numdat,egrid_dat,df0drho_dat,0.0,0.0,3,yp,temp,sigma,ierr)
          if (ierr .NE. 0) then
             write(*,*) "fitp_curv1 returned error code ", ierr
             stop 1
          end if

          do ie = 1,negrid
             ! Interpolate to get f0 at grid points
             dlnf0drho(ie,is)     = fitp_curv2(egrid(ie,is),numdat,egrid_dat, &
                                 df0drho_dat,yp,sigma) 
          end do
 
          if (minval(f0_values(:,is)) .GT. epsilon(0.0)) then 
             acceptable_spline = .true.
          else
             sigma = sigma * sigma_fac
          end if

       end do

       dlnf0drho(:,is) = - spec(is)%fprim

    else
       write(*,*) "ERROR. First line in f0 input file should be mode" 
       stop 1
    end if

    ! Now calculate moments of F0 from trapezoidal rule.
    ! Should probably do this by quadrature, 
    ! but le_grids is not defined yet

    moment0 = 0.0
    do ie = 1,numdat-1
       moment0  = moment0 + 0.5*(egrid_dat(ie+1)-egrid_dat(ie)) * &
                           ( sqrt(egrid_dat(ie))*f0_values_dat(ie) + &
                             sqrt(egrid_dat(ie+1))*f0_values_dat(ie+1) )
    end do
    moment0 = moment0*4.0*pi

    
    deallocate(egrid_dat,f0_values_dat,df0dE_dat,df0drho_dat,yp,temp)

  end subroutine calculate_f0_arrays_tabulated

  !> Determine if any of the known species are electron like
  !>
  !> @note This should be effectively constant, could cache during init
  pure logical function has_electron_species (spec)
    implicit none
    type (specie), dimension (:), intent (in) :: spec
    has_electron_species = any(is_electron_species(spec))
  end function has_electron_species

  !> Determine if any of the known species are hybrid electrons
  !>
  !> @note This should be effectively constant, could cache during init
  pure logical function has_hybrid_electron_species (spec)
    implicit none
    type (specie), dimension (:), intent (in) :: spec
    has_hybrid_electron_species = any(is_hybrid_electron_species(spec))
  end function has_hybrid_electron_species

  !> Determine if any of the known species are ion like
  !>
  !> @note This should be effectively constant, could cache during init
  pure logical function has_ion_species (spec)
    implicit none
    type (specie), dimension (:), intent (in) :: spec
    has_ion_species = any(is_ion_species(spec))
  end function has_ion_species

  !> Determine if the passed species instance is considered electron
  !> like. In other words if spec%type matches any of the recognised
  !> electron like types (electron, hybrid_electron)
  elemental logical function is_electron_species(spec)
    implicit none
    type(specie), intent(in) :: spec
    is_electron_species = any(spec%type == electron_like_species)
  end function is_electron_species

  !> Determine if the passed species instance corresponds to hybrid electrons
  elemental logical function is_hybrid_electron_species(spec)
    implicit none
    type(specie), intent(in) :: spec
    is_hybrid_electron_species = spec%type == hybrid_electron_species
  end function is_hybrid_electron_species

  !> Determine if the passed species instance corresponds to ions
  elemental logical function is_ion_species(spec)
    implicit none
    type(specie), intent(in) :: spec
    is_ion_species = spec%type == ion_species
  end function is_ion_species

  !> Determine if any of the known species are non-Maxwellian
  pure logical function has_nonmaxw_species (spec)
    implicit none
    type (specie), dimension (:), intent (in) :: spec
    has_nonmaxw_species = .NOT. all(spec%f0type == f0_maxwellian)
  end function has_nonmaxw_species

  !> FIXME : Add documentation
  subroutine finish_species
    implicit none
    if(allocated(spec)) deallocate (spec)
    !    if(allocated(f0_values)) deallocate (f0_values,nonmaxw_corr,dlnf0drho)
    call species_config%reset()
    if (allocated(species_element_config)) deallocate(species_element_config)
    initialized = .false.
  end subroutine finish_species

  !> FIXME : Add documentation   
  subroutine set_overrides(prof_ov)
    use overrides, only: profiles_overrides_type
    use unit_tests, only: debug_message
    type(profiles_overrides_type), intent(in) :: prof_ov
    integer :: is
    integer, parameter :: verbosity=3
    character(len=800) :: message
    call debug_message(verbosity, 'species::set_overrides starting')
    do is = 1,nspec
      message = ''
      write(message,*) 'species::set_overrides setting species ', is
      call debug_message(verbosity, trim(message))
      if (prof_ov%override_tprim(is)) spec(is)%tprim = prof_ov%tprim(is)
      if (prof_ov%override_fprim(is)) spec(is)%fprim = prof_ov%fprim(is)
      if (prof_ov%override_temp(is)) spec(is)%temp = prof_ov%temp(is)
      if (prof_ov%override_dens(is)) spec(is)%dens = prof_ov%dens(is)
      if (prof_ov%override_vnewk(is)) spec(is)%vnewk = prof_ov%vnewk(is)
    end do
    call debug_message(verbosity,&
    'species::set_overrides calling calculate_and_broadcast_species_properties')
    call calculate_and_broadcast_species_properties
    call debug_message(verbosity, 'species::set_overrides finishing')
  end subroutine set_overrides

  !> FIXME : Add documentation 
  subroutine calculate_and_broadcast_species_properties
    use mp, only: proc0, broadcast
    use kt_grids, only: gryfx
    integer :: is
    if (proc0) then 
      do is = 1, nspec
         spec(is)%stm = sqrt(spec(is)%temp/spec(is)%mass)
         spec(is)%zstm = spec(is)%z/sqrt(spec(is)%temp*spec(is)%mass)
         spec(is)%tz = spec(is)%temp/spec(is)%z
         spec(is)%zt = spec(is)%z/spec(is)%temp
         spec(is)%smz = abs(sqrt(spec(is)%temp*spec(is)%mass)/spec(is)%z)
         if (gryfx()) then
            !WARNING THIS NORMALISATION IS A GUESS!!!
            !FEEL FREE TO QUESTION IT AND WORRY ABOUT IT
           spec(is)%vnewk = 0.01 ! spec(is)%vnewk / sqrt(2.0)
         end if
      end do
    end if
    call broadcast (nspec)

    do is = 1, nspec
       call broadcast (spec(is)%dens)
       call broadcast (spec(is)%temp)
       call broadcast (spec(is)%fprim)
       call broadcast (spec(is)%tprim)
       call broadcast (spec(is)%vnewk)
       call broadcast (spec(is)%stm)
       call broadcast (spec(is)%zstm)
       call broadcast (spec(is)%tz)
       call broadcast (spec(is)%zt)
       call broadcast (spec(is)%smz)
    end do
  end subroutine calculate_and_broadcast_species_properties

  !> FIXME : Add documentation   
  subroutine write_trinity_parameters(trinpars_unit)
    integer, intent(in) :: trinpars_unit
    integer :: is
    do is = 1,nspec
      if (is<10) then
        write(trinpars_unit, "(A20,I1)") '&species_parameters_', is
      else
        write(trinpars_unit, "(A20,I2)") '&species_parameters_', is
      end if
      write (trinpars_unit, *) ' temp = ', spec(is)%temp
      write (trinpars_unit, *) ' dens = ', spec(is)%dens
      write (trinpars_unit, *) ' tprim = ', spec(is)%tprim
      write (trinpars_unit, *) ' fprim = ', spec(is)%fprim
      write (trinpars_unit, *) ' vnewk = ', spec(is)%vnewk
      write (trinpars_unit, "(A1)") '/'
    end do
  end subroutine write_trinity_parameters

  !> Set the module level config types
  !> Will abort if the module has already been initialised to avoid
  !> inconsistencies.
  subroutine set_species_config(species_config_in, species_element_config_in)
    use mp, only: mp_abort
    type(species_config_type), intent(in), optional :: species_config_in
    type(species_element_config_type), intent(in), dimension(:), allocatable, optional :: species_element_config_in

    if (initialized) then
       call mp_abort("Trying to set species_config when already initialized.", to_screen = .true.)
    end if
    if (present(species_config_in)) then
       species_config = species_config_in
    end if
    if (present(species_element_config_in)) then
       species_element_config = species_element_config_in
    end if    
  end subroutine set_species_config   

  !> Get the module level config instance
  function get_species_config()
    type(species_config_type) :: get_species_config
    get_species_config = species_config
  end function get_species_config

  !> Get the module level config instance
  function get_species_element_config()
    type(species_element_config_type), allocatable :: get_species_element_config(:)
    if (allocated(species_element_config)) then
      get_species_element_config = species_element_config
    else
      allocate(get_species_element_config(0))
    end if
  end function get_species_element_config

  !---------------------------------------
  ! Following is for the config_type
  !---------------------------------------

  !> Reads in the species_knobs namelist and populates the member variables
  subroutine read_species_config(self)
    use file_utils, only: input_unit_exist, get_indexed_namelist_unit
    use mp, only: proc0
    implicit none
    class(species_config_type), intent(in out) :: self
    logical :: exist
    integer :: in_file

    ! Note: When this routine is in the module where these variables live
    ! we shadow the module level variables here. This is intentional to provide
    ! isolation and ensure we can move this routine to another module easily.    
    integer :: nspec
    real :: me, zi_fac

    namelist /species_knobs/ me, nspec, zi_fac

    ! Only proc0 reads from file
    if (.not. proc0) return

    ! First set local variables from current values
    me = self%me
    nspec = self%nspec
    zi_fac = self%zi_fac

    ! Now read in the main namelist
    in_file = input_unit_exist(self%get_name(), exist)
    if (exist) read(in_file, nml = species_knobs)

    ! Now copy from local variables into type members
    self%me = me
    self%nspec = nspec
    self%zi_fac = zi_fac

    self%exist = exist
  end subroutine read_species_config

  !> Writes out a namelist representing the current state of the config object
  subroutine write_species_config(self, unit)
    implicit none
    class(species_config_type), intent(in) :: self
    integer, intent(in) , optional:: unit
    integer :: unit_internal

    unit_internal = 6 ! @todo -- get stdout from file_utils
    if (present(unit)) then
       unit_internal = unit
    endif

    call self%write_namelist_header(unit_internal)
    call self%write_key_val("me", self%me, unit_internal)
    call self%write_key_val("nspec", self%nspec, unit_internal)
    call self%write_key_val("zi_fac", self%zi_fac, unit_internal)
    call self%write_namelist_footer(unit_internal)
  end subroutine write_species_config

  !> Resets the config object to the initial empty state
  subroutine reset_species_config(self)
    class(species_config_type), intent(in out) :: self
    type(species_config_type) :: empty
    select type (self)
    type is (species_config_type)
       self = empty
    end select
  end subroutine reset_species_config

  !> Broadcasts all config parameters so object is populated identically on
  !! all processors
  subroutine broadcast_species_config(self)
    use mp, only: broadcast
    implicit none
    class(species_config_type), intent(in out) :: self
    call broadcast(self%me)
    call broadcast(self%nspec)
    call broadcast(self%zi_fac)

    call broadcast(self%exist)
  end subroutine broadcast_species_config

  !> Gets the default name for this namelist
  function get_default_name_species_config()
    implicit none
    character(len = CONFIG_MAX_NAME_LEN) :: get_default_name_species_config
    get_default_name_species_config = "species_knobs"
  end function get_default_name_species_config

  !> Gets the default requires index for this namelist
  function get_default_requires_index_species_config()
    implicit none
    logical :: get_default_requires_index_species_config
    get_default_requires_index_species_config = .false.
  end function get_default_requires_index_species_config
  
  !---------------------------------------
  ! Following is for the config_type
  !---------------------------------------
  
  !> Reads in the species_parameters namelist and populates the member variables
  subroutine read_species_element_config(self)
    use file_utils, only: input_unit_exist, get_indexed_namelist_unit
    use mp, only: proc0
    implicit none
    class(species_element_config_type), intent(in out) :: self
    logical :: exist
    integer :: in_file

    ! Note: When this routine is in the module where these variables live
    ! we shadow the module level variables here. This is intentional to provide
    ! isolation and ensure we can move this routine to another module easily.    
    character(len = 20) :: f0type, type
    real :: bess_fac, dens, dens0, fprim, mass, nu_h, temp, tpar0, tperp0, tprim, u0, uprim, uprim2, vcprim, vcrit, vnewk, z

    namelist /species_parameters/ bess_fac, dens, dens0, f0type, fprim, mass, nu_h, temp, tpar0, tperp0, tprim, type, u0, uprim, uprim2, vcprim, vcrit, vnewk, z

    ! Only proc0 reads from file
    if (.not. proc0) return

    ! First set local variables from current values
    bess_fac = self%bess_fac
    dens = self%dens
    dens0 = self%dens0
    f0type = self%f0type
    fprim = self%fprim
    mass = self%mass
    nu_h = self%nu_h
    temp = self%temp
    tpar0 = self%tpar0
    tperp0 = self%tperp0
    tprim = self%tprim
    type = self%type
    u0 = self%u0
    uprim = self%uprim
    uprim2 = self%uprim2
    vcprim = self%vcprim
    vcrit = self%vcrit
    vnewk = self%vnewk
    z = self%z

    ! Now read in the main namelist
    call get_indexed_namelist_unit(in_file, trim(self%get_name()), self%index, exist)
    if (exist) read(in_file, nml = species_parameters)
    close(unit = in_file)

    ! Now copy from local variables into type members
    self%bess_fac = bess_fac
    self%dens = dens
    self%dens0 = dens0
    self%f0type = f0type
    self%fprim = fprim
    self%mass = mass
    self%nu_h = nu_h
    self%temp = temp
    self%tpar0 = tpar0
    self%tperp0 = tperp0
    self%tprim = tprim
    self%type = type
    self%u0 = u0
    self%uprim = uprim
    self%uprim2 = uprim2
    self%vcprim = vcprim
    self%vcrit = vcrit
    self%vnewk = vnewk
    self%z = z

    self%exist = exist
  end subroutine read_species_element_config

  !> Writes out a namelist representing the current state of the config object
  subroutine write_species_element_config(self, unit)
    implicit none
    class(species_element_config_type), intent(in) :: self
    integer, intent(in) , optional:: unit
    integer :: unit_internal

    unit_internal = 6 ! @todo -- get stdout from file_utils
    if (present(unit)) then
       unit_internal = unit
    endif

    call self%write_namelist_header(unit_internal)
    call self%write_key_val("bess_fac", self%bess_fac, unit_internal)
    call self%write_key_val("dens", self%dens, unit_internal)
    call self%write_key_val("dens0", self%dens0, unit_internal)
    call self%write_key_val("f0type", self%f0type, unit_internal)
    call self%write_key_val("fprim", self%fprim, unit_internal)
    call self%write_key_val("mass", self%mass, unit_internal)
    call self%write_key_val("nu_h", self%nu_h, unit_internal)
    call self%write_key_val("temp", self%temp, unit_internal)
    call self%write_key_val("tpar0", self%tpar0, unit_internal)
    call self%write_key_val("tperp0", self%tperp0, unit_internal)
    call self%write_key_val("tprim", self%tprim, unit_internal)
    call self%write_key_val("type", self%type, unit_internal)
    call self%write_key_val("u0", self%u0, unit_internal)
    call self%write_key_val("uprim", self%uprim, unit_internal)
    call self%write_key_val("uprim2", self%uprim2, unit_internal)
    call self%write_key_val("vcprim", self%vcprim, unit_internal)
    call self%write_key_val("vcrit", self%vcrit, unit_internal)
    call self%write_key_val("vnewk", self%vnewk, unit_internal)
    call self%write_key_val("z", self%z, unit_internal)
    call self%write_namelist_footer(unit_internal)
  end subroutine write_species_element_config

  !> Resets the config object to the initial empty state
  subroutine reset_species_element_config(self)
    class(species_element_config_type), intent(in out) :: self
    type(species_element_config_type) :: empty
    select type (self)
    type is (species_element_config_type)
       self = empty
    end select
  end subroutine reset_species_element_config

  !> Broadcasts all config parameters so object is populated identically on
  !! all processors
  subroutine broadcast_species_element_config(self)
    use mp, only: broadcast
    implicit none
    class(species_element_config_type), intent(in out) :: self
    call broadcast(self%bess_fac)
    call broadcast(self%dens)
    call broadcast(self%dens0)
    call broadcast(self%f0type)
    call broadcast(self%fprim)
    call broadcast(self%mass)
    call broadcast(self%nu_h)
    call broadcast(self%temp)
    call broadcast(self%tpar0)
    call broadcast(self%tperp0)
    call broadcast(self%tprim)
    call broadcast(self%type)
    call broadcast(self%u0)
    call broadcast(self%uprim)
    call broadcast(self%uprim2)
    call broadcast(self%vcprim)
    call broadcast(self%vcrit)
    call broadcast(self%vnewk)
    call broadcast(self%z)

    call broadcast(self%exist)
  end subroutine broadcast_species_element_config 

  !> Gets the default name for this namelist
  function get_default_name_species_element_config()
    implicit none
    character(len = CONFIG_MAX_NAME_LEN) :: get_default_name_species_element_config
    get_default_name_species_element_config = "species_parameters"
  end function get_default_name_species_element_config

  !> Gets the default requires index for this namelist
  function get_default_requires_index_species_element_config()
    implicit none
    logical :: get_default_requires_index_species_element_config
    get_default_requires_index_species_element_config = .true.
  end function get_default_requires_index_species_element_config
end module species