program species_gen

! GWH: generate various species-related quantities, including the 
! collisionality of each species.
!
! to compile and run, do:
!
! f90 -o species_gen species_gen.f90
! species_gen < species_gen.in
!
! could modify this to read a gs2 namelist file, to check the species 
! parameter lists for consistency and to calculate the proper collisionality.

  implicit none
  integer nspec
  ! allocatable arrays can't be in a namelist, so hardwire max size:
  integer, parameter :: maxspec=10
  real, dimension(maxspec) :: Temp, dens, tprim, fprim, mass, Z
  real :: a_ref, R_major, n_e, pi, T_e, T_i, T_ref, beta_ref
  real :: coulog, vnewk_e, vnewk_i, vnewk_Z, sum, gnu, mass_ref
  integer i, ib, ie

  namelist /species_knobs/ nspec
  namelist /species_parameters/ mass, Z, dens, Temp, fprim, tprim
  namelist /plasma_parameters/ a_ref, R_major, T_ref, mass_ref, n_e, beta_ref

  pi=abs(acos(-1.0))

  read(5,nml=species_knobs)
  read(5,nml=species_parameters)
  read(5,nml=plasma_parameters)

  ! find the electron index
  do i=1, nspec
     if(Z(i) .lt. 0.0) then
        ie=i
        go to 30
     endif
  enddo
  write(*,*) "No electrons detected."
  stop
30 T_e=T_ref*temp(ie)

  T_i=T_ref*temp(1)    ! Assumes that the first species is ions


  coulog=15.94-0.5*alog(n_e/1.e19/T_e**2)
  write(*,*) "coulog =", coulog
  write(*,*) "n_e, T_e, a_ref, mass_ref, T_ref =", &
       n_e, T_e, a_ref, mass_ref, T_ref
  vnewk_e=0.00279*coulog*n_e/1.e19/T_e**1.5*a_ref*(mass_ref/T_ref)**0.5
  write(*,*) "vnewk_e =", vnewk_e
  
  gnu=2.5e-7*n_e/1.e6/T_e**1.5/T_i**0.5/1.e6*R_major
  write(*,*) "IFS-PPPL subroutine gnu=", gnu

  sum=0.0
  do i=1,nspec
     if(Z(i) .gt. 0.0) sum=sum+dens(i)*Z(i)**2
  enddo
  write(*,*) "Z_eff = ", sum

  sum=0.0
  do i=1,nspec
     sum=sum+dens(i)*Z(i)
  enddo
  write(*,*) "Sum_i dens(i)*Z(i) = ", sum, &
       " (=0 for quasineutrality.)"

  sum=0.0
  do i=1,nspec
     sum=sum+dens(i)*Z(i)*fprim(i)
  enddo
  write(*,*) "Sum_i dens(i)*Z(i)*fprim(i) = ", sum, &
       " (=0 for radial quasineutrality.)"

  sum=0.0
  do i=1,nspec
     sum=sum+beta_ref*dens(i)*temp(i)*(fprim(i)+tprim(i))
  enddo
  write(*,*) "beta_prime_input  = ", sum

  do i=1, nspec
     sum=0.0
     do ib=1,nspec
        sum=sum + &
             dens(ib)*Z(ib)**2*2/(1+(mass(i)/mass(ib)*Temp(ib)/Temp(i))**0.5)
     enddo
     vnewk_i = vnewk_e * &
          sqrt(1.0/1836/mass(i))*(Temp(ie)/Temp(i))**1.5 *Z(i)**2 *sum 
     if(Z(i) .gt. 0) then
        write(*,*) "i, vnewk_i      =", i, vnewk_i
     else
        write(*,*) "electron vnewk_e=  ", vnewk_e
     endif
  enddo

end