C 
C Program MTT to calculate Multivalued Travel Times either on a 3-D grid
C of points, or in specified individual points.
C
C Version: 7.30
C Date: 2016, February 10
C                                                  
C Coded by Petr Bulant
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     http://sw3d.cz/staff/bulant.htm
C
C----------------------------------------------------------------------
C
C The post processing program to the CRT program performs interpolation
C of travel times and other quantities either to the gridpoints of
C arbitrary regular rectangular 3-D grid of points, or to the individual
C points specified by the input file 'PTS'. 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
C 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
C of the program 'crt.for'.
C 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
C PRM0(3) of the program
C 'crt.for'. Such a behaviour is useful for two-point ray tracing but
C has 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)
C of the program 'crt.for' may thus considerably influence the results
C of the 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
C     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 file
C             writ.for.
C             Description specific to this program:
C               Just the first set of CRT output files is read from
C               file CRTOUT.  Only files
C               'CRT-R',
C               'CRT-I' and
C               'CRT-T'
C               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 for the interpolation to the target 3-D grid:
C     The interpolation to the 3-D grid of points is performed if the
C     name of the output file MTTPTS is not specified (default).
C   Data specifying the parameters of the target 3-D 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             Note that positive values of D1, D2, D3 are recommended.
C             Negative values of D1, D2, D3 are allowed, but
C             the target grid will be treated as if Di were positive,
C             and the value of corresponding Oi were Oi+(Ni-1)*Di.
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                                                  
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     MTI='string'... Name of the output file with the array of
C             values, containing for each gridpoint the imaginary
C             travel times of all determined arrivals.
C             For MTI=' ' the file is not generated.
C             Description of output files M*.
C             Default:  MTI=' ' (the file is not generated)
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', MP2='string', MP3='string'... Names of the output
C             files with the array of the first, second and third
C             component of the slowness vector (i.e. with the
C             derivatives of travel time with respect to the coordinates
C             of the point).
C             Description of output files M*.
C             Defaults:  MP1=' ', MP2=' ', MP3=' ' (no file generated)
C     MTT11='string', MTT12='string', MTT22='string', MTT13='string',
C     MTT23='string', MTT33='string'... Names of the output files
C             with the array of the second derivatives of travel time
C             with respect to the coordinates of the point.
C             Description of output files M*.
C             Defaults:  MTT11=' ' to MTT33=' ' (no file 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 file 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             If you choose this option, do not forget to specify the
C             values of AMPKI1, AMPKI2 or AMPKI3.
C             Default:  AMPKI=' ' (no file generated)
C     AMPR11='string',AMPI11='string',AMPR21='string',AMPI21='string',
C     AMPR31='string',AMPI31='string',AMPR12='string',AMPI12='string',
C     AMPR22='string',AMPI22='string',AMPR32='string',AMPI32='string',
C     AMPR13='string',AMPI13='string',AMPR23='string',AMPI23='string',
C     AMPR33='string',AMPI33='string'... Names of the output files
C             with the grid values of the real or imaginary part of the
C             vectorial amplitude of the Green function.
C             Default: AMPRij=' ', AMPIij=' ' (no file generated)
C     MX4='string', MX5='string', MX6='string'... Names of the output
C             files with the X1, X2 and X3 coordinates of the initial
C             point of the ray corresponding to the arrival. (In case of
C             point source the coordinates of the source.)
C             Default:  MX4=' ', MX5=' ', MX6=' ' (no file generated)
C     MP4='string', MP5='string', MP6='string'... Names of the output
C             files with the derivatives of travel time with respect
C             to the coordinates of the source (i.e. three components of
C             the vector opposite to the slowness vector in the source).
C             Default:  MP4=' ', MP5=' ', MP6=' ' (no file generated)
C                                                   
C   Numerical parameters specifying the interpolated quantities:
C     MTTLIN=integer... Switch between bicubic and bilinear
C             interpolation of travel times:
C             MTTLIN=0 ... Bicubic interpolation
C             MTTLIN=1 ... Bilinear interpolation
C             Default: MTTLIN=0
C     TTDMAX=real... Maximum value of the second travel-time derivatives
C             MTT11 to MTT33. If the derivative is greater, it is
C             replaced by TTDMAX.  Note that TTDMAX is applied during
C             the calculation of the values of the derivatives in the
C             points along the rays, which are then used during the
C             interpolation.
C             Default:  TTDMAX=999999.
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 allowed only at the gridpoints
C             not used during interpolation within ray cells.
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 allowed only at the gridpoints
C             not used during interpolation within ray cells.
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. If the amplitude is greater, it is replaced
C             by AMPMAX.  Note that AMPMAX is applied during the
C             calculation of the values of the Green function in the
C             points along the rays, which are then used during the
C             interpolation.
C             Similarly, AMPMAX and -AMPMAX are used as the limits for
C             the individual components AMPR11 to AMPI33 of the
C             vectorial amplitude.
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             Default: AMP2D1=0., AMP2D2=0., AMP2D3=0.
C                      (amplitudes of the 3-D Green function)
C     AMPKI1=real, AMPKI2=real, AMPKI3=real... Unit vector perpendicular
C             to the surface of integration, multiplied by the surface
C             area or line length corresponding to one surface point.
C             Modification of amplitudes for Kirchhoff integral.
C             Default: AMPKI1=0., AMPKI2=0., AMPKI3=0.
C     QIPV=integer... Optional quasi-isotropic approximation of the
C             polarization vectors.
C             QIPV=0: Correct polarization vectors read from the input
C               file are used.
C             QIPV=1: Quasi-isotropic approximation of the polarization
C               vectors: The basis vectors of the ray-centred coordinate
C               system are used in place of the polarization vectors.
C               This option corresponds to versions 6.60 and older of
C               this program.
C             Default: QIPV=0
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                                                  
C Data for alternative interpolation to the individual points:
C     The interpolation to the individual points is performed if the
C     name of the output file MTTPTS is specified.
C   Name of the input file with coordinates of the individual points:
C     PTS='string'... Name of the file with the points.  If MTTPTS
C             is specified, PTS must be given too.
C             Description of file PTS
C             Default: PTS=' '
C   Name of the output file with interpolated quantities:
C     MTTPTS='string'... Name of the output file with the interpolated
C             values.  If the name is left blank, the program performs
C             the interpolation to the gridpoints of the 3-D grid.
C             If MTTPTS is specified, PTS must be specified too,
C             and the program interpolates to the individual points
C             given in file PTS.
C             Description of file MTTPTS.
C             Default: MTTPTS=' '  (interpolation to the 3-D grid)
C                                                  
C   Parameters specifying the quantities to be written into the output
C   file MTTPTS.
C     COLUMN01='string' to COLUMN69='string'... Strings which specifie
C             the quantities to be written to the columns 1 to 69 of the
C             file MTTPTS.  First three columns usually contain
C             coordinates of the point.  Column zero is reserved for the
C             name of the point.
C             Following strings are allowed:
C             COLUMNii=' ' (a space)... Nothing is to be written to the
C               column ii and to all the following columns.
C             COLUMNii='X1'... First coordinate of the point.
C             COLUMNii='X2'... Second coordinate of the point.
C             COLUMNii='X3'... Third coordinate of the point.
C             COLUMNii='NAMEOUT'... The string specifiing the quantity
C               to be written to column ii may equal to any of
C               the names of SEP parameters used for the names
C               of output formatted files with interpolated
C               quantities when interpolating to the output 3-D
C               grid.
C               Here 'NAMEOUT' stays for any of 'NUM',
C               'MTT', 'MTI', 'MHIST', 'MP1', 'MP2', 'MP3',
C               'MTT11', 'MTT12', 'MTT22', 'MTT13', 'MTT23',
C               'MTT33',
C               'MPQ11', 'MPQ21', 'MPQ31', 'MPQ41', 'MPQ12',
C               'MPQ22', 'MPQ32', 'MPQ42', 'MPQ13', 'MPQ23',
C               'MPQ33', 'MPQ43', 'MPQ14', 'MPQ24', 'MPQ34',
C               'MPQ44', 'GBW11', 'GBW1', 'GBW22', 'GBW2',
C               'AMP', 'AMPKI', 'AMPR11' to 'AMPI33',
C               'MX4', 'MX5', 'MX6', 'MP4', 'MP5', 'MP6'.
C               Description of NAMEOUT.
C             All strings may be entered either in uppercase or in
C             lowercase letters.
C             Note that for COLUMNii='NUM', the sequential number of
C             the arrival is written to column ii.
C             Also note that all quantities, including number of arrival
C             or ray history, are written to MTTPTS as real values.
C             Defaults: COLUMN01='X1', COLUMN02='X2', COLUMN03='X3',
C               COLUMN04 to COLUMN69=' '.
C Projection vector for 2-D calculations:
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 or of the
C             individual points to the plane defined by rays.
C             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 in case of
C             interpolation to the grid.
C             For the case of 2-D ray field and interpolation to the
C             individual points, projection vector must be specified.
C             Default: PROJ1=0., PROJ2=0., PROJ3=0.
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.
C             Description of sorting parameters.
C             Default: MTTSORT='MTT' (sorting according to travel 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 Optional parameters specifying the form of the quantities
C written in the output formatted files:
C     MINDIG,MAXDIG=positive integers... See the description in file
C             forms.for.
C     NUMLIN=positive integer... Number of reals or integers in one
C             line of the output files NUM or M*. See the description in
C             file 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 Input file PTS with the coordinates of points:
C (1) None to several strings terminated by / (a slash)
C (2) For each point data:
C (2.1) 'NAME',X1,X2,X3,V1,...,VN,/
C     'NAME'... Name of the point.
C     X1,X2,X3... Coordinates of the point
C     V1,...,VN... Optional values, not read by MTT.
C     /...    Values must be terminated by a slash.
C (3) / or end of file.
C
C                                                  
C Output file MTTPTS with the values interpolated in the points:
C (1) /... a slash.
C (2) For each point with at least one arrival data (2.1).  In case of
C     multivalued arrivals data (2.1) are repeated for each arrival.
C (2.1) 'NAME',R1,R2,... /
C     'NAME'... Name of the point.
C     R1,R2,... /... None to several real values terminated by a slash.
C             See parameters COLUMN01 to COLUMN69 above.
C (3) /...    A slash at the end of file.
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,LOWER,FORM2,LENGTH,
     *ERROR,WARN,RSEP1,RSEP3T,RSEP3R,RSEP3I,WARRAY,WARRAI,INDEXI,AP00
      INTEGER LENGTH
      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     LOWER,LENGTH... File length.for.
C     ERROR,WARN... File error.for.
C     RSEP1,RSEP3I,RSEP3T,RSEP3R... File sep.for.
C     WARRAI,WARRAI,FORM2... File forms.for.
C     INDEXI... File indexi.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,I4,I6,INDX,INDX1,INDX2,IT1,IIT1,JWAVES
      INTEGER I
      INTEGER LU0,LU1,LU2,LU3,LU5
      PARAMETER (LU0=10,LU1=1,LU2=2,LU3=3,LU5=5)
C     LU4=4 used in CIPTS, LU=13 used in MTTQ
      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 parameters of target grid or coordinates of target points,
C     reading filename(s) of the output file(s),
C     recording which quantities are to be written into the file(s):
      CALL CIQD
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
          IF (LPTS) THEN
C           MTT-06
            CALL ERROR('MTT-06: Projection vector not specified.')
C           For the case of 2-D ray field and interpolation to the
C           individual points, the projection vector must be specified,
C           there is no default.
          ENDIF
          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
        IRAM(I1)=IRAM((I1-NRAMP-1)*3+1)
  30  CONTINUE
      CALL INDEXI(NT,IRAM(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 LOWER(SORT)
      CALL RSEP3R('MTTORDER',ORDER,1.)
      DO 500, I4=1,NOUT
        IF (CHOUT(I4).EQ.SORT) THEN
          IF (.NOT.LPTS) THEN
            I6=N1*N2*N3
          ELSE
            I6=NPTS
          ENDIF
          DO 310, I1=0,I6-1
  300       CONTINUE
            INDX=MAXOUT-I1
            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
        ENDIF
  500 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
C             bottom of the ray cell and C denotes top of the ray cell.
C             Each integer contains the address just before
C             the parameters 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...   Displays 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
C       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
C             the parameters 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
C             the parameters 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
      REAL CIPPP
      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
      REAL DWB,DWI,HEIG,WID
      PARAMETER (DWB=0.001)
      REAL X1,X2,X3
      INTEGER INDX
      REAL ERRMAX
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 range
C             [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
C             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-----------------------------------------------------------------------
      IF (.NOT.LPTS) THEN
C       Interpolation to the grid:
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
      ENDIF
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
      IF (.NOT.LPTS) THEN
C       Interpolation to the grid:
C       Loop for all the selected gridpoints:
        DO 19, I1=K1MI,K1MA
        DO  18, I2=K2MI,K2MA
        DO  17, 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=MAXOUT-INDX+1
          CALL CIINT2(I1B,I2B,I3B,I1C,I2C,I3C,X1,X2,X3,INDX,ERRMAX,DWI)
  17    CONTINUE
  18    CONTINUE
  19    CONTINUE
      ELSE
C       Interpolation to the individual points:
        DO 25, I1=1,NPTS
          I2=MAXRAM-(3*I1)
          X1=RAM(I2+1)
          X2=RAM(I2+2)
          X3=RAM(I2+3)
          INDX=MAXOUT-I1+1
          IF (.NOT.(L2D.AND.(PROJ1.NE.0.))) THEN
            IF ((X1.LT.AMIN1(RAM(I1B+1),RAM(I2B+1),RAM(I3B+1),
     *                       RAM(I1C+1),RAM(I2C+1),RAM(I3C+1))).OR.
     *          (X1.GT.AMAX1(RAM(I1B+1),RAM(I2B+1),RAM(I3B+1),
     *                       RAM(I1C+1),RAM(I2C+1),RAM(I3C+1)))) GOTO 25
          ENDIF
          IF (.NOT.(L2D.AND.(PROJ2.NE.0.))) THEN
            IF ((X2.LT.AMIN1(RAM(I1B+2),RAM(I2B+2),RAM(I3B+2),
     *                       RAM(I1C+2),RAM(I2C+2),RAM(I3C+2))).OR.
     *          (X2.GT.AMAX1(RAM(I1B+2),RAM(I2B+2),RAM(I3B+2),
     *                       RAM(I1C+2),RAM(I2C+2),RAM(I3C+2)))) GOTO 25
          ENDIF
          IF (.NOT.(L2D.AND.(PROJ3.NE.0.))) THEN
            IF ((X3.LT.AMIN1(RAM(I1B+3),RAM(I2B+3),RAM(I3B+3),
     *                       RAM(I1C+3),RAM(I2C+3),RAM(I3C+3))).OR.
     *          (X3.GT.AMAX1(RAM(I1B+3),RAM(I2B+3),RAM(I3B+3),
     *                       RAM(I1C+3),RAM(I2C+3),RAM(I3C+3)))) GOTO 25
          ENDIF
          CALL CIINT2(I1B,I2B,I3B,I1C,I2C,I3C,X1,X2,X3,INDX,ERRMAX,DWI)
  25    CONTINUE
      ENDIF
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE
     *  CIINT2(I1B,I2B,I3B,I1C,I2C,I3C,X1,X2,X3,INDX,ERRMAX,DWI)
C
C----------------------------------------------------------------------
C Subroutine designed to interpolate to the given point.
C
      INTEGER I1B,I2B,I3B,I1C,I2C,I3C,INDX
      REAL X1,X2,X3,ERRMAX,DWI
C
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
C             the parameters of the corresponding point:
C             the first parameter of the first point: RAM(I1B+1)
C     X1,X2,X3... Coordinates of the point for interpolation.
C     INDX... Index in IRAM of the address of the point for writing
C             the interpolated values.
C     ERRMAX... Upper error bound.
C     DWI...  Range [0-DWI,1+DWI] is used for local coordinates
C             w1,w2,w3.
C No output.
C Subroutines and external functions required:
      EXTERNAL CIPPP,CILSIG
      REAL CIPPP
      LOGICAL CILSIG
C
C ...........................
C Common block /MTTC/ to store the information about triangles and rays:
      INCLUDE 'mtt.inc'
C mtt.inc.
C.......................................................................
      INTEGER I4
      REAL DWB
      PARAMETER (DWB=0.001)
      REAL Y1,Y2,Y3
      INTEGER IR
      REAL ROOT(3)
      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-----------------------------------------------------------------------
        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
      RETURN
      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
C             the parameters 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
C             coordinate of point X in the real cell, WBB is the local
C             coordinate of X 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,WW,UU,VV
      INTEGER IR,I,J
      REAL ROOT(3)
      REAL Y11,Y12,Y13,Y21,Y22,Y23
      CHARACTER*72 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,UU,VV... Coordinates of the vectors
C             at the side of the ray cell, U and V are the edges
C             at the bottom and at the top of the ray cell, W is
C             connecting points in the middles of the edges, UU and VV
C             are lengths of the vectors U and V.
C             W1,W2,W3 are later used as local coordinates.
C             Optionally vectors U and V are later used for edges
C             at the sides of the ray cell.
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
        WW=SQRT(W1*W1+W2*W2+W3*W3)
        WW=(SQRT(3.)/2.)/WW
        W1=W1*WW
        W2=W2*WW
        W3=W3*WW
        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)
        WW=SQRT(PROJ1*PROJ1+PROJ2*PROJ2+PROJ3*PROJ3)
        Y11=(PROJ2/WW)*W3-W2*(PROJ3/WW)
        Y12=(PROJ3/WW)*W1-W3*(PROJ1/WW)
        Y13=(PROJ1/WW)*W2-W1*(PROJ2/WW)
        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
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 (A.GT.0.01) THEN
        U1=RAM(I1C+1)-RAM(I1B+1)
        U2=RAM(I1C+2)-RAM(I1B+2)
        U3=RAM(I1C+3)-RAM(I1B+3)
        V1=RAM(I2C+1)-RAM(I2B+1)
        V2=RAM(I2C+2)-RAM(I2B+2)
        V3=RAM(I2C+3)-RAM(I2B+3)
        UU=U1*U1+U2*U2+U3*U3
        IF (UU.GT.DWARF) UU=SQRT(UU)
        VV=V1*V1+V2*V2+V3*V3
        IF (VV.GT.DWARF) VV=SQRT(VV)
        A=A*ABS((UU*W1+VV*W2)/AMAX1(UU,VV))
        IF (A.GT.0.01) THEN
C         MTT-30
          WRITE(TXTERR,'(A,F8.5)')
     *    'MTT-30: WB in virtual and real cells too different. A=',A
          CALL WARN(TXTERR)
C         WB was computed for the point X in the real cell. Now we are
C         looking for WBB in the virtual cell, which corresponds to the
C         above mentioned WB in the real cell. Thus WB and WBB should
C         not be too different. It is possible, that the value 0.01
C         should be enlarged.
C         This warning should not appear, or at least not too often.
        ENDIF
      ENDIF
      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
C             the parameters 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,IP)
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,IP
      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(IP+1)
      P2=RAM(IP+2)
      P3=RAM(IP+3)
      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
C             functional values of cubic equation have different signes
C             in RR(i) and 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)/(ABS(RA+RB)+DWARF)).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:
C     R1,R2,R3,R4... Real quantities.
C Output:
C     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:
C     I1,I2... Addresses of the two points in array RAM.
C Output:
C     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
      MAXOUT=MAXRAM
      MRAMP=MAXOUT-N1*N2*N3
      IF (MRAMP.LE.1) CALL CIEROR(1,TXTERR)
      NRAMT=MRAMP+1
      DO 5, I1=NRAMT,MAXOUT
        IRAM(I1)=0
   5  CONTINUE
      RETURN
C
  10  CONTINUE
C     MTT-10
      CALL ERROR('MTT-10: This specification of the interpolation grid
     *may cause problems. Please, specify D1,D2,D3 and N1,N2,N3 greater
     * than 0. Di may equal 0 in case that corresponding Ni equals 1.')
C
      END
C
C
C=======================================================================
C
      SUBROUTINE CIPTS
C
C----------------------------------------------------------------------
C Subroutine to read in the file with coordinates of individual points
C for interpolation. The read quantities are then written into common
C 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 LU4,I1
      PARAMETER (LU4=4)
      CHARACTER*52 TXTERR
C-----------------------------------------------------------------------
C
C     Reading the names and the coordinates of points:
      OPEN(LU4,FILE=FILPTS,STATUS='OLD')
      READ(LU4,*) (TXTPTS(I1),I1=1,20)
      NPTS=0
      MAXOUT=MAXRAM
   20 CONTINUE
        IF (MAXOUT-2.LT.1) CALL CIEROR(1,TXTERR)
        IF (NPTS+1.GT.MPTS) THEN
C         MTT-07
          CALL ERROR('MTT-07: Too small array TXTPTS')
        ENDIF
        TXTPTS(NPTS+1)='$'
        RAM(MAXOUT-2)=0.
        RAM(MAXOUT-1)=0.
        RAM(MAXOUT  )=0.
        READ(LU4,*,END=29) TXTPTS(NPTS+1),(RAM(I1),I1=MAXOUT-2,MAXOUT)
        IF (TXTPTS(NPTS+1).EQ.'$') THEN
          GOTO 29
        ENDIF
        MAXOUT=MAXOUT-3
        NPTS=NPTS+1
      GOTO 20
   29 CONTINUE
      CLOSE(LU4)
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.)
      MRAMP=MAXOUT-NPTS
      IF (MRAMP.LE.1) CALL CIEROR(1,TXTERR)
      NRAMT=MRAMP+1
      DO 30, I1=NRAMT,MAXOUT
        IRAM(I1)=0
  30  CONTINUE
      RETURN
      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
C             the parameters 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 Subroutines and external functions required:
      EXTERNAL LENGTH
      INTEGER LENGTH
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
      REAL OUTMIN(MOUT),OUTMAX(MOUT)
      CHARACTER*(4+8*MOUT) FORMAT
      CHARACTER*20 FORMA1,TEXT
      CHARACTER*52 TXTERR
C.......................................................................
C
C     Reading filename of the output file with points:
      CALL RSEP3T('MTTPTS',FILOUT,' ')
      IF (FILOUT.NE.' ') THEN
        LPTS=.TRUE.
      ELSE
        LPTS=.FALSE.
      ENDIF
      CALL RSEP3I('MTTLIN',MTTLIN,0)
C
      IF (.NOT.LPTS) THEN
C       Interpolation to the 3-D grid:
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
        CALL CIGRID
      ELSE
C       Interpolation to the individual points:
        NOUT=0
C       Reading the list QNAMES of the possible interpolated quantities:
        CALL MTTQ
C       Reading the input filename:
        CALL RSEP3T('PTS',FILPTS,' ')
        IF (FILPTS.EQ.' ') THEN
C         MTT-19
          CALL ERROR('MTT-19: File PTS not specified.')
C         If MTTPTS is given, PTS must be specified too.
C         See description of input data.
        ENDIF
C       Reading the data specifying the quantities to be written into
C       the individual columns of the output file:
        FORMA1='COLUMN00'
        I1=0
   7    CONTINUE
          I1=I1+1
          FORMA1(8:8)=CHAR(ICHAR('0')+MOD(I1,10))
          FORMA1(7:7)=CHAR(ICHAR('0')+I1/10)
          IF     (I1.EQ.1) THEN
            CALL RSEP3T(FORMA1,TEXT,'x1')
          ELSEIF (I1.EQ.2) THEN
            CALL RSEP3T(FORMA1,TEXT,'x2')
          ELSEIF (I1.EQ.3) THEN
            CALL RSEP3T(FORMA1,TEXT,'x3')
          ELSE
            CALL RSEP3T(FORMA1,TEXT,' ')
          ENDIF
          IF (TEXT.NE.' ') THEN
            CALL LOWER(TEXT)
            IF (TEXT.EQ.'x1'   ) GOTO 9
            IF (TEXT.EQ.'x2'   ) GOTO 9
            IF (TEXT.EQ.'x3'   ) GOTO 9
            IF (TEXT.EQ.'num'  ) GOTO 9
            IF (TEXT.EQ.'mtt'  ) GOTO 9
            IF (TEXT.EQ.'mhist') GOTO 9
            IF (TEXT.EQ.'mp1'  ) GOTO 9
            IF (TEXT.EQ.'mp2'  ) GOTO 9
            IF (TEXT.EQ.'mp3'  ) GOTO 9
            DO 8, I2=1,MOUT
              IF (TEXT.EQ.QNAMES(I2)) GOTO 9
  8         CONTINUE
C           MTT-20
            CALL ERROR('MTT-20: Wrong value of COLUMN.')
C           See allowed values of COLUMNii in the
C           description of file SEP.
  9         CONTINUE
            IF (I1.EQ.69) THEN
              CALL RSEP3T('COLUMN70',TEXT,' ')
              IF (TEXT.NE.' ') THEN
C               MTT-21
                CALL ERROR('MTT-21: More than 69 COLUMNs.')
C               Currently up to 69 values of COLUMNii may be specified.
C               See allowed values of COLUMNii in the
C               description of file SEP.
              ENDIF
            ENDIF
C
            NOUT=NOUT+1
            IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR)
            CHOUT(NOUT)=TEXT
            FOUT(NOUT)=' '
