rmn2sp Subroutine

   subroutine rmn2sp(m,n,c,x,Cv,Df,Kd,R2f,R2d)

      integer,intent(in) :: m
      integer,intent(in) :: n
      real(wp) :: c
      real(wp) :: x
      real(wp) :: Cv
      real(wp) :: Df(200)
      integer :: Kd
      real(wp),intent(out) :: R2f
      real(wp),intent(out) :: R2d

      real(wp),dimension(0:251)  :: pd , pm , qd , qm
      real(wp),dimension(200)  :: dn
      real(wp) :: ck1 , ck2 , ga , gb , gc , r1 , r2 , &
                  r3 , r4, sd , sd0 , sd1 , sd2 , sdm , sf , spd1 , &
                  spd2 , spl , su0 , su1 , su2 , sum , sw
      integer :: ip , j , j1 , j2 , k , ki , l1 , nm , nm1 ,  &
                 nm2 , nm3

      real(wp),parameter :: eps = 1.0e-14_wp

      if ( abs(Df(1))<1.0e-280_wp ) then
         R2f = 1.0e+300_wp
         R2d = 1.0e+300_wp
         return
      endif

      ip = 1
      nm1 = int((n-m)/2)
      if ( n-m==2*nm1 ) ip = 0
      nm = 25 + nm1 + int(c)
      nm2 = 2*nm + m
      call kmn(m,n,c,Cv,Kd,Df,dn,ck1,ck2)
      call lpmns(m,nm2,x,pm,pd)
      call lqmns(m,nm2,x,qm,qd)
      su0 = 0.0_wp
      sw = 0.0_wp
      do k = 1 , nm
         j = 2*k - 2 + m + ip
         su0 = su0 + Df(k)*qm(j)
         if ( k>nm1 .and. abs(su0-sw)<abs(su0)*eps ) exit
         sw = su0
      enddo
      sd0 = 0.0_wp
      do k = 1 , nm
         j = 2*k - 2 + m + ip
         sd0 = sd0 + Df(k)*qd(j)
         if ( k>nm1 .and. abs(sd0-sw)<abs(sd0)*eps ) exit
         sw = sd0
      enddo
      su1 = 0.0_wp
      sd1 = 0.0_wp
      do k = 1 , m
         j = m - 2*k + ip
         if ( j<0 ) j = -j - 1
         su1 = su1 + dn(k)*qm(j)
         sd1 = sd1 + dn(k)*qd(j)
      enddo
      ga = ((x-1.0_wp)/(x+1.0_wp))**(0.5_wp*m)
      do k = 1 , m
         j = m - 2*k + ip
         if ( j<0 ) then
            if ( j<0 ) j = -j - 1
            r1 = 1.0_wp
            do j1 = 1 , j
               r1 = (m+j1)*r1
            enddo
            r2 = 1.0_wp
            do j2 = 1 , m - j - 2
               r2 = j2*r2
            enddo
            r3 = 1.0_wp
            sf = 1.0_wp
            do l1 = 1 , j
               r3 = 0.5_wp*r3*(-j+l1-1.0_wp)*(j+l1)/((m+l1)*l1)*(1.0_wp-x)
               sf = sf + r3
            enddo
            if ( m-j>=2 ) gb = (m-j-1.0_wp)*r2
            if ( m-j<=1 ) gb = 1.0_wp
            spl = r1*ga*gb*sf
            su1 = su1 + (-1)**(j+m)*dn(k)*spl
            spd1 = m/(x*x-1.0_wp)*spl
            gc = 0.5_wp*j*(j+1.0_wp)/(m+1.0_wp)
            sd = 1.0_wp
            r4 = 1.0_wp
            do l1 = 1 , j - 1
               r4 = 0.5_wp*r4*(-j+l1)*(j+l1+1.0_wp)/((m+l1+1.0_wp)*l1)*(1.0_wp-x)
               sd = sd + r4
            enddo
            spd2 = r1*ga*gb*gc*sd
            sd1 = sd1 + (-1)**(j+m)*dn(k)*(spd1+spd2)
         endif
      enddo
      su2 = 0.0_wp
      ki = (2*m+1+ip)/2
      nm3 = nm + ki
      do k = ki , nm3
         j = 2*k - 1 - m - ip
         su2 = su2 + dn(k)*pm(j)
         if ( j>m .and. abs(su2-sw)<abs(su2)*eps ) exit
         sw = su2
      enddo
      sd2 = 0.0_wp
      do k = ki , nm3
         j = 2*k - 1 - m - ip
         sd2 = sd2 + dn(k)*pd(j)
         if ( j>m .and. abs(sd2-sw)<abs(sd2)*eps ) exit
         sw = sd2
      enddo
      sum = su0 + su1 + su2
      sdm = sd0 + sd1 + sd2
      R2f = sum/ck2
      R2d = sdm/ck2

   end subroutine rmn2sp