SUBROUTINE DAF(sigma, delta, mval, dafx, kmax, hermit, ndaf, nxval, nd_daf)
USE Numeric_Kinds_Module
USE FileUnits_Module
!-----------------------------------------------------------------------
! 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
!-----------------------------------------------------------------------
IMPLICIT NONE
INTEGER lval, mval, kmax, mval2, kval, kmaxx, kdaf, ndaf, nxval, khalf, idaf, nd_daf
REAL(KIND=WP_Kind) sigma, delta, twopi, dafx(0:nd_daf,0:ndaf), hermit(0:mval+ndaf+1), harg, expfac, tol, fac
PARAMETER (tol=1.d-50)

WRITE(DAF_Unit,*)
WRITE(DAF_Unit,*)'   Begin(DAF)'
!----------------------------------------------------------------------- 
! calculate the daf deltas.
!-----------------------------------------------------------------------
twopi=8.d0*atan(1.d0)

mval2=mval/2
kmaxx=kmax
DO kval=0,kmax
   harg=(kval*delta)/(SQRT(2.d0)*sigma)
   expfac=EXP(-harg**2)
   CALL hermit_ex(harg, mval+ndaf, hermit, expfac)
   DO kdaf=0,ndaf
      khalf=kdaf/2
      dafx(kval,kdaf)=0.d0
      fac=1.d0
      DO idaf=1,khalf
         fac=fac*idaf
      ENDDO
      DO lval=0,mval2
         IF(lval.ne.0)fac=(lval+khalf)*fac/lval
         dafx(kval,kdaf)=dafx(kval,kdaf)+hermit(2*lval+kdaf)*fac
      ENDDO
      dafx(kval,kdaf)=(-1.d0/(SQRT(2.d0)*sigma))**kdaf*dafx(kval,kdaf)*delta/(SQRT(twopi)*sigma)*(-4.d0)**khalf
   ENDDO

   IF(ABS(dafx(kval,ndaf)).gt.tol)THEN
      kmaxx=kval
   ELSE IF(ABS(dafx(kval,ndaf)).lt.tol) THEN
      GO TO 5
   ENDIF

ENDDO

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

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

WRITE(DAF_Unit,*)'   kmaxx=',kmaxx
kmax=kmaxx

WRITE(DAF_Unit,*)'   END(DAF)'
WRITE(DAF_Unit,*)

RETURN
END SUBROUTINE DAF