C           Reading the input data for the interpolated quantity:
            CALL MTTQD(TEXT)
            IF (TEXT.EQ.'mtt') THEN
              IF (NOUT+3.GT.MOUT) CALL CIEROR(27,TXTERR)
              NOUT=NOUT+1
              CHOUT(NOUT)=' '
              FOUT(NOUT)=' '
              NOUT=NOUT+1
              CHOUT(NOUT)=' '
              FOUT(NOUT)=' '
              NOUT=NOUT+1
              CHOUT(NOUT)=' '
              FOUT(NOUT)=' '
            ENDIF
            GOTO 7
          ENDIF
C       End of the loop over future columns of the output file.
        NQ=5+NOUT
        CALL CIPTS
      ENDIF
      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.' ') .OR.(CHOUT(I1).EQ.'x1').OR.
     *          (CHOUT(I1).EQ.'x2').OR.(CHOUT(I1).EQ.'x3')) 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.' ') .OR.(CHOUT(I1).EQ.'x1').OR.
     *          (CHOUT(I1).EQ.'x2').OR.(CHOUT(I1).EQ.'x3')) 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
          IF (MTTLIN.EQ.1) THEN
C           ............................................................
C           Travel time - 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) THEN
              RAM(NRAMT+I5)=RAM(NRAMT+I5)+
     *           WB*W3*RAM(I3B+(5+I5))+
     *           WC*W3*RAM(I3C+(5+I5))
            ENDIF
