C 
C Program INVTT to calculate the derivatives of travel times
C with respect to the model parameters (B-spline coefficients)
C
C Version: 6.60
C Date: 2011, November 29
C
C Coded by: Ludek Klimes
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     http://sw3d.cz/staff/klimes.htm
C
C.......................................................................
C
C Program INVTT assumes all model parameters (B-spline coefficients)
C stored in the common block /VALC/ as in the submitted versions of
C user-defined model specification FORTRAN77 source code files
C 'srfc.for', 'parm.for' and 'val.for'.  Thus, unlike the other parts of
C the complete ray tracing, the INVTT program cannot work with user's
C modifications of subroutines SRFC1, SRFC2, PARM1, and PARM2.
C
C In the input data file
C DCRT of ray tracing
C program 'crt.for', step STORE of independent variable for storing the
C computed quantities along rays should be sufficiently small to render
C the numerical quadrature along rays accurate.  The relative error of
C the numerical quadrature is proportional to the fourth power of the
C step STORE along the rays.  When integrating a B-spline in a regular
C grid, the relative error is about 0.01 for the step of the size of the
C grid interval.  For example, in a regular grid, the relative error is
C about 0.01*0.1**4=0.000001 for the step of 0.1 the size of the grid
C interval.  Note that the numerical quadrature is performed by
C subroutine AP28, called
C by subroutine AP29.
C
C.......................................................................
C                                                    
C Description of data files:
C
C Input data read from the standard input device (*):
C     The data are read by the list directed input (free format) and
C     consist of a single string 'SEP':
C     'SEP'...String in apostrophes containing the name of the input
C             SEP parameter or history file with the input data.
C     No default, 'SEP' must be specified and cannot be blank.
C
C                                                     
C Input data file 'SEP':
C     File 'SEP' has the form of the SEP
C     parameter file.  The parameters, which do not differ from their
C     defaults, need not be specified in file 'SEP'.
C Data specifying the model:
C     MODEL='string'... Name of the input data file describing the
C             model.  For the description of the data refer to
C             file 
C             model.for.
C             Default: MODEL='model.dat'
C     NEGPAR=integer... Flag whether the negative values of material
C             parameters are allowed.  The default value is recommended
C             for this program.
C             NEGPAR=0: Negative values of material parameters or zero
C               P-wave velocity are reported as errors.
C             NEGPAR=1: Negative values of material parameters or zero
C               P-wave velocity are not reported as errors.
C             Default: NEGPAR=0
C Data specifying other input files:
C     CRTOUT='string'...Name of the file with the names of the output
C             files of program CRT. If blank, default names are
C             considered.
C             Only the files 'CRT-R', 'CRT-S' and 'CRT-I' are read
C             by INVTT.
C             For general description of file CRTOUT refer to file
C             writ.for.
C             Default: CRTOUT=' ' (usually sufficient for the first
C               elementary wave)
C     PTS='string'... String with the name of the input data file
C             containing the coordinates of shot and receiver points.
C             Description of file PTS
C             No default, PTS must be specified and cannot be blank.
C     PTS1='string',...,PTS9='string',... Optional strings with the
C             names of the input data files containing the coordinates
C             of additional shot or receiver points, supplementing file
C             PTS.
C             Default: PTS1=' ',...,PTS9=' ' (no additional points).
C     FTT='string'... Input data file containing the travel times from
C             the field measurements.  Its structure is described below.
C             This program may be executed several times to take into
C             account several files with travel times.
C             Description of file FTT
C             No default, FTT must be specified and cannot be blank.
C     M1IN='string'... Name of the input file containing number M1IN of
C             model and source parameters to be inverted.
C             The corresponding M1IN*M2IN input values of file GM1
C             are supplemented by zeros to get M1*M2IN input values.
C             Default: M1IN=' ' means M1IN=NM (number of model
C               parameters).
C     M2IN='string'... Name of the input file containing number M2IN of
C             values already processed by programs 'invpts.for' or
C             'invtt.for', i.e., the name of output file M2 from the
C             previous execution of 'invpts.for' or 'invtt.for'.
C             The corresponding output values of this program will be
C             appended after M1IN*M2IN input values of file GM1 and
C             after M2IN input values of files GM2, GM3 and DM1,
C             respectively.
C             Default: M2IN=' ' means M2IN=0 (no previous values).
C Data specifying output files:
C     M1='string'... Name of the output file containing number M1=NM of
C             model parameters (a single integer).
C             If simultaneous source location is enabled (INVSRC is not
C             blank), number M1=NM+4*NSRC is the number of model
C             parameters plus the number of parameters (coordinates and
C             time) of NSRC sources to be located.  If INVSRC=' ',
C             an equal file may be generated by program 'invsoft.for'.
C             The file is not generated if the value of M1 is blank.
C             Note that the model parameters describing structural
C             interfaces only, or describing material parameters only,
C             may be selected by means of input parameter ICLASS, see
C             below.
C             Default: M1=' '
C             Note: Default of 'invsoft.for' is M1='m1.out',
C                   default of 'invpts.for' is M1=' '.
C     M2='string'... Name of the output file containing number M2 of
C             accumulated values to be inverted (a single integer).
C             M2 is M2IN increased by the number of given travel times
C             matched by two-point rays.
C             The file is not generated if the value of M2 is blank.
C             Default: M2='m2.out'
C     INVSRC='string'... Name of the optional output file containing the
C             names of sources which coordinates are to be simultaneosly
C             inverted.
C             INVSRC=' ' (default): No simultaneous source location.
C             INVSRC.NE.' ': The output file contains the list of the
C               names of sources, which space-time hypocentral
C               coordinates are to be simultaneously inverted.
C               In this case, the travel times in file FTT are deemed
C               to be arrival times at the receivers.
C     INVLOG='string'... Output log file.  Just for your information.
C             Not opened and generated if INVLOG=' '.
C             Description of file INVLOG
C             Default: INVLOG='invlog.out'
C Data specifying names of header files of the input/output files with
C matrix and vector components:
C     GM1='string'... String with the name of the input/output file
C             to accumulate the matrix of the derivatives of M2 travel
C             times with respect to M1 model parameters.
C             The matrix has M1*M2IN components on input and M1*M2
C             components on output.  The formatted file is composed of
C             real-valued matrix components to be read at once by the
C             list-directed input.
C             If the filename is blank, no file is read nor written.
C             The file is not read, just written if M2IN=0.
C             Default: GM1=' '
C     GM2='string'... String with the name of the input/output file
C             containing the vector composed of differences between the
C             given travel times and the travel times calculated in
C             the model.  In the case of multiple two-point rays, the
C             ray of the smallest travel time is considered.
C             The vector has M2IN components on input and M2 components
C             on output.  The formatted file is composed of real-valued
C             vector components to be read at once by list-directed
C             input.
C             If the filename is blank, no file is read nor written.
C             The file is not read, just written if M2IN=0.
C             Default: GM2=' '
C     GM3='string'... String with the name of the input/output file
C             containing the vector composed of travel times calculated
C             in the model.  The vector has M2IN components on input and
C             M2 components on output.  The formatted file is composed
C             of real-valued vector components to be read at once by
C             list-directed input.
C             If the filename is blank, no file is read nor written.
C             The file is not read, just written if M2IN=0.
C             Default: GM3=' '
C     DM1='string'... String with the name of the input/output file
C             containing the diagonal matrix composed of the variances
C             (squares of standard deviations) of the travel times.
C             The diagonal has M2IN components on input and M2
C             components on output.  The formatted file is composed
C             of real-valued diagonal components to be read at once by
C             list-directed input.
C             If the filename is blank, no file is read nor written.
C             The file is not read, just written if M2IN=0.
C             Default: DM1=' '
C     Recent version of the program cannot deal with sparse or
C     unformatted matrices.
C     For general description of the files with matrices refer to file
C     forms.htm.
C Other numerical parameters:
C     ICLASS=integer... Class of model parameters to be inverted:
C             ICLASS=0: All model parameters are inverted.
C             ICLASS=1: Only model parameters describing interfaces are
C               inverted.
C             ICLASS=2: Only model parameters describing material
C               parameters are inverted.
C             Default: ICLASS=0
C     DIST=real... Maximum distance between the source or receiver point
C             and the initial or end point of a synthetic ray.
C             No default, DIST must be specified and must be positive.
C     VPOWER=real... Velocity power for the computation of the
C             travel-time check sum.  If the VPOWER-th velocity power is
C             expressed, in all blocks of the model, in terms of
C             explicit functions of model coordinates, linearly
C             homogeneous in their parameters (e.g., in B-spline
C             coefficients), the travel time minus its check sum (see
C             the log output file INVLOG) should be zero within rounding
C             errors.  Otherwise, the check sum may have no sense.
C             VPOWER=0: no check sum is evaluated and printed.
C             Default: VPOWER=0.
C
C                                                     
C Input data file PTS:
C     This data file contains the coordinates of shot and receiver
C     points.  The data are read in by the list directed input
C     (free format).  In the list of input data below, each numbered
C     paragraph indicates the beginning of a new input operation (new
C     READ statement).  The CHARACTER strings are explicitly mentioned
C     in this description.  Otherwise, if the first letter of the
C     symbolic name of the input variable is I-N, the corresponding
C     value in input data must be of the type INTEGER.  Otherwise, the
C     input parameter is of the type REAL.
C (1) Several strings terminated by / (a slash).
C (2) List of the sources and receivers: Any times the following data:
C (2.1) POINT,X1,X2,X3,X4
C     POINT...CHARACTER*11 string containing the name of the source or
C             receiver point.
C     X1,X2,X3... Coordinates of the source or receiver point.
C     X4...   Hypocentral time.  Considered only in the case of
C             simultaneous location (INVSRC.NE.' ') and only if POINT is
C             a source listed in file FTT.
C (3) / (a slash) or the end of file.
C Example of data
C PTS.
C
C                                                     
C Input data file FTT:
C     This data file contains the travel times from the field
C     measurements.  The data are read in by the list directed input
C     (free format).  In the list of input data below, each numbered
C     paragraph indicates the beginning of a new input operation (new
C     READ statement).  The CHARACTER strings are explicitly mentioned
C     in this description.  Otherwise, if the first letter of the
C     symbolic name of the input variable is I-N, the corresponding
C     value in input data must be of the type INTEGER.  Otherwise, the
C     input parameter is of the type REAL.
C (1) Several strings terminated by / (a slash).
C (2) List of the travel times: Any times the following data (2.1):
C (2.1) 'SOURCE','RECEIVER',TFIELD,TERR
C     'SOURCE'... String up to 11 characters enclosed in apostrophes,
C             containing the name of the source point.
C     'RECEIVER'... String up to 11 characters enclosed in apostrophes,
C             containing the name of the receiver point.
C             INVSRC=' ' (no simultaneous location), the source and
C             receiver points may be mutually interchanged.
C     TFIELD..Travel time from a field measurement if INVSRC=' '.
C             Arrival time at the receiver otherwise.
C     TERR... Error of the travel time from a field measurement.
C (3) / (a slash) or the end of file.
C Example of data
C FTT.
C
C                                                  
C Output log file INVLOG:
C (1) For each considered ray:
C (1.1) SOURCE,RECEIVER,TFIELD,TDIF,SDIST,RDIST,CHECK
C     SOURCE..Name of the source point.
C     RECEIVER... Name of the receiver point.
C     TFIELD..Travel time from a field measurement.
C     TDIF... Field travel time minus the minimum synthetic travel time.
C     SDIST...Distance between the source and the initial point of the
C             synthetic ray.
C     RDIST...Distance between the receiver and the end point of the
C             synthetic ray.
C     CHECK...Synthetic travel time minus the travel time resulting from
C             the derivatives of the theoretical travel time with
C             respect to the model parameters.  This quantity should
C             not exceed in order the numerical error of the synthetic
C             travel time.
C             In this version defined just for the models described in
C             the terms of velocity.
C
C-----------------------------------------------------------------------
C
C Common block /VALC/:
      INCLUDE 'val.inc'
