C 
C=======================================================================
C
      PROGRAM MTT
C to calculate Multivalued Travel Times on a 3-D grid of points.
C
C Version: 5.30
C Date: 1999, June 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)-st 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 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 The interpolation is not performed within the cells whose bottom
C or top is not formed by different points. This happens mainly
C in the case of point source. To obtain the results of interpolation
C also around the point source, CRT should be run with parameter
C ADVANC.
C.......................................................................
C
C                                                    
C Main input data read from external interactive device (*):
C     The data consist of one character string, read by list
C     directed (free format) input.  The string has thus to be
C     enclosed in apostrophes.
C     The interactive * external unit may be redirected to
C     the file containing the data.
C (1) SEP /
C     Name of the formatted data file with input parameters. If blank,
C     default values of the corresponding data are considered. The file
C     has the form of SEP file, for the description of the SEP format
C     refer to the file  'sep.for'.
C     Only the parameters, which should differ from the default, must be
C     specified in file SEP.
C Default: SEP=' '.
C
C                                                     
C Input parameters in the input data file SEP:
C Names of input files related to ray tracing (output of CRT):
C     CRTOUT='string'... File with the names of the output files of
C             program CRT. If blank, default names of the output files
C             of program CRT are considered.
C             Description of file CRTOUT.
C             Default: CRTOUT=' '.
C Data specifying the parameters of the target grid:
C     O1=real, O2=real, O3=real ... Coordinates of the origin of the
C             grid.
C             Default: O1=0., O2=0., O3=0..
C     D1=real... Grid interval along the X1 axis.
C             Default: D1=1..
C     D2=real... Grid interval along the X2 axis.
C             Default: D2=1..
C     D3=real... Grid interval along the X3 axis.
C             Default: D3=1..
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     PROJ1=real, PROJ2=real, PROJ3=real ... Three components of the
C             projection vector which is used in 2-D computations to
C             project the coordinates of gridpoints to the plane defined
C             by rays. If not specified (default), a vector in the
C             direction of the lowest dimension of the interpolation
C             grid is used as the projection vector.
C             Default: PROJ1=0., PROJ2=0., PROJ3=0..
C Names of output formatted files with interpolated quantities:
C     NUM='string' ...  Name of the output file with N1*N2*N3 array
C             of integer values, containing the number of arrivals at
C             each gridpoint of the target grid.
C             Description of file NUM.
C             Default:  NUM='mtt-num.out'.
C     MTT='string' ... Name of the output file with the array of
C          values, containing for each gridpoint the travel times
C          of all determined arrivals.
C          Description of output files M*.
C          Default:  MTT='mtt-tt.out'.
C     MHIST='string' ... Name of the output file with the array of
C          the ray histories of all determined arrivals.
C          Description of output files M*.
C          Default:  MHIST=' ' (the file is not generated).
C     MP1='string' ... Name of the output file with the array of the
C          first component of the slowness vector.
C     MP2='string' ... Name of the output file with the array of the
C          second component of the slowness vector.
C     MP3='string' ... Name of the output file with the array of the
C          third component of the slowness vector.
C          Description of output files M*.
C          Defaults:  MP1=' ', MP2=' ', MP3=' ' (no files
C          generated).
C     MPQ11='string' ... Name of the output file with the array of the
C          component 1,1 of the 4x4 ray propagator matrix.
C     MPQ21='string' ... Name of the output file with the array of the
C          component 2,1 of the 4x4 ray propagator matrix.
C     MPQ31='string' ... Name of the output file with the array of the
C          component 3,1 of the 4x4 ray propagator matrix.
C     MPQ41='string' ... Name of the output file with the array of the
C          component 4,1 of the 4x4 ray propagator matrix.
C     MPQ12='string' ... Name of the output file with the array of the
C          component 1,2 of the 4x4 ray propagator matrix.
C     MPQ22='string' ... Name of the output file with the array of the
C          component 2,2 of the 4x4 ray propagator matrix.
C     MPQ32='string' ... Name of the output file with the array of the
C          component 3,2 of the 4x4 ray propagator matrix.
C     MPQ42='string' ... Name of the output file with the array of the
C          component 4,2 of the 4x4 ray propagator matrix.
C     MPQ13='string' ... Name of the output file with the array of the
C          component 1,3 of the 4x4 ray propagator matrix.
C     MPQ23='string' ... Name of the output file with the array of the
C          component 2,3 of the 4x4 ray propagator matrix.
C     MPQ33='string' ... Name of the output file with the array of the
C          component 3,3 of the 4x4 ray propagator matrix.
C     MPQ43='string' ... Name of the output file with the array of the
C          component 4,3 of the 4x4 ray propagator matrix.
C     MPQ14='string' ... Name of the output file with the array of the
C          component 1,4 of the 4x4 ray propagator matrix.
C     MPQ24='string' ... Name of the output file with the array of the
C          component 2,4 of the 4x4 ray propagator matrix.
C     MPQ34='string' ... Name of the output file with the array of the
C          component 3,4 of the 4x4 ray propagator matrix.
C     MPQ44='string' ... Name of the output file with the array of the
C             component 4,4 of the 4x4 ray propagator matrix.
C             Description of output files M*.
C             Defaults:  MPQ*=' ' (no files generated).
C
C                                                  
C Input formatted file CRTOUT:
C     The data consist of character strings read by list directed (free
C     format) input. The strings have thus to be enclosed in apostrophes.
C (1) CRT-R, CRT-S, 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-S .. Not used here.
C     CRT-I .. Name of the file with the quantities at the initial
C              points of rays, (see C.R.T.6.1).
C              Description of file CRT-I.
C     CRT-T .. Name of the file with the indices of rays forming
C              the homogeneous triangles in the ray-parameter domain.
C              Description of file CRT-T.
C Default: CRT-R='r01.out', CRT-I='r01i.out', CRT-T='t01.out'.
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 files M*:
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. The loop over gridpoints
C     is the same as in file NUM.
C     If M* is blank, the output file is not generated.
C
C.......................................................................
C
C                                                  
C Coded by Petr Bulant
C     bulant@seis.karlov.mff.cuni.cz
C
C-----------------------------------------------------------------------
C Subroutines and external functions required:
      EXTERNAL LDWARF
      LOGICAL LDWARF
C
C Common block /MTTC/ to store the information about triangles and rays:
      INCLUDE 'mtt.inc'
C mtt.inc.
C ...........................
C Common block /CIGRD/ to store the parameters of the
C interpolation grid.
      REAL O1,O2,O3,D1,D2,D3
      REAL PROJ1,PROJ2,PROJ3
      INTEGER N1,N2,N3
      COMMON/CIGRD/O1,O2,O3,D1,D2,D3,N1,N2,N3,PROJ1,PROJ2,PROJ3
      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     PROJ1,PROJ2,PROJ3 ... Three components of the projection vector,
C                  which is used in 2-D computations to project the
C                  coordinates of gridpoints to the plane defined by
C                  rays.
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)
      REAL AUX
      CHARACTER*80 FILSEP,FILCRT,FILE1,FILE2,FILE3,FILE5,CH
      CHARACTER*52 TXTERR
C     'FOUT(i)'.Name of the output formatted file with the array of
C               interpolated values.
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               'p1' ... First component of the slowness vector.
C               'p2' ... Second component of the slowness vector.
C               'p3' ... Third component of the slowness vector.
C               'hi' ... Ray history.
C               'pq11'...Component 1,1 of the 4x4 ray propagator matrix.
C               'pq21'...Component 2,1 of the 4x4 ray propagator matrix.
C               'pq31'...Component 3,1 of the 4x4 ray propagator matrix.
C               'pq41'...Component 4,1 of the 4x4 ray propagator matrix.
C               'pq12'...Component 1,2 of the 4x4 ray propagator matrix.
C               'pq22'...Component 2,2 of the 4x4 ray propagator matrix.
C               'pq32'...Component 3,2 of the 4x4 ray propagator matrix.
C               'pq42'...Component 4,2 of the 4x4 ray propagator matrix.
C               'pq13'...Component 1,3 of the 4x4 ray propagator matrix.
C               'pq23'...Component 2,3 of the 4x4 ray propagator matrix.
C               'pq33'...Component 3,3 of the 4x4 ray propagator matrix.
C               'pq43'...Component 4,3 of the 4x4 ray propagator matrix.
C               'pq14'...Component 1,4 of the 4x4 ray propagator matrix.
C               'pq24'...Component 2,4 of the 4x4 ray propagator matrix.
C               'pq34'...Component 3,4 of the 4x4 ray propagator matrix.
C               'pq44'...Component 4,4 of the 4x4 ray propagator matrix.
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:
      FILSEP=' '
      WRITE(*,'(A)') ' Enter the name of the input data file SEP: '
      READ(*,*) FILSEP
      WRITE(*,'(''+'',79('' ''))')
      WRITE(*,'(A)') '+MTT: Reading input data.'
C
C     Reading all the data from the SEP file into the memory:
      CALL RSEP1(LU1,FILSEP)
C
C     Reading in filenames of the files with computed rays
C     and triangles:
      CALL RSEP3T('CRTOUT',FILCRT,' ')
      FILE1='r01.out'
      FILE2='r01i.out'
      FILE3='t01.out'
      IF (FILCRT.NE.' ') THEN
        OPEN (LU2,FILE=FILCRT,FORM='FORMATTED',STATUS='OLD')
        READ (LU2,*) FILE1,CH,FILE2,FILE3
        CLOSE (LU2)
      ENDIF
C
C     Opening the CRT output files:
      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
