C                                                            <HTML><PRE>
C=======================================================================
C
      PROGRAM MTT
C to calculate Multivalued Travel Times on a 3-D grid of points.
C
C Version: 5.20
C Date: 1998, November 11
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)-st 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 performs bicubic interpolation of
C travel times within ray cells (in the coordinate system connected
C with each cell), other quantities are interpolated bilinearly.
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 <A HREF="rpar.for#data2">AERR</A> 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 <A HREF="rpar.for#data2">PRM0(3)</A> 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 <A HREF="rpar.for#data2">PRM0(3)</A> 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 Subroutines and external functions required: see the end of the file.
C
C.......................................................................
C
C                                                    <A NAME="data"></A>
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               <A HREF="#SEP">Description of file SEP</A>
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               <A HREF="writ.for#out1">Description of file CRT-R</A>
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               <A HREF="writ.for#outi">Description of file CRT-I</A>
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               <A HREF="writ.for#outt">Description of file CRT-T</A>
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               <A HREF="#NUM">Description of file NUM</A>
C     Default:  'NUM'='num.out'
C (3) Zero to MOUT times following line (for MOUT refer to 'mtt.inc'):
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               <A HREF="#FOUT">Description of files 'FOUT(i)'</A>
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 time.
C               'hi'   .......  Ray history.
C               'p1'   .......  First component of the slowness vector.
C               'p2'   .......  Second component of the slowness vector.
C               'p3'   .......  Third component of the slowness vector.
C     Defaults: 'FOUT(1)'   ='mtt.out',  'OUT(1)'  ='tt',
C               'FOUT(2)'   =' '      ,  'OUT(2)'  =' ' ,
C                   .                      .
C                   .                      .
C               'FOUT(MOUT)'=' '      ,  'OUT(MOUT)=' ' .
C
C
C                                                     <A NAME="SEP"></A>
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 <A HREF="../forms/sep.for#data"> 'sep.for'</A>.
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 <A HREF="len/leni-mtt.h">Example of data SEP</A>
C
C                                                     <A NAME="NUM"></A>
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                                                    <A NAME="FOUT"></A>
C Output formatted file 'FOUT(i)':
C     N numbers, where N is the sum of integers in file NUM.
C     For each gridpoint, all values of a single interpolated quantity
C     determined at the gridpoint.  The values are stored in descending
C     order according to the travel time.
C     The loop over gridpoints is the same as in file 'NUM'.
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                                                  <A NAME="author"></A>
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     Auxiliary storage locations:
      INTEGER IRAY1,IRAY2,IRAY3
      INTEGER I1,I2,I3,I4,I5,INDX,INDX1,INDX2,IT1,IIT1
      INTEGER I
      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
      CHARACTER*80 OUTI(MOUT),FOUTI(MOUT)
C
C-----------------------------------------------------------------------
C
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
      WRITE(*,'(''+'',79('' ''))')
      WRITE(*,'(A)') '+Reading input data.'
      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:
      OUTI(1)='tt'
      OUTI(2)=' '
      OUTI(3)=' '
      OUTI(4)=' '
      OUTI(5)=' '
      OUTI(6)=' '
      FILE5='num.out'
      FOUTI(1)='mtt.out'
      FOUTI(2)=' '
      FOUTI(3)=' '
      FOUTI(4)=' '
      FOUTI(5)=' '
      FOUTI(6)=' '
      READ(LU0,*) FILE5
      I1=1
   2  CONTINUE
        READ(LU0,*,END=4) FOUTI(I1),OUTI(I1)
        I1=I1+1
        IF (I1.GT.MOUT) THEN
C         <A NAME="error19"></A><FONT COLOR="RED">MTT-19</FONT>
          CALL ERROR('MTT-19: Small array FOUT for the output files.')
C         This error may be caused by wrongly prescribed
C         <A HREF="#data"> input data FOUT and OUT</A>, 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:
      DO 5, I1=1,MOUT
        CALL LOWER(OUTI(I1))
        IF ((OUTI(I1)(1:2).NE.' ').AND.(OUTI(I1)(1:2).NE.'tt').AND.
     *      (OUTI(I1)(1:2).NE.'hi').AND.(OUTI(I1)(1:2).NE.'p1').AND.
     *      (OUTI(I1)(1:2).NE.'p2').AND.(OUTI(I1)(1:2).NE.'p3')) THEN
C         <A NAME="error20"></A><FONT COLOR="RED">MTT-20</FONT>
          CALL ERROR('MTT-20: Wrong value of OUT.')
C         This error may be caused by wrongly prescribed value of
C         <A HREF="#data"> input data OUT</A>.
        ELSEIF (((OUTI(I1)(1:2).EQ.' ').AND.(FOUTI(I1).NE.' ')) .OR.
     *          ((OUTI(I1)(1:2).NE.' ').AND.(FOUTI(I1).EQ.' '))) THEN
C         <A NAME="error21"></A><FONT COLOR="RED">MTT-21</FONT>
          CALL ERROR('MTT-21: Inconsistency between FOUT and OUT.')
C         Either both of FOUT(i) and OUT(i) or none of them must be
C         prescribed. See the description of
C         <A HREF="#data"> input data FOUT and OUT</A>.
        ENDIF
   5  CONTINUE
C     Looking, whether travel time is to be interpolated:
      NOUT=0
      DO 11, I1=1,MOUT
        IF (OUTI(I1)(1:2).EQ.'tt') THEN
          IF (NOUT+4.GT.MOUT) THEN
C           <A NAME="error26"></A><FONT COLOR="RED">MTT-26</FONT>
            CALL ERROR('MTT-26: Small array FOUT for the output files.')
C           This error may be caused by wrongly prescribed
C           <A HREF="#data"> input data FOUT and OUT</A>, or by
C           small parameter MOUT defined at the top of this file.
          ENDIF
          FOUT(1)=FOUTI(I1)
          OUT(1)='tt'
          OUT(2)='  '
          OUT(3)='  '
          OUT(4)='  '
          NOUT=4
        ENDIF
  11  CONTINUE
C     Looking, whether any component of the slowness vector
C     is to be interpolated:
      DO 12, I1=1,MOUT
        IF (OUTI(I1)(1:2).EQ.'p1') THEN
          FOUT(2)=FOUTI(I1)
          OUT(2)='p1'
        ELSEIF (OUTI(I1)(1:2).EQ.'p2') THEN
          FOUT(3)=FOUTI(I1)
          OUT(3)='p2'
        ELSEIF (OUTI(I1)(1:2).EQ.'p3') THEN
          FOUT(4)=FOUTI(I1)
          OUT(4)='p3'
        ENDIF
  12  CONTINUE
C     Looking, whether any other quantity is to be interpolated:
      DO 13, I1=1,MOUT
        IF (OUTI(I1)(1:2).EQ.'hi') THEN
          IF (NOUT+1.GT.MOUT) THEN
C           <A NAME="error27"></A><FONT COLOR="RED">MTT-27</FONT>
            CALL ERROR('MTT-27: Small array FOUT for the output files.')
