C Grid travel time tracing: second-order method for the first arrivals
C in smooth models.  Testing 2-D version based on the NET program.
C
C Version: 3.00
C Date: 1996, June 18
C
C Author:
C     Ludek Klimes
C             Department of Geophysics, Charles University Prague
C             Ke Karlovu 3
C             121 16  Praha 2,  Czech Republic
C             E-mail: klimes@seis.karlov.mff.cuni.cz
C
C To be linked with 'net.for' version 3.00 or later.
C
C Reference:
C     Klimes, L. (1996):  Grid travel-time tracing: second-order method
C             for the first arrivals in smooth media.
C             PAGEOPH, 148, 539-563.
C
C=======================================================================
C
      SUBROUTINE TTTSRC()
C
C Auxiliary subroutine to SOURCE.
C
C Date: 1996, June 18
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Common blocks /GRID/ /FS/ /NODE/ /SRCC/ are required here:
      INCLUDE 'net.inc'
      INCLUDE 'netnode.inc'
C
C-----------------------------------------------------------------------
C
      INTEGER I1,I2,I
      REAL U2,UG1,UG2,AUX0,AUX1,AUX2
C
C.......................................................................
C
C     Annulling second slowness-vector components:
      DO 10 I=1,NL1*NL2*NL3
        RAM(IP20+I)=0.
   10 CONTINUE
C
C     Slowness gradient:
      IF(0.LT.IPOS1) THEN
        IF(IPOS1.LT.NL1-1) THEN
          UG1=(RAM(IP10+IPOS+1)-RAM(IP10+IPOS-1))/(DSX1+DSX1)
        ELSE
          UG1=(RAM(IP10+IPOS)  -RAM(IP10+IPOS-1))/ DSX1
        END IF
      ELSE
        IF(IPOS1.LT.NL1-1) THEN
          UG1=(RAM(IP10+IPOS+1)-RAM(IP10+IPOS)  )/ DSX1
        ELSE
          UG1=0.
        END IF
      END IF
      IF(0.LT.IPOS2) THEN
        IF(IPOS2.LT.NL2-1) THEN
          UG2=(RAM(IP10+IPOS+NL1)-RAM(IP10+IPOS-NL1))/(DSX2+DSX2)
        ELSE
          UG2=(RAM(IP10+IPOS)    -RAM(IP10+IPOS-NL1))/ DSX2
        END IF
      ELSE
        IF(IPOS1.LT.NL1-1) THEN
          UG2=(RAM(IP10+IPOS+NL1)-RAM(IP10+IPOS)    )/ DSX2
        ELSE
          UG2=0.
        END IF
      END IF
C
C     Loop over forward star of size 1:
      DO 22 I2=MAX0(-IPOS2,-1),MIN0(NL2-1-IPOS2,1)
        DO 21 I1=MAX0(-IPOS1,-1),MIN0(NL1-1-IPOS1,1)
            I=IPOS+I1+NL1*I2
            U2=RAM(IP10+I)
            AUX1=FLOAT(I1)*DSX1-DPOS1
            AUX2=FLOAT(I2)*DSX2-DPOS2
            AUX0=SQRT(AUX1*AUX1+AUX2*AUX2)
            RAM(ITT0+I)=TTI+0.5*(PI+U2)*AUX0
            IF(AUX0.NE.0.) THEN
C             AUX1=0.5*((PI+U2)*AUX1/AUX0+UG1*AUX0)
C             AUX2=0.5*((PI+U2)*AUX2/AUX0+UG2*AUX0)
              AUX0=AUX0*AUX0/(PI+U2)
              AUX1=AUX1+UG1*AUX0
              AUX2=AUX2+UG2*AUX0
              AUX0=SQRT(AUX1*AUX1+AUX2*AUX2)
              RAM(IP10+I)=U2*AUX1/AUX0
              RAM(IP20+I)=U2*AUX2/AUX0
            ELSE
              RAM(IP10+I)=U2
              RAM(IP20+I)=0.
            END IF
C           Updating queue by adding new node as Q(MAXQ)
            MAXQ=MAXQ+1
            IRAM(IPOSQ0+MAXQ)=I
          IF(I1.EQ.0.AND.I2.EQ.0) THEN
            RAM(ITT0+I)=-RAM(ITT0+I)
          END IF
   21   CONTINUE
   22 CONTINUE
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TTTUPD()
C
C This subroutine performs the loop over the grid-line-neighbour nodes.
C
C Auxiliary subroutine to GENER.
C
C Date: 1996, June 18
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Common blocks /GRID/ /FS/ /NODE/ are required here:
      INCLUDE 'net.inc'
      INCLUDE 'netnode.inc'
