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: 2.00
C DATE: 1994, AUGUST 29
C
C CODED BY: MICHAL KVASNICKA
C     DEPARTMENT OF GEOPHYSICS, CHARLES UNIVERSITY PRAGUE,
C     KE KARLOVU 3, 121 16 PRAHA 2, CZECH REPUBLIC,
C     E-MAIL: GEOFMK@ACI.CVUT.CZ
C AND LUDEK KLIMES
C     DEPARTMENT OF GEOPHYSICS, CHARLES UNIVERSITY PRAGUE,
C     KE KARLOVU 3, 121 16 PRAHA 2, CZECH REPUBLIC,
C     E-MAIL: PSENCIK@EARN.CVUT.CZ, OR  PSENCIK@CSEARN.BITNET
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.
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
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 INEXED 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 (ADRESSES 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 BICK DEVIDED 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 EXPLICITELY 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 CONUS 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.......................................................................
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) 'NET'
C     THE STRING CONTAINING THE NAME OF THE MAIN INPUT DATA FILE.
C     ONLY THE FIRST 80 CHARACTERS OF THE STRING ARE SIGNIFICANT.
C DEFAULT: 'NET'='NET.DAT'.
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) N1,N2,N3,L1,L2,L3,/
C     N1,N2,N3... NUMBERS OF BIG BRICKS IN THE X1,X2,X3-AXIS DIRECTIONS.
C             MUST BE POSITIVE.
C             N3 NEED NOT BE SPECIFIED FOR A 2-D MODEL.
C     L1,L2,L3... NUMBERS OF SMALL BRICKS IN ONE BIG BRICK IN THE
C             X1,X2,X3-AXIS DIRECTIONS.  IF SPECIFIED, MUST BE POSITIVE.
C             L1,L2,L3 MAY BE LEFT OUT AND SHOULD ONLY BE SPECIFIED BY
C             EXPERIENCED USERS.
C     DEFAULT: N3=1, L1=1, L2=1, 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 (2) X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,/
C     X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX ... BOUNDARIES OF THE MODEL
C             VOLUME.
C             BECAUSE THE MODEL VOLUME IS DEVIDED 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             THUS, X3MIN,X3MAX NEED NOT BE SPECIFIED IF N3 AT (1) IS
C             NOT SPECIFIED FOR A 2-D MODEL.
C (3) NFSMAX,/
C     NFSMAX..MAXIMUM SIZE OF A FORWARD STAR.  SIMULTANEOUSLY, THE
C             DEFAULT SIZE OF A FORWARD STAR IF THE SIZE IS NOT
C             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 (4) /
C     LINE RESERVED FOR FUTURE EXETENSIONS LIKE FACTORIZED ANISOTROPIC
C     MEDIA.
C (5) 'SRC','REC','RAYS','END',/
C     'SRC'...NAME OF THE INPUT FILE WITH SOURCE COORDINATES.
C     'REC'...NAME OF THE INPUT FILE WITH RECEIVER COORDINATES.
C             IF BLANK, NO RAYS ARE STORED WITHIN THE FILE 'RAYS'.
C     'RAYS'..NAME OF THE OUTPUT FILE WITH RAYS.
C             IF BLANK, NO RAYS ARE STORED WITHIN THE FILE '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     DEFAULT: 'REC'=' ', 'RAYS'=' ', 'END'=' '.
C (6) 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 (7.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 (7) ONCE (7.1), THEN NREFL-TIMES (7.2) AND (7.1):
C (7.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     '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     '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     '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     '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     '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 LEFT OUT FOR MOST APPLICATIONS.
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 (3).
C     DEFAULT: 'IND(I)'=' ', 'VEL(I)'=' ', 'ICB(I)'=' ',
C             'TT(I)'=' ','ERR(I)'=' ', 'PRED(I)'=' ', 'NFS(I)'=' '.
C (7.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
C INPUT FILE 'SRC' CONTAINING SOURCE COORDINATES:
C (1) SEVERAL STRINGS TERMINATED BY / (A SLASH).
C             THE SIMLEST 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 (3) / (A SLASH).
C
C INPUT FILE 'REC' CONTAINING RECEIVER COORDINATES:
C (1) SEVERAL STRINGS TERMINATED BY / (A SLASH).
C             THE SIMLEST 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
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 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 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 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
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 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 OUTPUT FILES 'ERR(I)' (TIME ERROR FIELD):
C (1) (ERR(I),I=1,L1234)
C     ERR(I)..ESTIMATED MAXIMUM TRAVEL TIME ERROR 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 ERRORS IN A FREE SPACE.
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 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
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     INPUT...READS THE WHOLE MAIN INPUT DATA FILE 'NET', SOURCE AND
C             RECEIVER COORDINATES.
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     GENER...PERFORMS SHORTEST PATH RAY TRACING (ONE ITERATION).
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     ERRORS..GENERATES THE FIELD OF ARRIVAL-TIME ERRORS BY MEANS OF
C             BACKWARD RAY TRACING - WRITES 'ERR(I)' (ONE ITERATION).
C     RECVRS..UPDATES ARRIVAL TIMES AT THE RECEIVERS.
C     TRACER..PERFORMS BACKWARD RAY TRACING FROM THE RECEIVERS,
C             INCLUDING THE ESTIMATION OF ARRIVAL-TIME ERRORS.
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 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 TEMPLATE FORWARD STARS:
C     READFS..CALLED BY SOURCE, TO READ AND CONSTRUCT OPTIMIZED TEMPLATE
C             FORWARD STARS.
C     MAKEFS..CALLED BY SOURCE AND RECVRS, TO CONSTRUCT FULL SPHERICAL
C             TEMPLATE FORWARD STARS.
C     SETFS...AUXILIARY SUBROUTINE TO READFS AND MAKEFS.
C     STORFS..AUXILIARY SUBROUTINE TO SETFS.
C BACKWARD RAY TRACING:
C     ONERAY..SUBROUTINE CALLED BY ERRORS AND TRACER, TO PERFORM
C             BACKWARD RAY TRACING.
C ACCURACY ESTIMATION:
C     SETERR..SUBROUTINE CALLED BY ONERAY, TO ESTIMATE AND ACCUMULATE
C             THE ARRIVAL-TIME ERRORS DURING BACKWARD RAY TRACING.
C     OPTNFS..SUBROUTINE CALLED BY SOURCE AND SETERR, TO SET UP
C             OPTIMUM SIZES OF FORWARD STARS.
C     OPTMAT..AUXILIARY SUBROUTINE TO SETERR AND OPTNFS, TO CALCULATE
C             THE MATRIX COMPOSED OF THE FIRST AND SECOND SLOWNESS
C             DERIVATIVES.
C     TRYMAT..AUXILIARY SUBROUTINE TO OPTMAT, TO TRY THE CALCULATION OF
C             THE MATRIX.
C     MIXDER..AUXILIARY SUBROUTINE TO TRYMAT, TO CALCULATE MIXED SECOND
C             PARTIAL SLOWNESS DERIVATIVES.
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     POS...  AUXILIARY SUBROUTINE TO SLOW, TO FIND THE NEAREST
C             GRIDPONT.
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
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 'ARRAY.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 ERROR MESSAGES:
C     ERROR 10: NUMBER OF GRIDPOINTS IS NOT POSITIVE:
C               N1,N2,N3, AND, IF SPECIFIED, ALSO L1,L2,L3, IN THE MAIN
C               INPUT FILE 'NET', (1), MUST BE POSITIVE.
C     ERROR 11: MAXIMUM SIZE OF A FORWARD STAR TOO GREAT:
C               NFSMAX IN THE MAIN INPUT DATA FILE 'NET', (3), SHOULD BE
C               DECREASED,
C               OR THE DIMENSION MFSMAX IN THE COMMON BLOCK /FS/ SHOULD
C               BE INCREASED.
C     ERROR 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', (5).
C     ERROR 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', (5).
C     ERROR 14: TOO MANY REFLECTIONS:
C               NUMBER NREFL (MAIN INPUT FILE 'NET', (5)) SHOULD BE
C               DECREASED TO SATISFY THE INEQUALITY NREFL.LE.MREFL,
C               OR THE DIMENSION MREFL IN THE COMMON BLOCK /FILES/
C               SHOULD BE INCREASED.
C     ERROR 15: REFLECTOR FILE NOT SPECIFIED:
C               IF YOU WANT TO COMPUTE REFLECTED TIME FIELD, YOU
C               MUST SPECIFY NAME OF REFLECTOR FILE.
C               SEE THE MAIN INPUT FILE 'NET', (7.2).
C     ERROR 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 THE MAIN INPUT FILE 'NET', (7.1).
C     ERROR 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.
C     ERROR 18: DIFFERENT 1-ST SOURCE COORDINATES:
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     ERROR 19: DIFFERENT 2-ND SOURCE COORDINATES:
C               SEE ERROR 18.
C     ERROR 20: DIFFERENT 3-RD SOURCE COORDINATES:
C               SEE ERROR 18.
C     ERROR 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), AND INPUT FILE
C               'SRC'.
C     ERROR 22: TOO MANY SOURCE POINTS:
C               NUMBER OF RECEIVERS IN THE INPUT FILE 'REC' SHOULD BE
C               LESS THAN MREC,
C               OR THE DIMENSION MREC IN THE COMMON BLOCK /POINTS/
C               SHOULD BE INCREASED.
C     ERROR 23: NO SOURCE POINT GIVEN:
C               YOU MUST SPECIFY AT LEAST ONE SOURCE POINT AT THE
C               INPUT SOURCE FILE 'SRC'.
C     ERROR 24: DIFFERENT 1-ST SOURCE AND RECEIVER COORDINATE:
C               SEE ERROR 18.
C     ERROR 25: DIFFERENT 2-ND SOURCE AND RECEIVER COORDINATE:
C               SEE ERROR 18.
C     ERROR 26: DIFFERENT 3-RD SOURCE AND RECEIVER COORDINATE:
C               SEE ERROR 18.
C     ERROR 27: RECEIVER IS NOT INSIDE THE MODEL:
C               RECEIVER COORDINATES X1R,X2R,X3R MUST SATISFY 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), AND INPUT FILE
C               'REC'.
C     ERROR 28: TOO MANY RECEIVERS:
C               NUMBER OF RECEIVERS IN THE INPUT FILE 'REC' SHOULD BE
C               LESS THAN MREC,
C               OR THE DIMENSION MREC IN THE COMMON BLOCK /POINTS/
C               SHOULD BE INCREASED.
C     ERROR 29: NO INDEX FILE SPECIFIED:
C               IF L1,L2,L3 (MAIN INPUT DATA FILE 'NET', LINE (1))
C               ARE SPECIFIED THE INDEX FILE 'IND' HAS TO BE SUBMITTED.
C               SEE THE MAIN INPUT FILE 'NET', (7.1).
C     ERROR 30: TOO MANY BIG BRICKS:
C               NUMBERS N1,N2,N3 (MAIN INPUT DATA 'NET', (1)) SHOULD BE
C               DECREASED TO SATISFY THE INEQUALITY N1*N2*N3.LE.MIND, OR
C               THE DIMENSION MIND OF THE ARRAY IND IN THE COMMON BLOCK
C               /GRID/ SHOULD BE INCREASED.
C     ERROR 31: TOO MANY NETWORK NODES:
C               THE NUMBER OF NETWORK NODES, SEE THE MAIN INPUT FILE
C               'NET', (1), SHOULD NOT EXCEED MGRID,
C               OR THE DIMENSION MGRID IN THE COMMON BLOCK /GRID/ SHOULD
C               BE INCREASED.
C     ERROR 32: INSUFFICIENT MEMORY FOR INDICES OF BLOCKS:
C               THE NUMBER OF NETWORK NODES, SEE THE MAIN INPUT FILE
C               'NET', (1), SHOULD NOT EXCEED MICB,
C               OR THE DIMENSION MICB IN THE COMMON BLOCK /GRID/ SHOULD
C               BE INCREASED.
C     ERROR 33: INSUFFICIENT MEMORY FOR PREDECESSORS:
C               THE NUMBER OF NETWORK NODES, SEE THE MAIN INPUT FILE
C               'NET', (1), SHOULD NOT EXCEED MPRED,
C               OR THE DIMENSION MPRED IN THE COMMON BLOCK /GRID/ SHOULD
C               BE INCREASED.
C     ERROR 34: INSUFFICIENT MEMORY FOR FORWARD STAR SIZES:
C               THE NUMBER OF NETWORK NODES, SEE THE MAIN INPUT FILE
C               'NET', (1), SHOULD NOT EXCEED MNFS,
C               OR THE DIMENSION MNFS IN THE COMMON BLOCK /GRID/ SHOULD
C               BE INCREASED.
C     ERROR 35: ONE-POINT MODEL:
C               N1*L1=N2*L2=N3*L3=1 IN THE MAIN INPUT FILE 'NET', (1).
C     ERROR 36: LINE MODEL:
C               N2*L2=N3*L3=1 IN THE MAIN INPUT FILE 'NET', (1).
C     ERROR 37: LINE MODEL:
C               N1*L1=N3*L3=1 IN THE MAIN INPUT FILE 'NET', (1).
C     ERROR 38: LINE MODEL:
C               N1*L1=N2*L2=1 IN THE MAIN INPUT FILE 'NET', (1).
C     ERROR 39: TOO MANY POINTS OF REFLECTOR:
C               NUMBER OF REFLECTOR POINTS SHOULD BE LESS THAN MGRID,
C               OR THE DIMENSION MGRID IN THE COMMON BLOCK /GRID/ SHOULD
C               BE INCREASED.
C     ERROR 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 THE
C               MAIN INPUT DATA FILE 'NET', (3), SHOULD BE DECREASED.
C     ERROR 41: TOO MANY INPUT NODES OF A FORWARD STAR:
C               NFSMAX IN THE MAIN INPUT DATA FILE 'NET', (3), SHOULD BE
C               DECREASED,
C               OR THE DIMENSION MD IN THE SUBROUTINE READFS SHOULD BE
C               INCREASED.
C     ERROR 42: TOO BIG FORWARD STAR:
C               NFSMAX IN THE MAIN INPUT DATA FILE 'NET', (3), SHOULD BE
C               DECREASED,
C               OR THE DIMENSION MFS IN THE COMMON BLOCK /FS/ SHOULD
C               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) REMOVE ASTERISKS FROM THE FIRST COLUMN OF SECTION
C                     'PRINTING THE TABLE...' AT THE END OF THE
C                     SUBROUTINE READFS.
C                 (C) COMPILE AND RUN THE PROGRAM.
C                 (D) LOOK AT THE FILE 'NET.TAB' 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) RESTORE THE REMOVED ASTERISKS.
C     ERROR 43: TIME ERRORS CANNOT BE COMPUTED AFTER REFLECTION:
C               IF THE VELOCITY FIELD OF THE REFLECTED WAVE IS DIFFERENT
C               FROM THE VELOCITY FIELD BEFORE THE REFLECTION, THE
C               ERRORS OF THE ARRIVAL TIMES CANNOT BE COMPUTED CORRECTLY
C               BY THIS VERSION.  WHEN THE RAYS ARE TRACED, ONLY TIME
C               FIELDS AND PREDECESSORS ARE ALTERNATED IN THE MOMORY,
C               UNLIKE SLOWNESS FIELDS, SEE SUBROUTINE SETERR.
C               THUS, IF THE FILENAMES 'VEL(I)' AT THE LINES (7.1) OF
C               THE MAIN INPUT DATA 'NET' DIFFER, THE FILENAMES 'ERR(I)'
C               SHOULD NOT BE SPECIFIED AT THE LINES (7.1) FOR THE
C               REFLECTED WAVES.  SIMULTANEOUSLY, IF THE FILENAME 'END'
C               AT THE LINE (5) IS SPECIFIED, ERRORS WRITTEN IN THE
C               OUTPUT FILE 'END' ARE INCORRECT.
C     WARNING 44: TIME ERRORS CANNOT BE COMPUTED AFTER REFLECTION:
C               IF, IN A CASE OF REFLECTION, FILENAMES 'VEL(I)' AT THE
C               LINES (7.1) OF THE MAIN INPUT DATA 'NET' DIFFER AND THE
C               FILENAME 'END' AT THE LINE (5) IS SPECIFIED, ERRORS
C               WRITTEN IN THE OUTPUT FILE 'END' ARE INCORRECT.
C               FOR MORE DETAILS SEE ALSO ERROR 43.
C     ERROR 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).
C     ERROR 46: SOURCE OR RECEIVER POINT OUTSIDE THE MODEL:
C               SOURCE OR RECEIVER POINT IS SITUATED OUTSIDE THE MODEL
C               VOLUME.
C     ERROR 47: NO SOURCE POINT:
C               THIS ERROR SHOULD NOT APPEAR.  CONTACT THE AUTHORS.
C     ERROR 48: NFSMAX NOT POSITIVE:
C               MAXIMUM FORWARD-STAR SIZE NFSMAX IN THE MAIN INPUT DATA
C               FILE MUST BE POSITIVE.
C     ERROR 49: INSUFFICIENT MEMORY FOR P2:
C               REFER TO ERROR 48.
C     ERROR 50: INSUFFICIENT MEMORY FOR P3:
C               REFER TO ERROR 48.
C
C-----------------------------------------------------------------------
C
      PROGRAM NET
