C
C Program INVTT to calculate the derivatives of travel times
C with respect to the model B-spline coefficients
C
C Version: 5.50
C Date: 2001, January 20
C
C Coded by: Ludek Klimes
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     E-mail: klimes@seis.karlov.mff.cuni.cz
C
C.......................................................................
C
C Program INVTT assumes all model parameters (coefficients) stored in
C the common block /VALC/ as in the submitted versions of user-defined
C model specification FORTRAN77 source code files 'srfc.for', 'parm.for'
C and 'val.for'.  Thus, unlike the other parts of the complete ray
C tracing, the INVTT program cannot work with user's modifications of
C subroutines SRFC1, SRFC2, PARM1, and PARM2.
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 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
C             file 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     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, PTS must be specified and cannot be blank.
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 M1*M2IN input values of file GM1 and after
C             M2IN input values of files GM2, GM3 and DM1, 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 of
C             model coefficients to be determined during inversion.
C             The same file, containing a single integer, may be
C             generated by programs 'invsoft.for' and 'invpts.for'.
C             The file is not generated if the value of M1 is blank.
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     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 input/output files with matrix and vector components:
C     GM1='string'... String with the name of the input/output data file
C             to accumulate the matrix of the derivatives of M2 travel
C             times with respect to M1 model coefficients.
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 data 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 data 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 data 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     For general description of the files with matrices refer to file
C     forms.htm.
C Other numerical parameters:
C     DIST=real... Maximum distance between the source or receiver point
C             and the initial or end point of a synthetic ray.
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 coefficients (e.g. B-splines), the
C             travel time minus its check sum (see the log output file)
C             should be zero within rounding errors.  Otherwise, the
C             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
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 (3) / (a slash) or the end of file.
C Example of data 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.  The source and
C             receiver points may be mutually interchanged.
C     TFIELD..Travel time from a field measurement.
C     TERR... Error of the travel time from a field measurement.
C (3) / (a slash) or the end of file.
C Example of data 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 coefficients.  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(3,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+3*MP+1))
      EQUIVALENCE (INDM  ,IRAM(4*MT+3*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
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*1  TEXT(10)
      CHARACTER*80 ERRTXT
      CHARACTER*11 SRC,REC
      INTEGER NP,NT
      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
C     NM...   Number of the unknown model parameters.
C
C     Auxiliary storage locations:
      INTEGER IS,IR,IT,ND,IRAYTT,I,J,K
      REAL TT,TTAUX,TDIF,XI1,XI2,XI3,XE1,XE2,XE3,PI(6),PE(6)
      REAL SI,SE,RI,RE,SDIST,RDIST
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(G13.7,1X))'
C
      IF(4*MT+3*MP+2*NPAR.GT.MRAM) THEN
C       INVTT-06
        CALL ERROR('INVTT-06: Too small dimension MRAM of array RAM')
C       Dimension MRAM of array RAM in include file
C       ram.inc should 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-01
        CALL ERROR('INVTT-01: 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:
      DO 10 I=1,47
        W(I,1)=0.
        W(I,2)=0.
   10 CONTINUE
      CALL SOFT(2,0,0,0,0,0,0,47,W,NM,INDM,CS,1,CS)
      WRITE(*,'(A,I4,A)') '+INVTT:',NM,' model parameters'
C
C     Reading numerical parameters:
      CALL RSEP3R('DIST'  ,DIST  ,0.)
      CALL RSEP3R('VPOWER',VPOWER,0.)
      DIST2=DIST*DIST
C
C.......................................................................
C
C     Reading source and receiver points:
      CALL RSEP3T('PTS',FILE1 ,' ')
      OPEN(LU1,FILE=FILE1,STATUS='OLD')
      NP=-1
      READ(LU1,*,END=12) TEXT
   11 CONTINUE
        NP=NP+1
        IF(NP.GE.MP) THEN
C         INVTT-04
          CALL ERROR('INVTT-04: Too many source and receiver points')
        END IF
        POINT(NP+1)=' '
        READ(LU1,*,END=12) POINT(NP+1),(COOR(I,NP+1),I=1,3)
      IF(POINT(NP+1).NE.' ') GO TO 11
   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=12) TEXT
   13 CONTINUE
        NT=NT+1
        IF(NT.GT.MT) THEN
C         INVTT-05
          CALL ERROR('INVTT-05: 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-02
          ERRTXT(1:38)='INVTT-02: 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-03
        ERRTXT(1:40)='INVTT-03: 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.dat')
      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+3*MP+NPAR+NM
      LU5=IT0+NT
      IF(LU5.GT.MRAM) THEN
C       
C       INVTT-08
        CALL ERROR('INVTT-08: Too small array RAM')
C       Dimension MRAM of array RAM in include file
C       ram.inc should be increased.
      END IF
      DO 21 I=IT0+1,LU5
        IRAM(I)=0
   21 CONTINUE
C
C     Loop for the files with computed rays
   30 CONTINUE
C       Reading output filenames of the CRT program:
        IF(FILE4.NE.' ') THEN
          READ(LU4,*,END=70) FILE1,FILE2,FILE3
        END IF
        IF(FILE1.EQ.' ') THEN
          GO TO 70
        END IF
        I=INDEX(FILE1,'          ')
        J=INDEX(FILE2,'          ')
        K=INDEX(FILE3,'          ')
        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.GT.MRAM) THEN
