C
C=======================================================================
C
      PROGRAM MTT
C     to calculate Multivalued Travel Times on a 3-D grid of points
C
C----------------------------------------------------------------------
C
C The post processing program to the CRT program performs interpolation
C of travel times to the gridpoints of arbitrary regular rectangular
C 3-D grid of points. The interpolation inside ray tubes formed by
C vertices of homogeneous triangles created during 3-D two-point
C ray tracing is realized using the controlled initial-value ray tracing
C algorithm.
C If the I-th and (I+1)-th points of the three rays forming the ray tube
C are situated inside the same complex block, the program performes
C interpolation within ray cells whose bottom is formed by the I-th
C points and top by (I+1)-th points of the rays.  The best results can
C thus be obtained if the travel time is the independent variable for
C numerical integration in CRT program (default).
C
C This version of the program interpolates travel times within ray cells
C trilinearly (in the coordinate system connected with each cell),
C causing large travel-time errors inside wide ray cells.  This should
C be accounted for when interpolating travel times in the areas of large
C geometrical spreading.
C
C The boundaries between ray cells are determined numerically, with
C different rounding errors in different cells.  The rounding errors
C produce narrow strips along the boundaries between the cells, where
C the travel times of the corresponding history are randomly missing
C or doubled.  If a gridpoint falls, by a chance, into the strip,
C the corresponding arrival at the gridpoint may be missing or doubled.
C A small percentage of gridpoints with missing or doubled arrivals is
C thus inevitable with the present interpolation algorithm.
C
C The interpolation is not performed in demarcation belts between
C different ray histories.  The maximum width of the belts is controlled
C by input parameter AERR of the CRT
C program.  The typical width of the belts, measured in take-off angles,
C is 0.0001 radians in order of magnitude.  The belts may become wide
C in areas of large geometrical spreading.  The division of rays into
C different histories is, by default, very detailed and is controlled
C by input parameter PRM0(3) of the CRT
C program.  Such a behaviour is useful for two-point ray tracing but has
C some awkward consequences on the interpolation.  For example, the
C demarcation belt between the rays incident at different sides of the
C computational volume for ray tracing extends across the whole model,
C causing an artificial gap within the continuous travel times.
C If the ray history does not account for the termination of a ray at
C different sides of the computational volume, the gaps corresponding
C to the edges of the computational volume are removed, but the corners
C along the edges are not filled with the ray cells any more.  Then the
C computational volume for ray tracing should sufficiently exceed the
C extent of target volume covered by the grid for interpolation.
C Input parameter PRM0(3) of the CRT
C program may thus considerably influence the results of the
C interpolation.  However, gaps due to the demarcation belts
C corresponding to structural interfaces cannot be removed within
C the present interpolation algorithm.
C
C.......................................................................
C
C                                                    
C Main input data read from external interactive device (*):
C     The data consist of character strings, read by list
C     directed (free format) input.  The strings have thus to be
C     enclosed in apostrophes.
C 'SEP','INPUTFILE',/
C     'SEP'...  Name of the file describing the configuration of
C               the target grid, where the output quantities are
C               to be interpolated.
C               Description of file SEP
C     'INPUTFILE' ...  Name of file with input data.
C Defaults: 'SEP'='grd.h', 'INPUTFILE'='mtt.dat'
C
C Structure of the input data file 'INPUTFILE':
C     The data consist of following lines of
C     character strings, read by list directed (free format) input:
C (1) 'CRT-R','CRT-I','CRT-T',/
C     'CRT-R'...Name of the input unformatted file with the quantities
C               stored along rays (see C.R.T.5.5.1).
C               Description of file CRT-R
C     'CRT-I'...Name of the input unformatted file with the quantities
C               at the initial points of rays, corresponding to file
C               CRT-R (see C.R.T.6.1).
C               Description of file CRT-I
C     'CRT-T'...Name of the input formatted file with the indices of
C               rays forming the homogeneous triangles in the
C               ray-parameter domain.
C               Description of file CRT-T
C     Defaults: 'CRT-R'='r01.out', 'CRT-I'='r01i.out',
C               'CRT-T'='t01.out'.
C (2) 'NUM',/
C     'NUM'...  Name of the output formatted file with N1*N2*N3 array of
C               integer values, containing the number of arrivals at
C               each gridpoint of the target grid.
C               Description of file NUM
C     Default:  'NUM'='num.out'
C (3) Zero to MOUT times following line:
C     'FOUT(i)','OUT(i)',/
C     'FOUT(i)'.Name of the output formatted file with the array of
C               values, containing for each gridpoint the quantities
C               of all determined arrivals.  The number of values
C               thus equals the sum of the integers in file 'NUM'.
C               The kind of computed quantities written into this
C               file is given by value of 'OUT(i)'.
C               Description of files 'FOUT(i)'
C     'OUT(i)'..Character string, describing which quantities are to be
C               written into corresponding output file 'FOUT(i)'.
C               Following is a list of possible values of 'OUT(i)' and
C               corresponding quantities written into output files:
C               'tt'   .......  Travel times
C               'hi'   .......  Ray histories
C     Defaults: 'FOUT(1)'   ='mtt.out',  'OUT(1)'  ='tt',
C               'FOUT(2)'   =' '      ,  'OUT(2)'  =' ' ,
C                   .                      .
C                   .                      .
C               'FOUT(MOUT)'=' '      ,  'OUT(MOUT)=' ' .
C
C
C                                                     
C Input formatted input SEP parameter file:
C     The file with coordinates of the origin of the target grid
C     O1, O2, O3, numbers of gridpoints N1, N2, N3 and spatial step
C     between gridpoints D1, D2, D3. The file has the form of SEP file,
C     for the description of the SEP format refer
C     to the file  'sep.for'.
C     N1=positive integer... Number of gridpoints along the X1 axis.
C             Default: N1=1
C     N2=positive integer... Number of gridpoints along the X2 axis.
C             Default: N2=1
C     N3=positive integer... Number of gridpoints along the X3 axis.
C             Default: N3=1
C     O1=real, O2=real, O3=real... Coordinates of the first point of the
C             grid.  Default: O1=0, O2=0, O3=0.
C     D1=real... Grid interval along the X1 axis.  Default: D1=1.
C     D2=real... Grid interval along the X2 axis.  Default: D2=1.
C     D3=real... Grid interval along the X3 axis.  Default: D3=1.
C Example of data SEP
C
C                                                     
C Output formatted file 'NUM':
C     N1*N2*N3 integers:  For each gridpoint, the number of travel
C     times determined at the gridpoint.  The innermost loop corresponds
C     to N1, the outermost loop corresponds to N3.
C
C                                                   
C Output formatted file 'FOUT(i)':
C     N numbers, where N is the sum of integers in file NUM.
C     For each gridpoint, all quantities determined at the gridpoint.
C     The kind of quantities is determined by corresponding
C     parameter OUT(i), for example if OUT(1)='TT', then travel times
C     are written into output file named 'FOUT(1)'.
C
C.......................................................................
C
C Date: 1997, October 25
C                                                  
C Coded by Petr Bulant
C     bulant@seis.karlov.mff.cuni.cz
C
C-----------------------------------------------------------------------
C
C Common block /RAM/ to store the information about triangles and rays:
      INCLUDE 'mtt.inc'
C ...........................
C Common block /CIGRD/ to store the parameters of the
C interpolation grid.
      REAL O1,O2,O3,D1,D2,D3
      INTEGER N1,N2,N3
      COMMON/CIGRD/O1,O2,O3,D1,D2,D3,N1,N2,N3
      SAVE /CIGRD/
C     O1,O2,O3 ... Coordinates of the origin of the grid.
C     D1,D2,D3 ... Steps per gridpoint.
C     N1,N2,N3 ... Numbers of gridpoints.
C.......................................................................
C     External functions required:
C
C     Auxiliary storage locations:
      INTEGER IRAY1,IRAY2,IRAY3
      INTEGER I1,I2,I3,I4,I5,INDX,IT1
      INTEGER I,J,K,L,M,N
      INTEGER LU0,LU1,LU2,LU3,LU4,LU5
      PARAMETER (LU0=10,LU1=1,LU2=2,LU3=3,LU4=4,LU5=5)
      CHARACTER*80 FILE0,FILE1,FILE2,FILE3,FILE4,FILE5
C
C-----------------------------------------------------------------------
C
      WRITE(*,'(''+'',79('' ''))')
      WRITE(*,'(A)') '+Reading input data.'
C     Reading in a name of the file with the data for interpolation
C     grid, and a name of the main input data file:
      FILE4='grd.h'
      FILE0='mtt.dat'
      WRITE(*,'(A)')
     *  ' Enter the names of the input data files (SEP,INPUTFILE): '
      READ(*,*) FILE4,FILE0
      OPEN(LU0,FILE=FILE0,FORM='FORMATTED',STATUS='OLD')
C     Reading in filenames of the files with
C     computed triangles and rays:
      FILE1='r01.out'
      FILE2='r01i.out'
      FILE3='t01.out'
      READ(LU0,*) FILE1,FILE2,FILE3
      OPEN(LU1,FILE=FILE1,FORM='UNFORMATTED',STATUS='OLD')
      OPEN(LU2,FILE=FILE2,FORM='UNFORMATTED',STATUS='OLD')
      OPEN(LU3,FILE=FILE3,FORM='FORMATTED',STATUS='OLD')