C     val.inc
C None of the storage locations of the common block are altered.
C
C Common block /POINTC/:
      INCLUDE 'pointc.inc'
C     pointc.inc
C
C-----------------------------------------------------------------------
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
      INTEGER IRAM(MRAM)
      EQUIVALENCE (IRAM,RAM)
C
C     Allocation of working arrays:
C
      INTEGER MP,MT
      PARAMETER (MP=1000)
      PARAMETER (MT=4000)
      INTEGER KS(MT),KR(MT),INDM(NPAR)
      REAL COOR(4,MP),TFIELD(MT),TERR(MT),SUM(NPAR)
      EQUIVALENCE (KS    ,IRAM(    1))
      EQUIVALENCE (KR    ,IRAM( MT+1))
      EQUIVALENCE (TFIELD,RAM(2*MT+1))
      EQUIVALENCE (TERR  ,RAM(3*MT+1))
      EQUIVALENCE (COOR  ,RAM(4*MT+1))
      EQUIVALENCE (SUM   ,RAM(4*MT+4*MP+1))
      EQUIVALENCE (INDM  ,IRAM(4*MT+4*MP+NPAR+1))
      CHARACTER*11 POINT(MP)
      COMMON/PTSC/ POINT
C     INDM... Indices of the unknown model parameters.
C
C     Arrays to call subroutine SOFT
      REAL W(47,2),CS(1)
