C 
C Program WATABLE to read the results of common anisotropic ray tracing,
C isotropic ray tracing, and anisotropic ray tracing, and to compare
C the results.
C
C Version: 5.90
C Date: 2005, June 12
C
C Coded by Petr Bulant
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     E-mail: bulant@seis.karlov.mff.cuni.cz
C
C.......................................................................
C
C                                                    
C Description of data files:
C
C Main input data read from external interactive device (*):
C     The data consist of character strings read by list directed (free
C     format) input.  The strings have thus to be enclosed in
C     apostrophes.  The interactive * external unit may be redirected to
C     the file containing the data.
C (1) 'SEP'/
C     'SEP'...Name of the file with input parameters.
C             Description of file SEP
C             If blank, default values of the corresponding data are
C             considered.
C     No default, 'SEP' must be specified and cannot be blank.
C
C                                                     
C Input file 'SEP' in the SEP format:
C     The file has the form of a SEP parameter file.
C     For the description of the SEP format refer to file
C     'sep.for'.
C Names of optional input/output files related to ray tracing using
C program 'crt.for':
C     QILSTI='string' ... Name of the optional input/output file
C             containing second-order perturbations calculated along
C             isotropic reference rays.
C             Default: QILSTI=' ' (the file is not read/written)
C     QILSTA='string' ... Name of the optional input/output file
C             containing second-order perturbations calculated along
C             common anisotropic reference rays.
C             Default: QILSTA=' ' (the file is not read/written)
C             Description of QILSTI and QILSTA.
C Names of optional input files related to ray tracing using
C program 'anray.for':
C     ANRTIM1='string' ... Name of the optional input file containing
C             travel times of the faster anisotropic S wave.
C             Default: ANRTIM1=' ' (the file is not read)
C     ANRTIM2='string' ... Name of the optional input file containing
C             travel times of the slower anisotropic S wave.
C             Default: ANRTIM2=' ' (the file is not read)
C             Description of
C             ANRTIM1 and ANRTIM2.
C
C                                                   
C Input/output formatted file QILSTI:
C (1) / (a slash).
C (2) For each two-point ray one line (2.1):
C     On input, the lines have the following form:
C (2.1) IREC,TI,DTLI1,DTLI2,DTQI1,DTQI2,T1,T1,T2,T2,
C       DTLIA,DTQIA,TA,TA,DTQA1,DTQA2,DTQ12 /
C       For the description of the quantities refer to file
C       'green.for'.
C     On output, the lines have the following form:
C (2.1) IREC,TI,DTLI1,DTLI2,DTQI1,DTQI2,T1,DT1,T2,DT2,
C       DTLIA,DTQIA,TA,DTA,DTQA1,DTQA2,DTQ12 /
C     DT1,DT2 ... If the input files ANRTIM1 and ANRTIM2 with the
C             anisotropic travel times are given, the difference
C             between the anisotropic travel times and their
C             approximation by the second-order perturbation expansion
C             calculated along the isotropic S-wave reference ray.
C             Otherwise the copy of T1 and T2.
C     DTA ... If the input file QILSTA with the perturbations along
C             common anisotropic rays is given, the difference
C             between the common anisotropic travel time and its
C             approximation by the second-order perturbation expansion
C             calculated along the isotropic S-wave reference ray.
C             Otherwise the copy of TA.
C (3) / (a slash).
C
C Input/output formatted file QILSTA:
C (1) / (a slash).
C (2) For each two-point ray one line (2.1):
C     On input, the lines have the following form:
C (2.1) IREC,TA,DTLA1,DTLA2,DTQA1,DTQA2,T1,T1,T2,T2,
C       DTLAI,DTQAI,TI,TI,DTQI1,DTQI2,DTQ12 /
C       For the description of the quantities refer to file
C       'green.for'.
C     On output, the lines have the following form:
C (2.1) IREC,TA,DTLA1,DTLA2,DTQA1,DTQA2,T1,DT1,T2,DT2,
C       DTLAI,DTQAI,TI,DTI,DTQI1,DTQI2,DTQ12 /
C     DT1,DT2 ... If the input files ANRTIM1 and ANRTIM2 with the
C             anisotropic travel times are given, the difference
C             between the anisotropic travel times and their
C             approximation by the second-order perturbation expansion
C             calculated along the common anisotropic reference ray.
C             Otherwise the copy of T1 and T2.
C     DTI ... If the input file QILSTI with the perturbations along
C             isotropic rays is given, the difference
C             between the isotropic travel time and its
C             approximation by the second-order perturbation expansion
C             calculated along the common anisotropic reference ray.
C             Otherwise the copy of TI.
C (3) / (a slash).
C
C-----------------------------------------------------------------------
C
      INTEGER MTPR
      PARAMETER (MTPR=30)
      INTEGER IREC(MTPR),I,J
      REAL TT1(MTPR),TT2(MTPR),TTI(MTPR),DLI1(MTPR),DLI2(MTPR),
     *     DQI1(MTPR),DQI2(MTPR),TTI1(MTPR),TTI2(MTPR),DLIA(MTPR),
     *     DQIA(MTPR),TTIA(MTPR),TTA(MTPR),DLA1(MTPR),
     *     DLA2(MTPR),DQA1(MTPR),DQA2(MTPR),TTA1(MTPR),TTA2(MTPR),
     *     DTIA(MTPR),DTI1(MTPR),DTI2(MTPR),
     *     DTAI(MTPR),DTA1(MTPR),DTA2(MTPR),
     *     DLAI(MTPR),DQAI(MTPR),TTAI(MTPR),DQA1I(MTPR),DQA2I(MTPR),
     *     DQ12I(MTPR),DQI1A(MTPR),DQI2A(MTPR),DQ12A(MTPR)
      REAL QICOR(16)
      REAL T,TI,DTLI1,DTLI2,DTQI1,DTQI2,T1,TI1,T2,TI2,DTLIA,DTQIA,
     *     TA,TIA,DTQA1,DTQA2,DTQ12,DTLA1,DTLA2,TA1,TA2,
     *     DTLAI,DTQAI,TAI
      CHARACTER*135 FORMQI
      CHARACTER*80 FQI,FQA,FA1,FA2,FILSEP,TEXT
      DATA IREC/MTPR*0/
      DATA TT1/MTPR*0./,TT2/MTPR*0./,TTI/MTPR*0./,DLI1/MTPR*0./,
     *     DLI2/MTPR*0./,
     *     DQI1/MTPR*0./,DQI2/MTPR*0./,TTI1/MTPR*0./,TTI2/MTPR*0./,
     *     DLIA/MTPR*0./,
     *     DQIA/MTPR*0./,TTIA/MTPR*0./,TTA/MTPR*0./,
     *     DLA1/MTPR*0./,
     *     DLA2/MTPR*0./,DQA1/MTPR*0./,DQA2/MTPR*0./,TTA1/MTPR*0./,
     *     TTA2/MTPR*0./,
     *     DTIA/MTPR*0./,DTI1/MTPR*0./,DTI2/MTPR*0./,
     *     DTAI/MTPR*0./,DTA1/MTPR*0./,DTA2/MTPR*0./,
     *     DLAI/MTPR*0./,DQAI/MTPR*0./,TTAI/MTPR*0./,
     *     DQA1I/MTPR*0./,DQA2I/MTPR*0./,DQ12I/MTPR*0./,
     *     DQI1A/MTPR*0./,DQI2A/MTPR*0./,DQ12A/MTPR*0./
