kmn Subroutine

   subroutine kmn(m,n,c,Cv,Kd,Df,Dn,Ck1,Ck2)

      integer :: m
      integer :: n
      real(wp) :: c
      integer :: Kd
      real(wp),dimension(200) :: Df
      real(wp),dimension(200) :: Dn
      real(wp) :: Ck1
      real(wp) :: Ck2

      real(wp) :: cs , Cv , dnp , g0 , gk0 , gk1 , gk2 , gk3 , &
                  r , r1 , r2 , r3 , r4 , r5 , rk(200) , sa0 , sb0 , &
                  su0 , sw , t , tp(200) , u(200) , v(200) , w(200)
      integer :: i , ip , j , k , l , nm , nm1 , nn

      nm = 25 + int(0.5_wp*(n-m)+c)
      nn = nm + m
      cs = c*c*Kd
      ip = 1
      if ( n-m==2*int((n-m)/2) ) ip = 0
      k = 0
      do i = 1 , nn + 3
         if ( ip==0 ) k = -2*(i-1)
         if ( ip==1 ) k = -(2*i-3)
         gk0 = 2.0_wp*m + k
         gk1 = (m+k)*(m+k+1.0_wp)
         gk2 = 2.0_wp*(m+k) - 1.0_wp
         gk3 = 2.0_wp*(m+k) + 3.0_wp
         u(i) = gk0*(gk0-1.0_wp)*cs/(gk2*(gk2+2.0_wp))
         v(i) = gk1 - Cv + (2.0_wp*(gk1-m*m)-1.0_wp)*cs/(gk2*gk3)
         w(i) = (k+1.0_wp)*(k+2.0_wp)*cs/((gk2+2.0_wp)*gk3)
      enddo
      do k = 1 , m
         t = v(m+1)
         do l = 0 , m - k - 1
            t = v(m-l) - w(m-l+1)*u(m-l)/t
         enddo
         rk(k) = -u(k)/t
      enddo
      r = 1.0_wp
      do k = 1 , m
         r = r*rk(k)
         Dn(k) = Df(1)*r
      enddo
      tp(nn) = v(nn+1)
      do k = nn - 1 , m + 1 , -1
         tp(k) = v(k+1) - w(k+2)*u(k+1)/tp(k+1)
         if ( k>m+1 ) rk(k) = -u(k)/tp(k)
      enddo
      if ( m==0 ) dnp = Df(1)
      if ( m/=0 ) dnp = Dn(m)
      Dn(m+1) = (-1)**ip*dnp*cs/((2.0_wp*m-1.0_wp)*(2.0_wp*m+1.0-4.0_wp*ip)*tp(m+1))
      do k = m + 2 , nn
         Dn(k) = rk(k)*Dn(k-1)
      enddo
      r1 = 1.0_wp
      do j = 1 , (n+m+ip)/2
         r1 = r1*(j+0.5_wp*(n+m+ip))
      enddo
      nm1 = (n-m)/2
      r = 1.0_wp
      do j = 1 , 2*m + ip
         r = r*j
      enddo
      su0 = r*Df(1)
      sw = 0.0_wp
      do k = 2 , nm
         r = r*(m+k-1.0_wp)*(m+k+ip-1.5_wp)/(k-1.0_wp)/(k+ip-1.5_wp)
         su0 = su0 + r*Df(k)
         if ( k>nm1 .and. abs((su0-sw)/su0)<1.0e-14_wp ) exit
         sw = su0
      enddo
      if ( Kd/=1 ) then
         r2 = 1.0_wp
         do j = 1 , m
            r2 = 2.0_wp*c*r2*j
         enddo
         r3 = 1.0_wp
         do j = 1 , (n-m-ip)/2
            r3 = r3*j
         enddo
         sa0 = (2.0_wp*(m+ip)+1.0_wp)*r1/(2.0_wp**n*c**ip*r2*r3*Df(1))
         Ck1 = sa0*su0
         if ( Kd==-1 ) return
      endif
      r4 = 1.0_wp
      do j = 1 , (n-m-ip)/2
         r4 = 4.0_wp*r4*j
      enddo
      r5 = 1.0_wp
      do j = 1 , m
         r5 = r5*(j+m)/c
      enddo
      g0 = Dn(m)
      if ( m==0 ) g0 = Df(1)
      sb0 = (ip+1.0_wp)*c**(ip+1)/(2.0_wp*ip*(m-2.0_wp)+1.0_wp)/(2.0_wp*m-1.0_wp)
      Ck2 = (-1)**ip*sb0*r4*r5*g0/r1*su0

   end subroutine kmn