SUBROUTINE LagZro ( n, x )
USE Numeric_Kinds_Module
!
!-------------------- LagZro -----------------------------------------
!
! LagZro -- zeroes of Laguere polynomials
!
! === Latest Revision:    | 900524 at 15:38 | fixed up the convergence
!
!  input --
!
!     n -- order of the Laguere 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
!
!---------------------------------------------------------------------
!
USE FileUnits_Module

IMPLICIT NONE
!
!
LOGICAL   LAST
INTEGER(KIND=IW_Kind)   n, np, it, i, ibisect
!
REAL(KIND=WP_Kind) xg(300), hp(300), x(n), lower(300), upper(300)
REAL(KIND=WP_Kind) beta, beta8, beta8sq, bj, max, delta, xl, xu, xz, cor, abs, xzn, sqrt, hl, hu

IF( n .lt. 1 ) RETURN
!
x(1) = One
IF( n .eq. 1 ) RETURN
!
np = n + 1
!
DO i = 1, n
   beta = (i-Fourth)*pi
   beta8 = Eight*beta
   beta8sq = beta8**2
   bj = beta + (One + (-124.e0_WP_Kind/Three + (120928.e0_WP_Kind/15.e0_WP_Kind    &
       -401743168.e0_WP_Kind/105.e0_WP_Kind/beta8sq)/beta8sq)/beta8sq)/ beta8
   xg(i) = bj**2*(One+(-Two+bj**2)/(48.e0_WP_kind*(n+Half)**2))/ (Four*(n+Half))
   lower(i) = bj**2/(Four*(n+Half))
   upper(i) = (i+Half)/(n+Half)*(Two*(i+Half)+sqrt(Four*(i+Half)**2+Fourth))
ENDDO
!
DO 10 i = 1, n
!
! Make lower bound greater than previous root.
!
   IF(i.gt.2) THEN
      IF(lower(i).lt.x(i-1)+Half*(x(i-1)-x(i-2))/2) THEN
         lower(i) = x(i-1)+Half*(x(i-1)-x(i-2))/2
         IF(xg(i).lt.lower(i))xg(i)=(lower(i)+upper(i))/2
      ENDIF
   ENDIF
   xz = xg(i)
   CALL Laguer_fun ( xz, n, hp, np )
   it = 0
   LAST = .false.
!
 11      CONTINUE
      it = it + 1
      cor = xz/n*hp(np)/(hp(np)-hp(n))

      IF(it.eq.1.or.ABS(cor).gt.Half*(upper(i)-lower(i)).or.xz-cor.lt.lower(i).or.xz-cor.gt.upper(i))THEN
         delta = (upper(i)-lower(i))/10
         DO 30 ibisect = 1, 10
            xl = lower(i) + delta*(ibisect-1)
            CALL Laguer_fun (xl, n, hp, np)
            hl = hp(np)
            xu = lower(i) + delta*ibisect
            CALL Laguer_fun (xu, n, hp, np)
            hu = hp(np)
            IF(hl*hu.lt.Zero) THEN
               xzn = (xl*hu-xu*hl)/(hu-hl)
               cor = xz - xzn
               lower(i) = xl
               upper(i) = xu
               GOTO 31
            ENDIF
   30          CONTINUE
   31          CONTINUE
      ENDIF

      xzn = xz - cor

!
      IF( LAST ) GOTO 12
      LAST = One + .01e0_WP_Kind*ABS(cor/max(One,xzn)) .eq. One
      IF( it .gt. 40 ) GOTO 13
!
      CALL Laguer_fun ( xzn, n, hp, np )
      xz = xzn
   GOTO 11
!
 13      CONTINUE
   WRITE (Msg_Unit, 1000 ) it, i, n, xzn, xz
   STOP 'lagzro'
!
 12      CONTINUE
   x(i) = xzn
 10   CONTINUE
!
 1000 FORMAT( ' poor convergence after ', i2, ' iterations ***' ' LagZro -- root(order) = ', i3,"(",i3,")",/,   &
               ' best zero = ', e20.10, ', previous guess = ', e20.10 )
!
RETURN
END