C           ............................................................
          ELSE
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+(5+I5))+AWC*RAM(I1C+(5+I5))+
     *          CIBB(WB,I1B,I1C,I1B+5+I5)+CIBB(WC,I1C,I1B,I1C+5+I5)
            TA2=AWB*RAM(I2B+(5+I5))+AWC*RAM(I2C+(5+I5))+
     *          CIBB(WB,I2B,I2C,I2B+5+I5)+CIBB(WC,I2C,I2B,I2C+5+I5)
            PA(1,1)=AWB*RAM(I1B+(5+I5+1))+AWC*RAM(I1C+(5+I5+1))
            PA(2,1)=AWB*RAM(I1B+(5+I5+2))+AWC*RAM(I1C+(5+I5+2))
            PA(3,1)=AWB*RAM(I1B+(5+I5+3))+AWC*RAM(I1C+(5+I5+3))
            PA(1,2)=AWB*RAM(I2B+(5+I5+1))+AWC*RAM(I2C+(5+I5+1))
            PA(2,2)=AWB*RAM(I2B+(5+I5+2))+AWC*RAM(I2C+(5+I5+2))
            PA(3,2)=AWB*RAM(I2B+(5+I5+3))+AWC*RAM(I2C+(5+I5+3))
            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+(5+I5))+AWC*RAM(I3C+(5+I5))+
     *          CIBB(WB,I3B,I3C,I3B+5+I5)+CIBB(WC,I3C,I3B,I3C+5+I5)
              PA(1,3)=AWB*RAM(I3B+(5+I5+1))+AWC*RAM(I3C+(5+I5+1))
              PA(2,3)=AWB*RAM(I3B+(5+I5+2))+AWC*RAM(I3C+(5+I5+2))
              PA(3,3)=AWB*RAM(I3B+(5+I5+3))+AWC*RAM(I3C+(5+I5+3))
              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
          ENDIF
