C 
C Program MTT to calculate Multivalued Travel Times on a 3-D grid
C of points
C
C Version: 5.90
C Date: 2005, June 19
C                                                  
C Coded by Petr Bulant
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     E-mail: bulant@seis.karlov.mff.cuni.cz
C
C----------------------------------------------------------------------
C
C The post processing program to the CRT program performs interpolation
C of travel times and other quantities to the gridpoints of arbitrary
C regular rectangular 3-D grid of points. The interpolation inside ray
C tubes formed by vertices of homogeneous triangles created during 3-D
C two-point ray tracing is realized using the controlled initial-value
C ray tracing algorithm. The interpolation is performed within all of
C the ray tubes stored in the specified CRT output file 'CRT-T'.
C
C The description of the quantities which may be interpolated can be
C found below. The user may include many other
C quantities by editting the file mttq.for.
C
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 Description of data files:
C
C Input data read from the standard input device (*):
C     The data are read by the list directed input (free format) and
C     consist of a single string 'SEP':
C     'SEP'...String in apostrophes containing the name of the input
C             SEP parameter or history file with the input data.
C     No default, 'SEP' must be specified and cannot be blank.
C
C                                                     
C Input data file 'SEP':
C     File 'SEP' has the form of the SEP
C     parameter file.  The parameters, which do not differ from their
C     defaults, need not be specified in file 'SEP'.
C 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.
C             For general description of file CRTOUT refer to
C             file writ.for.
C             Description specific to this program:
C               Just the first set of CRT output files is read from
C               file CRTOUT.
C               Only files 'CRT-R',
C               'CRT-I' and
C               'CRT-T' are read by MTT.
C               File 'CRT-R' must contain all rays traced by CRT, not
C               only two-point rays.
C             Default: CRTOUT=' ' which means 'CRT-R'='r01.out',
C             'CRT-I'='r01i.out', 'CRT-T'='t01.out'
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             Defaults: 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             For NUM=' ' the file is not generated.
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             For MTT=' ' the file is not generated.
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 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     GBW11='string' ... Name of the output file containing the grid
C             values of the sum of squares of Gaussian beam widths
C             corresponding to GBY11.
C             The sum of squares of Gaussian beam widths is defined here
C             as the trace of the inverse imaginary part of the matrix
C             of second derivatives of the travel time of the Gaussian
C             beam.  The sum of squares of Gaussian beam widths depends
C             on parameters GBEij, GBR11, GBR12, GBR22, GBY11 and GBY22
C             described below.  It may be expressed as sum W11+W22,
C             where W11 is dependent on GBY11 but independent on GBY22,
C             and W22 is dependent on GBY22 but independent on GBY11.
C             Values of W11 are written to the file given by parameter
C             GBW11.
C             Note that positive GBY11 must be specified in order to
C             generate file GBW11.
C             Default:  GBW11=' ' (no file generated)
C     GBW1='string' ... Name of the output file containing the grid
C             values of SQRT(W11).  See the description of GBW11.
C             Note that positive GBY11 must be specified in order to
C             generate file GBW1.
C             Default:  GBW1=' ' (no file generated)
C     GBW22='string' ... Name of the output file containing the grid
C             values of the sum W22 of squares of Gaussian beam widths
C             corresponding to GBY22.  See the description of GBW11.
C             Note that positive GBY22 must be specified in order to
C             generate file GBW22.
C             Default:  GBW22=' ' (no file generated)
C     GBW2='string' ... Name of the output file containing the grid
C             values of SQRT(W22).  See the description of GBW11.
C             Note that positive GBY22 must be specified in order to
C             generate file GBW2.
C             Default:  GBW2=' ' (no file generated)
C     AMP='string' ... Name of the output file with the grid values of
C             the norm of the vectorial amplitude of the Green function.
C             Default:  AMP=' ' (no file generated)
C     AMPKI='string' ... Name of the output file with the grid values of
C             the amplitude modified for use in the Kirchoff integral.
C             Default:  AMPKI=' ' (no file generated)
C                                                   
C Order of the multivalued interpolated quantities:
C     MTTSORT='string'... Quantity according which the multivalued
C             interpolated quantities are sorted.  The value of the
C             'string' should equal the name of one of the above SEP
C             parameters specifying the names of output formatted files
C             with interpolated quantities, typed in uppercase.
C             Default:  MTTSORT='MTT' (sorting according to travel
C                                     times)
C     MTTORDER=real... Order of sorting:
C             MTTORDER=1 for sorting in the ascending order,
C             MTTORDER=-1 for sorting in the descending order.
C             Default:  MTTORDER=1.
C                                                   
C Numerical parameters specifying the interpolated quantities:
C     GBE11=real, GBE21=real, GBE31=real... Components of the first of
C             two vectors specifying the plane in which the initial
C             second derivatives of the complex-valued travel times of
C             Gaussian beams are given.  The vectors are basis vectors
C             of the linear coordinate system for the specification
C             of the second derivatives.
C             Defaults:  GBE11=0., GBE21=1., GBE31=0.
C     GBE12=real, GBE22=real, GBE32=real... Components of the second of
C             two vectors specifying the plane and coordinates for the
C             initial second derivatives of the complex-valued travel
C             times of Gaussian beams.
C             Defaults:  GBE12=1., GBE22=0., GBE32=0.
C     GBR11=real, GBR12=real, GBR22=real... Real part of the second
C             derivatives of the complex-valued travel times of Gaussian
C             beams in the directions of vectors GBEi1 and GBEi2.
C             Defaults:  GBR11=0., GBR12=0., GBR22=0.
C     GBY11=real, GBY22=real... Imaginary part of the second derivatives
C             of the complex-valued travel times of Gaussian beams in
C             the directions of vectors GBEi1 and GBEi2.
C             Vectors GBEi1 and GBEi2 are assumed to be specified in
C             such a way that mixed second derivatives GBY12=0.
C             Defaults:  GBY11=0., GBY22=0.
C     FGBE11=string, FGBE21=string, FGBE31=string, FGBE12=string,
C     FGBE22=string, FGBE32=string, FGBR11=string, FGBR12=string,
C     FGBR22=string, FGBY11=string, FGBY22=string... Filenames
C             corresponding to the above quantities, if the quantities
C             depend on two-way travel time and ray parameters and
C             are specified on an optional ray-parameter grid of
C             NTIME*NPAR1*NPAR2*NWAVES*NCRT values, see the grid
C             parameters below.
C             Undefined grid values are not allowed in this version.
C             If the filename is blank, a constant value given by the
C             corresponding real-valued parameter is used.
C             Defaults:  All blank.
C     FTIME=string... Name of the file containing two-way travel time,
C             gridded in dependence on travel times (inner loop) and ray
C             parameters.  NTT*NPAR1*NPAR2*NWAVES*NCRT values.
C             Undefined grid values are not allowed in this version.
C             If FTIME is blank, two-way travel time is deemed to
C             coincide with the travel time and parameters NTT, OTT, DTT
C             need not be specified.
C             Default:  FTIME=' '
C     AMPMAX=real... Maximum value of the amplitude AMP of the Green
C             function to be stored in the output file.
C             If the amplitude of the Green function is greater, it is
C             replaced by AMPMAX.
C             Default:  AMPMAX=999999.
C     AMP2D1=real, AMP2D2=real, AMP2D3=real... Unit vector perpendicular
C             to a 2-D section in the 3-D model.
C             Writing the amplitudes of the 2-D Green function instead
C             of those of the 3-D Green function.
C     AMPKI1=real, AMPKI2=real, AMPKI3=real... Unit vector perpendicular
C             to the surface of integration.
C             Modification of amplitudes for Kirchhoff integral
C                                                   
C Numerical parameters specifying optional ray-parameter grid:
C             Number of gridpoints is NTT*NPAR1*NPAR2*NWAVES*NCRT or
C             NTIME*NPAR1*NPAR2*NWAVES*NCRT, inner loop corresponding to
C             NTT or NTIME, outer loop to NCRT.
C     NTT=positive integer... Number of travel time samples (inner loop)
C             in the file given by parameter FTIME.
C             Not used if FTIME=' '.
C             Default:  NTT=1
C     OTT=real... Travel time of the first sample in the file given by
C             parameter FTIME.
C             Not used if FTIME=' '.
C             Default:  OTT=0.
C     DTT=real... Travel-time increment in the file given by parameter
C             FTIME.
C             Not used if FTIME=' '.
C             Default:  DTT=1.
C     NTIME=positive integer... Number of two-way travel times (inner
C             loop) in the files given by parameters FGBE11, FGBE21,
C             FGBE31, FGBE12, FGBE22, FGBE32, FGBR11, FGBR12, FGBR22,
C             FGBY11 and FGBY22.  If all these filenames are blank,
C             parameters NTIME, OTIME, DTIME, NPAR1, OPAR1, DPAR1,
C             NPAR2, OPAR2, DPAR2, NWAVES, NCRT and ICRT are not used.
C             Default:  NTIME=1
C     OTIME=real... Two-way travel time of the first sample.
C             Default:  OTIME=0.
C     DTIME=real... Two-way travel-time increment.
C             Default:  DTIME=1.
C     NPAR1=positive integer... Number of gridpoints corresponding to
C             the first ray take-off parameter (inner but one loop) in
C             the files given by parameters FTIME, FGBE11, FGBE21,
C             FGBE31, FGBE12, FGBE22, FGBE32, FGBR11, FGBR12, FGBR22,
C             FGBY11 and FGBY22.
C             Default:  NPAR1=1
C     OPAR1=real... The first ray take-off parameter corresponding to
C             the first gridpoint.
C             Default:  OPAR1=0.
C     DPAR1=real... Increment of the first ray take-off parameter.
C             Default:  DPAR1=1.
C     NPAR2=positive integer... Number of gridpoints corresponding to
C             the second ray take-off parameter in the files given by
C             parameters FTIME, FGBE11, ..., FGBY22.
C             Default:  NPAR1=1
C     OPAR2=real... The second ray take-off parameter corresponding to
C             the first gridpoint.
C             Default:  OPAR2=0.
C     DPAR2=real... Increment of the second ray take-off parameter.
C             Default:  DPAR2=1.
C     NWAVES=positive integer... Maximum number of elementary waves
C             or sources per a single execution of program 'crt.for'.
C             The waves are indexed 1,2,...,NWAVES.
C             Default:  NWAVES=1
C     NCRT=positive integer... Total number of executions of program
C             'crt.for'.
C             Default:  NCRT=1
C     ICRT=positive integer... Sequential index of the execution of
C             program 'crt.for'.
C             Default:  ICRT=1
C Optional parameters specifying the form of the real quantities
C written in the output formatted files:
C     MINDIG,MAXDIG=positive integers ... See the description in file
C             forms.for.
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     For general description of the files with gridded data refer to
C     file forms.htm.
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 ascending
C     order according to the travel time. The loop over gridpoints
C     is the same as in file NUM.
C     For general description of the files with multivalued gridded data
C     refer to file forms.htm.
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      EXTERNAL CIEROR,CIGRID,CIREAD,CIERAS,CITUBE,CIINTP,CIWB,CIWI,
     *CIQD,CIQW,CIQRI,CIQRA,CIQI,LDWARF,CILSID,CIPPP,CILSIG,CIL2P,
     *CIDET3,CISUBD,CIQUAR,CICUBR,CIQUA,CICUB,CIAA,CIBB,
     *MTTQ,
     *ERROR,WARN,RSEP1,RSEP3T,RSEP3R,RSEP3I,WARRAY,WARRAI,INDEXX,AP00
      LOGICAL LDWARF
