      SUBROUTINE HYTRT
C--GIVEN DEPTH AND DISTANCE, THIS ROUTINE CALCULATES TRAVEL TIME AND ITS
C--DERIVATIVES, AND EMERGENCE ANGLE AT THE SOURCE FOR ALL ARRIVALS.
C--HYTRT USES A CONDENSED AND REDUCED TRAVEL TIME TABLE GENERATED 
C--BY THE PROGRAM TTGEN.
C--IND AND KINDX ARE INDEXING ARRAYS NECESSARY TO REFERENCE THE CORRECT
C--ELEMENTS OF THE STATION ARRAYS.

      INCLUDE 'common.inc'
C--ALTMOD CONTROLS THE BEHAVIOR OF HYTRT. IF F, HYTRT CALCS TT DATA FOR ALL
C  STATIONS USING MODEL MOD. IF T, HYTRT CALCS TT'S ONLY FOR STATIONS WITH
C  A DESIGNATED ALTERNATE MODEL MODALT(MOD) (IF JLMOD(J)=T)
      LOGICAL ALTMOD

C--THESE VALUES ARE PASSED BY THE ARRAY A.
C DTDR      !TT DERIVATIVE WRT DISTANCE
C DTDZ      !TT DERIVATIVE WRT DEPTH
C T      !TRAVEL TIME
C AIN      !ANGLE OF EMERGENCE AT THE SOURCE, MEASURED IN DEGREES FROM NADIR.

C--WE DO NOT NEED THESE FOR SEVERAL MODELS AT ONCE, SINCE THE SAME
C  MODEL IS USED FOR ALL STATIONS.
C      DIMENSION VH(LN)      !VELOCITY AT HYPOCENTER
C      DIMENSION DHRZ(LN)      !DIST AT WHICH HORIZ RAY EMERGES
C      DIMENSION TZ(LN,42)      !TT AT EACH DISTANCE POINT
C      DIMENSION DTZ(LN,42)      !TT DEPTH DERIV AT EACH DISTANCE POINT

      DIMENSION TZ(42)      !TT AT EACH DISTANCE POINT
      DIMENSION DTZ(42)      !TT DEPTH DERIV AT EACH DISTANCE POINT
      SAVE SCFAC
      DATA SCFAC/2000./

C--ALLOW 2 PASSES OVER THE CALCULATIONS FOR MAIN & OPTIONAL ALTERNATE MODEL
      ALTMOD=.FALSE.
      MD=MOD
      DO 70 IPASS=1,2
      IF (IPASS.EQ.2) THEN
        MD=MODALT(MOD)
        IF (MD.EQ.0) THEN
          RETURN
        ELSE
          ALTMOD=.TRUE.
        END IF
      END IF

C--PERFORM SOME PRELIMINARY CALCS AND INTERPOLATIONS WHICH DEPEND ONLY ON
C  DEPTH. DO THIS FOR EACH LINEAR GRAD MODEL ONCE.
      NZM=NZ(MD)
      NZ1M=NZ1(MD)
      NZ2M=NZ2(MD)
      DZ1M=DZ1(MD)
      DZ2M=DZ2(MD)
      LAYM=LAY(MD)
      ND1M=ND1(MD)

C--PERFORM DEPTH INTERPOLATION FIRST
C--FIND DEPTH INDEX = INDEX OF TABLE ENTRY NEAREST HYPOCENTER
C--ALSO DEPTH SPACING H AND FRACTION OF INTERVAL X
C--USE I (NEAREST POINT) FOR 3-POINT INTERPOLATION
C--USE I1 (POINT ABOVE) FOR 2-POINT INTERPOLATION
      TEMP=NZ1M*DZ1M
      IF ((.NOT.GZ1(MD) .OR. Z1.GT.TEMP) .AND. GD2(MD)) THEN

C--HYPO IN LOWER PART OF TABLE
        TMP2=(Z1-TEMP)/DZ2M
        I=NZ1M+TMP2+1.5001
        IF (I.GE.NZM) I=NZM-1
        IF (I.LT.NZ1M+2) I=NZ1M+2
        H=DZ2M
        X=TMP2-(I-NZ1M-1)
        I1=NZ1M+TMP2+1
        IF (I1.GE.NZM) I1=NZM-1
        Y=TMP2-(I1-NZ1M-1)
      ELSE

C--HYPO IN UPPER PART OF TABLE
        TMP2=Z1/DZ1M
        I=TMP2+1.5
        IF (I.GT.NZ1M) I=NZ1M
        IF (I.LT.2) I=2
        H=DZ1M
        X=TMP2-(I-1)
        I1=TMP2+1
        Y=TMP2-(I1-1)
      END IF

C--FIND EXACT VELOCITY AT HYPOCENTER. HYPO IS IN LAYER L.
      DO 20 K=1,LAYM
        IF (D(K,MD).LT.Z1) L=K
20    CONTINUE
      VH=VEL(LAYM,MD)
      IF (L.LT.LAYM) VH=VEL(L,MD)+(VEL(L+1,MD)-VEL(L,MD))
     2 *(Z1-D(L,MD))/(D(L+1,MD)-D(L,MD))