C
C.......................................................................
C
C     Filenames:
      CHARACTER*80 FILE1,FILE2,FILE3,FILE4,FILED
C
C     Logical unit numbers:
      INTEGER LU1,LU2,LU3,LU4,LU5,LU6,IU1,IU2
      PARAMETER (LU1=11)
      PARAMETER (LU2=12)
      PARAMETER (LU3=13)
      PARAMETER (LU4=14)
      PARAMETER (LU6=16)
C
C     Input data:
      CHARACTER*4  FILPTS
      CHARACTER*1  TEXT(10)
      CHARACTER*80 ERRTXT
      CHARACTER*11 SRC,REC
      INTEGER NP,NT,ICLASS
      REAL DIST,DIST2,VPOWER
C     POINT...Names of the source and receiver points.
C     NP...   Number of source and receiver points.
C     NT...   Number of field travel times.
C     KS(I)...Index of the source point corresponding to the I-th field
C             travel time.
C     KR(I)...Index of the receiver point corresponding to the I-th
C             field travel time.
C     DIST... Maximum distance between the source or receiver point and
C             the initial or end point of a synthetic ray.
C     DIST2...DIST**2
C     VPOWER... Velocity power for the computation of the travel-time
C             check SUM.
C     COOR... Coordinates of the source or receiver points.
C     TFIELD..Field travel times.
C     TERR... Field travel time errors.
C
C     Output data:  variations of the synthetic travel time:
      INTEGER NSUM,NM,NSRCIN,NSRC
C     NM...   Number of the unknown model parameters.
C     NSRC... Number of sources to be simultaneously located.
C             NSRC=0 for NLOC=0, see below.
C
C     Auxiliary storage locations:
      INTEGER NLOC,IS,IR,IT,ND,IRAYTT,I,J,K,IPTS,IT0,M2IN,M2,LINLEN,M1
      REAL TT,TTAUX,TDIF,XI1,XI2,XI3,XE1,XE2,XE3,PI(6),PE(6)
      REAL SI,SE,RI,RE,SDIST,RDIST
C     NLOC... Number of derivatives of travel time with respect to the
C             source space-time coordinates for simultaneous location.
C             NLOC=0: No simultaneous location.
C             NLOC=4: Simultaneous location.
C     IS...   Index of the source point.
C     IR...   Index of the receiver point.
C     IT...   Index of the field travel time.
C     ND...   Number of synthetic travel times corresponding to the
C             field travel times.
C     IRAYTT..Index of the last processed ray.
C     I,J,K...Temporary storage locations.
C     TT...   Synthetic travel time.
C     TTAUX...Temporary storage location.
C     TDIF... Field travel time minus the minimum synthetic travel time.
C     XI1,XI2,XI3,XE1,XE2,XE3... Coordinates of the initial and end
C             points of the last processed ray.
C     PI,PE...Slowness vectors at the initial and end points of the
C             last processed ray.
C     SI,SE,RI,RE... Squares of the distances between the source or
C             receiver points and the initial or end points of the ray.
C     SDIST...Distance between the source and the initial point of the
C             synthetic ray.
C     RDIST...Distance between the receiver and the end point of the
C             synthetic ray.
C
      CHARACTER*13 FORMAT
