SUBROUTINE GlgnZro ( n, x )
USE Numeric_Kinds_Module
USE Precision_Module
USE FileUnits_Module
!
!-------------------- GlgnZro -----------------------------------------
!
! GlgnZro -- zeroes of Guass Legendre polynomials
!
!  input --
!
!   n -- order of the Guass Legendre polynomial whose zeroes are desired
!        (h0 has NONE, h1 has one zero, h25 has 25 zeroes)
!
!  output --
!
!     x -- array listing the non-negative zeroes of nth polynomial
!          x will be returned with 13 numbers in it IF n=25
!           note-- dimension statements limit n=200
!
!-----------------------------------------------------------------------
IMPLICIT NONE
SAVE
CHARACTER(LEN=21), PARAMETER    :: ProcName='GlgnZro'
INTEGER(KIND=IW_Kind), PARAMETER:: MaxDim=201
LOGICAL                         :: ODD, LAST
INTEGER(KIND=IW_Kind)           :: i, ifpr, iodd, iroot, it, n, nnr, np, npr
REAL(KIND=WP_Kind)              :: xg(MaxDim), hp(MaxDim), x(*)
REAL(KIND=WP_Kind)              :: deriv, cor, theta, xz, xzn
!
!------------------- execution begins here ----------------------------
!
IF( n .lt. 1 ) RETURN
!
x(1) = Zero
IF( n .eq. 1 ) RETURN
!
IF(n+1.gt.MaxDim)THEN
   WRITE(Out_Unit,*)'Error in GlgnZro :',n+1, MaxDim
   WRITE(Msg_Unit,*)'Error in GlgnZro :',n+1, MaxDim
   STOP 'GlgnZro'
ENDIF

npr = n/2
np = n + 1
iodd = mod(n,2)
ODD = iodd .eq. 1
nnr = iodd + npr
ifpr = iodd + 1
IF( ODD ) x(1) = Zero
!
DO iroot = ifpr, nnr
   theta = Pi*(Four*(npr+iroot+1)-Five)/(Four*n+Two)
   theta = theta + One/(Eight*n**2*TAN(theta))
   xg(iroot) = -COS(theta)
ENDDO
!
DO i = ifpr, nnr
   xz = xg(i)
   CALL Pn ( xz, n, hp)
   
   it = 0
   LAST = .false.
!
 11      CONTINUE
      it = it + 1
      deriv = (n*(-xz*hp(np)+hp(n)))/(One-xz*xz)
      cor = hp(np)/deriv
      xzn = xz - cor
!
      IF( LAST ) GOTO 12
      !LAST = One+ Half*ABS(cor/xzn) .eq. One
      LAST=ABS(Half*ABS(cor/xzn)).lt.Ten*REpsilon
      IF( it .gt. 30 ) GOTO 13
!
      CALL Pn ( xzn, n, hp)
      xz = xzn
   GOTO 11
!
 13      CONTINUE
   WRITE (Out_Unit, 1000 ) it, iroot, n, xzn, xz
   WRITE (Msg_Unit, 1000 ) it, iroot, n, xzn, xz
   STOP 'GlgnZro'
!
 12      CONTINUE
   x(i) = xzn
ENDDO
!
 1000 FORMAT( ' poor convergence after ', i2, ' iterations *** GlgnZro -- root(order) = '  &
             , i3,"(",i3,")",/,' best zero = ', e20.10, ', previous guess = ', e20.10 )
!
RETURN
END SUBROUTINE GlgnZro