SUBROUTINE vibrational(s_dim, s_vals, s_min, s_max, infunction, vibdim, vib, io)

USE Numeric_Kinds_Module
USE CommonInfo_Module
USE Jacobi_Module
USE APH_Module
USE SurfaceJacobi_Module
IMPLICIT NONE
!=========================================================================================
!    Written by: Jeff Crawford
!
! This routine maps the vibrational wave functions to the APH grid via interpolation.
! Routine is called for initial wavepacket and for final state eigenfunctions. 
!
! Variables:
!   infunction = array containing the vibrational eigenfunctions before interpolation
!   io         = set to 'I' if infunction is the initial vbrational state, and set to 'O' 
!                if its an output product state
!   s_dim      = dimension of Jacobi small s grid
!   s_max      = max Jacobi small s grid value
!   s_min      = min Jacobi small s grid value
!   s_vals     = array containing the Jacobi small s grid points
!   vib        = array containing mapped vibrational function
!   vibdim     = dimension of APH grid onto which infunction will be mapped
!=========================================================================================
!     I N P U T
CHARACTER(LEN=1)    :: io
INTEGER,INTENT(IN)  :: vibdim, s_dim
REAL(dp),INTENT(IN) :: s_min, s_max, infunction(s_dim), s_vals(s_dim)

!=========================================================================================
!     O U T P U T
REAL(dp),INTENT(OUT) :: vib(vibdim)

!=========================================================================================
!     I N T E R N A L S
INTEGER  :: irho, itheta, ichi, ispt, ipt, rhostart, rhoend
REAL(dp) :: vcoord, vibmag, interp_vib

!=========================================================================================
!     A L L O C A T A B L E
REAL(dp),ALLOCATABLE :: y2_vib(:)

!=========================================================================================
!     F U N C T I O N S
INTEGER :: aphindex

  IF (io.eq.'I') THEN
    rhostart=1
    rhoend=nrho
  ELSEIF (io.eq.'O') THEN
    rhostart=ind_rho_infty
    rhoend=rhostart
  ELSE
    PRINT*,'Incorrect I/O designation in vibrational'
    STOP'vibrational'
  ENDIF


  ALLOCATE(y2_vib(s_dim))

  y2_vib=0.d0
  ipt=0

  CALL spline(s_vals,infunction,s_dim,0.d0,0.d0,y2_vib)

  DO irho=rhostart,rhoend 
    DO itheta=1,ntheta 
      DO ichi=1,nchi 

        ispt=aphindex(ntheta,nchi,irho,itheta,ichi,.false.)
        ipt=ipt+1
 
        vcoord=smalls(ispt)

!        IF (ichi.le.nchieven.and.vcoord.le.smax.and.vcoord.ge.smin) THEN 
        IF (vcoord.le.s_max.and.vcoord.ge.s_min) THEN 
          CALL splint( s_dim, s_vals, infunction, y2_vib, vcoord, interp_vib )
          vib(ipt)=interp_vib
          vibmag=vib(ipt)
          !vib(ipt)=CMPLX(interp_vib,0.d0,dp)
          !vibmag=(REAL(vib(ipt)))**2+(AIMAG(vib(ipt))**2)
        ELSE 
          vib(ipt)=0.d0
          !vib(ipt)=Cmplx(0.d0,0.d0,dp)
          vibmag=0.d0
        ENDIF

      ENDDO 
    ENDDO 
  ENDDO 

DEALLOCATE(y2_vib)

 Return
 End Subroutine vibrational