C     FORMAT..String containing the output format.
C
C-----------------------------------------------------------------------
C
C     Setting output format:
      FORMAT='(5(E13.7,1X))'
C
      IF(4*MT+4*MP+2*NPAR.GT.MRAM) THEN
C       INVTT-01
        CALL ERROR('INVTT-01: Too small dimension MRAM of array RAM')
C       Dimension MRAM of array RAM in include file
C       ram.inc should
C       be increased.
      END IF
C
C.......................................................................
C
C     Opening data files and reading the input data:
C
C     Reading main input data:
      WRITE(*,'(A)') '+INVTT: Enter input filename: '
      FILE1=' '
      READ (*,*) FILE1
      IF(FILE1.EQ.' ') THEN
C       INVTT-02
        CALL ERROR('INVTT-02: No input file specified')
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.
      END IF
      WRITE(*,'(A)') '+INVTT: Working...            '
      CALL RSEP1(LU1,FILE1)
C
C     Reading input data for the model:
      CALL RSEP3T('MODEL',FILE1 ,'model.dat')
      OPEN(LU1,FILE=FILE1,STATUS='OLD')
      CALL MODEL1(LU1)
      CLOSE(LU1)
C
C     Number of unknown model parameters:
      CALL RSEP3I('ICLASS',ICLASS,0)
      IF(ICLASS.LT.0.OR.2.LT.ICLASS) THEN
C       INVTT-03
        CALL ERROR('INVTT-03: Incorrect class index ICLASS')
C       The value of ICLASS must be 0, 1 or 2.
C       Check the input data.
      END IF
      DO 1 I=1,47
        W(I,1)=0.
        W(I,2)=0.
    1 CONTINUE
      NM=0
      IF(ICLASS.LE.1) THEN
        CALL SOFT(1,0,0,0,0,0,0,47,W,NM,INDM,CS,1,CS)
      END IF
      IF(ICLASS.EQ.0.OR.ICLASS.EQ.2) THEN
        CALL SOFT(2,0,0,0,0,0,0,47,W,NM,INDM,CS,1,CS)
      END IF
      WRITE(*,'(A,I4,A)') '+INVTT:',NM,' model parameters'
C
C     Simultaneous source location:
      CALL RSEP3T('INVSRC',FILE1 ,' ')
      IF(FILE1.EQ.' ') THEN
C       No simultaneous source location
        NLOC=0
      ELSE
C       Simultaneous source location
        NLOC=4
      END IF
C
C     Reading numerical parameters:
      CALL RSEP3R('DIST'  ,DIST  ,0.)
      IF(DIST.LE.0.) THEN
C       INVTT-04
        CALL ERROR('INVTT-04: Input parameter DIST not specified')
C       Input parameter DIST must be specified in the input SEP
C       parameter or history file, and must be positive.
C       There is no default value.
      END IF
      CALL RSEP3R('VPOWER',VPOWER,0.)
      DIST2=DIST*DIST
C
C.......................................................................
C
C     Reading source and receiver points:
      NP=0
      FILPTS='PTS'
      DO 12 IPTS=0,9
        IF(IPTS.GT.0) THEN
          FILPTS(4:4)=CHAR(ICHAR('0')+IPTS)
        ENDIF
        CALL RSEP3T(FILPTS,FILE1 ,' ')
        IF(FILE1.EQ.' '.AND.IPTS.EQ.0) THEN
C         INVTT-05
          CALL ERROR ('INVTT-05: File PTS not specified')
        ENDIF
        IF(FILE1.NE.' ') THEN
          OPEN(LU1,FILE=FILE1,STATUS='OLD')
          READ(LU1,*,END=11) TEXT
   10     CONTINUE
            IF(NP+1.GE.MP) THEN
C             INVTT-06
              CALL ERROR
     *                 ('INVTT-06: Too many source and receiver points')
            END IF
            POINT(NP+1)=' '
            READ(LU1,*,END=11) POINT(NP+1),(COOR(I,NP+1),I=1,4)
            IF(POINT(NP+1).EQ.' ') THEN
              GO TO 11
            END IF
            NP=NP+1
          GO TO 10
   11     CONTINUE
        ENDIF
   12 CONTINUE
      CLOSE(LU1)
C
C     Reading field travel times:
      CALL RSEP3T('FTT',FILE1 ,' ')
      OPEN(LU1,FILE=FILE1,STATUS='OLD')
      NT=0
      READ(LU1,*,END=19) TEXT
   13 CONTINUE
        NT=NT+1
        IF(NT.GT.MT) THEN
C         INVTT-07
          CALL ERROR('INVTT-07: Too many field travel times')
        END IF
        SRC=' '
        READ(LU1,*,END=19) SRC,REC,TFIELD(NT),TERR(NT)
        IF(SRC.EQ.' ') THEN
          GO TO 19
        END IF
        DO 14 I=1,NP
          IF(SRC.EQ.POINT(I)) THEN
            KS(NT)=I
            GO TO 15
          END IF
   14   CONTINUE
C         INVTT-08
          ERRTXT(1:38)='INVTT-08: Source name not recognized: '
          ERRTXT(39:49)=SRC
          CALL ERROR(ERRTXT(1:49))
