lqmn Subroutine

   subroutine lqmn(Mm,m,n,x,Qm,Qd)

      real(wp),intent(in) :: x !! Argument of Qmn(x)
      integer,intent(in) :: m !! Order of Qmn(x)  ( m = 0,1,2,… )
      integer,intent(in) :: n !! Degree of Qmn(x) ( n = 0,1,2,… )
      integer,intent(in) :: mm !! Physical dimension of QM and QD
      real(wp),intent(out) :: Qm(0:Mm,0:n) !! Qmn(x)
      real(wp),intent(out) :: Qd(0:Mm,0:n) !! Qmn'(x)

      integer :: i , j , k , km , ls
      real(wp) :: q0 , q1 , q10 , qf , qf0 , qf1 , qf2 , xq , xs

      if ( abs(x)==1.0_wp ) then
         do i = 0 , m
            do j = 0 , n
               Qm(i,j) = 1.0e+300_wp
               Qd(i,j) = 1.0e+300_wp
            enddo
         enddo
         return
      endif
      ls = 1
      if ( abs(x)>1.0_wp ) ls = -1
      xs = ls*(1.0_wp-x*x)
      xq = sqrt(xs)
      q0 = 0.5_wp*log(abs((x+1.0_wp)/(x-1.0_wp)))
      if ( abs(x)<1.0001_wp ) then
         Qm(0,0) = q0
         Qm(0,1) = x*q0 - 1.0_wp
         Qm(1,0) = -1.0_wp/xq
         Qm(1,1) = -ls*xq*(q0+x/(1.0_wp-x*x))
         do i = 0 , 1
            do j = 2 , n
               Qm(i,j) = ((2.0_wp*j-1.0_wp)*x*Qm(i,j-1)-(j+i-1.0_wp) &
                          *Qm(i,j-2))/(j-i)
            enddo
         enddo
         do j = 0 , n
            do i = 2 , m
               Qm(i,j) = -2.0_wp*(i-1.0_wp)*x/xq*Qm(i-1,j) &
                         - ls*(j+i-1.0_wp)*(j-i+2.0_wp)*Qm(i-2,j)
            enddo
         enddo
      else
         if ( abs(x)>1.1_wp ) then
            km = 40 + m + n
         else
            km = (40+m+n)*int(-1.0_wp-1.8_wp*log(x-1.0_wp))
         endif
         qf2 = 0.0_wp
         qf1 = 1.0_wp
         qf0 = 0.0_wp
         do k = km , 0 , -1
            qf0 = ((2*k+3.0_wp)*x*qf1-(k+2.0_wp)*qf2)/(k+1.0_wp)
            if ( k<=n ) Qm(0,k) = qf0
            qf2 = qf1
            qf1 = qf0
         enddo
         do k = 0 , n
            Qm(0,k) = q0*Qm(0,k)/qf0
         enddo
         qf2 = 0.0_wp
         qf1 = 1.0_wp
         do k = km , 0 , -1
            qf0 = ((2*k+3.0_wp)*x*qf1-(k+1.0_wp)*qf2)/(k+2.0_wp)
            if ( k<=n ) Qm(1,k) = qf0
            qf2 = qf1
            qf1 = qf0
         enddo
         q10 = -1.0_wp/xq
         do k = 0 , n
            Qm(1,k) = q10*Qm(1,k)/qf0
         enddo
         do j = 0 , n
            q0 = Qm(0,j)
            q1 = Qm(1,j)
            do i = 0 , m - 2
               qf = -2.0_wp*(i+1)*x/xq*q1 + (j-i)*(j+i+1.0_wp)*q0
               Qm(i+2,j) = qf
               q0 = q1
               q1 = qf
            enddo
         enddo
      endif
      Qd(0,0) = ls/xs
      do j = 1 , n
         Qd(0,j) = ls*j*(Qm(0,j-1)-x*Qm(0,j))/xs
      enddo
      do j = 0 , n
         do i = 1 , m
            Qd(i,j) = ls*i*x/xs*Qm(i,j) + (i+j)*(j-i+1.0_wp)/xq*Qm(i-1,j)
         enddo
      enddo
   end subroutine lqmn