C     Reading in filenames of the output files together
C     with quantities describing what is to be written into the files:
      OUT(1)='tt'
      OUT(2)=' '
      OUT(3)=' '
      OUT(4)=' '
      OUT(5)=' '
      FILE5='num.out'
      FOUT(1)='mtt.out'
      FOUT(2)=' '
      FOUT(3)=' '
      FOUT(4)=' '
      FOUT(5)=' '
      READ(LU0,*) FILE5
      I1=1
   2  CONTINUE
        READ(LU0,*,END=4) FOUT(I1),OUT(I1)
        I1=I1+1
        IF (I1.GT.MOUT) THEN
C         MTT-19
          PAUSE 'Error MTT-19: Small array FOUT for the output files.'
          STOP
C         This error may be caused by wrongly prescribed
C          input data FOUT and OUT, or by
C         small parameter MOUT defined at the top of this file.
        ENDIF
      GOTO 2
   4  CONTINUE
C     Checking the requests for output files:
      NOUT=0
      DO 6, I1=1,MOUT
        CALL LOWER(OUT(I1))
        IF ((OUT(I1)(1:2).NE.' ').AND.(OUT(I1)(1:2).NE.'tt').AND.
     *      (OUT(I1)(1:2).NE.'hi')) THEN
C         MTT-20
          PAUSE 'Error MTT-20: Wrong value of OUT.'
          STOP
C         This error may be caused by wrongly prescribed value of
C          input data OUT.
        ELSEIF (((OUT(I1)(1:2).EQ.' ').AND.(FOUT(I1).NE.' ')) .OR.
     *          ((OUT(I1)(1:2).NE.' ').AND.(FOUT(I1).EQ.' '))) THEN
C         MTT-21
          PAUSE 'Error MTT-21: Inconsistency between FOUT and OUT.'
          STOP
C         Either both of FOUT(i) and OUT(i) or none of them must be
C         prescribed. See the description of
C          input data FOUT and OUT.
        ELSEIF ((OUT(I1)(1:2).NE.' ').AND.(FOUT(I1).NE.' ')) THEN
          NOUT=NOUT+1
          FOUT(NOUT)=FOUT(I1)
          OUT(NOUT)=OUT(I1)(1:2)
        ENDIF
   6  CONTINUE
      NQ=5+NOUT
C
C     Reading in the file with parameters of target grid:
      CALL CIGRID(LU4,FILE4)
C
C
C     Reading in file with computed triangles,
C     sorting the rays in triangles:
      NT=0
      NRAMP=0
      IRAYMI=999999
      IRAYMA=0
  10  CONTINUE
        READ(LU3,*,END=20) IRAY1,IRAY2,IRAY3
        IF (IRAY1.GT.IRAYMA) IRAYMA=IRAY1
        IF (IRAY2.GT.IRAYMA) IRAYMA=IRAY2
        IF (IRAY3.GT.IRAYMA) IRAYMA=IRAY3
        IF (IRAY1.LT.IRAYMI) IRAYMI=IRAY1
        IF (IRAY2.LT.IRAYMI) IRAYMI=IRAY2
        IF (IRAY3.LT.IRAYMI) IRAYMI=IRAY3
        IF (IRAY1.LT.IRAY2) THEN
          I=IRAY1
          IRAY1=IRAY2
          IRAY2=I
        ENDIF
        IF (IRAY2.LT.IRAY3) THEN
          I=IRAY2
          IRAY2=IRAY3
          IRAY3=I
        ENDIF
        IF (IRAY1.LT.IRAY2) THEN
          I=IRAY1
          IRAY1=IRAY2
          IRAY2=I
        ENDIF
        NRAMP=NRAMP+1
        IF (NRAMP.GT.MRAMP) CALL CIEROR(1)
        IRAM(NRAMP)=IRAY1
        NRAMP=NRAMP+1
        IF (NRAMP.GT.MRAMP) CALL CIEROR(1)
        IRAM(NRAMP)=IRAY2
        NRAMP=NRAMP+1
        IF (NRAMP.GT.MRAMP) CALL CIEROR(1)
        IRAM(NRAMP)=IRAY3
        NT=NT+1
      GOTO 10
  20  CONTINUE
      CLOSE(LU3)
      NR=IRAYMA-IRAYMI+1
C
C     Sorting the triangles:
      DO 40, I1=NRAMP-5,1,-3
        DO 30, I2=1,I1,3
          L=I2+3
          IF (IRAM(I2).GT.IRAM(L)) THEN
            J=I2+1
            K=I2+2
            M=I2+4
            N=I2+5
            I       =IRAM(I2)
            IRAM(I2)=IRAM(L)
            IRAM(L) =I
            I       =IRAM(J)
            IRAM(J) =IRAM(M)
            IRAM(M) =I
            I       =IRAM(K)
            IRAM(K) =IRAM(N)
            IRAM(N) =I
          ENDIF
  30    CONTINUE
  40  CONTINUE
C
C
C     Forming an auxiliary array with information about when can be rays
C     erased from the memory ("deleting array"):
      IF (NRAMP+NR.GT.MRAMP) CALL CIEROR(1)
      DO 50, I1=NRAMP+1,NRAMP+NR
        IRAM(I1)=0
  50  CONTINUE
      NRAMP=NRAMP+NR
      ORAYE=IRAYMI-1-3*NT
      DO 60, I1=1,3*NT,3
        IRAM(IRAM(I1  )-ORAYE)=IRAM(I1)
        IRAM(IRAM(I1+1)-ORAYE)=IRAM(I1)
        IRAM(IRAM(I1+2)-ORAYE)=IRAM(I1)
  60  CONTINUE
C
C
C     Forming an auxiliary array with information about addresses
C     of the ends of records for rays in array RAM ("addressing array"):
C     "Ray" IRAYMI-1:
      NRAMP=NRAMP+1
      IF (NRAMP.GT.MRAMP) CALL CIEROR(1)
      IRAM(NRAMP)=NRAMP+NR
C     All other rays:
      IF (NRAMP+NR.GT.MRAMP) CALL CIEROR(1)
      DO 70, I1=NRAMP+1,NRAMP+NR
        IRAM(I1)=0
  70  CONTINUE
      NRAMP=NRAMP+NR
      ORAYA=IRAYMI-1-3*NT-NR-1
C
C
C     Loop for all the triangles:
      IT1=1
      I1=-1
  100 CONTINUE
        I2=NINT(100*IT1/3/NT)
        IF (I2.NE.I1) THEN
          WRITE(*,'(''+'',79('' ''))')
          WRITE(*,'(A,1I4,A)') '+Interpolating. ',I2,'% completed.'
          I1=I2
        ENDIF
C
C       If necessary, reading in new rays:
        IF
     *  (((IRAM(IRAM(IT1)-ORAYA+1).EQ.0).AND.(IRAM(IT1).NE.IRAYMA)) .OR.
     *   ((IRAM(IRAM(IT1)-ORAYA  ).EQ.0).AND.(IRAM(IT1).EQ.IRAYMA)))
     *  CALL CIREAD(LU1,LU2,IT1)
C
C       Empting the array RAM to enable writing of possible interpolated
C       quantities:
        IF ((MRAMP-NRAMP).LT.(NINT(MRAM/10))) CALL CIERAS
C
C       Interpolation along the rays of the triangle:
        CALL CITUBE(IT1)
C
      IT1=IT1+3
      IF (IT1.LT.3*NT) GOTO 100
C     End of the loop for all the triangles.
      CLOSE(LU1)
      CLOSE(LU2)
C
C
C     Sorting the results according to the travel times:
      DO 340, I4=1,NOUT
        IF (OUT(I4)(1:2).EQ.'tt') THEN
          DO 330, I3=0,N3-1
          DO 320, I2=0,N2-1
          DO 310, I1=0,N1-1
  300       CONTINUE
            INDX=I3*N1*N2+I2*N1+I1
            INDX=MRAM-INDX
            INDX1=IRAM(INDX)
            IF (INDX1.EQ.0) GOTO 306
  302       CONTINUE
              INDX2=IRAM(INDX1)
              IF (INDX2.EQ.0) GOTO 306
              IF (RAM(INDX2+I4).LT.RAM(INDX1+I4)) THEN
                IRAM(INDX1)=IRAM(INDX2)
                IRAM(INDX2)=INDX1
                IRAM(INDX)=INDX2
                GOTO 300
              ENDIF
              INDX=INDX1
              INDX1=INDX2
            GOTO 302
  306       CONTINUE
  310     CONTINUE
  320     CONTINUE
  330     CONTINUE
        ENDIF
  340 CONTINUE
C
C     Writing the results of interpolation:
C     Writing numbers of computed arrivals:
      IF (N1*N2*N3.GT.MRAMP) CALL CIEROR(1)
      NRAMP=0
      I5=0
      DO 130, I3=0,N3-1
      DO 120, I2=0,N2-1
      DO 110, I1=0,N1-1
        NRAMP=NRAMP+1
        I4=0
        INDX=I3*N1*N2+I2*N1+I1
        INDX=MRAM-INDX
  105   CONTINUE