C         Source name used in file FTT with field travel times has not
C         been found in file PTS containing the list of source and
C         receiver points.
   15   CONTINUE
        DO 16 I=1,NP
          IF(REC.EQ.POINT(I)) THEN
            KR(NT)=I
            GO TO 17
          END IF
   16   CONTINUE
C       INVTT-09
        ERRTXT(1:40)='INVTT-09: Receiver name not recognized: '
        ERRTXT(41:51)=REC
        CALL ERROR(ERRTXT(1:51))
C       Receiver name used in file FTT with field travel times has not
C       been found in file PTS containing the list of source and
C       receiver points.
   17   CONTINUE
      GO TO 13
   19 CONTINUE
      NT=NT-1
      CLOSE(LU1)
C
C.......................................................................
C
C     Computing quantities describing objective prior information:
C
C     Opening output log file:
      CALL RSEP3T('INVLOG',FILE1 ,'invlog.out')
      OPEN(LU6,FILE=FILE1)
C
C     Opening input file with output filenames of the CRT program:
      CALL RSEP3T('CRTOUT',FILE4,' ')
      FILE1='r01.out'
      FILE2='s01.out'
      FILE3='s01i.out'
      IF(FILE4.NE.' ') THEN
        OPEN(LU4,FILE=FILE4,STATUS='OLD')
      END IF
C
      KS(NT+1)=NP+1
      KR(NT+1)=NP+1
      TFIELD(NT+1)=0.
      IRAY=0
      IWAVE=0
      NSUM=IPAR(IPAR(IPAR(2)))
      IT0=4*MT+4*MP+NPAR+NM
      LU5=IT0+NT
      IF(LU5.GT.MRAM) THEN
C       INVTT-10
        CALL ERROR('INVTT-10: Too small array RAM')
C       Dimension MRAM of array RAM in include file
C       ram.inc should
C       be increased.
      END IF
      DO 21 I=IT0+1,LU5
        IRAM(I)=0
   21 CONTINUE
C
C     Loop for the files with computed rays
      LINLEN=0
   30 CONTINUE
C       Reading output filenames of the CRT program:
        IF(FILE4.NE.' ') THEN
          READ(LU4,*,END=51) FILE1,FILE2,FILE3
        END IF
        IF(FILE1.EQ.' ') THEN
          GO TO 51
        END IF
        I=INDEX(FILE1,'          ')
        J=INDEX(FILE2,'          ')
        K=INDEX(FILE3,'          ')
        LINLEN=MAX0(LINLEN,I+J+K)
        WRITE(*,'(4A)')
     *    '+INVTT:                               ,  CRT files ',
     *    FILE1(1:I),FILE2(1:J),FILE3(1:K)
        OPEN(LU1,FILE=FILE1,FORM='UNFORMATTED',STATUS='OLD')
        IF(FILE2.EQ.' ') THEN
          IU1=0
          IU2=LU1
        ELSE
          IU1=LU1
          IU2=LU2
          OPEN(LU2,FILE=FILE2,FORM='UNFORMATTED',STATUS='OLD')
        END IF
        OPEN(LU3,FILE=FILE3,FORM='UNFORMATTED',STATUS='OLD')
        IRAYTT=0
C
C       Loop for the points of intersection of rays with the surface
   40   CONTINUE
C         Reading the results of the complete ray tracing
          CALL AP00(IU1,IU2,LU3)
          IF(IPT.LE.1.AND.IRAYTT.NE.0)THEN
C           New ray - recording results for the last ray IRAYTT
C           loop for field travel times - searching for two-point ray
            DO 45 IT=1,NT
              IS=KS(IT)
              IR=KR(IT)
              SI=(COOR(1,IS)-XI1)**2+(COOR(2,IS)-XI2)**2
     *                              +(COOR(3,IS)-XI3)**2
              RI=(COOR(1,IR)-XI1)**2+(COOR(2,IR)-XI2)**2
     *                              +(COOR(3,IR)-XI3)**2
              SE=(COOR(1,IS)-XE1)**2+(COOR(2,IS)-XE2)**2
     *                              +(COOR(3,IS)-XE3)**2
              RE=(COOR(1,IR)-XE1)**2+(COOR(2,IR)-XE2)**2
     *                              +(COOR(3,IR)-XE3)**2
              IF(SE.LE.SI.AND.RI.LE.RE) THEN
C               Interchanging source and receiver points
                SI=SE
                RE=RI
                IS=KR(IT)
                IR=KS(IT)
              END IF
              IF((SI.LE.DIST2.AND.RE.LE.DIST2)) THEN
C               Synthetic ray may correspond to the field travel time
C               check for ray directions near the source
C               (allowable angle deviation +-30 deg: cosine**2=0.750)
C               (allowable angle deviation +-15 deg: cosine**2=0.933)
                TTAUX=(COOR(1,IR)-COOR(1,IS))*(XE1-XI1)
     *               +(COOR(2,IR)-COOR(2,IS))*(XE2-XI2)
     *               +(COOR(3,IR)-COOR(3,IS))*(XE3-XI3)
                IF(TTAUX.GT.0..AND.TTAUX*TTAUX.GT.
     *                         0.933*((COOR(1,IR)-COOR(1,IS))**2
     *                               +(COOR(2,IR)-COOR(2,IS))**2
     *                               +(COOR(3,IR)-COOR(3,IS))**2)
     *                  *((XE1-XI1)**2+(XE2-XI2)**2+(XE3-XI3)**2) ) THEN
