subroutine besks (xnu, x, nin, bk)
  c july 1980 edition.   w. fullerton, c3, los alamos scientific lab.
        dimension bk(1)
        external  r1mach
        data xmax / 0.0 /
  c
        if (xmax.ne.0.0) go to 10
        xmax = -alog (r1mach(1))
        xmax = xmax + 0.5*alog(3.14*0.5/xmax)
  c
   10   if (x.gt.xmax) call seteru (
       1  36hbesks   x so big bessel k underflows, 36, 1, 2)
  c
        call beskes (xnu, x, nin, bk)
  c
        expxi = exp (-x)
        n = iabs (nin)
        do 20 i=1,n
          bk(i) = expxi * bk(i)
   20   continue
  c
        return
        end