C
C-----------------------------------------------------------------------
C
      REAL GIANT
      PARAMETER (GIANT=1.E+10)
C
      INTEGER I,I1
      REAL TT0,P10,P20,TT1,P11,P21,PR2
C
C.......................................................................
C
C     Loop over gridline segments starting at node IPOS:
      TT0=RAM(ITT0+IPOS)
      P10=RAM(IP10+IPOS)
      P20=RAM(IP20+IPOS)
      IF(IPOS1.GT.0) THEN
        I1=IPOS-1
        TT1=RAM(ITT0+I1)
        P11=RAM(IP10+I1)
        P21=RAM(IP20+I1)
        IF(ABS(TT1).LE.TT0) THEN
          PR2=P11+P11-P10
          IF(IPOS1.GT.1) THEN
            IF(ABS(RAM(ITT0+IPOS-2)).LE.TT0) THEN
              PR2=RAM(IP10+IPOS-2)
            END IF
          END IF
          IF(IPOS2.GT.0) THEN
            I=I1  -NL1
            IF(RAM(ITT0+I).GE.GIANT) THEN
              MAXQ=MAXQ+1
              IRAM(IPOSQ0+MAXQ)=I
            END IF
            CALL TTT(TT1,P11,P21, DSX1,TT0,P10,P20,
     *         -DSX2,RAM(ITT0+I),RAM(IP10+I),RAM(IP20+I),-P10,-P11,-PR2)
            I=IPOS-NL1
            IF(RAM(ITT0+I).GE.GIANT) THEN
              MAXQ=MAXQ+1
              IRAM(IPOSQ0+MAXQ)=I
            END IF
            CALL TTT(TT0,P10,P20,-DSX1,TT1,P11,P21,
     *         -DSX2,RAM(ITT0+I),RAM(IP10+I),RAM(IP20+I),-P10,-P11,-PR2)
          END IF
          IF(IPOS2.LT.NL2-1) THEN
            I=I1  +NL1
            IF(RAM(ITT0+I).GE.GIANT) THEN
              MAXQ=MAXQ+1
              IRAM(IPOSQ0+MAXQ)=I
            END IF
            CALL TTT(TT1,P11,P21, DSX1,TT0,P10,P20,
     *          DSX2,RAM(ITT0+I),RAM(IP10+I),RAM(IP20+I),-P10,-P11,-PR2)
            I=IPOS+NL1
            IF(RAM(ITT0+I).GE.GIANT) THEN
              MAXQ=MAXQ+1
              IRAM(IPOSQ0+MAXQ)=I
            END IF
            CALL TTT(TT0,P10,P20,-DSX1,TT1,P11,P21,
     *          DSX2,RAM(ITT0+I),RAM(IP10+I),RAM(IP20+I),-P10,-P11,-PR2)
          END IF
        END IF
      END IF
      IF(IPOS1.LT.NL1-1) THEN
        I1=IPOS+1
        TT1=RAM(ITT0+I1)
        P11=RAM(IP10+I1)
        P21=RAM(IP20+I1)
        IF(ABS(TT1).LE.TT0) THEN
          PR2=P11+P11-P10
          IF(IPOS1.LT.NL1-2) THEN
            IF(ABS(RAM(ITT0+IPOS+2)).LE.TT0) THEN
              PR2=RAM(IP10+IPOS+2)
            END IF
          END IF
          IF(IPOS2.GT.0) THEN
            I=I1  -NL1
            IF(RAM(ITT0+I).GE.GIANT) THEN
              MAXQ=MAXQ+1
              IRAM(IPOSQ0+MAXQ)=I
            END IF
            CALL TTT(TT1,P11,P21,-DSX1,TT0,P10,P20,
     *         -DSX2,RAM(ITT0+I),RAM(IP10+I),RAM(IP20+I), P10, P11, PR2)
            I=IPOS-NL1
            IF(RAM(ITT0+I).GE.GIANT) THEN
              MAXQ=MAXQ+1
              IRAM(IPOSQ0+MAXQ)=I
            END IF
            CALL TTT(TT0,P10,P20, DSX1,TT1,P11,P21,
     *         -DSX2,RAM(ITT0+I),RAM(IP10+I),RAM(IP20+I), P10, P11, PR2)
          END IF
          IF(IPOS2.LT.NL2-1) THEN
            I=I1  +NL1
            IF(RAM(ITT0+I).GE.GIANT) THEN
              MAXQ=MAXQ+1
              IRAM(IPOSQ0+MAXQ)=I
            END IF
            CALL TTT(TT1,P11,P21,-DSX1,TT0,P10,P20,
     *          DSX2,RAM(ITT0+I),RAM(IP10+I),RAM(IP20+I), P10, P11, PR2)
            I=IPOS+NL1
            IF(RAM(ITT0+I).GE.GIANT) THEN
              MAXQ=MAXQ+1
              IRAM(IPOSQ0+MAXQ)=I
            END IF
            CALL TTT(TT0,P10,P20, DSX1,TT1,P11,P21,
     *          DSX2,RAM(ITT0+I),RAM(IP10+I),RAM(IP20+I), P10, P11, PR2)
          END IF
        END IF
      END IF
      IF(IPOS2.GT.0) THEN
        I1=IPOS-NL1
        TT1=RAM(ITT0+I1)
        P11=RAM(IP10+I1)
        P21=RAM(IP20+I1)
        IF(ABS(TT1).LE.TT0) THEN
          PR2=P21+P21-P20
          IF(IPOS2.GT.1) THEN
            IF(ABS(RAM(ITT0+IPOS-2*NL1)).LE.TT0) THEN
              PR2=RAM(IP20+IPOS-2*NL1)
            END IF
          END IF
          IF(IPOS1.GT.0) THEN
            I=I1  -1
            IF(RAM(ITT0+I).GE.GIANT) THEN
              MAXQ=MAXQ+1
              IRAM(IPOSQ0+MAXQ)=I
            END IF
            CALL TTT(TT1,P21,P11, DSX2,TT0,P20,P10,
     *         -DSX1,RAM(ITT0+I),RAM(IP20+I),RAM(IP10+I),-P20,-P21,-PR2)
            I=IPOS-1
            IF(RAM(ITT0+I).GE.GIANT) THEN
              MAXQ=MAXQ+1
              IRAM(IPOSQ0+MAXQ)=I
            END IF
            CALL TTT(TT0,P20,P10,-DSX2,TT1,P21,P11,
     *         -DSX1,RAM(ITT0+I),RAM(IP20+I),RAM(IP10+I),-P20,-P21,-PR2)
          END IF
          IF(IPOS1.LT.NL1-1) THEN
            I=I1  +1
            IF(RAM(ITT0+I).GE.GIANT) THEN
              MAXQ=MAXQ+1
              IRAM(IPOSQ0+MAXQ)=I
            END IF
            CALL TTT(TT1,P21,P11, DSX2,TT0,P20,P10,
     *          DSX1,RAM(ITT0+I),RAM(IP20+I),RAM(IP10+I),-P20,-P21,-PR2)
            I=IPOS+1
            IF(RAM(ITT0+I).GE.GIANT) THEN
              MAXQ=MAXQ+1
              IRAM(IPOSQ0+MAXQ)=I
            END IF
            CALL TTT(TT0,P20,P10,-DSX2,TT1,P21,P11,
     *          DSX1,RAM(ITT0+I),RAM(IP20+I),RAM(IP10+I),-P20,-P21,-PR2)
          END IF
        END IF
      END IF
      IF(IPOS2.LT.NL2-1) THEN
        I1=IPOS+NL1
        TT1=RAM(ITT0+I1)
        P11=RAM(IP10+I1)
        P21=RAM(IP20+I1)
        IF(ABS(TT1).LE.TT0) THEN
          PR2=P21+P21-P20
          IF(IPOS2.LT.NL2-2) THEN
            IF(ABS(RAM(ITT0+IPOS+2*NL1)).LE.TT0) THEN
              PR2=RAM(IP20+IPOS+2*NL1)
            END IF
          END IF
          IF(IPOS1.GT.0) THEN
            I=I1  -1
            IF(RAM(ITT0+I).GE.GIANT) THEN
              MAXQ=MAXQ+1
              IRAM(IPOSQ0+MAXQ)=I
            END IF
            CALL TTT(TT1,P21,P11,-DSX2,TT0,P20,P10,
     *         -DSX1,RAM(ITT0+I),RAM(IP20+I),RAM(IP10+I), P20, P21, PR2)
            I=IPOS-1
            IF(RAM(ITT0+I).GE.GIANT) THEN
              MAXQ=MAXQ+1
              IRAM(IPOSQ0+MAXQ)=I
            END IF
            CALL TTT(TT0,P20,P10, DSX2,TT1,P21,P11,
     *         -DSX1,RAM(ITT0+I),RAM(IP20+I),RAM(IP10+I), P20, P21, PR2)
          END IF
          IF(IPOS1.LT.NL1-1) THEN
            I=I1  +1
            IF(RAM(ITT0+I).GE.GIANT) THEN
              MAXQ=MAXQ+1
              IRAM(IPOSQ0+MAXQ)=I
            END IF
            CALL TTT(TT1,P21,P11,-DSX2,TT0,P20,P10,
     *          DSX1,RAM(ITT0+I),RAM(IP20+I),RAM(IP10+I), P20, P21, PR2)
            I=IPOS+1
            IF(RAM(ITT0+I).GE.GIANT) THEN
              MAXQ=MAXQ+1
              IRAM(IPOSQ0+MAXQ)=I
            END IF
            CALL TTT(TT0,P20,P10, DSX2,TT1,P21,P11,
     *          DSX1,RAM(ITT0+I),RAM(IP20+I),RAM(IP10+I), P20, P21, PR2)
          END IF
        END IF
      END IF
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TTT
     *         (TT0,P10,P20,DX1,TT1,P11,P21,DX2,TT2,P12,P22,PR0,PR1,PR2)
      REAL      TT0,P10,P20,DX1,TT1,P11,P21,DX2,TT2,P12,P22,PR0,PR1,PR2
