subroutine LegendreP(l, m, n, x, plm)

c***********************************************************************
c     LegendreP computes the associated legendre polynomial p_l^m
c     for the array x(0:n), and return the results in plm(0:n)
c***********************************************************************

      integer l, m, n, i
      double precision x(0:n), plm(0:n)
      double precision half, one, three

      parameter (half = 0.5d0, one = 1.0d0, three = 3.0d0)
c*********************************************************

      if (m .eq. 0) then
         if (l .eq. 0) then
            do 10000 i = 0, n
               plm(i) = one
10000       continue
         else if (l .eq. 1) then
            do 10010 i = 0, n
               plm(i) = x(i)
10010       continue
         else if (l .eq. 2) then
            do 10020 i = 0, n
               plm(i) = half*(three*x(i)*x(i)-one)
10020       continue
         else
            write(8, *)
     &         'LegendreP: l can only be 0, 1, 2 when m = 0, stop'
            stop
         endif
      else if (m .eq. 1) then
         if (l .eq. 1) then
            do 10110 i = 0, n
               plm(i) = -sqrt(one-x(i)*x(i))
10110       continue
         else if (l .eq. 2) then
            do 10120 i = 0, n
               plm(i) = -three*x(i)*sqrt(one-x(i)*x(i))
10120       continue
         else
            write(8, *)
     &         'LegendreP: l can only be 1, 2 when m = 1, stop'
            stop
         endif
      endif

      return
      end