C           This error may be caused by wrongly prescribed
C           <A HREF="#data"> input data FOUT and OUT</A>, or by
C           small parameter MOUT defined at the top of this file.
          ENDIF
          NOUT=NOUT+1
          FOUT(NOUT)=FOUTI(I1)
          OUT(NOUT)=OUTI(I1)(1:2)
        ENDIF
  13  CONTINUE
      NQ=5+NOUT
C
C     Reading in the file with parameters of target grid:
      CALL CIGRID(LU4,FILE4)
C
C     Computing parameter DWARF:
      DWARF=1.
   7  CONTINUE
        DWARF=DWARF*0.1
      IF (1.+DWARF.NE.1.) GOTO 7
      DWARF=DWARF*100.
C     Computing parameter GIANT:
      I1=0
   8  CONTINUE
        I1=I1+100000000
      IF (I1.GT.0) GOTO 8
      GIANT=FLOAT(IABS(I1))
C
C     Reading in file with computed triangles,
C     sorting the rays in triangles:
      WRITE(*,'(''+'',79('' ''))')
      WRITE(*,'(A)') '+Reading the file with 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:
      WRITE(*,'(''+'',79('' ''))')
      WRITE(*,'(A)') '+Sorting the triangles.'
      IF (NRAMP+2*NT.GT.MRAMP) CALL CIEROR(1)
      DO 30, I1=NRAMP+1,NRAMP+NT
        RAM(I1)=FLOAT(IRAM((I1-NRAMP-1)*3+1))
  30  CONTINUE
      CALL INDEXX(NT,RAM(NRAMP+1),IRAM(NRAMP+NT+1))
      DO 40, I1=NRAMP+1,NRAMP+NT
        IRAM(I1)=IRAM(I1+NT)
  40  CONTINUE
      NRAMP=NRAMP+NT
C
C
C     Forming an auxiliary array with information about when can be rays
C     erased from the memory ("deleting array"):
      WRITE(*,'(''+'',79('' ''))')
      WRITE(*,'(A)') '+Forming the ''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-4*NT
      DO 60, I2=3*NT+1,4*NT
        I1=(IRAM(I2)-1)*3+1
        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"):
      WRITE(*,'(''+'',79('' ''))')
      WRITE(*,'(A)') '+Forming the ''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-4*NT-NR-1
