#include "fpreproc/library_names.h"
#include "fpreproc/fp_const.h"
*DECK DBKISR
      SUBROUTINE DBKISR (X, N, SUM, IERR)
C***BEGIN PROLOGUE  DBKISR
C***SUBSIDIARY
C***PURPOSE  Subsidiary to DBSKIN
C***LIBRARY   SLATEC
C***TYPE      REAL*8 (BKISR-S, DBKISR-D)
C***AUTHOR  Amos, D. E., (SNLA)
C***DESCRIPTION
C
C     DBKISR computes repeated integrals of the K0 Bessel function
C     by the series for N=0,1, and 2.
C
C***SEE ALSO  DBSKIN
C***ROUTINES CALLED  D1MACH, DPSIXN
C***REVISION HISTORY  (YYMMDD)
C   820601  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890911  Removed unnecessary intrinsics.  (WRB)
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900328  Added TYPE section.  (WRB)
C   910722  Updated AUTHOR section.  (ALS)
C***END PROLOGUE  DBKISR
      IMPLICIT NONE
      INTEGER I, IERR, K, KK, KKN, K1, N, NP
      REAL*8 AK, ATOL, BK, C, FK, FN, HX, HXS, POL, PR, SUM,
     * TKP, TOL, TRM, X, XLN
      REAL*8 DPSIXN, D1MACH
      DIMENSION C(2)
      SAVE C
C
      DATA C(1), C(2) /CONST(1.57079632679489662,+00),CONST(1.0,0)/
C***FIRST EXECUTABLE STATEMENT  DBKISR
      IERR=0
      TOL = MAX(D1MACH(4),CONST(1.0,-18))
      IF (X.LT.TOL) GO TO 50
      PR = CONST(1.0,0)
      POL = CONST(0.0,0)
      IF (N.EQ.0) GO TO 20
      DO 10 I=1,N
        POL = -POL*X + C(I)
        PR = PR*X/I
   10 CONTINUE
   20 CONTINUE
      HX = X*CONST(0.5,0)
      HXS = HX*HX
      XLN = LOG(HX)
      NP = N + 1
      TKP = CONST(3.0,0)
      FK = CONST(2.0,0)
      FN = N
      BK = CONST(4.0,0)
      AK = CONST(2.0,0)/((FN+CONST(1.0,0))*(FN+CONST(2.0,0)))
      SUM = AK*(DPSIXN(N+3)-DPSIXN(3)+DPSIXN(2)-XLN)
      ATOL = SUM*TOL*CONST(0.75,0)
      DO 30 K=2,20
        AK = AK*(HXS/BK)*((TKP+CONST(1.0,0))/(TKP+FN+CONST(1.0,0)))*
     +      (TKP/(TKP+FN))
        K1 = K + 1
        KK = K1 + K
        KKN = KK + N
        TRM = (DPSIXN(K1)+DPSIXN(KKN)-DPSIXN(KK)-XLN)*AK
        SUM = SUM + TRM
        IF (ABS(TRM).LE.ATOL) GO TO 40
        TKP = TKP + CONST(2.0,0)
        BK = BK + TKP
        FK = FK + CONST(1.0,0)
   30 CONTINUE
      GO TO 80
   40 CONTINUE
      SUM = (SUM*HXS+DPSIXN(NP)-XLN)*PR
      IF (N.EQ.1) SUM = -SUM
      SUM = POL + SUM
      RETURN
C-----------------------------------------------------------------------
C     SMALL X CASE, X.LT.WORD TOLERANCE
C-----------------------------------------------------------------------
   50 CONTINUE
      IF (N.GT.0) GO TO 60
      HX = X*CONST(0.5,0)
      SUM = DPSIXN(1) - LOG(HX)
      RETURN
   60 CONTINUE
      SUM = C(N)
      RETURN
   80 CONTINUE
      IERR=2
      RETURN
      END
! 22Jun2000 fgtok -s r8_precision.sub "r8con.csh conversion"
! 22Jun2000 fgtok
! 22Jun2000 fgtok