C                 Synthetic ray corresponds to the field travel time
                  TTAUX=TT-PI(1)*(COOR(1,IS)-XI1)
     *                    -PI(2)*(COOR(2,IS)-XI2)
     *                    -PI(3)*(COOR(3,IS)-XI3)
     *                    +PE(1)*(COOR(1,IR)-XE1)
     *                    +PE(2)*(COOR(2,IR)-XE2)
     *                    +PE(3)*(COOR(3,IR)-XE3)
C                 Possible minimum synthetic travel time
                  SDIST=SQRT(SI)
                  RDIST=SQRT(RE)
                  J=IRAM(IT0+IT)
                  IF(J.EQ.0) THEN
C                   First two-point ray
                    IF(LU5+4+NM+NLOC.GT.MRAM) THEN
C                     
C                     INVTT-11
                      CALL ERROR('INVTT-11: Too small array RAM')
C                     Dimension MRAM of array RAM in include file
C                     
C                     ram.inc should be
C                     increased.
                    END IF
                    IRAM(IT0+IT)=LU5
                    J=LU5
                    LU5=LU5+4+NM+NLOC
                  ELSE IF(TTAUX.GE.RAM(J+2)) THEN
C                   Synthetic travel time is not minimal
                    GO TO 42
                  END IF
                    RAM(J+1)=TT
                    RAM(J+2)=TTAUX
                    RAM(J+3)=SDIST
                    RAM(J+4)=RDIST
                    J=J+4
                    DO 41 I=1,NM
                      RAM(J+I)=SUM(INDM(I))
   41               CONTINUE
                    IF(NLOC.EQ.4) THEN
                      RAM(J-2)=TTAUX+COOR(4,KS(IT))
                      J=J+NM+1
                      IRAM(J)=KS(IT)
                      IF(KS(IT).EQ.IS) THEN
                        RAM(J+1)=-PI(1)
                        RAM(J+2)=-PI(2)
                        RAM(J+3)=-PI(3)
                      ELSE
                        RAM(J+1)=PE(1)
                        RAM(J+2)=PE(2)
                        RAM(J+3)=PE(3)
                      END IF
                    END IF
   42             CONTINUE
                END IF
              END IF
   45       CONTINUE
            IRAYTT=0
          END IF
          IF(IWAVE.EQ.0) THEN
            GO TO 50
          END IF
C         *** for future extensions (selection of two-point rays):
C         IF(IU1.NE.0) THEN
C           CALL AP30(IREC)
C           IF(IREC.EQ.0) THEN
C             IF(IPT.LE.1.) THEN
C               WRITE(*,'(''+WAVE:'',I3,''    RAY:'',I4,''    POINT:'',
C    *                                              I4)') IWAVE,IRAY,IPT
C             END IF
C             GO TO 40
C           END IF
C         END IF
C         ***
CPTS      IF(IPT.EQ.1.OR.MOD(IPT,10).EQ.0) THEN
CPTS        WRITE(*,'(''+INVTT: +WAVE:'',I3,''    RAY:'',I4,
CPTS *                               ''    POINT:'',I4)') IWAVE,IRAY,IPT
CPTS      END IF
          IF(IPT.EQ.1) THEN
            WRITE(*,'(A,I3,A,I3,A,I5)')
     *               '+INVTT: Source',ISRC,',  Wave',IWAVE,',  Ray',IRAY
C                   38 characters
          END IF
          IRAYTT=IRAY
          XI1=YI(3)
          XI2=YI(4)
          XI3=YI(5)
          XE1=Y(3)
          XE2=Y(4)
          XE3=Y(5)
          CALL AP01(TT,TTAUX)
          CALL AP02(PI,PE)
          CALL AP29(NSUM,SUM)
        GO TO 40
C       End of the loop for points of intersection of rays with surface
C
   50   CONTINUE
        CLOSE(LU1)
        IF(FILE2.NE.' ') THEN
          CLOSE(LU2)
        END IF
        CLOSE(LU3)
        FILE1=' '
      GO TO 30
C
   51 CONTINUE
      IF(FILE4.NE.' ') THEN
        CLOSE(LU4)
      END IF
C
C     All minimum travel times and their derivatives are stored in RAM
C
C.......................................................................
C
C     Writing objective prior information:
      WRITE(*,'(A)') '+INVTT: Writing matrices.             '
C
C     List of sources to be located:
      NSRC=0
      IF(NLOC.EQ.4) THEN
C       Old source points
        M1IN=NM
        CALL RSEP3T('M1IN',FILE1,' ')
        IF(FILE1.NE.' ') THEN
          OPEN(LU1,FILE=FILE1,STATUS='OLD')
          READ(LU1,*) M1IN
          CLOSE(LU1)
        END IF
        NSRCIN=(M1IN-NM)/NLOC
        IF(LU5+NSRCIN.GT.MRAM) THEN
C         INVTT-12
          CALL ERROR('INVTT-12: Too small dimension of array RAM')
C         Dimension MRAM of array RAM in include file
C         ram.inc
C         should be increased.
        END IF
        CALL RSEP3T('INVSRC',FILE1 ,' ')
        OPEN(LU1,FILE=FILE1)
        DO 60 J=1,NSRCIN
          READ(LU1,*) SRC
          DO 54 I=1,NP
            IF(SRC.EQ.POINT(I)) THEN
              IRAM(LU5+J)=I
              GO TO 55
            END IF
   54     CONTINUE
C           INVTT-13
            CALL ERROR('INVTT-13: Wrong source to be located')
   55     CONTINUE
   60   CONTINUE