C     CIEROR,CIGRID,CIREAD,CIERAS,CITUBE,CIINTP,CIWB,CIWI,
C     CIQD,CIQW,CIQRI,CIQRA,CIQI,LDWARF,CILSID,CIPPP,CILSIG,CIL2P,
C     CIDET3,CISUBD,CIQUAR,CICUBR,CIQUA,CICUB,CIAA,CIBB ... This file.
C     MTTQ ... File mttq.for.
C     ERROR,WARN ... File error.for.
C     RSEP1,RSEP3I,RSEP3T,RSEP3R ...
C               File sep.for.
C     WARRAI,WARRAI ... File forms.for.
C     INDEXX ... File indexx.for.
C     AP00 ... File ap.for.
C
C Common block /MTTC/ to store the information about triangles and rays:
      INCLUDE 'mtt.inc'
C mtt.inc.
C
C.......................................................................
C
C     Auxiliary storage locations:
      INTEGER IRAY1,IRAY2,IRAY3
      INTEGER I1,I2,I3,I4,INDX,INDX1,INDX2,IT1,IIT1,JWAVES
      INTEGER I
      INTEGER LU0,LU1,LU2,LU3,LU4,LU5
      PARAMETER (LU0=10,LU1=1,LU2=2,LU3=3,LU4=4,LU5=5)
      REAL ORDER,AUX
      CHARACTER*80 FILSEP,FILCRT,FILE1,FILE2,FILE3,CH
      CHARACTER*52 TXTERR
      CHARACTER*6 SORT
      LOGICAL LENDT
