o prevent overflows, the scaling approach of Holmes and Featherstone
       (2002) is utilized. These functions are accurate to about degree 2800.
       The index of the array corresponding to a given degree L and angular
       order M can be computed by a call to <B>PlmIndex</B>.

       The integral of the squared Legendre functions over the interval [-1,
       1] is 2 * (2 - delta(0,m))/(2L+1), where delta is the Kronecker delta
       function. If the optional parameter CNORM is set equal to 1, the
 