text: function plm(l,m,x)
  !----------------------------------------------------
  ! Subroutine:
  !   plm
  !
  ! Purpose:
  !  Calculate assiociated Legendre polynomials      
  !----------------------------------------------------
    implicit none
    integer,parameter::dp=selected_real_kind(8) !double precision
    integer::l,m
    integer::i,ll
    real(kind=dp)::x,fact,pll,pmm,pmmp1,somx2
    real(kind=dp)::plm
    IF((m.lt.0).or.(m.gt.l).or.(abs(x).gt.1.D0))THEN
        WRITE(*,*) 'bad arguments in plm'
    ENDIF
    pmm=1.D0
    if(m.gt.0)then
       somx2=sqrt((1.D0-x)*(1.D0+x))
       fact=1.D0
       do 11 i=1,m
          pmm=-pmm*fact*somx2
          fact=fact+2.D0
  11   continue
    endif
    if(l.eq.m)then
       plm=pmm
    else
       pmmp1=x*dble(2*m+1)*pmm
       if(l.eq.m+1)then
          plm=pmmp1
       else
          do 12 ll=m+2,l
             pll=(x*dble(2*ll-1)*pmmp1-dble(ll+m-1)*pmm)/dble(ll-m)
             pmm=pmmp1
             pmmp1=pll
  12      continue
          plm=pll
       endif
    endif
    return
  end function plm
  