subroutine dbesks (xnu, x, nin, bk)
  c july 1980 edition.   w. fullerton, c3, los alamos scientific lab.
        double precision xnu, x, bk(1), expxi, xmax, d1mach, dexp, dlog
        external d1mach
        data xmax / 0.d0 /
  c
        if (xmax.ne.0.0d0) go to 10
        xmax = -dlog(d1mach(1))
        xmax = xmax + 0.5d0*dlog(3.14d0*0.5d0/xmax)
  c
   10   if (x.gt.xmax) call seteru (
       1  36hdbesks  x so big bessel k underflows, 36, 1, 2)
  c
        call dbskes (xnu, x, nin, bk)
  c
        expxi = dexp (-x)
        n = iabs (nin)
        do 20 i=1,n
          bk(i) = expxi * bk(i)
   20   continue
  c
        return
        end