legendre_roots Subroutine

      subroutine legendre_roots(l,root,w)
  !-----------------------------------------------------------------------------
  ! Subroutine finds the l roots (in theta) and gaussian weights associated
  ! with the legendre polynomial of degree l > 1. (Source from SCRIP create
  ! gaussian grid routine
  !-----------------------------------------------------------------------------
  implicit none

  !-----------------------------------------------------------------------------
  !   arguments
  !-----------------------------------------------------------------------------
  integer, intent(in) :: l

  real(ESMF_KIND_R8), dimension(l), intent(out) :: root, w

  !-----------------------------------------------------------------------------
  !   local variables
  !-----------------------------------------------------------------------------
  integer :: l1, l2, l22, l3, k, i, j
  real(ESMF_KIND_R8) :: del,co,p1,p2,p3,t1,t2,slope,s,c,pp1,pp2,p00

  !-----------------------------------------------------------------------------
  !   Define useful constants.
  !-----------------------------------------------------------------------------
  del= pi/float(4*l)
  l1 = l+1
  co = float(2*l+3)/float(l1**2)
  p2 = 1.0
  t2 = -del
  l2 = l/2
  k = 1
  p00 = one/sqrt(two)

  !-----------------------------------------------------------------------------
  !   Start search for each root by looking for crossing point.
  !-----------------------------------------------------------------------------
  do i=1,l2
 10  t1 = t2
     t2 = t1+del
     p1 = p2
     s = sin(t2)
     c = cos(t2)
     pp1 = 1.0
     p3 = p00
     do j=1,l1
        pp2 = pp1
        pp1 = p3
        p3 = 2.0*sqrt((float(j**2)-0.250)/float(j**2))*c*pp1-                  &
             sqrt(float((2*j+1)*(j-1)*(j-1))/float((2*j-3)*j*j))*pp2
     enddo
     p2 = pp1
     if ((k*p2).gt.0) goto 10

  !-----------------------------------------------------------------------------
  !      Now converge using Newton-Raphson.
  !-----------------------------------------------------------------------------
     k = -k
 20  continue
     slope = (t2-t1)/(p2-p1)
     t1 = t2
     t2 = t2-slope*p2
     p1 = p2
     s = sin(t2)
     c = cos(t2)
     pp1 = 1.0
     p3 = p00
     do j=1,l1
        pp2 = pp1
        pp1 = p3
        p3 = 2.0*sqrt((float(j**2)-0.250)/float(j**2))*c*pp1-                  &
             sqrt(float((2*j+1)*(j-1)*(j-1))/float((2*j-3)*j*j))*pp2
     enddo
     p2 = pp1
     if (abs(p2).gt.1.e-10) goto 20
     root(i) = t2
     w(i) = co*(sin(t2)/p3)**2
  enddo

  !-----------------------------------------------------------------------------
  !   If l is odd, take care of odd point.
  !-----------------------------------------------------------------------------
  l22 = 2*l2
  if (l22 .ne. l) then
     l2 = l2+1
     t2 = pi/2.0
     root(l2) = t2
     s = sin(t2)
     c = cos(t2)
     pp1 = 1.0
     p3 = p00
     do j=1,l1
        pp2 = pp1
        pp1 = p3
        p3 = 2.0*sqrt((float(j**2)-0.250)/float(j**2))*c*pp1-                  &
             sqrt(float((2*j+1)*(j-1)*(j-1))/float((2*j-3)*j*j))*pp2
     enddo
     p2 = pp1
     w(l2) = co/p3**2
  endif

  !-----------------------------------------------------------------------------
  !   Use symmetry to compute remaining roots and weights.
  !-----------------------------------------------------------------------------
  l3 = l2+1
  do i=l3,l
     root(i) = pi-root(l-i+1)
     w(i) = w(l-i+1)
  enddo

  return

  !-----------------------------------------------------------------------------
      end subroutine legendre_roots