C       NSRC=NSRCIN
C       New source points
        DO 63 IT=1,NT
          J=IRAM(IT0+IT)
          IF(J.GT.0) THEN
            J=J+4+NM+1
            IF(J.LT.1.OR.J.GT.MRAM) THEN
C             INVTT-14
              CALL ERROR('INVTT-14: Wrong index in array RAM')
C             This error should not appear.  Contact the author.
            END IF
            DO 61 I=1,NSRC
              IF(IRAM(J).EQ.IRAM(LU5+I)) THEN
                IRAM(J)=I
                GO TO 62
              END IF
   61       CONTINUE
C             New source point
              NSRC=NSRC+1
              IF(LU5+NSRC.GT.MRAM) THEN
C               INVTT-15
                CALL ERROR('INVTT-15: Too small dimension of array RAM')
C               Dimension MRAM of array RAM in include file
C               ram.inc
C               should be increased.
              END IF
              IRAM(LU5+NSRC)=IRAM(J)
              IRAM(J)=NSRC
   62       CONTINUE
          END IF
   63   CONTINUE
        DO 64 I=NSRCIN+1,NSRC
              IF(IRAM(LU5+I).LT.1.OR.IRAM(LU5+I).GT.MP) THEN
                write(*,'(6I12)') NSRCIN,I
                write(*,'(6I12)') (IRAM(LU5+Ii),Ii=1,NSRC)
C               INVTT-16
                CALL ERROR('INVTT-16: Wrong source index')
C               This error should not appear.  Contact the author.
              END IF
          WRITE(LU1,'(3A)') '''',POINT(IRAM(LU5+I)),''''
   64   CONTINUE
        CLOSE(LU1)
      END IF
C
C     Writing the number of model parameters:
      CALL RSEP3T('M1',FILE1,' ')
      IF(FILE1.NE.' ') THEN
        OPEN(LU1,FILE=FILE1)
        WRITE(LU1,'(I10)') NM+4*NSRC
        M1=NM+4*NSRC
        CLOSE(LU1)
      END IF
C
C     Reading number of travel times processed previously:
      M2IN=0
      CALL RSEP3T('M2IN',FILE1,' ')
      IF(FILE1.NE.' ') THEN
        OPEN(LU1,FILE=FILE1,STATUS='OLD')
        READ(LU1,*) M2IN
        CLOSE(LU1)
      END IF
      IF(LU5+(NM+4*NSRCIN)*M2IN.GT.MRAM) THEN
C       INVTT-17
        CALL ERROR('INVTT-17: Too small dimension MRAM of array RAM')
C       Dimension MRAM of array RAM in include file
C       ram.inc
C       should be increased.
      END IF
C
C     Determining number ND of equations:
      ND=0
      DO 65 I=IT0+1,IT0+NT
        IF(IRAM(I).GT.0) THEN
          ND=ND+1
        END IF
   65 CONTINUE
C
C     Writing the total number of equations:
      M2=M2IN+ND
      CALL RSEP3T('M2',FILE1,'m2.out')
      IF(FILE1.NE.' ') THEN
        OPEN(LU1,FILE=FILE1)
        WRITE(LU1,'(I10)') M2
        CLOSE(LU1)
      END IF
C
C     Writing matrix header files:
      CALL RSEP3T('GM1',FILE1 ,' ')
      CALL RSEP3T('GM2',FILE2 ,' ')
      CALL RSEP3T('GM3',FILE3 ,' ')
      CALL RSEP3T('DM1',FILE4 ,' ')
      IF (FILE1.NE.' ') THEN
        FILED=' '
        CALL WMATH(LU1,FILE1,FILED,M1,M2,' ',0,' ','formatted')
        FILE1=FILED
      ENDIF
      IF (FILE2.NE.' ') THEN
        FILED=' '
        CALL WMATH(LU1,FILE2,FILED,M2,1,' ',0,' ','formatted')
        FILE2=FILED
      ENDIF
      IF (FILE3.NE.' ') THEN
        FILED=' '
        CALL WMATH(LU1,FILE3,FILED,M2,1,' ',0,' ','formatted')
        FILE3=FILED
      ENDIF
      IF (FILE4.NE.' ') THEN
        FILED=' '
        CALL WMATH(LU1,FILE4,FILED,M2,M2,' ',0,'diag','formatted')
        FILE4=FILED
      ENDIF
C
C     Opening input/output files with matrix components:
c     CALL RSEP3T('GM1',FILE1 ,' ')
      IF(FILE1.NE.' ') THEN
        IF(M2IN.GT.0) THEN
          CALL RARRAY(LU1,FILE1,'FORMATTED',.TRUE.,0.,
     *                                    (NM+4*NSRCIN)*M2IN,RAM(LU5+1))
        END IF
        OPEN(LU1,FILE=FILE1)
        IF(M2IN.GT.0) THEN
          DO 71 J=LU5,LU5+(NM+4*NSRCIN)*(M2IN-1),NM+4*NSRCIN
            WRITE(LU1,FORMAT) (RAM(I),I=J+1,J+NM+4*NSRCIN),
     *                        (0.,I=4*NSRCIN+1,4*NSRC)
   71     CONTINUE
        END IF
      END IF
c     CALL RSEP3T('GM2',FILE2 ,' ')
      IF(FILE2.NE.' ') THEN
        IF(M2IN.GT.0) THEN
          CALL RARRAY(LU2,FILE2,'FORMATTED',.TRUE.,0.,
     *                                                  M2IN,RAM(LU5+1))
        END IF
        OPEN(LU2,FILE=FILE2)
        IF(M2IN.GT.0) THEN
          DO 72 J=LU5+1,LU5+M2IN
            WRITE(LU2,FORMAT) RAM(J)
   72     CONTINUE
        END IF
      END IF
