SUBROUTINE Simpson (npt, xpt, wht, a, b)
USE Numeric_Kinds_Module
USE Precision_Module
USE FileUnits_Module
IMPLICIT NONE
INTEGER(KIND=IW_Kind) :: i, npt
REAL(KIND=WP_Kind)    :: xpt(npt), wht(npt), a, b, h
!---------------------------------------------------------------------
!
!  purpose
!     supplies the Simpson abscissas and weights
!
!  description of parameters
!     npt-number of points desired
!     xpt-resultant array containing the gauss legenedre abscissas.
!     wht-resultant array containing the gauss legenedre weights
!---------------------------------------------------------------------

IF(2*(Npt/2).eq.Npt)THEN
  WRITE(Out_Unit,*)'Error in Simpson Number of Points=',Npt
  WRITE(Out_Unit,*)'Number of points must be odd'
  STOP 'Error in Simpson'
ENDIF


h=(b-a)/(npt-1)
DO i = 1, npt
   xpt(i)=a+(i-1)*h
   IF(2*(i/2).eq.i)THEN
     wht(i)=Four*h/Three
   ELSE
     wht(i)=two*h/Three
   ENDIF
ENDDO
wht(1)=h/Three
wht(npt)=h/Three

IF(ABS(SUM(wht)-(Npt-1)*h).gt.Ten*REpsilon)THEN
  WRITE(Out_Unit,*)'Error in Simpson weights'
  WRITE(Out_Unit,*)ABS(SUM(wht))
  WRITE(Out_Unit,*)(Npt-1)*h
  STOP 'Error in Simpson'
ENDIF

RETURN
END SUBROUTINE Simpson