C         ..............................................................
        ELSEIF (CHOUT(I5).EQ.'mhist') THEN
C         Ray history - no interpolation:
          IRAM(NRAMT+I5)=IRAM(I1B+(5+I5))
        ELSEIF ((CHOUT(I5).EQ.' ') .OR.(CHOUT(I5).EQ.'x1').OR.
     *          (CHOUT(I5).EQ.'x2').OR.(CHOUT(I5).EQ.'x3')) 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
  29    CONTINUE
  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:
      IF (.NOT.LPTS) THEN
        I6=N1*N2*N3
      ELSE
        I6=NPTS
      ENDIF
      IF (I6.GT.MRAMP) THEN
        I1=MRAMP-I6
        WRITE(TXTERR,'(A,I9,A)')
     *  'MTT-31: Array RAM too small,',I1,' units missing.'
        CALL CIEROR(31,TXTERR)
      ENDIF
      NRAMP=0
      I5=0
      DO 110, I1=0,I6-1
        NRAMP=NRAMP+1
        I4=0
        INDX=MAXOUT-I1
  105   CONTINUE
C         Check for the consistentcy of the results
          IF (INDX.LE.MRAMP.OR.INDX.GT.MAXOUT) 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
            IF (.NOT.LPTS) IRAM(NRAMP)=I4
            I5=I5+I4
            GOTO 110
          ENDIF
          I4=I4+1
        GOTO 105
  110 CONTINUE
      IF (NRAMP.NE.I6) 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 ((.NOT.LPTS).AND.(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).AND.(.NOT.LPTS)) THEN
