SUBROUTINE daf (sigma, delta, mval, dafx, kmax, hermit, ndaf, nxval, nd_daf)
USE Numeric_Kinds_Module
USE FileUnits_Module
IMPLICIT NONE
!==============================================================================
! This routine was written by G. A. Parker
! IF you find an error or have an improvement please send a messge to
! parker@phyast.nhn.uoknor.edu
!-----------------------------------------------------------------------
INTEGER :: lval, mval, kmax, mval2, kval, kmaxx, kdaf, ndaf, nxval, khalf, idaf, nd_daf
INTEGER :: wp
REAL (dp) :: sigma, delta, twopi, harg, expfac, fac
REAL (dp) :: dafx (0:nd_daf, 0:ndaf), hermit ( 0:mval + ndaf + 1) 
REAL (dp),PARAMETER :: tol = 1.d-50
!==============================================================================
! calculate the daf deltas.

twopi = 8.d0 * atan (1.d0)
mval2 = mval / 2
kmaxx = kmax
DO 2 kval = 0, kmax
   harg = (kval * delta) / (sqrt (2.d0) * sigma)
   expfac = exp ( - harg**2)
   CALL hermit_ex (harg, mval + ndaf, hermit, expfac)
   DO 33 kdaf = 0, ndaf
      khalf = kdaf / 2
      dafx (kval, kdaf) = 0.d0
      fac = 1.d0
      DO 34 idaf = 1, khalf
         fac = fac * idaf
   34       END DO
      DO 1 lval = 0, mval2
         IF (lval.ne.0) fac = (lval + khalf) * fac / lval
         dafx (kval, kdaf) = dafx (kval, kdaf) + hermit (2 * lval + kdaf) * fac
    1       END DO
      dafx (kval, kdaf) = ( - 1.d0 / (sqrt (2.d0) * sigma) ) ** &
       kdaf * dafx (kval, kdaf) * delta / (sqrt (twopi) * sigma)  * ( - 4.d0) **khalf

   33    END DO
   IF (abs (dafx (kval, ndaf) ) .gt.tol) THEN
      kmaxx = kval
   ELSEIF (abs (dafx (kval, ndaf) ) .lt.tol) THEN
      goto 5
   ENDIF

    2 END DO

WRITE(Output_Unit, * ) 'Program must be modified to increase dimension'
WRITE(Output_Unit, * ) 'on array dafx'
WRITE(Output_Unit, * ) 'kmaxx=', kmaxx
WRITE(Output_Unit, * ) 'mval=', mval
WRITE(Output_Unit,  * ) '1  nd_daf=', nd_daf

STOP 'daf'
    5 IF (kmaxx.gt.nd_daf) THEN
   WRITE(Output_Unit, * ) 'Program must be modified to increase dimension'
   WRITE(Output_Unit, * ) 'on array dafx'
   WRITE(Output_Unit, * ) 'kmaxx=', kmaxx
   WRITE(Output_Unit, * ) 'mval=', mval
WRITE(Output_Unit,  * ) '2  nd_daf=', nd_daf
   STOP 'daf'

ENDIF
WRITE(Output_Unit, * ) 'kmaxx=', kmaxx

kmax = kmaxx
WRITE(Output_Unit, * ) 'END(DAF)'

WRITE(Output_Unit, * )
RETURN
END SUBROUTINE daf