C
C-----------------------------------------------------------------------
C
C COMMON BLOCKS:
C
      INTEGER MIND,MGRID,MICB,MNFS,MPRED
      PARAMETER (MIND=1,MGRID=100000)
      PARAMETER (MICB=MGRID,MNFS=MGRID,MPRED=MGRID)
      LOGICAL LN1,LICB,LNFS,LPRED
      INTEGER N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,NREFL
      INTEGER IND(MIND),ICB(MICB),NFS(MNFS),IPRED(MPRED)
      INTEGER IPOSQ(MGRID),MINQ,MAXQ
      REAL X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3
      REAL P(MGRID),TT(MGRID),TTMIN,DMIN,ASX1,ASX2,ASX3
      COMMON/GRID/ N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,
     *             X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3,
     *             NREFL,LN1,LICB,LNFS,LPRED,IND,P,ICB,NFS,TT,IPRED,
     *             IPOSQ,MINQ,MAXQ,TTMIN,DMIN,ASX1,ASX2,ASX3
C
      INTEGER MFSMAX,MFS
      PARAMETER (MFSMAX=100,MFS=6626)
      INTEGER NFSMAX,KFS0(0:MFSMAX)
      INTEGER KFS1(MFS),KFS2(MFS),KFS3(MFS),KFS4(MFS),KFS5(MFS)
      REAL DFS(MFS)
      COMMON/FS/ NFSMAX,KFS0,KFS1,KFS2,KFS3,KFS4,KFS5,DFS
C
      INTEGER MSRC,MREC
      PARAMETER (MSRC=100,MREC=100)
      CHARACTER*6 NAMES,NAMER(MREC)
      INTEGER NSRC,NREC,IPREDR(MREC)
      REAL X1S(MSRC),X2S(MSRC),X3S(MSRC),TTS(MSRC),TTSERR(MSRC)
      REAL X1R(MREC),X2R(MREC),X3R(MREC),TTR(MREC)
      COMMON/NAMESR/ NAMES,NAMER
      COMMON/POINTS/ NSRC,X1S,X2S,X3S,TTS,TTSERR,NREC,X1R,X2R,X3R,TTR,
     *               IPREDR
C
      INTEGER MREFL
      PARAMETER (MREFL=2)
      CHARACTER*80 FSRC,FREC,FRAYS,FEND
      CHARACTER*80 FIND(0:MREFL),FVEL(0:MREFL),FICB(0:MREFL)
      CHARACTER*80 FTT(0:MREFL),FERR(0:MREFL),FPRED(0:MREFL)
      CHARACTER*80 FNFS(0:MREFL),FINTF(1:MREFL)
      COMMON/FILES/ FSRC,FREC,FRAYS,FEND,
     *              FIND,FVEL,FICB,FTT,FERR,FPRED,FINTF,FNFS
C
C     NREFL.. NUMBER OF REFLECTIONS.
C     P...    ARRAY OF THE SLOWNESSES AT NETWORK NODES.
C     ASX1,ASX2,ASX3... ABS(DSX1),ABS(DSX2),ABS(DSX3), RESPECTIVELY.
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
      SUBROUTINE INPUT()
C
C PROCEDURE TO READ THE WHOLE MAIN INPUT DATA FILE 'NET', SOURCE AND
C RECEIVER COORDINATES.
C
C-----------------------------------------------------------------------
C
      INTEGER MIND,MGRID,MICB,MNFS,MPRED
      PARAMETER (MIND=1,MGRID=100000)
      PARAMETER (MICB=MGRID,MNFS=MGRID,MPRED=MGRID)
      LOGICAL LN1,LICB,LNFS,LPRED
      INTEGER N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,NREFL
      INTEGER IND(MIND),ICB(MICB),NFS(MNFS),IPRED(MPRED)
      INTEGER IPOSQ(MGRID),MINQ,MAXQ
      REAL X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3
      REAL P(MGRID),TT(MGRID),TTMIN,DMIN,ASX1,ASX2,ASX3
      COMMON/GRID/ N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,
     *             X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3,
     *             NREFL,LN1,LICB,LNFS,LPRED,IND,P,ICB,NFS,TT,IPRED,
     *             IPOSQ,MINQ,MAXQ,TTMIN,DMIN,ASX1,ASX2,ASX3
C
      INTEGER MFSMAX,MFS
      PARAMETER (MFSMAX=100,MFS=6626)
      INTEGER NFSMAX,KFS0(0:MFSMAX)
      INTEGER KFS1(MFS),KFS2(MFS),KFS3(MFS),KFS4(MFS),KFS5(MFS)
      REAL DFS(MFS)
      COMMON/FS/ NFSMAX,KFS0,KFS1,KFS2,KFS3,KFS4,KFS5,DFS
C
      INTEGER MSRC,MREC
      PARAMETER (MSRC=100,MREC=100)
      CHARACTER*6 NAMES,NAMER(MREC)
      INTEGER NSRC,NREC,IPREDR(MREC)
      REAL X1S(MSRC),X2S(MSRC),X3S(MSRC),TTS(MSRC),TTSERR(MSRC)
      REAL X1R(MREC),X2R(MREC),X3R(MREC),TTR(MREC)
      COMMON/NAMESR/ NAMES,NAMER
      COMMON/POINTS/ NSRC,X1S,X2S,X3S,TTS,TTSERR,NREC,X1R,X2R,X3R,TTR,
     *               IPREDR
C
      INTEGER MREFL
      PARAMETER (MREFL=2)
      CHARACTER*80 FSRC,FREC,FRAYS,FEND
      CHARACTER*80 FIND(0:MREFL),FVEL(0:MREFL),FICB(0:MREFL)
      CHARACTER*80 FTT(0:MREFL),FERR(0:MREFL),FPRED(0:MREFL)
      CHARACTER*80 FNFS(0:MREFL),FINTF(1:MREFL)
      COMMON/FILES/ FSRC,FREC,FRAYS,FEND,
     *              FIND,FVEL,FICB,FTT,FERR,FPRED,FINTF,FNFS
C
C-----------------------------------------------------------------------
C
      INTEGER LU1
      REAL DWARF
      PARAMETER(LU1=1)
      PARAMETER (DWARF=1.E-10)
C
      CHARACTER*80 TEXT
      INTEGER ISRC,IREC,I,L
      REAL A1111,A1122,A2222,A1133,A2233,A3333,A1123
      REAL A2223,A3323,A2323,A1113,A2213,A3313,A2313
      REAL A1313,A1112,A2212,A3312,A2312,A1312,A1212
C
C     TEXT...NAME OF THE MAIN INPUT FILE 'NET', 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:
      WRITE(*,'(A)')
     * ' ENTER NAME OF THE MAIN INPUT DATA FILE (DEFAULT ''NET.DAT''): '
      TEXT='NET.DAT'
      READ(*,*) TEXT
C     END OF READING FROM THE * UNIT.
      OPEN(LU1,FILE=TEXT,STATUS='OLD')
      WRITE(*,'(2A)')
     *  '+READING THE MAIN INPUT DATA FILE: ',TEXT(1:27)
C
C     (1) NUMBERS OF GRIDPOINTS:
      N3=1
      L1=1
      L2=1
      L3=1
      READ(LU1,*) N1,N2,N3,L1,L2,L3
      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
        PAUSE 'ERROR 10: NUMBER OF GRIDPOINTS IS NOT POSITIVE'
        STOP
      END IF
      NL1=N1*L1
      NL2=N2*L2
      NL3=N3*L3
C
C     (2) BOUNDARIES OF MODEL VOLUME:
      READ(LU1,*) X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX
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     (3) SIZE OF THE FORWARD STAR:
      READ(LU1,*) NFSMAX
      IF(NFSMAX.GT.MFSMAX) THEN
        PAUSE 'ERROR 11: MAXIMUM SIZE OF A FORWARD STAR TOO GREAT'
        STOP
      END IF
