LegendreQandL Subroutine

  subroutine LegendreQandL(N,x,q,qprime,lN)
    implicit none
    integer    :: N
    real(real64) :: x
    real(real64) :: lN,q,qprime
    ! Local
    real(real64) :: lNm1,lNm2,dlNm1,dlNm2,dlN,lN1,dlN1
    integer    :: i

    lNm2 = 1.0_real64
    lNm1 = x
    dlNm2 = 0.0_real64
    dlNm1 = 1.0_real64

    do i = 2,N

      lN = (2.0_real64*i-1.0_real64)/(real(i,real64))*x*lNm1-(real(i,real64)-1.0_real64)/(real(i,real64))*lNm2
      dlN = dlNm2+(2.0_real64*real(i,real64)-1.0_real64)*lNm1
      lNm2 = lNm1
      lNm1 = lN
      dlNm2 = dlNm1
      dlNm1 = dlN

    enddo

    i = N+1
    lN1 = (2.0_real64*i-1.0_real64)/(real(i,real64))*x*lN-(real(i,real64)-1.0_real64)/(real(i,real64))*lNm2
    dlN1 = dlNm2+(2.0_real64*real(i,real64)-1.0_real64)*lNm1
    q = lN1-lNm2
    qprime = dlN1-dlNm2

  endsubroutine LegendreQandL