C
C
C     Loop for all the triangles:
      IIT1=1
      I1=-1
  100 CONTINUE
        IT1=(IRAM(3*NT+IIT1)-1)*3+1
        I2=NINT((100.*IIT1)/(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
      IIT1=IIT1+1
      IF (IIT1.LT.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           <A NAME="error14"></A><FONT COLOR="RED">MTT-14</FONT>
            CALL ERROR('MTT-14: Disorder in array RAM')
C           This error should not appear.
C           Please <A HREF="#author">contact the author</A>.
          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       <A NAME="error15"></A><FONT COLOR="RED">MTT-15</FONT>
        CALL ERROR('MTT-15: Disorder in array RAM')
C       This error should not appear.
C       Please <A HREF="#author">contact the author</A>.
      ENDIF
      CALL WARRAI(LU5,FILE5,'FORMATTED',.FALSE.,0,.FALSE.,0,N1*N2*N3,
     *                                                    IRAM)
C
C
C     Writing interpolated quantities:
      IF (3*I5.GT.MRAMP) CALL CIEROR(1)
      I4=1
  240 CONTINUE
        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             <A NAME="error16"></A><FONT COLOR="RED">MTT-16</FONT>
              CALL ERROR('MTT-16: Disorder in array RAM')
C             This error should not appear.
C             Please <A HREF="#author">contact the author</A>.
            ENDIF
            NRAMP=NRAMP+1
            IF ((OUT(I4)(1:2).EQ.'tt').OR.(OUT(I4)(1:2).EQ.'p1').OR.
     *          (OUT(I4)(1:2).EQ.'p2').OR.(OUT(I4)(1:2).EQ.'p3')) THEN
C             Travel time or slowness vector:
              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         <A NAME="error17"></A><FONT COLOR="RED">MTT-17</FONT>
          CALL ERROR('MTT-17: Disorder in array RAM')
C         This error should not appear.
C         Please <A HREF="#author">contact the author</A>.
        ENDIF
        IF ((OUT(I4)(1:2).EQ.'tt').OR.(OUT(I4)(1:2).EQ.'p1').OR.
     *      (OUT(I4)(1:2).EQ.'p2').OR.(OUT(I4)(1:2).EQ.'p3')) THEN
C         Travel time or slowness vector:
          CALL WARRAY(LU5,FOUT(I4),'FORMATTED',.FALSE.,0.,.FALSE.,0.,I5,
     *                                                    RAM)
        ELSEIF (OUT(I4)(1:2).EQ.'hi') THEN
C         Ray history:
          CALL WARRAI(LU5,FOUT(I4),'FORMATTED',.FALSE.,0,.FALSE.,0,I5,
     *                                                    IRAM)
        ENDIF
        I4=I4+1
      IF (I4.LE.NOUT) GOTO 240
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           <A NAME="error02"></A><FONT COLOR="RED">MTT-02</FONT>
            CALL ERROR
     *           ('MTT-02: end of file with rays, wanted ray not found')
C           This error may be caused by <A HREF="writ.for#data">K2P </A>
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             <A NAME="error18"></A><FONT COLOR="RED">MTT-18</FONT>
              CALL ERROR
     *              ('MTT-18: a ray with history equal to 0 in CIREAD.')
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             <A HREF="rpar.for#data5">rpar.for.</A>
            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:
            I1=1
  20        CONTINUE
              IF (OUT(I1)(1:2).EQ.'tt') THEN
C               Travel time:
                RAM(NRAMP+(5+I1))=YI(1)
C               Three components of the slowness vector:
                RAM(NRAMP+(5+I1+1))=YI(6)
                RAM(NRAMP+(5+I1+2))=YI(7)
                RAM(NRAMP+(5+I1+3))=YI(8)
              ELSEIF (OUT(I1)(1:2).EQ.'p1') THEN
C               First component of the slowness vector:
                RAM(NRAMP+(5+I1))=YI(6)
              ELSEIF (OUT(I1)(1:2).EQ.'p2') THEN
C               Second component of the slowness vector:
                RAM(NRAMP+(5+I1))=YI(7)
              ELSEIF (OUT(I1)(1:2).EQ.'p3') THEN
C               Third component of the slowness vector:
                RAM(NRAMP+(5+I1))=YI(8)
              ELSEIF (OUT(I1)(1:2).EQ.'hi') THEN
C               Ray history:
                IRAM(NRAMP+(5+I1))=ISHEET
              ENDIF
              I1=I1+1
            IF (I1.LE.NOUT) GOTO 20
            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:
          I1=1
  30      CONTINUE
            IF (OUT(I1)(1:2).EQ.'tt') THEN
C             Travel time:
              RAM(NRAMP+(5+I1))=Y(1)
C             Three components of the slowness vector:
              RAM(NRAMP+(5+I1+1))=Y(6)
              RAM(NRAMP+(5+I1+2))=Y(7)
              RAM(NRAMP+(5+I1+3))=Y(8)
            ELSEIF (OUT(I1)(1:2).EQ.'hi') THEN
C             Ray history:
              IRAM(NRAMP+(5+I1))=ISHEET
            ELSEIF (OUT(I1)(1:2).EQ.'p1') THEN
C             First component of the slowness vector:
              RAM(NRAMP+(5+I1))=Y(6)
            ELSEIF (OUT(I1)(1:2).EQ.'p2') THEN
C             Second component of the slowness vector:
              RAM(NRAMP+(5+I1))=Y(7)
            ELSEIF (OUT(I1)(1:2).EQ.'p3') THEN
C             Third component of the slowness vector:
              RAM(NRAMP+(5+I1))=Y(8)
            ENDIF
            I1=I1+1
          IF (I1.LE.NOUT) GOTO 30
          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 Common block /CIPOIN/:
      INTEGER MPOINT
      PARAMETER (MPOINT=10000)
      INTEGER KPOINT(5,MPOINT),NPOINT
      COMMON/CIPOIN/KPOINT,NPOINT
      SAVE/CIPOIN/
C     KPOINT(1,I1)   ... Index of a point of receiver grid.
C     KPOINT(2:4,I1) ... Indices of four points forming
C                        a side of a ray cell.
C     NPOINT         ... Number of records in KPOINT.
C.......................................................................
C     Auxiliary storage locations:
      INTEGER I1,I2,I3,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)
            DO 5, I3=1,NPOINT
              IF (KPOINT(2,I3).EQ.I2) KPOINT(2,I3)=J1
              IF (KPOINT(3,I3).EQ.I2) KPOINT(3,I3)=J1
              IF (KPOINT(4,I3).EQ.I2) KPOINT(4,I3)=J1
              IF (KPOINT(5,I3).EQ.I2) KPOINT(5,I3)=J1
    5       CONTINUE
   10     CONTINUE
        ELSE
C         This ray is to be erased:
          IRAM(I1-ORAYE)=0
        ENDIF
        IADDRP=IRAM(I1-ORAYA)
        IRAM(I1-ORAYA)=J1
   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       <A NAME="error03"></A><FONT COLOR="RED">MTT-03</FONT>
        CALL ERROR
     *     ('MTT-03: parameters of a ray not found in memory in CITUBE')
C       This error may be caused by <A HREF="writ.for#data">K2P </A>
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             <A NAME="error04"></A><FONT COLOR="RED">MTT-04</FONT>
              CALL ERROR
     *              ('MTT-04: the points reached different interfaces.')
C             This error should not appear.
C             Please <A HREF="#author">contact the author</A>.
            ENDIF
          ENDIF
        ENDIF
        IF ((I1B.GT.I1MA).OR.(I2B.GT.I2MA).OR.(I3B.GT.I3MA)) THEN
C         <A NAME="error05"></A><FONT COLOR="RED">MTT-05</FONT>
          CALL ERROR('MTT-05: point exceded the ray.')
C         This error should not appear.
C         Please <A HREF="#author">contact the author</A>.
        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 CIWB(I1A,I2A,I3A,I1B,I2B,I3B,X1,X2,X3,DWB,IR,ROOT)
C
C----------------------------------------------------------------------
C Subroutine for the computation of the local coordinate wb.
      INTEGER I1A,I2A,I3A,I1B,I2B,I3B,IR
      REAL X1,X2,X3,ROOT(3),DWB
C Input:
C     I1A,I2A,I3A,I1B,I2B,I3B ... Integers specifying the six
C       points on the rays, which 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     X1,X2,X3 ...Coordinates of the examined point.
C     DWB .. The points with local coordinate WB from
C            range [0-DWB,1+DWB] are taken into account.
C Output:
C     IR ... Number of computed values of wb from range [0-DWB,1+DWB].
C     ROOT . Computed values of wb from range [0-DWB,1+DWB].
C
C Subroutines and external functions required:
      EXTERNAL CICUB,CICUBR,CIKVAR,CIDET3,CISUBD
      REAL CICUB,CIDET3,CISUBD
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 I4,I5
      REAL XX(3,3),ZZ(3,3)
      REAL AAA,BBB,CCC,DDD,EEE,FFF
      INTEGER NRR
      REAL RR(4),CRR(4)
      REAL R1,R2
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
      IF (DDD.EQ.0.) THEN
        IF(((ZZ(1,1).EQ.0.).AND.(ZZ(1,2).EQ.0.).AND.(ZZ(1,3).EQ.0.)).OR.
     *     ((ZZ(2,1).EQ.0.).AND.(ZZ(2,2).EQ.0.).AND.(ZZ(2,3).EQ.0.)).OR.
     *     ((ZZ(3,1).EQ.0.).AND.(ZZ(3,2).EQ.0.).AND.(ZZ(3,3).EQ.0.)))
     *     THEN
C         Special case, point X coincides with one of the points forming
C         the top of the ray cell:
          IR=1
          ROOT(1)=1.
          RETURN
        ENDIF
      ENDIF
C
C     Solving the equation according to their coefficients:
      IF (AAA.EQ.0.) THEN
        IF((XX(1,1).EQ.0.).AND.(XX(1,2).EQ.0.).AND.(XX(1,3).EQ.0.).AND.
     *     (XX(2,1).EQ.0.).AND.(XX(2,2).EQ.0.).AND.(XX(2,3).EQ.0.).AND.
     *     (XX(3,1).EQ.0.).AND.(XX(3,2).EQ.0.).AND.(XX(3,3).EQ.0.))
     *     THEN
C         The top of the ray cell coincides with the bottom of the
C         ray cell, point X cannot lie inside the cell:
          IR=0
          RETURN
        ENDIF
        IF (BBB.EQ.0.) THEN
C         Linear equation:
          IF (CCC.EQ.0.) THEN
C           <A NAME="error06"></A><FONT COLOR="RED">MTT-06</FONT>
            CALL ERROR
     *                ('MTT-06: wrong coefficients of equation in CIWB')
C           This error should not appear.
C           Please <A HREF="#author">contact the author</A>.
          ENDIF
          R1=-DDD/CCC
          IF (R1.GE.0.-DWB.AND.R1.LE.1.+DWB) THEN
            IR=1
            ROOT(1)=R1
          ELSE
            IR=0
          ENDIF
        ELSE
C         Kvadratic equation:
          CALL CIKVAR(BBB,CCC,DDD,DWB,IR,R1,R2)
          IF (IR.GE.1) THEN
            IF ((R1.LT.0.-DWB).OR.(R1.GT.1.+DWB)) THEN
C             <A NAME="error11"></A><FONT COLOR="RED">MTT-11</FONT>
              CALL ERROR('MTT-11: First root of quadratic eq. wrong')
C             This error should not appear.
C             Please <A HREF="#author">contact the author</A>.
            ENDIF
            ROOT(1)=R1
          ENDIF
          IF (IR.GE.2) THEN
            IF ((R2.LT.0.-DWB).OR.(R2.GT.1.+DWB)) THEN
C             <A NAME="error12"></A><FONT COLOR="RED">MTT-12</FONT>
              CALL ERROR('MTT-12: Second root of quadratic eq. wrong')
C             This error should not appear.
C             Please <A HREF="#author">contact the author</A>.
            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,DWB,IR,R1,R2)
C       Deciding, whether the cubic equation is to be solved:
        RR(1)=0.-DWB
        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.+DWB
        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.GT.0) THEN