C
C     (4) FACTORIZED ANISOTROPY:
      A1111=1.
      A1122=0.
      A2222=1.
      A1133=0.
      A2233=0.
      A3333=1.
      A1123=0.
      A2223=0.
      A3323=0.
      A2323=1.
      A1113=0.
      A2213=0.
      A3313=0.
      A2313=0.
      A1313=1.
      A1112=0.
      A2212=0.
      A3312=0.
      A2312=0.
      A1312=0.
      A1212=1.
C                 A 1 1,A 1 2,A 2 2,A 1 3,A 2 3,A 3 3,A 1 4,
C                 A 2 4,A 3 4,A 4 4,A 1 5,A 2 5,A 3 5,A 4 5,
C                 A 5 5,A 1 6,A 2 6,A 3 6,A 4 6,A 5 6,A 6 6
      READ(LU1,*) A1111,A1122,A2222,A1133,A2233,A3333,A1123,
     *            A2223,A3323,A2323,A1113,A2213,A3313,A2313,
     *            A1313,A1112,A2212,A3312,A2312,A1312,A1212
      A1111     = A1111*A1122*A2222*A1133*A2233*A3333*A1123*
     *            A2223*A3323*A2323*A1113*A2213*A3313*A2313*
     *            A1313*A1112*A2212*A3312*A2312*A1312*A1212
C
C     (5) 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
        PAUSE'ERROR 12: NAME OF INPUT FILE WITH RECEIVERS NOT SPECIFIED'
        STOP
      END IF
      IF(FREC.NE.' '.AND.FRAYS.EQ.' '.AND.FEND.EQ.' ') THEN
        PAUSE 'ERROR 13: NAME OF OUTPUT FILE WITH RAYS NOT SPECIFIED'
        STOP
      END IF
C
C     (6) NUMBER OF REFLECTIONS:
      NREFL=0
      READ(LU1,*) NREFL
      IF(NREFL.GT.MREFL) THEN
        PAUSE 'ERROR 14: TOO MANY REFLECTIONS'
        STOP
      END IF
C
C     (7) 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
            PAUSE 'ERROR 15: REFLECTOR FILE NOT SPECIFIED'
            STOP
          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
          PAUSE 'ERROR 16: VELOCITY FILE NOT SPECIFIED'
          STOP
        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)')
     *  '+READING SOURCE COORDINATES                                   '
      OPEN(LU1,FILE=FSRC,STATUS='OLD')
      READ(LU1,*) (TEXT,I=1,20)
      DO 48 ISRC=1,MSRC
        X1S(ISRC)=-999999.
        X2S(ISRC)=-999999.
        X3S(ISRC)=-999999.
        TTS(ISRC)=0.
        TTSERR(ISRC)=0.
        READ(LU1,*)
     *        NAMES,X1S(ISRC),X2S(ISRC),X3S(ISRC),TTS(ISRC),TTSERR(ISRC)
        IF(X1S(ISRC).EQ.-999999..AND.X2S(ISRC).EQ.-999999..AND.
     *                               X3S(ISRC).EQ.-999999.) THEN
          NSRC=ISRC-1
          GO TO 49
        END IF
C       CHECKING THE INITIAL ARRIVAL TIME AT SOURCE POINTS
        IF(TTS(ISRC).LT.0.) THEN
          PAUSE 'ERROR 17: NEGATIVE TIME AT THE SOURCE NOT ACCEPTABLE'
          STOP
        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
              PAUSE 'ERROR 18: DIFFERENT 1-ST SOURCE COORDINATES'
              STOP
            END IF
          END IF
          IF(NL2.EQ.1) THEN
            IF(X2S(L).NE.X2MIN) THEN
              PAUSE 'ERROR 19: DIFFERENT 2-ND SOURCE COORDINATES'
              STOP
            END IF
          END IF
          IF(NL3.EQ.1) THEN
            IF(X3S(L).NE.X3MIN) THEN
              PAUSE 'ERROR 20: DIFFERENT 3-RD SOURCE COORDINATES'
              STOP
            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
          PAUSE 'ERROR 21: SOURCE IS NOT INSIDE THE MODEL'
          STOP
        ENDIF
   48 CONTINUE
      PAUSE 'ERROR 22: TOO MANY SOURCE POINTS'
      STOP
   49 CONTINUE
      CLOSE(LU1)
      IF(NSRC.LE.0) THEN
        PAUSE 'ERROR 23: NO SOURCE POINT GIVEN'
        STOP
      END IF
C
C.......................................................................
C
C     READING RECEIVER COORDINATES:
C
      WRITE(*,'(A)')
     *  '+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)=-999999.
          X2R(IREC)=-999999.
          X3R(IREC)=-999999.
          READ(LU1,*) NAMER(IREC),X1R(IREC),X2R(IREC),X3R(IREC)
          IF(X1R(IREC).EQ.-999999..AND.X2R(IREC).EQ.-999999..AND.
     *                                 X3R(IREC).EQ.-999999.) THEN
            NREC=IREC-1
            GO TO 59
          END IF
C         CHECKING THE RECEIVER POSITIONS
          IF(NL1.EQ.1.AND..NOT.LN1) THEN
            IF(X1R(IREC).NE.X1MIN) THEN
              PAUSE
     *        'ERROR 24: DIFFERENT 1-ST SOURCE AND RECEIVER COORDINATE.'
              STOP
            END IF
          END IF
          IF(NL2.EQ.1) THEN
            IF(X2R(IREC).NE.X2MIN) THEN
              PAUSE
     *        'ERROR 25: DIFFERENT 2-ND SOURCE AND RECEIVER COORDINATE.'
              STOP
            END IF
          END IF
          IF(NL3.EQ.1) THEN
            IF(X3R(IREC).NE.X3MIN) THEN
              PAUSE
     *        'ERROR 26: DIFFERENT 3-RD SOURCE AND RECEIVER COORDINATE.'
              STOP
            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
            PAUSE 'ERROR 27: RECEIVER IS NOT INSIDE THE MODEL'
            STOP
          ENDIF
   58   CONTINUE
        PAUSE 'ERROR 28: TOO MANY RECEIVERS'
        STOP
   59   CONTINUE
        CLOSE(LU1)
      END IF
C
C.......................................................................
C
C     WRITING PARAMETERS ON THE SCREEN:
      WRITE(*,'(A,I2,A,3(I4,A),I7,A,I2,A,I4,A)')
     *  '+',NREFL,' REFLECTIONS,  ',NL1,'*',NL2,'*',NL3,'=',NL1*NL2*NL3,
     *  ' GRIDPOINTS, F.S.SIZE:',NFSMAX,',',NREC,' RECEIVERS'
      WRITE(*,*)
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE SOURCE(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
      INTEGER MIND,MGRID,MICB,MNFS,MPRED
      PARAMETER (MIND=1,MGRID=100000)
      PARAMETER (MICB=MGRID,MNFS=MGRID,MPRED=MGRID)
      LOGICAL LN1,LICB,LNFS,LPRED
      INTEGER N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,NREFL
      INTEGER IND(MIND),ICB(MICB),NFS(MNFS),IPRED(MPRED)
      INTEGER IPOSQ(MGRID),MINQ,MAXQ
      REAL X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3
      REAL P(MGRID),TT(MGRID),TTMIN,DMIN,ASX1,ASX2,ASX3
      COMMON/GRID/ N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,
     *             X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3,
     *             NREFL,LN1,LICB,LNFS,LPRED,IND,P,ICB,NFS,TT,IPRED,
     *             IPOSQ,MINQ,MAXQ,TTMIN,DMIN,ASX1,ASX2,ASX3
C
      INTEGER MFSMAX,MFS
      PARAMETER (MFSMAX=100,MFS=6626)
      INTEGER NFSMAX,KFS0(0:MFSMAX)
      INTEGER KFS1(MFS),KFS2(MFS),KFS3(MFS),KFS4(MFS),KFS5(MFS)
      REAL DFS(MFS)
      COMMON/FS/ NFSMAX,KFS0,KFS1,KFS2,KFS3,KFS4,KFS5,DFS
C
      INTEGER MSRC,MREC
      PARAMETER (MSRC=100,MREC=100)
      CHARACTER*6 NAMES,NAMER(MREC)
      INTEGER NSRC,NREC,IPREDR(MREC)
      REAL X1S(MSRC),X2S(MSRC),X3S(MSRC),TTS(MSRC),TTSERR(MSRC)
      REAL X1R(MREC),X2R(MREC),X3R(MREC),TTR(MREC)
      COMMON/NAMESR/ NAMES,NAMER
      COMMON/POINTS/ NSRC,X1S,X2S,X3S,TTS,TTSERR,NREC,X1R,X2R,X3R,TTR,
     *               IPREDR
C
      INTEGER MREFL
      PARAMETER (MREFL=2)
      CHARACTER*80 FSRC,FREC,FRAYS,FEND
      CHARACTER*80 FIND(0:MREFL),FVEL(0:MREFL),FICB(0:MREFL)
      CHARACTER*80 FTT(0:MREFL),FERR(0:MREFL),FPRED(0:MREFL)
      CHARACTER*80 FNFS(0:MREFL),FINTF(1:MREFL)
      COMMON/FILES/ FSRC,FREC,FRAYS,FEND,
     *              FIND,FVEL,FICB,FTT,FERR,FPRED,FINTF,FNFS
C
      INTEGER IPOS1,IPOS2,IPOS3,IPOS,IADR,ICBI,IFS
      REAL PI,TTI
      COMMON /NODE/ IPOS1,IPOS2,IPOS3,IPOS,IADR,PI,TTI,ICBI,IFS
C
      INTEGER ISRCI
      REAL DPOS1,DPOS2,DPOS3
      COMMON /SRCC/ ISRCI,DPOS1,DPOS2,DPOS3
C
C-----------------------------------------------------------------------
C
      EXTERNAL SRCREC,INDX
      INTEGER INDX
C
      INTEGER LU1
      PARAMETER (LU1=1)
      REAL GIANT
      PARAMETER (GIANT=1.0E+10)
C
      INTEGER N123,L123,ISRC,INDQ,NFSMAK,I,L,M
      INTEGER JPOS,JPOS1,JPOS2,JPOS3,JADR,NFSOPT,JCOUNT
      REAL VMAX,V,C1,C2MIN,C2MAX,AUX
      REAL P2(MPRED),P3(MNFS)
      EQUIVALENCE (IPRED,P2),(NFS,P3)
C
C     P2,P3... FOR FUTURE EXTENSIONS.
C
C.......................................................................
C
      WRITE(*,'(A)') '+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
      IF(FIND(IREFL).EQ.' ') THEN
        L1234=N1*N2*N3
        IF(L1.GT.1.OR.L2.GT.1.OR.L3.GT.1) THEN
          PAUSE 'ERROR 29: NO INDEX FILE SPECIFIED'
          STOP
        END IF
      ELSE
        N123=N1*N2*N3
        L123=L1*L2*L3
        IF(N123.GT.MIND) THEN
          PAUSE 'ERROR 30: TOO MANY BIG BRICKS'
          STOP
        END IF
        DO 11 L=1,N123
          IND(L)=L
   11   CONTINUE
        OPEN(LU1,FILE=FIND(IREFL),STATUS='OLD')
        READ(LU1,*) (IND(L),L=1,N123)
        CLOSE(LU1)
        DO 12 L=1,N123
          L4=MAX0(IND(L),L4)
          IND(L)=(IND(L)-1)*L123+1
   12   CONTINUE
        L1234=L1*L2*L3*L4
      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
        PAUSE 'ERROR 31: TOO MANY NETWORK NODES'
        STOP
      END IF
      IF(IREFL.GT.0) THEN
        IF(FERR(IREFL).NE.' ') THEN
          PAUSE
     *      'ERROR 43: TIME ERRORS CANNOT BE COMPUTED AFTER REFLECTION.'
          STOP
        END IF
        IF(FEND.NE.' ') THEN
          PAUSE
     *    'WARNING 44: TIME ERRORS CANNOT BE COMPUTED AFTER REFLECTION.'
        END IF
      END IF
      DO 21 L=1,L1234
        P(L)=1.
   21 CONTINUE
      OPEN(LU1,FILE=FVEL(IREFL),STATUS='OLD')
      READ(LU1,*) (P(L),L=1,L1234)
      VMAX=0.
      AUX=1.001/GIANT
      DO 22 L=1,L1234
        V=P(L)
        IF(ABS(V).LT.AUX) THEN
          P(L)=GIANT
        ELSE
          VMAX=AMAX1(V,VMAX)
          P(L)=1.0/V
        ENDIF
   22 CONTINUE
      CLOSE(LU1)
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
          PAUSE 'ERROR 32: INSUFFICIENT MEMORY FOR INDICES OF BLOCKS'
          STOP
        END IF
        DO 31 L=1,L1234
          ICB(L)=1
   31   CONTINUE
        OPEN(LU1,FILE=FICB(IREFL),STATUS='OLD')
        READ(LU1,*) (ICB(L),L=1,L1234)
        CLOSE(LU1)
      END IF
C
C     SWITCH FOR STORING PREDECESSORS:
   40 CONTINUE
      IF(NFSMAX.LE.0) THEN
C       PRIVATE BRANCH FOR FUTURE EXTENSIONS.
        LPRED=.FALSE.
        IF(L1234.GT.MPRED) THEN
          PAUSE 'ERROR 49: INSUFFICIENT MEMORY FOR P2'
          STOP
        END IF
        DO 42 L=1,L1234
          P2(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
          PAUSE 'ERROR 33: INSUFFICIENT MEMORY FOR PREDECESSORS'
          STOP
        END IF
      END IF
C
C     COMPUTING THE FORWARD STAR SIZES:
      IF(NFSMAX.LE.0) THEN
C       PRIVATE BRANCH FOR FUTURE EXTENSIONS.
        LNFS=.FALSE.
        IF(NL3.GT.1.AND.L1234.GT.MNFS) THEN
          PAUSE 'ERROR 50: INSUFFICIENT MEMORY FOR P3'
          STOP
        END IF
        DO 43 L=1,L1234
          P3(L)=0.
   43   CONTINUE
      ELSE IF(FNFS(IREFL)(1:1).EQ.'*') THEN
        LNFS=.FALSE.
      ELSE
        LNFS=.TRUE.
        IF(L1234.GT.MNFS) THEN
          PAUSE 'ERROR 34: INSUFFICIENT MEMORY FOR FORWARD STAR SIZES'
          STOP
        END IF
        IF(FNFS(IREFL)(1:1).EQ.'+') THEN
          DO 51 L=1,L1234
            NFS(L)=NFSMAX
   51     CONTINUE
          OPEN(LU1,FILE=FNFS(IREFL)(2:LEN(FNFS(IREFL))),STATUS='OLD')
          READ(LU1,*) (NFS(L),L=1,L1234)
          CLOSE(LU1)
          DO 52 L=1,L1234
            IF(NFS(L).GT.NFSMAX) THEN
              NFS(L)=NFSMAX
            ENDIF
   52     CONTINUE
        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)
                  NFS(JADR)=NFSOPT
                  JCOUNT=JCOUNT+1
                  IF(JCOUNT/1000*1000.EQ.JCOUNT) THEN
                    WRITE(*,'(A,I7,A,I7)')
     *                    '+GENERATING OPTIMUM SIZES OF FORWARD STARS:',
     *                    JCOUNT,' OF',L1234
                  END IF
                END IF
   56         CONTINUE
   57       CONTINUE
   58     CONTINUE
          WRITE(*,'(A,I7,A,I7)')
     *                    '+GENERATING OPTIMUM SIZES OF FORWARD STARS:',
     *                    JCOUNT,' OF',L1234
          IF(FNFS(IREFL).NE.' ') THEN
            OPEN(LU1,FILE=FNFS(IREFL)(2:LEN(FNFS(IREFL))))
            WRITE(LU1,'(20I3)') (NFS(L),L=1,L1234)
            CLOSE(LU1)
          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
            PAUSE 'ERROR 35: ONE-POINT MODEL'
            STOP
          ELSE
            PAUSE 'ERROR 36: LINE MODEL'
            STOP
          END IF
        ELSE
          IF(NL1.EQ.1) THEN
            PAUSE 'ERROR 37: LINE MODEL'
            STOP
          ELSE
            DMIN=AMIN1(ASX1,ASX2)/VMAX
          END IF
        END IF
      ELSE
        IF(NL2.EQ.1) THEN
          IF(NL1.EQ.1) THEN
            PAUSE 'ERROR 38: LINE MODEL'
            STOP
          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.LE.0) THEN
