      SUBROUTINE SPBESL (LMAXA,XX,FJ)
C *** SEE R H LANDAU, PROGRAM LPOTT, 1981
C     SPBESL CALCULATES SPHERICAL BESSEL FUCTIONS OF 1ST      KIND FOR A
C     FOR ALL ORDER FROM ZERO TO LMAXA, AT ARGUMRNT XX
      IMPLICIT REAL*8  (A-H,O-Z)
      DIMENSION FJ(50)
      LMAX = LMAXA
      R = XX
      X = XX
      LX = 2.*R
      IF (49..LT.LX) LX = 49.
C     CALCULATE SMALL ARGUMENT LIMIT
      IF (R.GT..01) GO TO 10
      FJ0 = 1.-X**2/6.+X**4/120.-X**6/42./120.
      IF (LMAX.GT.1) GO TO 90
      FJ(1) = FJ0
      FJ(2) = X*(1.-X**2/10.+X**4/280.-X**6/30./56./9.)/3.
      GO TO 130
C     REGULAR VALUE OF X, START WITH LARGE L,ITERATE DOWN,RENORMALIZE
 10   FJ0 = SIN(X)/X
 20   L2 = MAX0(LMAX+5,LX)
      XL2 = 2.*L2
      FJ(L2+1) = 1.E-10
      FJ(L2) = 1.E-10*((XL2+1.)/X-X/(XL2+3.))
      L3 = L2-1
      DO 40 LL=1,L3
         L1 = L2-LL
         FL1 = L1
         FJ(L1) = (2.*FL1+1.)*FJ(L1+1)/X-FJ(L1+2)
         RHL = ABS(FJ(L1))
         IF (RHL.LE.1.E30) GO TO 40
C     FJ TOO LARGE, SCALE DOWN AND/OR SET TOO SMALL FJ*S TO 0
         DO 30 L4=L1,L2
            IF ((ABS(FJ(L4)).LT.1.E-26).AND.(L4.NE.1)) FJ(L4) = 0.
C     NB CHANGED 50 TO 26 FOR VAX
            FJ(L4) = 1.E-10*FJ(L4)
 30      CONTINUE
 40   CONTINUE
C     APPLAY RENORMALIZ FACTOR
      ZZJ = FJ0/FJ(1)
      DO 60 L1=1,2
         IF (FJ(L1).EQ.0.) GO TO 50
         IF (LOG10(ABS(ZZJ))+LOG10(ABS(FJ(L1))).LE.-60.) FJ(L1) = 0.
 50      FJ(L1) = ZZJ*FJ(L1)
 60   CONTINUE
      IF (LMAX.LE.1) GO TO 130
      L2 = L2-4
      DO 80 L1=3,L2
         IF (FJ(L1).EQ.0.) GO TO 70
         IF (LOG10(ABS(ZZJ))+LOG10(ABS(FJ(L1))).LE.-60.) FJ(L1) = 0.
 70      FJ(L1) = ZZJ*FJ(L1)
 80   CONTINUE
      GO TO 130
 90   IF (R) 120,100,20
C     X=0 VALUES
 100  WRITE (6,140)
      L2 = LMAX+1
      DO 110 LL=2,L2
         FJ(LL) = 0.
 110  CONTINUE
      FJ(1) = FJ0
      GO TO 130
 120  WRITE (6,150)
 130  CONTINUE
      RETURN
C
 140  FORMAT (26H SPBESL(X) CALLED FOR X=0,)
 150  FORMAT (31H ERROR IN SPBESL, NEG. ARGUMENT)
      END