C--INTERPOLATE DIST AT WHICH A HORIZONTAL RAY EMERGES
C      DHRZ=(KDHR(MD,I)+X*(KDHR(MD,I+1)-KDHR(MD,I)))*.1
      DHRZ=(KDHR(MD,I1)+Y*(KDHR(MD,I1+1)-KDHR(MD,I1)))*.1

C--DEPTH INTERPOLATION
      IF (Z1.LE.TEMP+DZ2M*NZ2M) THEN
C--USE 3 POINT INTERPOLATION
        CA=X*.5*(X-1.)
        CB=1.-X**2
        CC=X*.5*(X+1.)
        DA=(X-.5)/H
        DB=-2.*X/H
        DC=(X+.5)/H
      ELSE
C--USE LINEAR EXTRAPOLATION
        CA=0.
        CB=1.-X
        CC=X
        DA=0.
        DC=1./H
        DB=-DC
        DHRZ=1000.
      END IF

C--INTERPOLATE TT AND ITS DEPTH DERIVATIVE FOR ALL DISTANCE GRID POINTS
      TEMP=ND1M*DD1(MD)
      DO 30 J=1,ND(MD)
        DX=(J-1)*DD1(MD)
        IF (J.GT.ND1M+1) DX=TEMP+(J-ND1M-1)*DD2(MD)
        TZ(J)=(CA*KT(MD,I-1,J)+CB*KT(MD,I,J)
     2  +CC*KT(MD,I+1,J)+32000.)/SCFAC+DX*REDV(MD)
30    DTZ(J)=(DA*KT(MD,I-1,J)+DB*KT(MD,I,J)+DC*KT(MD,I+1,J))/SCFAC

C--NOW FIND TRAVEL TIME, DTDR, DTDZ, AND ANGLE OF EMERGENCE FOR EACH ARRIVAL
C--LOOP OVER ALL ARRIVALS
      DO 60 I=1,M

C--FIND STATION INDEX 
      KI=IND(I)
      KPS=KI/10000
      KI=KI-10000*KPS
      J=KINDX(KI)

C--DECIDE WHETHER TO USE AN ALTERNATE MODEL FOR THIS STATION
C  (SKIP STATION IF WE NEED BUT DONT HAVE AN ALTERNATE CALCULATION)
      IF (ALTMOD .AND. .NOT.JLMOD(J)) GOTO 60

C--START DISTANCE INTERPOLATION
      NDM=ND(MD)
      ND1M=ND1(MD)
      ND2M=ND2(MD)
      DD1M=DD1(MD)
      DD2M=DD2(MD)
      TEMP=ND1M*DD1M

C--STATION DISTANCE
      DX=DIS(KI)

C--FIND DISTANCE INDEX = INDEX OF TABLE ENTRY NEAREST THE STATION DIST
C--ALSO DISTANCE SPACING H AND FRACTION OF INTERVAL X
      IF ((.NOT.GD1(MD) .OR. DX.GT.TEMP) .AND. GD2(MD)) THEN
C--HYPO IN FAR PART OF TABLE
        J=ND1M+(DX-TEMP)/DD2M+1.5001
        IF (J.GE.NDM) J=NDM-1
        IF (J.LT.ND1M+2) J=ND1M+2
        H=DD2M
        X=(DX-TEMP)/DD2M-(J-ND1M-1)
      ELSE
C--HYPO IN NEAR PART OF TABLE
        J=DX/DD1M+1.5
        IF (J.GT.ND1M) J=ND1M
        IF (J.LT.2) J=2
        H=DD1M
        X=DX/DD1M-(J-1)
      END IF

C--DISTANCE INTERPOLATION
      IF (DX.LE.TEMP+DD2M*ND2M) THEN
C--USE 3 POINT INTERPOLATION
        CA=X*.5*(X-1.)
        CB=1.-X**2
        CC=X*.5*(X+1.)
        DA=X-.5
        DB=-2.*X
        DC=X+.5
      ELSE
C--USE LINEAR EXTRAPOLATION
        CA=0.
        CB=1.-X
        CC=X
        DA=0.
        DC=1.
        DB=-1.
      END IF

C--INTERPOLATE TT AND ITS 2 DERIVATIVES
      A(I,2)=CA*TZ(J-1)+CB*TZ(J)+CC*TZ(J+1)
      A(I,4)=CA*DTZ(J-1)+CB*DTZ(J)+CC*DTZ(J+1)
      DTDR=(DA*TZ(J-1)+DB*TZ(J)+DC*TZ(J+1))/H
      IF (DTDR.LT.0.) DTDR=0.
      A(I,3)=DTDR

C--CALC EMERGENCE ANGLE
      CA=VH*DTDR
      IF (CA.GT..99999) CA=.99999
      AIN=RDEG*ASIN(CA)
      IF (DX.LT.DHRZ) AIN=180.-AIN
      A(I,1)=AIN

C--END OF LOOP OVER ARRIVALS
60    CONTINUE
70    CONTINUE
      RETURN
      END