C         Check for the consistentcy of the results
          IF (INDX.LE.MRAMP.OR.INDX.GT.MRAM) THEN
C           MTT-14
            PAUSE 'Error MTT-14: Disorder in array RAM'
            STOP
C           This error should not appear.
C           Please contact the author.
          ENDIF
          INDX=IRAM(INDX)
          IF (INDX.EQ.0) THEN
            IRAM(NRAMP)=I4
            I5=I5+I4
            GOTO 110
          ENDIF
          I4=I4+1
        GOTO 105
  110 CONTINUE
  120 CONTINUE
  130 CONTINUE
      IF (NRAMP.NE.N1*N2*N3) THEN
C       MTT-15
        PAUSE 'Error MTT-15: Disorder in array RAM'
        STOP
C       This error should not appear.
C       Please contact the author.
      ENDIF
      CALL WARRAI(LU5,FILE5,'FORMATTED',.FALSE.,0,.FALSE.,0,N1*N2*N3,
     *                                                    IRAM)
C
C
C     Writing interpolated quantities:
      IF (I5.GT.MRAMP) CALL CIEROR(1)
      DO 240, I4=1,NOUT
        NRAMP=0
        DO 230, I3=0,N3-1
        DO 220, I2=0,N2-1
        DO 210, I1=0,N1-1
          INDX=I3*N1*N2+I2*N1+I1
          INDX=MRAM-INDX
          IF (IRAM(INDX).EQ.0) GOTO 206
          INDX=IRAM(INDX)
  205     CONTINUE
C           Check for the consistentcy of the results
            IF (INDX.LE.MRAMP.OR.INDX.GT.MRAM) THEN
C             MTT-16
              PAUSE 'Error MTT-16: Disorder in array RAM'
              STOP
C             This error should not appear.
C             Please contact the author.
            ENDIF
            NRAMP=NRAMP+1
            IF (OUT(I4)(1:2).EQ.'tt') THEN
C             Travel time:
              RAM(NRAMP)=RAM(INDX+I4)
            ELSEIF (OUT(I4)(1:2).EQ.'hi') THEN
C             Ray history:
              IRAM(NRAMP)=IRAM(INDX+I4)
            ENDIF
            INDX=IRAM(INDX)
          IF (INDX.NE.0) GOTO 205
  206     CONTINUE
  210   CONTINUE
  220   CONTINUE
  230   CONTINUE
        IF (NRAMP.NE.I5) THEN
C         MTT-17
          PAUSE 'Error MTT-17: Disorder in array RAM'
          STOP
C         This error should not appear.
C         Please contact the author.
        ENDIF
        IF (OUT(I4)(1:2).EQ.'tt') THEN
          CALL WARRAY(LU5,FOUT(I4),'FORMATTED',.FALSE.,0.,.FALSE.,0.,I5,
     *                                                    RAM)
        ELSE
          CALL WARRAI(LU5,FOUT(I4),'FORMATTED',.FALSE.,0,.FALSE.,0,I5,
     *                                                    IRAM)
        ENDIF
  240 CONTINUE
C
      STOP
      END
C
C
C=======================================================================
C
      SUBROUTINE CIREAD(LU1,LU2,IT1)
C
C----------------------------------------------------------------------
C Subroutine to read the unformatted output of program CRT and
C to write it into array (I)RAM used in program MTT.
C Reading the output files is completed by a simple invocation of
C subroutine AP00 of file 'ap.for'.
C
      INTEGER LU1,LU2,IT1
C Input:
C     LU1 ... Number of logical unit corresponding to the file with
C             the quantities stored along rays (see C.R.T.5.5.1).
C     LU2 ... Number of logical unit corresponding to the file with
C             the quantities at the initial points of rays,
C             corresponding to file 'RAYPOINTS' (see C.R.T.6.1).
C     IT1 ... Position of the first ray of the actually processed
C             triangle in the array IRAM.
C No output.
C
C Coded by Petr Bulant
C
C ...........................
C Common block /POINTC/ to store the results of complete ray tracing:
      INCLUDE 'pointc.inc'
C None of the storage locations of the common block are altered.
C ...........................
C Common block /RAM/ to store the information about triangles and rays:
      INCLUDE 'mtt.inc'
C.......................................................................
C
C     Auxiliary storage locations:
      INTEGER IRAY0,I1
C     IRAY0.. Index of the last ray read in into the array RAM.
C
C-----------------------------------------------------------------------
C     Loop for the points of rays:
   10 CONTINUE
        IF ((NRAMP+2*NQ).GT.MRAMP) THEN
C         Empting the memory:
          CALL CIERAS
          IF ((NRAMP+2*NQ).GT.MRAMP) CALL CIEROR(1)
        ENDIF
C       Reading the results of the complete ray tracing:
        IRAY0=IRAY
        CALL AP00(0,LU1,LU2)
        IF (IWAVE.LT.1) THEN
C         End of rays:
          IF (IRAY0.LT.IRAM(IT1)) THEN
C           MTT-02
            PAUSE
     *       'Error MTT-02: end of file with rays, wanted ray not found'
            STOP
C           This error may be caused by K2P 
C           not equal to zero, then only two-point rays are stored in
C           output files of CRT.
          ENDIF
          RETURN
        ENDIF
        IF (IRAY.LT.IRAYMI) GOTO 10
        IF (IRAM(IRAY-ORAYE).NE.0) THEN
C         Writing the results of the complete ray tracing - recording
C         new point on a ray:
          IF (IPT.LE.1) THEN
            IF (ISHEET.EQ.0) THEN
C             MTT-18
              PAUSE
     *        'Error MTT-18: a ray with history equal to 0 in CIREAD.'
              STOP
C             All the rays are of the same history, the interpolation
C             in such a case makes no sense. Check the value of the
C             parameter IPOINT in file
C             rpar.for.
            ENDIF
C           New ray - recording the initial point:
C           Coordinates:
            RAM(NRAMP+1)=YI(3)
            RAM(NRAMP+2)=YI(4)
            RAM(NRAMP+3)=YI(5)
C           Index of surface:
            IRAM(NRAMP+4)=0
C           Sequential index of point:
            IRAM(NRAMP+5)=0
C           Quantities to be interpolated:
            DO 20, I1=1,NOUT
              IF (OUT(I1)(1:2).EQ.'tt') THEN
C               Travel time:
                RAM(NRAMP+(5+I1))=YI(1)
              ELSEIF (OUT(I1)(1:2).EQ.'hi') THEN
C               Ray history:
                IRAM(NRAMP+(5+I1))=ISHEET
              ENDIF
  20        CONTINUE
            NRAMP=NRAMP+NQ
          ENDIF
C         Recording the new point:
C         Coordinates:
          RAM(NRAMP+1)=Y(3)
          RAM(NRAMP+2)=Y(4)
          RAM(NRAMP+3)=Y(5)
C         Index of surface:
          IRAM(NRAMP+4)=ISRF
C         Sequential index of point:
          IRAM(NRAMP+5)=IPT
C         Quantities to be interpolated:
          DO 30, I1=1,NOUT
            IF (OUT(I1)(1:2).EQ.'tt') THEN
C             Travel time:
              RAM(NRAMP+(5+I1))=Y(1)
            ELSEIF (OUT(I1)(1:2).EQ.'hi') THEN
C             Ray history:
              IRAM(NRAMP+(5+I1))=ISHEET
            ENDIF
  30      CONTINUE
          NRAMP=NRAMP+NQ
        ENDIF
        IRAM(IRAY-ORAYA)=NRAMP
        IF ((IPT.LE.1).AND.(IRAY.GT.IRAM(IT1)).AND.(IRAY.NE.IRAYMA))
     *     RETURN
      GOTO 10
C
      END
C
C
C=======================================================================
C
      SUBROUTINE CIERAS
C
C----------------------------------------------------------------------
C Subroutine for empting the array (I)RAM. All the parameters
C of all the rays, which will no more be used, are erased.
C
C No input.
C No output.
C
C Coded by Petr Bulant
C
C ...........................
C Common block /POINTC/ to store the results of complete ray tracing:
      INCLUDE 'pointc.inc'
C     IRAY .. Index of the ray being actually read in by CIREAD.
C             This procedure supposes, that any ray with higher
C             index than IRAY was not read in.
C None of the storage locations of the common block are altered.
C ...........................
C Common block /RAM/ to store the information about triangles and rays:
      INCLUDE 'mtt.inc'
C.......................................................................
C     Auxiliary storage locations:
      INTEGER I1,I2,J1
      INTEGER IADDRP
C     I1 ... Controls main loop over rays.
C     I2 ... Controls the loop over parameters of ray I1.
C     J1 ... address of the last used record of array RAM.
C
C-----------------------------------------------------------------------
      J1=IRAM(IRAYMI-1-ORAYA)
      IADDRP=J1
C     Loop for the rays:
      DO 20, I1=IRAYMI,IRAY
        IF (IRAM(I1-ORAYE).GE.(IRAY-1)) THEN
C         This ray is not to be erased:
          DO 10, I2=IADDRP+1,IRAM(I1-ORAYA)
            J1=J1+1
            IRAM(J1)=IRAM(I2)
   10     CONTINUE
          IADDRP=IRAM(I1-ORAYA)
          IRAM(I1-ORAYA)=J1
        ELSE