C
C Subroutine calculating, in 2-D rectangular grid, travel time of the
C wave propagating from the segment between 2 given gridpoints Y0,Y1 to
C the third given gridpoint Y2.  Diffracted waves from the gridpoints
C Y0,Y1 are also calculated at Y2 and the first arrival is considered.
C Travel time is between points y0 and y1 interpolated by means of the
C Aitken-Hermite polynomial interpolation of the third order.
C
C Input:
C     DX1,DX2... Grid steps in temporary Cartesian coordinates X1,X2
C             defined as follows:
C             (0,0)...   Point Y0
C             (DX1,0)... Point Y1
C             (0,DX2)... Point Y2
C     TT0,P10,P20... Travel time and slowness vector at point Y0.
C     TT1,P11,P21... Travel time and slowness vector at point Y1.
C                    TT1 may have the minus sign which is ignored.
C     TT2,P12,P22... Travel time and slowness vector at point Y2.
C     PR0,PR1,PR2... First slowness components at points Y0,Y1,2*Y1-Y0
C                    for Y0.GT.Y1, or Y1,Y0,2*Y0-Y1 for Y1.GT.Y0.
C                    The sign corresponds to the coordinate increasing
C                    from PR0 to PR2.
C
C Output:
C     TT2,P12,P22... Travel time and slowness vector at point Y2
C             calculated from points Y0,Y1, if it is less than input
C             value of TT2.  Otherwise unchanged input.
C
C No subroutines and external functions required.
C
C Date: 1996, December 22
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Common block /TTTPAR/ is required here:
      INCLUDE 'ttt.inc'
