lamv Subroutine

   subroutine lamv(v,xi,Vm,Vl,Dl)

      real(wp),intent(in) :: v !! Order of lambda function
      real(wp),intent(in) :: xi !! Argument of lambda function
      real(wp),intent(out) :: Vm !! Highest order computed
      real(wp),dimension(0:*),intent(out) :: Vl !! Lambda function of order `n+v0`
      real(wp),dimension(0:*),intent(out) :: Dl !! Derivative of lambda function

      real(wp) :: a0 , bjv0 , bjv1 , bk , ck , cs , f , f0 ,  &
                  f1 , f2 , fac , ga , px , qx , r , r0 , rc ,&
                  rp, x, rq , sk , uk , v0 , vk , vv , x2 , xk
      integer :: i , j , k , k0 , m , n

      real(wp),parameter :: rp2 = 2.0_wp / pi ! 0.63661977236758d0

      x = abs(xi)
      x2 = x*x
      n = int(v)
      v0 = v - n
      Vm = v
      if ( x<=12.0_wp ) then
         do k = 0 , n
            vk = v0 + k
            bk = 1.0_wp
            r = 1.0_wp
            do i = 1 , 50
               r = -0.25_wp*r*x2/(i*(i+vk))
               bk = bk + r
               if ( abs(r)<abs(bk)*1.0e-15_wp ) exit
            enddo
            Vl(k) = bk
            uk = 1.0_wp
            r = 1.0_wp
            do i = 1 , 50
               r = -0.25_wp*r*x2/(i*(i+vk+1.0_wp))
               uk = uk + r
               if ( abs(r)<abs(uk)*1.0e-15_wp ) exit
            enddo
            Dl(k) = -0.5_wp*x/(vk+1.0_wp)*uk
         enddo
         return
      endif
      k0 = 11
      if ( x>=35.0_wp ) k0 = 10
      if ( x>=50.0_wp ) k0 = 8
      bjv0 = 0.0_wp
      bjv1 = 0.0_wp
      do j = 0 , 1
         vv = 4.0_wp*(j+v0)*(j+v0)
         px = 1.0_wp
         rp = 1.0_wp
         do k = 1 , k0
            rp = -0.78125e-2_wp*rp*(vv-(4.0_wp*k-3.0_wp)**2.0_wp) &
                 *(vv-(4.0_wp*k-1.0_wp)**2.0_wp)/(k*(2.0_wp*k-1.0_wp)*x2)
            px = px + rp
         enddo
         qx = 1.0_wp
         rq = 1.0_wp
         do k = 1 , k0
            rq = -0.78125e-2_wp*rq*(vv-(4.0_wp*k-1.0_wp)**2.0_wp) &
                 *(vv-(4.0_wp*k+1.0_wp)**2.0_wp)/(k*(2.0_wp*k+1.0_wp)*x2)
            qx = qx + rq
         enddo
         qx = 0.125_wp*(vv-1.0_wp)*qx/x
         xk = x - (0.5_wp*(j+v0)+0.25_wp)*pi
         a0 = sqrt(rp2/x)
         ck = cos(xk)
         sk = sin(xk)
         if ( j==0 ) bjv0 = a0*(px*ck-qx*sk)
         if ( j==1 ) bjv1 = a0*(px*ck-qx*sk)
      enddo
      if ( v0==0.0_wp ) then
         ga = 1.0_wp
      else
         call gam0(v0,ga)
         ga = v0*ga
      endif
      fac = (2.0_wp/x)**v0*ga
      Vl(0) = bjv0
      Dl(0) = -bjv1 + v0/x*bjv0
      Vl(1) = bjv1
      Dl(1) = bjv0 - (1.0_wp+v0)/x*bjv1
      r0 = 2.0_wp*(1.0_wp+v0)/x
      if ( n<=1 ) then
         Vl(0) = fac*Vl(0)
         Dl(0) = fac*Dl(0) - v0/x*Vl(0)
         Vl(1) = fac*r0*Vl(1)
         Dl(1) = fac*r0*Dl(1) - (1.0_wp+v0)/x*Vl(1)
         return
      endif
      if ( n>=2 .and. n<=int(0.9_wp*x) ) then
         f0 = bjv0
         f1 = bjv1
         do k = 2 , n
            f = 2.0_wp*(k+v0-1.0_wp)/x*f1 - f0
            f0 = f1
            f1 = f
            Vl(k) = f
         enddo
      elseif ( n>=2 ) then
         m = msta1(x,200)
         if ( m<n ) then
            n = m
         else
            m = msta2(x,n,15)
         endif
         f = 0.0_wp
         f2 = 0.0_wp
         f1 = 1.0e-100_wp
         do k = m , 0 , -1
            f = 2.0_wp*(v0+k+1.0_wp)/x*f1 - f2
            if ( k<=n ) Vl(k) = f
            f2 = f1
            f1 = f
         enddo
         cs = 0.0_wp
         if ( abs(bjv0)>abs(bjv1) ) then
            cs = bjv0/f
         else
            cs = bjv1/f2
         endif
         do k = 0 , n
            Vl(k) = cs*Vl(k)
         enddo
      endif
      Vl(0) = fac*Vl(0)
      do j = 1 , n
         rc = fac*r0
         Vl(j) = rc*Vl(j)
         Dl(j-1) = -0.5_wp*x/(j+v0)*Vl(j)
         r0 = 2.0_wp*(j+v0+1)/x*r0
      enddo
      Dl(n) = 2.0_wp*(v0+n)*(Vl(n-1)-Vl(n))/x
      Vm = n + v0

   end subroutine lamv