C         Solving the cubic equation:
          CALL CICUBR(AAA,BBB,CCC,DDD,RR,IR,ROOT)
        ENDIF
      ENDIF
      RETURN
      END
C=======================================================================
C
      SUBROUTINE CIWI(I1A,I2A,I3A,I1B,I2B,I3B,X1,X2,X3,WB,WC,
     *                W1,W2,W3)
C
C----------------------------------------------------------------------
C Subroutine for the computation of the local coordinates w1,w2,w3.
      INTEGER I1A,I2A,I3A,I1B,I2B,I3B
      REAL X1,X2,X3,WB,WC,W1,W2,W3
C Input:
C     I1A,I2A,I3A,I1B,I2B,I3B ... Integers specifying the six
C       points on the rays, which 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     X1,X2,X3 ...Coordinates of the examined point.
C     WB,WC ... Already computed local coordinates.
C Output:
C     W1,W2,W3 ... Computed values of the local coordinates W1,W2,W3.
C
C Coded by Petr Bulant
C
C ...........................
C Common block /RAM/ to store the information about triangles and rays:
      INCLUDE 'mtt.inc'
C.......................................................................
      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,AUX
C-----------------------------------------------------------------------
      W1=-1.
      W2=-1.
      W3=-1.
      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
        AUX=(UU22*VV1-UU12*VV2)/DETUU
        W1=AUX
        AUX=(UU11*VV2-UU21*VV1)/DETUU
        W2=AUX
        W3=1.-W1-W2
      ENDIF
      RETURN
      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, which 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,CIAA,CIBB
      REAL CIPPP,CIAA,CIBB
      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.......................................................................
      INTEGER K1MI,K1MA,K2MI,K2MA,K3MI,K3MA
      REAL AUX
      INTEGER I1,I2,I3,I4,I5
      REAL DWB,DWI,HEIG,WID
      PARAMETER (DWB=0.001)
      REAL X1,X2,X3
      REAL Y1,Y2,Y3
      INTEGER IR
      REAL ROOT(3)
      INTEGER INDX
      REAL ERRMAX
      REAL WB,WC,W1,W2,W3
      INTEGER I1AA,I2AA,I1BB,I2BB,ISID
      REAL SID
C     Bicubic interpolation of travel time
      INTEGER I6,I7,I8
      REAL AWB,AWC,TA1,TA2,TA3,PA(3,3),AA(3,3),KSI
C     DWB ... The points with local coordinates wb,wc from
C            range [0-DWB,1+DWB] are taken into account.
C            The points with local coordinates wb,wc from ran-
C            ge [0+DWB,1-DWB] are considered to lie within the ray cell.
C            For the other points further check is done.
C     DWI ... Range [0-DWI,1+DWI] is used for local coordinates w1,w2,w3.
C     HEIG,WID ... Average height and average width of the ray cell.
C            We consider averaged heights of the triangles forming the top
C            and the bottom of the ray cell as the average cell's width.
C-----------------------------------------------------------------------
C     Choosing gridpoints, which may be contained in the ray cell:
      AUX=(AMIN1(RAM(I1A+1),RAM(I2A+1),RAM(I3A+1),
     *           RAM(I1B+1),RAM(I2B+1),RAM(I3B+1)) - O1 ) / D1
      IF (AUX.GT. GIANT) AUX= GIANT
      IF (AUX.LT.-GIANT) AUX=-GIANT
      K1MI=MAX0(0,NINT(AUX))
      AUX=(AMAX1(RAM(I1A+1),RAM(I2A+1),RAM(I3A+1),
     *           RAM(I1B+1),RAM(I2B+1),RAM(I3B+1)) - O1 ) / D1
      IF (AUX.GT. GIANT) AUX= GIANT
      IF (AUX.LT.-GIANT) AUX=-GIANT
      K1MA=MIN0(N1-1,NINT(AUX))
      IF (K1MA.LT.K1MI) RETURN
      AUX=(AMIN1(RAM(I1A+2),RAM(I2A+2),RAM(I3A+2),
     *           RAM(I1B+2),RAM(I2B+2),RAM(I3B+2)) - O2 ) / D2
      IF (AUX.GT. GIANT) AUX= GIANT
      IF (AUX.LT.-GIANT) AUX=-GIANT
      K2MI=MAX0(0,NINT(AUX))
      AUX=(AMAX1(RAM(I1A+2),RAM(I2A+2),RAM(I3A+2),
     *           RAM(I1B+2),RAM(I2B+2),RAM(I3B+2)) - O2 ) / D2
      IF (AUX.GT. GIANT) AUX= GIANT
      IF (AUX.LT.-GIANT) AUX=-GIANT
      K2MA=MIN0(N2-1,NINT(AUX))
      IF (K2MA.LT.K2MI) RETURN
      AUX=(AMIN1(RAM(I1A+3),RAM(I2A+3),RAM(I3A+3),
     *           RAM(I1B+3),RAM(I2B+3),RAM(I3B+3)) - O3 ) / D3
      IF (AUX.GT. GIANT) AUX= GIANT
      IF (AUX.LT.-GIANT) AUX=-GIANT
      K3MI=MAX0(0,NINT(AUX))
      AUX=(AMAX1(RAM(I1A+3),RAM(I2A+3),RAM(I3A+3),
     *           RAM(I1B+3),RAM(I2B+3),RAM(I3B+3)) - O3 ) / D3
      IF (AUX.GT. GIANT) AUX= GIANT
      IF (AUX.LT.-GIANT) AUX=-GIANT
      K3MA=MIN0(N3-1,NINT(AUX))
      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.*DWARF**2
      HEIG=SQRT(((RAM(I1B+1)-RAM(I1A+1))**2+(RAM(I1B+2)-RAM(I1A+2))**2+
     *           (RAM(I1B+3)-RAM(I1A+3))**2+
     *           (RAM(I2B+1)-RAM(I2A+1))**2+(RAM(I2B+2)-RAM(I2A+2))**2+
     *           (RAM(I2B+3)-RAM(I2A+3))**2+
     *           (RAM(I3B+1)-RAM(I3A+1))**2+(RAM(I3B+2)-RAM(I3A+2))**2+
     *           (RAM(I3B+3)-RAM(I3A+3))**2)/3.)
      WID=SQRT(((RAM(I1A+1)-RAM(I2A+1))**2+(RAM(I1A+2)-RAM(I2A+2))**2+
     *          (RAM(I1A+3)-RAM(I2A+3))**2+
     *          (RAM(I2A+1)-RAM(I3A+1))**2+(RAM(I2A+2)-RAM(I3A+2))**2+
     *          (RAM(I2A+3)-RAM(I3A+3))**2+
     *          (RAM(I3A+1)-RAM(I1A+1))**2+(RAM(I3A+2)-RAM(I1A+2))**2+
     *          (RAM(I3A+3)-RAM(I1A+3))**2+
     *          (RAM(I1B+1)-RAM(I2B+1))**2+(RAM(I1B+2)-RAM(I2B+2))**2+
     *          (RAM(I1B+3)-RAM(I2B+3))**2+
     *          (RAM(I2B+1)-RAM(I3B+1))**2+(RAM(I2B+2)-RAM(I3B+2))**2+
     *          (RAM(I2B+3)-RAM(I3B+3))**2+
     *          (RAM(I3B+1)-RAM(I1B+1))**2+(RAM(I3B+2)-RAM(I1B+2))**2+
     *          (RAM(I3B+3)-RAM(I1B+3))**2)/6.)*SQRT(3.)/2.
      DWI=DWB*HEIG/WID
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
        Y1=(RAM(I1A+1)+RAM(I2A+1)+RAM(I3A+1))/3.
        Y2=(RAM(I1A+2)+RAM(I2A+2)+RAM(I3A+2))/3.
        Y3=(RAM(I1A+3)+RAM(I2A+3)+RAM(I3A+3))/3.
        IF (.NOT.CILSIG(
     *    CIPPP(I1B,I2B,I1B,I3B,I1B,Y1,Y2,Y3),
     *    CIPPP(I1B,I2B,I1B,I3B,I1B,X1,X2,X3) , 0. , 0.             ))
     *    GOTO 71