c     CALL RSEP3T('GM3',FILE3 ,' ')
      IF(FILE3.NE.' ') THEN
        IF(M2IN.GT.0) THEN
          CALL RARRAY(LU3,FILE3,'FORMATTED',.TRUE.,0.,
     *                                                  M2IN,RAM(LU5+1))
        END IF
        OPEN(LU3,FILE=FILE3)
        IF(M2IN.GT.0) THEN
          DO 73 J=LU5+1,LU5+M2IN
            WRITE(LU3,FORMAT) RAM(J)
   73     CONTINUE
        END IF
      END IF
c     CALL RSEP3T('DM1',FILE4 ,' ')
      IF(FILE4.NE.' ') THEN
        IF(M2IN.GT.0) THEN
          CALL RARRAY(LU4,FILE4,'FORMATTED',.TRUE.,0.,
     *                                                  M2IN,RAM(LU5+1))
        END IF
        OPEN(LU4,FILE=FILE4)
        IF(M2IN.GT.0) THEN
          DO 74 J=LU5+1,LU5+M2IN
            WRITE(LU4,FORMAT) RAM(J)
   74     CONTINUE
        END IF
      END IF
C
C     Writing input/output files with matrix components:
      WRITE(LU6,'(2A)') ' SOURCE      RECEIVER       TFIELD    ',
     *                  'TFIELD-TT     SDIST       RDIST    TT-CHECKSUM'
      DO 79 IT=1,NT
        J=IRAM(IT0+IT)
        IF(J.GT.0) THEN
          TT   =RAM(J+1)
          TTAUX=RAM(J+2)
          SDIST=RAM(J+3)
          RDIST=RAM(J+4)
          J=J+4
          DO 75 I=1,NM
            SUM(I)=RAM(J+I)
   75     CONTINUE
          J=J+NM+1
C
C         System of linear equations:
          TDIF=TFIELD(IT)-TTAUX
          IF((FILE1.NE.' ').AND.(NM+4*NSRC.GT.0)) THEN
            WRITE(LU1,FORMAT) (SUM(I),I=1,NM)
            IF(NLOC.EQ.4) THEN
              DO 76 I=1,IRAM(J)-1
                WRITE(LU1,FORMAT) 0.,0.,0.,0.
   76         CONTINUE
              WRITE(LU1,FORMAT) (RAM(I),I=J+1,J+3),1.
              DO 77 I=IRAM(J)+1,NSRC
                WRITE(LU1,FORMAT) 0.,0.,0.,0.
   77         CONTINUE
            END IF
          END IF
          IF(FILE2.NE.' ') THEN
            WRITE(LU2,FORMAT) TDIF
          END IF
          IF(FILE3.NE.' ') THEN
            WRITE(LU3,FORMAT) TTAUX
          END IF
          IF(FILE4.NE.' ') THEN
            WRITE(LU4,FORMAT) TERR(IT)**2
          END IF
C
C         Check sums and log output:
          IS=KS(IT)
          IR=KR(IT)
          IF(VPOWER.NE.0.) THEN
            TTAUX=0.
            DO 78 I=1,NM
              J=INDM(I)
              IF(IPAR(IPAR(IPAR(1))).LT.J) THEN
                IF(SUM(I).NE.0.) THEN
                  TTAUX=TTAUX+RPAR(J)*SUM(I)
                END IF
              END IF
   78       CONTINUE
            TTAUX=TT+VPOWER*TTAUX
            WRITE(LU6,'(2(1X,A),5F12.6)') POINT(IS),POINT(IR),
     *           TFIELD(IT),TDIF,SDIST,RDIST,TTAUX
          ELSE
            WRITE(LU6,'(2(1X,A),5F12.6)') POINT(IS),POINT(IR),
     *           TFIELD(IT),TDIF,SDIST,RDIST
          END IF
        END IF
   79 CONTINUE
C
      IF(FILE1.NE.' ') THEN
        CLOSE(LU1)
      END IF
      IF(FILE2.NE.' ') THEN
        CLOSE(LU2)
      END IF
      IF(FILE3.NE.' ') THEN
        CLOSE(LU3)
      END IF
      IF(FILE4.NE.' ') THEN
        CLOSE(LU4)
      END IF
      CLOSE(LU6)
      LINLEN=MIN0(LINLEN,81)
      FILE1=' '
      WRITE(*,'(2A)')
     *    '+INVTT: Done.                                      ',
     *    FILE1(1:LINLEN-1)
      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
      INCLUDE 'mat.for'
C     mat.for
      INCLUDE 'modelv.for'
C     modelv.for
      INCLUDE 'metric.for'
C     metric.for
      INCLUDE 'srfc.for'
C     srfc.for
      INCLUDE 'parmv.for'
C     parmv.for
      INCLUDE 'valv.for'
C     valv.for
      INCLUDE 'fitv.for'
C     fitv.for
      INCLUDE 'var.for'
C     var.for
      INCLUDE 'spsp.for'
C     spsp.for
      INCLUDE 'soft.for'
C     soft.for
      INCLUDE 'means.for'
C     means.for
      INCLUDE 'ap.for'
C     ap.for
      INCLUDE 'apvar.for'
C     apvar.for
C
C=======================================================================
C