sk_ku Subroutine

   subroutine sk_ku(xx, sk, ku)
   implicit none
   real(kind=R8), intent(in ), dimension(1:2) :: xx
   real(kind=R8), intent(out)                 :: sk
   real(kind=R8), intent(out)                 :: ku

      real(kind=R8)    :: xa, xb, mu, si, a, b, un
      real(kind=R8)    :: h, hh, b1, b2, alp, bet
      integer(kind=I4) :: ia, ib, deb, fin, npts, long, i, k

      real(kind=R8), dimension(1:2) :: x
      long  = nn

      do k = 1, 2
        x(k) = max( xx(k), 1.e-4_R8 )
      enddo

      ia = long
      ib = long

      npts = long*long

      deb = 1    +ia
      fin = npts -ib

      a = x(1)
      b = x(2)

      hh = (2._R8 -a -b)/(npts -1)
      h  = (pi_R8/2)*hh

      xa = a +ia*hh
      xb = b +ib*hh

      b1 = -pi_R8/2 *(1.-a)
      b2 = +pi_R8/2 *(1.-b)

      alp = -(b2-npts*b1)/(b2-b1)
      bet =      (npts-1)/(b2-b1)

      un = 1._R8
      !. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
      mu = log(1._R8/sin((pi_R8*xb)/2.)*sin((pi_R8*xa)/2.))
      mu = (un/h)*mu +add_tang(1, deb, fin, alp, bet, mu=0._R8, si=1._R8)
      do i = 1, ia-1
         mu = mu +tang(i*un, 1, alp, bet, mu=0._R8, si=1._R8)
      enddo
      do i = npts, npts -(ib -2), -1
         mu = mu +tang(i*un, 1, alp, bet, mu=0._R8, si=1._R8)
      enddo
      mu = mu/npts
      !. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
      si = (pi_R8*(-2 + xa + xb))/2. - (mu**2*pi_R8*(-2 + xa + xb))/2. + 1._R8/tan((pi_R8*xa)/2.) + 1._R8/tan((pi_R8*xb)/2.) - 2*mu*Log(1._R8/sin((pi_R8*xb)/2.)*Sin((pi_R8*xa)/2.))
      si = (un/h)*si +add_tang(2, deb, fin, alp, bet, mu, si=1._R8)
      do i = 1, ia-1
         si = si+ tang(i*un, 2, alp, bet, mu, si=1._R8)
      enddo
      do i = npts, npts -(ib -2), -1
         si = si +tang(i*un, 2, alp, bet, mu, si=1._R8)
      enddo
      si = si/npts
      si = sqrt(si)
      !. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
      sk = (6*mu*pi_R8 - 2*mu**3*pi_R8 - 3*mu*pi_R8*xa + mu**3*pi_R8*xa - 3*mu*pi_R8*xb + mu**3*pi_R8*xb - 6*mu*1._R8/tan((pi_R8*xa)/2.) - 6*mu*1._R8/tan((pi_R8*xb)/2.) - 1._R8/sin((pi_R8*xa)/2.)**2 + 1._R8/sin((pi_R8*xb)/2.)**2 - 2*Log(Sin((pi_R8*xa)/2.)) + 6*mu**2*Log(Sin((pi_R8*xa)/2.)) + 2*Log(Sin((pi_R8*xb)/2.)) - 6*mu**2*Log(Sin((pi_R8*xb)/2.)))/(2.*si**3)
      sk = (un/h)*sk +add_tang(3, deb, fin, alp, bet, mu, si)
      do i = 1, ia-1
         sk = sk +tang(i*un, 3, alp, bet, mu, si)
      enddo
      do i = npts, npts -(ib -2), -1
         sk = sk +tang(i*un, 3, alp, bet, mu, si)
      enddo
      sk = sk/npts
      !. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
      ku = (6*pi_R8 - 36*mu**2*pi_R8 + 6*mu**4*pi_R8 - 3*pi_R8*xa + 18*mu**2*pi_R8*xa - 3*mu**4*pi_R8*xa - 3*pi_R8*xb + 18*mu**2*pi_R8*xb - 3*mu**4*pi_R8*xb + 4*(-2 + 9*mu**2)*1._R8/tan((pi_R8*xa)/2.) + 4*(-2 + 9*mu**2)*1._R8/tan((pi_R8*xb)/2.) + 12*mu*1._R8/sin((pi_R8*xa)/2.)**2 - 12*mu*1._R8/sin((pi_R8*xb)/2.)**2 + 24*mu*Log(Sin((pi_R8*xa)/2.)) - 24*mu**3*Log(Sin((pi_R8*xa)/2.)) - 24*mu*Log(Sin((pi_R8*xb)/2.)) + 24*mu**3*Log(Sin((pi_R8*xb)/2.)) + 1._R8/sin((pi_R8*xa)/2.)**4*Sin(pi_R8*xa) + 1._R8/sin((pi_R8*xb)/2.)**4*Sin(pi_R8*xb))/(6.*si**4)
      ku = (un/h)*ku +add_tang(4, deb, fin, alp, bet, mu, si)
      do i = 1, ia-1
         ku = ku +tang(i*un, 4, alp, bet, mu, si)
      enddo
      do i = npts, npts -(ib -2), -1
         ku = ku +tang(i*un, 4, alp, bet, mu, si)
      enddo
      ku = ku/npts
      !. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

   return
   endsubroutine sk_ku