C
        Y1=(RAM(I1B+1)+RAM(I2B+1)+RAM(I3B+1))/3.
        Y2=(RAM(I1B+2)+RAM(I2B+2)+RAM(I3B+2))/3.
        Y3=(RAM(I1B+3)+RAM(I2B+3)+RAM(I3B+3))/3.
        IF (.NOT.CILSIG(
     *    CIPPP(I1A,I2A,I1A,I3A,I1A,Y1,Y2,Y3),
     *    CIPPP(I1A,I2A,I1A,I3A,I1A,X1,X2,X3) , 0. , 0.             ))
     *    GOTO 71
C
        IF (I1A.NE.I1B.AND.I2A.NE.I2B.AND.I3A.NE.I3B) THEN
          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 71
          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 71
          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 71
          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 71
          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 71
          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 71
          ENDIF
        ENDIF
C
C       Computation of the local coordinate wb:
        CALL CIWB(I1A,I2A,I3A,I1B,I2B,I3B,X1,X2,X3,DWB,IR,ROOT)
C
        DO 70, I4=1,IR
          WB=ROOT(I4)
          IF (WB.EQ.1.) GOTO 71
          WC=1.-WB
          IF ((WB.LT.0.-DWB).OR.(WC.LT.0.-DWB)) THEN
C           <A NAME="error13"></A><FONT COLOR="RED">MTT-13</FONT>
            CALL ERROR('MTT-13: Root outside the cell')
C           This error should not appear.
C           Please <A HREF="#author">contact the author</A>.
          ENDIF
          CALL CIWI(I1A,I2A,I3A,I1B,I2B,I3B,X1,X2,X3,WB,WC,W1,W2,W3)
          IF ((W1.LT.0.-DWI).OR.(W1.GT.1.+DWI).OR.
     *        (W2.LT.0.-DWI).OR.(W2.GT.1.+DWI).OR.
     *        (W3.LT.0.-DWI).OR.(W3.GT.1.+DWI))
C           For the considered WB the gridpoint is outside the ray cell.
C           Continuing with the next WB:
     *      GOTO 69
          IF ((W1.GE.0.-DWI).AND.(W1.LE.0.+DWI).AND.
     *        ((I2A.NE.I2B).OR.(I3A.NE.I3B))) THEN
C           The gridpoint is very close the side of the ray cell,
C           checking, whether it is inside or outside the cell:
            I1AA=MIN0(I2A,I3A)
            I2AA=MAX0(I2A,I3A)
            I1BB=MIN0(I2B,I3B)
            I2BB=MAX0(I2B,I3B)
            CALL CILSID
     *               (I1AA,I2AA,I1BB,I2BB,X1,X2,X3,DWB,WB,Y1,Y2,Y3,ISID)
            SID=Y1*(WB*RAM(I1A+1)+WC*RAM(I1B+1)-X1)
     *         +Y2*(WB*RAM(I1A+2)+WC*RAM(I1B+2)-X2)
     *         +Y3*(WB*RAM(I1A+3)+WC*RAM(I1B+3)-X3)
            IF (SID*FLOAT(ISID).LT.0.) THEN
C             For this WB the gridpoint is outside the ray cell:
              GOTO 69
            ENDIF
          ENDIF
          IF ((W2.GE.0.-DWI).AND.(W2.LE.0.+DWI).AND.
     *        ((I1A.NE.I1B).OR.(I3A.NE.I3B))) THEN
C           The gridpoint is very close the side of the ray cell,
C           checking, whether it is inside or outside the cell:
            I1AA=MIN0(I1A,I3A)
            I2AA=MAX0(I1A,I3A)
            I1BB=MIN0(I1B,I3B)
            I2BB=MAX0(I1B,I3B)
            CALL CILSID
     *               (I1AA,I2AA,I1BB,I2BB,X1,X2,X3,DWB,WB,Y1,Y2,Y3,ISID)
            SID=Y1*(WB*RAM(I2A+1)+WC*RAM(I2B+1)-X1)
     *         +Y2*(WB*RAM(I2A+2)+WC*RAM(I2B+2)-X2)
     *         +Y3*(WB*RAM(I2A+3)+WC*RAM(I2B+3)-X3)
            IF (SID*FLOAT(ISID).LT.0.) THEN
C             For this WB the gridpoint is outside the ray cell:
              GOTO 69
            ENDIF
          ENDIF
          IF ((W3.GE.0.-DWI).AND.(W3.LE.0.+DWI).AND.
     *        ((I1A.NE.I1B).OR.(I2A.NE.I2B))) THEN
C           The gridpoint is very close the side of the ray cell,
C           checking, whether it is inside or outside the cell:
            I1AA=MIN0(I2A,I1A)
            I2AA=MAX0(I2A,I1A)
            I1BB=MIN0(I2B,I1B)
            I2BB=MAX0(I2B,I1B)
            CALL CILSID
     *               (I1AA,I2AA,I1BB,I2BB,X1,X2,X3,DWB,WB,Y1,Y2,Y3,ISID)
            SID=Y1*(WB*RAM(I3A+1)+WC*RAM(I3B+1)-X1)
     *         +Y2*(WB*RAM(I3A+2)+WC*RAM(I3B+2)-X2)
     *         +Y3*(WB*RAM(I3A+3)+WC*RAM(I3B+3)-X3)
            IF (SID*FLOAT(ISID).LT.0.) THEN
