module species
  implicit none

  private
  public :: init_species, finish_species !, reinit_species, init_trin_species, finish_trin_species

  public :: write_trinity_parameters
  public :: wnml_species, check_species
  public :: nspec, specie, spec
  public :: ion_species, electron_species, tracer_species
  public :: f0_maxwellian, f0_tabulated, f0_sdanalytic
  public :: has_ion_species, has_electron_species, has_nonmaxw_species
  public :: calculate_slowingdown_parameters
  public :: ions, electrons, impurity
  public :: set_overrides
  public :: determine_species_order
  public :: calculate_f0_arrays, set_current_f0_species, eval_f0
  public :: nonmaxw_corr, f0_values, dlnf0drho

  type :: specie
     real :: z
     real :: mass
     real :: dens, dens0, u0
     real :: tpar0,tperp0
     real :: temp
     real :: tprim
     real :: fprim
     real :: uprim, uprim2
     real :: vnewk, nustar
     real :: nu, nu_h  ! nu will be the preferred collisionality parameter moving forward
                       ! as it will have a consistent v_t scaling
                       ! 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 :: 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

  integer :: ions, electrons, impurity
  integer :: ntspec_trin
  real, dimension (:), allocatable :: dens_trin, temp_trin, fprim_trin, tprim_trin, nu_trin
  real, dimension (:,:), allocatable :: nonmaxw_corr, f0_values, dlnf0drho

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

contains
  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 == 1) write (report_unit, fmt="('    Type:             Ion')")
       if (spec(is)%type == 2) write (report_unit, fmt="('    Type:             Electron')")
       if (spec(is)%type == 3) write (report_unit, fmt="('    Type:             Slowing-down')")
       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
    
    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

  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 == 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
  
  subroutine init_species
    !use mp, only: trin_flag
    implicit none

    if (initialized) return
    initialized = .true.

    call read_parameters
    !if (trin_flag) then
       !call reinit_species (ntspec_trin, dens_trin, &
         !temp_trin, fprim_trin, tprim_trin, nu_trin)
    !endif
  end subroutine init_species

  subroutine read_parameters
    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, broadcast, mp_abort
    implicit none
    real :: z, mass, dens, dens0, u0, temp, tprim, fprim, uprim, uprim2, vnewk, nustar, nu, nu_h
    real :: tperp0, tpar0, bess_fac, vcrit, vcprim
    real, dimension(:), allocatable :: tmp_bcast
    character(20) :: type, f0type
    integer :: unit
    integer :: is, iostat
    namelist /species_knobs/ nspec, ZI_fac, me
    namelist /species_parameters/ z, mass, dens, dens0, u0, temp, &
         tprim, fprim, uprim, uprim2, vnewk, nustar, type, nu, nu_h, &
         tperp0, tpar0, bess_fac, f0type, vcrit, vcprim
    integer :: ierr, in_file

    type (text_option), dimension (5), 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) /)

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

    if (proc0) then
       nspec = 2
       ZI_fac = -1.0
       me = -1.0
       in_file = input_unit_exist("species_knobs", exist)
