SUBROUTINE get_asymp_rho

USE Numeric_Kinds_Module
USE CommonInfo_Module
USE APH_Module
USE NIP_Module

IMPLICIT NONE

!=========================================================================================
! Written by: Jeff Crawford
!
! This routine obtains user-defined rho values that define where final state
! analysis occurs and where the negative imaginary potential inner boundary
! is located. These values are matched to the closest rho grid point.
!
! Variables:
!   nip_rho        = actual rho value for inner boundary of negative imaginary
!                    potential (matched to grid)
!   nip_rho_user   = user-defined rho value for inner boundary of the negative 
!                    imaginary potential (not matched to grid)
!   order          = order of negative imaginary potential function
!   rho_infity     = actual rho value for final analysis (matched to grid)
!   rho_infty_user = user-defined rho value for final wave packet analysis (not
!                    matched to grid)
!   ufac           = amplitude of negative imaginary potential (in eV)


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

LOGICAL  :: foundrhoinf
INTEGER  :: istat, irho
REAL(dp) :: rho_infty_user, nip_rho_user

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

NAMELIST / asymp_rho / rho_infty_user
NAMELIST / nip_input / nip_rho_user, order, ufac

!=========================================================================================
! Read in NAMELIST parameters

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

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

!=========================================================================================
! Test that rho_infty_user is within existing grid
 
 IF (rho_infty_user.ge.rho_val(nrho)) THEN
   PRINT*,'RHO_INFTY_USER is not on existing grid'
   PRINT*,'Maximum rho value = ', rho_val(nrho)
   PRINT*,'rho_infty_user = ', rho_infty_user
   STOP'GET_ASYMP_RHO'
 ENDIF

!=========================================================================================
! Test that nip_rho_user is within existing grid

 IF (nip_rho_user.ge.rho_val(nrho)) THEN
   PRINT*,'NIP_RHO_USER is not on existing grid'
   PRINT*,'Maximum rho value = ', rho_val(nrho)
   PRINT*,'nip_rho_user = ', nip_rho_user
   STOP'GET_ASYMP_RHO'
 ENDIF

!=========================================================================================
! Check that rho_infty is not inside negative imaginary potential region.

 IF (rho_infty_user.ge.nip_rho_user) THEN
   PRINT*,'RHO_INFTY_USER is inside the negative imaginary potential'
   PRINT*,'Minimum NIP boundary = ', nip_rho_user
   PRINT*,'rho_infty_user = ', rho_infty_user
   STOP'GET_ASYMP_RHO'
 ENDIF

!=========================================================================================
!  Pick first value of rho that is greater than or equal to rho_infty_user

 foundrhoinf=.false.

 DO irho=1,nrho
   IF (.not.foundrhoinf) THEN
     IF (rho_val(irho).ge.rho_infty_user) THEN
       rho_infty=rho_val(irho)
       ind_rho_infty=irho
       foundrhoinf=.true.
     ENDIF
   ENDIF
 ENDDO

!=========================================================================================
!  Pick first value of rho that is greater than or equal to nip_rho_user

 foundrhoinf=.false.

 DO irho=1,nrho
   IF (.not.foundrhoinf) THEN
     IF (rho_val(irho).ge.nip_rho_user) THEN
       nip_rho=rho_val(irho)
       ind_nip_rho=irho
       foundrhoinf=.true.
     ENDIF
   ENDIF
 ENDDO

!=========================================================================================
! Write info to log file and screen

 WRITE(Out_Unit,21) 'Rho Infty and NIP Info:','au'
 WRITE(Out_Unit,22) 'rho_infty_user', rho_infty_user
 WRITE(Out_Unit,22) 'rho_infty     ', rho_infty
 WRITE(Out_Unit,22) 'nip_rho_user  ', nip_rho_user
 WRITE(Out_Unit,22) 'nip_rho       ', nip_rho

 WRITE(*,23) 'rho_infty', rho_infty
 WRITE(*,23) 'nip_rho  ', nip_rho

 21 FORMAT(/1X,a,T29,a,/1X,50('-'))
 22 FORMAT(1X,a,10('.'),f7.3)
 23 FORMAT(1X,a,6('.'),f7.3)

END SUBROUTINE get_asymp_rho