C             For this WB the gridpoint is outside the ray cell:
              GOTO 69
            ENDIF
          ENDIF
          IF ((W1.GE.0.-DWI).AND.(W1.LE.1.+DWI).AND.
     *        (W2.GE.0.-DWI).AND.(W2.LE.1.+DWI).AND.
     *        (W3.GE.0.-DWI).AND.(W3.LE.1.+DWI)) 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             <A NAME="error23"></A><FONT COLOR="RED">MTT-23</FONT>
cc            CALL ERROR('MTT-23: Inexact interpolation coefficients')
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, check, that
C             you are computing with all possible accuracy.
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
            I5=1
  50        CONTINUE
              IF (OUT(I5)(1:2).EQ.'tt') THEN
C               Travel time - bilinear interpolation:
cc              RAM(NRAMT+I5)=
cc   *             WB*W1*RAM(I1A+(5+I5))+WB*W2*RAM(I2A+(5+I5))+
cc   *             WB*W3*RAM(I3A+(5+I5))+
cc   *             WC*W1*RAM(I1B+(5+I5))+WC*W2*RAM(I2B+(5+I5))+
cc   *             WC*W3*RAM(I3B+(5+I5))
C               Travel time - bicubic interpolation:
C For the interpolation scheme used below refer to:
C Bulant, P. & Klimes, L.: Interpolation of ray-theory travel times
C within ray cells, Seismic Waves in Complex 3-D Structures, Report 7
C Department of Geophysics, Charles University, Prague, 1998.
                AWB=CIAA(WB)
                AWC=CIAA(WC)
                TA1=AWB*RAM(I1A+6)+AWC*RAM(I1B+6)+
     *              CIBB(WB,I1A,I1B)+CIBB(WC,I1B,I1A)
                TA2=AWB*RAM(I2A+6)+AWC*RAM(I2B+6)+
     *              CIBB(WB,I2A,I2B)+CIBB(WC,I2B,I2A)
                TA3=AWB*RAM(I3A+6)+AWC*RAM(I3B+6)+
     *              CIBB(WB,I3A,I3B)+CIBB(WC,I3B,I3A)
                PA(1,1)=AWB*RAM(I1A+7)+AWC*RAM(I1B+7)
                PA(2,1)=AWB*RAM(I1A+8)+AWC*RAM(I1B+8)
                PA(3,1)=AWB*RAM(I1A+9)+AWC*RAM(I1B+9)
                PA(1,2)=AWB*RAM(I2A+7)+AWC*RAM(I2B+7)
                PA(2,2)=AWB*RAM(I2A+8)+AWC*RAM(I2B+8)
                PA(3,2)=AWB*RAM(I2A+9)+AWC*RAM(I2B+9)
                PA(1,3)=AWB*RAM(I3A+7)+AWC*RAM(I3B+7)
                PA(2,3)=AWB*RAM(I3A+8)+AWC*RAM(I3B+8)
                PA(3,3)=AWB*RAM(I3A+9)+AWC*RAM(I3B+9)
                AA(1,1)=WB*RAM(I1A+1)+WC*RAM(I1B+1)
                AA(2,1)=WB*RAM(I1A+2)+WC*RAM(I1B+2)
                AA(3,1)=WB*RAM(I1A+3)+WC*RAM(I1B+3)
                AA(1,2)=WB*RAM(I2A+1)+WC*RAM(I2B+1)
                AA(2,2)=WB*RAM(I2A+2)+WC*RAM(I2B+2)
                AA(3,2)=WB*RAM(I2A+3)+WC*RAM(I2B+3)
                AA(1,3)=WB*RAM(I3A+1)+WC*RAM(I3B+1)
                AA(2,3)=WB*RAM(I3A+2)+WC*RAM(I3B+2)
                AA(3,3)=WB*RAM(I3A+3)+WC*RAM(I3B+3)
                KSI=0
                DO 54, I6=1,3
                  DO 53, I7=1,3
                    DO 52, I8=1,3
                      KSI=KSI+PA(I8,I7)*(AA(I8,I6)-AA(I8,I7))
  52                CONTINUE
  53              CONTINUE
  54            CONTINUE
                KSI=KSI*0.5 + 2*(TA1+TA2+TA3)
                RAM(NRAMT+I5)=
     *          (CIAA(W1)*TA1+CIAA(W2)*TA2+CIAA(W3)*TA3) +
     *          W1**2*(PA(1,1)*(X1-AA(1,1))+PA(2,1)*(X2-AA(2,1))+
     *                 PA(3,1)*(X3-AA(3,1))) +
     *          W2**2*(PA(1,2)*(X1-AA(1,2))+PA(2,2)*(X2-AA(2,2))+
     *                 PA(3,2)*(X3-AA(3,2))) +
     *          W3**2*(PA(1,3)*(X1-AA(1,3))+PA(2,3)*(X2-AA(2,3))+
     *                 PA(3,3)*(X3-AA(3,3))) +
     *          W1*W2*W3*KSI
              ELSEIF ((OUT(I5)(1:2).EQ.'p1').OR.
     *            (OUT(I5)(1:2).EQ.'p2').OR.(OUT(I5)(1:2).EQ.'p3')) THEN
C               Slowness vector - bilinear interpolation:
                RAM(NRAMT+I5)=
     *             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))
              ELSEIF (OUT(I5)(1:2).EQ.'hi') THEN
C               Ray history:
                IRAM(NRAMT+I5)=IRAM(I1A+(5+I5))
              ENDIF
              I5=I5+1
            IF (I5.LE.NOUT) GOTO 50
          ENDIF
  69      CONTINUE
  70    CONTINUE
  71    CONTINUE
  80  CONTINUE
  90  CONTINUE
 100  CONTINUE
      END
C=======================================================================
C
      SUBROUTINE CILSID(I1A,I2A,I1B,I2B,X1,X2,X3,DWB,WB,Y1,Y2,Y3,ISID)
C
C----------------------------------------------------------------------
C Subroutine for the decision on which side of the side of the ray
C cell defined by points I1A,I2A,I1B,I2B is located point X.
C The side of the ray cell defined by points I1A,I2A,I1B,I2B is
C completed with points I3A,I3B (for details see the code bellow) to
C create virtual ray cell. Local coordinates are then computed for
C point X. If (0.-DWB)<WB<(1.+DWB) and W3>0. hold, point X is
C considered to be located at positive side of the side of the ray cell.
C Finaly vector Y with origin in X and direction towards the virtual
C cell is computed.
      INTEGER I1A,I2A,I1B,I2B
      REAL X1,X2,X3,DWB,WB,Y1,Y2,Y3
      INTEGER ISID