C     Reading in filenames of the output files, recording
C     which quantities are to be written into the files:
      CALL RSEP3T('NUM',FILE5,'mtt-num.out')
      CALL RSEP3T('MTT',FOUT(1),'mtt-tt.out')
      IF (FOUT(1).NE.' ') THEN
        OUT(1)='tt'
        OUT(2)='  '
        OUT(3)='  '
        OUT(4)='  '
        NOUT=4
        CALL RSEP3T('MP1',FOUT(2),' ')
        IF (FOUT(2).NE.' ') OUT(2)='p1'
        CALL RSEP3T('MP2',FOUT(3),' ')
        IF (FOUT(3).NE.' ') OUT(3)='p2'
        CALL RSEP3T('MP3',FOUT(4),' ')
        IF (FOUT(4).NE.' ') OUT(4)='p3'
      ELSE
        NOUT=0
        CALL RSEP3T('MP1',FOUT(NOUT+1),' ')
        IF (FOUT(NOUT+1).NE.' ') THEN
          OUT(NOUT+1)='p1'
          NOUT=NOUT+1
        ENDIF
        CALL RSEP3T('MP2',FOUT(NOUT+1),' ')
        IF (FOUT(NOUT+1).NE.' ') THEN
          OUT(NOUT+1)='p2'
          NOUT=NOUT+1
        ENDIF
        CALL RSEP3T('MP3',FOUT(NOUT+1),' ')
        IF (FOUT(NOUT+1).NE.' ') THEN
          OUT(NOUT+1)='p3'
          NOUT=NOUT+1
        ENDIF
      ENDIF
      CALL RSEP3T('MHIST',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        OUT(NOUT+1)=  'hi'
        NOUT=NOUT+1
      ENDIF
      CALL RSEP3T('MPQ11',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        OUT(NOUT+1)=  'pq11'
        NOUT=NOUT+1
      ENDIF
      CALL RSEP3T('MPQ21',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        OUT(NOUT+1)=  'pq21'
        NOUT=NOUT+1
      ENDIF
      CALL RSEP3T('MPQ31',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        OUT(NOUT+1)=  'pq31'
        NOUT=NOUT+1
      ENDIF
      CALL RSEP3T('MPQ41',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        OUT(NOUT+1)=  'pq41'
        NOUT=NOUT+1
      ENDIF
      CALL RSEP3T('MPQ12',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        OUT(NOUT+1)=  'pq12'
        NOUT=NOUT+1
      ENDIF
      CALL RSEP3T('MPQ22',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        OUT(NOUT+1)=  'pq22'
        NOUT=NOUT+1
      ENDIF
      CALL RSEP3T('MPQ32',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        OUT(NOUT+1)=  'pq32'
        NOUT=NOUT+1
      ENDIF
      CALL RSEP3T('MPQ42',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        OUT(NOUT+1)=  'pq42'
        NOUT=NOUT+1
      ENDIF
      CALL RSEP3T('MPQ13',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        OUT(NOUT+1)=  'pq13'
        NOUT=NOUT+1
      ENDIF
      CALL RSEP3T('MPQ23',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        OUT(NOUT+1)=  'pq23'
        NOUT=NOUT+1
      ENDIF
      CALL RSEP3T('MPQ33',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        OUT(NOUT+1)=  'pq33'
        NOUT=NOUT+1
      ENDIF
      CALL RSEP3T('MPQ43',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        OUT(NOUT+1)=  'pq43'
        NOUT=NOUT+1
      ENDIF
      CALL RSEP3T('MPQ14',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        OUT(NOUT+1)=  'pq14'
        NOUT=NOUT+1
      ENDIF
      CALL RSEP3T('MPQ24',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        OUT(NOUT+1)=  'pq24'
        NOUT=NOUT+1
      ENDIF
      CALL RSEP3T('MPQ34',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        OUT(NOUT+1)=  'pq34'
        NOUT=NOUT+1
      ENDIF
      CALL RSEP3T('MPQ44',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        OUT(NOUT+1)=  'pq44'
        NOUT=NOUT+1
      ENDIF
C
      NQ=5+NOUT
C
C     Reading in the file with parameters of target grid:
      CALL CIGRID
C
C     Computing parameter DWARF:
      DWARF=1.
   7  CONTINUE
        DWARF=DWARF*0.1
        AUX=1.+DWARF
      IF (LDWARF(AUX)) GOTO 7
      DWARF=DWARF*100.
C
C     Reading in file with computed triangles,
C     sorting the rays in triangles:
      WRITE(*,'(''+'',79('' ''))')
      WRITE(*,'(A)') '+MTT: 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.NE.0).AND.(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
        IF (NRAMP+3.GT.MRAMP) CALL CIEROR(1,TXTERR)
        NRAMP=NRAMP+1
        IRAM(NRAMP)=IRAY1
        NRAMP=NRAMP+1
        IRAM(NRAMP)=IRAY2
        NRAMP=NRAMP+1
        IRAM(NRAMP)=IRAY3
        NT=NT+1
      GOTO 10
  20  CONTINUE
      CLOSE(LU3)
      NR=IRAYMA-IRAYMI+1
C
C     Check for the 'triangles' in case of 2D computation:
      IF (IRAM(3).EQ.0) THEN
        DO 25, I1=6,3*NT,3
          IF (IRAM(I1).NE.0) THEN
C           MTT-32
            CALL ERROR('MTT-32: Disorder in triangles.')
C           The first 'triangle' is formed by two rays with zero
C           instead of the third ray. This identifies 2-D calculation
C           and all the 'triangles' should be formed by two rays and
C           zero instead of the third ray. This is not the case of the
C           current triangle. Either file CRT-T
C           is wrong, or there is a bug in the code. In the latter case,
C           please, contact the author.
          ENDIF
  25    CONTINUE
        IF ((PROJ1.EQ.0.).AND.(PROJ2.EQ.0.).AND.(PROJ3.EQ.0.)) THEN
          I1=MIN0(N1,N2,N3)
          I2=0
          IF (N1.EQ.I1) I2=I2+1
          IF (N2.EQ.I1) I2=I2+1
          IF (N3.EQ.I1) I2=I2+1
          IF (I2.GE.2) THEN
C           MTT-33
            CALL ERROR('MTT-33: No projection vector.')
C           The projection vector PROJ1(2,3), described in the
C           description of file SEP, was not
C           specified and the specification of the grid does not enable
C           the use of default projection vector.
          ENDIF
          IF (N1.EQ.I1) PROJ1=1.
          IF (N2.EQ.I1) PROJ2=1.
          IF (N3.EQ.I1) PROJ3=1.
        ENDIF
        L3D=.FALSE.
        L2D=.TRUE.
      ELSE
        L3D=.TRUE.
        L2D=.FALSE.
      ENDIF
C
C     Sorting the triangles:
      WRITE(*,'(''+'',79('' ''))')
      WRITE(*,'(A)') '+MTT: Sorting the triangles.'
      IF (NRAMP+2*NT.GT.MRAMP) CALL CIEROR(1,TXTERR)
      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)') '+MTT: Forming the ''deleting array''.'
      IF (NRAMP+NR.GT.MRAMP) CALL CIEROR(1,TXTERR)
      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)
        IF (IRAM(I1+2).NE.0) 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)') '+MTT: Forming the ''addressing array''.'
C     "Ray" IRAYMI-1:
      NRAMP=NRAMP+1
      IF (NRAMP.GT.MRAMP) CALL CIEROR(1,TXTERR)
      IRAM(NRAMP)=NRAMP+NR
C     All other rays:
      IF (NRAMP+NR.GT.MRAMP) CALL CIEROR(1,TXTERR)
      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)')
     *    '+MTT: Interpolating. ',I2,'% of ray tubes 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.LE.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) THEN
        I1=MRAMP-I5
        WRITE(TXTERR,'(A,I9,A)')
     *  'MTT-31: Array RAM too small,',I1,' units missing.'
        CALL CIEROR(31,TXTERR)
      ENDIF
      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
            CALL ERROR('MTT-14: Disorder in array RAM')
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
        CALL ERROR('MTT-15: Disorder in array RAM')
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) THEN
        I1=MRAMP-I5
        WRITE(TXTERR,'(A,I9,A)')
     *  'MTT-31: Array RAM too small,',I1,' units missing.'
        CALL CIEROR(31,TXTERR)
      ENDIF
      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             MTT-16
              CALL ERROR('MTT-16: Disorder in array RAM')
C             This error should not appear.
C             Please contact the author.
            ENDIF
            NRAMP=NRAMP+1
            IF ((OUT(I4)(1:1).EQ.'t').OR.(OUT(I4)(1:1).EQ.'p')) THEN
C             Travel time, slowness vector or ray propagator matrix:
              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
          CALL ERROR('MTT-17: Disorder in array RAM')
C         This error should not appear.
C         Please contact the author.
        ENDIF
        IF ((OUT(I4)(1:1).EQ.'t').OR.(OUT(I4)(1:1).EQ.'p')) THEN
C         Travel time, slowness vector or ray propagator matrix:
          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 ...........................
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 /MTTC/ to store the information about triangles and rays:
      INCLUDE 'mtt.inc'
C mtt.inc.
C.......................................................................
C
C     Auxiliary storage locations:
      INTEGER IRAY0,I1
      CHARACTER*52 TXTERR
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,TXTERR)
        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.IRAYMA) THEN
C           Some rays at the end of the CRT output file missing:
            DO 12, I1=IRAY0+1,IRAYMA
              IF (I1.GE.IRAYMI) THEN
                IRAM(I1-ORAYA)=NRAMP
              ENDIF
  12        CONTINUE
          ENDIF
          RETURN
        ENDIF
        IF (IRAY.LT.IRAYMI) GOTO 10
        IF (IRAY.GT.IRAYMA) GOTO 10
        IF ((IRAY-IRAY0).GE.2) THEN
C         Some rays skipped by AP00:
          DO 15, I1=IRAY0+1,IRAY-1
            IF (I1.GE.IRAYMI) THEN
              IRAM(I1-ORAYA)=IRAM(IRAY0-ORAYA)
            ENDIF
  15      CONTINUE
        ENDIF
        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
              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             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:
            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:4).EQ.'pq11') THEN
C               Ray propagator matrix:
                RAM(NRAMP+(5+I1))=1.
              ELSEIF (OUT(I1)(1:4).EQ.'pq21') THEN
C               Ray propagator matrix:
                RAM(NRAMP+(5+I1))=0.
              ELSEIF (OUT(I1)(1:4).EQ.'pq31') THEN
C               Ray propagator matrix:
                RAM(NRAMP+(5+I1))=0.
              ELSEIF (OUT(I1)(1:4).EQ.'pq41') THEN
C               Ray propagator matrix:
                RAM(NRAMP+(5+I1))=0.
              ELSEIF (OUT(I1)(1:4).EQ.'pq12') THEN
C               Ray propagator matrix:
                RAM(NRAMP+(5+I1))=0.
              ELSEIF (OUT(I1)(1:4).EQ.'pq22') THEN
C               Ray propagator matrix:
                RAM(NRAMP+(5+I1))=1.
              ELSEIF (OUT(I1)(1:4).EQ.'pq32') THEN
C               Ray propagator matrix:
                RAM(NRAMP+(5+I1))=0.
              ELSEIF (OUT(I1)(1:4).EQ.'pq42') THEN
C               Ray propagator matrix:
                RAM(NRAMP+(5+I1))=0.
              ELSEIF (OUT(I1)(1:4).EQ.'pq13') THEN
C               Ray propagator matrix:
                RAM(NRAMP+(5+I1))=0.
              ELSEIF (OUT(I1)(1:4).EQ.'pq23') THEN
C               Ray propagator matrix:
                RAM(NRAMP+(5+I1))=0.
              ELSEIF (OUT(I1)(1:4).EQ.'pq33') THEN
C               Ray propagator matrix:
                RAM(NRAMP+(5+I1))=1.
              ELSEIF (OUT(I1)(1:4).EQ.'pq43') THEN
C               Ray propagator matrix:
                RAM(NRAMP+(5+I1))=0.
              ELSEIF (OUT(I1)(1:4).EQ.'pq14') THEN
C               Ray propagator matrix:
                RAM(NRAMP+(5+I1))=0.
              ELSEIF (OUT(I1)(1:4).EQ.'pq24') THEN
C               Ray propagator matrix:
                RAM(NRAMP+(5+I1))=0.
              ELSEIF (OUT(I1)(1:4).EQ.'pq34') THEN
C               Ray propagator matrix:
                RAM(NRAMP+(5+I1))=0.
              ELSEIF (OUT(I1)(1:4).EQ.'pq44') THEN
C               Ray propagator matrix:
                RAM(NRAMP+(5+I1))=1.
              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)
            ELSEIF (OUT(I1)(1:4).EQ.'pq11') THEN
C             Ray propagator matrix:
              RAM(NRAMP+(5+I1))=Y(12)
            ELSEIF (OUT(I1)(1:4).EQ.'pq21') THEN
C             Ray propagator matrix:
              RAM(NRAMP+(5+I1))=Y(13)
            ELSEIF (OUT(I1)(1:4).EQ.'pq31') THEN
C             Ray propagator matrix:
              RAM(NRAMP+(5+I1))=Y(14)
            ELSEIF (OUT(I1)(1:4).EQ.'pq41') THEN
C             Ray propagator matrix:
              RAM(NRAMP+(5+I1))=Y(15)
            ELSEIF (OUT(I1)(1:4).EQ.'pq12') THEN
C             Ray propagator matrix:
              RAM(NRAMP+(5+I1))=Y(16)
            ELSEIF (OUT(I1)(1:4).EQ.'pq22') THEN
C             Ray propagator matrix:
              RAM(NRAMP+(5+I1))=Y(17)
            ELSEIF (OUT(I1)(1:4).EQ.'pq32') THEN
C             Ray propagator matrix:
              RAM(NRAMP+(5+I1))=Y(18)
            ELSEIF (OUT(I1)(1:4).EQ.'pq42') THEN
C             Ray propagator matrix:
              RAM(NRAMP+(5+I1))=Y(19)
            ELSEIF (OUT(I1)(1:4).EQ.'pq13') THEN
C             Ray propagator matrix:
              RAM(NRAMP+(5+I1))=Y(20)
            ELSEIF (OUT(I1)(1:4).EQ.'pq23') THEN
C             Ray propagator matrix:
              RAM(NRAMP+(5+I1))=Y(21)
            ELSEIF (OUT(I1)(1:4).EQ.'pq33') THEN
C             Ray propagator matrix:
              RAM(NRAMP+(5+I1))=Y(22)
            ELSEIF (OUT(I1)(1:4).EQ.'pq43') THEN
C             Ray propagator matrix:
              RAM(NRAMP+(5+I1))=Y(23)
            ELSEIF (OUT(I1)(1:4).EQ.'pq14') THEN
C             Ray propagator matrix:
              RAM(NRAMP+(5+I1))=Y(24)
            ELSEIF (OUT(I1)(1:4).EQ.'pq24') THEN
C             Ray propagator matrix:
              RAM(NRAMP+(5+I1))=Y(25)
            ELSEIF (OUT(I1)(1:4).EQ.'pq34') THEN
C             Ray propagator matrix:
              RAM(NRAMP+(5+I1))=Y(26)
            ELSEIF (OUT(I1)(1:4).EQ.'pq44') THEN
C             Ray propagator matrix:
              RAM(NRAMP+(5+I1))=Y(27)
            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 ...........................
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 /MTTC/ to store the information about triangles and rays:
      INCLUDE 'mtt.inc'
C 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
        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 ...........................
      EXTERNAL CIL2P
      LOGICAL CIL2P
C ...........................
C Common block /MTTC/ to store the information about triangles and rays:
      INCLUDE 'mtt.inc'
C mtt.inc.
C.......................................................................
      INTEGER I1B,I2B,I3B,I1C,I2C,I3C
      INTEGER I1MA,I2MA,I3MA,J1,J2
C     I1B,I2B,I3B,I1C,I2C,I3C ... 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 B denotes bottom of
C       the ray cell and 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(I1B+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.
     *    (L3D.AND.(IRAM(IRAM(IT1+2)-ORAYA).EQ.0))) THEN
C       MTT-03
        CALL ERROR
     *     ('MTT-03: Parameters of a ray not found in memory in CITUBE')
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. We reccomend to run CRT with file
C       writall.dat.
      ENDIF
C
      I1MA=IRAM(IRAM(IT1  )-ORAYA)-NQ
      I2MA=IRAM(IRAM(IT1+1)-ORAYA)-NQ
      IF (L3D) I3MA=IRAM(IRAM(IT1+2)-ORAYA)-NQ
C     The first ray cell:
      I1B=IRAM(IRAM(IT1  )-ORAYA-1)
      I2B=IRAM(IRAM(IT1+1)-ORAYA-1)
      IF (L3D) I3B=IRAM(IRAM(IT1+2)-ORAYA-1)
      IF ((I1B.GT.I1MA).OR.(I2B.GT.I2MA).OR.(I3B.GT.I3MA))
     *  RETURN
      I1C=I1B
      I2C=I2B
      IF (L3D) I3C=I3B
C     Loop for points on the rays (loop for ray cells):
C     J1 counts the number of points, which were not shifted.
C     J2 displays the sequential number of points, which is the
C        maximum allowed number for current cell.
  10  CONTINUE
        IF (L3D) THEN
          J1=3
          J2=MAX0(IRAM(I1B+5),IRAM(I2B+5),IRAM(I3B+5))
          IF ((IRAM(I1B+5).EQ.IRAM(I2B+5)).AND.
     *        (IRAM(I1B+5).EQ.IRAM(I3B+5)))  J2=J2+1
        ELSE
          J1=2
          J2=MAX0(IRAM(I1B+5),IRAM(I2B+5))
          IF (IRAM(I1B+5).EQ.IRAM(I2B+5))  J2=J2+1
        ENDIF
C       Forming standard ray cells: ----------------------------------
        IF ((IRAM(I1B+4).EQ.0).AND.(IRAM(I1B+5).LT.J2)) THEN
          I1C=I1B+NQ
          J1=J1-1
        ENDIF
        IF ((IRAM(I2B+4).EQ.0).AND.(IRAM(I2B+5).LT.J2)) THEN
          I2C=I2B+NQ
          J1=J1-1
        ENDIF
        IF (L3D) THEN
          IF ((IRAM(I3B+4).EQ.0).AND.(IRAM(I3B+5).LT.J2)) THEN
            I3C=I3B+NQ
            J1=J1-1
          ENDIF
        ENDIF
ccC       Forming degenerate ray cells (tetrahedrons): -----------------
cc        IF ((IRAM(I1B+4).EQ.0).AND.(IRAM(I1B+5).LT.J2)) THEN
cc          I1C=I1B+NQ
cc          J1=J1-1
cc        ELSE
cc          IF ((IRAM(I2B+4).EQ.0).AND.(IRAM(I2B+5).LT.J2)) THEN
cc            I2C=I2B+NQ
cc            J1=J1-1
cc          ELSE
cc            IF (L3D) THEN
cc              IF ((IRAM(I3B+4).EQ.0).AND.(IRAM(I3B+5).LT.J2)) THEN
cc                I3C=I3B+NQ
cc                J1=J1-1
cc              ENDIF
cc            ENDIF
cc          ENDIF
cc        ENDIF
C       ----------------------------------------------------------------
        IF ((L3D.AND.(J1.EQ.3)).OR.(L2D.AND.(J1.EQ.2))) THEN
          IF (L3D.AND.(IRAM(I1B+4).EQ.IRAM(I2B+4)).AND.
     *                (IRAM(I1B+4).EQ.IRAM(I3B+4))) THEN
C           Crossing the interface in 3D:
            I1B=I1B+NQ
            I2B=I2B+NQ
            I3B=I3B+NQ
            J1=0
            IF ((I1B.GT.I1MA).OR.(I2B.GT.I2MA).OR.(I3B.GT.I3MA))
     *        RETURN
            IRAM(I1B+4)=0
            IRAM(I2B+4)=0
            IRAM(I3B+4)=0
            I1C=I1B
            I2C=I2B
            I3C=I3B
            GOTO 10
          ELSEIF (L2D.AND.(IRAM(I1B+4).EQ.IRAM(I2B+4)).AND.
     *                    (IRAM(I1B+4).EQ.IRAM(I3B+4))) THEN
C           Crossing the interface in 2D:
            I1B=I1B+NQ
            I2B=I2B+NQ
            J1=0
            IF ((I1B.GT.I1MA).OR.(I2B.GT.I2MA))
     *        RETURN
            IRAM(I1B+4)=0
            IRAM(I2B+4)=0
            I1C=I1B
            I2C=I2B
            GOTO 10
          ELSE
C           Moving the points in a complex block:
C           Forming standard ray cells: ------------------------------
            IF (IRAM(I1B+4).EQ.0) THEN
              I1C=I1B+NQ
              J1=J1-1
            ENDIF
            IF (IRAM(I2B+4).EQ.0) THEN
              I2C=I2B+NQ
              J1=J1-1
            ENDIF
            IF (L3D) THEN
              IF (IRAM(I3B+4).EQ.0) THEN
                I3C=I3B+NQ
                J1=J1-1
              ENDIF
            ENDIF
ccC           Forming degenerate ray cells (tetrahedrons): -------------
cc            IF (IRAM(I1B+4).EQ.0) THEN
cc              I1C=I1B+NQ
cc              J1=J1-1
cc            ELSE
cc              IF (IRAM(I2B+4).EQ.0) THEN
cc                I2C=I2B+NQ
cc                J1=J1-1
cc              ELSE
cc                IF (L3D) THEN
cc                  IF (IRAM(I3B+4).EQ.0) THEN
cc                    I3C=I3B+NQ
cc                    J1=J1-1
cc                  ENDIF
cc                ENDIF
cc              ENDIF
cc            ENDIF
C             ----------------------------------------------------------
            IF ((L3D.AND.(J1.EQ.3)).OR.(L2D.AND.(J1.EQ.2))) THEN
              IF ((I1B.GE.I1MA).AND.(I2B.GE.I2MA).AND.(I3B.GE.I3MA).AND.
     *            L3D) RETURN
              IF ((I1B.GE.I1MA).AND.(I2B.GE.I2MA).AND.L2D) RETURN
C             MTT-04
              CALL ERROR
     *              ('MTT-04: The points reached different interfaces.')
C             This error should not appear.
C             Please contact the author.
            ENDIF
          ENDIF
        ENDIF
        IF ((I1C.GT.I1MA).OR.(I2C.GT.I2MA).OR.
     *      (L3D.AND.(I3C.GT.I3MA))) THEN
C         MTT-05
          CALL ERROR('MTT-05: point exceeded the ray.')
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 I1B,I2B,(I3B,)I1C,I2C,(I3C,)
C       is prepared for interpolation:
        IF (L3D) THEN
C         3-D case:
          IF(
     *     ((CIL2P(I1B,I2B).AND.CIL2P(I2B,I3B).AND.CIL2P(I1B,I3B)).AND.
     *      (CIL2P(I1C,I2C).AND.CIL2P(I2C,I3C).AND.CIL2P(I1C,I3C))).AND.
     *      (CIL2P(I1B,I1C).OR.CIL2P(I2B,I2C).OR.CIL2P(I3B,I3C))   )THEN
C           Ray cell formed by 4, 5 or 6 points. The
C           bottom and the top of the ray cell are triangles
C           (i.e. they are formed by different points), the three sides
C           of the ray cell are formed by lines, triangles or tetragons.
            CALL CIINTP(I1B,I2B,I3B,I1C,I2C,I3C)
          ENDIF
        ELSE
C         2-D case:
C         To prevent from calling undefined values:
                                                   I3B=I2B
                                                   I3C=I2C
          IF ((CIL2P(I1B,I2B).OR.CIL2P(I1C,I2C)).AND.
     *        (CIL2P(I1B,I1C).OR.CIL2P(I2B,I2C))) THEN
C           Ray cell formed by 3 or 4 points.
            CALL CIINTP(I1B,I2B,I3B,I1C,I2C,I3C)
          ELSE
C           Strange ray cell:
C           MTT-35
            CALL ERROR('MTT-26: Strange ray cell.')
C           This error should not appear.
C           Please contact the author.
          ENDIF
        ENDIF
        I1B=I1C
        I2B=I2C
        IF (L3D) I3B=I3C
      GOTO 10
C     End of the loop along the ray tube.
      END
C=======================================================================
C
      SUBROUTINE CIWB(I1B,I2B,I3B,I1C,I2C,I3C,X1,X2,X3,DWB,IR,ROOT)
C
C----------------------------------------------------------------------
C Subroutine for the computation of the local coordinate wb.
      INTEGER I1B,I2B,I3B,I1C,I2C,I3C,IR
      REAL X1,X2,X3,ROOT(3),DWB
C Input:
C     I1B,I2B,I3B,I1C,I2C,I3C ... 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 B denotes the bottom of the ray cell and C
C       denotes the 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(I1B+1)
C       In 2-D computations I3B and I3C need not be specified.
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,CIQUAR,CIDET3,CISUBD
      REAL CICUB,CIDET3,CISUBD
C
C ...........................
C Common block /MTTC/ to store the information about triangles and rays:
      INCLUDE 'mtt.inc'
C mtt.inc.
C ...........................
C Common block /CIGRD/ to store the parameters of the
C interpolation grid.
      REAL O1,O2,O3,D1,D2,D3
      REAL PROJ1,PROJ2,PROJ3
      INTEGER N1,N2,N3
      COMMON/CIGRD/O1,O2,O3,D1,D2,D3,N1,N2,N3,PROJ1,PROJ2,PROJ3
      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     PROJ1,PROJ2,PROJ3 ... Three components of the projection vector,
C                  which is used in 2-D computations to project the
C                  coordinates of gridpoints to the plane defined by
C                  rays.
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(I1B+1)-RAM(I1C+1)
      XX(2,1)=RAM(I1B+2)-RAM(I1C+2)
      XX(3,1)=RAM(I1B+3)-RAM(I1C+3)
      XX(1,2)=RAM(I2B+1)-RAM(I2C+1)
      XX(2,2)=RAM(I2B+2)-RAM(I2C+2)
      XX(3,2)=RAM(I2B+3)-RAM(I2C+3)
      IF (L3D) THEN
        XX(1,3)=RAM(I3B+1)-RAM(I3C+1)
        XX(2,3)=RAM(I3B+2)-RAM(I3C+2)
        XX(3,3)=RAM(I3B+3)-RAM(I3C+3)
      ELSE
        XX(1,3)=0.
        XX(2,3)=0.
        XX(3,3)=0.
      ENDIF
C
      ZZ(1,1)=RAM(I1C+1)-X1
      ZZ(2,1)=RAM(I1C+2)-X2
      ZZ(3,1)=RAM(I1C+3)-X3
      ZZ(1,2)=RAM(I2C+1)-X1
      ZZ(2,2)=RAM(I2C+2)-X2
      ZZ(3,2)=RAM(I2C+3)-X3
      IF (L3D) THEN
        ZZ(1,3)=RAM(I3C+1)-X1
        ZZ(2,3)=RAM(I3C+2)-X2
        ZZ(3,3)=RAM(I3C+3)-X3
      ELSE
        ZZ(1,3)=PROJ1
        ZZ(2,3)=PROJ2
        ZZ(3,3)=PROJ3
      ENDIF
C
      AAA=CIDET3(XX)
      IF (L2D.AND.AAA.NE.0.) THEN
C       MTT-35
        CALL ERROR('MTT-35: Cubic equation in 2-D.')
C       This error should not appear.
C       Please contact the author.
      ENDIF
      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(2,1).EQ.0.).AND.(ZZ(3,1).EQ.0.)).OR.
     *     ((ZZ(1,2).EQ.0.).AND.(ZZ(2,2).EQ.0.).AND.(ZZ(3,2).EQ.0.)).OR.
     *     ((ZZ(1,3).EQ.0.).AND.(ZZ(2,3).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(2,1).EQ.0.).AND.(XX(3,1).EQ.0.).AND.
     *     (XX(1,2).EQ.0.).AND.(XX(2,2).EQ.0.).AND.(XX(3,2).EQ.0.).AND.
     *     (XX(1,3).EQ.0.).AND.(XX(2,3).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           The equation has no roots:
            IR=0
            RETURN
          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         Quadratic equation:
          CALL CIQUAR(BBB,CCC,DDD,DWB,IR,R1,R2)
          IF (IR.GE.1) THEN
            IF ((R1.LT.0.-DWB).OR.(R1.GT.1.+DWB)) THEN
C             MTT-11
              CALL ERROR('MTT-11: First root of quadratic eq. wrong')
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.-DWB).OR.(R2.GT.1.+DWB)) THEN
C             MTT-12
              CALL ERROR('MTT-12: Second root of quadratic eq. wrong')
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 CIQUAR(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(I1B,I2B,I3B,I1C,I2C,I3C,X1,X2,X3,WB,WC,
     *                W1,W2,W3)
C
C----------------------------------------------------------------------
C Subroutine for the computation of the local coordinates w1,w2,w3.
      INTEGER I1B,I2B,I3B,I1C,I2C,I3C
      REAL X1,X2,X3,WB,WC,W1,W2,W3
C Input:
C     I1B,I2B,I3B,I1C,I2C,I3C ... 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 B denotes the bottom of the ray cell and C
C       denotes the 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(I1B+1)
C       In 2-D computations I3B and I3C need not be specified.
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                  In 2-D computations W3 is not computed.
C
C ...........................
C Common block /MTTC/ to store the information about triangles and rays:
      INCLUDE 'mtt.inc'
C mtt.inc.
C ...........................
C Common block /CIGRD/ to store the parameters of the
C interpolation grid.
      REAL O1,O2,O3,D1,D2,D3
      REAL PROJ1,PROJ2,PROJ3
      INTEGER N1,N2,N3
      COMMON/CIGRD/O1,O2,O3,D1,D2,D3,N1,N2,N3,PROJ1,PROJ2,PROJ3
      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     PROJ1,PROJ2,PROJ3 ... Three components of the projection vector,
C                  which is used in 2-D computations to project the
C                  coordinates of gridpoints to the plane defined by
C                  rays.
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
      REAL XA1,XA2,XA3,AA1,AA2,AA3,E1,E2,E3,EE
      EQUIVALENCE (E1,PROJ1),(E2,PROJ2),(E3,PROJ3)
C-----------------------------------------------------------------------
      W1=-1.
      W2=-1.
      W3=-1.
      A11=WB*RAM(I1B+1) + WC*RAM(I1C+1)
      A21=WB*RAM(I1B+2) + WC*RAM(I1C+2)
      A31=WB*RAM(I1B+3) + WC*RAM(I1C+3)
      A12=WB*RAM(I2B+1) + WC*RAM(I2C+1)
      A22=WB*RAM(I2B+2) + WC*RAM(I2C+2)
      A32=WB*RAM(I2B+3) + WC*RAM(I2C+3)
      IF (L2D) THEN
C       2-D case:
        XA1=X1-A12
        XA2=X2-A22
        XA3=X3-A32
        AA1=A11-A12
        AA2=A21-A22
        AA3=A31-A32
        EE=E1*E1+E2*E2+E3*E3
        W1=(XA1*AA1+XA2*AA2+XA3*AA3)*EE -
     *     (XA1*E1+XA2*E2+XA3*E3)*(AA1*E1+AA2*E2+AA3*E3)
        W1=W1/((AA1*AA1+AA2*AA2+AA3*AA3)*EE-(AA1*E1+AA2*E2+AA3*E3)**2)
        W2=1. - W1
        IF (ABS(W1-1.).LT.DWARF) THEN
          XA1=X1-A11
          XA2=X2-A21
          XA3=X3-A31
          AA1=A12-A11
          AA2=A22-A21
          AA3=A32-A31
          EE=E1*E1+E2*E2+E3*E3
          W2=(XA1*AA1+XA2*AA2+XA3*AA3)*EE -
     *       (XA1*E1+XA2*E2+XA3*E3)*(AA1*E1+AA2*E2+AA3*E3)
          W2=W2/((AA1*AA1+AA2*AA2+AA3*AA3)*EE-(AA1*E1+AA2*E2+AA3*E3)**2)
          W1=1. - W2
        ENDIF
      ELSE
C       3-D case:
        A13=WB*RAM(I3B+1) + WC*RAM(I3C+1)
        A23=WB*RAM(I3B+2) + WC*RAM(I3C+2)
        A33=WB*RAM(I3B+3) + WC*RAM(I3C+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
      ENDIF
      RETURN
      END
C=======================================================================
C
      SUBROUTINE CIINTP(I1B,I2B,I3B,I1C,I2C,I3C)
C
C----------------------------------------------------------------------
C Subroutine for interpolation within standard ray cell formed by the
C points I1B,I2B,I3B,I1C,I2C,I3C in 3-D or by points I1B,I2B,I1C,I2C
C in 2-D.
C
      INTEGER I1B,I2B,I3B,I1C,I2C,I3C
C Input:
C     I1B,I2B,I3B,I1C,I2C,I3C ... 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 B denotes the bottom of the ray cell and C
C       denotes the 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(I1B+1)
C No output.
C
C Subroutines and external functions required:
      EXTERNAL CIPPP,CILSIG,CIAA,CIBB
      REAL CIPPP,CIAA,CIBB
      LOGICAL CILSIG
C
C ...........................
C Common block /MTTC/ to store the information about triangles and rays:
      INCLUDE 'mtt.inc'
C mtt.inc.
C ...........................
C Common block /CIGRD/ to store the parameters of the
C interpolation grid.
      REAL O1,O2,O3,D1,D2,D3
      REAL PROJ1,PROJ2,PROJ3
      INTEGER N1,N2,N3
      COMMON/CIGRD/O1,O2,O3,D1,D2,D3,N1,N2,N3,PROJ1,PROJ2,PROJ3
      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     PROJ1,PROJ2,PROJ3 ... Three components of the projection vector,
C                  which is used in 2-D computations to project the
C                  coordinates of gridpoints to the plane defined by
C                  rays.
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 I1BA,I2BA,I1CB,I2CB,ISID
      REAL SID
      CHARACTER*52 TXTERR
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            Averaged heights of the triangles forming the top
C            and the bottom of the ray cell are considered as the average
C            cell's width.
C-----------------------------------------------------------------------
C     Choosing gridpoints, which may be contained in the ray cell:
      IF (L2D) THEN
        I1=0
        IF (PROJ1.EQ.0.) I1=I1+1
        IF (PROJ2.EQ.0.) I1=I1+1
        IF (PROJ3.EQ.0.) I1=I1+1
        IF (I1.EQ.3) THEN
C         In this case all the gridpoints may be inside the cell:
          K1MI=0
          K1MA=N1-1
          K2MI=0
          K2MA=N2-1
          K3MI=0
          K3MA=N3-1
          GOTO 10
        ENDIF
      ENDIF
      AUX=(AMIN1(RAM(I1B+1),RAM(I2B+1),RAM(I3B+1),
     *           RAM(I1C+1),RAM(I2C+1),RAM(I3C+1)) - O1 ) / D1
      IF (AUX.GT. GIANT) AUX= GIANT
      IF (AUX.LT.-GIANT) AUX=-GIANT
      K1MI=MAX0(0,NINT(AUX))
      IF (L2D.AND.(PROJ1.NE.0.)) K1MI=0
      AUX=(AMAX1(RAM(I1B+1),RAM(I2B+1),RAM(I3B+1),
     *           RAM(I1C+1),RAM(I2C+1),RAM(I3C+1)) - O1 ) / D1
      IF (AUX.GT. GIANT) AUX= GIANT
      IF (AUX.LT.-GIANT) AUX=-GIANT
      K1MA=MIN0(N1-1,NINT(AUX))
      IF (L2D.AND.(PROJ1.NE.0.)) K1MA=N1-1
      IF (K1MA.LT.K1MI) RETURN
      AUX=(AMIN1(RAM(I1B+2),RAM(I2B+2),RAM(I3B+2),
     *           RAM(I1C+2),RAM(I2C+2),RAM(I3C+2)) - O2 ) / D2
      IF (AUX.GT. GIANT) AUX= GIANT
      IF (AUX.LT.-GIANT) AUX=-GIANT
      K2MI=MAX0(0,NINT(AUX))
      IF (L2D.AND.(PROJ2.NE.0.)) K2MI=0
      AUX=(AMAX1(RAM(I1B+2),RAM(I2B+2),RAM(I3B+2),
     *           RAM(I1C+2),RAM(I2C+2),RAM(I3C+2)) - O2 ) / D2
      IF (AUX.GT. GIANT) AUX= GIANT
      IF (AUX.LT.-GIANT) AUX=-GIANT
      K2MA=MIN0(N2-1,NINT(AUX))
      IF (L2D.AND.(PROJ2.NE.0.)) K2MA=N2-1
      IF (K2MA.LT.K2MI) RETURN
      AUX=(AMIN1(RAM(I1B+3),RAM(I2B+3),RAM(I3B+3),
     *           RAM(I1C+3),RAM(I2C+3),RAM(I3C+3)) - O3 ) / D3
      IF (AUX.GT. GIANT) AUX= GIANT
      IF (AUX.LT.-GIANT) AUX=-GIANT
      K3MI=MAX0(0,NINT(AUX))
      IF (L2D.AND.(PROJ3.NE.0.)) K3MI=0
      AUX=(AMAX1(RAM(I1B+3),RAM(I2B+3),RAM(I3B+3),
     *           RAM(I1C+3),RAM(I2C+3),RAM(I3C+3)) - O3 ) / D3
      IF (AUX.GT. GIANT) AUX= GIANT
      IF (AUX.LT.-GIANT) AUX=-GIANT
      K3MA=MIN0(N3-1,NINT(AUX))
      IF (L2D.AND.(PROJ3.NE.0.)) K3MA=N3-1
      IF (K3MA.LT.K3MI) RETURN
C
  10  CONTINUE
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 (L3D.AND.I1B.NE.I1C.AND.I2B.NE.I2C.AND.I3B.NE.I3C) THEN
        IF (.NOT.CILSIG(
     *    CIPPP(I1B,I2B,I1B,I3B,I1B,RAM(I1C+1),RAM(I1C+2),RAM(I1C+3)),
     *    CIPPP(I1B,I2B,I1B,I3B,I1B,RAM(I2C+1),RAM(I2C+2),RAM(I2C+3)),
     *    CIPPP(I1B,I2B,I1B,I3B,I1B,RAM(I3C+1),RAM(I3C+2),RAM(I3C+3)),
     *          0.)) RETURN
        IF (.NOT.CILSIG(
     *    CIPPP(I1C,I2C,I1C,I3C,I1C,RAM(I1B+1),RAM(I1B+2),RAM(I1B+3)),
     *    CIPPP(I1C,I2C,I1C,I3C,I1C,RAM(I2B+1),RAM(I2B+2),RAM(I2B+3)),
     *    CIPPP(I1C,I2C,I1C,I3C,I1C,RAM(I3B+1),RAM(I3B+2),RAM(I3B+3)),
     *          0.)) RETURN
      ENDIF
C
      AUX=4.
      ERRMAX=RAM(I1B+1)**2+RAM(I2B+1)**2+RAM(I3B+1)**2
     *      +RAM(I1C+1)**2+RAM(I2C+1)**2+RAM(I3C+1)**2
     *      +RAM(I1B+2)**2+RAM(I2B+2)**2+RAM(I3B+2)**2
     *      +RAM(I1C+2)**2+RAM(I2C+2)**2+RAM(I3C+2)**2
      IF (L3D) THEN
        ERRMAX=ERRMAX
     *      +RAM(I1B+3)**2+RAM(I2B+3)**2+RAM(I3B+3)**2
     *      +RAM(I1C+3)**2+RAM(I2C+3)**2+RAM(I3C+3)**2
        AUX=6.
      ENDIF
      ERRMAX=ERRMAX/AUX*DWARF**2
C
      AUX=2.
      HEIG=      (RAM(I1C+1)-RAM(I1B+1))**2+(RAM(I1C+2)-RAM(I1B+2))**2+
     *           (RAM(I1C+3)-RAM(I1B+3))**2+
     *           (RAM(I2C+1)-RAM(I2B+1))**2+(RAM(I2C+2)-RAM(I2B+2))**2+
     *           (RAM(I2C+3)-RAM(I2B+3))**2
      IF (L3D)THEN
        HEIG=HEIG+(RAM(I3C+1)-RAM(I3B+1))**2+(RAM(I3C+2)-RAM(I3B+2))**2
     *           +(RAM(I3C+3)-RAM(I3B+3))**2
        AUX=3.
      ENDIF
      HEIG=SQRT(HEIG/AUX)
C
      AUX=2.
      WID=      (RAM(I1B+1)-RAM(I2B+1))**2+(RAM(I1B+2)-RAM(I2B+2))**2+
     *          (RAM(I1B+3)-RAM(I2B+3))**2+
     *          (RAM(I1C+1)-RAM(I2C+1))**2+(RAM(I1C+2)-RAM(I2C+2))**2+
     *          (RAM(I1C+3)-RAM(I2C+3))**2
      IF (L3D) THEN
        WID=WID+(RAM(I2B+1)-RAM(I3B+1))**2+(RAM(I2B+2)-RAM(I3B+2))**2+
     *          (RAM(I2B+3)-RAM(I3B+3))**2+
     *          (RAM(I2C+1)-RAM(I3C+1))**2+(RAM(I2C+2)-RAM(I3C+2))**2+
     *          (RAM(I2C+3)-RAM(I3C+3))**2+
     *          (RAM(I3B+1)-RAM(I1B+1))**2+(RAM(I3B+2)-RAM(I1B+2))**2+
     *          (RAM(I3B+3)-RAM(I1B+3))**2+
     *          (RAM(I3C+1)-RAM(I1C+1))**2+(RAM(I3C+2)-RAM(I1C+2))**2+
     *          (RAM(I3C+3)-RAM(I1C+3))**2
        AUX=6./(SQRT(3.)/2.)
      ENDIF
      WID=SQRT(WID/AUX)
C
      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
        IF (L3D) THEN
C         Checking the position of the gridpoint with respect to
C         the planes bounding the ray cell:
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(I1C,I2C,I1C,I3C,I1C,Y1,Y2,Y3),
     *      CIPPP(I1C,I2C,I1C,I3C,I1C,X1,X2,X3) , 0. , 0.             ))
     *      GOTO 71
C
          Y1=(RAM(I1C+1)+RAM(I2C+1)+RAM(I3C+1))/3.
          Y2=(RAM(I1C+2)+RAM(I2C+2)+RAM(I3C+2))/3.
          Y3=(RAM(I1C+3)+RAM(I2C+3)+RAM(I3C+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
          IF (I1B.NE.I1C.AND.I2B.NE.I2C.AND.I3B.NE.I3C) THEN
            IF (CILSIG(
     *      CIPPP(I1B,I2C,I2B,I1C,I1B,RAM(I2B+1),RAM(I2B+2),RAM(I2B+3)),
     *      CIPPP(I1B,I2C,I2B,I1C,I1B,RAM(I3B+1),RAM(I3B+2),RAM(I3B+3)),
     *      CIPPP(I1B,I2C,I2B,I1C,I1B,RAM(I3C+1),RAM(I3C+2),RAM(I3C+3)),
     *            0.)) THEN
              IF (.NOT.CILSIG(
     *      CIPPP(I1B,I2C,I2B,I1C,I1B,RAM(I2B+1),RAM(I2B+2),RAM(I2B+3)),
     *      CIPPP(I1B,I2C,I2B,I1C,I1B,X1,X2,X3),0.,0.))          GOTO 71
            ELSEIF (CILSIG(
     *      CIPPP(I1B,I2C,I2B,I1C,I2B,RAM(I1B+1),RAM(I1B+2),RAM(I1B+3)),
     *      CIPPP(I1B,I2C,I2B,I1C,I2B,RAM(I3B+1),RAM(I3B+2),RAM(I3B+3)),
     *      CIPPP(I1B,I2C,I2B,I1C,I2B,RAM(I3C+1),RAM(I3C+2),RAM(I3C+3)),
     *            0.)) THEN
              IF (.NOT.CILSIG(
     *      CIPPP(I1B,I2C,I2B,I1C,I2B,RAM(I1B+1),RAM(I1B+2),RAM(I1B+3)),
     *      CIPPP(I1B,I2C,I2B,I1C,I2B,X1,X2,X3),0.,0.))          GOTO 71
            ENDIF
C
            IF (CILSIG(
     *      CIPPP(I2B,I3C,I3B,I2C,I2B,RAM(I3B+1),RAM(I3B+2),RAM(I3B+3)),
     *      CIPPP(I2B,I3C,I3B,I2C,I2B,RAM(I1B+1),RAM(I1B+2),RAM(I1B+3)),
     *      CIPPP(I2B,I3C,I3B,I2C,I2B,RAM(I1C+1),RAM(I1C+2),RAM(I1C+3)),
     *            0.)) THEN
              IF (.NOT.CILSIG(
     *      CIPPP(I2B,I3C,I3B,I2C,I2B,RAM(I3B+1),RAM(I3B+2),RAM(I3B+3)),
     *      CIPPP(I2B,I3C,I3B,I2C,I2B,X1,X2,X3),0.,0.))          GOTO 71
            ELSEIF (CILSIG(
     *      CIPPP(I2B,I3C,I3B,I2C,I3B,RAM(I2B+1),RAM(I2B+2),RAM(I2B+3)),
     *      CIPPP(I2B,I3C,I3B,I2C,I3B,RAM(I1B+1),RAM(I1B+2),RAM(I1B+3)),
     *      CIPPP(I2B,I3C,I3B,I2C,I3B,RAM(I1C+1),RAM(I1C+2),RAM(I1C+3)),
     *            0.)) THEN
              IF (.NOT.CILSIG(
     *      CIPPP(I2B,I3C,I3B,I2C,I3B,RAM(I2B+1),RAM(I2B+2),RAM(I2B+3)),
     *      CIPPP(I2B,I3C,I3B,I2C,I3B,X1,X2,X3),0.,0.))          GOTO 71
            ENDIF
C
            IF (CILSIG(
     *      CIPPP(I3B,I1C,I1B,I3C,I1B,RAM(I3B+1),RAM(I3B+2),RAM(I3B+3)),
     *      CIPPP(I3B,I1C,I1B,I3C,I1B,RAM(I2B+1),RAM(I2B+2),RAM(I2B+3)),
     *      CIPPP(I3B,I1C,I1B,I3C,I1B,RAM(I2C+1),RAM(I2C+2),RAM(I2C+3)),
     *            0.)) THEN
              IF (.NOT.CILSIG(
     *      CIPPP(I3B,I1C,I1B,I3C,I1B,RAM(I3B+1),RAM(I3B+2),RAM(I3B+3)),
     *      CIPPP(I3B,I1C,I1B,I3C,I1B,X1,X2,X3),0.,0.))          GOTO 71
            ELSEIF (CILSIG(
     *      CIPPP(I3B,I1C,I1B,I3C,I3B,RAM(I1B+1),RAM(I1B+2),RAM(I1B+3)),
     *      CIPPP(I3B,I1C,I1B,I3C,I3B,RAM(I2B+1),RAM(I2B+2),RAM(I2B+3)),
     *      CIPPP(I3B,I1C,I1B,I3C,I3B,RAM(I2C+1),RAM(I2C+2),RAM(I2C+3)),
     *            0.)) THEN
              IF (.NOT.CILSIG(
     *      CIPPP(I3B,I1C,I1B,I3C,I3B,RAM(I1B+1),RAM(I1B+2),RAM(I1B+3)),
     *      CIPPP(I3B,I1C,I1B,I3C,I3B,X1,X2,X3),0.,0.))          GOTO 71
            ENDIF
          ENDIF
        ENDIF
C
C       Computation of the local coordinate wb:
        CALL CIWB(I1B,I2B,I3B,I1C,I2C,I3C,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           MTT-13
            CALL ERROR('MTT-13: Root outside the cell')
C           This error should not appear.
C           Please contact the author.
          ENDIF
          CALL CIWI(I1B,I2B,I3B,I1C,I2C,I3C,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.
     *        (L3D.AND.((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 (L3D) THEN
C           Tests for the points close to the sides of the cell:
            IF ((W1.GE.0.-DWI).AND.(W1.LE.0.+DWI).AND.
     *          ((I2B.NE.I2C).OR.(I3B.NE.I3C))) THEN
C             The gridpoint is very close the side of the ray cell,
C             checking, whether it is inside or outside the cell:
              I1BA=MIN0(I2B,I3B)
              I2BA=MAX0(I2B,I3B)
              I1CB=MIN0(I2C,I3C)
              I2CB=MAX0(I2C,I3C)
              CALL CILSID
     *             (I1BA,I2BA,I1CB,I2CB,X1,X2,X3,DWB,WB,Y1,Y2,Y3,ISID)
              SID=Y1*(WB*RAM(I1B+1)+WC*RAM(I1C+1)-X1)
     *           +Y2*(WB*RAM(I1B+2)+WC*RAM(I1C+2)-X2)
     *           +Y3*(WB*RAM(I1B+3)+WC*RAM(I1C+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.
     *          ((I1B.NE.I1C).OR.(I3B.NE.I3C))) THEN
C             The gridpoint is very close the side of the ray cell,
C             checking, whether it is inside or outside the cell:
              I1BA=MIN0(I1B,I3B)
              I2BA=MAX0(I1B,I3B)
              I1CB=MIN0(I1C,I3C)
              I2CB=MAX0(I1C,I3C)
              CALL CILSID
     *             (I1BA,I2BA,I1CB,I2CB,X1,X2,X3,DWB,WB,Y1,Y2,Y3,ISID)
              SID=Y1*(WB*RAM(I2B+1)+WC*RAM(I2C+1)-X1)
     *           +Y2*(WB*RAM(I2B+2)+WC*RAM(I2C+2)-X2)
     *           +Y3*(WB*RAM(I2B+3)+WC*RAM(I2C+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.
     *          ((I1B.NE.I1C).OR.(I2B.NE.I2C))) THEN
C             The gridpoint is very close the side of the ray cell,
C             checking, whether it is inside or outside the cell:
              I1BA=MIN0(I2B,I1B)
              I2BA=MAX0(I2B,I1B)
              I1CB=MIN0(I2C,I1C)
              I2CB=MAX0(I2C,I1C)
              CALL CILSID
     *             (I1BA,I2BA,I1CB,I2CB,X1,X2,X3,DWB,WB,Y1,Y2,Y3,ISID)
              SID=Y1*(WB*RAM(I3B+1)+WC*RAM(I3C+1)-X1)
     *           +Y2*(WB*RAM(I3B+2)+WC*RAM(I3C+2)-X2)
     *           +Y3*(WB*RAM(I3B+3)+WC*RAM(I3C+3)-X3)
              IF (SID*FLOAT(ISID).LT.0.) THEN
C               For this WB the gridpoint is outside the ray cell:
                GOTO 69
              ENDIF
            ENDIF
          ENDIF
          IF (L2D) THEN
C           ------------------------------------------------------------
ccC           No tests for the points close to the sides of the cell:
cc            IF ((W1.LT.0.).OR.(W1.GE.1.).OR.
cc     *          (W2.LE.0.).OR.(W2.GT.1.).OR.
cc     *          (WB.LT.0.).OR.(WB.GE.1.))
ccC               The gridpoint is situated outside the ray cell.
cc     *          GOTO 69
C           ------------------------------------------------------------
C           Tests for the points close to the sides of the cell:
            IF ((WB.LT.0.).OR.(WB.GE.1.))
C               The gridpoint is situated outside the ray cell.
     *          GOTO 69
            IF (ABS(W1-1.).LT.DWI) THEN
C             The gridpoint is very close the side of the ray cell,
C             checking, whether it is inside or outside the cell:
              I1BA=I1B
              I1CB=I1C
              CALL CILSID
     *             (I1BA,I2BA,I1CB,I2CB,X1,X2,X3,DWB,WB,Y1,Y2,Y3,ISID)
              SID=Y1*(WB*RAM(I2B+1)+WC*RAM(I2C+1)-X1)
     *           +Y2*(WB*RAM(I2B+2)+WC*RAM(I2C+2)-X2)
     *           +Y3*(WB*RAM(I2B+3)+WC*RAM(I2C+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 (ABS(W2-1.).LT.DWI) THEN
C             The gridpoint is very close the side of the ray cell,
C             checking, whether it is inside or outside the cell:
              I1BA=I2B
              I1CB=I2C
              CALL CILSID
     *             (I1BA,I2BA,I1CB,I2CB,X1,X2,X3,DWB,WB,Y1,Y2,Y3,ISID)
              SID=Y1*(WB*RAM(I1B+1)+WC*RAM(I1C+1)-X1)
     *           +Y2*(WB*RAM(I1B+2)+WC*RAM(I1C+2)-X2)
     *           +Y3*(WB*RAM(I1B+3)+WC*RAM(I1C+3)-X3)
              IF (SID*FLOAT(ISID).LT.0.) THEN
C               For this WB the gridpoint is outside the ray cell:
                GOTO 69
              ENDIF
            ENDIF
C           ------------------------------------------------------------
          ENDIF
C
C         The gridpoint is situated inside the ray cell:
C
          IF (L3D) THEN
C           Checking accuracy of interpolation coefficients:
            Y1=WB*(W1*RAM(I1B+1)+W2*RAM(I2B+1)+W3*RAM(I3B+1))+
     *         WC*(W1*RAM(I1C+1)+W2*RAM(I2C+1)+W3*RAM(I3C+1))
            Y2=WB*(W1*RAM(I1B+2)+W2*RAM(I2B+2)+W3*RAM(I3B+2))+
     *         WC*(W1*RAM(I1C+2)+W2*RAM(I2C+2)+W3*RAM(I3C+2))
            Y3=WB*(W1*RAM(I1B+3)+W2*RAM(I2B+3)+W3*RAM(I3B+3))+
     *         WC*(W1*RAM(I1C+3)+W2*RAM(I2C+3)+W3*RAM(I3C+3))
            IF (((Y1-X1)**2 + (Y2-X2)**2 + (Y3-X3)**2).GT.ERRMAX) THEN
C             MTT-23
              CALL WARN('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.
            ENDIF
          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,TXTERR)
          IRAM(INDX)=NRAMT
          IRAM(NRAMT)=0
          MRAMP=NRAMT-1
          I5=1
  50      CONTINUE
            IF (OUT(I5)(1:2).EQ.'tt') THEN
C             ----------------------------------------------------------
C             Travel time - bilinear interpolation:
cc            RAM(NRAMT+I5)=
cc   *           WB*W1*RAM(I1B+(5+I5))+WB*W2*RAM(I2B+(5+I5))+
cc   *           WC*W1*RAM(I1C+(5+I5))+WC*W2*RAM(I2C+(5+I5))
cc            IF (L3D) THEN
cc              RAM(NRAMT+I5)=RAM(NRAMT+I5)+
cc   *             WB*W3*RAM(I3B+(5+I5))+
cc   *             WC*W3*RAM(I3C+(5+I5))
cc            ENDIF
C             ----------------------------------------------------------
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(I1B+6)+AWC*RAM(I1C+6)+
     *            CIBB(WB,I1B,I1C)+CIBB(WC,I1C,I1B)
              TA2=AWB*RAM(I2B+6)+AWC*RAM(I2C+6)+
     *            CIBB(WB,I2B,I2C)+CIBB(WC,I2C,I2B)
              PA(1,1)=AWB*RAM(I1B+7)+AWC*RAM(I1C+7)
              PA(2,1)=AWB*RAM(I1B+8)+AWC*RAM(I1C+8)
              PA(3,1)=AWB*RAM(I1B+9)+AWC*RAM(I1C+9)
              PA(1,2)=AWB*RAM(I2B+7)+AWC*RAM(I2C+7)
              PA(2,2)=AWB*RAM(I2B+8)+AWC*RAM(I2C+8)
              PA(3,2)=AWB*RAM(I2B+9)+AWC*RAM(I2C+9)
              AA(1,1)=WB*RAM(I1B+1)+WC*RAM(I1C+1)
              AA(2,1)=WB*RAM(I1B+2)+WC*RAM(I1C+2)
              AA(3,1)=WB*RAM(I1B+3)+WC*RAM(I1C+3)
              AA(1,2)=WB*RAM(I2B+1)+WC*RAM(I2C+1)
              AA(2,2)=WB*RAM(I2B+2)+WC*RAM(I2C+2)
              AA(3,2)=WB*RAM(I2B+3)+WC*RAM(I2C+3)
              RAM(NRAMT+I5)=
     *          CIAA(W1)*TA1+CIAA(W2)*TA2 +
     *          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)))
              IF (L3D) THEN
                TA3=AWB*RAM(I3B+6)+AWC*RAM(I3C+6)+
     *            CIBB(WB,I3B,I3C)+CIBB(WC,I3C,I3B)
                PA(1,3)=AWB*RAM(I3B+7)+AWC*RAM(I3C+7)
                PA(2,3)=AWB*RAM(I3B+8)+AWC*RAM(I3C+8)
                PA(3,3)=AWB*RAM(I3B+9)+AWC*RAM(I3C+9)
                AA(1,3)=WB*RAM(I3B+1)+WC*RAM(I3C+1)
                AA(2,3)=WB*RAM(I3B+2)+WC*RAM(I3C+2)
                AA(3,3)=WB*RAM(I3B+3)+WC*RAM(I3C+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)=RAM(NRAMT+I5) + CIAA(W3)*TA3 +
     *            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
              ENDIF
C             ----------------------------------------------------------
            ELSEIF (OUT(I5)(1:1).EQ.'p') THEN
C             Slowness vector or ray propagator matrix -
C             - bilinear interpolation:
              RAM(NRAMT+I5)=
     *           WB*W1*RAM(I1B+(5+I5))+WB*W2*RAM(I2B+(5+I5))+
     *           WC*W1*RAM(I1C+(5+I5))+WC*W2*RAM(I2C+(5+I5))
              IF (L3D) RAM(NRAMT+I5)=RAM(NRAMT+I5)+
     *           WB*W3*RAM(I3B+(5+I5))+
     *           WC*W3*RAM(I3C+(5+I5))
            ELSEIF (OUT(I5)(1:2).EQ.'hi') THEN
C             Ray history:
              IRAM(NRAMT+I5)=IRAM(I1B+(5+I5))
            ENDIF
            I5=I5+1
          IF (I5.LE.NOUT) GOTO 50
  69      CONTINUE
  70    CONTINUE
  71    CONTINUE
  80  CONTINUE
  90  CONTINUE
 100  CONTINUE
      END
C=======================================================================
C
      SUBROUTINE CILSID(I1B,I2B,I1C,I2C,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 I1B,I2B,I1C,I2C is located point X.
C The side of the ray cell defined by points I1B,I2B,I1C,I2C is
C completed with points I3B,I3C (for details see the code bellow) to
C create virtual ray cell. Local coordinates are then computed for
C point X. If (0.-DWB)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 I1B,I2B,I1C,I2C
      REAL X1,X2,X3,DWB,WB,Y1,Y2,Y3
      INTEGER ISID
C Input:
C     I1B,I2B,I1C,I2C ... 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 B denotes the bottom of the ray cell and C
C       denotes the 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(I1B+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:
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 ...........................
C Common block /MTTC/ to store the information about triangles and rays:
      INCLUDE 'mtt.inc'
C mtt.inc.
C ...........................
C Common block /CIGRD/ to store the parameters of the
C interpolation grid.
      REAL O1,O2,O3,D1,D2,D3
      REAL PROJ1,PROJ2,PROJ3
      INTEGER N1,N2,N3
      COMMON/CIGRD/O1,O2,O3,D1,D2,D3,N1,N2,N3,PROJ1,PROJ2,PROJ3
      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     PROJ1,PROJ2,PROJ3 ... Three components of the projection vector,
C                  which is used in 2-D computations to project the
C                  coordinates of gridpoints to the plane defined by
C                  rays.
C.......................................................................
      INTEGER I3B,I3C
      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
      CHARACTER*52 TXTERR
C     I3B,I3C ... 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-----------------------------------------------------------------------
      IF (NRAMP+6.GT.MRAMP) CALL CIEROR(1,TXTERR)
      IF (L3D) THEN
C       Computation of the points I3B, I3C:
        A1=(RAM(I1B+1)+RAM(I2B+1))/2.
        A2=(RAM(I1B+2)+RAM(I2B+2))/2.
        A3=(RAM(I1B+3)+RAM(I2B+3))/2.
        B1=(RAM(I1C+1)+RAM(I2C+1))/2.
        B2=(RAM(I1C+2)+RAM(I2C+2))/2.
        B3=(RAM(I1C+3)+RAM(I2C+3))/2.
        U1=RAM(I1B+1)-RAM(I2B+1)
        U2=RAM(I1B+2)-RAM(I2B+2)
        U3=RAM(I1B+3)-RAM(I2B+3)
        V1=RAM(I1C+1)-RAM(I2C+1)
        V2=RAM(I1C+2)-RAM(I2C+2)
        V3=RAM(I1C+3)-RAM(I2C+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
        I3B=NRAMP
        RAM(I3B+1)=A1+Y11
        RAM(I3B+2)=A2+Y12
        RAM(I3B+3)=A3+Y13
        I3C=NRAMP+3
        RAM(I3C+1)=B1+Y21
        RAM(I3C+2)=B2+Y22
        RAM(I3C+3)=B3+Y23
      ENDIF
      IF (L2D) THEN
C       Computation of the points I2B, I2C:
        W1=RAM(I1C+1)-RAM(I1B+1)
        W2=RAM(I1C+2)-RAM(I1B+2)
        W3=RAM(I1C+3)-RAM(I1B+3)
        Y11=PROJ2*W3-W2*PROJ3
        Y12=PROJ3*W1-W3*PROJ1
        Y13=PROJ1*W2-W1*PROJ2
        Y21=Y11
        Y22=Y12
        Y23=Y13
        I2B=NRAMP
        I3B=I2B
        RAM(I2B+1)=RAM(I1B+1)+Y11
        RAM(I2B+2)=RAM(I1B+2)+Y12
        RAM(I2B+3)=RAM(I1B+3)+Y13
        I2C=NRAMP+3
        I3C=I2C
        RAM(I2C+1)=RAM(I1C+1)+Y21
        RAM(I2C+2)=RAM(I1C+2)+Y22
        RAM(I2C+3)=RAM(I1C+3)+Y23
      ENDIF
C
C     Computation of the local coordinate wb:
      CALL CIWB(I1B,I2B,I3B,I1C,I2C,I3C,X1,X2,X3,DWB,IR,ROOT)
C
      IF(IR.LE.0) THEN
C       MTT-29
        CALL ERROR('MTT-29: No WB found in the virtual cell')
C       This error should not appear.
C       Please contact the author.
      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       MTT-30
        CALL WARN('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.
      ENDIF
C
      WBB=ROOT(J)
      WCC=1.-WBB
      IF ((WBB.LT.0.-DWB).OR.(WCC.LT.0.-DWB)) THEN
C       MTT-28
        CALL ERROR('MTT-28: Root outside the cell')
C       This error should not appear.
C       Please contact the author.
      ENDIF
      CALL CIWI(I1B,I2B,I3B,I1C,I2C,I3C,X1,X2,X3,WBB,WCC,W1,W2,W3)
      IF ((L3D.AND.W3.GE.0.).OR.(L2D.AND.W2.GE.0.)) THEN
        ISID=+1
      ELSE
        ISID=-1
      ENDIF
C
      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 ...........................
C Common block /MTTC/ to store the information about triangles and rays:
      INCLUDE 'mtt.inc'
C 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 /MTTC/ to store the information about triangles and rays:
      INCLUDE 'mtt.inc'
C 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 CIQUA(AAA,BBB,CCC,XX)
C
C----------------------------------------------------------------------
      REAL AAA,BBB,CCC,XX
C     Subroutine designed to calculate the value of quadratic function
C     given by coefficients AAA,BBB,CCC in point XX.
C......................................................................
      CIQUA=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......................................................................
      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......................................................................
      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 CIQUAR(AAA,BBB,CCC,DWB,IR,R1,R2)
C
C----------------------------------------------------------------------
C Subroutine to solve the quadratic 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.
      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 
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,CIQUA
      REAL CICUB,CIQUA
C.......................................................................
C Common block /MTTC/ to store the information about triangles and rays:
      INCLUDE 'mtt.inc'
C mtt.inc.
C.......................................................................
      REAL TRIAA0,DVEBB0,POM1,POM2,POM3,DER,RA,RB,DA,DB
      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
C           MTT-08
            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 contact the author.
          ENDIF
C
          POM1=CICUB(AA0,BB0,CC0,DD0,RB)
          POM2=CICUB(AA0,BB0,CC0,DD0,RA)
          IF (((POM1.LE.0.).AND.(POM2.LE.0.)).OR.
     *        ((POM1.GE.0.).AND.(POM2.GE.0.))) THEN
C          MTT-24
            CALL ERROR('MTT-24: values in RA and RB have the same sign')
C           This error should not appear.
C           Please contact the author.
          ENDIF
          DER=CIQUA(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.LT.0.).AND.(POM3.GT.0.)).OR.
     *        ((POM1.GT.0.).AND.(POM3.LT.0.))) THEN
            RA=ROOT(I2)
            POM2=POM3
          ELSE
            RB=ROOT(I2)
            POM1=POM3
          ENDIF
C
          IF (((POM1.LE.0.).AND.(POM2.LE.0.)).OR.
     *        ((POM1.GE.0.).AND.(POM2.GE.0.))) THEN
C          MTT-25
            CALL ERROR('MTT-25: values in RA and RB have the same sign')
C           This error should not appear.
C           Please contact the author.
          ENDIF
          POM3=POM1-POM2
          IF (POM3.EQ.0.) THEN
C           MTT-09
            CALL ERROR('MTT-09: division by zero in CICUBR.')
C           This error should not appear.
C           Please contact the author.
          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.LT.0.).AND.(POM3.GT.0.)).OR.
     *        ((POM1.GT.0.).AND.(POM3.LT.0.))) THEN
            RA=ROOT(I2)
            POM2=POM3
          ELSE
            RB=ROOT(I2)
            POM1=POM3
          ENDIF
        IF (ABS(RA-RB).GT.DWARF) GOTO 10
        IF (RA.NE.RB) THEN
          IF (POM1.EQ.0.) THEN
            ROOT(I2)=RB
          ELSEIF (POM2.EQ.0.) THEN
            ROOT(I2)=RA
          ELSEIF (POM1-POM2.NE.0.) THEN
C           Linear interpolation of the root:
            DA=RA-ROOT(I2)
            DB=RB-ROOT(I2)
            POM3=(DA*POM1-DB*POM2)/(POM1-POM2)
            ROOT(I2)=ROOT(I2)+POM3
            IF((ROOT(I2).LE.RA).OR.(ROOT(I2).GE.RB)) ROOT(I2)=(RA+RB)/2.
          ELSE
C           MTT-22
            CALL ERROR('MTT-22: Equal functional values in CICUBR.')
C           This error should not appear.
C           Please contact the author.
          END IF
        END IF
  19    CONTINUE
  20  CONTINUE
C
      END
C
C=======================================================================
C
      LOGICAL FUNCTION LDWARF(AUX)
C
C----------------------------------------------------------------------
      REAL AUX
C     Subroutine to force the computer to take the value of DWARF from
C     the memory.
      IF (AUX.NE.1.) THEN
        LDWARF=.TRUE.
      ELSE
        LDWARF=.FALSE.
      ENDIF
      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
      LOGICAL FUNCTION CIL2P(I1,I2)
C
C----------------------------------------------------------------------
      INTEGER I1,I2
C     Subroutine to compare coordinates of two points from array RAM.
C
C Input: I1,I2 ... Addresses of the two points in array RAM.
C Output: CIL2P... .TRUE. if the coordinates of the points are different
C                  .FALSE. if they are the same.
C ...........................
C Common block /MTTC/ to store the information about triangles and rays:
      INCLUDE 'mtt.inc'
C mtt.inc.
C-----------------------------------------------------------------------
      IF ((RAM(I1+1).NE.RAM(I2+1)).OR.(RAM(I1+2).NE.RAM(I2+2))
     *       .OR.(RAM(I1+3).NE.RAM(I2+3))) THEN
        CIL2P=.TRUE.
      ELSE
        CIL2P=.FALSE.
      ENDIF
      END
C
C
C=======================================================================
C
      SUBROUTINE CIGRID
C
C----------------------------------------------------------------------
C Subroutine to read in the file with parameters of an interpolation
C grid. The read quantities are then written into common block CIGRD.
C No input:
C No output.
C
C ...........................
C Common block /MTTC/ to store the information about triangles and rays:
      INCLUDE 'mtt.inc'
C mtt.inc.
C ...........................
C Common block /CIGRD/ to store the parameters of the
C interpolation grid.
      REAL O1,O2,O3,D1,D2,D3
      REAL PROJ1,PROJ2,PROJ3
      INTEGER N1,N2,N3
      COMMON/CIGRD/O1,O2,O3,D1,D2,D3,N1,N2,N3,PROJ1,PROJ2,PROJ3
      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     PROJ1,PROJ2,PROJ3 ... Three components of the projection vector,
C                  which is used in 2-D computations to project the
C                  coordinates of gridpoints to the plane defined by
C                  rays.
C.......................................................................
      INTEGER I1
      CHARACTER TXTERR
C-----------------------------------------------------------------------
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     Recalling the data specifying projection vector, which is used
C     in 2-D calculations for projection of the gridpoint coordinates
C     to the plane defined by rays.
      CALL RSEP3R('PROJ1',PROJ1,0.)
      CALL RSEP3R('PROJ2',PROJ2,0.)
      CALL RSEP3R('PROJ3',PROJ3,0.)
C
      IF (D1.LT.0.) THEN
        O1=O1+(N1-1)*D1
        D1=-D1
      ENDIF
      IF (D2.LT.0.) THEN
        O2=O2+(N2-1)*D2
        D2=-D2
      ENDIF
      IF (D3.LT.0.) THEN
        O3=O3+(N3-1)*D3
        D3=-D3
      ENDIF
      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,TXTERR)
      NRAMT=MRAMP+1
      DO 5, I1=NRAMT,MRAM
        IRAM(I1)=0
   5  CONTINUE
      RETURN
C
  10  CONTINUE
C     MTT-10
      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,TXTERR)
C
C-----------------------------------------------------------------------
      INTEGER IERR
      CHARACTER*52 TXTERR
C
C Subroutine designed to print error messages of different
C subroutines of the program MTT using subroutine ERROR.
C
C Input:
C     IERR ... Index of the error.
C     TXTERR ... Description of the error.
C No output.
C-----------------------------------------------------------------------
C
      IF     (IERR.EQ.01) THEN
C       MTT-01
        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       ram.inc.
      ELSEIF (IERR.EQ.27) THEN
C       MTT-27
        CALL ERROR('MTT-27: Small array FOUT for the output files.')
C       This error may be caused by wrongly prescribed
C        input data FOUT and OUT, or by
C       small parameter MOUT defined in the file
C       mtt.inc.
      ELSEIF (IERR.EQ.31) THEN
C       MTT-31
        CALL ERROR(TXTERR)
C       This error may be caused by too small dimension of array
C       RAM. Try to enlarge dimension MRAM in include file
C       ram.inc. Number of
C       needed aditional units should have been written on
C       standard output by subroutine ERROR.
      ENDIF
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'ap.for'
C     ap.for
      INCLUDE 'forms.for'
C     forms.for
      INCLUDE 'length.for'
C     length.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'indexx.for'
C     indexx.for
C
C=======================================================================
C