C 
C Program NET for a network shortest path calculation of the first
C arrival travel time and the corresponding rays, on a rectangular grid
C of points.
C
C Version: 3.10
C Date: 1998, September 26
C
C Authors:
C     Ludek Klimes
C             Department of Geophysics, Charles University Prague
C             Ke Karlovu 3
C             121 16  Praha 2,  Czech Republic
C             E-mail: klimes@seis.karlov.mff.cuni.cz
C     Michal Kvasnicka
C             Department of Geophysics, Charles University Prague
C             Ke Karlovu 3
C             121 16  Praha 2,  Czech Republic
C             E-mail: qasnicka@seis.karlov.mff.cuni.cz
C
C.......................................................................
C
C Program 'NET' is able to:
C (A) Perform both 2-D and 3-D ray tracing in a rectangular grid
C     discretization of an arbitrarily complex model according to
C     the paper by Klimes and Kvasnicka (1994).
C (B) Compute first arrivals of both direct and multiply reflected
C     waves.
C (C) Compute first arrivals of both P-waves, S-waves and converted
C     waves.
C (D) Employ forward stars of any kind (template forward stars are
C     stored in a disk file, and may be modified and optimized in
C     various ways).
C (E) Trace rays both in the whole rectangular model volume and in a
C     Fresnel volume corresponding to a two-point ray, or in a volume
C     of any other shape.  This may increase accuracy of two-point ray
C     tracing by one order, or to decrease computational complexity
C     by about two orders.
C (F) Adjust the size of a forward star at a structural interface in
C     order to maximize the accuracy.
C (G) Adjust the size of a forward star according to a local degree of
C     heterogeneity in order to maximize the accuracy.
C (H) Estimate maximum errors of all computed arrival times.
C (I) In a special 2-D mode referred hereinafter as TTT perform 2-D
C     grid travel-time tracing of the second order according to the
C     paper by Klimes (1996).
C
C References:
C     Klimes L. and Kvasnicka M. (1994):  3-D network ray tracing.
C             Geophys.J.int., 116, 726-738.
C     Klimes L. (1996):  Grid travel-time tracing: second-order method
C             for the first arrivals in smooth media.
C             PAGEOPH, 148, 539-563.
C
C.......................................................................
C
C Notes concerning the method:
C
C Although the program deals with arbitrarily complex models, the
C accuracy of the shortest path approximation of the travel time and
C rays, of course, depends on the model discretization step and on the
C model complexity.  Fortunately, the relative travel-time error can
C be estimated.
C
C Because of the local optimization of the sizes of forward stars, the
C network becomes oriented (i.e. backward stars do not equal forward
C stars) and the resulting shortest path travel time need not be
C reciprocal.
C
C This program employs two-point (trapezoidal) travel time approximation
C along a network edge.
C
C Hereinafter the way of building the rectangular grid and the network
C is briefly described.
C
C Big bricks, small bricks, gridpoints:
C     The rectangular volume bounded by coordinate limits X1MIN,X1MAX,
C     X2MIN,X2MAX, and X3MIN,X3MAX is divided into N1*N2*N3 big bricks.
C     If the numbers L1,L2,L3 are specified in addition to N1,N2,N3,
C             each big brick is subdivided into L1*L2*L3 small bricks.
C     If the numbers L1,L2,L3 are not specified, each big brick contains
C             just one small brick, as large as big one.
C     Gridpoints are defined as centres of the small bricks.
C
C Gridpoint indexing (positions of network nodes):
C     Gridpoints are indexed by integers 1,2,...,N1*L1*N2*L2*N3*L3.
C     Indices of gridpoints describe the their locations within the
C     model volume.  The gridpoint indexing is similar to 1-D indexing
C     of 3-D Fortran array.  A gridpoint index define the location of
C     the network node situated at the gridpoint.
C
C Network node indexing (addresses of network nodes):
C     If L1,L2,L3 are not specified, gridpoint indices equal to the
C     corresponding gridpoint indices.
C     Description for experienced users:
C     Network nodes are indexed just within the computational volume
C     composed of big bricks.  Nodes are indexed consecutively within
C     each big brick.  The node indexing within a big brick is similar
C     to 1-D indexing of 3-D Fortran array.  The number of indices
C     within each big brick is L1*L2*L3.  If the number of big bricks
C     within the computational volume is L4, then the node indices take
C     the values 1,2,...,L1*L2*L3*L4.
C
C Index file:
C     Index file enables to specify the computational volume of a
C     complex shape.  The network ray tracing is then performed only
C     in the specified computational volume instead of the whole
C     rectangular model volume.  This enables, e.g., to perform the
C     network ray tracing just in the estimated Fresnel volume.
C     The option of the index file is designed for experienced users.
C     The index file to specifies mapping of gridpoints onto the network
C     nodes.  If it exists, it contains N1*N2*N3 items (integer
C     indices), each corresponding to one big brick.  Each item contains
C     the index of the first node of the big brick divided by L1*L2*L3.
C
C Network nodes:
C     The set of network nodes if composed of:
C     (1) Given source nodes.
C     (2) The subset of gridpoints:
C             A gridpoint belongs to the network if contained in the big
C             brick belonging to the network.  The relation of big
C             bricks with respect to the network is defined by means of
C             the index file:
C             If the index file is not specified, all big bricks belong
C               to the network.
C             If the index file is specified, each its item (index)
C               corresponds to one big brick:
C               If the index is zero, the big brick does not belong to
C                 the network.
C               If the index is not zero, the big brick belongs to the
C                 network and the index denotes the set of nodes within
C                 the brick.
C             Big bricks belonging to the network are assumed to be
C             indexed by positive integers.
C     (3) Given receiver nodes.
C
C Generation of network edges:
C     Template forward star:
C             The concept of a template forward star is based on the
C             idea of a set of nodes (gridpoints) periodic in the space:
C             Template forward star is of the same geometry for all
C             network nodes.  The template forward star is assumed
C             symmetric with respect to axis reflections and
C             interchangings.  The edges corresponding to the template
C             forward star are read from the file 'net.fs*'.
C             Template forward stars are considered to contain the edge
C             edge connecting the central node with itself, although
C             this edge is not explicitly referred in the file
C             'net.fs*'.
C             Size of a template forward star roughly describes the
C               extent of the forward star in grid intervals.
C             Many different kinds of template forward stars may be
C             considered:
C             Full 3-D cubic forward star, size NFS:
C               Contains all edges connecting the central node with all
C               nodes within the cube of the linear size of 2*NFS+1
C               nodes, centred at the central node.
C             Full 2-D square forward star, size NFS:
C               Contains all edges connecting the central node with all
C               nodes within the square of the linear size of 2*NFS+1
C               nodes, centred at the central node.
C             Full 3-D spherical forward star, size NFS:
C               Contains all edges connecting the central node with all
C               nodes within the sphere of the radius of SQRT(NFS*NFS+2)
C               grid intervals, centred at the central node.
C             Full 3-D circular forward star, size NFS:
C               Contains all edges connecting the central node with all
C               nodes within the circle of the radius of SQRT(NFS*NFS+1)
C               grid intervals, centred at the central node.
C             Moser-Saito 3-D cubic forward star, size NFS:
C               Full 3-D cubic forward star, size NFS, without the edges
C               that would be parallel in a homogeneous medium.
C             Moser-Saito 2-D square forward star, size NFS:
C               Full 2-D square forward star, size NFS, without the
C               edges that would be parallel in a homogeneous medium.
C             Optimized 3-D spherical forward star, size NFS:
C               Full 3-D spherical forward star, size NFS, without the
C               edges that can be removed while keeping,
C               in a homogeneous medium, each angular cone of the
C               radius grater than 1/(SQRT(2)*NFS) radians covered by
C               edges.
C             Optimized 2-D circular forward star, size NFS:
C               Full 2-D circular forward star, size NFS, without the
C               edges that can be removed while keeping,
C               in a homogeneous medium, each angle of the size greater
C               than 1/NFS radians covered by edges.
C               Note that this forward star is not a subset of the
C               optimized 3-D spherical forward star.
C
C     Actual forward stars (without receiver nodes):
C             All forward stars corresponding to nodes situated within
C             a same small brick are the same, equal to the forward
C             star corresponding to its central gridpoint.
C             A forward star corresponding to a gridpoint is the
C             intersection of the template forward star with the set of
C             network nodes.
C             At source nodes, full spherical or circular template
C               forward stars are considered.
C             At gridpoint nodes, optimized spherical or circular
C               template forward stars from the files 'net.fs3' and
C               'net.fs2' are considered.
C     Receiver backward stars:
C             The edges generated by the backward stars of the receiver
C             nodes are added to the edges generated by the above
C             forward stars.  The backward stars are the same as forward
C             stars from the central nodes of the corresponding small
C             bricks.
C             At receivers, full spherical or circular template
C               backward stars are considered.
C
C.......................................................................
C
C Files 'net.fs2' and 'net.fs3' with template forward stars:
C     The files describing forward stars of various sizes must have the
C     names 'net.fs2' and 'net.fs3'.  Forward stars submitted with this
C     version:
C     'net.fs2'... Optimized spherical 2-D forward stars, sizes 1-100.
C     'net.fs3'... Optimized circular 3-D forward stars, sizes 1-27.
C
C Note:
C     Files 'net.fs2' and 'net.fs3' may be reduced to approximately 2/3
C     of their size if replacing all multiple spaces by a single space.
C     This may slightly save disk space and speed up program starting.
C
C TTT mode:
C     Template forward stars 'net.fs2' and 'net.fs3' are not required.
C
C.......................................................................
C
C                                                    
C Data files:
C
C Input data read from the * external unit:
C     The data consist of a single character string, read by list
C     directed (free format) input.  Thus the string has to be enclosed
C     in apostrophes.  The interactive * external unit may be redirected
C     to the file containing the string.
C (1) 'SEP',/
C     'SEP'...String in apostrophes containing the name of the input SEP
C             parameter file containing the data specifying grid
C             dimensions, references to other data files and optionally
C             some numerical parameters.
C             Description of file SEP
C     Default: 'SEP'='net.h'
C
C                                                     
C Data file SEP has the form of the SEP (Stanford Exploration Project)
C parameter file:
C     All the data are specified in the form of PARAMETER=VALUE, e.g.
C     N1=50, with PARAMETER directly preceding = without intervening
C     spaces and with VALUE directly following = without intervening
C     spaces.  The PARAMETER=VALUE couple must be delimited by a space
C     or comma from both sides.
C     The PARAMETER string is not case-sensitive.
C     PARAMETER= followed by a space resets the default parameter value.
C     All other text in the input files is ignored.  The file thus may
C     contain unused data or comments without leading comment character.
C     Everything between comment character # and the end of the
C     respective line is ignored, too.
C     The PARAMETER=VALUE couples may be specified in any order.
C     The last appearance takes precedence.
C Other data filenames:
C     NET=string... String with the name of the input data file
C             containing the references to other input and output files.
C             Description of input data NET
C             Default: NET='net.dat'
C     FSTAB=string... String with the name of the optional output file
C             containing the information on memory required to store
C             forward stars.  Useful if array dimension MFS declared in
C             net.inc has to be updated.
C             The file is not created by default.
C             Default: FSTAB=' '
C Data specifying grid dimensions:
C             Boundaries of the model volume:
C               X1MIN=O1-0.5*D1,  X1MAX=X1MIN+FLOAT(N1)*D1,
C               X2MIN=O2-0.5*D2,  X2MAX=X2MIN+FLOAT(N2)*D2,
C               X3MIN=O3-0.5*D3,  X3MAX=X3MIN+FLOAT(N3)*D3.
C             Because the model volume is divided into the small bricks
C             and the gridpoints are situated in the centres of the
C             small bricks, the left-hand model boundary is situated
C             half a grid interval to the left from the leftmost
C             gridpoints.  Similarly for other model boundaries.
C             In a 2-D model, boundaries parallel with the model plane
C             are ignored.
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             N3 need not be specified for a 2-D model.
C             Default: N3=1
C     D1=positive real... Grid interval in the direction of the first
C             coordinate axis.
C             Default: D1=1.
C     D2=positive real... Grid interval in the direction of the second
C             coordinate axis.
C             Default: D2=1.
C     D3=positive real... Grid interval in the direction of the third
C             coordinate axis.
C             Default: D3=1.
C     O1=real... First coordinate of the grid origin (first point of the
C             grid).
C             Default: O1=0.
C     O2=real... Second coordinate of the grid origin.
C             Default: O2=0.
C     O3=real... Third coordinate of the grid origin.
C             Default: O3=0.
C Additional parameters for network ray tracing in Fresnel volumes:
C             L1,L2,L3 may be left out and should only be specified by
C             experienced users.
C             (If the numbers L1,L2,L3 are not specified, each big brick
C             contains just one small brick, as large as big one.)
C     L1=positive integer... Number of small bricks in one big brick in
C             the direction of axis X1.  If specified, must be positive.
C             Default: L1=1
C     L2=positive integer... Number of small bricks in one big brick in
C             the direction of axis X2.  If specified, must be positive.
C             Default: L2=1
C     L3=positive integer... Number of small bricks in one big brick in
C             the direction of axis X3.  If specified, must be positive.
C             Default: L3=1
C Note: Only one of N1*L1, N2*L2, N3*L3 may equal 1.
C             In such a case, two-dimensional ray tracing is performed:
C             Two-dimensional forward stars are considered,
C             model boundaries in the proper direction are ignored,
C             source and receiver coordinates perpendicular to the plane
C             of computation have to equal.
C Numerical parameters:
C     NFSMAX=integer... Options: NFSMAX=positive integer, 0, -1.
C             Default: NFSMAX=0.
C     NFSMAX=positive integer:
C             NFSMAX is the maximum size of a forward star.
C             Simultaneously, the default size of a forward star if the
C             size is not adjusted.
C             The size of a forward star is measured in grid intervals,
C             and describes the maximum length of the edge between two
C             consecutive points (nodes) of a ray.
C             The maximum size of a forward star has to be estimated by
C             the user according the model, in order to prevent the ray
C             to skip over some details in the model.  These details may
C             be, e.g., small low or high velocity regions or bumps on
C             the structural interfaces.
C             NFSMAX has to be positive and not exceeding dimension
C             MFSMAX in the common block /FS/.
C     NFSMAX=0: The maximum size of forward stars is determined
C             automatically: New NFSMAX**(-2) is average of NFS**(-2),
C             where NFS is optimum size at each node.  This is default.
C TTT mode (specified by NFSMAX=-1):
C     N1,N2,N3... N3=1 is obligatory in the input data.
C     L1,L2,L3... L1=L2=L3=1 is strongly recommended.
C     NFSMAX=-1: 2-D 'Second-order' grid travel-time tracing by Klimes
C             (1996) is used instead of network ray tracing.
C             Default: NFSMAX=0.
C     RIDGE1=positive real, RIDGE2=positive real...
C             Numerical parameters controlling the check for
C             the ridges of the first-arrival travel time.  At the
C             ridges, the rays going from different directions meet.
C             The interpolation across the ridges should be avoided
C             because it would introduce the first-order errors.
C             Instead, paraxial approximations from gridpoints at each
C             side of the ridge should be applied.  For the paraxial
C             approximations, the second travel-time derivative along
C             the gridline crossing the ridge is taken zero.
C             The serious problem is to indicate the ridge.
C             The slowness-vector difference across the ridge is
C             negative and lower than would correspond to the
C             slowness-vector difference along the neighbouring gridline
C             segment, i.e. to the second travel-time derivative at the
C             corresponding side of the ridge.
C             The following algorithm is now applied:
C               The slowness-vector difference along the neighbouring
C               gridline is multiplied by RIDGE1 if negative or by 0 if
C               positive, and decreased by RIDGE2 times the norm of the
C               slowness gradient.  If the result is greater than the
C               slowness-vector difference along the gridline being
C               checked, the ridge is indicated.
C             RIDGE1 should thus aways be greater than 1, whereas
C             RIDGE2 should be positive, small with respect to 1.
C             Values like RIDGE1=1.1, 1.2, 1.3, 1.4, 1.5, 3.0, 6.0,
C             and RIDGE2=0.05, 0.10, 0.20, 1.00 have been tested with
C             very unstable results.  The questionable defaults of
C             RIDGE1=1.5 and RIDGE2=0.10 have finally been chosen.
C             Very large values of RIDGE1 and RIDGE2 should disable
C             the test for ridges and thus correspond to program TTT-2D
C             version 0.17 (slightly fixed version 0.07).
C             Default: RIDGE1=1.5, RIDGE2=0.1.
C     VER1=real... If the distance of the centre of wavefront curvature
C             from the gridline segment is greater than VER1 grid
C             intervals:
C               Cubic interpolation of travel time along the gridline
C               segment,
C             else:  cubic interpolation of travel time residual along
C               the gridline segment.  The travel time residual is taken
C               with respect to spherical wavefronts having the common
C               centre at the point of intersection of two given
C               slowness vectors.
C             For more details refer to Klimes (1996), Sec.3.1.4, where
C             VER1 is denoted by N.
C             The results are not very sensitive to this parameter and
C             there has never arisen a need to change the default value.
C             Default: VER1=9.999.
C     VER2=real... Smoothing and stabilizing slowness vectors:
C             Residual travel times with respect to quadratic or
C             spherical interpolation are approximated by
C               TTRES=(1-VER2)*TTCUB+VER2*TTLIN,
C             where TTCUB is cubic interpolation, and TTLIN is linear
C             approximation adjusting the resulting slowness vector
C             according to travel time residuals.
C             VER1=0 corresponds to interpolation described in Sections
C               3.1.3.1 and 3.1.4.1 and is default in TTT-2D ver. 0.06.
C             VER1=1 corresponds to interpolation described in Sections
C               3.1.3.2 and 3.1.4.2 and is default since TTT-2D version
C               0.07.
C             Default: VER2=1.0.
C Example of data set SEP
C
C                                                     
C Main input file 'NET':
C     Sequential file, read by list directed (free format) input,
C     containing model parameters, source/receiver coordinates, and
C     names of other input and output files.
C     In the list of input data below, each numbered paragraph indicates
C     the beginning of a new input operation (new read statement).
C     'ITEMS' in the list of input variables enclosed in apostrophes
C     represent character strings enclosed in apostrophes.  Otherwise,
C     if the first letter of the symbolic name in the list of input
C     variables is I-N, the corresponding value in input data is
C     integer, otherwise, the input parameter is of the type real.
C     / in the list of input variables indicates an obligatory slash.
C     The slash may also be used instead of default values.
C (1) 'SRC','REC','RAYS','END',/
C     'SRC'...Name of the input file with source coordinates.
C             Description of input data SRC
C     'REC'...Name of the input file with receiver coordinates.
C             If blank, no rays are stored within the file 'RAYS'.
C             Description of input data REC
C     'RAYS'..Name of the output file with rays.
C             If blank, no rays are stored within the file 'RAYS'.
C             Description of output data RAYS
C     'END'...Name of the output file with endpoints of rays (receiver
C             coordinates, receiver arrival times, and estimates of the
C             corresponding maximum travel-time errors.
C             If blank, no file 'END' is generated.
C             Description of output data END
C     Default: 'REC'=' ', 'RAYS'=' ', 'END'=' '.
C TTT mode:
C     'REC'...File read if specified but its data ignored.  No receivers
C             are considered in the present version.
C     'RAYS'..Ignored - output file not generated.
C     'END'...Ignored - output file not generated.
C (2) NREFL,/
C     NREFL...Number of reflections.
C             Attention concerning this version:
C               If, in the case of reflections (NREFL.GT.0), some
C               'VEL(I)' differs from 'VEL(0)' at lines (3.2) below, the
C               errors cannot be evaluated correctly and the filename
C               'END' above should be left blank.
C             Note: To calculate a reflected wave, it is recommended to
C               keep NREFL=0 and to submit the reflector points as a set
C               'REC' of the receiver points.  Then to submit the output
C               'END' file as the source file 'SRC' for the subsequent
C               calculation of the reflected wave.  In this way, the
C               arrival time errors connected with the discretization of
C               the reflector can be removed.  On the other hand, the
C               whole rays are split into several files 'RAYS' and have
C               to be put together after finishing the network ray
C               tracing.
C     Default: NREFL=0.
C (3) Once (3.1), then NREFL-times (3.2) and (3.1):
C TTT mode:   NREFL=0.
C (3.1) 'IND(I)','VEL(I)','ICB(I)','TT(I)','ERR(I)','PRED(I)','NFS(I)',/
C     'IND(I)'... Name of the input index file, specifying for each
C             big brick if its gridpoints belong to the network.
C             May be blank - then the default indexing is assumed.
C             Must not be blank if (L1.GT.1.OR.L2.GT.1.OR.L3.GT.1) at
C             input (1).
C             This file should only be used by experienced users, others
C             should leave it blank.
C             Description of input data IND(I)
C     'VEL(I)'... Name of the input file containing velocities at all
C             network nodes, for I-times reflected wave.
C             Has to be specified.
C             Description of input data VEL(I)
C     'ICB(I)'... Name of the input file containing indices of
C             (geological) blocks.  Only network edges corresponding to
C             a forward star of the size NFS=1 are allowed to cross an
C             interface between two different blocks (i.e. blocks with
C             different indices).  This considerably increases the
C             accuracy in presence of structural interfaces (often 5
C             times).  Note that the limitation to size 1 is considered
C             to be the optimum one, at least, for velocity contrasts of
C             1.18/1 or greater.
C             'ICB(I)' may be blank - especially in the case of a single
C             block (no structural interfaces).
C             In the presence of structural interfaces, this file is
C             recommended to be submitted in order to increase the
C             accuracy.
C             Description of input data ICB(I)
C     'TT(I)'... Name of the output file containing travel-times at all
C             network nodes after I reflections.
C             If blank, the file is not generated.
C             Description of output data TT(I)
C     'ERR(I)'... Name of the output file containing estimated upper
C             bounds for the errors of the computed travel-times at all
C             network nodes after I reflections.
C             If blank, the file is not generated.
C             Attentions concerning this version:
C             (a) in the presence of structural interfaces, the errors
C               are evaluated correctly only if the file 'ICB' is
C               submitted.
C             (b) if, in the case of reflections, 'VEL(I)' differs from
C               'VEL(I-1)', the errors cannot be evaluated correctly and
C               the corresponding 'ERR(I)' should be left blank.
C             Description of output data ERR(I)
C     'PRED(I)'... Name of the output file containing predecessors of
C             all network nodes after I reflections.
C             If blank, the file is not generated.
C             May be blank for most applications.
C             Description of output data PRED(I)
C     'NFS(I)'... Just for debugging.  The parameter should be left out
C             or blank in the input data.
C             String 'NFS(I)' may be blank, '*' or the name of the file
C             containing optimum sizes of forward stars, preceded by +
C             for input or by - for output.
C             'NFS(I)'=' ':  optimum sizes of forward stars at the
C               network nodes are estimated and used for calculation.
C             'NFS(I)'='+...' (input):  sizes of forward stars at the
C               network nodes are read from the file.
C             'NFS(I)'='-...' (output):  optimum sizes of forward stars
C               at the network nodes are estimated, used for calculation
C               and written to the file.
C             'NFS(I)'='*': sizes of all forward stars equal NFSMAX from
C               the input (2).
C             Description of data NFS(I)
C     Default: 'IND(I)'=' ', 'VEL(I)'=' ', 'ICB(I)'=' ',
C             'TT(I)'=' ','ERR(I)'=' ', 'PRED(I)'=' ', 'NFS(I)'=' '.
C TTT mode:
C     'IND(I)'... Should be blank.
C             No index file may be used in the presend version.
C     'ICB(I)'... Indices of geological blocks are read but
C             not used.  Should be blank.
C     'ERR(I)'='P1(I)'... Name of the output file containing the first
C             slowness vector component.
C             If blank, the file is not generated.
C     'PRED(I)'='P2(I)'... Name of the output file containing the second
C             slowness vector component.
C             If blank, the file is not generated.
C     'NFS(I)'='P3(I)'... Name of the output file reserved for the third
C             slowness vector component.  Presently should be blank.
C             No file is generated.
C (3.2) 'INTF(I)',/
C     'INTF(I)'... Name of the input file containing refractor points.
C             Must not be blank if the reflection is considered.
C             Description of input data INTF(I)
C Example of data set NET
C
C                                                     
C Input file 'SRC' containing source coordinates:
C (1) Several strings terminated by / (a slash).
C             The simplest way is to submit just the /.
C (2) Several times (2.1):
C (2.1) 'NAMSRC',X1S,X2S,X3S,TTS,TTSERR/
C     'NAMSRC'... Name of the source point.  Truncated to the first six
C             characters.
C     X1S,X2S,X3S... Coordinates of a point of the source.
C             In a 2-D model, coordinates perpendicular to the model
C             plane have to be the same for all source points and all
C             receivers.
C     TTS...  Initial arrival time at the source point.  Must not be
C             negative.  Zero initial time is o.k., but is changed to
C             a very small positive value.
C     TTSERR..Initial value of the arrival time error at the source
C             point.  It is likely zero at the actual source and may
C             thus be omitted.  It is introduced especially for the
C             purposes of tracing reflected waves.  In such a case,
C             TTS is the arrival time of incident wave at the reflector
C             and TTSERR is its error resulting from network ray tracing
C             between the actual source and the reflector.
C     Default: TTS=0, TTSERR=0.
C TTT mode:
C     Only a point source is considered in the present version.
C (3) / (a slash).
C Example of data set SRC
C
C                                                     
C Input file 'REC' containing receiver coordinates:
C (1) Several strings terminated by / (a slash).
C             The simplest way is to submit just the /.
C (2) Several times (2.1):
C (2.1) 'NAMREC(IR)',X1R(IR),X2R(IR),X3R(IR),/
C     'NAMREC(IR)'... Name of the receiver point.  Truncated to the
C             first six characters.
C     X1R(IR),X2R(IR),X3R(IR)... Coordinates of the IR-th receiver.
C (3) / (a slash).
C Example of data set REC
C
C                                                    
C Output file 'RAYS':
C (1) / (a slash).
C (2) For each receiver (2.1), (2.2), and (2.3):
C (2.1) 'RAYnnnnn FROM NAMSRC TO NAMREC',/
C     'RAYnnnnn FROM NAMSRC TO NAMREC'... String in apostrophes,
C             composed of:
C             (a1) substring 'RAY',
C             (a2) the sequential index of the receiver written using
C               format (i5),
C             (b1) substring ' FROM ',
C             (b2) name of the source point truncated to the first 6
C               characters,
C             (c1) substring ' TO ',
C             (c2) name of the receiver point truncated to the first 6
C               characters.
C             If the name of the source point is blank, (b1) is replaced
C               by 6 spaces.
C             If the name of the receiver point is blank, (c1) and (c2)
C               are left out.
C             If both the names of the source and receiver point are
C               blank, (b1), (b2), (c1) and (c2) are left out.
C             Examples of the string:
C               'RAY  112'
C               'RAY  112             TO NAMREC'
C               'RAY  112 FROM NAMSRC'
C               'RAY  112 FROM NAMSRC TO NAMREC'
C (2.2) For each point of the ray (2.1.1):
C (2.2.1) X1,X2,X3,TT,/
C     X1,X2,X3... Coordinates of the point of the ray.
C     TT...   Arrival time at the point.
C     Note: Rays are stored backwards, starting with the receiver,
C             ending with the source point.
C (2.3) / (a slash).
C (3) / (a slash).
C
C                                                     
C Output file 'END':
C (1) / (a slash).
C (2) For each receiver (2.1):
C (2.1) 'NAMREC',X1R,X2R,X3R,TT,TTERR,/
C     'NAMREC'... Name of the receiver point truncated to the first six
C             characters.
C     X1R,X2R,X3R... Coordinates of the receiver.
C     TT...   Arrival time at the receiver point.
C     TTERR...Estimated maximum error of the arrival time.
C             Attentions concerning this version:
C             (a) in the presence of structural interfaces, the errors
C               are evaluated correctly only if the file 'ICB' is
C               submitted.
C             (b) if, in the case of reflections, 'VEL(I)' differs from
C               'VEL(I-1)', the errors cannot be evaluated correctly and
C               the corresponding 'ERR(I)' should be left blank.
C (3) / (a slash).
C
C                                                     
C Input index files 'IND(I)':
C (1) (IND(I),I=1,N1*N2*N3)
C     IND(I)..Zero if the I-th big brick does not belong to the
C               network,
C             otherwise the index the big brick.  The nodes within the
C               big brick are indexed by integers from
C               L1*L2*L3*(IND(I)-1)+1 to L1*L2*L3*IND(I).
C     Default: IND(I)=I.
C
C                                                     
C Input files 'VEL(I)' containing velocities at network nodes:
C (1) (V(I),I=1,L1234)
C     V(I)... Velocity at network node no.I.
C     L1234.. Maximum index of a network node.  L1234=L1*L2*L3*L4, where
C             L4 is the maximum index of a big brick belonging to the
C             network, see the file 'IND' below.  If the file 'IND'
C             is not specified, L4=N1*N2*N3 by the default.
C             Free space is indicated by a zero velocity V(I)=0.
C             Free space should not be indicated by a small positive
C             velocity.  A small velocity considerably increases time of
C             computation.
C     Default: V(I)=1.
C TTT mode:
C     At structural interfaces, it is recommended to specify average
C     slownesses over bricks.  However, the method is unstable in block
C     models.
C Example of data set VEL
C
C                                                     
C Input index files 'ICB(I)':
C (1) (ICB(I),I=1,L1234)
C     ICB(I)..Index of (geological) block in which the network node no.I
C             is situated.
C     Default: ICB(I)=1.
C
C                                                      
C Output files 'TT(I)' (time field):
C (1) (TT(I),I=1,L1234)
C     TT(I)... Travel time at network node no.I.
C     L1234.. Maximum index of a network node.
C The output file is designed to be read by the list-directed input
C (free format).  The null values are generated in place of undefined
C travel times in a free space.  The null values are treated as default
C values when read by list-directed input (free format).
C Example: 124 null values are written as '    124*'.
C
C                                                     
C Output files 'ERR(I)' (time error field):
C (1) (ERR(I),I=1,L1234)
C     ERR(I)..Estimated maximum absolute value of the travel-time error
C             at network node no.I.
C             Please, notice that the absolute (i.e. not relative)
C             travel-time errors are estimated and written.
C     L1234.. Maximum index of a network node.
C The output file is designed to be read by the list-directed input
C (free format).  The null values are generated in place of undefined
C errors in a free space.
C
C                                                    
C Output files 'PRED(I)' (predecessors):
C (1) (IPRED(I),I=1,L1234)
C     IPRED(I)... Predecessor to network node no.I.
C             IPRED(I)=I: node has no predecessor (is in a free space).
C             IPRED(I)=0: node belongs to a reflector, preceding node
C               has the same position.
C             IPRED(I).LT.0: preceding node is (-IPRED(I))-th source
C               point.
C             Otherwise, IPRED(I).GT.0.AND.IPRED(I).NE.I.: preceding
C               node is a gridpoint, IPRED(I) being its position index.
C     L1234.. Maximum index of a network node.
C
C                                                     
C Output files 'NFS(I)':
C (1) (NFS(I),I=1,L1234)
C     NFS(I)..Optimum forward star size at the network node no.I.
C             If NFS.GT.NFSMAX, it is automatically decreased to NFSMAX,
C             where NFSMAX takes its positive input value or is
C             determined automatically if NFSMAX=0 on input.
C
C                                                    
C Input files 'INTF(I)' containing indices of gridpoints forming the
C reflectors:
C (1) (INTF(I),I=1,N),/
C     INTF(I)... Index of the I-th gridpoint forming the reflector.
C     N...    Number of points forming the reflector.
C
C.......................................................................
C
C List of routines:
C
C High-level:
C     NET...  Main program calling subroutines:
C             INPUT, SOURCE, GENER, RESFIL, ERRORS, RECVRS, TRACER, and
C             OUTPUT.
C             NET
C     INPUT...Reads the whole main input data file 'NET', source and
C             receiver coordinates.
C             INPUT
C     SOURCE..Reads input data for each iteration, and initializes the
C             arrival-time field and other arrays.  Under iteration we
C             understand here the computation between two reflections.
C             If reflectors are not specified, there is just one
C             iteration.
C             SOURCE
C     GENER...Performs shortest path ray tracing (one iteration).
C             GENER
C     RESFIL..In a case of reflections, stores the results of an
C             iteration in an unformatted direct-access scratch file
C             (one iteration).
C             RESFIL
C     ERRORS..Generates the field of arrival-time errors by means of
C             backward ray tracing - writes 'ERR(I)' (one iteration).
C             ERRORS
C     RECVRS..Updates arrival times at the receivers.
C             RECVRS
C     TRACER..Performs backward ray tracing from the receivers,
C             including the estimation of arrival-time errors.
C             TRACER
C     OUTPUT..Rewrites arrival-time fields and predecessors from memory
C             or scratch files to formatted output files 'TT(I)' and
C             'IPRED(I)'.
C             OUTPUT
C Updating:
C     LOOP,FREVOL,UPDATE,SRCREC... Auxiliary routines to SOURCE, GENER,
C             and RECVRS, updating a node of a network.
C             Loop with FREVOL or UPDATE: employed by GENER.
C             Loop with SRCREC: employed by SOURCE and RECVRS.
C             LOOP
C             FREVOL
C             UPDATE
C             SRCREC
C Template forward stars:
C     READFS..Called by SOURCE, to read and construct optimized template
C             forward stars.
C             READFS
C     MAKEFS..Called by SOURCE and RECVRS, to construct full spherical
C             template forward stars.
C             MAKEFS
C     SETFS...Auxiliary subroutine to READFS and MAKEFS.
C             SETFS
C     STORFS..Auxiliary subroutine to SETFS.
C             STORFS
C Backward ray tracing:
C     ONERAY..Subroutine called by ERRORS and TRACER, to perform
C             backward ray tracing.
C             ONERAY
C Accuracy estimation:
C     SETERR..Subroutine called by ONERAY, to estimate and accumulate
C             the arrival-time errors during backward ray tracing.
C             SETERR
C     OPTNFS..Subroutine called by SOURCE and SETERR, to set up
C             optimum sizes of forward stars.
C             OPTNFS
C     OPTMAT..Auxiliary subroutine to SETERR and OPTNFS, to calculate
C             the matrix composed of the first and second slowness
C             derivatives.
C             OPTMAT
C     TRYMAT..Auxiliary subroutine to OPTMAT, to try the calculation of
C             the matrix.
C             TRYMAT
C     MIXDER..Auxiliary subroutine to TRYMAT, to calculate mixed second
C             partial slowness derivatives.
C             MIXDER
C Slowness interpolation:
C     SLOW... Subroutine called by SOURCE, RECVRS, TRACER, and ONERAY,
C             to find a close gridpoint, and to interpolate the slowness
C             within the same geological block.
C             SLOW
C     POS...  Auxiliary subroutine to SLOW, to find the nearest
C             gridpoint.
C             POS
C Indexing of nodes:
C     INDX... Integer function, called by many subroutines, to return
C             the index of the node at the given gridpoint.
C             INDX
C
C.......................................................................
C
C External subroutines required and not contained within this file:
C
C Formatted output:
C     WARRAY..Auxiliary subroutine to ERRORS and OUTPUT, designed to
C             write the given array into the file.
C             Source code file 'forms.for'.
C Eigenvalues:
C     EIGEN...Subroutine from the IBM scientific subroutine package,
C             called by OPTNFS.
C             Source code file 'eigen.for'.
C
C-----------------------------------------------------------------------
C
C     
C
      PROGRAM NET
C
C-----------------------------------------------------------------------
C
C Common blocks:
      INCLUDE 'net.inc'
C     net.inc
C All common blocks are declared here.
C
C-----------------------------------------------------------------------
C
      INTEGER IREFL
*     INTEGER IT1,IT2
*     REAL TIME
C
C     IREFL.. Number of reflections.  IREFL=0: direct wave.
C     IT1,IT2,TIME... Auxiliary variables for time watching.
C
C.......................................................................
C
C     Reading the input parameters
      CALL INPUT()
C
C     Loop over reflections:
      DO 10 IREFL=0,NREFL
C
C       Generating time field of the IREFL-times reflected wave:
        CALL SOURCE(IREFL)
C
C       Example timer for Lahey Fortran77 compiler:
*       IT1=0
*       CALL TIMER(IT1)
C
C       Calculation of shortest paths and arrival times:
        CALL GENER(IREFL)
C
C       Example timer for Lahey Fortran77 compiler:
*       IT2=0
*       CALL TIMER(IT2)
*       TIME=0.01*(FLOAT(IT2)-FLOAT(IT1))/60.
*       WRITE(*,'(A,F10.3,A)') ' Time=',TIME,' min.'
C
C       Post-processing:
        CALL RESFIL(IREFL)
        CALL ERRORS(IREFL)
C
   10 CONTINUE
C     End of loop over reflections.
C
C     End of ray computation and writing the results:
      CALL RECVRS()
      CALL TRACER()
      CALL OUTPUT()
      STOP
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE INPUT()
C
C Procedure to read the whole main input data file 'net', source and
C receiver coordinates.
C
C-----------------------------------------------------------------------
C
C     Common blocks /GRID/ /FS/ /NAMESR/ /POINTS/ /FILES/ /TTTPAR/ are
C     required here:
      INCLUDE 'net.inc'
C     net.inc
      INCLUDE 'ttt.inc'
C     ttt.inc
C
C-----------------------------------------------------------------------
C
      INTEGER LU1,LU2
      REAL DWARF,UNDEF
      PARAMETER(LU1=1)
      PARAMETER (DWARF=1.E-10,UNDEF=-999999.)
C
      CHARACTER*80 TEXT
      INTEGER ISRC,IREC,I,L
      REAL D1,D2,D3,O1,O2,O3
C     REAL A1111,A1122,A2222,A1133,A2233,A3333,A1123
C     REAL A2223,A3323,A2323,A1113,A2213,A3313,A2313
C     REAL A1313,A1112,A2212,A3312,A2312,A1312,A1212
C
C     TEXT...Names of input files NET and SEP, then a dummy storage
C            location for a text.
C     ID(II),JD(II),KD(II)... Vector specifying II-th node of the
C            forward star, in dimensions of a small brick.
C     A1111,A1122,A2222,A1133,A2233,A3333,A1123,A2223,A3323,A2323,A1113,
C            A2213,A3313,A2313,A1313,A1112,A2212,A3312,A2312,A1312,
C            A1212... Intended for future extension (anisotropy).
C
C.......................................................................
C
C     Name of the main input data file is read from the * external unit:
      TEXT='net.h'
      WRITE(*,'(A)')
     *          '+NET: Enter name of main input data file [''net.h'']: '
      READ(*,*) TEXT
C     End of reading from the * unit.
C
C     Reading SEP parameter file:
      CALL RSEP1(LU1,TEXT)
C
C     Reading input data file NET:
      CALL RSEP3T('NET',TEXT,'net.dat')
      OPEN(LU1,FILE=TEXT,STATUS='OLD')
      WRITE(*,'(2A)') '+NET: Reading input data file: ',TEXT(1:36)
C
C     Numbers of gridpoints:
      CALL RSEP3I('N1',N1,1)
      CALL RSEP3I('N2',N2,1)
      CALL RSEP3I('N3',N3,1)
      CALL RSEP3I('L1',L1,1)
      CALL RSEP3I('L2',L2,1)
      CALL RSEP3I('L3',L3,1)
      IF(N1.LE.0) THEN
        LN1=.TRUE.
        N1=1
        L1=MAX0(1,L1)
      ELSE
        LN1=.FALSE.
      END IF
      IF(N1.LT.1.OR.N2.LT.1.OR.N3.LT.1.OR.
     *   L1.LT.1.OR.L2.LT.1.OR.L3.LT.1) THEN
C       NET-10
        CALL ERROR('NET-10: Number of gridpoints is not positive')
C       NET-10: Number of gridpoints is not positive:
C               N1,N2,N3, and, if specified, also L1,L2,L3, in input
C               file SEP must be positive.
      END IF
      NL1=N1*L1
      NL2=N2*L2
      NL3=N3*L3
C
C     Boundaries of the model volume:
      CALL RSEP3R('D1',D1,1.)
      CALL RSEP3R('D2',D2,1.)
      CALL RSEP3R('D3',D3,1.)
      CALL RSEP3R('O1',O1,0.)
      CALL RSEP3R('O2',O2,0.)
      CALL RSEP3R('O3',O3,0.)
      IF(D1.LE.0..OR.D2.LE.0..OR.D3.LE.0.) THEN
C       NET-56
        CALL ERROR('NET-56: Grid interval D1, D2 or D3 is not positive')
C       NET-56: Grid interval D1, D2 or D3 is not positive:
C               Grid interval D1, D2 and D3 must be positive in this
C               version of the NET program, see the
C               description of input file SEP.
      END IF
      X1MIN=O1-0.5*D1
      X2MIN=O2-0.5*D2
      X3MIN=O3-0.5*D3
      X1MAX=X1MIN+FLOAT(N1)*D1
      X2MAX=X2MIN+FLOAT(N2)*D2
      X3MAX=X3MIN+FLOAT(N3)*D3
C     Space steps in the rectangular grid
      IF(LN1) THEN
        DSX1=(X1MAX-X1MIN)/FLOAT(NL2)
      ELSE IF(NL1.EQ.1) THEN
        DSX1=0.
      ELSE
        DSX1=(X1MAX-X1MIN)/FLOAT(NL1)
      END IF
      IF(NL2.EQ.1) THEN
        DSX2=0.
      ELSE
        DSX2=(X2MAX-X2MIN)/FLOAT(NL2)
      END IF
      IF(NL3.EQ.1) THEN
        DSX3=0.
      ELSE
        DSX3=(X3MAX-X3MIN)/FLOAT(NL3)
      END IF
      IF(LN1) THEN
        ASX1=0.
        ASX2=SQRT(DSX1**2+DSX2**2)
      ELSE
        ASX1=ABS(DSX1)
        ASX2=ABS(DSX2)
      ENDIF
      ASX3=ABS(DSX3)
C
C     Size of the forward star and the interpolation method for TTT:
C     Size of the forward star:
      CALL RSEP3I('NFSMAX',NFSMAX,0)
      IF(NFSMAX.GT.MFSMAX) THEN
C       NET-07
        CALL ERROR('NET-07: Maximum size of a forward star too great')
C       NET-07: Maximum size of a forward star too great:
C               NFSMAX in input data file SEP should
C               be decreased if positive, adjusted to MFSMAX if zero,
C               or the dimension MFSMAX in the common block /FS/
C               in net.inc should be increased.
      END IF
C     Interpolation method:
      CALL RSEP3R('RIDGE1',RIDGE1,1.5  )
      CALL RSEP3R('RIDGE2',RIDGE2,0.10 )
      CALL RSEP3R('VER1'  ,VER1  ,9.999)
      CALL RSEP3R('VER2'  ,VER2  ,1.000)
C
C     (1) Names of the files with source, receivers, rays, and errors:
      FREC=' '
      FRAYS=' '
      FEND=' '
      READ(LU1,*) FSRC,FREC,FRAYS,FEND
      IF(FREC.EQ.' '.AND.(FRAYS.NE.' '.OR.FEND.NE.' ')) THEN
C       NET-12
        CALL ERROR
     *       ('NET-12: Name of input file with receivers not specified')
C       NET-12: Name of input file with receivers not specified:
C               If you want to generate output file with rays, or with
C               arrival times and their errors at the receivers,
C               you have to specify the name of input file with
C               receivers.
C       See the main input file NET, (2).
      END IF
      IF(FREC.NE.' '.AND.FRAYS.EQ.' '.AND.FEND.EQ.' ') THEN
C       NET-13
        CALL ERROR
     *           ('NET-13: Name of output file with rays not specified')
C       NET-13: Name of output file with rays not specified:
C               If you specify receiver positions at the receiver
C               file, you should specify the name of the output file
C               with rays, or with arrival times and their errors at
C               the receivers.
C               See the main input file NET, (2).
      END IF
C
C     (2) Number of reflections:
      NREFL=0
      READ(LU1,*) NREFL
      IF(NREFL.GT.MREFL) THEN
C       NET-14
        CALL ERROR('NET-14: Too many reflections')
C       NET-14: Too many reflections:
C               Number NREFL (main input file NET,
C               (2)) should be decreased to satisfy the inequality
C               NREFL.LE.MREFL,
C               or the dimension MREFL in the common block /FILES/
C               in net.inc should be increased.
      END IF
C
C     (3) Names of the output travel-time and predecessor files,
C     input velocity and index files, and input refractor-point files:
      DO 32 L=0,NREFL
        IF(L.GT.0) THEN
          FINTF(L)=' '
          READ(LU1,*) FINTF(L)
          IF(FINTF(L).EQ.' ') THEN
C           NET-15
            CALL ERROR('NET-15: Reflector file not specified')
C           NET-15: Reflector file not specified:
C                   If you want to compute reflected time field, you
C                   must specify name of reflector file.
C                   See main input file NET, (4.2).
          END IF
        END IF
        FIND(L)=' '
        FVEL(L)=' '
        FICB(L)=' '
        FTT(L)=' '
        FERR(L)=' '
        FPRED(L)=' '
        FNFS(L)=' '
        READ(LU1,*)
     *           FIND(L),FVEL(L),FICB(L),FTT(L),FERR(L),FPRED(L),FNFS(L)
        IF(FVEL(L).EQ.' ') THEN
C         NET-16
          CALL ERROR('NET-16: Velocity file not specified')
C         NET-16: Velocity file not specified:
C                 Velocity file has a fundamental importance for time
C                 field computations. You must always specify name of
C                 velocity file.
C                 See main input file NET, (4.1).
        END IF
   32 CONTINUE
C
C     End of reading main input data file 'NET'
      CLOSE(LU1)
C
C.......................................................................
C
C     Reading source coordinates:
C
      WRITE(*,'(A)')
     *  '+NET: Reading source coordinates                              '
      OPEN(LU1,FILE=FSRC,STATUS='OLD')
      READ(LU1,*) (TEXT,I=1,20)
      DO 48 ISRC=1,MSRC
        X1S(ISRC)=UNDEF
        X2S(ISRC)=UNDEF
        X3S(ISRC)=UNDEF
        TTS(ISRC)=0.
        TTSERR(ISRC)=0.
        READ(LU1,*)
     *        NAMES,X1S(ISRC),X2S(ISRC),X3S(ISRC),TTS(ISRC),TTSERR(ISRC)
        IF(X1S(ISRC).EQ.UNDEF.AND.X2S(ISRC).EQ.UNDEF.AND.
     *                            X3S(ISRC).EQ.UNDEF) THEN
          NSRC=ISRC-1
          GO TO 49
        END IF
        IF(X1S(ISRC).EQ.UNDEF) X1S(ISRC)=0.
        IF(X2S(ISRC).EQ.UNDEF) X2S(ISRC)=0.
        IF(X3S(ISRC).EQ.UNDEF) X3S(ISRC)=0.
C       Checking the initial arrival time at source points
        IF(TTS(ISRC).LT.0.) THEN
C         NET-17
          CALL ERROR
     *            ('NET-17: Negative time at the source not acceptable')
C         NET-17: Negative time at the source not acceptable:
C                 Initial time at the source points must not be negative
C                 for a correct function of this code.
        ELSE IF(TTS(ISRC).EQ.0.) THEN
          TTS(ISRC)=DWARF
        END IF
C       Checking source positions
        IF(ISRC.EQ.1) THEN
          IF(NL1.EQ.1.AND..NOT.LN1) THEN
            X1MIN=X1S(1)
            X1MAX=X1S(1)
          END IF
          IF(NL2.EQ.1) THEN
            X2MIN=X2S(1)
            X2MAX=X2S(1)
          END IF
          IF(NL3.EQ.1) THEN
            X3MIN=X3S(1)
            X3MAX=X3S(1)
          END IF
        ELSE
          IF(NL1.EQ.1.AND..NOT.LN1) THEN
            IF(X1S(L).NE.X1MIN) THEN
C             NET-18
              CALL ERROR('NET-18: Different 1-st source coordinates')
C             NET-18: Different 1-st source coordinates:
C                     In a 2-D model, coordinates perpendicular to the
C                     model plane have to be the same for all source
C                     points and all receivers.
            END IF
          END IF
          IF(NL2.EQ.1) THEN
            IF(X2S(L).NE.X2MIN) THEN
C             NET-19
              CALL ERROR('NET-19: Different 2-nd source coordinates')
C             NET-19: Different 2-nd source coordinates:
C                     In a 2-D model, coordinates perpendicular to the
C                     model plane have to be the same for all source
C                     points and all receivers.
            END IF
          END IF
          IF(NL3.EQ.1) THEN
            IF(X3S(L).NE.X3MIN) THEN
C             NET-20
              CALL ERROR('NET-20: Different 3-rd source coordinates')
C             NET-20: Different 3-rd source coordinates:
C                     In a 2-D model, coordinates perpendicular to the
C                     model plane have to be the same for all source
C                     points and all receivers.
            END IF
          END IF
        END IF
        IF(X1S(ISRC).LT.X1MIN.OR.X1MAX.LT.X1S(ISRC).OR.
     *     X2S(ISRC).LT.X2MIN.OR.X2MAX.LT.X2S(ISRC).OR.
     *     X3S(ISRC).LT.X3MIN.OR.X3MAX.LT.X3S(ISRC)) THEN
C         NET-21
          CALL ERROR('NET-21: Source is not inside the model')
C         NET-21: Source is not inside the model:
C                 Source coordinates X1S,X2S,X3S must satisfy conditions
C                 X1MIN.LE.X1S.AND.X1S.LE.X1MAX,
C                 X2MIN.LE.X2S.AND.X2S.LE.X2MAX,
C                 X3MIN.LE.X3S.AND.X3S.LE.X3MAX.
C                 See the main input file NET, (2),
C                 and input file SRC.
        ENDIF
   48 CONTINUE
C     NET-22
      CALL ERROR('NET-22: Too many source points')
C     NET-22: Too many source points:
C             Number of receivers in input file REC
C             should be less than MREC,
C             or the dimension MREC in the common block /POINTS/
C             should be increased.
   49 CONTINUE
      CLOSE(LU1)
      IF(NSRC.LE.0) THEN
C       NET-23
        CALL ERROR('NET-23: No source point given')
C       NET-23: No source point given:
C               You must specify at least one source point at the
C               input source file SRC.
      END IF
C
C.......................................................................
C
C     Reading receiver coordinates:
C
      WRITE(*,'(A)')
     *  '+NET: Reading receiver coordinates                            '
      IF(FREC.EQ.' '.OR.FRAYS.EQ.' ') THEN
        NREC=0
      ELSE
        OPEN(LU1,FILE=FREC,STATUS='OLD')
        READ(LU1,*) (TEXT,I=1,20)
        DO 58 IREC=1,MREC
          X1R(IREC)=UNDEF
          X2R(IREC)=UNDEF
          X3R(IREC)=UNDEF
          READ(LU1,*) NAMER(IREC),X1R(IREC),X2R(IREC),X3R(IREC)
          IF(X1R(IREC).EQ.UNDEF.AND.X2R(IREC).EQ.UNDEF.AND.
     *                              X3R(IREC).EQ.UNDEF) THEN
            NREC=IREC-1
            GO TO 59
          END IF
          IF(X1R(IREC).EQ.UNDEF) X1R(IREC)=0.
          IF(X2R(IREC).EQ.UNDEF) X2R(IREC)=0.
          IF(X3R(IREC).EQ.UNDEF) X3R(IREC)=0.
C         Checking the receiver positions
          IF(NL1.EQ.1.AND..NOT.LN1) THEN
            IF(X1R(IREC).NE.X1MIN) THEN
C             NET-24
              CALL ERROR
     *         ('NET-24: Different 1-st source and receiver coordinate')
C             NET-24: Different 1-st source and receiver coordinate:
C                     In a 2-D model, coordinates perpendicular to the
C                     model plane have to be the same for all source
C                     points and all receivers.
            END IF
          END IF
          IF(NL2.EQ.1) THEN
            IF(X2R(IREC).NE.X2MIN) THEN
C             NET-25
              CALL ERROR
     *         ('NET-25: Different 2-nd source and receiver coordinate')
C             NET-25: Different 2-nd source and receiver coordinate:
C                     In a 2-D model, coordinates perpendicular to the
C                     model plane have to be the same for all source
C                     points and all receivers.
            END IF
          END IF
          IF(NL3.EQ.1) THEN
            IF(X3R(IREC).NE.X3MIN) THEN
C             NET-26
              CALL ERROR
     *         ('NET-26: Different 3-rd source and receiver coordinate')
C             NET-26: Different 3-rd source and receiver coordinate:
C                     In a 2-D model, coordinates perpendicular to the
C                     model plane have to be the same for all source
C                     points and all receivers.
            END IF
          END IF
          IF(X1R(IREC).LT.X1MIN.OR.X1MAX.LT.X1R(IREC).OR.
     *       X2R(IREC).LT.X2MIN.OR.X2MAX.LT.X2R(IREC).OR.
     *       X3R(IREC).LT.X3MIN.OR.X3MAX.LT.X3R(IREC)) THEN
C           NET-27
            CALL ERROR('NET-27: Receiver is not inside the model')
C           NET-27: Receiver is not inside the model:
C                   Receiver coordinates X1R,X2R,X3R must satisfy
C                   conditions
C                   X1MIN.LE.X1R.AND.X1R.LE.X1MAX,
C                   X2MIN.LE.X2R.AND.X2R.LE.X2MAX,
C                   X3MIN.LE.X3R.AND.X3R.LE.X3MAX.
C                   See the main input file NET, (2),
C                   and input file REC.
          ENDIF
   58   CONTINUE
C       NET-28
        CALL ERROR('NET-28: Too many receivers')
C       NET-28: Too many receivers:
C               Number of receivers in the input file
c               REC should be less than MREC,
C               or the dimension MREC in the common block /POINTS/
C               in net.inc should be increased.
   59   CONTINUE
        CLOSE(LU1)
      END IF
C
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE SOURCE(IREFL)
      INTEGER IREFL
C
C Initialization procedure for starting the program from the source or
C from the IREFL-th interface.  Reads input data for the iteration,
C and initializes the arrival-time field and other arrays.
C
C Input:
C     IREFL.. Number of reflections.
C
C No output.
C
C-----------------------------------------------------------------------
C
C     Common blocks /GRID/ /FS/ /NAMESR/ /POINTS/ /FILES/ /NODE/ /SRCC/
C     are required here:
      INCLUDE 'net.inc'
C     net.inc
      INCLUDE 'netnode.inc'
C     netnode.inc
C
C-----------------------------------------------------------------------
C
      EXTERNAL SRCREC,TTTSRC,INDX
      INTEGER INDX
C
      INTEGER LU1
      PARAMETER (LU1=1)
      REAL GIANT
      PARAMETER (GIANT=1.0E+10)
C
      INTEGER MIND,MGRID,MP2,MP3,MPRED,MNFS,MICB
      INTEGER N123,L123,ISRC,INDQ,NFSMAK,I,L,M
      INTEGER JPOS,JPOS1,JPOS2,JPOS3,JADR,NFSOPT,JCOUNT,NUMOPT
      REAL VMAX,V,C1,C2MIN,C2MAX,AUX,SUMOPT
C
C.......................................................................
C
      WRITE(*,'(A)') '+NET: Reading velocities                         '
C
C     Reading the index array:
      IF(IREFL.GT.0) THEN
        IF(FIND(IREFL).EQ.FIND(IREFL-1)) THEN
C         The index array has not been changed
          GO TO 20
        END IF
      END IF
      L4=0
      IND0=0
C
      IF(FIND(IREFL).EQ.' ') THEN
        L1234=N1*N2*N3
        IF(L1.GT.1.OR.L2.GT.1.OR.L3.GT.1) THEN
C         NET-29
          CALL ERROR('NET-29: No index file specified')
C         NET-29: No index file specified:
C                 If L1,L2,L3 (input data file SEP)
C                 are specified, the index file IND has to be
C                 submitted. See the main input file
c                 NET, (4.1).
        END IF
        MIND=0
      ELSE
        MIND=MRAM
        N123=N1*N2*N3
        L123=L1*L2*L3
        IF(N123.GT.MIND) THEN
C         NET-30
          CALL ERROR('NET-30: Too many big bricks')
C         NET-30: Too many big bricks:
C                 Numbers N1,N2,N3 (input data SEP)
C                 should be decreased to satisfy the inequality
C                 N1*N2*N3.LE.MRAM,
C                 or the dimension MRAM of array RAM in include file
C                 ram.inc should be
C                 increased.
        END IF
        DO 11 L=1,N123
          IRAM(IND0+L)=L
   11   CONTINUE
        CALL RARRAI(LU1,FIND(IREFL),'FORMATTED',
     *              .FALSE.,0,N123,IRAM(IND0+1))
        DO 12 L=1,N123
          L4=MAX0(IRAM(IND0+L),L4)
          IRAM(IND0+L)=(IRAM(IND0+L)-1)*L123+1
   12   CONTINUE
        L1234=L1*L2*L3*L4
        MIND=N123
      END IF
C
C     Dynamic array allocation:
      MGRID=L1234
      IF(NFSMAX.GE.0) THEN
C       Network ray tracing:
        IF(FICB(IREFL).EQ.' ') THEN
          MICB=0
        ELSE
          MICB=L1234
        END IF
        IF(FNFS(IREFL).EQ.'*') THEN
          MNFS=0
        ELSE
          MNFS=L1234
        END IF
        IF(FPRED(IREFL).EQ.' '.AND.
     *     FERR(IREFL).EQ.' '.AND.FRAYS.EQ.' '.AND.FEND.EQ.' ') THEN
          MPRED=0
        ELSE
          MPRED=L1234
        END IF
        MP2=0
        MP3=0
      ELSE
C       Second-order grid travel-time tracing:
        MICB=0
        MNFS=0
        MPRED=0
        MP2=L1234
        MP3=0
      END IF
      ITT0  =IND0  +MIND
      IPOSQ0=ITT0  +MGRID
      IP0   =IPOSQ0+MGRID
      IP20  =IP0   +MGRID
      IP30  =IP20  +MP2
      IPRED0=IP30  +MP3
      NFS0  =IPRED0+MPRED
      ICB0  =NFS0  +MNFS
      I     =ICB0  +MICB
      IF(I.GT.MRAM) THEN
C       NET-01
        CALL ERROR('NET-01: Small array RAM')
C       NET-01: Small array RAM:
C               Dimension MRAM of array RAM allocated in
C               ram.inc
C               must be at least
C                 MIND+L1*L2*L3*L4*(4+M1PRED+M1NFS+M1ICB),
C               where:
C                 MIND=0  if 'IND'=' ', see 'index file' (usually
C                   applicable),
C                 MIND=N1*N2*N3  otherwise.
C                 L4      is the number of nonempty big bricks.
C                 L4=N1*N2*N3 and L1=1,L2=1,L3=1 without volume indexing
C                   ('IND'=' ').
C                 N1,N2,N3,L1,L2,L3 are given by input data file
C                   SEP.
C                 M1NFS=-1 without the optimization of sizes of forward
C                   stars,
C                 M1NFS=0 with the optimization of sizes of forward
C                   stars and also for TTT,
C                 M1PRED=0 if travel-time errors are not calculated and
C                         the predecessor file is not generated and also
C                         for TTT,
C                 M1PRED=1 otherwise.
C                 M1ICB=0 in models without structural interfaces and
C                   also for TTT,
C                 M1ICB=1 in models with structural interfaces.
C               N1,N2,N3 (input data file SEP) should
C               be decreased to satisfy the inequality N1*N2*N3.LE.MRAM,
C               or the dimension MRAM of array RAM in include file
C               ram.inc should be
C               increased.
      END IF
C
C     Reading the velocity field
   20 CONTINUE
      IF(IREFL.GT.0) THEN
        IF(FVEL(IREFL).EQ.FVEL(IREFL-1)) THEN
C         The velocity field has not been changed
          GO TO 30
        END IF
      END IF
      IF(L1234.GT.MGRID) THEN
C       NET-31
        CALL ERROR('NET-31: Too many network nodes')
C       NET-31: Too many network nodes:
C               The number of network nodes, see input file
C               SEP, should not exceed MGRID.
C               This error should not appear.  Contact the authors.
      END IF
      IF(IREFL.GT.0) THEN
        IF(FERR(IREFL).NE.' ') THEN
C         NET-43
          CALL ERROR
     *   ('NET-43: Time errors cannot be computed after reflection')
C         NET-43: Time errors cannot be computed after reflection:
C                 If the velocity field of the reflected wave is
C                 different from the velocity field before the
C                 reflection, the errors of the arrival times cannot be
C                 computed correctly by this version.  When the rays are
C                 traced, only time fields and predecessors are
C                 alternated in the memory, unlike slowness fields, see
C                 subroutine SETERR.
C                 Thus, if the filenames VEL(I) at the lines (4.1) of
C                 main input data NET differ, the
C                 filenames ERR(I)
C                 should not be specified at the lines (4.1) for the
C                 reflected waves.  Simultaneously, if the filename
C                 END
C                 at the line (2) is specified, errors written in the
C                 output file END are incorrect.
        END IF
        IF(FEND.NE.' ') THEN
C         NET-44
          CALL ERROR
     *  ('NET-44: Time errors cannot be computed after reflection.')
C         NET-44: Time errors cannot be computed after reflection:
C                 If, in a case of reflection, filenames VEL(I) at the
C                 lines (4.1) of main input data NET
c                 differ and the filename END at line (2) is specified,
C                 errors written in the output file END are incorrect.
C                 For more details see also error NET-43 above.
        END IF
      END IF
      DO 21 L=1,L1234
        RAM(IP0+L)=1.
   21 CONTINUE
      CALL RARRAY(LU1,FVEL(IREFL),'FORMATTED',
     *            .FALSE.,0.,L1234,RAM(IP0+1))
      VMAX=0.
      AUX=1.001/GIANT
      DO 22 L=1,L1234
        V=RAM(IP0+L)
        IF(ABS(V).LT.AUX) THEN
          RAM(IP0+L)=GIANT
        ELSE
          VMAX=AMAX1(V,VMAX)
          RAM(IP0+L)=1.0/V
        ENDIF
   22 CONTINUE
C
C     Reading the indices of blocks:
   30 CONTINUE
      IF(IREFL.GT.0) THEN
        IF(FICB(IREFL).EQ.FICB(IREFL-1)) THEN
C         The indices of blocks have not been changed
          GO TO 40
        END IF
      END IF
      IF(FICB(IREFL).EQ.' ') THEN
        LICB=.FALSE.
      ELSE
        LICB=.TRUE.
        IF(L1234.GT.MICB) THEN
C         NET-32
          CALL ERROR
     *      ('NET-32: Insufficient memory for indices of blocks')
C         NET-32: Insufficient memory for indices of blocks:
C                 The number of network nodes, see input file
C                 SEP, should not exceed MICB.
C                 This error should not appear.  Contact the authors.
        END IF
        DO 31 L=1,L1234
          IRAM(ICB0+L)=1
   31   CONTINUE
        CALL RARRAI(LU1,FICB(IREFL),'FORMATTED',
     *              .FALSE.,0,L1234,IRAM(ICB0+1))
      END IF
C
C     Switch for storing predecessors:
   40 CONTINUE
      IF(NFSMAX.LT.0) THEN
C       Second-order grid travel-time tracing (array for P2):
        LPRED=.FALSE.
        IF(L1234.GT.MP2) THEN
C         NET-52
          CALL ERROR('NET-52: Insufficient memory for P2')
C         NET-52: Insufficient memory for P2:
C                 This error should not appear.  Contact the authors.
        END IF
        DO 42 L=1,L1234
          RAM(IP20+L)=0.
   42   CONTINUE
      ELSE IF(FREC.EQ.' '.AND.FPRED(IREFL).EQ.' '
     *                   .AND.FERR (IREFL).EQ.' ') THEN
        LPRED=.FALSE.
      ELSE
        LPRED=.TRUE.
        IF(L1234.GT.MPRED) THEN
C         NET-33
          CALL ERROR('NET-33: Insufficient memory for predecessors')
C         NET-33: Insufficient memory for predecessors:
C                 The number of network nodes, see input file
C                 SEP, should not exceed MPRED.
C                 This error should not appear.  Contact the authors.
        END IF
      END IF
C
C     Computing the forward star sizes:
      IF(NFSMAX.LT.0) THEN
C       Second-order grid travel-time tracing (array for P3):
        LNFS=.FALSE.
        IF(NL3.GT.1.AND.L1234.GT.MP3) THEN
C         NET-53
          CALL ERROR('NET-53: Insufficient memory for P3')
C         NET-53: Insufficient memory for P3:
C                 This error should not appear.  Contact the authors.
        END IF
C       DO 43 L=1,L1234
C         RAM(IP30+L)=0.
C  43   CONTINUE
      ELSE IF(FNFS(IREFL)(1:1).EQ.'*') THEN
        LNFS=.FALSE.
        IF(NFSMAX.LE.0) THEN
C         NET-08
          CALL ERROR('NET-08: Zero maximum size of a forward star')
C         NET-08: Zero maximum size of a forward star:
C                 For NFSMAX=0 in input data SEP,
C                 option 'NFS(I)'='*' cannot be used in input data
C                 NET in record (4.1).
        END IF
      ELSE
        LNFS=.TRUE.
        IF(L1234.GT.MNFS) THEN
C         NET-34
          CALL ERROR
     *     ('NET-34: Insufficient memory for forward star sizes')
C         NET-34: Insufficient memory for forward star sizes:
C                 The number of network nodes, see input file
C                 SEP, should not exceed MNFS.
C                 This error should not appear.  Contact the authors.
        END IF
        IF(FNFS(IREFL)(1:1).EQ.'+') THEN
          DO 51 L=1,L1234
            IRAM(NFS0+L)=NFSMAX
   51     CONTINUE
          CALL RARRAI(LU1,FNFS(IREFL)(2:LEN(FNFS(IREFL))),'FORMATTED',
     *                .FALSE.,0,L1234,IRAM(NFS0+1))
          IF(NFSMAX.EQ.0) THEN
            DO 52 L=1,L1234
              IF(IRAM(NFS0+L).GT.NFSMAX) THEN
                NFSMAX=IRAM(NFS0+L)
              ENDIF
   52       CONTINUE
          ELSE
            DO 53 L=1,L1234
              IF(IRAM(NFS0+L).GT.NFSMAX) THEN
                IRAM(NFS0+L)=NFSMAX
              ENDIF
   53       CONTINUE
          END IF
        ELSE
          JCOUNT=0
          DO 58 JPOS3=0,NL3-1
            DO 57 JPOS2=0,NL2-1
              JPOS=NL1*(JPOS2+NL2*JPOS3)
              DO 56 JPOS1=0,NL1-1
                JPOS=JPOS+1
C               JPOS=1+JPOS1+NL1*(JPOS2+NL2*JPOS3)
                JADR=INDX(JPOS)
                IF(JADR.GT.0) THEN
                  CALL OPTNFS
     *               (JPOS,JPOS1,JPOS2,JPOS3,JADR,NFSOPT,C1,C2MIN,C2MAX)
                  IRAM(NFS0+JADR)=NFSOPT
                  JCOUNT=JCOUNT+1
                  IF(JCOUNT/1000*1000.EQ.JCOUNT) THEN
                    WRITE(*,'(A,I7,A,I7)')
     *               '+NET: Generating optimum sizes of forward stars:',
     *                    JCOUNT,' of',L1234
                  END IF
                END IF
   56         CONTINUE
   57       CONTINUE
   58     CONTINUE
          IF(NFSMAX.EQ.0) THEN
C           Automatic estimation of NFSMAX:
            NUMOPT=0
            SUMOPT=0.
            DO 61 L=1,L1234
              I=IRAM(NFS0+L)
              IF(I.GT.0) THEN
                NUMOPT=NUMOPT+1
                SUMOPT=SUMOPT+1./FLOAT(I)**2
              END IF
   61       CONTINUE
            SUMOPT=SQRT(FLOAT(NUMOPT)/SUMOPT)
            NFSMAX=INT(SUMOPT+0.5)
            NFSMAX=MIN0(NFSMAX,MAX0(N1*L1,N2*L2,N3*L3)-1)
            IF(NFSMAX.LE.0) THEN
C             NET-09
              CALL ERROR('NET-09: Zero maximum size of a forward star')
C             NET-09: Zero maximum size of a forward star:
C                     Contact the authors.
            END IF
            IF(NFSMAX.GT.MFSMAX) THEN
C             NET-11
              CALL ERROR
     *              ('NET-11: Maximum size of a forward star too great')
C             NET-11: Maximum size of a forward star too great:
C                     NFSMAX in input data file SEP
C                     should be decreased if positive,
C                     adjusted to MFSMAX if zero,
C                     or the dimension MFSMAX in the common block /FS/
C                     in net.inc should be increased.
            END IF
            DO 62 L=1,L1234
              IF(IRAM(NFS0+L).GT.NFSMAX) THEN
                IRAM(NFS0+L)=NFSMAX
              END IF
   62       CONTINUE
          END IF
          WRITE(*,'(A,I7,A,I7)')
     *               '+NET: Generating optimum sizes of forward stars:',
     *                    JCOUNT,' of',L1234
          IF(FNFS(IREFL).NE.' ') THEN
            CALL WARRAI(LU1,FNFS(IREFL)(2:LEN(FNFS(IREFL))),'FORMATTED',
     *                  .FALSE.,0,.FALSE.,0,L1234,IRAM(NFS0+1))
          END IF
        END IF
      END IF
C
C     2-D or 3-D model:
      IF(NL3.EQ.1) THEN
        IF(NL2.EQ.1) THEN
          IF(NL1.EQ.1) THEN
C           NET-35
            CALL ERROR('NET-35: One-point model')
C           NET-35: One-point model:
C                   N1*L1=N2*L2=N3*L3=1 in input file
C                   SEP.
          ELSE
C           NET-36
            CALL ERROR('NET-36: Line model')
C           NET-36: Line model:
C                   N2*L2=N3*L3=1 in input file SEP.
          END IF
        ELSE
          IF(NL1.EQ.1) THEN
C           NET-37
            CALL ERROR('NET-37: Line model')
C           NET-37: Line model:
C                   N1*L1=N3*L3=1 in input file SEP.
          ELSE
            DMIN=AMIN1(ASX1,ASX2)/VMAX
          END IF
        END IF
      ELSE
        IF(NL2.EQ.1) THEN
          IF(NL1.EQ.1) THEN
C           NET-38
            CALL ERROR('NET-38: Line model')
C           NET-38: Line model:
C                   N1*L1=N2*L2=1 in input file SEP.
          ELSE
            DMIN=AMIN1(ASX1,ASX3)/VMAX
          END IF
        ELSE
          IF(NL1.EQ.1) THEN
            DMIN=AMIN1(ASX2,ASX3)/VMAX
          ELSE
            DMIN=AMIN1(ASX1,ASX2,ASX3)/VMAX
          END IF
        END IF
      END IF
      IF(NFSMAX.LT.0) THEN
C       Second-order grid travel-time tracing:
        AUX=ASX1*ASX1+ASX2*ASX2+ASX3*ASX3
        DMIN=SQRT(AUX)-SQRT(AUX-DMIN*DMIN)
      END IF
C
C.......................................................................
C
C     Writing parameters on the screen:
      IF(IREFL.EQ.0) THEN
        WRITE(*,'(A,I2,A,3(I4,A),I9,A,I2,A,I4,A)')
     *  '+',NREFL,' reflections,',NL1,'*',NL2,'*',NL3,'=',NL1*NL2*NL3,
     *  ' gridpoints, f.s.size:',NFSMAX,',',NREC,' receivers'
        WRITE(*,*)
      END IF
C
C.......................................................................
C
C     Initialization:
C
      WRITE(*,'(A)')
     *  '+NET: Initializing source nodes                               '
C
      IF(IREFL.LE.0) THEN
C
C       Given source points:
C
C       Initialization of arrays
        DO 71 I=1,NL1*NL2*NL3
          IADR=INDX(I)
          IF(IADR.GT.0) THEN
            RAM(ITT0+IADR)=GIANT
            IF(LPRED) THEN
              IRAM(IPRED0+IADR)=I
            END IF
          END IF
   71   CONTINUE
C
C       MINQ,MAXQ describe the extent of the queue.
        MINQ=1
        MAXQ=0
        TTMIN=TTS(1)
C
C       Extent of a dense forward star at source points:
        IF(LNFS) THEN
          NFSMAK=0
        ELSE IF(NFSMAX.GE.0) THEN
          NFSMAK=NFSMAX
          CALL MAKEFS(NFSMAK)
        END IF
C
C       Loop over source points
        DO 79 ISRC=1,NSRC
C
C         Minimum first-arrival time:
          TTMIN=AMIN1(TTS(ISRC),TTMIN)
C
C         Parameters of the source point:
          CALL SLOW(X1S(ISRC),X2S(ISRC),X3S(ISRC),DPOS1,DPOS2,DPOS3,
     *                                   IPOS1,IPOS2,IPOS3,IPOS,IADR,PI)
          TTI=TTS(ISRC)
          ISRCI=-ISRC
C
C         Updating:
          IF(NFSMAX.GE.0) THEN
C           Network ray tracing:
C           Adjusting extent of dense forward star at the source point
            IF(LNFS) THEN
              IF(IRAM(NFS0+IADR).NE.NFSMAK) THEN
                NFSMAK=IRAM(NFS0+IADR)
                CALL MAKEFS(NFSMAK)
              END IF
            END IF
C           Updating
            CALL LOOP(SRCREC)
          ELSE
C           Second-order grid travel-time tracing:
            CALL TTTSRC()
          END IF
C
   79   CONTINUE
C
C       Reading and assembling optimized forward star:
        CALL READFS()
C
C.......................................................................
C
      ELSE
C
C       Given reflector:
C
C       Reading the interface file:
C       Creating queue for travel-times on the interface
        DO 81 INDQ=1,MGRID
          IRAM(IPOSQ0+INDQ)=0
   81   CONTINUE
        CALL RARRAI(LU1,FINTF(IREFL),'FORMATTED',
     *              .FALSE.,0,L1234,IRAM(IPOSQ0+1))
        DO 82 INDQ=1,MGRID
          IF(IRAM(IPOSQ0+INDQ).LE.0) THEN
            MAXQ=INDQ-1
            GO TO 83
          END IF
   82   CONTINUE
C       NET-39
        CALL ERROR('NET-39: Too many points of reflector')
C       NET-39: Too many points of reflector:
C               Number of reflector points should be less than the
C               number of network nodes (gridpoints).
   83   CONTINUE
C
C       Labeling the time field in the queue (TT(I)=RAM(ITT0+I).LT.0):
        M=0
        TTMIN=GIANT
        DO 84 I=1,MAXQ
          IADR=INDX(IRAM(IPOSQ0+I))
C         Check for the computational volume and for a free space:
          IF(IADR.GT.0) THEN
            IF(RAM(IP0+IADR).LT.GIANT) THEN
              M=M+1
              IRAM(IPOSQ0+M)=IRAM(IPOSQ0+I)
              TTMIN=AMIN1(RAM(ITT0+IADR),TTMIN)
              RAM(ITT0+IADR)=-RAM(ITT0+IADR)
            END IF
          END IF
   84   CONTINUE
        MINQ=1
        MAXQ=M
C
C       Defining the 1st approximation of the s.p.t. And initially
C       updating the interface:
        DO 88 I=1,NL1*NL2*NL3
          IADR=INDX(I)
          IF(IADR.GT.0) THEN
            IF(RAM(ITT0+IADR).LT.0..AND.RAM(IP0+IADR).LT.GIANT) THEN
              RAM(ITT0+IADR)=-RAM(ITT0+IADR)
              IF(LPRED) THEN
                IRAM(IPRED0+IADR)=0
              ENDIF
            ELSE
              RAM(ITT0+IADR)=GIANT
              IF(LPRED) THEN
                IRAM(IPRED0+IADR)=I
              ENDIF
            ENDIF
          ENDIF
   88   CONTINUE
C
      ENDIF
      IF(MAXQ.LT.MINQ) THEN
C       NET-50
        CALL ERROR('NET-50: No source point')
C       NET-50: No source point:
C               This error should not appear.  Contact the authors.
      ENDIF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE GENER(IREFL)
      INTEGER IREFL
C
C Procedure generating travel-time field and predecessors by performing
C shortest path ray tracing (one iteration).
C
C Input:
C     IREFL.. Number of reflections.
C
C No output.
C
C-----------------------------------------------------------------------
C
C     Common blocks /GRID/ /FS/ /NODE/ are required here:
      INCLUDE 'net.inc'
C     net.inc
      INCLUDE 'netnode.inc'
C     netnode.inc
C
C-----------------------------------------------------------------------
C
      REAL GIANT
      PARAMETER (GIANT=1.E+10)
C
      EXTERNAL INDX,UPDATE,TTTUPD,FREVOL
      INTEGER  INDX
C
      INTEGER NUMNOD,NUMA,NUMINC,NUMOLD,INDQ,I
      REAL TTINDQ
      PARAMETER (NUMINC=1000)
C
C     NUMNOD..Number of network nodes not situated in a free space.
C     NUMA... Counter of finished nodes (nodes moved to set A).
C     NUMINC..Minimum increment of the number of finished nodes between
C             the reports on the screen.  If the increment is smaller,
C             screen output is suppressed.
C     NUMOLD..Counter of finished nodes reported on the screen.
C     INDQ... Loop variable (index in Q).
C     I...    Loop variable (index of a node).
C     TTINDQ..Travel time in the INDQ-th Q element.
C
C.......................................................................
C
C     Check for a free space:
      WRITE(*,'(A)')
     *  '+NET: Eliminating nodes in a free space.                      '
      NUMNOD=L1234
      DO 8 I=1,L1234
        IF(RAM(IP0+I).GE.GIANT) THEN
          RAM(ITT0+I)=-GIANT
          NUMNOD=NUMNOD-1
        ENDIF
    8 CONTINUE
C
C     Counter of finished nodes (nodes moved to set A):
      NUMA=0
      NUMOLD=-NUMINC
C
C     Loop for intervals:
   10 CONTINUE
C       New interval:
        TTMIN=TTMIN+DMIN
C
C       Determination of the first element MINQ of set B in Q:
        DO 21 INDQ=MINQ,MAXQ
          TTINDQ=RAM(ITT0+INDX(IRAM(IPOSQ0+INDQ)))
          IF(0..LT.TTINDQ) THEN
            MINQ=INDQ
            GO TO 22
          END IF
   21   CONTINUE
        MINQ=MAXQ+1
   22   CONTINUE
C
C       Screen output:
        IF(NUMA.GE.NUMOLD+NUMINC.OR.MINQ.GT.MAXQ) THEN
          WRITE(*,'(A,I2,4(A,I7))')
     *          '+',IREFL,'-th reflection',NUMA,' nodes of',NUMNOD,
     *          ' finished,  QMIN=',MINQ,',  QMAX=',MAXQ
          NUMOLD=NUMA
        END IF
C
C       (4) Iteration check: testing the end of time field generation
C       condition to finish the generation of the s.p.t.
        IF(MINQ.GT.MAXQ) THEN
          DO 31 IADR=1,L1234
            RAM(ITT0+IADR)=ABS(RAM(ITT0+IADR))
   31     CONTINUE
          RETURN
        ENDIF
C
        DO 40 INDQ=MINQ,MAXQ
C
C         (2) Selection:
          TTINDQ=RAM(ITT0+INDX(IRAM(IPOSQ0+INDQ)))
          IF(0..LT.TTINDQ.AND.TTINDQ.LT.TTMIN) THEN
C
C           New node 'I' (position IPOS, address IADR):
            IPOS=IRAM(IPOSQ0+INDQ)
            IPOS1=IPOS-1
            IPOS3=IPOS1/(NL1*NL2)
            IPOS2=IPOS1/NL1-IPOS3*NL2
            IPOS1=IPOS1-(IPOS2+IPOS3*NL2)*NL1
            IADR=INDX(IPOS)
            PI=RAM(IP0+IADR)
            TTI=RAM(ITT0+IADR)
C
C           (3) Updating nodes 'J' of the forward star FS(I), which are
C           the neighbours to the node 'I':
            IF(NFSMAX.GE.0) THEN
C             Network ray tracing:
              IF(L4.EQ.0) THEN
                CALL LOOP(UPDATE)
              ELSE
                CALL LOOP(FREVOL)
              END IF
            ELSE
C             Second-order grid travel-time tracing:
              IF(L4.EQ.0) THEN
                CALL TTTUPD()
              ELSE
C               NET-51
                CALL ERROR('NET-51: Fresnel volumes disabled')
C               NET-51: Fresnel volumes disabled:
C                       Second-order grid travel-time tracing cannot be
C                       performed in Fresnel or other volumes specified
C                       by means of the index file.
              END IF
            END IF
C
C           Moving 'I' from set B to set A:
            IADR=INDX(IPOS)
            RAM(ITT0+IADR)=-RAM(ITT0+IADR)
C           Increment of the counter (number of nodes in the set A):
            NUMA=NUMA+1
C
          END IF
   40   CONTINUE
C       End of loop for nodes 'I', which are secondary sources of the
C       time field.
C
      GO TO 10
C     End of loop for intervals.
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE RESFIL(IREFL)
      INTEGER IREFL
C
C In a case of reflections, stores the results of an iteration in an
C unformatted direct-access scratch file (one iteration).
C
C Input:
C     IREFL.. Number of reflections.
C
C No output.
C
C-----------------------------------------------------------------------
C
C     Common blocks /GRID/ /FILES/ are required here:
      INCLUDE 'net.inc'
C     net.inc
C
C-----------------------------------------------------------------------
C
      INTEGER LU0,LU
      PARAMETER(LU0=10)
C
      INTEGER I
C
C     I...    Temporary loop variable.
C
C.......................................................................
C
      WRITE(*,'(A)')
     *  ' NET: Writing travel time field                               '
C
      IF(IREFL.LT.NREFL) THEN
        LU=LU0+IREFL
        OPEN(LU,RECL=8,FORM='UNFORMATTED',ACCESS='DIRECT',
     *                                                 STATUS='SCRATCH')
        IF(LPRED) THEN
          DO 11 I=1,L1234
            WRITE(LU,REC=I) RAM(ITT0+I),IRAM(IPRED0+I)
11        CONTINUE
        ELSE
          DO 12 I=1,L1234
            WRITE(LU,REC=I) RAM(ITT0+I)
12        CONTINUE
        END IF
      END IF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE ERRORS(IREFL)
      INTEGER IREFL
C
C Procedure generating the travel-time errors at all network nodes
C coinciding with gridpoints by means of backward ray tracing.
C Writes 'ERR(I)' (one iteration).
C
C Input:
C     IREFL.. Number of reflections.
C
C No output.
C
C-----------------------------------------------------------------------
C
C     Common blocks /GRID/ /FS/ /FILES/ are required here:
      INCLUDE 'net.inc'
C     net.inc
C
C-----------------------------------------------------------------------
C
      INTEGER LUERR
      PARAMETER (LUERR=5)
C
      EXTERNAL INDX
      INTEGER  INDX
      LOGICAL LTRACE,LRAYS,LEND
      INTEGER IPREDE,IEND1,IEND2,IEND3,IEND,IENDA,IAUX,JCOUNT
      REAL DUMMY
C
C     ERR(I)=RAM(IERR0+I):
      INTEGER  IERR0
      EQUIVALENCE (IERR0,IPOSQ0)
C
C.......................................................................
C
      IF(NFSMAX.LT.0) THEN
C       Second-order grid travel-time tracing: Errors are not generated
        RETURN
      END IF
C
      WRITE(*,'(A)')
     *  '+NET: Generating travel-time errors.                          '
C
      IF(FERR(IREFL).NE.' ') THEN
        LTRACE=.FALSE.
        LRAYS=.FALSE.
        LEND=.FALSE.
C
        DO 10 IEND=1,L1234
          RAM(IERR0+IEND)=-1.
   10   CONTINUE
C
        JCOUNT=0
        DO 23 IEND3=0,NL3-1
          DO 22 IEND2=0,NL2-1
            IEND=NL1*(IEND2+NL2*IEND3)
            DO 21 IEND1=0,NL1-1
              IEND=IEND+1
C             IEND=1+IEND1+NL1*(IEND2+NL2*IEND3)
              IENDA=INDX(IEND)
              IF(IENDA.GT.0) THEN
                IF(RAM(IERR0+IENDA).LT.0.) THEN
                  IPREDE=IRAM(IPRED0+IENDA)
                  IF(IPREDE.NE.IEND) THEN
                    CALL ONERAY(LTRACE,LRAYS,LEND,IAUX,IREFL,IPREDE,
     *                          IEND1,IEND2,IEND3,IENDA,RAM(ITT0+IENDA),
     *                          RAM(IERR0+1),DUMMY,IAUX)
                  END IF
                END IF
                JCOUNT=JCOUNT+1
                IF(JCOUNT/1000*1000.EQ.JCOUNT) THEN
                  WRITE(*,'(A,I7,A,I7)')
     *              '+NET: Generating travel-time errors:',
     *              JCOUNT,' of',L1234
                END IF
              END IF
   21       CONTINUE
   22     CONTINUE
   23   CONTINUE
C
C.......................................................................
C
C       Writing the travel-time error estimates:
C
        WRITE(*,'(A)')
     *  '+NET: Writing travel-time errors.                             '
C
        CALL WARRAY(LUERR,FERR(IREFL),'FORMATTED',
     *              .TRUE.,-1.,.FALSE.,0.,L1234,RAM(IERR0+1))
C
      END IF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE RECVRS()
C
C Procedure updating arrival times at the receivers.
C
C No input, no output.
C
C-----------------------------------------------------------------------
C
C     Common blocks /GRID/ /FS/ /NAMESR/ /POINTS/ /FILES/ /NODE/ /SRCC/
C     are required here:
      INCLUDE 'net.inc'
C     net.inc
      INCLUDE 'netnode.inc'
C     netnode.inc
C
C-----------------------------------------------------------------------
C
      EXTERNAL SRCREC
C
      REAL GIANT
      PARAMETER (GIANT=1.E+10)
C
      INTEGER IR,ISRC,NFSREC,NFSSRC,NFSMAK,I1,I2,I3,I4,I5
      REAL DIM,PS,A1,A2,A3,DIST2,TTSR
C
C     IR...   Index of the receiver.
C     ISRC... Index of the source.
C     NFSREC..Size of the f.s. at the receiver.
C     NFSSRC..Size of the f.s. at the source.
C     NFSMAK..Size of recently generated full forward star at source
C             points.
C     DIM...  2 or 3 for 2-D or 3-D calculation, respectively.
C     PS...   Slowness at the source.
C     A1,A2,A3,I1,I2,I3,I4,I5... Other source parameters.
C     DIST2.. Square of the source-reiver distance.
C     TTSR... Arrival time at the receiver from the source under
C             consideration.
C
C.......................................................................
C
      IF(NFSMAX.LT.0) THEN
C       Second-order grid travel-time tracing:
        IF(NREC.GT.0) THEN
          WRITE(*,'(A)')
     *'+NET: Receivers are not supported by the 2-D second-order method'
          WRITE(*,'(A)') ' '
        END IF
        RETURN
      END IF
C
C     Computing approximate arrival time and rays at the receiver nodes:
      WRITE(*,'(A)')
     *  '+NET: Updating receivers...                                   '
C
C     Extent of a dense forward star at receiver points:
      IF(LNFS) THEN
        NFSMAK=0
      ELSE
        NFSMAK=NFSMAX
        CALL MAKEFS(NFSMAK)
      END IF
C
C     Loop over receivers:
      DO 90 IR=1,NREC
C
C       Initially, receiver is considered not connected:
        IPREDR(IR)=0
C
C       Parameters of the receiver point:
        CALL SLOW(X1R(IR),X2R(IR),X3R(IR),DPOS1,DPOS2,DPOS3,
     *                                   IPOS1,IPOS2,IPOS3,IPOS,IADR,PI)
        TTI=GIANT
        ISRCI=IR
C
C.......................................................................
C
C       For no reflections, testing the source-receiver position:
C
        IF(NREFL.EQ.0) THEN
C         Loop over source points:
          DO 39 ISRC=1,NSRC
C
C           Square of the source-reiver distance:
            DIST2=0.
            IF(LN1) THEN
              DIST2=DIST2+((X1S(ISRC)-X1R(IR))**2
     *                    +(X2S(ISRC)-X2R(IR))**2)/ASX2**2
            END IF
            IF(NL1.GT.1) THEN
              DIST2=DIST2+((X1S(ISRC)-X1R(IR))/ASX1)**2
            END IF
            IF(NL2.GT.1) THEN
              DIST2=DIST2+((X2S(ISRC)-X2R(IR))/ASX2)**2
            END IF
            IF(NL3.GT.1) THEN
              DIST2=DIST2+((X3S(ISRC)-X3R(IR))/ASX3)**2
            END IF
C
C           Size of the f.s. at the receiver point:
            IF(LNFS) THEN
              NFSREC=IRAM(NFS0+IADR)
            ELSE
              NFSREC=NFSMAX
            END IF
C
C           Dimension of the model and forward stars (2 or 3):
            IF(NL1.EQ.1.OR.NL2.EQ.1.OR.NL3.EQ.1) THEN
              DIM=2.
            ELSE
              DIM=3.
            END IF
C
            IF(DIST2.LE.FLOAT(NFSREC)**2+DIM-1.) THEN
C             Source is situated within the receiver f.s. ellipsoid
              CALL SLOW(X1S(ISRC),X2S(ISRC),X3S(ISRC),A1,A2,A3,
     *                                                I1,I2,I3,I4,I5,PS)
C
C             Check for crossing an interface:
              IF(LICB) THEN
                IF(IRAM(ICB0+IADR).NE.IRAM(ICB0+I5)) THEN
                  IF(DIST2.GT.DIM) GO TO 38
                END IF
              END IF
C
C             Size of the f.s. at the source point:
              IF(LNFS) THEN
                NFSSRC=IRAM(NFS0+I5)
              ELSE
                NFSSRC=NFSMAX
              END IF
C
C             Updating:
              IF(DIST2.LE.FLOAT(NFSSRC)**2+DIM-1.) THEN
C               Travel time from source to the receiver
                TTSR=SQRT((X1R(IR)-X1S(ISRC))**2+
     *                    (X2R(IR)-X2S(ISRC))**2+
     *                    (X3R(IR)-X3S(ISRC))**2)*0.5*(PI+PS)+TTS(ISRC)
C               Updating
                IF(IPREDR(IR).LT.0) THEN
                  IF(TTR(IR).GT.TTSR) THEN
                    IPREDR(IR)=-ISRC
                    TTR(IR)=TTSR
                  END IF
                ELSE
                  IPREDR(IR)=-ISRC
                  TTR(IR)=TTSR
                END IF
              END IF
C
   38         CONTINUE
            END IF
   39     CONTINUE
        END IF
C
C.......................................................................
C
C       Connecting the receiver still not connected:
C
        IF(IPREDR(IR).EQ.0) THEN
C         Adjusting extent of the dense forward star at the receiver
C         point:
          IF(LNFS) THEN
            IF(IRAM(NFS0+IADR).NE.NFSMAK) THEN
              NFSMAK=IRAM(NFS0+IADR)
              CALL MAKEFS(NFSMAK)
            END IF
          END IF
C
C         Updating:
          CALL LOOP(SRCREC)
          TTR(IR)=TTI
        ENDIF
C
   90 CONTINUE
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE TRACER()
C
C Procedure performing backward ray tracing from the receivers,
C including the estimation of arrival-time errors.
C
C No input, no output.
C
C-----------------------------------------------------------------------
C
C     Common blocks /GRID/ /FS/ /NAMESR/ /POINTS/ /FILES/ are required:
      INCLUDE 'net.inc'
C     net.inc
C
C-----------------------------------------------------------------------
C
      INTEGER LURAYS,LUEND
      PARAMETER (LURAYS=4)
      PARAMETER (LUEND=5)
C
      LOGICAL LTRACE,LRAYS,LEND
      INTEGER IPREDE,IEND1,IEND2,IEND3,IAUX,IENDA,IREC,ISRC
      REAL DUMMY1,DUMMY2,DUMMY3,PEND,TERR,DUMMY
C
C     IREC... Loop variable - index of a receiver.
C
C.......................................................................
C
      IF(NFSMAX.LT.0) THEN
C       Second-order grid travel-time tracing: No posterior ray tracing
        RETURN
      END IF
C
      LTRACE=.TRUE.
      WRITE(*,'(A)')
     *  '+NET: Backward ray tracing...                                 '
      IF(FRAYS.EQ.' ') THEN
        LRAYS=.FALSE.
      ELSE
        LRAYS=.TRUE.
      END IF
      IF(FEND.EQ.' ') THEN
        LEND=.FALSE.
      ELSE
        LEND=.TRUE.
      END IF
C
      IF(LRAYS) THEN
        OPEN(LURAYS,FILE=FRAYS)
        WRITE(LURAYS,'(A)') '/'
      END IF
      IF(LEND) THEN
        OPEN(LUEND,FILE=FEND)
        WRITE(LUEND,'(A)') '/'
      END IF
      DO 80 IREC=1,NREC
        IPREDE=IPREDR(IREC)
        IF(IPREDE.NE.0) THEN
          IF(LRAYS) THEN
            IF(NAMES.EQ.' ') THEN
              IF(NAMER(IREC).EQ.' ') THEN
                WRITE(LURAYS,'(2A,I5,5A)')
     *          '''','RAY',IREC,                                  ''' /'
              ELSE
                WRITE(LURAYS,'(2A,I5,5A)')
     *          '''','RAY',IREC,'      ',NAMES,' TO ',NAMER(IREC),''' /'
              END IF
            ELSE
              IF(NAMER(IREC).EQ.' ') THEN
                WRITE(LURAYS,'(2A,I5,5A)')
     *          '''','RAY',IREC,' FROM ',NAMES,                   ''' /'
              ELSE
                WRITE(LURAYS,'(2A,I5,5A)')
     *          '''','RAY',IREC,' FROM ',NAMES,' TO ',NAMER(IREC),''' /'
              END IF
            END IF
            WRITE(LURAYS,'(4F12.6,A)')
     *                      X1R(IREC),X2R(IREC),X3R(IREC),TTR(IREC),' /'
          END IF
          IF(LEND) THEN
            CALL SLOW(X1R(IREC),X2R(IREC),X3R(IREC),
     *           DUMMY1,DUMMY2,DUMMY3,IEND1,IEND2,IEND3,IAUX,IENDA,PEND)
          END IF
          CALL ONERAY(LTRACE,LRAYS,LEND,LURAYS,NREFL,
     *         IPREDE,IEND1,IEND2,IEND3,IENDA,TTR(IREC),DUMMY,TERR,ISRC)
          IF(LRAYS) THEN
            WRITE(LURAYS,'(4F12.6,A)')
     *                      X1S(ISRC),X2S(ISRC),X3S(ISRC),TTS(ISRC),' /'
            WRITE(LURAYS,'(''/'')')
          END IF
          IF(LEND) THEN
            WRITE(LUEND,'(3A,5F12.6,A)') '''',NAMER(IREC),'''',
     *                 X1R(IREC),X2R(IREC),X3R(IREC),TTR(IREC),TERR,' /'
          END IF
        END IF
   80 CONTINUE
      IF(LRAYS) THEN
        WRITE(LURAYS,'(''/'')')
        CLOSE(LURAYS)
      END IF
      IF(LEND) THEN
        WRITE(LUEND,'(''/'')')
        CLOSE(LUEND)
      END IF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE OUTPUT()
C
C Subroutine to rewrite arrival-time fields and predecessors from memory
C or scratch files to formatted output files 'TT(I)' and 'IPRED(I)'.
C
C No input, no output.
C
C-----------------------------------------------------------------------
C
C     Common blocks /GRID/ /FS/ /FILES/ are required here:
      INCLUDE 'net.inc'
C     net.inc
C
C-----------------------------------------------------------------------
C
      INTEGER LU0,LU1,LU
      PARAMETER(LU0=10)
      PARAMETER(LU1=1)
      REAL GIANT
      PARAMETER (GIANT=1.0E+10)
C
      INTEGER NF,I
C     NF,I... Loop variables.
C
C.......................................................................
C
      WRITE(*,'(A)')
     *  '+NET: Generating output files...                              '
C
      DO 30 NF=NREFL,0,-1
C
C       Reading scratch file:
        IF(NF.NE.NREFL) THEN
          LU=LU0+NF
          IF(FTT(NF).NE.' '.OR.FPRED(NF).NE.' ') THEN
            IF(LPRED) THEN
              DO 11 I=1,L1234
                READ(LU,REC=I) RAM(ITT0+I),IRAM(IPRED0+I)
   11         CONTINUE
            ELSE
              DO 12 I=1,L1234
                READ(LU,REC=I) RAM(ITT0+I)
   12         CONTINUE
            END IF
          END IF
C         Closing the scratch file
C         *********
          CLOSE(LU)
C         *********
        END IF
C
C       Writing arrival times:
        IF(FTT(NF).NE.' ') THEN
          CALL WARRAY(LU1,FTT(NF),'FORMATTED',
     *                .FALSE.,0.,.TRUE.,GIANT,L1234,RAM(ITT0+1))
        END IF
C
C       Writing predecessors:
        IF(NFSMAX.LT.0) THEN
C         Second-order grid travel-time tracing
C         (writing slowness vector):
          IF(N1.GT.1.AND.FERR(NF).NE.' ') THEN
            CALL WARRAY(LU1,FERR(NF),'FORMATTED',
     *                  .FALSE.,0.,.FALSE.,GIANT,L1234,RAM(IP0+1))
          END IF
          IF(N2.GT.1.AND.FPRED(NF).NE.' ') THEN
            CALL WARRAY(LU1,FPRED(NF),'FORMATTED',
     *                  .FALSE.,0.,.FALSE.,GIANT,L1234,RAM(IP20+1))
          END IF
C         IF(N3.GT.1.AND.FNFS(NF).NE.' ') THEN
C           CALL WARRAY(LU1,FNFS(NF),'FORMATTED',
C    *                  .FALSE.,0.,.FALSE.,GIANT,L1234,RAM(IP30+1))
C         END IF
        ELSE IF(FPRED(NF).NE.' ') THEN
          CALL WARRAI(LU1,FPRED(NF),'FORMATTED',
     *                .FALSE.,0,.FALSE.,0,L1234,IRAM(IPRED0+1))
        END IF
   30 CONTINUE
C
      WRITE(*,'(A)')
     *  '+NET: Done.                                                   '
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE LOOP(UPDATE)
      EXTERNAL UPDATE
C
C This subroutine performs the loop over the nodes 'J' of the forward
C star FS(I).  A clear but exact arrangement of the loop is expressed
C inside asterisks.  The code following the asterisk rectangular is
C equivalent, but much longer and somewhat faster.
C Auxiliary subroutine to SOURCE, GENER, and RECVRS.
C LOOP(FREVOL) or LOOP(UPDATE): employed by GENER,
C LOOP(SRCREC): employed by SOURCE and RECVRS.
C
C-----------------------------------------------------------------------
C
C     Common blocks /GRID/ /FS/ /NODE/ are required here:
      INCLUDE 'net.inc'
C     net.inc
      INCLUDE 'netnode.inc'
C     netnode.inc
C
C-----------------------------------------------------------------------
C
      INTEGER I1MIN,I1MAX,I2MIN,I2MAX,I3MIN,I3MAX,I1,I2,I3
      INTEGER NFSOLD,NLEX1,NLEX2,NLEX3,IFSMIN,IFSMAX
      SAVE    NFSOLD,NLEX1,NLEX2,NLEX3,IFSMIN,IFSMAX
C     Still no forward star used:
      DATA    NFSOLD/0/
C
C     I1MIN,I1MAX,I2MIN,I2MAX,I3MIN,I3MAX...
C     I1,I2,I3...
C     NFSOLD...
C     NLEX1,NLEX2,NLEX3...
C     IFSMIN,IFSMAX...
C
C.......................................................................
C
C     Additional check of input data
      IF(NFSMAX.LT.0) THEN
C       Second-order grid travel-time tracing -- inaccessible branch:
C       NET-55
        CALL ERROR('NET-55: LOOP should not be called')
C       NET-55: LOOP should not be called:
C               This error should not appear.  Contact the authors.
      END IF
C
C     Index of the block (if applicable):
      IF(LICB) THEN
        ICBI=IRAM(ICB0+IADR)
      END IF
C
C     Adjusting extent of the forward star:
      IF(NFSOLD.EQ.0) THEN
        NFSOLD=NFSMAX
        NLEX1=MIN0(NFSOLD,NL1-1)
        NLEX2=MIN0(NFSOLD,NL2-1)
        NLEX3=MIN0(NFSOLD,NL3-1)
      END IF
      IF(LNFS) THEN
        IF(IRAM(NFS0+IADR).NE.NFSOLD) THEN
          NFSOLD=IRAM(NFS0+IADR)
          NLEX1=MIN0(NFSOLD,NL1-1)
          NLEX2=MIN0(NFSOLD,NL2-1)
          NLEX3=MIN0(NFSOLD,NL3-1)
        END IF
      END IF
C
C     Location of the forward star in the memory:
      IFSMIN=KFS0(NFSOLD-1)+1
      IFSMAX=KFS0(NFSOLD)
C
C     Extent of the intersection of the forward star with the grid:
      I1MIN=     -IPOS1
      I1MAX=NL1-1-IPOS1
      I2MIN=     -IPOS2
      I2MAX=NL2-1-IPOS2
      I3MIN=     -IPOS3
      I3MAX=NL3-1-IPOS3
C
************************************************************************
*     DO 163 IFS=IFSMIN,IFSMAX                                         *
*       I1=KFS1(IFS)                                                   *
*       I2=KFS2(IFS)                                                   *
*       I3=KFS3(IFS)                                                   *
*       IF(I1MIN.LE.I1.AND.I1.LE.I1MAX.AND.                            *
*    *     I2MIN.LE.I2.AND.I2.LE.I2MAX.AND.                            *
*    *     I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN                           *
*         CALL UPDATE()                                                *
*       END IF                                                         *
* 163 CONTINUE                                                         *
************************************************************************
C
C     This is an optimized version of the above loop:
      IF(I3MIN.LE.-NLEX3) THEN
        IF(NLEX3.LE.I3MAX) THEN
          IF(I2MIN.LE.-NLEX2) THEN
            IF(NLEX2.LE.I2MAX) THEN
              IF(I1MIN.LE.-NLEX1) THEN
                IF(NLEX1.LE.I1MAX) THEN
                  DO 100 IFS=IFSMIN,IFSMAX
                      CALL UPDATE()
  100             CONTINUE
                ELSE
                  DO 101 IFS=IFSMIN,IFSMAX
                    IF(KFS1(IFS).LE.I1MAX) THEN
                      CALL UPDATE()
                    END IF
  101             CONTINUE
                END IF
              ELSE
                IF(NLEX1.LE.I1MAX) THEN
                  DO 102 IFS=IFSMIN,IFSMAX
                    IF(I1MIN.LE.KFS1(IFS)) THEN
                      CALL UPDATE()
                    END IF
  102             CONTINUE
                ELSE
                  DO 103 IFS=IFSMIN,IFSMAX
                    I1=KFS1(IFS)
                    IF(I1MIN.LE.I1.AND.I1.LE.I1MAX) THEN
                      CALL UPDATE()
                    END IF
  103             CONTINUE
                END IF
              END IF
            ELSE
              IF(I1MIN.LE.-NLEX1) THEN
                IF(NLEX1.LE.I1MAX) THEN
                  DO 104 IFS=IFSMIN,IFSMAX
                    IF(KFS2(IFS).LE.I2MAX) THEN
                      CALL UPDATE()
                    END IF
  104             CONTINUE
                ELSE
                  DO 105 IFS=IFSMIN,IFSMAX
                    IF(KFS1(IFS).LE.I1MAX.AND.
     *                 KFS2(IFS).LE.I2MAX) THEN
                      CALL UPDATE()
                    END IF
  105             CONTINUE
                END IF
              ELSE
                IF(NLEX1.LE.I1MAX) THEN
                  DO 106 IFS=IFSMIN,IFSMAX
                    IF(I1MIN.LE.KFS1(IFS).AND.
     *                          KFS2(IFS).LE.I2MAX) THEN
                      CALL UPDATE()
                    END IF
  106             CONTINUE
                ELSE
                  DO 107 IFS=IFSMIN,IFSMAX
                    I1=KFS1(IFS)
                    IF(I1MIN.LE.I1.AND.I1.LE.I1MAX.AND.
     *                          KFS2(IFS).LE.I2MAX) THEN
                      CALL UPDATE()
                    END IF
  107             CONTINUE
                END IF
              END IF
            END IF
          ELSE
            IF(NLEX2.LE.I2MAX) THEN
              IF(I1MIN.LE.-NLEX1) THEN
                IF(NLEX1.LE.I1MAX) THEN
                  DO 108 IFS=IFSMIN,IFSMAX
                    IF(I2MIN.LE.KFS2(IFS)) THEN
                      CALL UPDATE()
                    END IF
  108             CONTINUE
                ELSE
                  DO 109 IFS=IFSMIN,IFSMAX
                    IF(         KFS1(IFS).LE.I1MAX.AND.
     *                 I2MIN.LE.KFS2(IFS)) THEN
                      CALL UPDATE()
                    END IF
  109             CONTINUE
                END IF
              ELSE
                IF(NLEX1.LE.I1MAX) THEN
                  DO 110 IFS=IFSMIN,IFSMAX
                    IF(I1MIN.LE.KFS1(IFS).AND.
     *                 I2MIN.LE.KFS2(IFS)) THEN
                      CALL UPDATE()
                    END IF
  110             CONTINUE
                ELSE
                  DO 111 IFS=IFSMIN,IFSMAX
                    I1=KFS1(IFS)
                    IF(I1MIN.LE.I1.AND.I1.LE.I1MAX.AND.
     *                 I2MIN.LE.KFS2(IFS)) THEN
                      CALL UPDATE()
                    END IF
  111             CONTINUE
                END IF
              END IF
            ELSE
              IF(I1MIN.LE.-NLEX1) THEN
                IF(NLEX1.LE.I1MAX) THEN
                  DO 112 IFS=IFSMIN,IFSMAX
                    I2=KFS2(IFS)
                    IF(I2MIN.LE.I2.AND.I2.LE.I2MAX) THEN
                      CALL UPDATE()
                    END IF
  112             CONTINUE
                ELSE
                  DO 113 IFS=IFSMIN,IFSMAX
                    I2=KFS2(IFS)
                    IF(         KFS1(IFS).LE.I1MAX.AND.
     *                 I2MIN.LE.I2.AND.I2.LE.I2MAX) THEN
                      CALL UPDATE()
                    END IF
  113             CONTINUE
                END IF
              ELSE
                IF(NLEX1.LE.I1MAX) THEN
                  DO 114 IFS=IFSMIN,IFSMAX
                    I2=KFS2(IFS)
                    IF(I1MIN.LE.KFS1(IFS).AND.
     *                 I2MIN.LE.I2.AND.I2.LE.I2MAX) THEN
                      CALL UPDATE()
                    END IF
  114             CONTINUE
                ELSE
                  DO 115 IFS=IFSMIN,IFSMAX
                    I1=KFS1(IFS)
                    I2=KFS2(IFS)
                    IF(I1MIN.LE.I1.AND.I1.LE.I1MAX.AND.
     *                 I2MIN.LE.I2.AND.I2.LE.I2MAX) THEN
                      CALL UPDATE()
                    END IF
  115             CONTINUE
                END IF
              END IF
            END IF
          END IF
        ELSE
          IF(I2MIN.LE.-NLEX2) THEN
            IF(NLEX2.LE.I2MAX) THEN
              IF(I1MIN.LE.-NLEX1) THEN
                IF(NLEX1.LE.I1MAX) THEN
                  DO 116 IFS=IFSMIN,IFSMAX
                    IF(KFS3(IFS).LE.I3MAX) THEN
                      CALL UPDATE()
                    END IF
  116             CONTINUE
                ELSE
                  DO 117 IFS=IFSMIN,IFSMAX
                    IF(KFS1(IFS).LE.I1MAX.AND.
     *                 KFS3(IFS).LE.I3MAX) THEN
                      CALL UPDATE()
                    END IF
  117             CONTINUE
                END IF
              ELSE
                IF(NLEX1.LE.I1MAX) THEN
                  DO 118 IFS=IFSMIN,IFSMAX
                    IF(I1MIN.LE.KFS1(IFS).AND.
     *                          KFS3(IFS).LE.I3MAX) THEN
                      CALL UPDATE()
                    END IF
  118             CONTINUE
                ELSE
                  DO 119 IFS=IFSMIN,IFSMAX
                    I1=KFS1(IFS)
                    IF(I1MIN.LE.I1.AND.I1.LE.I1MAX.AND.
     *                          KFS3(IFS).LE.I3MAX) THEN
                      CALL UPDATE()
                    END IF
  119             CONTINUE
                END IF
              END IF
            ELSE
              IF(I1MIN.LE.-NLEX1) THEN
                IF(NLEX1.LE.I1MAX) THEN
                  DO 120 IFS=IFSMIN,IFSMAX
                    IF(KFS2(IFS).LE.I2MAX.AND.
     *                 KFS3(IFS).LE.I3MAX) THEN
                      CALL UPDATE()
                    END IF
  120             CONTINUE
                ELSE
                  DO 121 IFS=IFSMIN,IFSMAX
                    IF(KFS1(IFS).LE.I1MAX.AND.
     *                 KFS2(IFS).LE.I2MAX.AND.
     *                 KFS3(IFS).LE.I3MAX) THEN
                      CALL UPDATE()
                    END IF
  121             CONTINUE
                END IF
              ELSE
                IF(NLEX1.LE.I1MAX) THEN
                  DO 122 IFS=IFSMIN,IFSMAX
                    IF(I1MIN.LE.KFS1(IFS).AND.
     *                          KFS2(IFS).LE.I2MAX.AND.
     *                          KFS3(IFS).LE.I3MAX) THEN
                      CALL UPDATE()
                    END IF
  122             CONTINUE
                ELSE
                  DO 123 IFS=IFSMIN,IFSMAX
                    I1=KFS1(IFS)
                    IF(I1MIN.LE.I1.AND.I1.LE.I1MAX.AND.
     *                          KFS2(IFS).LE.I2MAX.AND.
     *                          KFS3(IFS).LE.I3MAX) THEN
                      CALL UPDATE()
                    END IF
  123             CONTINUE
                END IF
              END IF
            END IF
          ELSE
            IF(NLEX2.LE.I2MAX) THEN
              IF(I1MIN.LE.-NLEX1) THEN
                IF(NLEX1.LE.I1MAX) THEN
                  DO 124 IFS=IFSMIN,IFSMAX
                    IF(I2MIN.LE.KFS2(IFS).AND.
     *                          KFS3(IFS).LE.I3MAX) THEN
                      CALL UPDATE()
                    END IF
  124             CONTINUE
                ELSE
                  DO 125 IFS=IFSMIN,IFSMAX
                    IF(         KFS1(IFS).LE.I1MAX.AND.
     *                 I2MIN.LE.KFS2(IFS).AND.
     *                          KFS3(IFS).LE.I3MAX) THEN
                      CALL UPDATE()
                    END IF
  125             CONTINUE
                END IF
              ELSE
                IF(NLEX1.LE.I1MAX) THEN
                  DO 126 IFS=IFSMIN,IFSMAX
                    IF(I1MIN.LE.KFS1(IFS).AND.
     *                 I2MIN.LE.KFS2(IFS).AND.
     *                          KFS3(IFS).LE.I3MAX) THEN
                      CALL UPDATE()
                    END IF
  126             CONTINUE
                ELSE
                  DO 127 IFS=IFSMIN,IFSMAX
                    I1=KFS1(IFS)
                    IF(I1MIN.LE.I1.AND.I1.LE.I1MAX.AND.
     *                 I2MIN.LE.KFS2(IFS).AND.
     *                          KFS3(IFS).LE.I3MAX) THEN
                      CALL UPDATE()
                    END IF
  127             CONTINUE
                END IF
              END IF
            ELSE
              IF(I1MIN.LE.-NLEX1) THEN
                IF(NLEX1.LE.I1MAX) THEN
                  DO 128 IFS=IFSMIN,IFSMAX
                    I2=KFS2(IFS)
                    IF(I2MIN.LE.I2.AND.I2.LE.I2MAX.AND.
     *                          KFS3(IFS).LE.I3MAX) THEN
                      CALL UPDATE()
                    END IF
  128             CONTINUE
                ELSE
                  DO 129 IFS=IFSMIN,IFSMAX
                    I2=KFS2(IFS)
                    IF(         KFS1(IFS).LE.I1MAX.AND.
     *                 I2MIN.LE.I2.AND.I2.LE.I2MAX.AND.
     *                          KFS3(IFS).LE.I3MAX) THEN
                      CALL UPDATE()
                    END IF
  129             CONTINUE
                END IF
              ELSE
                IF(NLEX1.LE.I1MAX) THEN
                  DO 130 IFS=IFSMIN,IFSMAX
                    I2=KFS2(IFS)
                    IF(I1MIN.LE.KFS1(IFS).AND.
     *                 I2MIN.LE.I2.AND.I2.LE.I2MAX.AND.
     *                          KFS3(IFS).LE.I3MAX) THEN
                      CALL UPDATE()
                    END IF
  130             CONTINUE
                ELSE
                  DO 131 IFS=IFSMIN,IFSMAX
                    I1=KFS1(IFS)
                    I2=KFS2(IFS)
                    IF(I1MIN.LE.I1.AND.I1.LE.I1MAX.AND.
     *                 I2MIN.LE.I2.AND.I2.LE.I2MAX.AND.
     *                          KFS3(IFS).LE.I3MAX) THEN
                      CALL UPDATE()
                    END IF
  131             CONTINUE
                END IF
              END IF
            END IF
          END IF
        END IF
      ELSE
        IF(NLEX3.LE.I3MAX) THEN
          IF(I2MIN.LE.-NLEX2) THEN
            IF(NLEX2.LE.I2MAX) THEN
              IF(I1MIN.LE.-NLEX1) THEN
                IF(NLEX1.LE.I1MAX) THEN
                  DO 132 IFS=IFSMIN,IFSMAX
                    IF(I3MIN.LE.KFS3(IFS)) THEN
                      CALL UPDATE()
                    END IF
  132             CONTINUE
                ELSE
                  DO 133 IFS=IFSMIN,IFSMAX
                    IF(         KFS1(IFS).LE.I1MAX.AND.
     *                 I3MIN.LE.KFS3(IFS)) THEN
                      CALL UPDATE()
                    END IF
  133             CONTINUE
                END IF
              ELSE
                IF(NLEX1.LE.I1MAX) THEN
                  DO 134 IFS=IFSMIN,IFSMAX
                    IF(I1MIN.LE.KFS1(IFS).AND.
     *                 I3MIN.LE.KFS3(IFS)) THEN
                      CALL UPDATE()
                    END IF
  134             CONTINUE
                ELSE
                  DO 135 IFS=IFSMIN,IFSMAX
                    I1=KFS1(IFS)
                    IF(I1MIN.LE.I1.AND.I1.LE.I1MAX.AND.
     *                 I3MIN.LE.KFS3(IFS)) THEN
                      CALL UPDATE()
                    END IF
  135             CONTINUE
                END IF
              END IF
            ELSE
              IF(I1MIN.LE.-NLEX1) THEN
                IF(NLEX1.LE.I1MAX) THEN
                  DO 136 IFS=IFSMIN,IFSMAX
                    IF(         KFS2(IFS).LE.I2MAX.AND.
     *                 I3MIN.LE.KFS3(IFS)) THEN
                      CALL UPDATE()
                    END IF
  136             CONTINUE
                ELSE
                  DO 137 IFS=IFSMIN,IFSMAX
                    IF(         KFS1(IFS).LE.I1MAX.AND.
     *                          KFS2(IFS).LE.I2MAX.AND.
     *                 I3MIN.LE.KFS3(IFS)) THEN
                      CALL UPDATE()
                    END IF
  137             CONTINUE
                END IF
              ELSE
                IF(NLEX1.LE.I1MAX) THEN
                  DO 138 IFS=IFSMIN,IFSMAX
                    IF(I1MIN.LE.KFS1(IFS).AND.
     *                          KFS2(IFS).LE.I2MAX.AND.
     *                 I3MIN.LE.KFS3(IFS)) THEN
                      CALL UPDATE()
                    END IF
  138             CONTINUE
                ELSE
                  DO 139 IFS=IFSMIN,IFSMAX
                    I1=KFS1(IFS)
                    IF(I1MIN.LE.I1.AND.I1.LE.I1MAX.AND.
     *                          KFS2(IFS).LE.I2MAX.AND.
     *                 I3MIN.LE.KFS3(IFS)) THEN
                      CALL UPDATE()
                    END IF
  139             CONTINUE
                END IF
              END IF
            END IF
          ELSE
            IF(NLEX2.LE.I2MAX) THEN
              IF(I1MIN.LE.-NLEX1) THEN
                IF(NLEX1.LE.I1MAX) THEN
                  DO 140 IFS=IFSMIN,IFSMAX
                    IF(I2MIN.LE.KFS2(IFS).AND.
     *                 I3MIN.LE.KFS3(IFS)) THEN
                      CALL UPDATE()
                    END IF
  140             CONTINUE
                ELSE
                  DO 141 IFS=IFSMIN,IFSMAX
                    IF(         KFS1(IFS).LE.I1MAX.AND.
     *                 I2MIN.LE.KFS2(IFS).AND.
     *                 I3MIN.LE.KFS3(IFS)) THEN
                      CALL UPDATE()
                    END IF
  141             CONTINUE
                END IF
              ELSE
                IF(NLEX1.LE.I1MAX) THEN
                  DO 142 IFS=IFSMIN,IFSMAX
                    IF(I1MIN.LE.KFS1(IFS).AND.
     *                 I2MIN.LE.KFS2(IFS).AND.
     *                 I3MIN.LE.KFS3(IFS)) THEN
                      CALL UPDATE()
                    END IF
  142             CONTINUE
                ELSE
                  DO 143 IFS=IFSMIN,IFSMAX
                    I1=KFS1(IFS)
                    IF(I1MIN.LE.I1.AND.I1.LE.I1MAX.AND.
     *                 I2MIN.LE.KFS2(IFS).AND.
     *                 I3MIN.LE.KFS3(IFS)) THEN
                      CALL UPDATE()
                    END IF
  143             CONTINUE
                END IF
              END IF
            ELSE
              IF(I1MIN.LE.-NLEX1) THEN
                IF(NLEX1.LE.I1MAX) THEN
                  DO 144 IFS=IFSMIN,IFSMAX
                    I2=KFS2(IFS)
                    IF(I2MIN.LE.I2.AND.I2.LE.I2MAX.AND.
     *                 I3MIN.LE.KFS3(IFS)) THEN
                      CALL UPDATE()
                    END IF
  144             CONTINUE
                ELSE
                  DO 145 IFS=IFSMIN,IFSMAX
                    I2=KFS2(IFS)
                    IF(         KFS1(IFS).LE.I1MAX.AND.
     *                 I2MIN.LE.I2.AND.I2.LE.I2MAX.AND.
     *                 I3MIN.LE.KFS3(IFS)) THEN
                      CALL UPDATE()
                    END IF
  145             CONTINUE
                END IF
              ELSE
                IF(NLEX1.LE.I1MAX) THEN
                  DO 146 IFS=IFSMIN,IFSMAX
                    I2=KFS2(IFS)
                    IF(I1MIN.LE.KFS1(IFS).AND.
     *                 I2MIN.LE.I2.AND.I2.LE.I2MAX.AND.
     *                 I3MIN.LE.KFS3(IFS)) THEN
                      CALL UPDATE()
                    END IF
  146             CONTINUE
                ELSE
                  DO 147 IFS=IFSMIN,IFSMAX
                    I1=KFS1(IFS)
                    I2=KFS2(IFS)
                    IF(I1MIN.LE.I1.AND.I1.LE.I1MAX.AND.
     *                 I2MIN.LE.I2.AND.I2.LE.I2MAX.AND.
     *                 I3MIN.LE.KFS3(IFS)) THEN
                      CALL UPDATE()
                    END IF
  147             CONTINUE
                END IF
              END IF
            END IF
          END IF
        ELSE
          IF(I2MIN.LE.-NLEX2) THEN
            IF(NLEX2.LE.I2MAX) THEN
              IF(I1MIN.LE.-NLEX1) THEN
                IF(NLEX1.LE.I1MAX) THEN
                  DO 148 IFS=IFSMIN,IFSMAX
                    I3=KFS3(IFS)
                    IF(I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN
                      CALL UPDATE()
                    END IF
  148             CONTINUE
                ELSE
                  DO 149 IFS=IFSMIN,IFSMAX
                    I3=KFS3(IFS)
                    IF(         KFS1(IFS).LE.I1MAX.AND.
     *                 I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN
                      CALL UPDATE()
                    END IF
  149             CONTINUE
                END IF
              ELSE
                IF(NLEX1.LE.I1MAX) THEN
                  DO 150 IFS=IFSMIN,IFSMAX
                    I3=KFS3(IFS)
                    IF(I1MIN.LE.KFS1(IFS).AND.
     *                 I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN
                      CALL UPDATE()
                    END IF
  150             CONTINUE
                ELSE
                  DO 151 IFS=IFSMIN,IFSMAX
                    I1=KFS1(IFS)
                    I3=KFS3(IFS)
                    IF(I1MIN.LE.I1.AND.I1.LE.I1MAX.AND.
     *                 I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN
                      CALL UPDATE()
                    END IF
  151             CONTINUE
                END IF
              END IF
            ELSE
              IF(I1MIN.LE.-NLEX1) THEN
                IF(NLEX1.LE.I1MAX) THEN
                  DO 152 IFS=IFSMIN,IFSMAX
                    I3=KFS3(IFS)
                    IF(         KFS2(IFS).LE.I2MAX.AND.
     *                 I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN
                      CALL UPDATE()
                    END IF
  152             CONTINUE
                ELSE
                  DO 153 IFS=IFSMIN,IFSMAX
                    I3=KFS3(IFS)
                    IF(         KFS1(IFS).LE.I1MAX.AND.
     *                          KFS2(IFS).LE.I2MAX.AND.
     *                 I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN
                      CALL UPDATE()
                    END IF
  153             CONTINUE
                END IF
              ELSE
                IF(NLEX1.LE.I1MAX) THEN
                  DO 154 IFS=IFSMIN,IFSMAX
                    I3=KFS3(IFS)
                    IF(I1MIN.LE.KFS1(IFS).AND.
     *                          KFS2(IFS).LE.I2MAX.AND.
     *                 I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN
                      CALL UPDATE()
                    END IF
  154             CONTINUE
                ELSE
                  DO 155 IFS=IFSMIN,IFSMAX
                    I1=KFS1(IFS)
                    I3=KFS3(IFS)
                    IF(I1MIN.LE.I1.AND.I1.LE.I1MAX.AND.
     *                          KFS2(IFS).LE.I2MAX.AND.
     *                 I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN
                      CALL UPDATE()
                    END IF
  155             CONTINUE
                END IF
              END IF
            END IF
          ELSE
            IF(NLEX2.LE.I2MAX) THEN
              IF(I1MIN.LE.-NLEX1) THEN
                IF(NLEX1.LE.I1MAX) THEN
                  DO 156 IFS=IFSMIN,IFSMAX
                    I3=KFS3(IFS)
                    IF(I2MIN.LE.KFS2(IFS).AND.
     *                 I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN
                      CALL UPDATE()
                    END IF
  156             CONTINUE
                ELSE
                  DO 157 IFS=IFSMIN,IFSMAX
                    I3=KFS3(IFS)
                    IF(         KFS1(IFS).LE.I1MAX.AND.
     *                 I2MIN.LE.KFS2(IFS).AND.
     *                 I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN
                      CALL UPDATE()
                    END IF
  157             CONTINUE
                END IF
              ELSE
                IF(NLEX1.LE.I1MAX) THEN
                  DO 158 IFS=IFSMIN,IFSMAX
                    I3=KFS3(IFS)
                    IF(I1MIN.LE.KFS1(IFS).AND.
     *                 I2MIN.LE.KFS2(IFS).AND.
     *                 I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN
                      CALL UPDATE()
                    END IF
  158             CONTINUE
                ELSE
                  DO 159 IFS=IFSMIN,IFSMAX
                    I1=KFS1(IFS)
                    I3=KFS3(IFS)
                    IF(I1MIN.LE.I1.AND.I1.LE.I1MAX.AND.
     *                 I2MIN.LE.KFS2(IFS).AND.
     *                 I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN
                      CALL UPDATE()
                    END IF
  159             CONTINUE
                END IF
              END IF
            ELSE
              IF(I1MIN.LE.-NLEX1) THEN
                IF(NLEX1.LE.I1MAX) THEN
                  DO 160 IFS=IFSMIN,IFSMAX
                    I2=KFS2(IFS)
                    I3=KFS3(IFS)
                    IF(I2MIN.LE.I2.AND.I2.LE.I2MAX.AND.
     *                 I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN
                      CALL UPDATE()
                    END IF
  160             CONTINUE
                ELSE
                  DO 161 IFS=IFSMIN,IFSMAX
                    I2=KFS2(IFS)
                    I3=KFS3(IFS)
                    IF(         KFS1(IFS).LE.I1MAX.AND.
     *                 I2MIN.LE.I2.AND.I2.LE.I2MAX.AND.
     *                 I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN
                      CALL UPDATE()
                    END IF
  161             CONTINUE
                END IF
              ELSE
                IF(NLEX1.LE.I1MAX) THEN
                  DO 162 IFS=IFSMIN,IFSMAX
                    I2=KFS2(IFS)
                    I3=KFS3(IFS)
                    IF(I1MIN.LE.KFS1(IFS).AND.
     *                 I2MIN.LE.I2.AND.I2.LE.I2MAX.AND.
     *                 I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN
                      CALL UPDATE()
                    END IF
  162             CONTINUE
                ELSE
                  DO 163 IFS=IFSMIN,IFSMAX
                    I1=KFS1(IFS)
                    I2=KFS2(IFS)
                    I3=KFS3(IFS)
                    IF(I1MIN.LE.I1.AND.I1.LE.I1MAX.AND.
     *                 I2MIN.LE.I2.AND.I2.LE.I2MAX.AND.
     *                 I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN
                      CALL UPDATE()
                    END IF
  163             CONTINUE
                END IF
              END IF
            END IF
          END IF
        END IF
      END IF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE FREVOL()
C
C Procedure determining index of the node 'J', if its offset
C respectively to the node 'I' is given.
C This procedure has to be called if the indices of nodes differ from
C the indices of gridpoints.
C Called by LOOP if LOOP is called by GENER.
C
C-----------------------------------------------------------------------
C
C     Common blocks /GRID/ /FS/ /NODE/ are required here:
      INCLUDE 'net.inc'
C     net.inc
      INCLUDE 'netnode.inc'
C     netnode.inc
C
C-----------------------------------------------------------------------
C
      REAL GIANT
      PARAMETER (GIANT=1.E+10)
C
      INTEGER JPOS1,JPOS2,JPOS3,JN1,JN2,JN3,JL1,JL2,JL3,JADR,J
C
C     JPOS1,JPOS2,JPOS3... Position of the gridpoint 'J' within the
C             model.  JPOSI=0,1,2,...,NLI-1.
C     JN1,JN2,JN3... Position of the big brick with 'J' within the
C             model.   JNI=0,1,2,...,NI-1.
C     JL1,JL2,JL3... Position of the gridpoint 'J' within the big brick.
C             JLI=0,1,2,...,LI-1.
C     JADR... Index of the node 'J'.
C     J...    Index of the first node of the big brick.
C
      REAL TTJ,TTIJ
C
C     TTJ...  Arrival time at the node 'J'.
C     TTIJ... Arrival time through the node 'I' to the node 'J'.
C
C.......................................................................
C
C     22 integer operations
      JPOS1=IPOS1+KFS1(IFS)
      JPOS2=IPOS2+KFS2(IFS)
      JPOS3=IPOS3+KFS3(IFS)
      JN1=JPOS1/L1
      JN2=JPOS2/L2
      JN3=JPOS3/L3
      J=IRAM(IND0+1+JN1+N1*(JN2+N2*JN3))
      IF(J.LE.0) RETURN
        JL1=JPOS1-JN1*L1
        JL2=JPOS2-JN2*L2
        JL3=JPOS3-JN3*L3
        JADR=J+JL1+L1*(JL2+L2*JL3)
      GO TO 20
C
C=======================================================================
C
C     
C
      ENTRY UPDATE()
C
C     Procedure updating travel time at node 'J' of the network.
C     Called by LOOP if LOOP is called by GENER.
C
C-----------------------------------------------------------------------
C
C     Index of the node 'J':
      JADR=IPOS+KFS4(IFS)
C
   20 CONTINUE
      TTJ=RAM(ITT0+JADR)
C
      IF(TTJ.LT.0.) RETURN
C       Node 'J' is not in the set A, may be updated.
C
C       Check for crossing an interface:
        IF(LICB) THEN
          IF(IRAM(ICB0+JADR).NE.ICBI) THEN
            IF(KFS5(IFS).GT.1) RETURN
          END IF
        END IF
C
C***********************************************************************
C       Consolidating forward and backward stars (optional):
*       IF(LNFS) THEN
*         IF(KFS5(IFS).GT.IRAM(NFS0+JADR)) RETURN
*       END IF
C       To guarantee the coincidence of forward and backward stars, it
C       is necessary not only to enable the above three statements, but
C       also to submit such template forward stars that each template
C       forward star is a subset of the forward stars of greater sizes.
C***********************************************************************
C
C       Arrival time from the node 'I' to the node 'J':
        TTIJ=TTI+DFS(IFS)*(RAM(IP0+JADR)+PI)
C
C       Test of Bellman's condition for shortest path
        IF(TTIJ.GE.TTJ) RETURN
C
C         Testing if the node 'J' is in the queue
          IF(TTJ.GE.GIANT) THEN
C           Updating queue by adding 'J' as Q(MAXQ)
            MAXQ=MAXQ+1
            IRAM(IPOSQ0+MAXQ)=IPOS+KFS4(IFS)
          END IF
C
C         Updating
          RAM(ITT0+JADR)=TTIJ
          IF(LPRED) THEN
            IRAM(IPRED0+JADR)=IPOS
          END IF
C
C       END of IF
C     END of IF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE SRCREC()
C
C Auxiliary routine to LOOP updating a source or receiver node.
C Called by LOOP if LOOP is called by SOURCE or RECVRS.
C
C-----------------------------------------------------------------------
C
C     Common blocks /GRID/ /FS/ /NAMESR/ /POINTS/ /NODE/ /SRCC/ are
C     required here:
      INCLUDE 'net.inc'
C     net.inc
      INCLUDE 'netnode.inc'
C     netnode.inc
C
C-----------------------------------------------------------------------
C
      EXTERNAL INDX
      INTEGER  INDX
C
      REAL GIANT
      PARAMETER (GIANT=1.E+10)
C
      INTEGER JPOS,JADR
      REAL TTJ,TTIJ
C
C     JPOS... Position of the gridpoint 'J' within the model.
C             JPOS=1,2,3,...,NL1*NL2*NL3.
C     JADR... Index of the node 'J'.
C     TTJ...  Arrival time at the node 'J'.
C     TTIJ... Travel time from the node 'I' to the node 'J'.
C
C.......................................................................
C
C     Index of the node 'J':
      JPOS=IPOS+KFS4(IFS)
      JADR=INDX(JPOS)
C
C     Check for the Fresnel volume:
      IF(JADR.LE.0) RETURN
C
C       Check for a free space:
        IF(RAM(IP0+JADR).GE.GIANT) RETURN
C
          TTJ=RAM(ITT0+JADR)
C
C         Check for crossing an interface:
          IF(LICB) THEN
            IF(IRAM(ICB0+JADR).NE.ICBI) THEN
              IF(KFS5(IFS).GT.1) RETURN
            END IF
          END IF
C
C***********************************************************************
C         Consolidating forward and backward stars (optional):
*         IF(LNFS) THEN
*           IF(KFS5(IFS).GT.IRAM(NFS0+JADR)) RETURN
*         END IF
C***********************************************************************
C
C         Travel time from the node 'I' to the node 'J':
          IF(LN1) THEN
            TTIJ=SQRT((FLOAT(KFS2(IFS))*DSX1-DPOS1)**2
     *               +(FLOAT(KFS2(IFS))*DSX2-DPOS2)**2
     *               +(FLOAT(KFS3(IFS))*DSX3-DPOS3)**2)
     *           *0.5*(RAM(IP0+JADR)+PI)
          ELSE
            TTIJ=SQRT((FLOAT(KFS1(IFS))*DSX1-DPOS1)**2
     *               +(FLOAT(KFS2(IFS))*DSX2-DPOS2)**2
     *               +(FLOAT(KFS3(IFS))*DSX3-DPOS3)**2)
     *           *0.5*(RAM(IP0+JADR)+PI)
          END IF
C
          IF(ISRCI.LE.0) THEN
            IF(TTI+TTIJ.GE.TTJ) RETURN
              IF(TTJ.GE.GIANT) THEN
C               Updating queue by adding 'J' as Q(MAXQ)
                MAXQ=MAXQ+1
                IRAM(IPOSQ0+MAXQ)=IPOS+KFS4(IFS)
              END IF
              RAM(ITT0+JADR)=TTI+TTIJ
              IF(LPRED) THEN
                IRAM(IPRED0+JADR)=ISRCI
              END IF
C           END of IF
          ELSE
            IF(TTJ+TTIJ.GE.TTI) RETURN
              TTI=TTJ+TTIJ
              IPREDR(ISRCI)=JPOS
C           END of IF
          END IF
C
C       END of IF
C     END of IF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE READFS()
C
C Procedure called by SOURCE, to read and construct optimized template
C forward stars.
C
C-----------------------------------------------------------------------
C
C     Common blocks /GRID/ /FS/ are required here:
      INCLUDE 'net.inc'
C     net.inc
C
C-----------------------------------------------------------------------
C
      INTEGER LU1
      PARAMETER(LU1=1)
      CHARACTER*80 FSTAB
C
      INTEGER MD
      PARAMETER (MD=292)
C     MD=292 is good for optimized spherical 3-D forward stars up to
C     size MFSMAX=27.  For minimum values of MD refer to files
C     net.fs2 and net.fs3.
C
      INTEGER ID(MD),JD(MD),KD(MD)
      INTEGER JFS,ND,I1,I2,I3,I,L
C
C     ID(II),JD(II),KD(II)... Vector specifying II-th node of the
C            forward star, in dimensions of a small brick.
C     JFS...  Size of a template forward star currently being read.
C     ND...   Number of nodes stored in the file for a forward star.
C     I1,I2,I3,I,L... Temporary and loop variables.
C
C.......................................................................
C
C     Reading and assembling optimized forward star:
      WRITE(*,'(A)')
     *  '+NET: Reading and assembling the forward star.                '
C
C     Reading nodes of the forward stars:
      IF(NFSMAX.GT.0) THEN
        IF(NL1.EQ.1.OR.NL2.EQ.1.OR.NL3.EQ.1) THEN
          OPEN(LU1,FILE='net.fs2',STATUS='OLD')
        ELSE
          OPEN(LU1,FILE='net.fs3',STATUS='OLD')
        END IF
      END IF
      KFS0(0)=0
      DO 79 JFS=1,NFSMAX
        READ(LU1,*) I,ND
        IF(I.NE.JFS) THEN
C         NET-40
          CALL ERROR
     *           ('NET-40: Template forward star not found in the file')
C         NET-40: Template forward star not found in the file:
C                 The file 'net.fs3' or 'net.fs2' does not contain a
C                 template forward star of the size NFSMAX.  NFSMAX in
C                 main input data file NET, (3),
C                 should be decreased if positive,
C                 adjusted if zero,
C                 or larger forward stars should be added to file
C                 'net.fs3' or 'net.fs2'.
        END IF
        IF(ND.GT.MD) THEN
C         NET-41
          CALL ERROR('NET-41: Too many input nodes of a forward star')
C         NET-41: Too many input nodes of a forward star:
C                 NFSMAX in main input data file NET,
C                 (3), should be decreased if positive,
C                 adjusted if zero,
C                 or the dimension MD in the subroutine READFS should be
C                 increased.
        END IF
        READ(LU1,*) (ID(L),JD(L),KD(L),L=1,ND)
C
C       Assembling whole optimized forward star from the given nodes:
        KFS0(JFS)=KFS0(JFS-1)
        DO 69 L=1,ND
          I1=ID(L)
          I2=JD(L)
          I3=KD(L)
          CALL SETFS(JFS,I1,I2,I3)
   69   CONTINUE
C
   79 CONTINUE
      CLOSE(LU1)
C
C     Printing the table of the numbers of nodes corresponding to the
C     template forward stars of sizes 1,2,...,NFSMAX:
      CALL RSEP3T('FSTAB',FSTAB,' ')
      IF(FSTAB.NE.' ') THEN
        OPEN(LU1,FILE=FSTAB)
        WRITE(LU1,'(A)') '      Size    Number    Sum of'
        WRITE(LU1,'(A)') '     (NFS)  of nodes   numbers'
        WRITE(LU1,'(3I10)')
     *                (JFS,KFS0(JFS)-KFS0(JFS-1),KFS0(JFS),JFS=1,NFSMAX)
        CLOSE(LU1)
        WRITE(*,'(2A)') '+NET: FS table written to file ',FSTAB(1:48)
C       The central node of a forward star is not stored in the memory
C       and is not counted.
      END IF
C
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE MAKEFS(NFSMAK)
      INTEGER NFSMAK
C
C Procedure called by SOURCE and RECVRS, to construct full spherical
C template forward stars.
C
C-----------------------------------------------------------------------
C
C     Common blocks /GRID/ /FS/ are required here:
      INCLUDE 'net.inc'
C     net.inc
C
C-----------------------------------------------------------------------
C
      INTEGER NDIM,NH,IH,I1,I2,I3,I3I3,I
      REAL DH
C
C.......................................................................
C
      IF(NL1.EQ.1.OR.NL2.EQ.1.OR.NL3.EQ.1) THEN
        NDIM=2
      ELSE
        NDIM=3
      END IF
      NH=NFSMAK*NFSMAK+NDIM-1
      DO 11 I=1,NFSMAK
        KFS0(I)=0
   11 CONTINUE
C
      DH=1./SQRT(12.*FLOAT(NH))
      DH=DH/2.
      DO 69 IH=0,NH
        DO 58 I1=INT(SQRT(FLOAT(IH)/3.)-DH+1.),INT(SQRT(FLOAT(IH))+DH)
          I=IH-I1*I1
          DO 57 I2=INT(SQRT(FLOAT(I)/2.)-DH+1.),
     *                                   MIN0(INT(SQRT(FLOAT(I))+DH),I1)
            I3I3=I-I2*I2
            I3=INT(SQRT(FLOAT(I3I3))+0.500)
            IF(I3*I3.EQ.I3I3.AND.I3.LE.I2) THEN
C             Assembling the forward star:
              CALL SETFS(NFSMAK,I3,I2,I1)
            END IF
   57     CONTINUE
   58   CONTINUE
   69 CONTINUE
C
      DO 71 I=NFSMAK,NFSMAX
        KFS0(I)=KFS0(NFSMAK)
   71 CONTINUE
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE SETFS(JFS,I1,I2,I3)
      INTEGER JFS,I1,I2,I3
C
C Subroutine designed to set up all possible mirrorings of an edge
C of the forward star.  Auxiliary subroutine to READFS and MAKEFS.
C
C-----------------------------------------------------------------------
C
C     Common blocks /GRID/ /FS/ are required here:
      INCLUDE 'net.inc'
C     net.inc
C
C-----------------------------------------------------------------------
C
      INTEGER NLEX1,NLEX2,NLEX3,LFS,I
      REAL DIM
C
C     NLEX1,NLEX2,NLEX3...
C     LFS...  Number of template forward star nodes already stored.
C             NODES KFS0(JFS),...,LFS are being mirrored with respect to
C             the central node, while storing the results of mirroring.
C     I...    Loop variable.
C     DIM...  Dimension of the model (2-D or 3-D).
C
C.......................................................................
C
      LFS=KFS0(JFS)
      IF(NL1.EQ.1.OR.NL2.EQ.1.OR.NL3.EQ.1) THEN
        DIM=2.
      ELSE
        DIM=3.
      END IF
      NLEX1=MIN0(JFS,NL1-1)
      NLEX2=MIN0(JFS,NL2-1)
      NLEX3=MIN0(JFS,NL3-1)
C
      IF(I1.LE.NLEX1.AND.I2.LE.NLEX2.AND.I3.LE.NLEX3) THEN
        CALL STORFS(LFS,I1,I2,I3)
      END IF
      IF(I1.NE.I3) THEN
        IF(I2.LE.NLEX1.AND.I3.LE.NLEX2.AND.I1.LE.NLEX3) THEN
          CALL STORFS(LFS,I2,I3,I1)
        END IF
        IF(I3.LE.NLEX1.AND.I1.LE.NLEX2.AND.I2.LE.NLEX3) THEN
          CALL STORFS(LFS,I3,I1,I2)
        END IF
        IF(I1.NE.I2.AND.I2.NE.I3) THEN
          IF(I3.LE.NLEX1.AND.I2.LE.NLEX2.AND.I1.LE.NLEX3) THEN
            CALL STORFS(LFS,I3,I2,I1)
          END IF
          IF(I1.LE.NLEX1.AND.I3.LE.NLEX2.AND.I2.LE.NLEX3) THEN
            CALL STORFS(LFS,I1,I3,I2)
          END IF
          IF(I2.LE.NLEX1.AND.I1.LE.NLEX2.AND.I3.LE.NLEX3) THEN
            CALL STORFS(LFS,I2,I1,I3)
          END IF
        END IF
      END IF
      DO 61 I=KFS0(JFS)+1,LFS
        IF(KFS1(I).NE.0) THEN
          CALL STORFS(LFS,-KFS1(I), KFS2(I), KFS3(I))
        END IF
   61 CONTINUE
      DO 62 I=KFS0(JFS)+1,LFS
        IF(KFS2(I).NE.0) THEN
          CALL STORFS(LFS, KFS1(I),-KFS2(I), KFS3(I))
        END IF
   62 CONTINUE
      DO 63 I=KFS0(JFS)+1,LFS
        IF(KFS3(I).NE.0) THEN
          CALL STORFS(LFS, KFS1(I), KFS2(I),-KFS3(I))
        END IF
   63 CONTINUE
      DO 64 I=KFS0(JFS)+1,LFS
        KFS4(I)=KFS1(I)+NL1*(KFS2(I)+NL2*KFS3(I))
        KFS5(I)=INT(SQRT(AMAX1(FLOAT(KFS1(I))**2
     *                        +FLOAT(KFS2(I))**2
     *                        +FLOAT(KFS3(I))**2-DIM+1.,1.))+0.999)
        DFS(I)=SQRT((FLOAT(KFS1(I))*ASX1)**2
     *             +(FLOAT(KFS2(I))*ASX2)**2
     *             +(FLOAT(KFS3(I))*ASX3)**2)/2.
   64 CONTINUE
C
      KFS0(JFS)=LFS
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE STORFS(LFS,I1,I2,I3)
      INTEGER LFS,I1,I2,I3
C
C Subroutine designed to store one edge of the forward star in the
C memory.  Auxiliary subroutine to SETFS.
C
C-----------------------------------------------------------------------
C
C     Common block /FS/ is required here:
      INCLUDE 'net.inc'
C     net.inc
C
C-----------------------------------------------------------------------
C
      LFS=LFS+1
      IF(LFS.GT.MFS) THEN
C       NET-42
        CALL ERROR('NET-42: Too big forward star')
C       NET-42: Too big forward star:
C               NFSMAX in main input data file NET,
C               (3), should be decreased if positive, adjusted if zero,
C               or the dimension MFS in the common block /FS/
C               in net.inc should be increased.
C               The minimum value of the dimension MFS may be determined
C               as follows:
C                 (a) Choose MFS as large as possible.
C                 (b) Add parameter FSTAB='filename' into the input SEP
C                     parameter file and adjust NFSMAX.
C                 (c) Compile and run the program.
C                 (d) Look at file 'filename' just having been
C                     generated.  The last integer in the file is the
C                     minimum dimension MFS corresponding to the given
C                     value of NFSMAX.
C                 (e) Update MFS, restore NFSMAX and remove parameter
C                     FSTAB='filename' from the data.
      END IF
      KFS1(LFS)=I1
      KFS2(LFS)=I2
      KFS3(LFS)=I3
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE ONERAY(LTRACE,LRAYS,LEND,LURAYS,IREFL,
     *                IPREDE,IEND1,IEND2,IEND3,IENDA,TEND,ERR,TERR,ISRC)
      LOGICAL LTRACE,LRAYS,LEND
      INTEGER LURAYS,IREFL,IPREDE,IEND1,IEND2,IEND3,IENDA,ISRC
      REAL TEND,ERR(*),TERR
C
C Subroutine called by ERRORS and TRACER, to perform backward ray
C tracing.
C
C Input:
C     LTRACE,LRAYS,LEND,LURAYS,IREFL,IPREDE,IEND1,IEND2,IEND3,IENDA,TEND
C     IPREDE... Must be positive.
C
C Output:
C     TERR
C
C-----------------------------------------------------------------------
C
C     Common blocks /GRID/ /NAMESR/ /POINTS/ /FILES/ are required here:
      INCLUDE 'net.inc'
C     net.inc
C
C-----------------------------------------------------------------------
C
      EXTERNAL INDX
      INTEGER  INDX
C
      INTEGER LU0,LU1,LU
      PARAMETER(LU0=10)
      PARAMETER(LU1=1)
C
      INTEGER MTMP
      PARAMETER (MTMP=100)
      INTEGER KTMP(MTMP),NTMP,ITMP
      REAL TMPMIN(MTMP),TMPMAX(MTMP)
C
      INTEGER IPREDI,IPNEW,NF,N123,L123,IPOS1,IPOS2,IPOS3,IPOS,IADR
      INTEGER IOLD1,IOLD2,IOLD3,ICBOLD,NFOLD,K
      REAL DPOS1,DPOS2,DPOS3,X1,X2,X3,T,TOLD,POLD,PAUX,ERRMIN,ERRMAX
C
C     IPOS1,IPOS2,IPOS3,IPOS,IADR... Position and index (address) of the
C             nearest network node.
C     DPOS1,DPOS2,DPOS3... Offset from the nearest network node.
C     IPREDI..Predecessor of the current node.
C     IOLD1,IOLD2,IOLD3,ICBOLD,NFOLD,TOLD,POLD... Parameters of the last
C             network node along the ray.
C     X1,X2,X3,T,PAUX... Temporary storage locations.
C
C.......................................................................
C
      IPREDI=IPREDE
      NTMP=0
      IF(.NOT.LTRACE.OR.LEND) THEN
        IF(LICB) THEN
          ICBOLD=IRAM(ICB0+IENDA)
        ELSE
          ICBOLD=1
        END IF
        IOLD1=IEND1
        IOLD2=IEND2
        IOLD3=IEND3
        POLD=RAM(IP0+IENDA)
C       POLD is a constant approximation, not interpolation.
        TOLD=TEND
        NFOLD=IREFL
        ERRMIN=0.
        ERRMAX=0.
      END IF
      IF(IPREDI.GE.0) THEN
        DO 68 NF=IREFL,0,-1
          IF(IPREDI.EQ.0) THEN
C           Endpoint of the ray at the reflector:
C           This is not possible at the receiver.
C           Thus, this is possible only if called by subroutine ERRORS.
C           Then IOLD1,IOLD2,IOLD3 are defined.
            IPREDI=1+IOLD1+(NL1*IOLD2+NL2*IOLD3)
            GO TO 66
          ENDIF
          IF(NF.LT.IREFL) THEN
            LU=LU0+NF
            IF(L4.NE.0) THEN
C             Reading the index array:
              N123=N1*N2*N3
              L123=L1*L2*L3
              DO 11 K=1,N123
                IRAM(IND0+K)=K
   11         CONTINUE
              CALL RARRAI(LU1,FIND(NF),'FORMATTED',
     *                    .FALSE.,0,N123,IRAM(IND0+1))
              DO 12 K=1,N123
                IRAM(IND0+K)=(IRAM(IND0+K)-1)*L123+1
   12         CONTINUE
            END IF
          END IF
   50     CONTINUE
            IADR=INDX(IPREDI)
            IF(NF.LT.IREFL) THEN
              READ(LU,REC=IADR) T,IPNEW
            ELSE
              T=RAM(ITT0+IADR)
              IPNEW=IRAM(IPRED0+IADR)
            END IF
            IPOS1=IPREDI-1
            IPOS2=IPOS1/NL1
            IPOS3=IPOS2/NL2
            IPOS1=IPOS1-IPOS2*NL1
            IPOS2=IPOS2-IPOS3*NL2
            IF(LTRACE.AND.LRAYS) THEN
              IF(LN1) THEN
                X1=X1MIN+(FLOAT(IPOS2)+0.5)*DSX1
              ELSE
                X1=X1MIN+(FLOAT(IPOS1)+0.5)*DSX1
              END IF
              X2=X2MIN+(FLOAT(IPOS2)+0.5)*DSX2
              X3=X3MIN+(FLOAT(IPOS3)+0.5)*DSX3
              WRITE(LURAYS,'(4F12.6,A)') X1,X2,X3,T,' /'
            END IF
            IF(.NOT.LTRACE.OR.LEND) THEN
C             Travel time error estimate:
              CALL SETERR(NF,IPREDI,IPOS1,IPOS2,IPOS3,IADR,T,NFOLD,
     *                 IOLD1,IOLD2,IOLD3,ICBOLD,POLD,TOLD,ERRMIN,ERRMAX)
              IF(.NOT.LTRACE) THEN
                IF(NF.EQ.IREFL) THEN
                  IF(ERR(IADR).GE.0.) THEN
                    ERRMIN=ERRMIN-ERR(IADR)
                    ERRMAX=ERRMAX+ERR(IADR)
                    GO TO 70
                  ELSE
                    IF(NTMP.LT.MTMP) THEN
                      NTMP=NTMP+1
                      KTMP(NTMP)=IADR
                      TMPMIN(NTMP)=ERRMIN
                      TMPMAX(NTMP)=ERRMAX
                    END IF
                  END IF
                END IF
              END IF
            END IF
            IF(IPNEW.EQ.0) THEN
C             Reflector:
C             The same point at the reflector is considered twice,
C             i.e. the ray segment of zero length is situated at the
C             reflector.
              GO TO 66
            ELSE IF(IPNEW.LT.0) THEN
C             Source point:
              IPREDI=IPNEW
              GO TO 69
            ELSE
              IPREDI=IPNEW
            ENDIF
          GO TO 50
   66     CONTINUE
C         Node IPREDI is situated at a reflector.
*         AUX=POLD*AMAX1(ASX1,ASX2,ASX3)/2.
*         ERRMIN=ERRMIN-AUX
*         ERRMAX=ERRMIN+AUX
   68   CONTINUE
   69   CONTINUE
      END IF
C     Node IPREDI is a source point.
      ISRC=-IPREDI
C
C     Source point:
      IF(.NOT.LTRACE.OR.LEND) THEN
C       Travel time error estimate:
        CALL SLOW(X1S(ISRC),X2S(ISRC),X3S(ISRC),
     *               DPOS1,DPOS2,DPOS3,IPOS1,IPOS2,IPOS3,IPOS,IADR,PAUX)
        CALL SETERR(IREFL,IPOS,IPOS1,IPOS2,IPOS3,IADR,TTS(ISRC),NFOLD,
     *                 IOLD1,IOLD2,IOLD3,ICBOLD,POLD,TOLD,ERRMIN,ERRMAX)
        ERRMIN=ERRMIN-TTSERR(ISRC)
        ERRMAX=ERRMAX+TTSERR(ISRC)
      END IF
C
   70 CONTINUE
C
      IF(.NOT.LTRACE) THEN
C       Rays are not traced, the end node is a gridpoint,
C       thus IENDA is defined.
        ERR(IENDA)=AMAX1(ERRMAX,-ERRMIN)
        DO 71 ITMP=1,NTMP
          ERR(KTMP(ITMP))=
     *                   AMAX1(ERRMAX-TMPMAX(ITMP),-ERRMIN+TMPMIN(ITMP))
   71   CONTINUE
      ELSE IF(LEND) THEN
        TERR=AMAX1(ERRMAX,-ERRMIN)
      END IF
C
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE SETERR(NF,IPOS,IPOS1,IPOS2,IPOS3,IADR,T,NFOLD,
     *                 IOLD1,IOLD2,IOLD3,ICBOLD,POLD,TOLD,ERRMIN,ERRMAX)
      INTEGER NF,IPOS,IPOS1,IPOS2,IPOS3,IADR
      INTEGER NFOLD,IOLD1,IOLD2,IOLD3,ICBOLD
      REAL T,POLD,TOLD,ERRMIN,ERRMAX
C
C Subroutine called by ONERAY, to estimate and accumulate the
C arrival-time errors during backward ray tracing.
C
C Attention:
C     In this version, in a case of reflections, the slowness field is
C     correct only if it is the same for the direct wave and for all the
C     reflected waves.  When the rays are traced, only time fields and
C     predecessors are alternated in the memory, unlike slowness fields.
C
C Input:
C     NF,IPOS,IPOS1,IPOS2,IPOS3,IADR,T,
C     NFOLD,IOLD1,IOLD2,IOLD3,ICBOLD,POLD,TOLD,ERRMIN,ERRMAX
C
C Output:
C     NFOLD,IOLD1,IOLD2,IOLD3,ICBOLD,POLD,TOLD... Values at the old node
C             are replaced by the values at the new node.
C     ERRMIN,ERRMAX... Input values increased by the increment of the
C             error bounds between the old node and the new node.
C
C-----------------------------------------------------------------------
C
C     Common blocks /GRID/ /FS/ are required here:
      INCLUDE 'net.inc'
C     net.inc
C
C-----------------------------------------------------------------------
C
      INTEGER ICBNEW,NFS1,IERR
      REAL UIJ(6),C1,C2MIN,C2MAX,ERR3,PMEAN
      REAL AUX1,AUX2,AUX3
C
C.......................................................................
C
      CALL OPTMAT(IPOS,IPOS1,IPOS2,IPOS3,IADR,IERR,UIJ)
      IF(LNFS) THEN
        IF(NF.EQ.NREFL) THEN
          NFS1=IRAM(NFS0+IADR)
        ELSE
          CALL OPTNFS(IPOS,IPOS1,IPOS2,IPOS3,IADR,NFS1,C1,C2MIN,C2MAX)
        END IF
      ELSE
        NFS1=NFSMAX
      END IF
C
C     Calculating the coefficient of the network-geometry error:
      IF(NL1.EQ.1) THEN
        C1=((ASX2**2+ASX3**2)/AMIN1(ASX2,ASX3)**2-1.)/8.
      ELSE IF(NL2.EQ.1) THEN
        C1=((ASX1**2+ASX3**2)/AMIN1(ASX1,ASX3)**2-1.)/8.
      ELSE IF(NL3.EQ.1) THEN
        C1=((ASX1**2+ASX2**2)/AMIN1(ASX1,ASX2)**2-1.)/8.
      ELSE
        C1=((ASX1**2+ASX2**2+ASX3**2)/AMIN1(ASX1,ASX2,ASX3)**2-1.)/8.
      END IF
C
C     Calculating the relative heterogeneity error:
      AUX1=FLOAT(IPOS1-IOLD1)
      AUX2=FLOAT(IPOS2-IOLD2)
      AUX3=FLOAT(IPOS3-IOLD3)
      AUX1=AUX1*(UIJ(1)*AUX1+2.*(UIJ(2)*AUX2+UIJ(4)*AUX3))
     *    +AUX2*(UIJ(3)*AUX2+2.*UIJ(5)*AUX3)+AUX3*UIJ(6)*AUX3
      AUX1=AUX1/(12.*RAM(IP0+IADR))
C
C     Geometry and heterogeneity increments of the error bounds:
      ERRMIN=ERRMIN+(TOLD-T)* AUX1
      ERRMAX=ERRMAX+(TOLD-T)*(AUX1+C1/FLOAT(NFS1*NFS1+1))
C
C     Error due to structural interfaces:
      IF(LICB) THEN
        ICBNEW=IRAM(ICB0+IADR)
      ELSE
        ICBNEW=1
      END IF
      IF(ICBNEW.NE.ICBOLD) THEN
C       For ABS(DSX1)=ABS(DSX2)=ABS(DSX3)=DSX:
C       D(ERRMIN)=-DSX* D(P)*SQRT(N)/2
C       D(ERRMAX)=DSX*MIN(P,
C                (D(P)*SQRT(2)-SQRT(MIN(P)**2+((SQRT(N)-SQRT(2))*P)**2))
        AUX1=AMAX1(ASX1,ASX2,ASX3)
        AUX3=AMIN1(ASX1,ASX2,ASX3)
        AUX2=ASX1+ASX2+ASX3-AUX1-AUX3
        PMEAN=(RAM(IP0+IADR)+POLD)/2.
        ERRMIN=ERRMIN-SQRT(AUX1*AUX1+AUX2*AUX2+AUX3*AUX3)
     *                *ABS(RAM(IP0+IADR)-POLD)/2.
        ERR3=AMIN1(RAM(IP0+IADR),POLD)
C       In the worst case, interface is perpendicular to the greatest
C       grid step AUX1.
        IF(AUX3.LE.0.) THEN
C         2-D:
C         Now, ERR3 is the time derivative in the direction AUX2.
          ERR3=PMEAN*SQRT(AUX1*AUX1+AUX2*AUX2)-ERR3*AUX2
        ELSE
C         3-D:
          ERR3=ERR3*ERR3-(PMEAN*( SQRT(AUX1*AUX1+AUX2*AUX2+AUX3*AUX3)
     *                           -SQRT(AUX1*AUX1+AUX3*AUX3)   )/AUX2)**2
          IF(ERR3.GT.0.) THEN
            ERR3=SQRT(ERR3)
          ELSE
            ERR3=0.
          END IF
C         Now, ERR3 is the time derivative in the direction AUX3.
          ERR3=PMEAN*SQRT(AUX1*AUX1+AUX3*AUX3)-ERR3*AUX3
        END IF
        ERR3=AMIN1(ERR3,PMEAN*AUX1)
        ERRMAX=ERRMAX+ERR3
      END IF
C
C     Error due to reflections:
      IF(NFOLD.NE.NF) THEN
        AUX1=(RAM(IP0+IADR)+POLD)*AMAX1(ASX1,ASX2,ASX3)/2.
        ERRMIN=ERRMIN-AUX1
        ERRMAX=ERRMIN+AUX1
        NFOLD=NF
      END IF
C
      ICBOLD=ICBNEW
      IOLD1=IPOS1
      IOLD2=IPOS2
      IOLD3=IPOS3
      POLD=RAM(IP0+IADR)
      TOLD=T
C
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE OPTNFS
     *               (IPOS,IPOS1,IPOS2,IPOS3,IADR,NFSOPT,C1,C2MIN,C2MAX)
      INTEGER IPOS,IPOS1,IPOS2,IPOS3,IADR,NFSOPT
      REAL C1,C2MIN,C2MAX
C
C Subroutine called by SOURCE and SETERR, designed to estimate the
C optimum size of a forward star at the given node.
C
C Input:
C     IPOS,IPOS1,IPOS2,IPOS3... Indices of the given gridpoint,
C             describing the position of the network node situated at
C             the gridpoint.
C     IADR... Index of the network node.
C
C Output:
C     NFSOPT..Optimum size of the forward star at the given node.
C             Default of NFSOPT=1 is set if the optimum size cannot be
C             determined.
C     C1,C2MIN,C2MAX... Coefficients of the travel time error estimate.
C
C-----------------------------------------------------------------------
C
C     Common blocks /GRID/ /FS/ are required here:
      INCLUDE 'net.inc'
C     net.inc
C
C-----------------------------------------------------------------------
C
      REAL GIANT
      PARAMETER (GIANT=1.E+10)
C
      INTEGER IERR
      REAL C1SAVE,C2,UIJ(6),AUX0
      SAVE C1SAVE
      DATA C1SAVE/0./
C
C     C1SAVE..Coefficient C1.
C     C2...   Maximum absolute value of coefficients C2MIN, C2MAX.
C     UIJ...  Symmetric matrix describing the level of heterogeneity.
C     AUX0... Dummy storage location.
C
C.......................................................................
C
C     Calculating the coefficient of the network-geometry error:
      IF(C1SAVE.EQ.0.) THEN
        IF(NL1.EQ.1) THEN
          C1SAVE=((ASX2**2+ASX3**2)/AMIN1(ASX2,ASX3)**2-1.)/8.
        ELSE IF(NL2.EQ.1) THEN
          C1SAVE=((ASX1**2+ASX3**2)/AMIN1(ASX1,ASX3)**2-1.)/8.
        ELSE IF(NL3.EQ.1) THEN
          C1SAVE=((ASX1**2+ASX2**2)/AMIN1(ASX1,ASX2)**2-1.)/8.
        ELSE
          C1SAVE=
     *        ((ASX1**2+ASX2**2+ASX3**2)/AMIN1(ASX1,ASX2,ASX3)**2-1.)/8.
        END IF
      END IF
C
C     Coefficient of the network-geometry error (rel.error=C1/NFS**2):
      C1=C1SAVE
C
C     Check for a free space:
      IF(RAM(IP0+IADR).GE.GIANT) THEN
        NFSOPT=0
        C2MIN=0.
        C2MAX=0.
        RETURN
      END IF
C
      CALL OPTMAT(IPOS,IPOS1,IPOS2,IPOS3,IADR,IERR,UIJ)
C     IERR=0,1,...,16.  Free space: IERR=16.
      IF(IERR.GE.16) THEN
C       Optimum size of a forward star cannot be determined:
        NFSOPT=1
        RETURN
      END IF
C
C     Coefficient of the heterogeneity error (rel.error=C2*NFS**2):
      CALL EIGEN(UIJ,AUX0,3,1)
      C2=12.*RAM(IP0+IADR)
      C2MAX=AMAX1(UIJ(1),0.)/C2
      C2MIN=AMIN1(UIJ(6),0.)/C2
      C2=AMAX1(C2MAX,-C2MIN-C2MAX)
C
C     Maximum size of the forward star:
      IF(NFSMAX.GT.0) THEN
        NFSOPT=NFSMAX
      ELSE
        NFSOPT=999999
      END IF
C
C     Optimum size of the forward star:
      IF(C2.GT.0.) THEN
        NFSOPT=MAX0(1,MIN0(INT(SQRT(SQRT(C1/C2))+0.5),NFSOPT))
      END IF
C
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE OPTMAT(IPOS,IPOS1,IPOS2,IPOS3,IADR,IERR,UIJ)
      INTEGER IPOS,IPOS1,IPOS2,IPOS3,IADR,IERR
      REAL UIJ(6)
C
C Auxiliary subroutine to SETERR and OPTNFS, designed to calculate the
C symmetric matrix describing the level of heterogeneity at the given
C node.  The matrix is composed of the first and second slowness
C derivatives.
C
C Input:
C     IPOS,IPOS1,IPOS2,IPOS3... Indices of the given gridpoint,
C             describing the position of the network node situated at
C             the gridpoint.
C     IADR... Index of the network node.
C
C Output:
C     IERR... IERR=0 if the matrix is determined.
C     UIJ...  Symmetric matrix describing the level of heterogeneity.
C
C-----------------------------------------------------------------------
C
C     Common block /GRID/ is required here:
      INCLUDE 'net.inc'
C     net.inc
C
C-----------------------------------------------------------------------
C
      EXTERNAL INDX
      INTEGER  INDX
      INTEGER ISHIFT(3),ICBOPT,I1,I2,I3,IP,IP1,IP2,IP3,IAUX
      INTEGER IP1OLD,IP2OLD,IP3OLD,IEROLD
      DATA ISHIFT/0,1,-1/
C
C.......................................................................
C
      ICBOPT=-999
      CALL TRYMAT(IPOS,IPOS1,IPOS2,IPOS3,IADR,ICBOPT,IERR,UIJ)
C     ICBOPT unchanged: node not situated in the computational volume.
C     IERR.GT.0: matrix Uij has not been determined.
C
C     Check whether matrix Uij is determined:
      IF(ICBOPT.NE.-999.AND.IERR.GT.0) THEN
        IEROLD=IERR
        IP1OLD=IPOS1
        IP2OLD=IPOS2
        IP3OLD=IPOS3
C
C       Looking for matrix Uij around:
        DO 13 I3=1,3
          IP3=IPOS3+ISHIFT(I3)
C         Check for the model volume:
          IF(0.LE.IP3.AND.IP3.LT.NL3) THEN
            DO 12 I2=1,3
              IP2=IPOS2+ISHIFT(I2)
C             Check for the model volume:
              IF(0.LE.IP2.AND.IP2.LT.NL2) THEN
                DO 11 I1=1,3
                  IF(I1.NE.1.OR.I2.NE.1.OR.I3.NE.1) THEN
                    IP1=IPOS1+ISHIFT(I1)
C                   Check for the model volume:
                    IF(0.LE.IP1.AND.IP1.LT.NL1) THEN
                      IP=1+IP1+NL1*(IP2+NL2*IP3)
                      IAUX=INDX(IP)
                      CALL TRYMAT(IP,IP1,IP2,IP3,IAUX,ICBOPT,IERR,UIJ)
                      IF(IERR.EQ.0) THEN
C                       Matrix Uij is determined
                        RETURN
                      ELSE IF(IERR.LT.IEROLD) THEN
                        IEROLD=IERR
                        IP1OLD=IP1
                        IP2OLD=IP2
                        IP3OLD=IP3
                      END IF
                    END IF
                  END IF
   11           CONTINUE
              END IF
   12       CONTINUE
          END IF
   13   CONTINUE
C
        IP=1+IP1OLD+NL1*(IP2OLD+NL2*IP3OLD)
        IAUX=INDX(IP)
        CALL TRYMAT(IP,IP1OLD,IP2OLD,IP3OLD,IAUX,ICBOPT,IERR,UIJ)
      END IF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE TRYMAT(IPOS,IPOS1,IPOS2,IPOS3,IADR,ICBOPT,IERR,UIJ)
      INTEGER IPOS,IPOS1,IPOS2,IPOS3,IADR,ICBOPT,IERR
      REAL UIJ(6)
C
C Auxiliary subroutine to OPTMAT, to try the calculation of the matrix
C describing the level of heterogeneity at the given node.
C The matrix is composed of the first and second slowness derivatives.
C
C Input:
C     IPOS,IPOS1,IPOS2,IPOS3... Indices of the given gridpoint,
C             describing the position of the network node situated at
C             the gridpoint.
C     IADR... Index of the network node.
C
C Output:
C     ICBOPT..
C     IERR... IERR=0 if the matrix is determined.
C     UIJ...  Symmetric matrix describing the level of heterogeneity.
C
C-----------------------------------------------------------------------
C
C     Common block /GRID/ is required here:
      INCLUDE 'net.inc'
C     net.inc
C
C-----------------------------------------------------------------------
C
      REAL GIANT
      PARAMETER (GIANT=1.E+10)
C
      EXTERNAL INDX
      INTEGER  INDX
      INTEGER IA0,IA2,IERR12,IERR13,IERR23
      REAL U1,U2,U3,U11,U12,U22,U13,U23,U33,AUX0,AUX1,AUX2,AUX3
C
C     IERR12,IERR13,IERR23... Switches for evaluation of mixed second
C             derivatives.
C     U1,U2,U3,U11,U12,U22,U13,U23,U33... Slowness derivatives.
C     AUX0,AUX1,AUX2,AUX3... Temporary storage locations.
C
C.......................................................................
C
C     IERR=0,1,...,16.  IERR=16 if slowness is not defined at the point.
C     Otherwise, slowness derivatives cannot be determined along IERR/4
C     axes.  In addition, MOD(IERR,4) mixed derivatives cannot be
C     determined.  Initially:
      IERR=16
      IERR12=2
      IERR13=2
      IERR23=2
C
C     2-D model and default in 3-D:
      U1 =0.
      U2 =0.
      U3 =0.
      U11=0.
      U12=0.
      U22=0.
      U13=0.
      U23=0.
      U33=0.
C
C     Check for the computational volume at the given node:
      IF(IADR.EQ.0) THEN
        RETURN
      END IF
C
C     Determination of the geological block at the given node:
      IF(ICBOPT.EQ.-999) THEN
        IF(LICB) THEN
          ICBOPT=IRAM(ICB0+IADR)
        ELSE
          ICBOPT=1
        END IF
      END IF
C
C     Check for the geological block at the given node:
      IF(LICB) THEN
        IF(IRAM(ICB0+IADR).NE.ICBOPT) THEN
          RETURN
        END IF
      END IF
C     Check for a free space:
      AUX1=RAM(IP0+IADR)
      IF(AUX1.GE.GIANT) THEN
        RETURN
      END IF
      IERR=IERR-4
C
C     First and second partial slowness derivatives along the X1 axis:
      IF(NL1.GT.1) THEN
C       3-D model:
C       Check for the model volume:
        IF(IPOS1.LT.1.OR.NL1-2.LT.IPOS1) THEN
          GO TO 20
        END IF
C       Check for the computational volume:
        IA0=INDX(IPOS-1)
        IA2=INDX(IPOS+1)
        IF(IA0.EQ.0) THEN
          GO TO 20
        END IF
        IF(IA2.EQ.0) THEN
          GO TO 20
        END IF
C       Check for the geological block:
        IF(LICB) THEN
          IF(IRAM(ICB0+IA0).NE.ICBOPT) THEN
            GO TO 20
          END IF
          IF(IRAM(ICB0+IA2).NE.ICBOPT) THEN
            GO TO 20
          END IF
        END IF
C       Check for a free space:
        AUX0=RAM(IP0+IA0)
        AUX2=RAM(IP0+IA2)
        IF(AUX0.GE.GIANT) THEN
          GO TO 20
        END IF
        IF(AUX2.GE.GIANT) THEN
          GO TO 20
        END IF
C       Partial slowness derivatives scaled by the grid step:
        U1 =(AUX2-AUX0)/2.
        U11=AUX2-2.*AUX1+AUX0
        IERR12=IERR12-1
        IERR13=IERR13-1
      END IF
      IERR=IERR-4
C
C     First and second partial slowness derivatives along the X2 axis:
   20 CONTINUE
      IF(NL2.GT.1) THEN
C       3-D model:
C       Check for the model volume:
        IF(IPOS2.LT.1.OR.NL2-2.LT.IPOS2) THEN
          GO TO 30
        END IF
C       Check for the computational volume:
        IA0=INDX(IPOS-NL1)
        IA2=INDX(IPOS+NL1)
        IF(IA0.EQ.0) THEN
          GO TO 30
        END IF
        IF(IA2.EQ.0) THEN
          GO TO 30
        END IF
C       Check for the geological block:
        IF(LICB) THEN
          IF(IRAM(ICB0+IA0).NE.ICBOPT) THEN
            GO TO 30
          END IF
          IF(IRAM(ICB0+IA2).NE.ICBOPT) THEN
            GO TO 30
          END IF
        END IF
C       Check for a free space:
        AUX0=RAM(IP0+IA0)
        AUX2=RAM(IP0+IA2)
        IF(AUX0.GE.GIANT) THEN
          GO TO 30
        END IF
        IF(AUX2.GE.GIANT) THEN
          GO TO 30
        END IF
C       Partial slowness derivatives scaled by the grid step:
        U2 =(AUX2-AUX0)/2.
        U22=AUX2-2.*AUX1+AUX0
        IERR12=IERR12-1
        IERR23=IERR23-1
      END IF
      IERR=IERR-4
C
C     First and second partial slowness derivatives along the X3 axis:
   30 CONTINUE
      IF(NL3.GT.1) THEN
C       3-D model:
C       Check for the model volume:
        IF(IPOS3.LT.1.OR.NL3-2.LT.IPOS3) THEN
          GO TO 40
        END IF
C       Check for the computational volume:
        IA0=INDX(IPOS-NL1*NL2)
        IA2=INDX(IPOS+NL1*NL2)
        IF(IA0.EQ.0) THEN
          GO TO 40
        END IF
        IF(IA2.EQ.0) THEN
          GO TO 40
        END IF
C       Check for the geological block:
        IF(LICB) THEN
          IF(IRAM(ICB0+IA0).NE.ICBOPT) THEN
            GO TO 40
          END IF
          IF(IRAM(ICB0+IA2).NE.ICBOPT) THEN
            GO TO 40
          END IF
        END IF
C       Check for a free space:
        AUX0=RAM(IP0+IA0)
        AUX2=RAM(IP0+IA2)
        IF(AUX0.GE.GIANT) THEN
          GO TO 40
        END IF
        IF(AUX2.GE.GIANT) THEN
          GO TO 40
        END IF
C       Partial slowness derivatives scaled by the grid step:
        U3 =(AUX2-AUX0)/2.
        U33=AUX2-2.*AUX1+AUX0
        IERR13=IERR13-1
        IERR23=IERR23-1
      END IF
      IERR=IERR-4
C
C     Mixed second partial slowness derivatives:
   40 CONTINUE
      IF(NL1.GT.1.AND.NL2.GT.1.AND.IERR12.EQ.0) THEN
        CALL MIXDER(IPOS,  1,NL1    ,ICBOPT,U12,IERR)
      END IF
      IF(NL1.GT.1.AND.NL3.GT.1.AND.IERR13.EQ.0) THEN
        CALL MIXDER(IPOS,  1,NL1*NL2,ICBOPT,U13,IERR)
      END IF
      IF(NL2.GT.1.AND.NL3.GT.1.AND.IERR23.EQ.0) THEN
        CALL MIXDER(IPOS,NL1,NL1*NL2,ICBOPT,U23,IERR)
      END IF
C
C     Symmetric matrix describing the level of heterogeneity:
      AUX0=0.5/RAM(IP0+IADR)
      AUX1=U1*AUX0
      AUX2=U2*AUX0
      AUX3=U3*AUX0
      AUX0=AUX1*U1+AUX2*U2+AUX3*U3
      UIJ(1)=U11-AUX1*U1+AUX0
      UIJ(2)=U12-AUX1*U2
      UIJ(3)=U22-AUX2*U2+AUX0
      UIJ(4)=U13-AUX1*U3
      UIJ(5)=U23-AUX2*U3
      UIJ(6)=U33-AUX3*U3+AUX0
C
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE MIXDER(IPOS0,IPOS1,IPOS2,ICBREF,DERMIX,IERR)
      INTEGER IPOS0,IPOS1,IPOS2,ICBREF,IERR
      REAL DERMIX
C
C Auxiliary subroutine to TRYMAT, to calculate mixed second partial
C slowness derivatives.
C
C-----------------------------------------------------------------------
C
C     Common block /GRID/ is required here:
      INCLUDE 'net.inc'
C     net.inc
C
C-----------------------------------------------------------------------
C
      REAL GIANT
      PARAMETER (GIANT=1.E+10)
C
      EXTERNAL INDX
      INTEGER  INDX
      INTEGER IA00,IA20,IA02,IA22
C
C.......................................................................
C
C     Addresses of the corner gridpoints:
      IA00=INDX(IPOS0-IPOS1-IPOS2)
      IA20=INDX(IPOS0+IPOS1-IPOS2)
      IA02=INDX(IPOS0-IPOS1+IPOS2)
      IA22=INDX(IPOS0+IPOS1+IPOS2)
C
C     Check for the same geological block and for a free space:
      IF(IA00.GT.0) THEN
        IF(LICB) THEN
          IF(IRAM(ICB0+IA00).NE.ICBREF) THEN
            IA00=0
            GO TO 11
          END IF
        END IF
        IF(RAM(IP0+IA00).GE.GIANT) THEN
          IA00=0
        END IF
      END IF
   11 CONTINUE
      IF(IA20.GT.0) THEN
        IF(LICB) THEN
          IF(IRAM(ICB0+IA20).NE.ICBREF) THEN
            IA20=0
            GO TO 12
          END IF
        END IF
        IF(RAM(IP0+IA20).GE.GIANT) THEN
          IA20=0
        END IF
      END IF
   12 CONTINUE
      IF(IA02.GT.0) THEN
        IF(LICB) THEN
          IF(IRAM(ICB0+IA02).NE.ICBREF) THEN
            IA02=0
            GO TO 13
          END IF
        END IF
        IF(RAM(IP0+IA02).GE.GIANT) THEN
          IA02=0
        END IF
      END IF
   13 CONTINUE
      IF(IA22.GT.0) THEN
        IF(LICB) THEN
          IF(IRAM(ICB0+IA22).NE.ICBREF) THEN
            IA22=0
            GO TO 14
          END IF
        END IF
        IF(RAM(IP0+IA22).GE.GIANT) THEN
          IA22=0
        END IF
      END IF
   14 CONTINUE
C
C     Calculating mixed second partial slowness derivatives:
      IF(IA22.GT.0) THEN
        IF(IA02.GT.0) THEN
          IF(IA20.GT.0) THEN
            IF(IA00.GT.0) THEN
              DERMIX=(RAM(IP0+IA00)-RAM(IP0+IA20)
     *               -RAM(IP0+IA02)+RAM(IP0+IA22))/4.
            ELSE
              IA00=INDX(IPOS0-IPOS1)
              IA20=INDX(IPOS0+IPOS1)
              DERMIX=(RAM(IP0+IA00)-RAM(IP0+IA20)
     *               -RAM(IP0+IA02)+RAM(IP0+IA22))/2.
            END IF
          ELSE
            IF(IA00.GT.0) THEN
              IA20=INDX(IPOS0-IPOS2)
              IA22=INDX(IPOS0+IPOS2)
              DERMIX=(RAM(IP0+IA00)-RAM(IP0+IA20)
     *               -RAM(IP0+IA02)+RAM(IP0+IA22))/2.
            ELSE
              IA00=INDX(IPOS0-IPOS1)
              IA20=INDX(IPOS0+IPOS1)
              DERMIX=(RAM(IP0+IA00)-RAM(IP0+IA20)
     *               -RAM(IP0+IA02)+RAM(IP0+IA22))/2.
            END IF
          END IF
        ELSE
          IF(IA20.GT.0) THEN
            IF(IA00.GT.0) THEN
              IA02=INDX(IPOS0-IPOS1)
              IA22=INDX(IPOS0+IPOS1)
              DERMIX=(RAM(IP0+IA00)-RAM(IP0+IA20)
     *               -RAM(IP0+IA02)+RAM(IP0+IA22))/2.
            ELSE
              IA00=INDX(IPOS0-IPOS2)
              IA02=INDX(IPOS0+IPOS2)
              DERMIX=(RAM(IP0+IA00)-RAM(IP0+IA20)
     *               -RAM(IP0+IA02)+RAM(IP0+IA22))/2.
            END IF
          ELSE
            IF(IA00.GT.0) THEN
              IA02=INDX(IPOS0-IPOS1)
              IA20=INDX(IPOS0-IPOS2)
              IA22=INDX(IPOS0)
              DERMIX= RAM(IP0+IA00)-RAM(IP0+IA20)
     *               -RAM(IP0+IA02)+RAM(IP0+IA22)
            ELSE
              IA00=INDX(IPOS0)
              IA20=INDX(IPOS0+IPOS1)
              IA02=INDX(IPOS0+IPOS2)
              DERMIX= RAM(IP0+IA00)-RAM(IP0+IA20)
     *               -RAM(IP0+IA02)+RAM(IP0+IA22)
            END IF
          END IF
        END IF
      ELSE
        IF(IA02.GT.0) THEN
          IF(IA20.GT.0) THEN
            IF(IA00.GT.0) THEN
              IA20=INDX(IPOS0-IPOS2)
              IA22=INDX(IPOS0+IPOS2)
              DERMIX=(RAM(IP0+IA00)-RAM(IP0+IA20)
     *               -RAM(IP0+IA02)+RAM(IP0+IA22))/2.
            ELSE
              IA00=INDX(IPOS0-IPOS1)
              IA20=INDX(IPOS0)
              IA22=INDX(IPOS0+IPOS2)
              DERMIX= RAM(IP0+IA00)-RAM(IP0+IA20)
     *               -RAM(IP0+IA02)+RAM(IP0+IA22)
            END IF
          ELSE
            IF(IA00.GT.0) THEN
              IA20=INDX(IPOS0-IPOS2)
              IA22=INDX(IPOS0+IPOS2)
              DERMIX=(RAM(IP0+IA00)-RAM(IP0+IA20)
     *               -RAM(IP0+IA02)+RAM(IP0+IA22))/2.
            ELSE
              IA00=INDX(IPOS0-IPOS1)
              IA20=INDX(IPOS0)
              IA22=INDX(IPOS0+IPOS2)
              DERMIX= RAM(IP0+IA00)-RAM(IP0+IA20)
     *               -RAM(IP0+IA02)+RAM(IP0+IA22)
            END IF
          END IF
        ELSE
          IF(IA20.GT.0) THEN
            IF(IA00.GT.0) THEN
              IA02=INDX(IPOS0-IPOS1)
              IA22=INDX(IPOS0+IPOS1)
              DERMIX=(RAM(IP0+IA00)-RAM(IP0+IA20)
     *               -RAM(IP0+IA02)+RAM(IP0+IA22))/2.
            ELSE
              IA00=INDX(IPOS0-IPOS2)
              IA02=INDX(IPOS0)
              IA22=INDX(IPOS0+IPOS1)
              DERMIX= RAM(IP0+IA00)-RAM(IP0+IA20)
     *               -RAM(IP0+IA02)+RAM(IP0+IA22)
            END IF
          ELSE
            IF(IA00.GT.0) THEN
              IA02=INDX(IPOS0-IPOS1)
              IA20=INDX(IPOS0-IPOS2)
              IA22=INDX(IPOS0)
              DERMIX= RAM(IP0+IA00)-RAM(IP0+IA20)
     *               -RAM(IP0+IA02)+RAM(IP0+IA22)
            ELSE
              IERR=IERR+1
            END IF
          END IF
        END IF
      END IF
C
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE SLOW(X1,X2,X3,DPOS1,DPOS2,DPOS3,
     *                         IPOS1,IPOS2,IPOS3,IPOS,IADR,PS)
      REAL X1,X2,X3,DPOS1,DPOS2,DPOS3,PS
      INTEGER IPOS1,IPOS2,IPOS3,IPOS,IADR
C
C Subroutine called by SOURCE, RECVRS, TRACER, and ONERAY, to find
C a close gridpoint within the computational volume, and to interpolate
C the slowness within the same geological block.
C If the distance of the given point from the model volume is greater
C than a grid step, an error message is generated.
C If the point is outside the computational volume or in a free space,
C another error message is generated.
C
C Input:
C     X1,X2,X3... Coordinates of the point.
C
C Output:
C     DPOS1,DPOS2,DPOS3... Vectorial distance from the nearest
C             gridpoint situated in the same geological block.
C     IPOS1,IPOS2,IPOS3... Positional indices of the nearest gridpoint.
C     IPOS... Index of the nearest gridpoint in the index array.
C     IADR... Storage address of the nearest gridpoint if it is located
C             in the computational volume, otherwise zero.
C     PS...   Interpolated slowness.
C
C-----------------------------------------------------------------------
C
C     Common block /GRID/ is required here:
      INCLUDE 'net.inc'
C     net.inc
C
C-----------------------------------------------------------------------
C
      EXTERNAL INDX
      INTEGER  INDX
      REAL GIANT
      PARAMETER (GIANT=1.E+10)
C
      INTEGER I1,I2,I3,I1D,I2D,I3D,IAUX,ICBREF
      REAL DS,AUX
C
C     I1,I2,I3... Loop variables - position of a gridpoint.
C     I1D,I2D,I3D... Loop increments.
C     IAUX... Index of the node at a gridpoint.
C     DS...   Distance of the point from a gridpoint.
C     AUX...  Temporary variable.
C     ICBREF... Index of the geological block.
C
C.......................................................................
C
C     Determination of the nearest gridpoint:
      CALL POS(X1,X2,X3,DPOS1,DPOS2,DPOS3,IPOS1,IPOS2,IPOS3,IPOS,IADR)
      I1D=INT(SIGN(1.5,DSX1*DPOS1))
      I2D=INT(SIGN(1.5,DSX2*DPOS2))
      I3D=INT(SIGN(1.5,DSX3*DPOS3))
C
C     Slowness in the given point is determined by an interpolation
C     from neighbouring grid-point slownesses
      PS=0.0
      DS=0.0
      DO 23 I3=IPOS3,IPOS3+I3D,I3D
        DO 22 I2=IPOS2,IPOS2+I2D,I2D
          DO 21 I1=IPOS1,IPOS1+I1D,I1D
            IF(0.LE.I1.AND.I1.LE.NL1-1.AND.
     *         0.LE.I2.AND.I2.LE.NL2-1.AND.
     *         0.LE.I3.AND.I3.LE.NL3-1) THEN
              IAUX=1+I1+(I2+I3*NL2)*NL1
              IAUX=INDX(IAUX)
              IF(IAUX.GT.0) THEN
C
C               Test for nodes in a free space:
                IF(RAM(IP0+IAUX).GE.GIANT) THEN
                  GO TO 20
                ENDIF
C
C               Check for the nearest gridpoint in the indexed volume
                IF(DS.EQ.0.) THEN
C
C                 Redefinition of the nearest gridpoint
                  IF(IADR.NE.IAUX) THEN
                    IF(LN1) THEN
                      DPOS1=DPOS1+DSX1*FLOAT(IPOS2-I2)
                    ELSE
                      DPOS1=DPOS1+DSX1*FLOAT(IPOS1-I1)
                    END IF
                    DPOS2=DPOS2+DSX2*FLOAT(IPOS2-I2)
                    DPOS3=DPOS3+DSX3*FLOAT(IPOS3-I3)
                    IPOS1=I1
                    IPOS2=I2
                    IPOS3=I3
                    IPOS=1+I1+(I2+I3*NL2)*NL1
                    IADR=IAUX
                  END IF
C
C                 Setting reference geological block
                  IF(LICB) THEN
                    ICBREF=IRAM(ICB0+IADR)
                  END IF
                ELSE
C
C                 Check for the (geological) block:
                  IF(LICB) THEN
                    IF(IRAM(ICB0+IAUX).NE.ICBREF) THEN
                      GO TO 20
                    ENDIF
                  ENDIF
                ENDIF
C
C               Interpolation:
                IF(LN1) THEN
                  AUX=SQRT((FLOAT(I2-IPOS2)*DSX1-DPOS1)**2+
     *                     (FLOAT(I2-IPOS2)*DSX2-DPOS2)**2+
     *                     (FLOAT(I3-IPOS3)*DSX3-DPOS3)**2)
                ELSE
                  AUX=SQRT((FLOAT(I1-IPOS1)*DSX1-DPOS1)**2+
     *                     (FLOAT(I2-IPOS2)*DSX2-DPOS2)**2+
     *                     (FLOAT(I3-IPOS3)*DSX3-DPOS3)**2)
                END IF
                IF(AUX.EQ.0.) THEN
                  PS=RAM(IP0+IAUX)
                  RETURN
                END IF
                PS=PS+RAM(IP0+IAUX)/AUX
                DS=DS+1.0/AUX
C
              ENDIF
            ENDIF
   20       CONTINUE
   21     CONTINUE
   22   CONTINUE
   23 CONTINUE
      IF(DS.EQ.0.) THEN
C       NET-45
        CALL ERROR('NET-45: Source or receiver point in a free space')
C       NET-45: Source or receiver point in a free space:
C               Source or receiver point is situated outside the indexed
C               computational volume or in a free space (zero velocity).
      ENDIF
      PS=PS/DS
C
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE POS(X1,X2,X3,DPOS1,DPOS2,DPOS3,
     *                                      IPOS1,IPOS2,IPOS3,IPOS,IADR)
      REAL X1,X2,X3,DPOS1,DPOS2,DPOS3
      INTEGER IPOS1,IPOS2,IPOS3,IPOS,IADR
C
C Auxiliary subroutine to SLOW, to find the nearest gridpoint.
C This subroutine checks the position of the given point with respect to
C the boundaries of the model volume.  If the distance from the model
C volume is greater than a grid step, the error message is generated.
C
C Input:
C     X1,X2,X3... Coordinates of the point.
C
C Output:
C     DPOS1,DPOS2,DPOS3... Vectorial distance from the nearest
C             gridpoint.
C     IPOS1,IPOS2,IPOS3... Positional indices of the nearest gridpoint.
C     IPOS... Index of the nearest gridpoint in the index array.
C     IADR... Storage address of the nearest gridpoint if it is located
C             in the computational volume, otherwise zero.
C
C-----------------------------------------------------------------------
C
C     Common block /GRID/ is required here:
      INCLUDE 'net.inc'
C     net.inc
C
C-----------------------------------------------------------------------
C
      EXTERNAL INDX
      INTEGER  INDX
C
C.......................................................................
C
      IF(LN1) THEN
        IPOS1=0
        IPOS2=INT( ((X1-X1MIN)*DSX1+(X2-X2MIN)*DSX2)/(DSX1**2+DSX2**2) )
        IF(IPOS2.LT.0.OR.NL2.LT.IPOS2) THEN
C         NET-46
          CALL ERROR
     *            ('NET-46: Source or receiver point outside the model')
C         NET-46: Source or receiver point outside the model:
C                 Source or receiver point is situated outside the model
C                 volume.
        ELSE IF(IPOS2.GE.NL2) THEN
          IPOS2=NL2-1
        END IF
        DPOS1=X1-X1MIN-(FLOAT(IPOS2)+0.5)*DSX1
        DPOS2=X2-X2MIN-(FLOAT(IPOS2)+0.5)*DSX2
      ELSE
        IF(NL1.EQ.1) THEN
          IPOS1=0
          DPOS1=0.
        ELSE
          IPOS1=INT((X1-X1MIN)/DSX1)
          IF(IPOS1.LT.0.OR.NL1.LT.IPOS1) THEN
C           NET-47
            CALL ERROR
     *            ('NET-47: Source or receiver point outside the model')
C           NET-47: Source or receiver point outside the model:
C                   Source or receiver point is situated outside the
C                   model volume.
          ELSE IF(IPOS1.GE.NL1) THEN
            IPOS1=NL1-1
          END IF
          DPOS1=X1-X1MIN-(FLOAT(IPOS1)+0.5)*DSX1
        END IF
        IF(NL2.EQ.1) THEN
          IPOS2=0
          DPOS2=0.
        ELSE
          IPOS2=INT((X2-X2MIN)/DSX2)
          IF(IPOS2.LT.0.OR.NL2.LT.IPOS2) THEN
C           NET-48
            CALL ERROR
     *            ('NET-48: Source or receiver point outside the model')
C           NET-48: Source or receiver point outside the model:
C                   Source or receiver point is situated outside the
C                   model volume.
          ELSE IF(IPOS2.GE.NL2) THEN
            IPOS2=NL2-1
          END IF
          DPOS2=X2-X2MIN-(FLOAT(IPOS2)+0.5)*DSX2
        END IF
      END IF
      IF(NL3.EQ.1) THEN
        IPOS3=0
        DPOS3=0.
      ELSE
        IPOS3=INT((X3-X3MIN)/DSX3)
        IF(IPOS3.LT.0.OR.NL3.LT.IPOS3) THEN
C         NET-54
          CALL ERROR
     *            ('NET-54: Source or receiver point outside the model')
C         NET-54: Source or receiver point outside the model.
        ELSE IF(IPOS3.GE.NL3) THEN
          IPOS3=NL3-1
        END IF
        DPOS3=X3-X3MIN-(FLOAT(IPOS3)+0.5)*DSX3
      END IF
C
      IPOS =1+IPOS1+(IPOS2+IPOS3*NL2)*NL1
      IADR=INDX(IPOS)
      RETURN
      END
C
C=======================================================================
C
C     
C
      INTEGER FUNCTION INDX(IPOS)
      INTEGER IPOS
C
C Integer function, called by many subroutines, to return the index of
C the node at the given gridpoint.  The subroutine does not check the
C validity of the input argument.
C
C-----------------------------------------------------------------------
C
C     Common block /GRID/ is required here:
      INCLUDE 'net.inc'
C     net.inc
C
C-----------------------------------------------------------------------
C
      INTEGER IPOS1,IPOS2,IPOS3,I1,I2,I3,J1,J2,J3,J
C
C.......................................................................
C
      IF(L4.EQ.0) THEN
        INDX=IPOS
      ELSE
        IPOS1=IPOS-1
        IPOS2=IPOS1/NL1
        IPOS3=IPOS2/NL2
        IPOS1=IPOS1-IPOS2*NL1
        IPOS2=IPOS2-IPOS3*NL2
        J1=IPOS1/L1
        J2=IPOS2/L2
        J3=IPOS3/L3
        J=IRAM(IND0+1+J1+N1*(J2+N2*J3))
        IF(J.LE.0) THEN
          INDX=0
        ELSE
          I1=IPOS1-J1*L1
          I2=IPOS2-J2*L2
          I3=IPOS3-J3*L3
          INDX=J+I1+L1*(I2+L2*I3)
        END IF
      END IF
      RETURN
      END
C
C=======================================================================
C
      INCLUDE 'ttt.for'
C     ttt.for
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'sep.for'
c     sep.for
      INCLUDE 'length.for'
c     length.for
      INCLUDE 'forms.for'
c     forms.for
      INCLUDE 'eigen.for'
C     eigen.for
C
C=======================================================================
C