      FUNCTION YNTERP ( X, Y, XI, N, ND )
C LINEAR INTERPOLATION
C ND: DUMMY, TO BE COMPATIBLE WITH CIT LIBRARY ROUTINE
C FLAT OUTSIDE THE RANGE OF X(I), I.E. YNTERP=Y(1) OR Y(N)
C X(I) MUST BE MONOTONIC, EITHER DECREASING OR INCREASING
C X(I) MUST BE ALL DIFFERENT
C HIROO KANAMORI AUG.4 1978
       DIMENSION X(1), Y(1)
      I=1
      J=N
      K=1
      IF ( X(1) .GT. X(N) ) GO TO 11
1     IF( XI-X(K)) 2,2,3
2     J=K
      GO TO 4
3     I=K
4     K=(I+J)/2
      IF(K.NE.I) GO TO 1
      IF( XI.GT.X(N)) J=N+1
      GO TO 21
11    IF( X(K)-XI) 12,12,13
12    J=K
      GO TO 14
13    I=K
14    K=(I+J)/2
      IF(K.NE.I) GO TO 11
      IF( XI.LT. X(N)) J=N+1
21    IF ( J.EQ. 1) GO TO 23
      IF(J.EQ.N+1) GO TO 24
      YNTERP=Y(J-1)+(Y(J)-Y(J-1))*(XI-X(J-1))/(X(J)-X(J-1))
      RETURN
23    YNTERP=Y(1)
      RETURN
24    YNTERP=Y(N)
      RETURN
      END
      SUBROUTINE DTRD (X, N, Y, NDEG, AVEX)
      DIMENSION X(1),Y(1)
      AN = N
      S1 = 0.0
      S2 = 0.0
      IF(NDEG) 11, 11, 20
   11 DO 12 I=1,N
   12 S1 = S1+X(I)
      AVEX = S1/AN
      DO 15 I=1,N
   15 Y(I) = X(I) - AVEX
      RETURN
   20 DO 22 I=1,N
      S1 = S1+X(I)
   22 S2 = S2 + S1
      AVEI = 0.5 * (AN+1.0)
      AVEX = S1/AN
      SLOPE = -12.0*(S2-AVEI*S1)/(AN*(AN**2-1.0))
      DO 28 I=1,N
      AI = I
   28 Y(I) = X(I)-AVEX-SLOPE*(AI-AVEI)
      RETURN
      END
      SUBROUTINE TAPER (A, N, START, END, B)
      DIMENSION A(1),B(1)
      AN = N
      M1 = AN*START+0.5
      M2 = M1 + 1
      IF(M1) 10, 10, 11
   11 ANG = 3.1415926 / FLOAT(M1)
      DO 12 I=1, M1
      XI = I-1
      CS = (1.0-COS(XI*ANG))/2.0
   12 B(I) = A(I)*CS
   10 M3 = AN*END+0.5
      M5 = N-M3
      M4 = M5 + 1
      IF(M3) 13, 13, 14
   14 ANG = 3.1415926 / FLOAT (M3)
      DO 15 I=M4,N
      XI = I-N
      CS = (1.0-COS(XI*ANG))/2.0
   15 B(I) = A(I)*CS
   13 DO 16 I=M2, M5
   16 B(I) = A(I)
      RETURN
      END
      SUBROUTINE COOL(NN,DATAI,SIGNI)
      DIMENSION DATAI(1)
      N=2**(NN+1)
      J=1
      DO 5 I=1,N,2
      IF(I-J)1,2,2
    1 TEMPR=DATAI(J)
      TEMPI=DATAI(J+1)
      DATAI(J)=DATAI(I)
      DATAI(J+1)=DATAI(I+1)
      DATAI(I)=TEMPR
      DATAI(I+1)=TEMPI
    2 M=N/2
    3 IF(J-M)5,5,4
    4 J=J-M
      M=M/2
      IF(M-2)5,3,3
    5 J=J+M
      MMAX=2
    6 IF(MMAX-N)7,10,10
    7 ISTEP=2*MMAX
      THETA=SIGNI*6.28318531/FLOAT(MMAX)
      SINTH=SIN(THETA/2.)
      WSTPR=-2.0  *SINTH*SINTH
      WSTPI= SIN(THETA)
      WR=1.
      WI=0.
      DO 9 M=1,MMAX,2
      DO 8 I=M,N,ISTEP
      J=I+MMAX
      TEMPR=WR*DATAI(J)-WI*DATAI(J+1)
      TEMPI=WR*DATAI(J+1)+WI*DATAI(J)
      DATAI(J)=DATAI(I)-TEMPR
      DATAI(J+1)=DATAI(I+1)-TEMPI
      DATAI(I)=DATAI(I)+TEMPR
    8 DATAI(I+1)=DATAI(I+1)+TEMPI
      TEMPR=WR
      WR=WR*WSTPR-WI*WSTPI+WR
    9 WI=WI*WSTPR+TEMPR*WSTPI+WI
      MMAX=ISTEP
      GO TO 6
10    RETURN
      END
      SUBROUTINE  ORDRJX ( X, J, N )
C     ORDRJX REARRANGES ORIGINAL DATA IN ACCORDANCE WITH THEIR ARRIVAL TIMES OF
C     P WAVE
      DIMENSION X(200),J(200)
      J(1) = 1
      IF( N .EQ. 1 )  GO TO 320
      I = 2
  310 IF( I .GT. N )  GO  TO 320
      XTST = X(I)
      IM1 = I-1
      DO  300 JJ = 1, IM1
      K = I - JJ
      TST= X(K) - XTST
      IF( TST .GT. 0.0 )  GO TO  300
      KST = K + 1
      GO TO 330
  300 CONTINUE
      KST = K
  330 C=XTST
      IC = I
      DO 340 L=KST, I
      B = X(L)
      X(L) = C
      C = B
      IB = J(L)
      J(L) = IC
      IC = IB
  340 CONTINUE
      I = I + 1
      GO TO 310
  320 RETURN
      END