C       PRIVATE BRANCH FOR FUTURE EXTENSIONS.
        AUX=ASX1*ASX1+ASX2*ASX2+ASX3*ASX3
        DMIN=SQRT(AUX)-SQRT(AUX-DMIN*DMIN)
      END IF
C
C.......................................................................
C
C     INITIALIZATION:
C
      WRITE(*,'(A)')
     *  '+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
            TT(IADR)=GIANT
            IF(LPRED) THEN
              IPRED(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
          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         ADJUSTING EXTENT OF THE DENSE FORWARD STAR AT THE SOURCE POINT
          IF(LNFS) THEN
            IF(NFS(IADR).NE.NFSMAK) THEN
              NFSMAK=NFS(IADR)
              CALL MAKEFS(NFSMAK)
            END IF
          END IF
C
C         UPDATING:
          CALL LOOP(SRCREC)
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
          IPOSQ(INDQ)=0
   81   CONTINUE
        OPEN(LU1,FILE=FINTF(IREFL))
        READ(LU1,*) (IPOSQ(INDQ),INDQ=1,MGRID)
        CLOSE(LU1)
        DO 82 INDQ=1,MGRID
          IF(IPOSQ(INDQ).LE.0) THEN
            MAXQ=INDQ-1
            GO TO 83
          END IF
   82   CONTINUE
        PAUSE 'ERROR 39: TOO MANY POINTS OF REFLECTOR'
        STOP
   83   CONTINUE
C
C       LABELING THE TIME FIELD IN THE QUEUE (TT(I).LT.0):
        M=0
        TTMIN=GIANT
        DO 84 I=1,MAXQ
          IADR=INDX(IPOSQ(I))
C         CHECK FOR THE COMPUTATIONAL VOLUME AND FOR A FREE SPACE:
          IF(IADR.GT.0) THEN
            IF(P(IADR).LT.GIANT) THEN
              M=M+1
              IPOSQ(M)=IPOSQ(I)
              TTMIN=AMIN1(TT(IADR),TTMIN)
              TT(IADR)=-TT(IADR)
            END IF
          END IF
   84   CONTINUE
        MINQ=1
        MAXQ=M
C
C       DEFINING THE 1ST APPOXIMATION 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(TT(IADR).LT.0..AND.P(IADR).LT.GIANT) THEN
              TT(IADR)=-TT(IADR)
              IF(LPRED) THEN
                IPRED(IADR)=0
              ENDIF
            ELSE
              TT(IADR)=GIANT
              IF(LPRED) THEN
                IPRED(IADR)=I
              ENDIF
            ENDIF
          ENDIF
   88   CONTINUE
C
      ENDIF
      IF(MAXQ.LT.MINQ) THEN
        PAUSE 'ERROR 47: NO SOURCE POINT'
        STOP
      ENDIF
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE GENER(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
      INTEGER MIND,MGRID,MICB,MNFS,MPRED
      PARAMETER (MIND=1,MGRID=100000)
      PARAMETER (MICB=MGRID,MNFS=MGRID,MPRED=MGRID)
      LOGICAL LN1,LICB,LNFS,LPRED
      INTEGER N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,NREFL
      INTEGER IND(MIND),ICB(MICB),NFS(MNFS),IPRED(MPRED)
      INTEGER IPOSQ(MGRID),MINQ,MAXQ
      REAL X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3
      REAL P(MGRID),TT(MGRID),TTMIN,DMIN,ASX1,ASX2,ASX3
      COMMON/GRID/ N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,
     *             X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3,
     *             NREFL,LN1,LICB,LNFS,LPRED,IND,P,ICB,NFS,TT,IPRED,
     *             IPOSQ,MINQ,MAXQ,TTMIN,DMIN,ASX1,ASX2,ASX3
C
      INTEGER MFSMAX,MFS
      PARAMETER (MFSMAX=100,MFS=6626)
      INTEGER NFSMAX,KFS0(0:MFSMAX)
      INTEGER KFS1(MFS),KFS2(MFS),KFS3(MFS),KFS4(MFS),KFS5(MFS)
      REAL DFS(MFS)
      COMMON/FS/ NFSMAX,KFS0,KFS1,KFS2,KFS3,KFS4,KFS5,DFS
C
      INTEGER IPOS1,IPOS2,IPOS3,IPOS,IADR,ICBI,IFS
      REAL PI,TTI
      COMMON /NODE/ IPOS1,IPOS2,IPOS3,IPOS,IADR,PI,TTI,ICBI,IFS
C
C     IADR... INDEX OF THE NETWORK NODE (LOOP OR TEMPORARY VARIABLE).
C
C-----------------------------------------------------------------------
C
      REAL GIANT
      PARAMETER (GIANT=1.E+10)
C
      EXTERNAL INDX,UPDATE,FREVOL
      INTEGER  INDX
C
      INTEGER NUMNOD,NUMA,INDQ,I
      REAL TTINDQ
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     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)')
     *  '+ELIMINATING NODES IN A FREE SPACE.                           '
      NUMNOD=L1234
      DO 8 I=1,L1234
        IF(P(I).GE.GIANT) THEN
          TT(I)=-GIANT
          NUMNOD=NUMNOD-1
        ENDIF
    8 CONTINUE
C
C     COUNTER OF FINISHED NODES (NODES MOVED TO SET A):
      NUMA=0
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=TT(INDX(IPOSQ(INDQ)))
          IF(0..LT.TTINDQ) THEN
            MINQ=INDQ
            GO TO 22
          END IF
   21   CONTINUE
        MINQ=MAXQ+1
   22   CONTINUE
C
C       SCREEN OUTPUT:
        WRITE(*,'(A,I2,4(A,I7))')
     *        '+',IREFL,'-TH REFLECTION',NUMA,' NODES OF',NUMNOD,
     *        ' FINISHED,  QMIN=',MINQ,',  QMAX=',MAXQ
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
            TT(IADR)=ABS(TT(IADR))
   31     CONTINUE
          RETURN
        ENDIF
C
        DO 40 INDQ=MINQ,MAXQ
C
C         (2) SELECTION:
          TTINDQ=TT(INDX(IPOSQ(INDQ)))
          IF(0..LT.TTINDQ.AND.TTINDQ.LT.TTMIN) THEN
C
C           NEW NODE 'I' (POSITION IPOS, ADDRESS IADR):
            IPOS=IPOSQ(INDQ)
            IPOS1=IPOS-1
            IPOS3=IPOS1/(NL1*NL2)
            IPOS2=IPOS1/NL1-IPOS3*NL2
            IPOS1=IPOS1-(IPOS2+IPOS3*NL2)*NL1
            IADR=INDX(IPOS)
            PI=P(IADR)
            TTI=TT(IADR)
C
C           (3) UPDATING NODES 'J' OF THE FORWARD STAR FS(I), WHICH ARE
C           THE NEIGHBOURS TO THE NODE 'I':
            IF(L4.EQ.0) THEN
              CALL LOOP(UPDATE)
            ELSE
              CALL LOOP(FREVOL)
            END IF
C
C           MOVING 'I' FROM SET B TO SET A:
            IADR=INDX(IPOS)
            TT(IADR)=-TT(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
      SUBROUTINE RESFIL(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
      INTEGER MIND,MGRID,MICB,MNFS,MPRED
      PARAMETER (MIND=1,MGRID=100000)
      PARAMETER (MICB=MGRID,MNFS=MGRID,MPRED=MGRID)
      LOGICAL LN1,LICB,LNFS,LPRED
      INTEGER N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,NREFL
      INTEGER IND(MIND),ICB(MICB),NFS(MNFS),IPRED(MPRED)
      INTEGER IPOSQ(MGRID),MINQ,MAXQ
      REAL X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3
      REAL P(MGRID),TT(MGRID),TTMIN,DMIN,ASX1,ASX2,ASX3
      COMMON/GRID/ N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,
     *             X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3,
     *             NREFL,LN1,LICB,LNFS,LPRED,IND,P,ICB,NFS,TT,IPRED,
     *             IPOSQ,MINQ,MAXQ,TTMIN,DMIN,ASX1,ASX2,ASX3
C
      INTEGER MREFL
      PARAMETER (MREFL=2)
      CHARACTER*80 FSRC,FREC,FRAYS,FEND
      CHARACTER*80 FIND(0:MREFL),FVEL(0:MREFL),FICB(0:MREFL)
      CHARACTER*80 FTT(0:MREFL),FERR(0:MREFL),FPRED(0:MREFL)
      CHARACTER*80 FNFS(0:MREFL),FINTF(1:MREFL)
      COMMON/FILES/ FSRC,FREC,FRAYS,FEND,
     *              FIND,FVEL,FICB,FTT,FERR,FPRED,FINTF,FNFS
C
C-----------------------------------------------------------------------
C
      INTEGER LU0,LU
      PARAMETER(LU0=10)
C
      INTEGER I
C
C     I...    TEMPORARY LOOP VARIABLE.
C
C.......................................................................
C
      WRITE(*,'(A)')
     *  ' 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) TT(I),IPRED(I)
11        CONTINUE
        ELSE
          DO 12 I=1,L1234
            WRITE(LU,REC=I) TT(I)
12        CONTINUE
        END IF
      END IF
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE ERRORS(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
      INTEGER MIND,MGRID,MICB,MNFS,MPRED
      PARAMETER (MIND=1,MGRID=100000)
      PARAMETER (MICB=MGRID,MNFS=MGRID,MPRED=MGRID)
      LOGICAL LN1,LICB,LNFS,LPRED
      INTEGER N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,NREFL
      INTEGER IND(MIND),ICB(MICB),NFS(MNFS),IPRED(MPRED)
      INTEGER IPOSQ(MGRID),MINQ,MAXQ
      REAL X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3
      REAL P(MGRID),TT(MGRID),TTMIN,DMIN,ASX1,ASX2,ASX3
      COMMON/GRID/ N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,
     *             X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3,
     *             NREFL,LN1,LICB,LNFS,LPRED,IND,P,ICB,NFS,TT,IPRED,
     *             IPOSQ,MINQ,MAXQ,TTMIN,DMIN,ASX1,ASX2,ASX3
C
      INTEGER MFSMAX,MFS
      PARAMETER (MFSMAX=100,MFS=6626)
      INTEGER NFSMAX,KFS0(0:MFSMAX)
      INTEGER KFS1(MFS),KFS2(MFS),KFS3(MFS),KFS4(MFS),KFS5(MFS)
      REAL DFS(MFS)
      COMMON/FS/ NFSMAX,KFS0,KFS1,KFS2,KFS3,KFS4,KFS5,DFS
C
      INTEGER MREFL
      PARAMETER (MREFL=2)
      CHARACTER*80 FSRC,FREC,FRAYS,FEND
      CHARACTER*80 FIND(0:MREFL),FVEL(0:MREFL),FICB(0:MREFL)
      CHARACTER*80 FTT(0:MREFL),FERR(0:MREFL),FPRED(0:MREFL)
      CHARACTER*80 FNFS(0:MREFL),FINTF(1:MREFL)
      COMMON/FILES/ FSRC,FREC,FRAYS,FEND,
     *              FIND,FVEL,FICB,FTT,FERR,FPRED,FINTF,FNFS
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 ERR(MGRID),DUMMY
      EQUIVALENCE (ERR,IPOSQ)
C
C.......................................................................
C
      IF(NFSMAX.LE.0) THEN
C       THIS IS INACCESSIBLE PRIVATE BRANCH FOR FUTURE EXTENSIONS.
        RETURN
      END IF
C
      WRITE(*,'(A)')
     *  '+GENERATING TRAVEL-TIME ERRORS.                               '
C
      IF(FERR(IREFL).NE.' ') THEN
        LTRACE=.FALSE.
        LRAYS=.FALSE.
        LEND=.FALSE.
C
        DO 10 IEND=1,L1234
          ERR(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(ERR(IENDA).LT.0.) THEN
                  IPREDE=IPRED(IENDA)
                  IF(IPREDE.NE.IEND) THEN
                    CALL ONERAY(LTRACE,LRAYS,LEND,IAUX,IREFL,IPREDE,
     *                 IEND1,IEND2,IEND3,IENDA,TT(IENDA),ERR,DUMMY,IAUX)
                  END IF
                END IF
                JCOUNT=JCOUNT+1
                IF(JCOUNT/1000*1000.EQ.JCOUNT) THEN
                  WRITE(*,'(A,I7,A,I7)')
     *              '+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)')
     *  '+WRITING TRAVEL-TIME ERRORS.                                  '
C
        CALL WARRAY(LUERR,FERR(IREFL),'FORMATTED',
     *                                  .TRUE.,-1.,.FALSE.,0.,L1234,ERR)
C
      END IF
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE RECVRS()
C
C PROCEDURE UPDATING ARRIVAL TIMES AT THE RECEIVERS.
C
C NO INPUT, NO OUTPUT.
C
C-----------------------------------------------------------------------
C
      INTEGER MIND,MGRID,MICB,MNFS,MPRED
      PARAMETER (MIND=1,MGRID=100000)
      PARAMETER (MICB=MGRID,MNFS=MGRID,MPRED=MGRID)
      LOGICAL LN1,LICB,LNFS,LPRED
      INTEGER N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,NREFL
      INTEGER IND(MIND),ICB(MICB),NFS(MNFS),IPRED(MPRED)
      INTEGER IPOSQ(MGRID),MINQ,MAXQ
      REAL X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3
      REAL P(MGRID),TT(MGRID),TTMIN,DMIN,ASX1,ASX2,ASX3
      COMMON/GRID/ N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,
     *             X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3,
     *             NREFL,LN1,LICB,LNFS,LPRED,IND,P,ICB,NFS,TT,IPRED,
     *             IPOSQ,MINQ,MAXQ,TTMIN,DMIN,ASX1,ASX2,ASX3
C
      INTEGER MFSMAX,MFS
      PARAMETER (MFSMAX=100,MFS=6626)
      INTEGER NFSMAX,KFS0(0:MFSMAX)
      INTEGER KFS1(MFS),KFS2(MFS),KFS3(MFS),KFS4(MFS),KFS5(MFS)
      REAL DFS(MFS)
      COMMON/FS/ NFSMAX,KFS0,KFS1,KFS2,KFS3,KFS4,KFS5,DFS
C
      INTEGER MSRC,MREC
      PARAMETER (MSRC=100,MREC=100)
      CHARACTER*6 NAMES,NAMER(MREC)
      INTEGER NSRC,NREC,IPREDR(MREC)
      REAL X1S(MSRC),X2S(MSRC),X3S(MSRC),TTS(MSRC),TTSERR(MSRC)
      REAL X1R(MREC),X2R(MREC),X3R(MREC),TTR(MREC)
      COMMON/NAMESR/ NAMES,NAMER
      COMMON/POINTS/ NSRC,X1S,X2S,X3S,TTS,TTSERR,NREC,X1R,X2R,X3R,TTR,
     *               IPREDR
C
      INTEGER MREFL
      PARAMETER (MREFL=2)
      CHARACTER*80 FSRC,FREC,FRAYS,FEND
      CHARACTER*80 FIND(0:MREFL),FVEL(0:MREFL),FICB(0:MREFL)
      CHARACTER*80 FTT(0:MREFL),FERR(0:MREFL),FPRED(0:MREFL)
      CHARACTER*80 FNFS(0:MREFL),FINTF(1:MREFL)
      COMMON/FILES/ FSRC,FREC,FRAYS,FEND,
     *              FIND,FVEL,FICB,FTT,FERR,FPRED,FINTF,FNFS
C
      INTEGER IPOS1,IPOS2,IPOS3,IPOS,IADR,ICBI,IFS
      REAL PI,TTI
      COMMON /NODE/ IPOS1,IPOS2,IPOS3,IPOS,IADR,PI,TTI,ICBI,IFS
C
      INTEGER ISRCI
      REAL DPOS1,DPOS2,DPOS3
      COMMON /SRCC/ ISRCI,DPOS1,DPOS2,DPOS3
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 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...
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
C     COMPUTING APPROXIMATE ARRIVAL TIME AND RAYS IN THE RECEIVER NODES:
      WRITE(*,'(A)')
     *  '+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=NFS(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(ICB(IADR).NE.ICB(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=NFS(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(NFS(IADR).NE.NFSMAK) THEN
              NFSMAK=NFS(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
      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
      INTEGER MIND,MGRID,MICB,MNFS,MPRED
      PARAMETER (MIND=1,MGRID=100000)
      PARAMETER (MICB=MGRID,MNFS=MGRID,MPRED=MGRID)
      LOGICAL LN1,LICB,LNFS,LPRED
      INTEGER N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,NREFL
      INTEGER IND(MIND),ICB(MICB),NFS(MNFS),IPRED(MPRED)
      INTEGER IPOSQ(MGRID),MINQ,MAXQ
      REAL X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3
      REAL P(MGRID),TT(MGRID),TTMIN,DMIN,ASX1,ASX2,ASX3
      COMMON/GRID/ N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,
     *             X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3,
     *             NREFL,LN1,LICB,LNFS,LPRED,IND,P,ICB,NFS,TT,IPRED,
     *             IPOSQ,MINQ,MAXQ,TTMIN,DMIN,ASX1,ASX2,ASX3
C
      INTEGER MFSMAX,MFS
      PARAMETER (MFSMAX=100,MFS=6626)
      INTEGER NFSMAX,KFS0(0:MFSMAX)
      INTEGER KFS1(MFS),KFS2(MFS),KFS3(MFS),KFS4(MFS),KFS5(MFS)
      REAL DFS(MFS)
      COMMON/FS/ NFSMAX,KFS0,KFS1,KFS2,KFS3,KFS4,KFS5,DFS
C
      INTEGER MSRC,MREC
      PARAMETER (MSRC=100,MREC=100)
      CHARACTER*6 NAMES,NAMER(MREC)
      INTEGER NSRC,NREC,IPREDR(MREC)
      REAL X1S(MSRC),X2S(MSRC),X3S(MSRC),TTS(MSRC),TTSERR(MSRC)
      REAL X1R(MREC),X2R(MREC),X3R(MREC),TTR(MREC)
      COMMON/NAMESR/ NAMES,NAMER
      COMMON/POINTS/ NSRC,X1S,X2S,X3S,TTS,TTSERR,NREC,X1R,X2R,X3R,TTR,
     *               IPREDR
C
      INTEGER MREFL
      PARAMETER (MREFL=2)
      CHARACTER*80 FSRC,FREC,FRAYS,FEND
      CHARACTER*80 FIND(0:MREFL),FVEL(0:MREFL),FICB(0:MREFL)
      CHARACTER*80 FTT(0:MREFL),FERR(0:MREFL),FPRED(0:MREFL)
      CHARACTER*80 FNFS(0:MREFL),FINTF(1:MREFL)
      COMMON/FILES/ FSRC,FREC,FRAYS,FEND,
     *              FIND,FVEL,FICB,FTT,FERR,FPRED,FINTF,FNFS
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 DPOS1,DPOS2,DPOS3,PEND,TERR,DUMMY
C
C     IREC... LOOP VARIABLE - INDEX OF A RECEIVER.
C
C.......................................................................
C
      IF(NFSMAX.LE.0) THEN
C       THIS IS INACCESSIBLE PRIVATE BRANCH FOR FUTURE EXTENSIONS.
        RETURN
      END IF
C
      LTRACE=.TRUE.
      WRITE(*,'(A)')
     *  '+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),
     *              DPOS1,DPOS2,DPOS3,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
      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
      INTEGER MIND,MGRID,MICB,MNFS,MPRED
      PARAMETER (MIND=1,MGRID=100000)
      PARAMETER (MICB=MGRID,MNFS=MGRID,MPRED=MGRID)
      LOGICAL LN1,LICB,LNFS,LPRED
      INTEGER N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,NREFL
      INTEGER IND(MIND),ICB(MICB),NFS(MNFS),IPRED(MPRED)
      INTEGER IPOSQ(MGRID),MINQ,MAXQ
      REAL X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3
      REAL P(MGRID),TT(MGRID),TTMIN,DMIN,ASX1,ASX2,ASX3
      COMMON/GRID/ N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,
     *             X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3,
     *             NREFL,LN1,LICB,LNFS,LPRED,IND,P,ICB,NFS,TT,IPRED,
     *             IPOSQ,MINQ,MAXQ,TTMIN,DMIN,ASX1,ASX2,ASX3
C
      INTEGER MFSMAX,MFS
      PARAMETER (MFSMAX=100,MFS=6626)
      INTEGER NFSMAX,KFS0(0:MFSMAX)
      INTEGER KFS1(MFS),KFS2(MFS),KFS3(MFS),KFS4(MFS),KFS5(MFS)
      REAL DFS(MFS)
      COMMON/FS/ NFSMAX,KFS0,KFS1,KFS2,KFS3,KFS4,KFS5,DFS
C
      INTEGER MREFL
      PARAMETER (MREFL=2)
      CHARACTER*80 FSRC,FREC,FRAYS,FEND
      CHARACTER*80 FIND(0:MREFL),FVEL(0:MREFL),FICB(0:MREFL)
      CHARACTER*80 FTT(0:MREFL),FERR(0:MREFL),FPRED(0:MREFL)
      CHARACTER*80 FNFS(0:MREFL),FINTF(1:MREFL)
      COMMON/FILES/ FSRC,FREC,FRAYS,FEND,
     *              FIND,FVEL,FICB,FTT,FERR,FPRED,FINTF,FNFS
C
C-----------------------------------------------------------------------
C
      INTEGER LU0,LU1,LU
      PARAMETER(LU0=10)
      PARAMETER(LU1=1)
      REAL GIANT
      PARAMETER (GIANT=1.0E+10)
C
      INTEGER NF,IADR
      REAL P1(MGRID),P2(MPRED),P3(MNFS)
      EQUIVALENCE (P,P1),(IPRED,P2),(NFS,P3)
C
C     NF,IADR... LOOP VARIABLES.
C     P1,P2,P3... FOR FUTURE EXTENSIONS.
C
C.......................................................................
C
      WRITE(*,'(A)')
     *  '+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 IADR=1,L1234
                READ(LU,REC=IADR)TT(IADR),IPRED(IADR)
   11         CONTINUE
            ELSE
              DO 12 IADR=1,L1234
                READ(LU,REC=IADR)TT(IADR)
   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,TT)
        END IF
C
C       WRITING PREDECESSORS:
        IF(NFSMAX.LE.0) THEN
C         THIS IS INACCESSIBLE PRIVATE BRANCH FOR FUTURE EXTENSIONS.
          IF(FERR(NF).NE.' ') THEN
            CALL WARRAY(LU1,FERR(NF),'FORMATTED',
     *                                .FALSE.,0.,.FALSE.,GIANT,L1234,P1)
          END IF
          IF(FPRED(NF).NE.' ') THEN
            CALL WARRAY(LU1,FPRED(NF),'FORMATTED',
     *                                .FALSE.,0.,.FALSE.,GIANT,L1234,P2)
          END IF
          IF(FNFS(NF).NE.' ') THEN
            CALL WARRAY(LU1,FNFS(NF),'FORMATTED',
     *                                .FALSE.,0.,.FALSE.,GIANT,L1234,P3)
          END IF
        ELSE IF(FPRED(NF).NE.' ') THEN
          OPEN(LU1,FILE=FPRED(NF))
          WRITE(LU1,'(10I8)') (IPRED(IADR),IADR=1,L1234)
          CLOSE(LU1)
        END IF
   30 CONTINUE
C
      WRITE(*,'(A)')
     *  '+END                                                          '
      RETURN
      END
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 EQUIAVALENT, 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
      INTEGER MIND,MGRID,MICB,MNFS,MPRED
      PARAMETER (MIND=1,MGRID=100000)
      PARAMETER (MICB=MGRID,MNFS=MGRID,MPRED=MGRID)
      LOGICAL LN1,LICB,LNFS,LPRED
      INTEGER N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,NREFL
      INTEGER IND(MIND),ICB(MICB),NFS(MNFS),IPRED(MPRED)
      INTEGER IPOSQ(MGRID),MINQ,MAXQ
      REAL X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3
      REAL P(MGRID),TT(MGRID),TTMIN,DMIN,ASX1,ASX2,ASX3
      COMMON/GRID/ N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,
     *             X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3,
     *             NREFL,LN1,LICB,LNFS,LPRED,IND,P,ICB,NFS,TT,IPRED,
     *             IPOSQ,MINQ,MAXQ,TTMIN,DMIN,ASX1,ASX2,ASX3
C
      INTEGER MFSMAX,MFS
      PARAMETER (MFSMAX=100,MFS=6626)
      INTEGER NFSMAX,KFS0(0:MFSMAX)
      INTEGER KFS1(MFS),KFS2(MFS),KFS3(MFS),KFS4(MFS),KFS5(MFS)
      REAL DFS(MFS)
      COMMON/FS/ NFSMAX,KFS0,KFS1,KFS2,KFS3,KFS4,KFS5,DFS
C
      INTEGER IPOS1,IPOS2,IPOS3,IPOS,IADR,ICBI,IFS
      REAL PI,TTI
      COMMON /NODE/ IPOS1,IPOS2,IPOS3,IPOS,IADR,PI,TTI,ICBI,IFS
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.LE.0) THEN
        PAUSE 'ERROR 48: NFSMAX NOT POSITIVE'
        STOP
      END IF
C
C     INDEX OF THE BLOCK (IF APPLICABLE):
      IF(LICB) THEN
        ICBI=ICB(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(NFS(IADR).NE.NFSOLD) THEN
          NFSOLD=NFS(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
      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
      INTEGER MIND,MGRID,MICB,MNFS,MPRED
      PARAMETER (MIND=1,MGRID=100000)
      PARAMETER (MICB=MGRID,MNFS=MGRID,MPRED=MGRID)
      LOGICAL LN1,LICB,LNFS,LPRED
      INTEGER N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,NREFL
      INTEGER IND(MIND),ICB(MICB),NFS(MNFS),IPRED(MPRED)
      INTEGER IPOSQ(MGRID),MINQ,MAXQ
      REAL X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3
      REAL P(MGRID),TT(MGRID),TTMIN,DMIN,ASX1,ASX2,ASX3
      COMMON/GRID/ N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,
     *             X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3,
     *             NREFL,LN1,LICB,LNFS,LPRED,IND,P,ICB,NFS,TT,IPRED,
     *             IPOSQ,MINQ,MAXQ,TTMIN,DMIN,ASX1,ASX2,ASX3
C
      INTEGER MFSMAX,MFS
      PARAMETER (MFSMAX=100,MFS=6626)
      INTEGER NFSMAX,KFS0(0:MFSMAX)
      INTEGER KFS1(MFS),KFS2(MFS),KFS3(MFS),KFS4(MFS),KFS5(MFS)
      REAL DFS(MFS)
      COMMON/FS/ NFSMAX,KFS0,KFS1,KFS2,KFS3,KFS4,KFS5,DFS
C
      INTEGER IPOS1,IPOS2,IPOS3,IPOS,IADR,ICBI,IFS
      REAL PI,TTI
      COMMON /NODE/ IPOS1,IPOS2,IPOS3,IPOS,IADR,PI,TTI,ICBI,IFS
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=IND(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
      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=TT(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(ICB(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.NFS(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)*(P(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
            IPOSQ(MAXQ)=IPOS+KFS4(IFS)
          END IF
C
C         UPDATING
          TT(JADR)=TTIJ
          IF(LPRED) THEN
            IPRED(JADR)=IPOS
          END IF
C
C       END OF IF
C     END OF IF
      RETURN
      END
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
      INTEGER MIND,MGRID,MICB,MNFS,MPRED
      PARAMETER (MIND=1,MGRID=100000)
      PARAMETER (MICB=MGRID,MNFS=MGRID,MPRED=MGRID)
      LOGICAL LN1,LICB,LNFS,LPRED
      INTEGER N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,NREFL
      INTEGER IND(MIND),ICB(MICB),NFS(MNFS),IPRED(MPRED)
      INTEGER IPOSQ(MGRID),MINQ,MAXQ
      REAL X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3
      REAL P(MGRID),TT(MGRID),TTMIN,DMIN,ASX1,ASX2,ASX3
      COMMON/GRID/ N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,
     *             X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3,
     *             NREFL,LN1,LICB,LNFS,LPRED,IND,P,ICB,NFS,TT,IPRED,
     *             IPOSQ,MINQ,MAXQ,TTMIN,DMIN,ASX1,ASX2,ASX3
C
      INTEGER MSRC,MREC
      PARAMETER (MSRC=100,MREC=100)
      CHARACTER*6 NAMES,NAMER(MREC)
      INTEGER NSRC,NREC,IPREDR(MREC)
      REAL X1S(MSRC),X2S(MSRC),X3S(MSRC),TTS(MSRC),TTSERR(MSRC)
      REAL X1R(MREC),X2R(MREC),X3R(MREC),TTR(MREC)
      COMMON/NAMESR/ NAMES,NAMER
      COMMON/POINTS/ NSRC,X1S,X2S,X3S,TTS,TTSERR,NREC,X1R,X2R,X3R,TTR,
     *               IPREDR
C
      INTEGER MFSMAX,MFS
      PARAMETER (MFSMAX=100,MFS=6626)
      INTEGER NFSMAX,KFS0(0:MFSMAX)
      INTEGER KFS1(MFS),KFS2(MFS),KFS3(MFS),KFS4(MFS),KFS5(MFS)
      REAL DFS(MFS)
      COMMON/FS/ NFSMAX,KFS0,KFS1,KFS2,KFS3,KFS4,KFS5,DFS
C
      INTEGER IPOS1,IPOS2,IPOS3,IPOS,IADR,ICBI,IFS
      REAL PI,TTI
      COMMON /NODE/ IPOS1,IPOS2,IPOS3,IPOS,IADR,PI,TTI,ICBI,IFS
C
      INTEGER ISRCI
      REAL DPOS1,DPOS2,DPOS3
      COMMON /SRCC/ ISRCI,DPOS1,DPOS2,DPOS3
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(P(JADR).GE.GIANT) RETURN
C
          TTJ=TT(JADR)
C
C         CHECK FOR CROSSING AN INTERFACE:
          IF(LICB) THEN
            IF(ICB(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.NFS(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*(P(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*(P(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
                IPOSQ(MAXQ)=IPOS+KFS4(IFS)
              END IF
              TT(JADR)=TTI+TTIJ
              IF(LPRED) THEN
                IPRED(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
      SUBROUTINE READFS()
C
C PROCEDURE CALLED BY SOURCE, TO READ AND CONSTRUCT OPTIMIZED TEMPLATE
C FORWARD STARS.
C
C-----------------------------------------------------------------------
C
      INTEGER MIND,MGRID,MICB,MNFS,MPRED
      PARAMETER (MIND=1,MGRID=100000)
      PARAMETER (MICB=MGRID,MNFS=MGRID,MPRED=MGRID)
      LOGICAL LN1,LICB,LNFS,LPRED
      INTEGER N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,NREFL
      INTEGER IND(MIND),ICB(MICB),NFS(MNFS),IPRED(MPRED)
      INTEGER IPOSQ(MGRID),MINQ,MAXQ
      REAL X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3
      REAL P(MGRID),TT(MGRID),TTMIN,DMIN,ASX1,ASX2,ASX3
      COMMON/GRID/ N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,
     *             X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3,
     *             NREFL,LN1,LICB,LNFS,LPRED,IND,P,ICB,NFS,TT,IPRED,
     *             IPOSQ,MINQ,MAXQ,TTMIN,DMIN,ASX1,ASX2,ASX3
C
      INTEGER MFSMAX,MFS
      PARAMETER (MFSMAX=100,MFS=6626)
      INTEGER NFSMAX,KFS0(0:MFSMAX)
      INTEGER KFS1(MFS),KFS2(MFS),KFS3(MFS),KFS4(MFS),KFS5(MFS)
      REAL DFS(MFS)
      COMMON/FS/ NFSMAX,KFS0,KFS1,KFS2,KFS3,KFS4,KFS5,DFS
C
C-----------------------------------------------------------------------
C
      INTEGER LU1
      PARAMETER(LU1=1)
C
      INTEGER MD
      PARAMETER (MD=169)
C
      INTEGER ID(MD),JD(MD),KD(MD)
      INTEGER IFS,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     IFS...  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)')
     *  '+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 IFS=1,NFSMAX
        READ(LU1,*) I,ND
        IF(I.NE.IFS) THEN
          PAUSE 'ERROR 40: TEMPLATE FORWARD STAR NOT FOUND IN THE FILE'
          STOP
        END IF
        IF(ND.GT.MD) THEN
          PAUSE 'ERROR 41: TOO MANY INPUT NODES OF A FORWARD STAR'
          STOP
        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(IFS)=KFS0(IFS-1)
        DO 69 L=1,ND
          I1=ID(L)
          I2=JD(L)
          I3=KD(L)
          CALL SETFS(IFS,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:
*     OPEN(LU1,FILE='NET.TAB')
*     WRITE(LU1,'(A)') '      SIZE    NUMBER    SUM OF'
*     WRITE(LU1,'(A)') '     (NFS)  OF NODES   NUMBERS'
*     WRITE(LU1,'(3I10)')
*    *                (IFS,KFS0(IFS)-KFS0(IFS-1),KFS0(IFS),IFS=1,NFSMAX)
*     CLOSE(LU1)
*     PAUSE 'THE TABLE HAS BEEN WRITTEN TO THE FILE ''NET.TAB'''
*     STOP
C     THE CENTRAL NODE OF A FORWARD STAR IS NOT STORED IN THE MEMORY
C     AND IS NOT COUNTED.
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE MAKEFS(NFSMAK)
C
C PROCEDURE CALLED BY SOURCE AND RECVRS, TO CONSTRUCT FULL SPHERICAL
C TEMPLATE FORWARD STARS.
C
C-----------------------------------------------------------------------
C
      INTEGER MIND,MGRID,MICB,MNFS,MPRED
      PARAMETER (MIND=1,MGRID=100000)
      PARAMETER (MICB=MGRID,MNFS=MGRID,MPRED=MGRID)
      LOGICAL LN1,LICB,LNFS,LPRED
      INTEGER N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,NREFL
      INTEGER IND(MIND),ICB(MICB),NFS(MNFS),IPRED(MPRED)
      INTEGER IPOSQ(MGRID),MINQ,MAXQ
      REAL X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3
      REAL P(MGRID),TT(MGRID),TTMIN,DMIN,ASX1,ASX2,ASX3
      COMMON/GRID/ N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,
     *             X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3,
     *             NREFL,LN1,LICB,LNFS,LPRED,IND,P,ICB,NFS,TT,IPRED,
     *             IPOSQ,MINQ,MAXQ,TTMIN,DMIN,ASX1,ASX2,ASX3
C
      INTEGER MFSMAX,MFS
      PARAMETER (MFSMAX=100,MFS=6626)
      INTEGER NFSMAX,KFS0(0:MFSMAX)
      INTEGER KFS1(MFS),KFS2(MFS),KFS3(MFS),KFS4(MFS),KFS5(MFS)
      REAL DFS(MFS)
      COMMON/FS/ NFSMAX,KFS0,KFS1,KFS2,KFS3,KFS4,KFS5,DFS
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
      SUBROUTINE SETFS(IFS,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
      INTEGER MIND,MGRID,MICB,MNFS,MPRED
      PARAMETER (MIND=1,MGRID=100000)
      PARAMETER (MICB=MGRID,MNFS=MGRID,MPRED=MGRID)
      LOGICAL LN1,LICB,LNFS,LPRED
      INTEGER N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,NREFL
      INTEGER IND(MIND),ICB(MICB),NFS(MNFS),IPRED(MPRED)
      INTEGER IPOSQ(MGRID),MINQ,MAXQ
      REAL X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3
      REAL P(MGRID),TT(MGRID),TTMIN,DMIN,ASX1,ASX2,ASX3
      COMMON/GRID/ N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,
     *             X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3,
     *             NREFL,LN1,LICB,LNFS,LPRED,IND,P,ICB,NFS,TT,IPRED,
     *             IPOSQ,MINQ,MAXQ,TTMIN,DMIN,ASX1,ASX2,ASX3
C
      INTEGER MFSMAX,MFS
      PARAMETER (MFSMAX=100,MFS=6626)
      INTEGER NFSMAX,KFS0(0:MFSMAX)
      INTEGER KFS1(MFS),KFS2(MFS),KFS3(MFS),KFS4(MFS),KFS5(MFS)
      REAL DFS(MFS)
      COMMON/FS/ NFSMAX,KFS0,KFS1,KFS2,KFS3,KFS4,KFS5,DFS
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(IFS),...,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(IFS)
      IF(NL1.EQ.1.OR.NL2.EQ.1.OR.NL3.EQ.1) THEN
        DIM=2.
      ELSE
        DIM=3.
      END IF
      NLEX1=MIN0(IFS,NL1-1)
      NLEX2=MIN0(IFS,NL2-1)
      NLEX3=MIN0(IFS,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(IFS)+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(IFS)+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(IFS)+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(IFS)+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(IFS)=LFS
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE STORFS(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
      INTEGER MFSMAX,MFS
      PARAMETER (MFSMAX=100,MFS=6626)
      INTEGER NFSMAX,KFS0(0:MFSMAX)
      INTEGER KFS1(MFS),KFS2(MFS),KFS3(MFS),KFS4(MFS),KFS5(MFS)
      REAL DFS(MFS)
      COMMON/FS/ NFSMAX,KFS0,KFS1,KFS2,KFS3,KFS4,KFS5,DFS
C
C-----------------------------------------------------------------------
C
      LFS=LFS+1
      IF(LFS.GT.MFS) THEN
        PAUSE 'ERROR 42: TOO BIG FORWARD STAR'
        STOP
      END IF
      KFS1(LFS)=I1
      KFS2(LFS)=I2
      KFS3(LFS)=I3
      RETURN
      END
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
      INTEGER MIND,MGRID,MICB,MNFS,MPRED
      PARAMETER (MIND=1,MGRID=100000)
      PARAMETER (MICB=MGRID,MNFS=MGRID,MPRED=MGRID)
      LOGICAL LN1,LICB,LNFS,LPRED
      INTEGER N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,NREFL
      INTEGER IND(MIND),ICB(MICB),NFS(MNFS),IPRED(MPRED)
      INTEGER IPOSQ(MGRID),MINQ,MAXQ
      REAL X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3
      REAL P(MGRID),TT(MGRID),TTMIN,DMIN,ASX1,ASX2,ASX3
      COMMON/GRID/ N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,
     *             X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3,
     *             NREFL,LN1,LICB,LNFS,LPRED,IND,P,ICB,NFS,TT,IPRED,
     *             IPOSQ,MINQ,MAXQ,TTMIN,DMIN,ASX1,ASX2,ASX3
C
      INTEGER MSRC,MREC
      PARAMETER (MSRC=100,MREC=100)
      CHARACTER*6 NAMES,NAMER(MREC)
      INTEGER NSRC,NREC,IPREDR(MREC)
      REAL X1S(MSRC),X2S(MSRC),X3S(MSRC),TTS(MSRC),TTSERR(MSRC)
      REAL X1R(MREC),X2R(MREC),X3R(MREC),TTR(MREC)
      COMMON/NAMESR/ NAMES,NAMER
      COMMON/POINTS/ NSRC,X1S,X2S,X3S,TTS,TTSERR,NREC,X1R,X2R,X3R,TTR,
     *               IPREDR
C
      INTEGER MREFL
      PARAMETER (MREFL=2)
      CHARACTER*80 FSRC,FREC,FRAYS,FEND
      CHARACTER*80 FIND(0:MREFL),FVEL(0:MREFL),FICB(0:MREFL)
      CHARACTER*80 FTT(0:MREFL),FERR(0:MREFL),FPRED(0:MREFL)
      CHARACTER*80 FNFS(0:MREFL),FINTF(1:MREFL)
      COMMON/FILES/ FSRC,FREC,FRAYS,FEND,
     *              FIND,FVEL,FICB,FTT,FERR,FPRED,FINTF,FNFS
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,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=ICB(IENDA)
        ELSE
          ICBOLD=1
        END IF
        IOLD1=IEND1
        IOLD2=IEND2
        IOLD3=IEND3
        POLD=P(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
                IND(K)=K
   11         CONTINUE
              OPEN(LU1,FILE=FIND(NF),STATUS='OLD')
              READ(LU1,*) (IND(K),K=1,N123)
              CLOSE(LU1)
              DO 12 K=1,N123
                IND(K)=(IND(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=TT(IADR)
              IPNEW=IPRED(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 ERRROR 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 ERRROR 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
      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
      INTEGER MIND,MGRID,MICB,MNFS,MPRED
      PARAMETER (MIND=1,MGRID=100000)
      PARAMETER (MICB=MGRID,MNFS=MGRID,MPRED=MGRID)
      LOGICAL LN1,LICB,LNFS,LPRED
      INTEGER N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,NREFL
      INTEGER IND(MIND),ICB(MICB),NFS(MNFS),IPRED(MPRED)
      INTEGER IPOSQ(MGRID),MINQ,MAXQ
      REAL X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3
      REAL P(MGRID),TT(MGRID),TTMIN,DMIN,ASX1,ASX2,ASX3
      COMMON/GRID/ N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,
     *             X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3,
     *             NREFL,LN1,LICB,LNFS,LPRED,IND,P,ICB,NFS,TT,IPRED,
     *             IPOSQ,MINQ,MAXQ,TTMIN,DMIN,ASX1,ASX2,ASX3
C
      INTEGER MFSMAX,MFS
      PARAMETER (MFSMAX=100,MFS=6626)
      INTEGER NFSMAX,KFS0(0:MFSMAX)
      INTEGER KFS1(MFS),KFS2(MFS),KFS3(MFS),KFS4(MFS),KFS5(MFS)
      REAL DFS(MFS)
      COMMON/FS/ NFSMAX,KFS0,KFS1,KFS2,KFS3,KFS4,KFS5,DFS
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=NFS(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.*P(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=ICB(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=(P(IADR)+POLD)/2.
        ERRMIN=ERRMIN
     *         -SQRT(AUX1*AUX1+AUX2*AUX2+AUX3*AUX3)*ABS(P(IADR)-POLD)/2.
        ERR3=AMIN1(P(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=(P(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=P(IADR)
      TOLD=T
C
      RETURN
      END
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             DEFUALT OF NFS=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
      INTEGER MIND,MGRID,MICB,MNFS,MPRED
      PARAMETER (MIND=1,MGRID=100000)
      PARAMETER (MICB=MGRID,MNFS=MGRID,MPRED=MGRID)
      LOGICAL LN1,LICB,LNFS,LPRED
      INTEGER N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,NREFL
      INTEGER IND(MIND),ICB(MICB),NFS(MNFS),IPRED(MPRED)
      INTEGER IPOSQ(MGRID),MINQ,MAXQ
      REAL X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3
      REAL P(MGRID),TT(MGRID),TTMIN,DMIN,ASX1,ASX2,ASX3
      COMMON/GRID/ N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,
     *             X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3,
     *             NREFL,LN1,LICB,LNFS,LPRED,IND,P,ICB,NFS,TT,IPRED,
     *             IPOSQ,MINQ,MAXQ,TTMIN,DMIN,ASX1,ASX2,ASX3
C
      INTEGER MFSMAX,MFS
      PARAMETER (MFSMAX=100,MFS=6626)
      INTEGER NFSMAX,KFS0(0:MFSMAX)
      INTEGER KFS1(MFS),KFS2(MFS),KFS3(MFS),KFS4(MFS),KFS5(MFS)
      REAL DFS(MFS)
      COMMON/FS/ NFSMAX,KFS0,KFS1,KFS2,KFS3,KFS4,KFS5,DFS
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(P(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.*P(IADR)
      C2MAX=AMAX1(UIJ(1),0.)/C2
      C2MIN=AMIN1(UIJ(6),0.)/C2
      C2=AMAX1(C2MAX,-C2MIN-C2MAX)
C
C     OPTIMUM SIZE OF THE FORWARD STAR:
      IF(C2.GT.0.) THEN
        NFSOPT=MAX0(1,MIN0(INT(SQRT(SQRT(C1/C2))+0.5),NFSMAX))
      ELSE
        NFSOPT=NFSMAX
      END IF
C
      RETURN
      END
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
      INTEGER MIND,MGRID,MICB,MNFS,MPRED
      PARAMETER (MIND=1,MGRID=100000)
      PARAMETER (MICB=MGRID,MNFS=MGRID,MPRED=MGRID)
      LOGICAL LN1,LICB,LNFS,LPRED
      INTEGER N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,NREFL
      INTEGER IND(MIND),ICB(MICB),NFS(MNFS),IPRED(MPRED)
      INTEGER IPOSQ(MGRID),MINQ,MAXQ
      REAL X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3
      REAL P(MGRID),TT(MGRID),TTMIN,DMIN,ASX1,ASX2,ASX3
      COMMON/GRID/ N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,
     *             X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3,
     *             NREFL,LN1,LICB,LNFS,LPRED,IND,P,ICB,NFS,TT,IPRED,
     *             IPOSQ,MINQ,MAXQ,TTMIN,DMIN,ASX1,ASX2,ASX3
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
      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
      INTEGER MIND,MGRID,MICB,MNFS,MPRED
      PARAMETER (MIND=1,MGRID=100000)
      PARAMETER (MICB=MGRID,MNFS=MGRID,MPRED=MGRID)
      LOGICAL LN1,LICB,LNFS,LPRED
      INTEGER N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,NREFL
      INTEGER IND(MIND),ICB(MICB),NFS(MNFS),IPRED(MPRED)
      INTEGER IPOSQ(MGRID),MINQ,MAXQ
      REAL X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3
      REAL P(MGRID),TT(MGRID),TTMIN,DMIN,ASX1,ASX2,ASX3
      COMMON/GRID/ N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,
     *             X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3,
     *             NREFL,LN1,LICB,LNFS,LPRED,IND,P,ICB,NFS,TT,IPRED,
     *             IPOSQ,MINQ,MAXQ,TTMIN,DMIN,ASX1,ASX2,ASX3
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
C     IERR/4 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=ICB(IADR)
        ELSE
          ICBOPT=1
        END IF
      END IF
C
C     CHECK FOR THE GEOLOGICAL BLOCK AT THE GIVEN NODE:
      IF(LICB) THEN
        IF(ICB(IADR).NE.ICBOPT) THEN
          RETURN
        END IF
      END IF
C     CHECK FOR A FREE SPACE:
      AUX1=P(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(ICB(IA0).NE.ICBOPT) THEN
            GO TO 20
          END IF
          IF(ICB(IA2).NE.ICBOPT) THEN
            GO TO 20
          END IF
        END IF
C       CHECK FOR A FREE SPACE:
        AUX0=P(IA0)
        AUX2=P(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(ICB(IA0).NE.ICBOPT) THEN
            GO TO 30
          END IF
          IF(ICB(IA2).NE.ICBOPT) THEN
            GO TO 30
          END IF
        END IF
C       CHECK FOR A FREE SPACE:
        AUX0=P(IA0)
        AUX2=P(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(ICB(IA0).NE.ICBOPT) THEN
            GO TO 40
          END IF
          IF(ICB(IA2).NE.ICBOPT) THEN
            GO TO 40
          END IF
        END IF
C       CHECK FOR A FREE SPACE:
        AUX0=P(IA0)
        AUX2=P(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/P(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
      SUBROUTINE MIXDER(IPOS0,IPOS1,IPOS2,ICBREF,DERMIX,IERR)
C
C AUXILIARY SUBROUTINE TO TRYMAT, TO CALCULATE MIXED SECOND PARTIAL
C SLOWNESS DERIVATIVES.
C
C-----------------------------------------------------------------------
C
      INTEGER MIND,MGRID,MICB,MNFS,MPRED
      PARAMETER (MIND=1,MGRID=100000)
      PARAMETER (MICB=MGRID,MNFS=MGRID,MPRED=MGRID)
      LOGICAL LN1,LICB,LNFS,LPRED
      INTEGER N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,NREFL
      INTEGER IND(MIND),ICB(MICB),NFS(MNFS),IPRED(MPRED)
      INTEGER IPOSQ(MGRID),MINQ,MAXQ
      REAL X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3
      REAL P(MGRID),TT(MGRID),TTMIN,DMIN,ASX1,ASX2,ASX3
      COMMON/GRID/ N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,
     *             X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3,
     *             NREFL,LN1,LICB,LNFS,LPRED,IND,P,ICB,NFS,TT,IPRED,
     *             IPOSQ,MINQ,MAXQ,TTMIN,DMIN,ASX1,ASX2,ASX3
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(ICB(IA00).NE.ICBREF) THEN
            IA00=0
            GO TO 11
          END IF
        END IF
        IF(P(IA00).GE.GIANT) THEN
          IA00=0
        END IF
      END IF
   11 CONTINUE
      IF(IA20.GT.0) THEN
        IF(LICB) THEN
          IF(ICB(IA20).NE.ICBREF) THEN
            IA20=0
            GO TO 12
          END IF
        END IF
        IF(P(IA20).GE.GIANT) THEN
          IA20=0
        END IF
      END IF
   12 CONTINUE
      IF(IA02.GT.0) THEN
        IF(LICB) THEN
          IF(ICB(IA02).NE.ICBREF) THEN
            IA02=0
            GO TO 13
          END IF
        END IF
        IF(P(IA02).GE.GIANT) THEN
          IA02=0
        END IF
      END IF
   13 CONTINUE
      IF(IA22.GT.0) THEN
        IF(LICB) THEN
          IF(ICB(IA22).NE.ICBREF) THEN
            IA22=0
            GO TO 14
          END IF
        END IF
        IF(P(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=(P(IA00)-P(IA20)-P(IA02)+P(IA22))/4.
            ELSE
              IA00=INDX(IPOS0-IPOS1)
              IA20=INDX(IPOS0+IPOS1)
              DERMIX=(P(IA00)-P(IA20)-P(IA02)+P(IA22))/2.
            END IF
          ELSE
            IF(IA00.GT.0) THEN
              IA20=INDX(IPOS0-IPOS2)
              IA22=INDX(IPOS0+IPOS2)
              DERMIX=(P(IA00)-P(IA20)-P(IA02)+P(IA22))/2.
            ELSE
              IA00=INDX(IPOS0-IPOS1)
              IA20=INDX(IPOS0+IPOS1)
              DERMIX=(P(IA00)-P(IA20)-P(IA02)+P(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=(P(IA00)-P(IA20)-P(IA02)+P(IA22))/2.
            ELSE
              IA00=INDX(IPOS0-IPOS2)
              IA02=INDX(IPOS0+IPOS2)
              DERMIX=(P(IA00)-P(IA20)-P(IA02)+P(IA22))/2.
            END IF
          ELSE
            IF(IA00.GT.0) THEN
              IA02=INDX(IPOS0-IPOS1)
              IA20=INDX(IPOS0-IPOS2)
              IA22=INDX(IPOS0)
              DERMIX= P(IA00)-P(IA20)-P(IA02)+P(IA22)
            ELSE
              IA00=INDX(IPOS0)
              IA20=INDX(IPOS0+IPOS1)
              IA02=INDX(IPOS0+IPOS2)
              DERMIX= P(IA00)-P(IA20)-P(IA02)+P(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=(P(IA00)-P(IA20)-P(IA02)+P(IA22))/2.
            ELSE
              IA00=INDX(IPOS0-IPOS1)
              IA20=INDX(IPOS0)
              IA22=INDX(IPOS0+IPOS2)
              DERMIX= P(IA00)-P(IA20)-P(IA02)+P(IA22)
            END IF
          ELSE
            IF(IA00.GT.0) THEN
              IA20=INDX(IPOS0-IPOS2)
              IA22=INDX(IPOS0+IPOS2)
              DERMIX=(P(IA00)-P(IA20)-P(IA02)+P(IA22))/2.
            ELSE
              IA00=INDX(IPOS0-IPOS1)
              IA20=INDX(IPOS0)
              IA22=INDX(IPOS0+IPOS2)
              DERMIX= P(IA00)-P(IA20)-P(IA02)+P(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=(P(IA00)-P(IA20)-P(IA02)+P(IA22))/2.
            ELSE
              IA00=INDX(IPOS0-IPOS2)
              IA02=INDX(IPOS0)
              IA22=INDX(IPOS0+IPOS1)
              DERMIX= P(IA00)-P(IA20)-P(IA02)+P(IA22)
            END IF
          ELSE
            IF(IA00.GT.0) THEN
              IA02=INDX(IPOS0-IPOS1)
              IA20=INDX(IPOS0-IPOS2)
              IA22=INDX(IPOS0)
              DERMIX= P(IA00)-P(IA20)-P(IA02)+P(IA22)
            ELSE
              IERR=IERR+1
            END IF
          END IF
        END IF
      END IF
C
      RETURN
      END
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 LOCTED
C             IN THE COMPUTATIONAL VOLUME, OTHERWISE ZERO.
C     PS...   ITERPOLATED SLOWNESS.
C
C-----------------------------------------------------------------------
C
      INTEGER MIND,MGRID,MICB,MNFS,MPRED
      PARAMETER (MIND=1,MGRID=100000)
      PARAMETER (MICB=MGRID,MNFS=MGRID,MPRED=MGRID)
      LOGICAL LN1,LICB,LNFS,LPRED
      INTEGER N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,NREFL
      INTEGER IND(MIND),ICB(MICB),NFS(MNFS),IPRED(MPRED)
      INTEGER IPOSQ(MGRID),MINQ,MAXQ
      REAL X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3
      REAL P(MGRID),TT(MGRID),TTMIN,DMIN,ASX1,ASX2,ASX3
      COMMON/GRID/ N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,
     *             X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3,
     *             NREFL,LN1,LICB,LNFS,LPRED,IND,P,ICB,NFS,TT,IPRED,
     *             IPOSQ,MINQ,MAXQ,TTMIN,DMIN,ASX1,ASX2,ASX3
C
C-----------------------------------------------------------------------
C
      EXTERNAL 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(P(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=ICB(IADR)
                  END IF
                ELSE
C
C                 CHECK FOR THE (GEOLOGICAL) BLOCK:
                  IF(LICB) THEN
                    IF(ICB(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=P(IAUX)
                  RETURN
                END IF
                PS=PS+P(IAUX)/AUX
                DS=DS+1.0/AUX
C
              ENDIF
            ENDIF
   20       CONTINUE
   21     CONTINUE
   22   CONTINUE
   23 CONTINUE
      IF(DS.EQ.0.) THEN
        PAUSE 'ERROR 45: SOURCE OR RECEIVER POINT IN A FREE SPACE'
        STOP
      ENDIF
      PS=PS/DS
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE POS(X1,X2,X3,DPOS1,DPOS2,DPOS3,
     *                                      IPOS1,IPOS2,IPOS3,IPOS,IADR)
C
C AUXILIARY SUBROUTINE TO SLOW, TO FIND THE NEAREST GRIDPONT.
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 LOCTED
C             IN THE COMPUTATIONAL VOLUME, OTHERWISE ZERO.
C
C-----------------------------------------------------------------------
C
      INTEGER MIND,MGRID,MICB,MNFS,MPRED
      PARAMETER (MIND=1,MGRID=100000)
      PARAMETER (MICB=MGRID,MNFS=MGRID,MPRED=MGRID)
      LOGICAL LN1,LICB,LNFS,LPRED
      INTEGER N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,NREFL
      INTEGER IND(MIND),ICB(MICB),NFS(MNFS),IPRED(MPRED)
      INTEGER IPOSQ(MGRID),MINQ,MAXQ
      REAL X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3
      REAL P(MGRID),TT(MGRID),TTMIN,DMIN,ASX1,ASX2,ASX3
      COMMON/GRID/ N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,
     *             X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3,
     *             NREFL,LN1,LICB,LNFS,LPRED,IND,P,ICB,NFS,TT,IPRED,
     *             IPOSQ,MINQ,MAXQ,TTMIN,DMIN,ASX1,ASX2,ASX3
C
C-----------------------------------------------------------------------
C
      EXTERNAL 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
          PAUSE 'ERROR 46: SOURCE OR RECEIVER POINT OUTSIDE THE MODEL'
          STOP
        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
            PAUSE 'ERROR 46: SOURCE OR RECEIVER POINT OUTSIDE THE MODEL'
            STOP
          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
            PAUSE 'ERROR 46: SOURCE OR RECEIVER POINT OUTSIDE THE MODEL'
            STOP
          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
          PAUSE 'ERROR 46: SOURCE OR RECEIVER POINT OUTSIDE THE MODEL'
          STOP
        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
      INTEGER FUNCTION INDX(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
      INTEGER MIND,MGRID,MICB,MNFS,MPRED
      PARAMETER (MIND=1,MGRID=100000)
      PARAMETER (MICB=MGRID,MNFS=MGRID,MPRED=MGRID)
      LOGICAL LN1,LICB,LNFS,LPRED
      INTEGER N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,NREFL
      INTEGER IND(MIND),ICB(MICB),NFS(MNFS),IPRED(MPRED)
      INTEGER IPOSQ(MGRID),MINQ,MAXQ
      REAL X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3
      REAL P(MGRID),TT(MGRID),TTMIN,DMIN,ASX1,ASX2,ASX3
      COMMON/GRID/ N1,N2,N3,L1,L2,L3,L4,L1234,NL1,NL2,NL3,
     *             X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,DSX1,DSX2,DSX3,
     *             NREFL,LN1,LICB,LNFS,LPRED,IND,P,ICB,NFS,TT,IPRED,
     *             IPOSQ,MINQ,MAXQ,TTMIN,DMIN,ASX1,ASX2,ASX3
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=IND(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