C         This ray is to be erased:
          IRAM(I1-ORAYE)=0
          IADDRP=IRAM(I1-ORAYA)
          IRAM(I1-ORAYA)=J1
        ENDIF
   20 CONTINUE
      NRAMP=J1
      RETURN
      END
C=======================================================================
C
      SUBROUTINE CITUBE(IT1)
C
C----------------------------------------------------------------------
C Subroutine for interpolation within ray tube formed by the rays
C IRAM(IT1), IRAM(IT1+1), IRAM(IT1+2).
C
      INTEGER IT1
C Input:
C     IT1 ... The address of the index of the first ray of the ray tube,
C             in which the interpolation is to be performed.
C No output.
C
C Coded by Petr Bulant
C
C ...........................
C Common block /RAM/ to store the information about triangles and rays:
      INCLUDE 'mtt.inc'
C.......................................................................
      INTEGER I1A,I2A,I3A,I1B,I2B,I3B
      INTEGER I1MA,I2MA,I3MA,J1,J2
C     I1A,I2A,I3A,I1B,I2B,I3B ... Integers controling main loop over
C       points on the rays (along ray tube). Numbers 1,2,3 denote
C       the first, second and third ray, character A denotes bottom of
C       the ray cell and B denotes top of the ray cell.
C       Each integer contains the address just before the parameters
C       of the corresponding point:
C         the first parameter of the first point: RAM(I1A+1)
C     J1 ... Counts the number of identical points of the ray cell
C            ( J1=0,1,2,3 ).
C     J2 ... Displaies maximum sequential index of a point allowed
C            for actual ray cell.
C
C-----------------------------------------------------------------------
      IF ((IRAM(IRAM(IT1  )-ORAYA).EQ.0).OR.
     *    (IRAM(IRAM(IT1+1)-ORAYA).EQ.0).OR.
     *    (IRAM(IRAM(IT1+2)-ORAYA).EQ.0)) THEN
C       MTT-03
        PAUSE
     * 'Error MTT-03: parameters of a ray not found in memory in CITUBE'
        STOP
C       This error may be caused by K2P 
C       not equal to zero, then only two-point rays are stored in
C       output files of CRT.
      ENDIF
C
      I1MA=IRAM(IRAM(IT1  )-ORAYA)-NQ
      I2MA=IRAM(IRAM(IT1+1)-ORAYA)-NQ
      I3MA=IRAM(IRAM(IT1+2)-ORAYA)-NQ
C     Skipping the first ray cell:
      I1A=IRAM(IRAM(IT1  )-ORAYA-1)+NQ
      I2A=IRAM(IRAM(IT1+1)-ORAYA-1)+NQ
      I3A=IRAM(IRAM(IT1+2)-ORAYA-1)+NQ
C     Loop for points on the rays (loop for ray cells):
  10  CONTINUE
        J1=3
        J2=MAX0(IRAM(I1A+5),IRAM(I2A+5),IRAM(I3A+5))
        IF ((IRAM(I1A+5).EQ.IRAM(I2A+5)).AND.
     *      (IRAM(I1A+5).EQ.IRAM(I3A+5)))  J2=J2+1
        IF ((IRAM(I1A+4).EQ.0).AND.(IRAM(I1A+5).LT.J2)) THEN
          I1B=I1A+NQ
          J1=J1-1
        ELSE
          I1B=I1A
        ENDIF
        IF ((IRAM(I2A+4).EQ.0).AND.(IRAM(I2A+5).LT.J2)) THEN
          I2B=I2A+NQ
          J1=J1-1
        ELSE
          I2B=I2A
        ENDIF
        IF ((IRAM(I3A+4).EQ.0).AND.(IRAM(I3A+5).LT.J2)) THEN
          I3B=I3A+NQ
          J1=J1-1
        ELSE
          I3B=I3A
        ENDIF
        IF (J1.EQ.3) THEN
          IF ((IRAM(I1A+4).EQ.IRAM(I2A+4)).AND.
     *        (IRAM(I1A+4).EQ.IRAM(I3A+4))) THEN
C           Crossing the interface:
            I1A=I1A+NQ
            I2A=I2A+NQ
            I3A=I3A+NQ
            J1=0
            IF ((I1A.GT.I1MA).OR.(I2A.GT.I2MA).OR.(I3A.GT.I3MA))
     *        RETURN
            IRAM(I1A+4)=0
            IRAM(I2A+4)=0
            IRAM(I3A+4)=0
            GOTO 10
          ELSE
            IF (IRAM(I1A+4).EQ.0) THEN
              I1B=I1A+NQ
              J1=J1-1
            ELSE
              I1B=I1A
            ENDIF
            IF (IRAM(I2A+4).EQ.0) THEN
              I2B=I2A+NQ
              J1=J1-1
            ELSE
              I2B=I2A
            ENDIF
            IF (IRAM(I3A+4).EQ.0) THEN
              I3B=I3A+NQ
              J1=J1-1
            ELSE
              I3B=I3A
            ENDIF
            IF (J1.EQ.3) THEN
              IF ((I1A.GE.I1MA).AND.(I2A.GE.I2MA).AND.(I3A.GE.I3MA))
     *          RETURN
C             MTT-04
              PAUSE
     *          'Error MTT-04: the points reached different interfaces.'
              STOP
C             This error should not appear.
C             Please contact the author.
            ENDIF
          ENDIF
        ENDIF
        IF ((I1B.GT.I1MA).OR.(I2B.GT.I2MA).OR.(I3B.GT.I3MA)) THEN
C         MTT-05
          PAUSE 'Error MTT-05: point exceded the ray.'
          STOP
C         This error should not appear.
C         Please contact the author.
        ENDIF
C
C       The ray cell formed by points, whose parameters can be found
C       just behind the addresses I1A,I2A,I3A,I1B,I2B,I3B,
C       is prepared for interpolation:
        IF (((RAM(I1A+1).NE.RAM(I2A+1)).OR.(RAM(I1A+2).NE.RAM(I2A+2))
     *       .OR.(RAM(I1A+3).NE.RAM(I2A+3))) .AND.
     *      ((RAM(I1A+1).NE.RAM(I3A+1)).OR.(RAM(I1A+2).NE.RAM(I3A+2))
     *       .OR.(RAM(I1A+3).NE.RAM(I3A+3))) .AND.
     *      ((RAM(I2A+1).NE.RAM(I3A+1)).OR.(RAM(I2A+2).NE.RAM(I3A+2))
     *       .OR.(RAM(I2A+3).NE.RAM(I3A+3))) .AND.
     *      ((RAM(I1B+1).NE.RAM(I2B+1)).OR.(RAM(I1B+2).NE.RAM(I2B+2))
     *       .OR.(RAM(I1B+3).NE.RAM(I2B+3))) .AND.
     *      ((RAM(I1B+1).NE.RAM(I3B+1)).OR.(RAM(I1B+2).NE.RAM(I3B+2))
     *       .OR.(RAM(I1B+3).NE.RAM(I3B+3))) .AND.
     *      ((RAM(I2B+1).NE.RAM(I3B+1)).OR.(RAM(I2B+2).NE.RAM(I3B+2))
     *       .OR.(RAM(I2B+3).NE.RAM(I3B+3)))) THEN
C         Standard ray cell formed by 4, 5 or 6 points. The bottom
C         and the top of the ray cell are triangles (i.e. they are
C         formed by different points), the three sides of the ray cell
C         are formed by lines, triangles or tetragons.
          CALL CIINTP(I1A,I2A,I3A,I1B,I2B,I3B)
        ELSE
C TADY TO JESTE NENI HOTOVE !!!!!!!
          CONTINUE
        ENDIF
        I1A=I1B
        I2A=I2B
        I3A=I3B
      GOTO 10
C     End of the loop along the ray tube.
      END
C=======================================================================
C
      SUBROUTINE CIINTP(I1A,I2A,I3A,I1B,I2B,I3B)
C
C----------------------------------------------------------------------
C Subroutine for interpolation within standard ray cell formed by the
C points I1A,I2A,I3A,I1B,I2B,I3B.
C
      INTEGER I1A,I2A,I3A,I1B,I2B,I3B
C Input:
C     I1A,I2A,I3A,I1B,I2B,I3B ... Integers specifying the six
C       points on the rays, whch define the ray cell.
C       Numbers 1,2,3 denote the first, second and third ray,
C       character A denotes bottom of the ray cell and B
C       denotes top of the ray cell.
C       Each integer contains the address just before the parameters
C       of the corresponding point:
C         the first parameter of the first point: RAM(I1A+1)
C No output.
C
C Subroutines and external functions required:
      EXTERNAL CIPPP,CILSIG,CICUB,CIKVA,CICUBR,CIKVAR,CIDET3,CISUBD
      REAL CIPPP,CICUB,CIKVA,CIDET3,CISUBD
      LOGICAL CILSIG
C
C Coded by Petr Bulant
C
C ...........................
C Common block /RAM/ to store the information about triangles and rays:
      INCLUDE 'mtt.inc'
C ...........................
C Common block /CIGRD/ to store the parameters of the
C interpolation grid.
      REAL O1,O2,O3,D1,D2,D3
      INTEGER N1,N2,N3
      COMMON/CIGRD/O1,O2,O3,D1,D2,D3,N1,N2,N3
      SAVE /CIGRD/
