SUBROUTINE pleg( x,      &   ! Input
                     nleg,   &   ! Input
                     p       )   ! Output

    !# returns the l=0..(n-1) Pl(x) legendre polynomials (all evaluated at x).
    !# x is a scalar
    !# uses Bonnet's recursion equation for legendre polynomials

    !#--------------------------------------------------------------------------#
    !#                         -- Type declarations --                          #
    !#--------------------------------------------------------------------------#

    !------------
    !  Arguments
    !------------
    REAL    x
    INTEGER nleg
    REAL    p( nleg )
    
    INTEGER i

    P( 1 ) = 1.0
    P( 2 ) = x
    i_nleg_loop: DO i = 3, nleg
      p( i ) = ( ( 2 * i - 3 ) * x * p( i - 1 ) - ( i - 2 ) * p( i - 2 ) ) / ( i - 1 )
    END DO i_nleg_loop
    
    END SUBROUTINE pleg