!       if (exist) read (unit=input_unit("species_knobs"), nml=species_knobs)
       if (exist) then
          read (unit=in_file, nml=species_knobs)
       else 
          write(6,*) 'Error: species_knobs namelist file does not exist'
       endif
       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

    call broadcast (nspec)
    call broadcast (ZI_fac)
    call broadcast (me)
    allocate (spec(nspec))

    if (proc0) then
       do is = 1, nspec
          call get_indexed_namelist_unit (unit, "species_parameters", is)
          z = 1
          mass = 1.0
          dens = 1.0
          dens0 = 1.0
          u0 = 1.0
          tperp0 = 0.
          tpar0 = 0.
          temp = 1.0
          tprim = 6.9
          fprim = 2.2
          uprim = 0.0
          uprim2 = 0.0
          nustar = -1.0
          vnewk = 0.0
          nu = -1.0
          nu_h = 0.0
          bess_fac=1.0
          vcrit = -1.0
          vcprim = -999.9
          type = "default"
          f0type = "maxwellian"
          read (unit=unit, nml=species_parameters, iostat=iostat)
          if(iostat /= 0) write(6,*) 'Error ',iostat,'reading species parameters'
          close (unit=unit, iostat=iostat)
          if(iostat /= 0) write(6,*) 'Error ',iostat,'closing species parameters namelist'
          !write (*,*) 'type is ', type
          spec(is)%z = z
          spec(is)%mass = mass
          spec(is)%dens = dens
          spec(is)%dens0 = dens0
          spec(is)%u0 = u0
          spec(is)%tperp0 = tperp0
          spec(is)%tpar0 = tpar0
          spec(is)%temp = temp
          spec(is)%tprim = tprim
          spec(is)%fprim = fprim
          spec(is)%uprim = uprim
          spec(is)%uprim2 = uprim2
          spec(is)%vnewk = vnewk
          spec(is)%nustar = nustar
          spec(is)%nu = nu
          spec(is)%nu_h = nu_h

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

          spec(is)%bess_fac=bess_fac
          spec(is)%vcrit=vcrit
          spec(is)%vcprim=vcprim

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

    end if

    !We use a temporary array here to allow us to broadcast all species data at once
    !May be possible to avoid this by simply calling broadcast(spec%z) etc.
    allocate(tmp_bcast(nspec))
    tmp_bcast=spec%z ; call broadcast(tmp_bcast) ; spec%z=tmp_bcast
    tmp_bcast=spec%mass ; call broadcast(tmp_bcast) ; spec%mass=tmp_bcast
    tmp_bcast=spec%dens ; call broadcast(tmp_bcast) ; spec%dens=tmp_bcast
    tmp_bcast=spec%dens0 ; call broadcast(tmp_bcast) ; spec%dens0=tmp_bcast
    tmp_bcast=spec%u0 ; call broadcast(tmp_bcast) ; spec%u0=tmp_bcast
    tmp_bcast=spec%tperp0 ; call broadcast(tmp_bcast) ; spec%tperp0=tmp_bcast
    tmp_bcast=spec%tpar0 ; call broadcast(tmp_bcast) ; spec%tpar0=tmp_bcast
    tmp_bcast=spec%temp ; call broadcast(tmp_bcast) ; spec%temp=tmp_bcast
    tmp_bcast=spec%tprim ; call broadcast(tmp_bcast) ; spec%tprim=tmp_bcast
    tmp_bcast=spec%fprim ; call broadcast(tmp_bcast) ; spec%fprim=tmp_bcast
    tmp_bcast=spec%uprim ; call broadcast(tmp_bcast) ; spec%uprim=tmp_bcast
    tmp_bcast=spec%uprim2 ; call broadcast(tmp_bcast) ; spec%uprim2=tmp_bcast
    tmp_bcast=spec%vnewk ; call broadcast(tmp_bcast) ; spec%vnewk=tmp_bcast
    tmp_bcast=spec%nu ; call broadcast(tmp_bcast) ; spec%nu=tmp_bcast
    tmp_bcast=spec%nu_h ; call broadcast(tmp_bcast) ; spec%nu_h=tmp_bcast
    tmp_bcast=spec%nustar ; call broadcast(tmp_bcast) ; spec%nustar=tmp_bcast
    tmp_bcast=spec%stm ; call broadcast(tmp_bcast) ; spec%stm=tmp_bcast
    tmp_bcast=spec%zstm ; call broadcast(tmp_bcast) ; spec%zstm=tmp_bcast
    tmp_bcast=spec%tz ; call broadcast(tmp_bcast) ; spec%tz=tmp_bcast
    tmp_bcast=spec%zt ; call broadcast(tmp_bcast) ; spec%zt=tmp_bcast
    tmp_bcast=spec%smz ; call broadcast(tmp_bcast) ; spec%smz=tmp_bcast
    tmp_bcast=spec%type ; call broadcast(tmp_bcast) ; spec%type=int(tmp_bcast)
    tmp_bcast=spec%f0type ; call broadcast(tmp_bcast) ; spec%f0type=int(tmp_bcast)
    tmp_bcast=spec%bess_fac ; call broadcast(tmp_bcast) ; spec%bess_fac=tmp_bcast
    tmp_bcast=spec%vcrit ; call broadcast(tmp_bcast) ; spec%vcrit=tmp_bcast
    tmp_bcast=spec%vcprim ; call broadcast(tmp_bcast) ; spec%vcprim=tmp_bcast
    deallocate(tmp_bcast)

  end subroutine read_parameters

  pure function has_ion_species (spec)
    implicit none
    type (specie), dimension (:), intent (in) :: spec
    logical :: has_ion_species
    has_ion_species = any(spec%type == ion_species)
  end function has_ion_species

  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 (spec(is)%type .eq. electron_species)                            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

  subroutine calculate_f0_arrays(egrid)
     use mp, only: iproc
     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
           !write(*,*) iproc, is, ie, f0_values(ie,is), -nonmaxw_corr(ie,is)/spec(is)%temp,dlnf0drho(ie,is)
        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

  subroutine set_current_f0_species(is)
     implicit none
     integer, intent(in):: is
     is_global = is
  end subroutine set_current_f0_species


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

  subroutine eval_f0_sdanalytic(is,v,f0,df0dE_out,f0prim)
    use constants, only: pi
    use mp, only: mp_abort, proc0
    use spfunc, only: erf => erf_ext
    implicit none
    integer,intent(in):: is
    real,intent(in):: v
    real,intent(inout):: f0,df0dE_out, f0prim
    real:: A, df0dv,vcprim, vcrit
    integer:: ie, i, electron_spec

    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, il, it, negrid
    logical:: acceptable_spline
    real:: test, moment0, moment2, 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

  pure function has_electron_species (spec)
    implicit none
    type (specie), dimension (:), intent (in) :: spec
    logical :: has_electron_species
    has_electron_species = any(spec%type == electron_species)
  end function has_electron_species

  pure function has_nonmaxw_species (spec)
    implicit none
    type (specie), dimension (:), intent (in) :: spec
    logical :: has_nonmaxw_species
    has_nonmaxw_species = .NOT. all(spec%f0type == f0_maxwellian)
  end function has_nonmaxw_species

  subroutine finish_species
    implicit none
    if(allocated(spec)) deallocate (spec)
