      SUBROUTINE HYTRH
C--GIVEN DEPTH & DISTANCE, THIS ROUTINE CALCULATES TRAVEL TIME, ITS
C--DERIVATIVES AND EMERGENCE ANGLES AT THE SOURCE FOR ALL ARRIVALS.
C  USES HOMOGENEOUS LAYER VELOCITY MODELS & ALGORITHMS OF HYPOLAYR/HYPO71
      INCLUDE 'common.inc'
      DIMENSION TINJ(13),DIDJ(13),TR(12)
      LOGICAL ALTMOD

C--THE FOLLOWING ARE PASSED THRU THE ARRAY A:
C DTDR      !TT DERIV WRT DISTANCE
C DTDZ      !TT DERIV WRT DEPTH
C T      !TRAVEL TIME
C AIN      !ANGLE OF EMERGENCE AT SOURCE

      ZSQ=Z1*Z1
      ALTMOD=MODALT(MOD).GT.0

C--LOOP OVER ALL ARRIVALS
      DO 280 I=1,M
C--FIND STATION INDEX AND REMOVE KPS AS AN S FLAG
      KI=IND(I)
      KPS=KI/10000
      KI=KI-10000*KPS
      J=KINDX(KI)

C--DETERMINE THE MODEL NO. TO ACTUALLY USE FOR THIS STATION
      MD=MOD
      IF (ALTMOD .AND. JLMOD(J)) MD=MODALT(MOD)

C--STATION DISTANCE
      DX=DIS(KI)

C--JL IS THE LAYER IN WHICH THE HYPOCENTER LIES
      NL=LAY(MD)
      DO 10 L=1,NL
        IF (D(L,MD).GT.Z1) THEN
          JJ=L
          JL=L-1
          GOTO 30
        END IF
10    CONTINUE
      JL=NL

C--CALC SOME BASIC PARAMETERS FOR THIS MODEL
30    TKJ=Z1-D(JL,MD)
      TKJSQ=TKJ**2+0.000001
C--IF HYPO IS IN HALFSPACE (BOT LAYER), THEN ALL RAYS GO UP
      IF (JL.EQ.NL) GOTO 100
      TMIN=999.99

C--DETERMINE DIST OF CLOSEST CRITICALLY REFRACTED RAY XOVMAX
      DO 40 L=JJ,NL
        SQT=SQRT(VSQ(L,MD)-VSQ(JL,MD))
        TINJ(L)=HYTID(JL,L,MD)-TKJ*SQT/(VEL(L,MD)*VEL(JL,MD))
40    DIDJ(L)=HYDID(JL,L,MD)-TKJ*VEL(JL,MD)/SQT
      XOVMAX=VEL(JJ,MD)*VEL(JL,MD)*(TINJ(JJ)-HYTID(JL,JL,MD))/
     2 (VEL(JJ,MD)-VEL(JL,MD))

      DO 70 L=JJ,NL
        TR(L)=TINJ(L)+DX/VEL(L,MD)
        IF (TR(L).GT.TMIN .OR. DIDJ(L).GT.DX) GOTO 70
        K=L
        TMIN=TR(L)
70    CONTINUE
      IF (DX .LT. XOVMAX) GO TO 90

C--TRAVEL TIME & DERIVATIVES FOR REFRACTED WAVE
80    T=TR(K)
      DTDR=1.0/VEL(K,MD)
      DTDZ=-SQRT(VSQ(K,MD)-VSQ(JL,MD))/(VEL(K,MD)*VEL(JL,MD))
      ANIN=-VEL(JL,MD)/VEL(K,MD)
      GO TO 260

C--DIRECT WAVE WITH EARTHQ IN TOP LAYER
90    IF (JL.NE.1) GOTO 100
      SQT=SQRT(ZSQ+DX**2)
      TDJ1=SQT/VEL(1,MD)
      IF (TDJ1 .GE. TMIN) GO TO 80

C--TRAVEL TIME & DERIVATIVES FOR DIRECT WAVE IN FIRST LAYER
      T=TDJ1
      DTDR=DX/(VEL(1,MD)*SQT)
      DTDZ=Z1/(VEL(1,MD)*SQT)
      ANIN=DX/SQT
      GO TO 260

C--FIND A DIRECT WAVE THAT WILL EMERGE AT THE STATION
100   XBIG=DX
      XLIT=DX*TKJ/Z1
      UB=XBIG/SQRT(XBIG**2+TKJSQ)
      UL=XLIT/SQRT(XLIT**2+TKJSQ)
      UBSQ=UB**2
      ULSQ=UL**2
      DELBIG=TKJ*UB/SQRT(1.000001-UBSQ)
      DELLIT=TKJ*UL/SQRT(1.000001-ULSQ)
      J1=JL-1
      DO 110 L=1,J1
        DELBIG=DELBIG+(THK(L,MD)*UB)/SQRT(VSQ(JL,MD)/VSQ(L,MD)-UBSQ)