C
C-----------------------------------------------------------------------
C
      REAL ABSTT0,ABSTT1,U0,U1,U2,U,UG1,UG2,X1,TT,T1,T2,T11,T12,T22
      REAL P1,P2,A0,A1,A01,A001,A011,A0A1,AUX0,AUX1,AUX2,AUX3
      REAL P10P20,P11P21,X1C,X2C,X2CX2C,UG1X2C,UC,R10,R11
C
C     ABSTT0..ABS(TT0)
C     ABSTT1..ABS(TT1)
C     U0,U1,U2,U...Slownesses at points Y0,Y1,Y2,Y.
C     UG1,UG2... Slowness gradient.
C     X1...   X1 coordinate of point Y.
C     TT,T1,T2,T11,T12,T22... Travel time and its first and second
C             derivatives at point Y.
C     TT,P1,P2... Travel time and its first derivatives at point Y2.
C     A0,A1...Relative position of point Y with respect to points Y0,Y1.
C     A01,A001,A011,A0A1,AUX0,AUX1,AUX2,AUX3... Temporary storage
c             locations.
C
C.......................................................................
C
C     Travel times at points Y0,Y1:
      ABSTT0=ABS(TT0)
      ABSTT1=ABS(TT1)
C
C     Slownesses at points Y0,Y1,Y2:
      U0=SQRT(P10*P10+P20*P20)
      U1=SQRT(P11*P11+P21*P21)
      U2=SQRT(P12*P12+P22*P22)