C     O1,O2,O3 ... Coordinates of the origin of the grid.
C     D1,D2,D3 ... Grid steps.
C     N1,N2,N3 ... Numbers of gridpoints.
C.......................................................................
      REAL ZERO
      PARAMETER (ZERO=0.000001)
      INTEGER K1MI,K1MA,K2MI,K2MA,K3MI,K3MA
      INTEGER I1,I2,I3,I4,I5
      REAL X1,X2,X3
      REAL Y1,Y2,Y3
      REAL XX(3,3),ZZ(3,3)
      REAL AAA,BBB,CCC,DDD,EEE,FFF
      INTEGER NRR
      REAL RR(4),CRR(4),ROOT(3)
      INTEGER IR
      REAL R1,R2
      INTEGER INDX
      REAL ERRMAX
      REAL WB,WC,W1,W2,W3
      REAL A11,A21,A31,A12,A22,A32,A13,A23,A33,B11,B21,B31,B12,B22,B32
      REAL A1,A2,UU11,UU21,UU12,UU22,VV1,VV2,DETUU,TT
C-----------------------------------------------------------------------
C     Choosing gridpoints, which may be contained in the ray cell:
      K1MI=MAX0(0,NINT((
     *     ( AMIN1(RAM(I1A+1),RAM(I2A+1),RAM(I3A+1),
     *             RAM(I1B+1),RAM(I2B+1),RAM(I3B+1)) - O1 ) / D1)+0.5) )
      K1MA=MIN0(N1-1,NINT((
     *     ( AMAX1(RAM(I1A+1),RAM(I2A+1),RAM(I3A+1),
     *             RAM(I1B+1),RAM(I2B+1),RAM(I3B+1)) - O1 ) / D1)-0.5) )
      IF (K1MA.LT.K1MI) RETURN
      K2MI=MAX0(0,NINT((
     *     ( AMIN1(RAM(I1A+2),RAM(I2A+2),RAM(I3A+2),
     *             RAM(I1B+2),RAM(I2B+2),RAM(I3B+2)) - O2 ) / D2)+0.5) )
      K2MA=MIN0(N2-1,NINT((
     *     ( AMAX1(RAM(I1A+2),RAM(I2A+2),RAM(I3A+2),
     *             RAM(I1B+2),RAM(I2B+2),RAM(I3B+2)) - O2 ) / D2)-0.5) )
      IF (K2MA.LT.K2MI) RETURN
      K3MI=MAX0(0,NINT((
     *     ( AMIN1(RAM(I1A+3),RAM(I2A+3),RAM(I3A+3),
     *             RAM(I1B+3),RAM(I2B+3),RAM(I3B+3)) - O3 ) / D3)+0.5) )
      K3MA=MIN0(N3-1,NINT((
     *     ( AMAX1(RAM(I1A+3),RAM(I2A+3),RAM(I3A+3),
     *             RAM(I1B+3),RAM(I2B+3),RAM(I3B+3)) - O3 ) / D3)-0.5) )
      IF (K3MA.LT.K3MI) RETURN
C
C     Checking whether all the rays of the top (bottom) of the ray cell
C     are on the same side of the bottom (top) of the ray cell:
      IF (I1A.NE.I1B.AND.I2A.NE.I2B.AND.I3A.NE.I3B) THEN
        IF (.NOT.CILSIG(
     *    CIPPP(I1A,I2A,I1A,I3A,I1A,RAM(I1B+1),RAM(I1B+2),RAM(I1B+3)),
     *    CIPPP(I1A,I2A,I1A,I3A,I1A,RAM(I2B+1),RAM(I2B+2),RAM(I2B+3)),
     *    CIPPP(I1A,I2A,I1A,I3A,I1A,RAM(I3B+1),RAM(I3B+2),RAM(I3B+3)),
     *          0.)) RETURN
        IF (.NOT.CILSIG(
     *    CIPPP(I1B,I2B,I1B,I3B,I1B,RAM(I1A+1),RAM(I1A+2),RAM(I1A+3)),
     *    CIPPP(I1B,I2B,I1B,I3B,I1B,RAM(I2A+1),RAM(I2A+2),RAM(I2A+3)),
     *    CIPPP(I1B,I2B,I1B,I3B,I1B,RAM(I3A+1),RAM(I3A+2),RAM(I3A+3)),
     *          0.)) RETURN
      ENDIF
C
        ERRMAX=(RAM(I1A+1)**2+RAM(I2A+1)**2+RAM(I3A+1)**2
     *         +RAM(I1B+1)**2+RAM(I2B+1)**2+RAM(I3B+1)**2
     *         +RAM(I1A+2)**2+RAM(I2A+2)**2+RAM(I3A+2)**2
     *         +RAM(I1B+2)**2+RAM(I2B+2)**2+RAM(I3B+2)**2
     *         +RAM(I1A+3)**2+RAM(I2A+3)**2+RAM(I3A+3)**2
     *         +RAM(I1B+3)**2+RAM(I2B+3)**2+RAM(I3B+3)**2)/6.0E12
C
C     Loop for all the selected gridpoints:
      DO 100, I1=K1MI,K1MA
      DO  90, I2=K2MI,K2MA
      DO  80, I3=K3MI,K3MA
        X1=O1+D1*FLOAT(I1)
        X2=O2+D2*FLOAT(I2)
        X3=O3+D3*FLOAT(I3)
        INDX=I3*N1*N2+I2*N1+I1+1
C
C
C       Checking the position of the gridpoint with respect to
C       the planes bounding the ray cell:
C
        IF (I1A.NE.I1B.AND.I2A.NE.I2B.AND.I3A.NE.I3B) THEN
          IF (.NOT.CILSIG(
     *      CIPPP(I1B,I2B,I1B,I3B,I1B,RAM(I1A+1),RAM(I1A+2),RAM(I1A+3)),
     *      CIPPP(I1B,I2B,I1B,I3B,I1B,X1,X2,X3) , 0. , 0.             ))
     *      GOTO 79
C
          IF (.NOT.CILSIG(
     *      CIPPP(I1A,I2A,I1A,I3A,I1A,RAM(I1B+1),RAM(I1B+2),RAM(I1B+3)),
     *      CIPPP(I1A,I2A,I1A,I3A,I1A,X1,X2,X3) , 0. , 0.             ))
     *      GOTO 79
C
          IF (CILSIG(
     *    CIPPP(I1A,I2B,I2A,I1B,I1A,RAM(I2A+1),RAM(I2A+2),RAM(I2A+3)),
     *    CIPPP(I1A,I2B,I2A,I1B,I1A,RAM(I3A+1),RAM(I3A+2),RAM(I3A+3)),
     *    CIPPP(I1A,I2B,I2A,I1B,I1A,RAM(I3B+1),RAM(I3B+2),RAM(I3B+3)),
     *          0.)) THEN
            IF (.NOT.CILSIG(
     *      CIPPP(I1A,I2B,I2A,I1B,I1A,RAM(I2A+1),RAM(I2A+2),RAM(I2A+3)),
     *      CIPPP(I1A,I2B,I2A,I1B,I1A,X1,X2,X3),0.,0.))          GOTO 79
          ELSEIF (CILSIG(
     *    CIPPP(I1A,I2B,I2A,I1B,I2A,RAM(I1A+1),RAM(I1A+2),RAM(I1A+3)),
     *    CIPPP(I1A,I2B,I2A,I1B,I2A,RAM(I3A+1),RAM(I3A+2),RAM(I3A+3)),
     *    CIPPP(I1A,I2B,I2A,I1B,I2A,RAM(I3B+1),RAM(I3B+2),RAM(I3B+3)),
     *          0.)) THEN
            IF (.NOT.CILSIG(
     *      CIPPP(I1A,I2B,I2A,I1B,I2A,RAM(I1A+1),RAM(I1A+2),RAM(I1A+3)),
     *      CIPPP(I1A,I2B,I2A,I1B,I2A,X1,X2,X3),0.,0.))          GOTO 79
          ENDIF
C
          IF (CILSIG(
     *    CIPPP(I2A,I3B,I3A,I2B,I2A,RAM(I3A+1),RAM(I3A+2),RAM(I3A+3)),
     *    CIPPP(I2A,I3B,I3A,I2B,I2A,RAM(I1A+1),RAM(I1A+2),RAM(I1A+3)),
     *    CIPPP(I2A,I3B,I3A,I2B,I2A,RAM(I1B+1),RAM(I1B+2),RAM(I1B+3)),
     *          0.)) THEN
            IF (.NOT.CILSIG(
     *      CIPPP(I2A,I3B,I3A,I2B,I2A,RAM(I3A+1),RAM(I3A+2),RAM(I3A+3)),
     *      CIPPP(I2A,I3B,I3A,I2B,I2A,X1,X2,X3),0.,0.))          GOTO 79
          ELSEIF (CILSIG(
     *    CIPPP(I2A,I3B,I3A,I2B,I3A,RAM(I2A+1),RAM(I2A+2),RAM(I2A+3)),
     *    CIPPP(I2A,I3B,I3A,I2B,I3A,RAM(I1A+1),RAM(I1A+2),RAM(I1A+3)),
     *    CIPPP(I2A,I3B,I3A,I2B,I3A,RAM(I1B+1),RAM(I1B+2),RAM(I1B+3)),
     *          0.)) THEN
            IF (.NOT.CILSIG(
     *      CIPPP(I2A,I3B,I3A,I2B,I3A,RAM(I2A+1),RAM(I2A+2),RAM(I2A+3)),
     *      CIPPP(I2A,I3B,I3A,I2B,I3A,X1,X2,X3),0.,0.))          GOTO 79
          ENDIF