C       Interpolation to the grid:
        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
        DO 140, I4=1,NOUT
          IF (FOUT(I4).EQ.' ') GOTO 139
          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=MAXOUT-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.MAXOUT) 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
  139     CONTINUE
  140   CONTINUE
      ELSEIF ((NOUT.GE.1).AND.(LPTS).AND.(FILOUT.NE.' ')) THEN
C       Interpolation to individual points:
        IF (NOUT.GT.MRAMP) THEN
          I1=MRAMP-NOUT
          WRITE(TXTERR,'(A,I9,A)')
     *    'MTT-31: Array RAM too small,',I1,' units missing.'
          CALL CIEROR(31,TXTERR)
        ENDIF
        OPEN(LU5,FILE=FILOUT)
        WRITE(LU5,'(A)') '/'
C       Setting output format:
        FORMAT(1:4)='(4A,'
        DO 149, I1=1,NOUT
          OUTMIN(I1)=GIANT
          OUTMAX(I1)=DWARF
  149   CONTINUE
C       Loop over points:
        DO 170, I1=1,NPTS
C         I1 counts the points
          INDX=MAXOUT-I1+1
          I2=0
          I3=0
C         Loop over multivalued arrivals:
  150     CONTINUE
C           I2 counts the arrivals
            INDX=IRAM(INDX)
            IF (INDX.NE.0) THEN
              I3=0