C
C-----------------------------------------------------------------------
C
C     Reading a name of the file with the input data:
      FILSEP=' '
      WRITE(*,'(A)') '+MTT: Enter input filename: '
      READ(*,*) FILSEP
      IF (FILSEP.EQ.' ') THEN
C       MTT-02
        CALL ERROR('MTT-02: No input file specified')
C       Input file in the form of the SEP (Stanford Exploration Project)
C       parameter or history file must be specified.
C       There is no default filename.
      ENDIF
      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 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 filenames of the output files, recording
C     which quantities are to be written into the files:
      CALL CIQD
C
C     Reading 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 file with computed triangles,
C     sorting the rays in triangles:
      WRITE(*,'(A)') '+MTT: Reading the file with triangles.'
      JWAVES=0
      LENDT=.FALSE.
      READ(LU3,*,END=20) IRAY1,IRAY2,IRAY3
   8  CONTINUE
      NT=0
      NRAMP=0
      IRAYMI=999999
      IRAYMA=0
  10  CONTINUE
        IRAY1=0
        IRAY2=0
        IRAY3=0
        READ(LU3,*,END=20) IRAY1,IRAY2,IRAY3
        IF (IRAY1.EQ.0) GOTO 22
        IF ((IRAY1.LE.0).OR.(IRAY2.LE.0).OR.(IRAY3.LT.0)) THEN