C Input:
C     I1A,I2A,I1B,I2B ... Integers specifying the four
C       points on the rays, which define the side of the ray cell.
C       Numbers 1,2 denote the first and second 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     X1,X2,X3 ...Coordinates of the examined point.
C     DWB .. The points with local coordinate WBB from
C            range [0-DWB,1+DWB] are taken into account.
C     WB...  The point with the value of local coordinate WBB nearest
C       to input value WB is considered. (WB is the local coordinate
C       of point X in the real cell, WBB is the local coordinate of X
C       in the virtual ray cell.)
C Output:
CC    Y1,Y2,Y3... Three components of the vector Y from point X
CC      to the point with same WBB lying on the third edge of the virtual
CC      ray cell. The vector points from X towards the virtual cell.
C     Y1,Y2,Y3... Three components of the vector Y from point X
C       towards (pointing inside) the virtual cell.
C     ISID ... +1  if point X is on the positive side,
C              -1  otherwise.
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 I3A,I3B
      REAL A1,A2,A3,B1,B2,B3
      REAL U1,U2,U3,V1,V2,V3,W1,W2,W3,WBB,WCC,A
      INTEGER IR,I,J
      REAL ROOT(3)
      REAL Y11,Y12,Y13,Y21,Y22,Y23
C     I3A,I3B ... Addresses of the two points completing the ray cell.
C     A1,A2,A3,B1,B2,B3 ... Coordinates of two points in the middles of
C                 the edges of the ray cell.
C     U1,U2,U3,V1,V2,V3,W1,W2,W3 ... Coordinates of the vectors at the
C                 side of the ray cell, U and V are the edges, W is
C                 connecting points in the middles of the edges.
C                 W1,W2,W3 are later used as local coordinates.
C     Y11,Y12,Y13 .. Three components of the vector U x W.
C     Y21,Y22,Y23 .. Three components of the vector V x W.
C-----------------------------------------------------------------------
C     Computation of the points I3A, I3B:
      IF (NRAMP+6.GT.MRAMP) CALL CIEROR(1)
      A1=(RAM(I1A+1)+RAM(I2A+1))/2.
      A2=(RAM(I1A+2)+RAM(I2A+2))/2.
      A3=(RAM(I1A+3)+RAM(I2A+3))/2.
      B1=(RAM(I1B+1)+RAM(I2B+1))/2.
      B2=(RAM(I1B+2)+RAM(I2B+2))/2.
      B3=(RAM(I1B+3)+RAM(I2B+3))/2.
      U1=RAM(I1A+1)-RAM(I2A+1)
      U2=RAM(I1A+2)-RAM(I2A+2)
      U3=RAM(I1A+3)-RAM(I2A+3)
      V1=RAM(I1B+1)-RAM(I2B+1)
      V2=RAM(I1B+2)-RAM(I2B+2)
      V3=RAM(I1B+3)-RAM(I2B+3)
      W1=B1-A1
      W2=B2-A2
      W3=B3-A3
      Y11=U2*W3-W2*U3
      Y12=U3*W1-W3*U1
      Y13=U1*W2-W1*U2
      Y21=V2*W3-W2*V3
      Y22=V3*W1-W3*V1
      Y23=V1*W2-W1*V2
      I3A=NRAMP
      RAM(I3A+1)=A1+Y11
      RAM(I3A+2)=A2+Y12
      RAM(I3A+3)=A3+Y13
      I3B=NRAMP+3
      RAM(I3B+1)=B1+Y21
      RAM(I3B+2)=B2+Y22
      RAM(I3B+3)=B3+Y23
C
C     Computation of the local coordinate wb:
      CALL CIWB(I1A,I2A,I3A,I1B,I2B,I3B,X1,X2,X3,DWB,IR,ROOT)
C
      IF(IR.LE.0) THEN
C       <A NAME="error29"></A><FONT COLOR="RED">MTT-29</FONT>
        CALL ERROR('MTT-29: No WB found in the virtual cell')
C       This error should not appear.
C       Please <A HREF="#author">contact the author</A>.
      ENDIF
      J=1
      A=ABS(ROOT(1)-WB)
      DO 10 I=2,IR
        IF(ABS(ROOT(I)-WB).LT.A) THEN
          J=I
          A=ABS(ROOT(I)-WB)
        END IF
   10 CONTINUE
      IF (A.GT.0.01) THEN
C       <A NAME="error30"></A><FONT COLOR="RED">MTT-30</FONT>
        CALL ERROR('MTT-30: WB in virtual and real cells too different')
C       WB was computed for the point X in the real cell. Now we are
C       looking for WBB in the virtual cell, which corresponds to the
C       above mentioned WB in the real cell. Thus WB and WBB should
C       not be too different. It is possible, that the value 0.01
C       should be enlarged.
C
      ENDIF
C
      WBB=ROOT(J)
      WCC=1.-WBB
      IF ((WBB.LT.0.-DWB).OR.(WCC.LT.0.-DWB)) THEN
C       <A NAME="error28"></A><FONT COLOR="RED">MTT-28</FONT>
        CALL ERROR('MTT-28: Root outside the cell')
C       This error should not appear.
C       Please <A HREF="#author">contact the author</A>.
      ENDIF
      CALL CIWI(I1A,I2A,I3A,I1B,I2B,I3B,X1,X2,X3,WBB,WCC,W1,W2,W3)
      IF (W3.GE.0.) THEN
        ISID=+1
      ELSE
        ISID=-1
      ENDIF
C
cc    Y1=WBB*RAM(I3A+1)+WCC*RAM(I3B+1)-X1
cc    Y2=WBB*RAM(I3A+2)+WCC*RAM(I3B+2)-X2
cc    Y3=WBB*RAM(I3A+3)+WCC*RAM(I3B+3)-X3
      Y1=WBB*Y11+WCC*Y21
      Y2=WBB*Y12+WCC*Y22
      Y3=WBB*Y13+WCC*Y23
      RETURN
      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 coordinates.
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 CIAA(W)
C
C----------------------------------------------------------------------
C Subroutine to calculate the value of the cubic function A used for
C interpolation of travel time during bicubic travel time interpolation.
      REAL W
C......................................................................
      CIAA=(3-2*W)*W*W
      END
C
C=======================================================================
C
      REAL FUNCTION CIBB(W,IB,IC)
C
C----------------------------------------------------------------------
C Subroutine to calculate the value of the cubic function B used for
C interpolation of slowness during bicubic travel time interpolation.
      INTEGER IB,IC
      REAL W
C ...........................
C Common block /RAM/ to store the information about triangles and rays:
      INCLUDE 'mtt.inc'
C ...........................
      REAL B1,B2,B3,C1,C2,C3,P1,P2,P3,A1W,A2W,A3W
C......................................................................
      B1=RAM(IB+1)
      B2=RAM(IB+2)
      B3=RAM(IB+3)
      C1=RAM(IC+1)
      C2=RAM(IC+2)
      C3=RAM(IC+3)
      P1=RAM(IB+7)
      P2=RAM(IB+8)
      P3=RAM(IB+9)
      A1W=W*B1+(1-W)*C1
      A2W=W*B2+(1-W)*C2
      A3W=W*B3+(1-W)*C3
      CIBB=W*W*(P1*(A1W-B1)+P2*(A2W-B2)+P3*(A3W-B3))
      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,DWB,IR,R1,R2)