C
      WRITE(*,'(A)') '+WATABLE: Enter input filename: '
      FILSEP=' '
      READ(*,*) FILSEP
C
C     Reading all data from the SEP file into the memory:
      IF (FILSEP.NE.' ') THEN
        CALL RSEP1(5,FILSEP)
      ELSE
C       WATABLE-01
        CALL ERROR('WATABLE-01: SEP file not given')
C       Input file in the form of the SEP (Stanford Exploration Project)
C       parameter or history file must be specified.
C       There is no default filename.
      ENDIF
C
      WRITE(*,'(A)') '+WATABLE: Working ...           '
      CALL RSEP3T('QILSTI',  FQI,' ')
      CALL RSEP3T('QILSTA',  FQA,' ')
      CALL RSEP3T('ANRTIM1', FA1,' ')
      CALL RSEP3T('ANRTIM2', FA2,' ')
      IF (FQI.NE.' ') OPEN(1,FILE=FQI)
      IF (FQA.NE.' ') OPEN(2,FILE=FQA)
      IF (FA1.NE.' ') OPEN(3,FILE=FA1)
      IF (FA2.NE.' ') OPEN(4,FILE=FA2)
      IF (FQI.NE.' ') READ(1,*) TEXT
      IF (FQA.NE.' ') READ(2,*) TEXT
      FORMQI(1:7)='(1I6,A,'