C         MTT-36
          CALL ERROR('MTT-36: Disorder in file ''CRT-T''.')
C         The 'triangles' in the input file 'CRT-T' should be written
C         as three nonnegative integers (indices of rays), two of them
C         being positive (the third may equal zero in 2-D computations).
        ENDIF
        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
      LENDT=.TRUE.
      CLOSE(LU3)
  22  CONTINUE
      JWAVES=JWAVES+1
      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 input 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(*,'(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"):
      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"):
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,1I2,A,1I4,A)')
     *    '+MTT: Interpolating wave ',JWAVES,'. ',
     *    I2,'% of ray tubes completed.'
          I1=I2
        ENDIF
C
C       If necessary, reading 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(MAXRAM/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.
C     End of the loop for all the triangles of the elementary wave.
      if (.not.lendt) goto 8
      CLOSE(LU1)
      CLOSE(LU2)
C
C
C     Sorting the results according to the travel times:
      CALL RSEP3T('MTTSORT',SORT,'MTT')
      CALL RSEP3R('MTTORDER',ORDER,1.)
      DO 340, I4=1,NOUT
        IF (CHOUT(I4).EQ.SORT) 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=MAXRAM-INDX
            INDX1=IRAM(INDX)
            IF (INDX1.EQ.0) GOTO 306
  302       CONTINUE
              INDX2=IRAM(INDX1)
              IF (INDX2.EQ.0) GOTO 306
              IF (ORDER*RAM(INDX2+I4).LT.ORDER*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:
      CALL CIQW(LU5)
C
      WRITE(*,'(A)') '+MTT: Done.                           '
      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.
C     LU2 ... Number of logical unit corresponding to the file with
C             the quantities at the initial points of rays,
C             corresponding to file LU1.
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,IWAVE0,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
        IWAVE0=IWAVE
        IF ((IRAY0.EQ.0).AND.(IWAVE0.EQ.0)) THEN
C         Reading the first point on a ray of the first wave:
          CALL AP00(0,LU1,LU2)
          IF (IWAVE.LT.1) GOTO 20
        ELSEIF (IWAVE.EQ.IWAVE0) THEN
C         Reading the next point on a ray of the actual wave:
          CALL AP00(0,LU1,LU2)
          IF (IWAVE.NE.IWAVE0) GOTO 20
        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).AND.(IRAY.EQ.IRAYMI)) THEN