C             I3 counts the output values to be written to MTTPTS
              I2=I2+1
C             Looking for the minimum and maximum values:
              DO 160, I4=1,NOUT
                IF (CHOUT(I4).EQ.'num') THEN
                  I3=I3+1
                  OUTMIN(I3)=AMIN1(OUTMIN(I3),FLOAT(I2))
                  OUTMAX(I3)=AMAX1(OUTMAX(I3),FLOAT(I2))
                ELSEIF (CHOUT(I4).EQ.'mhist') THEN
                  I3=I3+1
                  OUTMIN(I3)=AMIN1(OUTMIN(I3),FLOAT(IRAM(INDX+I4)))
                  OUTMAX(I3)=AMAX1(OUTMAX(I3),FLOAT(IRAM(INDX+I4)))
                ELSEIF (CHOUT(I4).EQ.'x1') THEN
                  I3=I3+1
                  OUTMIN(I3)=AMIN1(OUTMIN(I3),RAM(MAXRAM-3*I1+1))
                  OUTMAX(I3)=AMAX1(OUTMAX(I3),RAM(MAXRAM-3*I1+1))
                ELSEIF (CHOUT(I4).EQ.'x2') THEN
                  I3=I3+1
                  OUTMIN(I3)=AMIN1(OUTMIN(I3),RAM(MAXRAM-3*I1+2))
                  OUTMAX(I3)=AMAX1(OUTMAX(I3),RAM(MAXRAM-3*I1+2))
                ELSEIF (CHOUT(I4).EQ.'x3') THEN
                  I3=I3+1
                  OUTMIN(I3)=AMIN1(OUTMIN(I3),RAM(MAXRAM-3*I1+3))
                  OUTMAX(I3)=AMAX1(OUTMAX(I3),RAM(MAXRAM-3*I1+3))
                ELSEIF (CHOUT(I4).NE.' ') THEN
                  I3=I3+1
                  OUTMIN(I3)=AMIN1(OUTMIN(I3),RAM(INDX+I4))
                  OUTMAX(I3)=AMAX1(OUTMAX(I3),RAM(INDX+I4))
                ENDIF
  160         CONTINUE
              GOTO 150
            ENDIF
