SUBROUTINE get_energy

USE Numeric_Kinds_Module
USE CommonInfo_Module
USE Convrsns_Module
USE EnergyGrid_Module
Use Cheby_Module
USE Gaussian_Module
USE masses_module
IMPLICIT NONE

!=========================================================================================
! Written by: Jeff Crawford
!
! This routine reads in the energy parameters and calculates the energy grid.
! The S-matrix will be calulated at each energy grid point. Thi routine also
! reads in the gaussian wave packet parameters.
!
! Variables:
!   deltae = energy grid spacing
!   edim   = Number of energy grid points
!   emax   = Maximum value of the energy grid
!   emin   = Minimum value of the energy grid
!   e_vals = array of energygrid values
!   k0     = wave number center of gaussian wave packet
!   sigma  = width parameter for the gaussian wave packet
!   s0     = coordinate center of gaussian wave packet (mass-scaled Jacobi coordinate cap S)
!   vcut   = Maximum allowed eigenvalue for the potential energy and the individual
!            rho, theta and chi kinetic energy matrices.

!=========================================================================================
!     I N T E R N A L S

INTEGER  :: ie, istat
REAL(dp) :: emax, emin,e0
!=========================================================================================
!   N A M E L I S T

NAMELIST / ek_dim / edim, emin, emax
NAMELIST / pot_cut / vcut
NAMELIST / aphgauss / sigma, k0, s0

!=========================================================================================
! Read in namelists

 OPEN(UNIT=Input_Unit,FILE=Trim(inputfile1),IOSTAT=istat,STATUS='old',ACTION='read')
 IF (istat.ne.0) STOP 'OPEN failed - get_energy'
 READ(Input_Unit,ek_dim)
 CLOSE(UNIT=Input_Unit,IOSTAT=istat,STATUS='keep') 
 IF (istat.ne.0) STOP 'CLOSE failed - get_energy'

 OPEN(UNIT=Input_Unit,FILE=Trim(inputfile1),IOSTAT=istat,STATUS='old',ACTION='read')
 IF (istat.ne.0) STOP 'OPEN failed - get_energy'
 READ(Input_Unit,pot_cut)
 CLOSE(UNIT=Input_Unit,IOSTAT=istat,STATUS='keep')
 IF (istat.ne.0) STOP 'CLOSE failed - get_energy'

 OPEN(UNIT=Input_Unit,FILE=Trim(inputfile1),IOSTAT=istat,STATUS='old',ACTION='read')
 IF (istat.ne.0) STOP 'OPEN failed - get_energy'
 READ(Input_Unit,aphgauss)
 CLOSE(UNIT=Input_Unit,IOSTAT=istat,STATUS='keep')
 IF (istat.ne.0) STOP 'CLOSE failed - get_energy'

!=========================================================================================
! Check for valid input

IF (emin.gt.emax) THEN
  PRINT*,'ERROR!! emin > emax'
  STOP 'get_energy_grid'
ENDIF

IF (edim.le.0) THEN
  PRINT*,'ERROR!! edim < 1'
  STOP 'get_energy_grid'
ENDIF

IF (s0.le.0.d0.or.sigma.le.0.or.k0.eq.0.d0) THEN
  PRINT*,'ERROR!! Gaussian parameters out-of-bounds'
  STOP 'get_energy_grid'
ENDIF

IF (vcut.le.emax) THEN
  PRINT*,'ERROR!! vcut <= emax'
  STOP 'get_energy_grid'
ENDIF

!=========================================================================================
! Allocate energy grid

 ALLOCATE( e_vals(edim) )

!========================================================
!  Obtain energy grid

 deltae=(emax-emin)/(edim-1)

 DO ie=1,edim
   e_vals(ie)=emin+(ie-1)*deltae
 ENDDO

!========================================================
!  Write energy info

 e0=k0**2/usys2

 WRITE(Out_Unit,21) 'Energy Grid Info:','au','eV'
 WRITE(Out_Unit,22) 'emin  ', emin, emin*autoev
 WRITE(Out_Unit,22) 'emax  ', emax, emax*autoev
 WRITE(Out_Unit,22) 'deltae', deltae, deltae*autoev
 WRITE(Out_Unit,23) 'edim',edim

 WRITE(*,24) 'edim',edim
 WRITE(*,25) 'emin',emin
 WRITE(*,25) 'emax',emax
 WRITE(*,25) 'vcut',vcut

 WRITE(Out_Unit,26) 'Gaussian Info:','au','eV'
 WRITE(Out_Unit,27) 'sigma',sigma 
 WRITE(Out_Unit,27) 's0   ',s0
 WRITE(Out_Unit,27) 'k0   ',k0
 WRITE(Out_Unit,28) 'E0   ',e0,e0*autoev

 WRITE(*,27) 'sigma',sigma
 WRITE(*,27) 's0   ',s0
 WRITE(*,27) 'k0   ',k0

 21 FORMAT(/1X,a,T25,a,T34,a,/1X,50('-'))
 22 FORMAT(1X,a,15('.'),f7.4,3X,f6.3)
 23 FORMAT(1X,a,4('.'),I4)
 24 FORMAT(1X,a,11('.'),I6)
 25 FORMAT(1X,a,11('.'),f7.4,' au')
 26 FORMAT(/1X,a,T19,a,7X,a,/1X,50('-'))
 27 FORMAT(1X,a,10('.'),f6.3)
 28 FORMAT(1X,a,10('.'),f6.3,3X,f6.3)

END SUBROUTINE get_energy