C             MTT-18
              CALL WARN
     *           ('MTT-18: a ray with history equal to 0 was observed.')
C              All the rays are probably of the same history 0. Only the
C             initial-value ray tracing was performed.  Width and shape
C             of ray tubes were not controlled.  The interpolation
C             in such a case makes sense only in smooth models.
C              Check the value of the parameter IPOINT in CRT input data
C             RPAR (4).
            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:
            CALL CIQRI
          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:
          CALL CIQRA
        ENDIF
        IRAM(IRAY-ORAYA)=NRAMP
        IF ((IPT.LE.1).AND.(IRAY.GT.IRAM(IT1)).AND.(IRAY.NE.IRAYMA))
     *     RETURN
      GOTO 10
C
  20  CONTINUE
C     End of rays:
      IF (IRAY0.LT.IRAYMA) THEN
C       Some rays at the end of the CRT output file missing:
        DO 22, I1=IRAY0+1,IRAYMA
          IF (I1.GE.IRAYMI) THEN
            IRAM(I1-ORAYA)=NRAMP
          ENDIF
  22    CONTINUE
      ENDIF
      RETURN
      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 recommend 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.(L3D.AND.(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))) 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).AND.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)
          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 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.......................................................................
      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=0.
      W2=0.
      W3=0.
      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=WB*(RAM(I1B+1)-RAM(I2B+1)) + WC*(RAM(I1C+1)-RAM(I2C+1))
        AA2=WB*(RAM(I1B+2)-RAM(I2B+2)) + WC*(RAM(I1C+2)-RAM(I2C+2))
        AA3=WB*(RAM(I1B+3)-RAM(I2B+3)) + WC*(RAM(I1C+3)-RAM(I2C+3))
        EE=E1*E1+E2*E2+E3*E3
        AUX=(AA1*AA1+AA2*AA2+AA3*AA3)*EE-(AA1*E1+AA2*E2+AA3*E3)**2
        IF (AUX.NE.0.) THEN
          W1=(XA1*AA1+XA2*AA2+XA3*AA3)*EE -
     *       (XA1*E1+XA2*E2+XA3*E3)*(AA1*E1+AA2*E2+AA3*E3)
          W1=W1/AUX
          W2=1. - W1
          IF (ABS(W1-1.).LT.DWARF) THEN
            XA1=X1-A11
            XA2=X2-A21
            XA3=X3-A31
            AA1=-AA1
            AA2=-AA2
            AA3=-AA3
            AUX=(AA1*AA1+AA2*AA2+AA3*AA3)*EE-(AA1*E1+AA2*E2+AA3*E3)**2
            IF (AUX.NE.0.) THEN
              W2=(XA1*AA1+XA2*AA2+XA3*AA3)*EE -
     *           (XA1*E1+XA2*E2+XA3*E3)*(AA1*E1+AA2*E2+AA3*E3)
              W2=W2/AUX
              W1=1. - W2
            ENDIF
          ENDIF
        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=WB*(RAM(I1B+1)-RAM(I3B+1)) + WC*(RAM(I1C+1)-RAM(I3C+1))
        A21=WB*(RAM(I1B+2)-RAM(I3B+2)) + WC*(RAM(I1C+2)-RAM(I3C+2))
        A31=WB*(RAM(I1B+3)-RAM(I3B+3)) + WC*(RAM(I1C+3)-RAM(I3C+3))
        A12=WB*(RAM(I2B+1)-RAM(I3B+1)) + WC*(RAM(I2C+1)-RAM(I3C+1))
        A22=WB*(RAM(I2B+2)-RAM(I3B+2)) + WC*(RAM(I2C+2)-RAM(I3C+2))
        A32=WB*(RAM(I2B+3)-RAM(I3B+3)) + WC*(RAM(I2C+3)-RAM(I3C+3))
        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.......................................................................
      INTEGER K1MI,K1MA,K2MI,K2MA,K3MI,K3MA
      REAL AUX
      INTEGER I1,I2,I3,I4
      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     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
