SUBROUTINE aph_kinetic !(t_rho,t_theta,t_chi,a,b,c,d,rhom2)
USE Numeric_Kinds_Module
USE FileUnits_Module
USE Convrsns_Module
USE APH_Module
USE Aph_Ham_Module
USE masses_module
!-----------------------------------------------------------------------------------------
! This routine was written by G. A. Parker and modified by J. Crawford
! IF you find an error or have an improvement please send a messge to
! parker@phyast.nhn.uoknor.edu
!-----------------------------------------------------------------------------------------
IMPLICIT NONE
LOGICAL :: debug
INTEGER :: kmax, itheta, irho
INTEGER,PARAMETER :: mval=32, nd_daf=40, ndaf=2 
REAL (dp),ALLOCATABLE :: dafx(:,:), hermit(:),toe(:)
!==============================================================================
! Allocate necessary arrays

ALLOCATE (t_rho (nrho , nrho), t_theta (ntheta , ntheta) )
ALLOCATE (t_chi (nchi , nchi))!, t_chir (nchi *nchi) )
ALLOCATE (a (ntheta), b (ntheta), c (ntheta), d (ntheta) )
!ALLOCATE (t_chi1 (nchi * nchi), t_chir1 (nchi * nchi) )
ALLOCATE (rhom2 (nrho) )
ALLOCATE(dafx(0:nd_daf, 0:ndaf), hermit(0:mval + ndaf + 1))
ALLOCATE(toe(0:nd_daf))

!-----------------------------------------------------------------------------------------
! Calculate 1/rho**2
!-----------------------------------------------------------------------------------------
DO irho = 1, nrho
   rhom2 (irho) = 1.d0 / rho_val (irho) **2
ENDDO
!-----------------------------------------------------------------------------------------
! Calculate rho DAF
!-----------------------------------------------------------------------------------------
 IF (nrho.gt.1) THEN
   CALL get_daf (nrho, rho_val, deltarho, mval, dafx, hermit, ndaf,  nd_daf, toe, kmax)
 ENDIF
!------------------------------------------------------------------------------
! Calculate rho Kinetic energy term.
!------------------------------------------------------------------------------
 CALL get_kinetic (nrho,toe,t_rho,rho_val,kmax,'rho     ',dafx,nd_daf,ndaf)

!-----------------------------------------------------------------------------------------
! Calculate theta DAF
!-----------------------------------------------------------------------------------------
 CALL DVR_kinetic(ntheta,t_theta,theta_val)
!-----------------------------------------------------------------------------------------
! Calculate A, B, C, and D
!----------------------------------------------------------------------------------------- 
 DO itheta = 1, ntheta
   a (itheta) = 1.d0 / (1.d0 + sin (theta_val (itheta) ) )
   c (itheta) = 1.d0 / (1.d0 - sin (theta_val (itheta) ) )
   d (itheta) = cos (theta_val (itheta) ) / sin (theta_val (itheta) ) **2
 ENDDO
 DO itheta = 1, ntheta
   b (itheta) = 0.5d0 / sin (theta_val (itheta)) **2
   IF(abs(theta_val(itheta)).lt.1.0d-7)b (itheta)=10000.0d4
 ENDDO
!-----------------------------------------------------------------------------------------
! Calculate chi DAF.
!-----------------------------------------------------------------------------------------
 IF (nchi.gt.1) THEN
   CALL get_daf (nchi, chi_val, deltachi, mval, dafx, hermit, ndaf, nd_daf, toe, kmax)
 ENDIF
!-----------------------------------------------------------------------------------------
! Calculate chi Kinetic energy term.
!-----------------------------------------------------------------------------------------
 CALL get_kinetic (nchi, toe, t_chi, chi_val, kmax, 'chi     ', dafx, nd_daf, ndaf)

 DEALLOCATE(dafx,hermit,toe)

RETURN
END SUBROUTINE aph_kinetic