C =============================================================================

C MULTIPLE SCATTERING SOURCE FUNCTION INTERPOLATION SET-UP ROUTINE.
C USING CUBIC SPLINE INTERPOLATION ON THE LOG IN THE RADIAL DIRECTION.
C NOTE THAT THE ARRAYS ARE PADDED OUT TO INCLUDE THE LOWERMOST AND
C UPPERMOST RADIAL BOUNDARIES.  THE MULTIPLE SCATTERING SOURCE FUNCTION
C AT THOSE POINTS IS ESTIMATED USING SIMPLE ``SCALE HEIGHT''
C EXTRAPOLATION (I.E., LINEAR EXTRAPOLATION ON THE LOGARITHMS).
C PASSES BACK:
C   XSMS  . . . PADDED RADIAL INTERPOLATION ARRAY.
C   D0SMS . . . 2-D ARRAY OF THE LOGARITHMS OF THE MULTIPLE SCATTERING
C               SOURCE FUNCTIONS AT THE POINTS DEFINED BY THE XSMS
C               AND CHIN ARRAYS.
C   D2SMS . . . SECOND (RADIAL) DERIVATIVES OF THE MULTIPLE SCATTERING
C               SOURCE FUNCTION LOGARITHMS AT THE POINTS DEFINED BY
C               THE XSMS AND CHIN ARRAYS.

      SUBROUTINE SMS_SPLINE(IKNT,JKNT,BRAD,RADN,SOL,SRC,
     &                      XSMS,D0SMS,D2SMS)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (NMAX = 100)      ! must be at least as large as IKNT+2
      DIMENSION BRAD(IKNT+1),RADN(IKNT),SOL(IKNT,JKNT),SRC(IKNT,JKNT)
      DIMENSION XSMS(IKNT+2),D0SMS(IKNT+2,JKNT),D2SMS(IKNT+2,JKNT)
      DIMENSION YSMS(NMAX),D2YSMS(NMAX)
      COMMON /NMBR/  PI,RTPI,PID2,OFFSET

c set up array of radial points
      DO 101 ISPLN = 1,IKNT
        XSMS(ISPLN+1) = RADN(ISPLN)
 101  END DO
c . . . lowermost point:  bottom boundary radius
      XSMS(1) = BRAD(1)
c . . . uppermost point:  top boundary radius
      XSMS(IKNT+2) = BRAD(IKNT+1)

c set up array of logarithms of multiple scattering source function
      DO 201 JJ = 1,JKNT
        DO 202 ISPLN = 1,IKNT
c currently have SRC & SOL archived to the LYAO_source.DAT file retaining
c five significant digits.  possibility arises of SRC-SOL=0 in optically
c thin situations, leading to NaN results.  protection trigger is set
c one decade down (current value of OFFSET is 1.0D-6). (JBishop 10Jul96)
          DELSMS = (SRC(ISPLN,JJ) - SOL(ISPLN,JJ)) / SRC(ISPLN,JJ)
          IF (DELSMS.GT.OFFSET) THEN
            D0SMS(ISPLN+1,JJ) = LOG(SRC(ISPLN,JJ) - SOL(ISPLN,JJ))
          ELSE
            D0SMS(ISPLN+1,JJ) = LOG(OFFSET * SRC(ISPLN,JJ))
          END IF
 202    END DO
c . . . lowermost point:  extrapolate to bottom boundary radius
        slope = (D0SMS(3,JJ)-D0SMS(2,JJ)) / (XSMS(3)-XSMS(2))
        D0SMS(1,JJ) = D0SMS(2,JJ) + (XSMS(1)-XSMS(2))*SLOPE
c . . . uppermost point:  extrapolate to top boundary radius
        slope = (D0SMS(IKNT+1,JJ)-D0SMS(IKNT,JJ)) /
     &          (XSMS(IKNT+1)-XSMS(IKNT))
        D0SMS(IKNT+2,JJ) = D0SMS(IKNT+1,JJ) +
     &                     (XSMS(IKNT+2)-XSMS(IKNT+1))*SLOPE
 201  END DO

c now to break into individual 1-D arrays, obtain second derivatives,
c then repack into 2-D array
      DO 301 JJ = 1,JKNT
        DO 302 ISPLN = 1,IKNT+2
          YSMS(ISPLN) = D0SMS(ISPLN,JJ)
 302    END DO
        CALL SPLINE_LYAO(XSMS,YSMS,IKNT+2,1.0D30,1.0D30,D2YSMS)
        DO 303 ISPLN = 1,IKNT+2
          D2SMS(ISPLN,JJ) = D2YSMS(ISPLN)
 303    END DO
 301  END DO

      RETURN
      END

C -----------------------------------------------------------------------------

C SPLINE SUBROUTINE FROM NUMERICAL RECIPES, MODIFIED TO DOUBLE PRECISION
      SUBROUTINE SPLINE_LYAO(X,Y,N,YP1,YPN,Y2)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (NMAX=100)
      DIMENSION X(N),Y(N),Y2(N),U(NMAX)
      IF (YP1.GT.0.99D30) THEN
        Y2(1)=0.0D0
        U(1) =0.0D0
      ELSE
        Y2(1)=-0.5D0
        U(1) =(3.0D0/(X(2)-X(1)))*((Y(2)-Y(1))/(X(2)-X(1))-YP1)
      ENDIF
      DO 11 I=2,N-1
        SIG=(X(I)-X(I-1))/(X(I+1)-X(I-1))
        P  =SIG*Y2(I-1)+2.0D0
        Y2(I)=(SIG-1.0D0)/P
        U(I) =(6.0D0*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1))
     *       /(X(I)-X(I-1)))/(X(I+1)-X(I-1))-SIG*U(I-1))/P
11    CONTINUE
      IF (YPN.GT.0.99D30) THEN
        QN=0.0D0
        UN=0.0D0
      ELSE
        QN=0.50D0
        UN=(3.0D0/(X(N)-X(N-1)))*(YPN-(Y(N)-Y(N-1))/(X(N)-X(N-1)))
      ENDIF
      Y2(N)=(UN-QN*U(N-1))/(QN*Y2(N-1)+1.0D0)
      DO 12 K=N-1,1,-1
        Y2(K)=Y2(K)*Y2(K+1)+U(K)
12    CONTINUE
      RETURN
      END

C -----------------------------------------------------------------------------

      SUBROUTINE SPLINT_LYAO(XA,YA,Y2A,N,X,Y)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION XA(N),YA(N),Y2A(N)
      KLO=1
      KHI=N
1     IF (KHI-KLO.GT.1) THEN
        K=(KHI+KLO)/2
        IF(XA(K).GT.X)THEN
          KHI=K
        ELSE
          KLO=K
        ENDIF
      GOTO 1
      ENDIF
      H=XA(KHI)-XA(KLO)
      IF (H.EQ.0.0D0) PAUSE 'Bad XA input.'
      A=(XA(KHI)-X)/H
      B=(X-XA(KLO))/H
      Y=A*YA(KLO)+B*YA(KHI)+
     *      ((A**3-A)*Y2A(KLO)+(B**3-B)*Y2A(KHI))*(H**2)/6.0D0
      RETURN
      END

C =============================================================================