110   DELLIT=DELLIT+(THK(L,MD)*UL)/SQRT(VSQ(JL,MD)/VSQ(L,MD)-ULSQ)

C--SHOOT UP TO 25 RAYS TO ITERATIVELY FIND THE STATION BY TRAPPING
C--IT BETWEEN TWO RAYS
      DO 170 LL=1,25
        XTEST=DELBIG-DELLIT
        IF (XTEST.LT..02) GOTO 180
        XTR=XLIT+(DX-DELLIT)*(XBIG-XLIT)/XTEST
        U=XTR/SQRT(XTR**2+TKJSQ)
        USQ=U**2
        DELXTR=TKJ*U/SQRT(1.000001-USQ)
        DO 120 L=1,J1
120     DELXTR=DELXTR+(THK(L,MD)*U)/SQRT(VSQ(JL,MD)/VSQ(L,MD)-USQ)
        XTEST=DX-DELXTR
C--STATION IS FOUND WHEN RAY IS WITHIN .02 KM
        IF (ABS(XTEST) .LE. 0.02) GOTO 190
        IF (XTEST.LT.0.) THEN
          XBIG=XTR
          DELBIG=DELXTR
        ELSE
          XLIT=XTR
          DELLIT=DELXTR
        END IF
        IF (LL.GT.10 .AND. U.GT..9999) GOTO 190
170   CONTINUE

180   XTR=0.5*(XBIG+XLIT)
      U=XTR/SQRT(XTR**2+TKJSQ)
      USQ=U**2

C--IF U IS TOO NEAR 1, COMPUTE TDIR AS WAVE ALONG THE TOP OF LAYER JL
190   IF (U.GT..9999) THEN
        TDC=HYTID(JL,JL,MD)+DX/VEL(JL,MD)
        IF (JL.NE.NL .AND. TDC.GE.TMIN) GO TO 80
        T=TDC
        DTDR=1./VEL(JL,MD)
        DTDZ=0.
        ANIN=0.99999
        GOTO 260
      END IF

C--TRAVEL TIME & DERIVATIVES FOR DIRECT WAVE BELOW FIRST LAYER
      TDIR=TKJ/(VEL(JL,MD)*SQRT(1.-USQ))
      DO 240 L=1,J1
240   TDIR=TDIR+(THK(L,MD)*VEL(JL,MD))/(VSQ(L,MD)*SQRT(VSQ(JL,MD)/
     2 VSQ(L,MD)-USQ))
      IF(JL.NE.NL .AND. TDIR.GE.TMIN) GOTO 80
      T=TDIR
      SRR=SQRT(1.-USQ)
      SRT=SRR**3
      ALFA=TKJ/SRT
      BETA=TKJ*U/(VEL(JL,MD)*SRT)
      DO 250 L=1,J1
        STK=(SQRT(VSQ(JL,MD)/VSQ(L,MD)-USQ))**3
        VTK=THK(L,MD)/(VSQ(L,MD)*STK)
        ALFA=ALFA+VTK*VSQ(JL,MD)
250   BETA=BETA+VTK*VEL(JL,MD)*U
      DTDR=BETA/ALFA
      DTDZ=(1.-VEL(JL,MD)*U*DTDR)/(VEL(JL,MD)*SRR)
      ANIN=U
260   AIN=RDEG*ASIN(ANIN)
      IF (AIN.LT.0.) AIN=AIN+180.
      AIN=180.-AIN

C--END OF STATION LOOP
270   A(I,1)=AIN
      A(I,2)=T
      A(I,3)=DTDR
280   A(I,4)=DTDZ

      RETURN
      END

C**************************************************
      FUNCTION HYDID (J,ML,MD)
C--CALLED BY HYTRH IN PROGRAM HYPOINVERSE
      INCLUDE 'common.inc'
      HYDID=0.
      IF (ML.EQ.1) GOTO 110
      DO 20 L=1,ML-1
        SQT=SQRT(VSQ(ML,MD)-VSQ(L,MD))
        DM=THK(L,MD)*VEL(L,MD)/SQT
        FLJ=1.
        IF (L.GE.J) FLJ=2.
        HYDID=HYDID+FLJ*DM
20    CONTINUE
110   RETURN
      END

C*********************************************************
      FUNCTION HYTID (J,ML,MD)
C--CALLED BY HYTRH IN PROGRAM HYPOINVERSE
      INCLUDE 'common.inc'
      HYTID=0.
      IF (ML.EQ.1) GOTO 110
      DO 20 L=1,ML-1
        SQT=SQRT(VSQ(ML,MD)-VSQ(L,MD))
        TIM=THK(L,MD)*SQT/(VEL(L,MD)*VEL(ML,MD))
        FLJ=1.
        IF (L.GE.J) FLJ=2.
        HYTID=HYTID+FLJ*TIM
20    CONTINUE
110   RETURN
      END