C            average 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
        INDX=MAXRAM-INDX+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.(WB.GT.1.+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.EQ.0.).AND.(W2.EQ.0.).AND.(W3.EQ.0.))
C           For the considered WB the gridpoint is outside the ray cell.
C           Continuing with the next WB:
     *      GOTO 69
          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).LE.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).LE.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).LE.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).LE.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).LE.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:
  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
          CALL CIQI(WB,WC,W1,W2,W3,I1B,I2B,I3B,I1C,I2C,I3C,X1,X2,X3)
  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               0  if the position with respect to the side
C                  was not found,
C              -1  otherwise.
C
C ...........................
C Common block /MTTC/ to store the information about triangles and rays:
      INCLUDE 'mtt.inc'
C mtt.inc.
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       No WB found in the virtual cell:
        ISID=0
        RETURN
      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.(WBB.GT.1.+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 ((W1.EQ.0.).AND.(W2.EQ.0.).AND.(W3.EQ.0.)) THEN
        ISID=0
      ELSEIF ((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
      RETURN
      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
      RETURN
      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
      RETURN
      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))
      RETURN
      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
      RETURN
      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))
      RETURN
      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 CICUBR(AA0,BB0,CC0,DD0,RR,IR,ROOT)
C
C----------------------------------------------------------------------
C Subroutine for numerical solution of the cubic equation
C on the interval .
      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)/(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
      RETURN
      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
      RETURN
      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 MTTC.
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.......................................................................
      INTEGER I1
      CHARACTER*52 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=MAXRAM-N1*N2*N3
      IF (MRAMP.LE.1) CALL CIEROR(1,TXTERR)
      NRAMT=MRAMP+1
      DO 5, I1=NRAMT,MAXRAM
        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
      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.......................................................................
      INTEGER I4,I5
      REAL YY(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
      YY(1,1)=RAM(I1B+1)-RAM(I1C+1)
      YY(2,1)=RAM(I1B+2)-RAM(I1C+2)
      YY(3,1)=RAM(I1B+3)-RAM(I1C+3)
      YY(1,2)=RAM(I2B+1)-RAM(I2C+1)
      YY(2,2)=RAM(I2B+2)-RAM(I2C+2)
      YY(3,2)=RAM(I2B+3)-RAM(I2C+3)
      IF (L3D) THEN
        YY(1,3)=RAM(I3B+1)-RAM(I3C+1)
        YY(2,3)=RAM(I3B+2)-RAM(I3C+2)
        YY(3,3)=RAM(I3B+3)-RAM(I3C+3)
      ELSE
        YY(1,3)=0.
        YY(2,3)=0.
        YY(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(YY)
      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,YY)
      CCC=CISUBD(YY,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((YY(1,1).EQ.0.).AND.(YY(2,1).EQ.0.).AND.(YY(3,1).EQ.0.).AND.
     *     (YY(1,2).EQ.0.).AND.(YY(2,2).EQ.0.).AND.(YY(3,2).EQ.0.).AND.
     *     (YY(1,3).EQ.0.).AND.(YY(2,3).EQ.0.).AND.(YY(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
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, R1mtt.inc.
C ...........................
C Common block /POINTC/ to store the results of complete ray tracing:
      INCLUDE 'pointc.inc'
C pointc.inc.
C None of the storage locations of the common block are altered.
C.......................................................................
C     Auxiliary storage locations for all the entries:
      INTEGER LU5
      INTEGER I1B,I2B,I3B,I1C,I2C,I3C
      REAL WB,WC,W1,W2,W3,X1,X2,X3
      INTEGER I1,I2,I3,I4,I5,I6,I7,I8,INDX
      REAL AWB,AWC,TA1,TA2,TA3,PA(3,3),AA(3,3),KSI
      CHARACTER*52 TXTERR
C.......................................................................
C
C     Reading filenames of the output files, recording
C     which quantities are to be written into the files:
      CALL RSEP3T('NUM',FOUT(0),'mtt-num.out')
      CALL RSEP3T('MTT',FOUT(1),'mtt-tt.out')
      IF (FOUT(1).NE.' ') THEN
        CHOUT(1)='MTT'
        CHOUT(2)=' '
        CHOUT(3)=' '
        CHOUT(4)=' '
        NOUT=4
        CALL RSEP3T('MP1',FOUT(2),' ')
        IF (FOUT(2).NE.' ') CHOUT(2)='MP1'
        CALL RSEP3T('MP2',FOUT(3),' ')
        IF (FOUT(3).NE.' ') CHOUT(3)='MP2'
        CALL RSEP3T('MP3',FOUT(4),' ')
        IF (FOUT(4).NE.' ') CHOUT(4)='MP3'
      ELSE
        NOUT=0
        CALL RSEP3T('MP1',FOUT(NOUT+1),' ')
        IF (FOUT(NOUT+1).NE.' ') THEN
          CHOUT(NOUT+1)='MP1'
          NOUT=NOUT+1
        ENDIF
        CALL RSEP3T('MP2',FOUT(NOUT+1),' ')
        IF (FOUT(NOUT+1).NE.' ') THEN
          CHOUT(NOUT+1)='MP2'
          NOUT=NOUT+1
        ENDIF
        CALL RSEP3T('MP3',FOUT(NOUT+1),' ')
        IF (FOUT(NOUT+1).NE.' ') THEN
          CHOUT(NOUT+1)='MP3'
          NOUT=NOUT+1
        ENDIF
      ENDIF
      CALL RSEP3T('MHIST',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        CHOUT(NOUT+1)=  'MHIST'
        NOUT=NOUT+1
      ENDIF
C
C     Reading the list QNAMES of the possible interpolated quantities:
      CALL MTTQ
C     Recording the quantities to be interpolated:
      DO 5, I1=1,MOUT
        IF (QNAMES(I1).EQ.' ') GOTO 6
        CALL RSEP3T(QNAMES(I1),FOUT(NOUT+1),' ')
        IF (FOUT(NOUT+1).NE.' ') THEN
          NOUT=NOUT+1
          IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR)
          CHOUT(NOUT)=QNAMES(I1)
C         Reading the input data for the interpolated quantity:
          CALL MTTQD(QNAMES(I1))
        ENDIF
  5   CONTINUE
  6   CONTINUE
C
      NQ=5+NOUT
      RETURN
C
C=======================================================================
C
      ENTRY CIQRI
C
C.......................................................................
C
C Entry designed to read the output of the CRT at an initial point of
C a ray.
C
      DO 10, I1=1,NOUT
        IF (CHOUT(I1).EQ.'MTT') 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 (CHOUT(I1).EQ.'MP1') THEN
C         First component of the slowness vector:
          RAM(NRAMP+(5+I1))=YI(6)
        ELSEIF (CHOUT(I1).EQ.'MP2') THEN
C         Second component of the slowness vector:
          RAM(NRAMP+(5+I1))=YI(7)
        ELSEIF (CHOUT(I1).EQ.'MP3') THEN
C         Third component of the slowness vector:
          RAM(NRAMP+(5+I1))=YI(8)
        ELSEIF (CHOUT(I1).EQ.'MHIST') THEN
C         Ray history:
          IRAM(NRAMP+(5+I1))=ISHEET
        ELSEIF (CHOUT(I1).EQ.' ') THEN
C         Nothing to do:
          CONTINUE
        ELSE
C         User-defined quantity:
          RAM(NRAMP+(5+I1))=0.
          CALL MTTQRI(CHOUT(I1),RAM(NRAMP+(5+I1)))
        ENDIF
  10  CONTINUE
      NRAMP=NRAMP+NQ
      RETURN
C=======================================================================
C
      ENTRY CIQRA
C
C.......................................................................
C
C Entry designed to read the output of the CRT at other than initial
C point of a ray.
C
      DO 20, I1=1,NOUT
        IF (CHOUT(I1).EQ.'MTT') 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 (CHOUT(I1).EQ.'MHIST') THEN
C         Ray history:
          IRAM(NRAMP+(5+I1))=ISHEET
        ELSEIF (CHOUT(I1).EQ.'MP1') THEN
C         First component of the slowness vector:
          RAM(NRAMP+(5+I1))=Y(6)
        ELSEIF (CHOUT(I1).EQ.'MP2') THEN
C         Second component of the slowness vector:
          RAM(NRAMP+(5+I1))=Y(7)
        ELSEIF (CHOUT(I1).EQ.'MP3') THEN
C         Third component of the slowness vector:
          RAM(NRAMP+(5+I1))=Y(8)
        ELSEIF (CHOUT(I1).EQ.' ') THEN
C         Nothing to do:
          CONTINUE
        ELSE
C         User-defined quantity:
          RAM(NRAMP+(5+I1))=0.
          CALL MTTQRA(CHOUT(I1),RAM(NRAMP+(5+I1)))
        ENDIF
  20  CONTINUE
      NRAMP=NRAMP+NQ
      RETURN
C
C=======================================================================
C
      ENTRY CIQI(WB,WC,W1,W2,W3,I1B,I2B,I3B,I1C,I2C,I3C,X1,X2,X3)
C
C.......................................................................
C
C Entry designed to perform the interpolation.
      DO 30, I5=1,NOUT
        IF (CHOUT(I5).EQ.'MTT') 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 24, I6=1,3
              DO 23, I7=1,3
                DO 22, I8=1,3
                  KSI=KSI+PA(I8,I7)*(AA(I8,I6)-AA(I8,I7))
  22            CONTINUE
  23          CONTINUE
  24        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 (CHOUT(I5).EQ.'MHIST') THEN
C         Ray history - no interpolation:
          IRAM(NRAMT+I5)=IRAM(I1B+(5+I5))
        ELSEIF (CHOUT(I5).EQ.' ') THEN
C         Nothing to do:
          CONTINUE
        ELSE
C         Slowness vector, ray propagator matrix or other real quantity
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))
        ENDIF
  30  CONTINUE
      RETURN
C
C=======================================================================
C
      ENTRY CIQW(LU5)
C
C.......................................................................
C
C Entry designed to write the output files.
C
C     Writing numbers of computed arrivals:
      I5=N1*N2*N3
      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
      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=MAXRAM-INDX
  105   CONTINUE
C         Check for the consistentcy of the results
          IF (INDX.LE.MRAMP.OR.INDX.GT.MAXRAM) 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
      IF (FOUT(0).NE.' ') THEN
        CALL WARRAI(LU5,FOUT(0),'FORMATTED',.FALSE.,0,.FALSE.,0,
     *                                                    N1*N2*N3,IRAM)
      ENDIF
C
C
C     Writing interpolated quantities:
      IF (NOUT.GE.1) THEN
        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
  31    CONTINUE
          NRAMP=0
          DO 37, I3=0,N3-1
          DO 36, I2=0,N2-1
          DO 35, I1=0,N1-1
            INDX=I3*N1*N2+I2*N1+I1
            INDX=MAXRAM-INDX
            IF (IRAM(INDX).EQ.0) GOTO 33
            INDX=IRAM(INDX)
  32        CONTINUE
C             Check for the consistentcy of the results
              IF (INDX.LE.MRAMP.OR.INDX.GT.MAXRAM) 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 (CHOUT(I4).EQ.'MHIST') THEN
C               Ray history:
                IRAM(NRAMP)=IRAM(INDX+I4)
              ELSEIF (CHOUT(I4).NE.' ') THEN
C               Travel time or other real quantity:
                RAM(NRAMP)=RAM(INDX+I4)
              ENDIF
              INDX=IRAM(INDX)
            IF (INDX.NE.0) GOTO 32
  33        CONTINUE
  35      CONTINUE
  36      CONTINUE
  37      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 (CHOUT(I4).EQ.'MHIST') THEN
C           Ray history:
            CALL WARRAI(LU5,FOUT(I4),'FORMATTED',.FALSE.,0,.FALSE.,0,
     *                                                   I5,IRAM)
          ELSEIF (CHOUT(I4).NE.' ') THEN
C           Travel time or other real quantity:
            CALL WARRAY(LU5,FOUT(I4),'FORMATTED',.FALSE.,0.,.FALSE.,0.,
     *                                                   I5,RAM)
          ENDIF
          I4=I4+1
        IF (I4.LE.NOUT) GOTO 31
      ENDIF
C
      RETURN
      END
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 or CHOUT.')
C       This error may be caused by wrongly prescribed
C       input data, 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 'mttq.for'
C     mttq.for
      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 of Numerical Recipes
C
C=======================================================================
C