C
C     Slowness gradient:
      UG1=(U1-U0)/DX1
      UG2=(U2-U0)/DX2
C
C     Diffraction at point Y0:
      AUX0=ABS(DX2)
      TT=ABSTT0+0.5*(U0+U2)*AUX0
      IF(TT.LT.TT2) THEN
        TT2=TT
C       AUX1= 0.
C       AUX2= DX2
C       AUX1=0.5*((U0+U2)*AUX1/AUX0+UG1*AUX0)
C       AUX2=0.5*((U0+U2)*AUX2/AUX0+UG2*AUX0)
C       AUX1=0.5*(                  UG1*AUX0)
C       AUX2=0.5*((U0+U2)* DX2/AUX0+UG2*AUX0)
        AUX0=AUX0*AUX0/(U0+U2)
        AUX1=     UG1*AUX0
        AUX2= DX2+UG2*AUX0
        AUX0=SQRT(AUX1*AUX1+AUX2*AUX2)
        P12=U2*AUX1/AUX0
        P22=U2*AUX2/AUX0
      END IF
C
C     Diffraction at point Y1:
      AUX0=SQRT(DX1*DX1+DX2*DX2)
      TT=ABSTT1+0.5*(U1+U2)*AUX0
      IF(TT.LT.TT2) THEN
        TT2=TT
        AUX1=-DX1
        AUX2= DX2
C       AUX1=0.5*((U1+U2)*AUX1/AUX0+UG1*AUX0)
C       AUX2=0.5*((U1+U2)*AUX2/AUX0+UG2*AUX0)
        AUX0=AUX0*AUX0/(U1+U2)
        AUX1=AUX1+UG1*AUX0
        AUX2=AUX2+UG2*AUX0
        AUX0=SQRT(AUX1*AUX1+AUX2*AUX2)
        P12=U2*AUX1/AUX0
        P22=U2*AUX2/AUX0
      END IF
C
C     Check for the vicinity of the point source:
      AUX1=0.20*(U0+U1)*ABS(DX1)
      IF(ABSTT0.LT.AUX1.OR.ABSTT1.LT.AUX1) THEN
C       Not updating inside very close vicinity of the point source
        RETURN
      END IF
C
C     Check for the travel-time ridge:
      AUX1=U1-U0
      AUX2=U2-U0
      IF(PR1-PR0.LT.
     *   RIDGE1*AMIN1(PR2-PR1,0.)-RIDGE2*SQRT(AUX1*AUX1+AUX2*AUX2)) THEN
        IF(P20/DX2.GT.0..AND.-P20/DX2.LE.P10/DX1.AND.
     *                                   P10/DX1.LE.0.) THEN
          T12= U*UG1         /P20
          T22=(U*UG2-P10*T12)/P20
          P1 =P10        +T12*DX2
          P2 =P20        +T22*DX2
          AUX0=SQRT(P1*P1+P2*P2)
          P1 =P1*U2/AUX0
          P2 =P2*U2/AUX0
          TT =ABSTT0+0.5* (P2+P20)*DX2
          IF(TT.LE.TT2) THEN
            TT2=TT
            P12=P1
            P22=P2
          END IF
        END IF
        IF(P21/DX2.GT.0..AND.-P21/DX2.LE.P11/DX1.AND.
     *                                   P11/DX1.LE.0.) THEN
          T12= U*UG1         /P21
          T22=(U*UG2-P11*T12)/P21
          P1 =P11        +T12*DX2
          P2 =P21-T12*DX1+T22*DX2
          AUX0=SQRT(P1*P1+P2*P2)
          P1 =P1*U2/AUX0
          P2 =P2*U2/AUX0
          TT =ABSTT1+0.5*((P2+P21)*DX2-(P1+P11)*DX1)
          IF(TT.LE.TT2) THEN
            TT2=TT
            P12=P1
            P22=P2
          END IF
        END IF
        RETURN
      END IF
C
C     Check for the direction of propagation (from Y0 and Y1 to Y2):
      IF(P20/DX2.LE..000001*U0/ABS(DX2).OR.
     *   P21/DX2.LE..000001*U1/ABS(DX2)) THEN
        RETURN
      END IF
C
C     Projecting points y0,y1 onto the parallel grid line through y2:
      P10P20=P10/P20
      P11P21=P11/P21
      A0=    DX2*P10P20
      A1=DX1+DX2*P11P21