C
C----------------------------------------------------------------------
C Subroutine to solve the kvadratic equation on the interval
C <0-DWB,1+DWB>.
      INTEGER IR
      REAL AAA,BBB,CCC,DWB,R1,R2
C Input:
C     AAA,BBB,CCC ... Coefficients of the equation.
C     DWB .. The roots from range [0-DWB,1+DWB] are taken into account.
C Output:
C     IR ... Number of real roots of the equation from the interval.
C     R1,R2. Real roots of the equation, R1<R2, or undefined.
C
C Coded by Petr Bulant
C
C.......................................................................
      REAL QQ
C----------------------------------------------------------------------
C
C     Solution of kvadratic equation following the algorithm
C     from Numerical Recipes:
C
      QQ=BBB*BBB-4.*AAA*CCC
      IF (QQ.LT.0.) THEN
C       Determinant negative, no real roots:
        IR=0
        RETURN
      ENDIF
      QQ=-0.5*(BBB+SIGN(1.,BBB)*SQRT(QQ))
      R1=QQ/AAA
      R2=CCC/QQ
      IF (R1.GT.R2) THEN
        QQ=R2
        R2=R1
        R1=QQ
      ENDIF
      IF ((R1.GE.0.-DWB).AND.(R1.LE.1.+DWB)) THEN
        IR=1
      ELSE
        IR=0
        R1=R2
      ENDIF
      IF ((R2.GE.0.-DWB).AND.(R2.LE.1.+DWB)) IR=IR+1
      RETURN
      END
C
C
C=======================================================================
C
      SUBROUTINE CICUBR(AA0,BB0,CC0,DD0,RR,IR,ROOT)
C
C----------------------------------------------------------------------
C Subroutine for numerical solution of the cubic equation
C on the interval <RR(1),RR(NRR)>.
      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.
C     RR ... Array of points, which divide the interval <RR(1),RR(NRR)>
C            into subintervals, each subinterval containing just one
C            root of the cubic equation. It is expected, that functional
C            values of cubic equation have different signes in RR(i) and
C            RR(i+1).
C         RR(1)   ... Lower bound of the interval.
C         RR(NRR) ... Upper bound of the interval.
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 the interval 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)
        POM1=CICUB(AA0,BB0,CC0,DD0,RB)
        POM2=CICUB(AA0,BB0,CC0,DD0,RA)
        IF (POM1.EQ.0.) THEN
          ROOT(I2)=RB
          GOTO 19
        ENDIF
        IF (POM2.EQ.0.) THEN
          ROOT(I2)=RA
          GOTO 19
        ENDIF
        I1=0
  10    CONTINUE
          I1=I1+1
          IF (I1.GT.100) THEN
            IF (ABS(RA-RB).LE.ZERO) RETURN
C           <A NAME="error08"></A><FONT COLOR="RED">MTT-08</FONT>
            CALL ERROR('MTT-08: infinite loop in CICUBR')
C           The root was not find even after 100 iterations.
C           This error should not appear.
C           Please <A HREF="#author">contact the author</A>.
          ENDIF
C
          POM1=CICUB(AA0,BB0,CC0,DD0,RB)
          POM2=CICUB(AA0,BB0,CC0,DD0,RA)
          IF (POM1*POM2.GE.0.) THEN
C          <A NAME="error24"></A><FONT COLOR="RED">MTT-24</FONT>
            CALL ERROR('MTT-24: values in RA and RB have the same sign')
C           This error should not appear.
C           Please <A HREF="#author">contact the author</A>.
          ENDIF
          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 (POM3.EQ.0.) GOTO 19
          IF (POM1*POM3.LT.0.) THEN
            RA=ROOT(I2)
            POM2=POM3
          ELSE
            RB=ROOT(I2)
            POM1=POM3
          ENDIF
C
          POM1=CICUB(AA0,BB0,CC0,DD0,RB)
          POM2=CICUB(AA0,BB0,CC0,DD0,RA)
          IF (POM1*POM2.GE.0.) THEN
C          <A NAME="error25"></A><FONT COLOR="RED">MTT-25</FONT>
            CALL ERROR('MTT-25: values in RA and RB have the same sign')
C           This error should not appear.
C           Please <A HREF="#author">contact the author</A>.
          ENDIF
          POM3=POM1-POM2
          IF (POM3.EQ.0.) THEN
C           <A NAME="error09"></A><FONT COLOR="RED">MTT-09</FONT>
            CALL ERROR('MTT-09: division by zero in CICUBR.')
C           This error should not appear.
C           Please <A HREF="#author">contact the author</A>.
          ENDIF
          ROOT(I2)=RA-(RB-RA)*POM2/POM3
          IF((ROOT(I2).LE.RA).OR.(ROOT(I2).GE.RB))ROOT(I2)=(RA+RB)/2.
          POM3=CICUB(AA0,BB0,CC0,DD0,ROOT(I2))
          IF (POM3.EQ.0.) GOTO 19
          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           <A NAME="error22"></A><FONT COLOR="RED">MTT-22</FONT>
            CALL ERROR('MTT-22: Equal functional values in CICUBR.')
C           This error should not appear.
C           Please <A HREF="#author">contact the author</A>.
          END IF
        END IF
  19    CONTINUE
  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.
C
C Input: R1,R2,R3,R4 ... Real quantities.
C Output: CILSIG ... .TRUE. for all the quantities nonnegative or
C                           all the quantities nonpositive.
C                    .FALSE. in other cases.
C-----------------------------------------------------------------------
      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     <A NAME="error10"></A><FONT COLOR="RED">MTT-10</FONT>
      CALL 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.')
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       <A NAME="error01"></A><FONT COLOR="RED">MTT-01</FONT>
        CALL ERROR('MTT-01: array (I)RAM too small.')
C       This error may be caused by too small dimension of array
C       RAM. Try to enlarge dimension MRAM in include file
C       <A HREF="../forms/ram.inc">ram.inc</A>.
      ENDIF
      END
C
C=======================================================================
C
      BLOCK DATA
C
C-----------------------------------------------------------------------
C Common block CIPOIN:
      INTEGER MPOINT
      PARAMETER (MPOINT=10000)
      INTEGER KPOINT(5,MPOINT),NPOINT
      COMMON/CIPOIN/KPOINT,NPOINT
      SAVE/CIPOIN/
      DATA NPOINT/0/
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     <A HREF="../forms/error.for">error.for</A>
      INCLUDE 'ap.for'
C     <A HREF="ap.for">ap.for</A>
      INCLUDE 'forms.for'
C     <A HREF="../forms/forms.for">forms.for</A>
      INCLUDE 'length.for'
C     <A HREF="../forms/length.for">length.for</A>
      INCLUDE 'sep.for'
C     <A HREF="../forms/sep.for">sep.for</A>
      INCLUDE 'indexx.for'
C     <A HREF="../nr/indexx.for">indexx.for</A>
C
C=======================================================================
C                                                                 </PRE>