C         End of the loop over multivalued arrivals.
          IF (I3.NE.0) THEN
            CALL FORM2(I3,OUTMIN,OUTMAX,FORMAT(5:4+8*I3))
          ENDIF
  170   CONTINUE
C
C       Writing the values at the points:
C       Loop over points:
        DO 200, I1=1,NPTS
C         I1 counts the points
          INDX=MAXOUT-I1+1
          I2=0
C         Loop over multivalued arrivals:
  180     CONTINUE
C           I2 counts the arrivals
            INDX=IRAM(INDX)
            IF (INDX.NE.0) THEN
              I2=I2+1
C             Preparing the output values for the point and the arrival:
              I3=0
C             I3 counts the output values to be written to MTTPTS
              DO 190, I4=1,NOUT
                IF (CHOUT(I4).EQ.'num') THEN
                  I3=I3+1
                  RAM(I3)=FLOAT(I2)
                ELSEIF (CHOUT(I4).EQ.'mhist') THEN
                  I3=I3+1
                  RAM(I3)=FLOAT(IRAM(INDX+I4))
                ELSEIF (CHOUT(I4).EQ.'x1') THEN
                  I3=I3+1
                  RAM(I3)=RAM(MAXRAM-3*I1+1)
                ELSEIF (CHOUT(I4).EQ.'x2') THEN
                  I3=I3+1
                  RAM(I3)=RAM(MAXRAM-3*I1+2)
                ELSEIF (CHOUT(I4).EQ.'x3') THEN
                  I3=I3+1
                  RAM(I3)=RAM(MAXRAM-3*I1+3)
                ELSEIF (CHOUT(I4).NE.' ') THEN
                  I3=I3+1
                  RAM(I3)=RAM(INDX+I4)
                ENDIF
  190         CONTINUE
              WRITE(LU5,FORMAT(1:4+8*I3))
     *          '''',TXTPTS(I1)(1:LENGTH(TXTPTS(I1))),'''',
     *          (' ',RAM(I4),I4=1,I3),' /'
              GOTO 180
            ENDIF
C         End of the loop over multivalued arrivals.
  200   CONTINUE
        WRITE(LU5,'(A)') '/'
        CLOSE(LU5)
      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.
C       Number of 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 'indexi.for'
C     indexi.for.
C
C=======================================================================
C