C
      IF (FA1.NE.' ') THEN
   10   CONTINUE
          I=0
          READ(3,*,END=11) I,T
          IF(I.EQ.0) GOTO 11
          IREC(I)=I
          TT1 (I)=T
        GOTO 10
   11   CONTINUE
      ENDIF
      IF (FA2.NE.' ') THEN
   20   CONTINUE
          I=0
          READ(4,*,END=21) I,T
          IF(I.EQ.0) GOTO 21
          IREC(I)=I
          TT2 (I)=T
        GOTO 20
   21   CONTINUE
      ENDIF
      IF (FQI.NE.' ') THEN
   30   CONTINUE
          I=0
          READ(1,*,END=31) I,TI,DTLI1,DTLI2,DTQI1,DTQI2,T1,TI1,T2,TI2,
     *                     DTLIA,DTQIA,TA,TIA,DTQA1,DTQA2,DTQ12
          IF(I.EQ.0) GOTO 31
          IREC (I)=I
          TTI  (I)=TI
          DLI1 (I)=DTLI1
          DLI2 (I)=DTLI2
          DQI1 (I)=DTQI1
          DQI2 (I)=DTQI2
          TTI1 (I)=T1
          DTI1 (I)=TI1
          TTI2 (I)=T2
          DTI2 (I)=TI2
          DLIA (I)=DTLIA
          DQIA (I)=DTQIA
          TTIA (I)=TA
          DTIA (I)=TIA
          DQA1I(I)=DTQA1
          DQA2I(I)=DTQA2
          DQ12I(I)=DTQ12
        GOTO 30
   31   CONTINUE
      ENDIF
      IF (FQA.NE.' ') THEN
   40   CONTINUE
          I=0
          READ(2,*,END=41) I,TA,DTLA1,DTLA2,DTQA1,DTQA2,T1,TA1,T2,TA2,
     *                     DTLAI,DTQAI,TI,TAI,DTQI1,DTQI2,DTQ12
          IF(I.EQ.0) GOTO 41
          IREC (I)=I
          TTA  (I)=TA
          DLA1 (I)=DTLA1
          DLA2 (I)=DTLA2
          DQA1 (I)=DTQA1
          DQA2 (I)=DTQA2
          TTA1 (I)=T1
          DTA1 (I)=TA1
          TTA2 (I)=T2
          DTA2 (I)=TA2
          DLAI (I)=DTLAI
          DQAI (I)=DTQAI
          TTAI (I)=TI
          DTAI (I)=TAI
          DQI1A(I)=DTQI1
          DQI2A(I)=DTQI2
          DQ12A(I)=DTQ12
        GOTO 40
   41   CONTINUE
      ENDIF
C
      IF (FQI.NE.' ') THEN
        REWIND(1)
        WRITE(1,'(A)') ' /'
        DO 50, I=1,MTPR
          IF (IREC(I).NE.0) THEN
            IF (TTA(I).NE.0..AND.TTIA(I).NE.0.) DTIA(I)=TTA(I)-TTIA(I)
            IF (TT1(I).NE.0..AND.TTI1(I).NE.0.) DTI1(I)=TT1(I)-TTI1(I)
            IF (TT2(I).NE.0..AND.TTI2(I).NE.0.) DTI2(I)=TT2(I)-TTI2(I)
            QICOR( 1)=TTI  (I)
            QICOR( 2)=DLI1 (I)
            QICOR( 3)=DLI2 (I)
            QICOR( 4)=DQI1 (I)
            QICOR( 5)=DQI2 (I)
            QICOR( 6)=TTI1 (I)
            QICOR( 7)=DTI1 (I)
            QICOR( 8)=TTI2 (I)
            QICOR( 9)=DTI2 (I)
            QICOR(10)=DLIA (I)
            QICOR(11)=DQIA (I)
            QICOR(12)=TTIA (I)
            QICOR(13)=DTIA (I)
            QICOR(14)=DQA1I(I)
            QICOR(15)=DQA2I(I)
            QICOR(16)=DQ12I(I)
            CALL FORM2(16,QICOR,QICOR,FORMQI(8:135))
            WRITE(1,FORMQI) IREC(I),(' ',QICOR(J),J=1,16),' /'
          ENDIF
   50   CONTINUE
        WRITE(1,'(A)') ' /'
      ENDIF
      IF (FQA.NE.' ') THEN
        REWIND(2)
        WRITE(2,'(A)') ' /'
        DO 60, I=1,MTPR
          IF (IREC(I).NE.0) THEN
            IF (TTI(I).NE.0..AND.TTAI(I).NE.0.) DTAI(I)=TTI(I)-TTAI(I)
            IF (TT1(I).NE.0..AND.TTA1(I).NE.0.) DTA1(I)=TT1(I)-TTA1(I)
            IF (TT2(I).NE.0..AND.TTA2(I).NE.0.) DTA2(I)=TT2(I)-TTA2(I)
            QICOR( 1)=TTA  (I)
            QICOR( 2)=DLA1 (I)
            QICOR( 3)=DLA2 (I)
            QICOR( 4)=DQA1 (I)
            QICOR( 5)=DQA2 (I)
            QICOR( 6)=TTA1 (I)
            QICOR( 7)=DTA1 (I)
            QICOR( 8)=TTA2 (I)
            QICOR( 9)=DTA2 (I)
            QICOR(10)=DLAI (I)
            QICOR(11)=DQAI (I)
            QICOR(12)=TTAI (I)
            QICOR(13)=DTAI (I)
            QICOR(14)=DQI1A(I)
            QICOR(15)=DQI2A(I)
            QICOR(16)=DQ12A(I)
            CALL FORM2(16,QICOR,QICOR,FORMQI(8:135))
            WRITE(2,FORMQI) IREC(I),(' ',QICOR(J),J=1,16),' /'
          ENDIF
   60   CONTINUE
        WRITE(2,'(A)') ' /'
      ENDIF
C
      IF (FQI.NE.' ') CLOSE(1)
      IF (FQA.NE.' ') CLOSE(2)
      IF (FA1.NE.' ') CLOSE(3)
      IF (FA2.NE.' ') CLOSE(4)
      WRITE(*,'(A)') '+WATABLE: Done.                 '
      STOP
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'forms.for'
C     forms.for
      INCLUDE 'length.for'
C     length.for
C
C=======================================================================
C