C
          IF (CILSIG(
     *    CIPPP(I3A,I1B,I1A,I3B,I1A,RAM(I3A+1),RAM(I3A+2),RAM(I3A+3)),
     *    CIPPP(I3A,I1B,I1A,I3B,I1A,RAM(I2A+1),RAM(I2A+2),RAM(I2A+3)),
     *    CIPPP(I3A,I1B,I1A,I3B,I1A,RAM(I2B+1),RAM(I2B+2),RAM(I2B+3)),
     *          0.)) THEN
            IF (.NOT.CILSIG(
     *      CIPPP(I3A,I1B,I1A,I3B,I1A,RAM(I3A+1),RAM(I3A+2),RAM(I3A+3)),
     *      CIPPP(I3A,I1B,I1A,I3B,I1A,X1,X2,X3),0.,0.))          GOTO 79
          ELSEIF (CILSIG(
     *    CIPPP(I3A,I1B,I1A,I3B,I3A,RAM(I1A+1),RAM(I1A+2),RAM(I1A+3)),
     *    CIPPP(I3A,I1B,I1A,I3B,I3A,RAM(I2A+1),RAM(I2A+2),RAM(I2A+3)),
     *    CIPPP(I3A,I1B,I1A,I3B,I3A,RAM(I2B+1),RAM(I2B+2),RAM(I2B+3)),
     *          0.)) THEN
            IF (.NOT.CILSIG(
     *      CIPPP(I3A,I1B,I1A,I3B,I3A,RAM(I1A+1),RAM(I1A+2),RAM(I1A+3)),
     *      CIPPP(I3A,I1B,I1A,I3B,I3A,X1,X2,X3),0.,0.))          GOTO 79
          ENDIF
        ENDIF
C
C
C       Computation of the coefficients AAA,BBB,CCC,DDD
C       of the cubic equation:
C
        XX(1,1)=RAM(I1A+1)-RAM(I1B+1)
        XX(1,2)=RAM(I1A+2)-RAM(I1B+2)
        XX(1,3)=RAM(I1A+3)-RAM(I1B+3)
        XX(2,1)=RAM(I2A+1)-RAM(I2B+1)
        XX(2,2)=RAM(I2A+2)-RAM(I2B+2)
        XX(2,3)=RAM(I2A+3)-RAM(I2B+3)
        XX(3,1)=RAM(I3A+1)-RAM(I3B+1)
        XX(3,2)=RAM(I3A+2)-RAM(I3B+2)
        XX(3,3)=RAM(I3A+3)-RAM(I3B+3)
C
        ZZ(1,1)=RAM(I1B+1)-X1
        ZZ(1,2)=RAM(I1B+2)-X2
        ZZ(1,3)=RAM(I1B+3)-X3
        ZZ(2,1)=RAM(I2B+1)-X1
        ZZ(2,2)=RAM(I2B+2)-X2
        ZZ(2,3)=RAM(I2B+3)-X3
        ZZ(3,1)=RAM(I3B+1)-X1
        ZZ(3,2)=RAM(I3B+2)-X2
        ZZ(3,3)=RAM(I3B+3)-X3
C
        AAA=CIDET3(XX)
        BBB=CISUBD(ZZ,XX)
        CCC=CISUBD(XX,ZZ)
        DDD=CIDET3(ZZ)
C
C
C       Solving the equation according to their coefficients:
        IF (AAA.EQ.0.) THEN
          IF (BBB.EQ.0.) THEN
C           Linear equation:
            IF (CCC.EQ.0.) THEN
C             MTT-06
              PAUSE
     *          'Error MTT-06: wrong coefficients of equation in CIINTP'
              STOP
C             This error should not appear.
C             Please contact the author.
            ENDIF
            R1=-DDD/CCC
            IF (R1.GE.0..AND.R1.LE.1.) THEN
              IR=1
              ROOT(1)=R1
            ELSE
              IR=0
            ENDIF
          ELSE
C           Kvadratic equation:
            CALL CIKVAR(BBB,CCC,DDD,IR,R1,R2)
            IF (IR.GE.1) THEN
              IF ((R1.LT.0.).OR.(R1.GT.1.)) THEN
C               MTT-11
                PAUSE 'Error MTT-11: First root of quadratic eq. wrong'
                STOP
C               This error should not appear.
C               Please contact the author.
              ENDIF
              ROOT(1)=R1
            ENDIF
            IF (IR.GE.2) THEN
              IF ((R2.LT.0.).OR.(R2.GT.1.)) THEN
C               MTT-12
                PAUSE 'Error MTT-12: Second root of quadratic eq. wrong'
                STOP
C               This error should not appear.
C               Please contact the author.
              ENDIF
              ROOT(2)=R2
            ENDIF
          ENDIF
        ELSE
C         Cubic equation:
C         Computation of the coefficients EEE,FFF,CCC of the derivative
C         of the cubic equation:
          EEE=3.*AAA
          FFF=2.*BBB
C
C         Solving the equation for the derivative:
          CALL CIKVAR(EEE,FFF,CCC,IR,R1,R2)
C         Deciding, whether the cubic equation is to be solved:
          RR(1)=0.
          IF (IR.GE.1) THEN
            RR(2)=R1
            NRR=3
          ELSE
            NRR=2
          ENDIF
          IF (IR.GE.2) THEN
            RR(NRR)=R2
            NRR=NRR+1
          ENDIF
          RR(NRR)=1.
          DO 10, I4=1,NRR
            CRR(I4)=CICUB(AAA,BBB,CCC,DDD,RR(I4))
  10      CONTINUE
          IR=0
          I4=1
  20      CONTINUE
            IF ((CRR(I4)*CRR(I4+1)).LE.0.) THEN
              IR=IR+1
              I4=I4+1
            ELSE
              DO 30, I5=I4,NRR-1
                CRR(I5)=CRR(I5+1)
                RR(I5)=RR(I5+1)
  30          CONTINUE
              NRR=NRR-1
            ENDIF
          IF (I4.LT.NRR) GOTO 20
          IF (IR.EQ.0) GOTO 79
C
C
C         Solving the cubic equation:
          CALL CICUBR(AAA,BBB,CCC,DDD,RR,IR,ROOT)
        ENDIF
C
        DO 70, I4=1,IR
          WB=ROOT(I4)
          IF (WB.EQ.0.) GOTO 79
          WC=1.-WB
          IF ((WB.LT.0.).OR.(WC.LT.0.)) THEN
C           MTT-13
            PAUSE 'Error MTT-13: Root outside the cell'
            STOP
C           This error should not appear.
C           Please contact the author.
          ENDIF
          A11=WB*RAM(I1A+1) + WC*RAM(I1B+1)
          A21=WB*RAM(I1A+2) + WC*RAM(I1B+2)
          A31=WB*RAM(I1A+3) + WC*RAM(I1B+3)
          A12=WB*RAM(I2A+1) + WC*RAM(I2B+1)
          A22=WB*RAM(I2A+2) + WC*RAM(I2B+2)
          A32=WB*RAM(I2A+3) + WC*RAM(I2B+3)
          A13=WB*RAM(I3A+1) + WC*RAM(I3B+1)
          A23=WB*RAM(I3A+2) + WC*RAM(I3B+2)
          A33=WB*RAM(I3A+3) + WC*RAM(I3B+3)
          A11=A11-A13
          A21=A21-A23
          A31=A31-A33
          A12=A12-A13
          A22=A22-A23
          A32=A32-A33
          A13=X1-A13
          A23=X2-A23
          A33=X3-A33
          A1=SQRT(A11*A11+A21*A21+A31*A31)
          A2=SQRT(A12*A12+A22*A22+A32*A32)
          B11=A11*A2+A12*A1
          B21=A21*A2+A22*A1
          B31=A31*A2+A32*A1
          B12=A11*A2-A12*A1
          B22=A21*A2-A22*A1
          B32=A31*A2-A32*A1
          UU11=B11*A11+B21*A21+B31*A31
          UU21=B12*A11+B22*A21+B32*A31
          UU12=B11*A12+B21*A22+B31*A32
          UU22=B12*A12+B22*A22+B32*A32
          VV1 =B11*A13+B21*A23+B31*A33
          VV2 =B12*A13+B22*A23+B32*A33
          DETUU=UU11*UU22-UU12*UU21
C         Determinant eq zero in case of infinitely thin ray cell,
C         in such a case gridpoint cannot lie inside the cell.
          IF (DETUU.NE.0.) THEN
            W1=(UU22*VV1-UU12*VV2)/DETUU
            W2=(UU11*VV2-UU21*VV1)/DETUU
            W3=1.-W1-W2
            IF ((W1.GE.0.).AND.(W1.LE.1.).AND.
     *          (W2.GE.0.).AND.(W2.LE.1.).AND.
     *          (W3.GE.0.).AND.(W3.LE.1.)) THEN