!    if(allocated(f0_values)) deallocate (f0_values,nonmaxw_corr,dlnf0drho)
    initialized = .false.
  end subroutine finish_species

  !> This routine is used only for Trinity
  subroutine determine_species_order
    use mp, only: broadcast, proc0, mp_abort

    implicit none

    integer is
     if (nspec == 1) then
        ions = 1
        electrons = 0
        impurity = 0
     else
        ! if 2 or more species in GS2 calculation, figure out which is main ion
        ! and which is electron via mass (main ion mass assumed to be one)
        do is = 1, nspec
           if (abs(spec(is)%mass-1.0) <= epsilon(0.0)) then
              ions = is
           else if (spec(is)%mass < 0.3) then
              ! for electrons, assuming electrons are at least a factor of 3 less massive
              ! than main ion and other ions are no less than 30% the mass of the main ion
              electrons = is
           else if (spec(is)%mass > 1.0 + epsilon(0.0)) then
              impurity = is
           else
              if (proc0) write (*,*) &
                   "Error: TRINITY requires the main ions to have mass 1", &
                   "and the secondary ions to be impurities (mass > 1)"
              call mp_abort("Error: TRINITY requires the main ions to have mass 1 and the secondary ions to be impurities (mass > 1)")
           end if
        end do
     end if
  end subroutine determine_species_order

  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


  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

  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

end module species