C                     
C                     INVTT-07
                      CALL ERROR('INVTT-07: Too small array RAM')
C                     Dimension MRAM of array RAM in include file
C                     ram.inc should be
C                     increased.
                    END IF
                    IRAM(IT0+IT)=LU5
                    J=LU5
                    LU5=LU5+4+NM
                  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
   42             CONTINUE
                END IF
              END IF
   45       CONTINUE
            IRAYTT=0
          END IF
          IF(IWAVE.EQ.0) THEN
            GO TO 60
          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
   60   CONTINUE
        CLOSE(LU1)
        IF(FILE2.NE.' ') THEN
          CLOSE(LU2)
        END IF
        CLOSE(LU3)
        FILE1=' '
      GO TO 30
C
   70 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     Writing the number of model coefficients:
      CALL RSEP3T('M1',FILE1,' ')
      IF(FILE1.NE.' ') THEN
        OPEN(LU1,FILE=FILE1)
        WRITE(LU1,'(I10)') NM
        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*M2IN.GT.MRAM) THEN
C       INVTT-09
        CALL ERROR('INVTT-09: Too small dimension MRAM of array RAM')
C       Dimension MRAM of array RAM in include file
C       ram.inc should be increased.
      END IF
C
C     Determining number ND of equations:
      ND=0
      DO 71 I=IT0+1,IT0+NT
        IF(IRAM(I).GT.0) THEN
          ND=ND+1
        END IF
   71 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     Opening input/output files with matrix components:
      CALL RSEP3T('GM1',FILE1 ,' ')
      IF(FILE1.NE.' ') THEN
        IF(M2IN.GT.0) THEN
          CALL RARRAY(LU1,FILE1,'FORMATTED',.TRUE.,0.,
     *                                               NM*M2IN,RAM(LU5+1))
        END IF
        OPEN(LU1,FILE=FILE1)
        IF(M2IN.GT.0) THEN
          DO 72 J=LU5,LU5+NM*(M2IN-1),NM
            WRITE(LU1,FORMAT) (RAM(I),I=J+1,J+NM)
   72     CONTINUE
        END IF
      END IF
      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 73 J=LU5+1,LU5+M2IN
            WRITE(LU2,FORMAT) RAM(J)
   73     CONTINUE
        END IF
      END IF
      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 74 J=LU5+1,LU5+M2IN
            WRITE(LU3,FORMAT) RAM(J)
   74     CONTINUE
        END IF
      END IF
      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 75 J=LU5+1,LU5+M2IN
            WRITE(LU4,FORMAT) RAM(J)
   75     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 77 I=1,NM
            SUM(I)=RAM(J+I)
   77     CONTINUE
C
C         System of linear equations:
          TDIF=TFIELD(IT)-TTAUX
          IF(FILE1.NE.' ') THEN
            WRITE(LU1,FORMAT) (SUM(I),I=1,NM)
          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:
          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
            IS=KS(IT)
            IR=KR(IT)
            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)
      WRITE(*,'(A)') '+INVTT: 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
      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