C             The gridpoint is situated inside the ray cell.
C
C             Checking accuracy of interpolation coefficients:
              Y1=WB*(W1*RAM(I1A+1)+W2*RAM(I2A+1)+W3*RAM(I3A+1))+
     *           WC*(W1*RAM(I1B+1)+W2*RAM(I2B+1)+W3*RAM(I3B+1))
              Y2=WB*(W1*RAM(I1A+2)+W2*RAM(I2A+2)+W3*RAM(I3A+2))+
     *           WC*(W1*RAM(I1B+2)+W2*RAM(I2B+2)+W3*RAM(I3B+2))
              Y3=WB*(W1*RAM(I1A+3)+W2*RAM(I2A+3)+W3*RAM(I3A+3))+
     *           WC*(W1*RAM(I1B+3)+W2*RAM(I2B+3)+W3*RAM(I3B+3))
              IF (((Y1-X1)**2 + (Y2-X2)**2 + (Y3-X3)**2).GT.ERRMAX) THEN
C               MTT-23
                PAUSE 'Error MTT-23: Inexact interpolation coefficients'
                STOP
C               Interpolation coefficients are not exact enough to
C               provide exact coordinates of gridpoints. Thus the error
C               of the same order will appear in interpolation of
C               travel times and other quantities.
C               In principle this error should not appear.
C               On the other hand, if you wish to continue in
C               interpolation even with inaccurate interpolation
C               coefficients, put the letter 'C' to the beginning of the
C               previous line with 'STOP'.
              ENDIF
C
C             Interpolation of output quantities:
              INDX=MRAM-INDX+1
  40          CONTINUE
              IF (IRAM(INDX).NE.0) THEN
                INDX=IRAM(INDX)
                GOTO 40
              ENDIF
              NRAMT=NRAMT-(NOUT+1)
              IF (NRAMT.LE.NRAMP) CALL CIEROR(1)
              IRAM(INDX)=NRAMT
              IRAM(NRAMT)=0
              MRAMP=NRAMT-1
              DO 50, I5=1,NOUT
                IF (OUT(I5)(1:2).EQ.'tt') THEN
C                 Travel time:
                  TT=WB*W1*RAM(I1A+(5+I5))+WB*W2*RAM(I2A+(5+I5))+
     *               WB*W3*RAM(I3A+(5+I5))+
     *               WC*W1*RAM(I1B+(5+I5))+WC*W2*RAM(I2B+(5+I5))+
     *               WC*W3*RAM(I3B+(5+I5))
                  RAM(NRAMT+I5)=TT
                ELSEIF (OUT(I5)(1:2).EQ.'hi') THEN
C                 Ray history:
                  IRAM(NRAMT+I5)=IRAM(I1A+(5+I5))
                ENDIF
  50          CONTINUE
              GOTO 79
            ENDIF
          ENDIF
  70    CONTINUE
  79    CONTINUE
  80  CONTINUE
  90  CONTINUE
 100  CONTINUE
      END
C
C=======================================================================
C
      REAL FUNCTION CIPPP(IU1,IU2,IV1,IV2,IA,X1,X2,X3)
C
C----------------------------------------------------------------------
      INTEGER IU1,IU2,IV1,IV2,IA
      REAL X1,X2,X3
C     Subroutine designed to indicate the position of
C     a point X: (X1,X2,X3) with respect to
C     a plane defined by two vectors u: IU1 --> IU2, v: IV1 --> IV2,
C                 and by a point IA.
C     The subroutine computes the scalar product a.w of vector
C     a: IA --> X with vector w: w = u x v being normal to the plane.
C     The sign of this product then indicates, on which side of the
C     plane is the point X located.
C     All the computation is performed in 3-D Cartesian space.
C
C Input:
C     IU1,IU2 .. Two points defining vector u.
C     IV1,IV2 .. Two points defining vector v.
C     IA ... Point of the plane.
C     Each integer contains the address just before the parameters
C     of the corresponding point:
C         the first parameter of the point IA: RAM(IA+1)
C     X1,X2,X3 ... Coordinates of the examined point.
C Output:
C     CIPPP ... The value of the scalar product.
C
C Coded by Petr Bulant
C
      INCLUDE 'mtt.inc'
C
C......................................................................
      REAL U1,U2,U3,V1,V2,V3,W1,W2,W3,A1,A2,A3
C-----------------------------------------------------------------------
      U1=RAM(IU2+1)-RAM(IU1+1)
      U2=RAM(IU2+2)-RAM(IU1+2)
      U3=RAM(IU2+3)-RAM(IU1+3)
      V1=RAM(IV2+1)-RAM(IV1+1)
      V2=RAM(IV2+2)-RAM(IV1+2)
      V3=RAM(IV2+3)-RAM(IV1+3)
      W1=U2*V3-U3*V2
      W2=U3*V1-U1*V3
      W3=U1*V2-U2*V1
      A1=X1-RAM(IA+1)
      A2=X2-RAM(IA+2)
      A3=X3-RAM(IA+3)
      CIPPP=W1*A1+W2*A2+W3*A3
      END
C
C=======================================================================
C
      REAL FUNCTION CICUB(AAA,BBB,CCC,DDD,XX)
C
C----------------------------------------------------------------------
      REAL AAA,BBB,CCC,DDD,XX
C     Subroutine designed to calculate the value of cubic function
C     given by coefficients AAA,BBB,CCC,DDD in point XX.
C......................................................................
      CICUB=AAA*XX*XX*XX+BBB*XX*XX+CCC*XX+DDD
      END
C
C=======================================================================
C
      REAL FUNCTION CIKVA(AAA,BBB,CCC,XX)
C
C----------------------------------------------------------------------
      REAL AAA,BBB,CCC,XX
C     Subroutine designed to calculate the value of kvadratic function
C     given by coefficients AAA,BBB,CCC in point XX.
C......................................................................
      CIKVA=AAA*XX*XX+BBB*XX+CCC
      END
C
C=======================================================================
C
      REAL FUNCTION CIDET3(A)
C
C----------------------------------------------------------------------
      REAL A(3,3)
C     Subroutine designed to calculate the determinant of 3x3 matrix AA.
C Input:
C     A  ... 3x3 real matrix.
C Output:
C     CIDET3 . Determinant.
C
C Coded by Petr Bulant
C......................................................................
      CIDET3=A(1,1)*(A(2,2)*A(3,3)-A(3,2)*A(2,3))+
     *       A(2,1)*(A(3,2)*A(1,3)-A(1,2)*A(3,3))+
     *       A(3,1)*(A(1,2)*A(2,3)-A(2,2)*A(1,3))
      END
C
C=======================================================================
C
      REAL FUNCTION CISUBD(A,B)
C
C----------------------------------------------------------------------
      REAL A(3,3),B(3,3)
C     Subroutine designed to calculate the sum
C
C           SUM  (  e  *  A * B * B   )
C           ijk      ijk   1i  2j  3k
C
C     where i,j,k=1..3 and e    is a Levi-Civita symbol.
C                           ijk
C Input:
C     A  ... 3x3 real matrix.
C     B  ... 3x3 real matrix.
C Output:
C     CISUBD . Value of the sum mentioned above.
C
C Coded by Petr Bulant
C......................................................................
      CISUBD=A(1,1)*(B(2,2)*B(3,3)-B(3,2)*B(2,3))+
     *       A(2,1)*(B(3,2)*B(1,3)-B(1,2)*B(3,3))+
     *       A(3,1)*(B(1,2)*B(2,3)-B(2,2)*B(1,3))+
     *       A(1,2)*(B(2,3)*B(3,1)-B(3,3)*B(2,1))+
     *       A(2,2)*(B(3,3)*B(1,1)-B(1,3)*B(3,1))+
     *       A(3,2)*(B(1,3)*B(2,1)-B(2,3)*B(1,1))+
     *       A(1,3)*(B(2,1)*B(3,2)-B(3,1)*B(2,2))+
     *       A(2,3)*(B(3,1)*B(1,2)-B(1,1)*B(3,2))+
     *       A(3,3)*(B(1,1)*B(2,2)-B(2,1)*B(1,2))
      END
C
C
C=======================================================================
C
      SUBROUTINE CIKVAR(AAA,BBB,CCC,IR,R1,R2)
C
C----------------------------------------------------------------------
C Subroutine to solve the kvadratic equation on the interval <0,1>.
      INTEGER IR
      REAL AAA,BBB,CCC,R1,R2
C Input:
C     AAA,BBB,CCC ... Coefficients of the equation.
C Output:
C     IR ... Number of real roots of the equation from <0,1>.
C     R1,R2. Real roots of the equation from <0,1>, R1.
      INTEGER IR
      REAL AA0,BB0,CC0,DD0,RR(4),ROOT(3)
C Input:
C     AA0,BB0,CC0,DD0 ... Coefficients of the equation.
C     IR ... Expected number of real roots from the interval <0,1>.
C     RR ... Array of points, which divide the interval <0,1> into
C            subintervals, each subinterval containing just one
C            root of the cubic equation.
C         RR(1)   ... 0.
C         RR(NRR) ... 1.
C         for IR=1 NRR=2
C         for IR=2 NRR=3, RR(2) is a point in which
C                               the cubic equation has zero derivative
C         for IR=3 NRR=4, RR(2) and RR(3) are points in which
C                               the cubic equation has zero derivative.
C Output:
C     ROOT  ... Real roots of the equation from <0,1> or undefined.
C
C Subroutines and external functions required:
      EXTERNAL CICUB,CIKVA
      REAL CICUB,CIKVA