C     Now:  A0,A1 are X1-coordinates of the projections of points Y0,Y1.
C
C     Check for the position of point Y2 between projections of Y0,Y1:
      IF(A0/DX1.GE.0..OR.A1/DX1.LE.0.) THEN
        RETURN
      END IF
C
C     Backprojecting point Y2 onto point Y at gridline through Y0,Y1:
      AUX0=A0-A1
      A0=A0/AUX0
      A1=A1/AUX0
      X1=A0*DX1
C     Now:  A0=(X1-X10)/(X11-X10),  A1=(X1-X11)/(X11-X10),  A0-A1=1,
C     where X1,X10,X11 are X1-coordinates of points Y,Y0,Y1.
C
C     Slowness at point Y:
C     U =A0*U1-A1*U0
      U =U0+UG1*X1
C
      IF(VER1*ABS(P10P20-P11P21).GE.1.) THEN
C
        X2C=DX1/(P10P20-P11P21)
        X1C=X2C*P10P20
C       Now:  X1C,X2C is the crosssection of two given slowness vectors.
C
C       Travel-time spherical correction
        UC=U0+UG1*X1C
        X2CX2C=X2C*X2C
        UG1X2C=UG1*X2CX2C
        AUX1=   -X1C
        AUX0=SQRT(AUX1*AUX1+X2CX2C)
        ABSTT0=ABSTT0-0.5*(AUX0*(U0+UC)-UG1X2C*ALOG(AUX0+AUX1))
        AUX1=DX1-X1C
        AUX0=SQRT(AUX1*AUX1+X2CX2C)
        ABSTT1=ABSTT1-0.5*(AUX0*(U1+UC)-UG1X2C*ALOG(AUX0+AUX1))
        AUX1= X1-X1C
        AUX0=SQRT(AUX1*AUX1+X2CX2C)
        TT    =       0.5*(AUX0*(U +UC)-UG1X2C*ALOG(AUX0+AUX1))
        T1    = U  *AUX1/AUX0
        T11   =(UG1*AUX1+U*X2CX2C/AUX0/AUX0)/AUX0
C
C       Travel-time cubic polynomial correction
        AUX0=ABSTT0-ABSTT1
        AUX1=6.*AUX0/DX1
        AUX2=1.-VER2
        TT =TT +AUX2*((A0+A0+1.)*A1*A1*AUX0 )+ABSTT1
        T1 =T1 +AUX2*(           A0*A1*AUX1 )
        T11=T11+AUX2*(          (A0+A1)*AUX1)/DX1
        TT =TT -VER2*A1*AUX0
        T1 =T1 -VER2*AUX1/6.
C
      ELSE
C
C       Interpolating travel-time and its derivatives at point Y by
C       means of the third-order interpolation from Y0,Y1:
        A01=A0+A1
        A001=A0+A01
        A011=A01+A1
        A0A1=A0*A1
        AUX0=ABSTT0-ABSTT1
        AUX1=6.*AUX0/DX1
        AUX3=VER2*((P10+P11)/2.+AUX1/6.)
        R10=     P10-AUX3
        R11=     P11-AUX3
        TT =(A0+A0+1.)*A1*A1*AUX0+ABSTT1+A0A1*(A1*R10+     A0*R11)*DX1
        T1 =           A0A1 *AUX1       + A001*A1*R10+A011*A0*R11
        T11=(          A01  *AUX1       +2.*(A011*R10+   A001*R11))/DX1
C
      END IF
C
C     Second-order Taylor expansion of travel time at point Y:
      T2=U*U-T1*T1
      IF(T2.LT.0.) THEN
        WRITE(*,'(1X,11F7.4)')
     *                       TT0,P10,P20,DX1,TT1,P11,P21,DX2,TT2,P12,P22
        WRITE(*,'(A)') ' TOO LARGE SLOWNESS COMPONENT'
        RETURN
      END IF
      T2=SIGN(SQRT(T2),DX2)
      T12=(U*UG1-T1*T11)/T2
      T22=(U*UG2-T1*T12)/T2
      P1 =T1-T11*X1+T12*DX2
      P2 =T2-T12*X1+T22*DX2
      AUX0=SQRT(P1*P1+P2*P2)
      P1=P1*U2/AUX0
      P2=P2*U2/AUX0
      TT =TT+0.5*((P2+T2)*DX2-(P1+T1)*X1)
      IF(TT.LE.TT2) THEN
        TT2=TT
        P12=P1
        P22=P2
      END IF
C
      RETURN
      END
C
C=======================================================================
C