C
C Coded by Petr Bulant
C
C.......................................................................
      REAL ZERO
      PARAMETER (ZERO=0.000001)
      REAL TRIAA0,DVEBB0,POM1,POM2,POM3,DER,RA,RB
      INTEGER I1,I2
C-----------------------------------------------------------------------
      TRIAA0=3.*AA0
      DVEBB0=2.*BB0
      DO 20, I2=1,IR
        RA=RR(I2)
        RB=RR(I2+1)
        I1=0
  10    CONTINUE
          I1=I1+1
          IF (I1.GT.100) THEN
            IF (ABS(RA-RB).LE.ZERO) RETURN
C           MTT-08
            PAUSE 'Error MTT-08: infinite loop in CICUBR'
c           STOP
C           This error should not appear.
C           Please contact the author.
          ENDIF
C
          POM1=CICUB(AA0,BB0,CC0,DD0,RB)
          POM2=CICUB(AA0,BB0,CC0,DD0,RA)
          DER=CIKVA(TRIAA0,DVEBB0,CC0,RA)
          IF (DER.NE.0.) THEN
            ROOT(I2)=RA-POM2/DER
            IF((ROOT(I2).LE.RA).OR.(ROOT(I2).GE.RB))ROOT(I2)=(RA+RB)/2.
          ELSE
            ROOT(I2)=(RA+RB)/2.
          ENDIF
          POM3=CICUB(AA0,BB0,CC0,DD0,ROOT(I2))
          IF (POM1*POM3.LT.0.) THEN
            RA=ROOT(I2)
          ELSE
            RB=ROOT(I2)
          ENDIF
C
          POM1=CICUB(AA0,BB0,CC0,DD0,RB)
          POM2=CICUB(AA0,BB0,CC0,DD0,RA)
          POM3=POM1-POM2
          IF (POM3.EQ.0.) THEN
C           MTT-09
            PAUSE 'Error MTT-09: division by zero in CICUBR.'
            STOP
C           This error should not appear.
C           Please contact the author.
          ENDIF
          ROOT(I2)=RA-(RB-RA)*POM2/POM3
          POM3=CICUB(AA0,BB0,CC0,DD0,ROOT(I2))
          IF (POM1*POM3.LT.0.) THEN
            RA=ROOT(I2)
            POM2=POM3
          ELSE
            RB=ROOT(I2)
            POM1=POM3
          ENDIF
        IF (ABS(RA-RB).GT.ZERO) GOTO 10
        IF (RA.NE.RB) THEN
          IF (POM1.EQ.0.) THEN
            ROOT(I2)=RB
          ELSE IF (POM2.EQ.0.) THEN
            ROOT(I2)=RA
          ELSE IF (POM1-POM2.NE.0.) THEN
C           Linear interpolation of the root:
            RA=RA-ROOT(I2)
            RB=RB-ROOT(I2)
            POM3=(RA*POM1-RB*POM2)/(POM1-POM2)
            ROOT(I2)=ROOT(I2)+POM3
          ELSE
C           MTT-22
            PAUSE 'Error MTT-22: Equal functional values in CICUBR.'
            STOP
C           This error should not appear.
C           Please contact the author.
            STOP
          END IF
        END IF
  20  CONTINUE
C
      END
C
C=======================================================================
C
      LOGICAL FUNCTION CILSIG(R1,R2,R3,R4)
C
C----------------------------------------------------------------------
      REAL R1,R2,R3,R4
C     Subroutine to compare signs of real quantities. All the quantities
C     with absolute value less then parameter ZERO are considered to
C     be equal 0.
C
C Input: R1,R2,R3,R4 ... Real quantities.
C Output: CILSIG ... .TRUE. for all the quantities nonnegative or
C                           all th quantities nonpositive.
C                    .FALSE. in other cases.
C.......................................................................
      REAL ZERO
      PARAMETER (ZERO=.000001)
C     ZERO ... Constant used to decide whether
C               the real quantiti .EQ. zero.
C-----------------------------------------------------------------------
      IF (ABS(R1).LT.ZERO) R1=0.
      IF (ABS(R2).LT.ZERO) R2=0.
      IF (ABS(R3).LT.ZERO) R3=0.
      IF (ABS(R4).LT.ZERO) R4=0.
      IF (((R1.LE.0.).AND.(R2.LE.0.).AND.(R3.LE.0.).AND.(R4.LE.0.)).OR.
     *    ((R1.GE.0.).AND.(R2.GE.0.).AND.(R3.GE.0.).AND.(R4.GE.0.)))
     *  THEN
        CILSIG=.TRUE.
      ELSE
        CILSIG=.FALSE.
      ENDIF
      END
C
C
C=======================================================================
C
      SUBROUTINE CIGRID(LU4,FILE4)
C
C----------------------------------------------------------------------
C Subroutine to read in the file with parameters of an interpolation
C grid. The file is assumed to be in form of a SEP-like parameter or
C header file. The read quantities are then written into common block
C CIGRD.
      INTEGER LU4
      CHARACTER*80 FILE4
C Input:
C     LU4 ... Number of logical unit corresponding to the file with
C             the input quantities for interpolation grid.
C     FILE4 . Name of the file.
C No output.
C
C Coded by Petr Bulant
C
C ...........................
C Common block /RAM/ to store the information about triangles and rays:
      INCLUDE 'mtt.inc'
C ...........................
C Common block /CIGRD/ to store the parameters of the
C interpolation grid.
      REAL O1,O2,O3,D1,D2,D3
      INTEGER N1,N2,N3
      COMMON/CIGRD/O1,O2,O3,D1,D2,D3,N1,N2,N3
      SAVE /CIGRD/
C     O1,O2,O3 ... Coordinates of the origin of the grid.
C     D1,D2,D3 ... Grid steps.
C     N1,N2,N3 ... Numbers of gridpoints.
C.......................................................................
      INTEGER I1
C-----------------------------------------------------------------------
C     Reading all the data from file FILE4 to the memory
C     (SEP parameter file form):
      CALL RSEP1(LU4,FILE4)
C
C     Recalling the data specifying grid dimensions
C     (arguments: Name of value in input data, Variable, Default):
      CALL RSEP3R('O1',O1,0.)
      CALL RSEP3R('O2',O2,0.)
      CALL RSEP3R('O3',O3,0.)
      CALL RSEP3R('D1',D1,1.)
      CALL RSEP3R('D2',D2,1.)
      CALL RSEP3R('D3',D3,1.)
      CALL RSEP3I('N1',N1,1)
      CALL RSEP3I('N2',N2,1)
      CALL RSEP3I('N3',N3,1)
C
      IF ((D1.LT.0.).OR.(D2.LT.0.).OR.(D3.LT.0.)) GOTO 10
      IF ((N1.LE.0).OR.(N2.LE.0).OR.(N3.LE.0)) GOTO 10
      IF (D1.EQ.0.) THEN
        IF (N1.EQ.1) THEN
          D1=1.
        ELSE
          GOTO 10
        ENDIF
      ENDIF
      IF (D2.EQ.0.) THEN
        IF (N2.EQ.1) THEN
          D2=1.
        ELSE
          GOTO 10
        ENDIF
      ENDIF
      IF (D3.EQ.0.) THEN
        IF (N3.EQ.1) THEN
          D3=1.
        ELSE
          GOTO 10
        ENDIF
      ENDIF
      MRAMP=MRAM-N1*N2*N3
      IF (MRAMP.LE.1) CALL CIEROR(1)
      NRAMT=MRAMP+1
      DO 5, I1=NRAMT,MRAM
        IRAM(I1)=0
   5  CONTINUE
      RETURN
C
  10  CONTINUE
C     MTT-10
      PAUSE 'Error MTT-10: This specification of the interpolation
     *       grid may cause problems. Please, specify
     *       D1,D2,D3 and N1,N2,N3 greater than 0.
     *       Di may equal 0 in case that corresponding Ni equals 1.    '
      STOP
C
      END
C
C
C=======================================================================
C
      SUBROUTINE CIEROR(IERR)
C
C-----------------------------------------------------------------------
      INTEGER IERR
C
C Subroutine designed to print error messages of different
C CI* subroutines using command 'PAUSE'.
C
C Input:
C     IERR ... Index of the error.
C No output.
C Coded by Petr Bulant
C-----------------------------------------------------------------------
C
      IF (IERR.EQ.001) THEN
C        MTT-01
        PAUSE  'Error MTT-01: array (I)RAM too small.'
        STOP
C       This error may be caused by too small dimension of array
C       RAM. Try to enlarge the parameter MRAM in common block
C       RAM.
      ENDIF
      END
C
C=======================================================================
C
      INCLUDE 'ap.for'
C     ap.for
      INCLUDE 'forms.for'
C     forms.for
      INCLUDE 'length.for'
C     length.for
      INCLUDE 'sep.for'
C     sep.for
C
C=======================================================================
C