C 
C Program NETIND to generate the index file mapping gridpoints onto the
C network nodes situated within the Fresnel volume.
C
C Version: 3.00
C Date: 1997, October 24
C
C Coded by: Ludek Klimes
C     Department of Geophysics, Charles University Prague
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic
C     E-mail: klimes@seis.karlov.mff.cuni.cz
C
C Note:
C     In most cases, it is reasonable to declare dimensions MPOS and MTT
C     of arrays IND(MPOS), TT1(MTT), and TT2(MTT) equal to dimensions
C     MIND and MGRID of program 'net.for', respectively.  Then program
C     'netind.for' is able to reasonably control the division of big
C     bricks into small bricks at all iterations of network two-point
C     ray tracing within the Fresnel volumes.
C
C.......................................................................
C
C                                                    
C Data files:
C
C Main input data read from the * device:
C     The data are read in by the list directed input (free format).
C     The strings have to be enclosed in apostrophes.
C (1) 'NET1','NET2','NET3',MIND,MGRID,/
C     'NET1'..String containing the name of the input file NET of
C             the 'net.for' program when it calculated the travel times
C             from the source point.
C     'NET2'..String containing the name of the input file NET of
C             the 'net.for' program when it calculated the travel times
C             from the receiver point.  Otherwise, file NET2 should be
C             similar to file NET1.
C     'NET3'..String containing the name of the input file NET of
C             the 'net.for' program to perform network ray tracing in
C             the Fresnel volume.  Program NETIND reads only the names
C             of file SEP, specifying the grid dimensions, and index
C             file IND from this file.  File SEP specifying the grid
C             dimensions will be updated by appending the dimensions
C             of the grid for the next network ray tracing.
C             Index file IND is the output of program NETIND.
C             Filename NET3 may coincide with NET1.
C     Description of input data NET1, NET2 and NET3
C     MIND... Zero, or the maximum number of big bricks within the
C             whole model volume.  See the array dimension MIND of the
C             program 'net.for'.  If nonzero, the division of big bricks
C             into small bricks is controlled in such a way that the
C             number of small bricks within the whole model volume (i.e.
C             The number of big bricks in the next iteration) does not
C             exceed MIND.  This limitation is not applied for the last
C             iteration (i.e. if the number of big bricks cannot be
C             increased within the limit of MIND).
C     MGRID.. Zero, or the maximum number of the small bricks within the
C             Fresnel volume.  See the array dimension MGRID of the
C             program 'net.for'.  If nonzero, the division of big bricks
C             into small bricks is controlled in such a way that the
C             number of small bricks within the Fresnel volume does not
C             exceed MGRID.
C     Default: 'NET1'='net1.dat', 'NET2'='net2.dat', 'NET3'='net3.dat',
C             MIND=MPOS, MGRID=MTT,
C     where:
C     MPOS... Is the dimension of array IND in this program, i.e. the
C             maximum number of input small bricks = output big bricks
C             in the whole model volume.  If using the default, it is
C             assumed that the arrays IND(MPOS) in 'netind.for' and
C             IND(MIND) in 'net.for' are of the same length.
C     MTT...  Is the dimension of arrays TT1 and TT2 in this program,
C             i.e. the maximum number of input small bricks in the input
C             computational (Fresnel) volume.  If using the default, it
C             is assumed that the arrays TT1(MTT),TT2(MTT) in
C             'netind.for' and TT(MGRID) in 'net.for' are of the same
C             length.
C
C                                                     
C Structure of input data files NET1, NET2, NET3:
C     Sequential files, read by list directed (free format) input,
C     containing model parameters, source/receiver coordinates, and
C     names of other input and output files for the 'net.for' program.
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) 'SEP1','SEP2',/
C     'SEP1','SEP2'... String(s) in apostrophes containing the name(s)
C             of one or two input file(s) with the data specifying grid
C             dimensions and optionally some numerical parameters.
C             Usually, only 'SEP1' is specified.
C             Description of file SEP
C     Default: 'SEP1'='net.h', 'SEP2'=' '.
C     NET1 and NET2 must refer to the same files SEP1 and SEP2.
C     NET3:   File SEP (SEP1 or SEP2), specifying the grid dimensions,
C             will be updated by appending the dimensions of the grid
C             for the next network ray tracing.  If both filenames
C             'SEP1' and 'SEP2' are nonblank, file SEP2 is updated.
C             If the same file SEP is used in all iterations, it
C             accumulates the history of the grid dimensions.
C (2) 'SRC','REC','RAYS','END',/
C     'SRC'...Name of the input file with source coordinates.
C             Description of input data SRC
C     'REC'...Name of the input file with receiver coordinates.
C             If blank, no rays are stored within the file 'RAYS'.
C             Description of input data REC
C     'RAYS'..Name of the output file with rays.
C             If blank, no rays are stored within the file 'RAYS'.
C             Description of data RAYS
C     'END'...Name of the output file with endpoints of rays (receiver
C             coordinates, receiver arrival times, and estimates of the
C             corresponding maximum travel-time errors.
C             If blank, no file 'END' is generated.
C             Description of data END
C     Default: 'REC'=' ', 'RAYS'=' ', 'END'=' '.
C     NET1:   Files SRC and END are input files for program NETIND.
C     NET2:   Files SRC and REC from NET1 must be swopped.
C     NET3:   This line is skipped.
C (3) NREFL,/
C     NREFL...Number of reflections.
C     Default: NREFL=0.
C     NET1, NET2 and NET3 must have the same NREFL.
C (4) Once (4.1), then NREFL-times (4.2) and (4.1):
C (4.1) 'IND(I)','VEL(I)','ICB(I)','TT(I)','ERR(I)','PRED(I)','NFS(I)',/
C     'IND(I)'... Name of the index file, specifying for each
C             big brick if its gridpoints belong to the network.
C             If it is blank, 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             Description of data IND(I)
C     'VEL(I)'... Name of the input file containing velocities at all
C             network nodes, for I-times reflected wave.
C             Has always to be specified.
C             Description of data VEL(I)
C     'ICB(I)'... Name of the input file containing indices of
C             (geological) blocks.  For more detail refer to the
C             description of this item in program
C             'net.for'.
C             Description of data ICB(I)
C     'TT(I)'... Name of the file containing travel-times at all
C             network nodes after I reflections.
C             Description of data TT(I)
C     'ERR(I)'... Name of the output file containing estimated upper
C             bounds for the errors of the computed travel-times at all
C             network nodes after I reflections.
C             Description of data ERR(I)
C     'PRED(I)'... Name of the file containing predecessors of
C             all network nodes after I reflections.
C             May be blank for most applications.
C             Description of data PRED(I)
C     'NFS(I)'... Unimportant file.  Refer to the description in
C             'net.for'.
C     Default: 'IND(I)'=' ', 'VEL(I)'=' ', 'ICB(I)'=' ',
C             'TT(I)'=' ','ERR(I)'=' ', 'PRED(I)'=' ', 'NFS(I)'=' '.
C     NET1 and NET2: Files IND(I), VEL(I) and ICB(I) must be the same
C             for each I.  Files IND(I) may be blank.
C     NET1:   Files TT(I) are the input files for program NETIND.
C     NET2:   Files TT(I) are the input files for program NETIND.
C     NET3:   Files IND(I) are the output files of program NETIND.
C (4.2) 'INTF(I)',/
C     'INTF(I)'... Name of the input file containing refractor points.
C             Description of data INTF(I)
C     NET1, NET2 and NET3: This line is skipped.
C Example of data set NET1
C Example of data set NET2
C Example of data set NET3
C
C-----------------------------------------------------------------------
C
      PROGRAM NETIND
C
C.......................................................................
C
      CHARACTER*80 FNET1,FNET2,FNET3,FSRC,FEND,FIND,FTT1,FTT2,FOUT
      CHARACTER*80 SEP1,SEP2,SEP3,SEP4
      CHARACTER*80 FVEL1,FICB1,FIND2,FVEL2,FICB2
      CHARACTER*60 LINE
      CHARACTER*1  FAUX
      INTEGER LU1,LU2,MPOS,MTT
      PARAMETER (LU1=1,LU2=2,MPOS=150000,MTT=100000)
      INTEGER MIND,MGRID,NREFL,IREFL
      INTEGER N1,N2,N3,L1,L2,L3,L4,L1234,NBIG,IBIG,NPOS,IPOS,IADR
      INTEGER L1MAX,L2MAX,L3MAX,IND(MPOS)
      INTEGER ISRC,ISRC1,ISRC2,ISRC3,IREC,IREC1,IREC2,IREC3
      INTEGER IN1,IN2,IN3,IL1,IL2,IL3,I,J
      REAL D1,D2,D3,O1,O2,O3
      REAL X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,TT1(MTT),TT2(MTT),TTMAX
      REAL AUX1,AUX2,AUX3,AUX4,AUX5,AUX6
C
C     FNET1,FNET2,FNET3... Main input and output files.
C     FSRC,FEND,FIND,FTT1,FTT2,FOUT... Other input and output files.
C     FVEL1,FICB1,FIND2,FVEL2,FICB2,FAUX... Temporary filenames or text
C            strings.
C     LU1,LU2... Input-output logical unit numbers used for different
C             files.
C     MPOS... Maximum number of input small bricks = output big bricks
C             in the whole model volume.
C     MTT...  Maximum number of input small bricks in the input
C             computational (Fresnel) volume.
C     MIND... Zero, or the maximum number of big bricks within the whole
C             model volume.
C     MGRID.. Zero, or the maximum number of the small bricks within the
C             Fresnel volume.
C     NREFL...Number of reflections.  NREFL=0 for a refracted wave.
C     IREFL...Loop variable over reflections.  IREFL=0 for a refracted
C             wave.
C     N1,N2,N3... Numbers of big bricks along gridlines.
C     L1,L2,L3... Numbers small bricks within a big brick.
C     L4...   Input:  Number of big bricks belonging to the network,
C               i.e. length of the travel-time files.
C             Output:  Number of small bricks belonging to the Fresnel
C               volume.
C     L1234...L1*L2*L3*L4 for input values.
C     NBIG... Number of big bricks, i.e. length of the input index file,
C             NBIG=N1*N2*N3.
C     IBIG... Index of a big brick (IBIG=1,2,...,NBIG).
C     NPOS... Number of small bricks, i.e. length of the output index
C             file:  NPOS=N1*N2*N3*L1*L2*L3.
C     IPOS... Index of a small brick (IPOS=1,2,...,NPOS).
C     IADR... Index within a travel time file or within a Fresnel volume
C             (IADR=1,2,3,...,L4 or IADR=0).
C     L1MAX,L2MAX,L3MAX... Maximum numbers of output small bricks in an
C             output big brick.
C     IND...  Array to store an index file.
C     ISRC,ISRC1,ISRC2,ISRC3,IREC,IREC1,IREC2,IREC3... Positions of the
C             source and receiver small bricks.
C     IN1,IN2,IN3,IL1,IL2,IL3,I... Loop and temporary variables.
C     X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX ... Boundaries of the model
C             volume.
C     TT1,TT2... Arrays to store travel times.
C     TTMAX...Maximum sum of travel times, limiting the Fresnel volume.
C     AUX1,AUX2,AUX3,AUX4,AUX5,AUX6... Temporary storage locations.
C
C.......................................................................
C
C     Reading main input data from the * external unit:
      FNET1='net1.dat'
      FNET2='net2.dat'
      FNET3='net3.dat'
      MIND=MPOS
      MGRID=MTT
      WRITE(*,'(2A)')  ' Enter 3 names of input files for ''NET'', and '
     *                                     ,'max.number of small bricks'
      WRITE(*,'(2A,2(I8,A))')    ' (default ''net1.dat'', ''net2.dat'','
     *              ,' ''net3.dat'',',MIND,',',MGRID,'): '
      READ(*,*) FNET1,FNET2,FNET3,MIND,MGRID
C     FNET1,FNET2,FNET3 are input/output data files.
C     MGRID is maximum number of output small bricks.
C
C.......................................................................
C
C     Loop over reflections
      IREFL=0
   10 CONTINUE
C
C.......................................................................
C
C     Reading the 1-st input file NET1 for the NET program:
      OPEN(LU1,FILE=FNET1,STATUS='OLD')
C     (1) SEP parameter files:
      SEP1='net.h'
      SEP2=' '
      READ(LU1,*) SEP1,SEP2
C     (2) names of the files with source, receivers, rays, and errors:
      FEND=' '
      READ(LU1,*) FSRC,FAUX,FAUX,FEND
      IF(FEND.EQ.' ') THEN
C       NETIND-01
        PAUSE
     *   'Error NETIND-01: Name of file with times at receivers missing'
        STOP
      END IF
C     (3) number of reflections:
      NREFL=0
      READ(LU1,*) NREFL
C     (4) names of the output travel-time and predecessor files,
C     input velocity and index files, and input refractor-point files:
      DO 11 I=0,IREFL
        IF(I.GT.0) THEN
          READ(LU1,*) FAUX
        END IF
        FIND=' '
        FVEL1=' '
        FICB1=' '
        FTT1=' '
        READ(LU1,*) FIND,FVEL1,FICB1,FTT1,FAUX,FAUX,FAUX
   11 CONTINUE
      IF(FIND.EQ.' ') THEN
        IF(L1.GT.1.OR.L2.GT.1.OR.L3.GT.1) THEN
C         NETIND-02
          PAUSE 'Error NETIND-02: No index file specified'
          STOP
        END IF
      END IF
      CLOSE(LU1)
C     End of reading the 1-st main input data file.
C     FSRC and FEND are the source and endpoint (receiver) filenames.
C     NREFL is the number of reflections.
C     FTT1 is the input travel-time file, with times from the source.
C     FIND is the input index file.
C
C     Reading the 2-nd input file NET2 for the NET program:
      OPEN(LU1,FILE=FNET2,STATUS='OLD')
C     (1) SEP parameter files:
      SEP3='net.h'
      SEP4=' '
      READ(LU1,*) SEP3,SEP4
      IF(SEP1.NE.SEP3.OR.SEP2.NE.SEP4) THEN
C       NETIND-03
        PAUSE
     *     'Error NETIND-03: Different files specifying the input grids'
        STOP
      END IF
C     (2) Names of the files with source, receivers, rays, and errors:
      READ(LU1,*) FAUX,FAUX,FAUX,FAUX
C     (3) Number of reflections:
      I=0
      READ(LU1,*) I
      IF(I.NE.NREFL) THEN
C       NETIND-04
        PAUSE 'Error NETIND-04: Different number of reflections in NET2'
        STOP
      END IF
C     (4) Names of the output travel-time and predecessor files,
C     input velocity and index files, and input refractor-point files:
      DO 12 I=0,NREFL-IREFL
        IF(I.GT.0) THEN
          READ(LU1,*) FAUX
        END IF
        FIND2=' '
        FVEL2=' '
        FICB2=' '
        FTT2=' '
        READ(LU1,*) FIND2,FVEL2,FICB2,FTT2,FAUX,FAUX,FAUX
   12 CONTINUE
      IF(FIND.NE.FIND2) THEN
C       NETIND-05
        PAUSE 'Error NETIND-05: Different input index files'
        STOP
      END IF
      IF(FVEL1.NE.FVEL2) THEN
C       NETIND-06
        PAUSE 'Error NETIND-06: Different velocity files'
        STOP
      END IF
      IF(FICB1.NE.FICB2) THEN
C       NETIND-07
        PAUSE 'Error NETIND-07: Different block files'
        STOP
      END IF
      CLOSE(LU1)
C     End of reading the 2-nd main input data file.
C     FTT2 is the input travel-time file, with times from the receiver.
C
C     Reading the 3-rd input file NET3 for the NET program:
      OPEN(LU1,FILE=FNET3,STATUS='OLD')
C     (1) SEP parameter files:
      SEP3='net.h'
      SEP4=' '
      READ(LU1,*) SEP3,SEP4
C     (2) Names of the files with source, receivers, rays, and errors:
      READ(LU1,*) FAUX,FAUX,FAUX,FAUX
C     (3) Number of reflections:
      I=0
      READ(LU1,*) I
      IF(I.NE.NREFL) THEN
C       NETIND-08
        PAUSE 'Error NETIND-08: Different number of reflections in NET3'
        STOP
      END IF
C     (4) Names of the output travel-time and predecessor files,
C     input velocity and index files, and input refractor-point files:
      DO 13 I=0,IREFL
        IF(I.GT.0) THEN
          READ(LU1,*) FAUX
        END IF
        FOUT=' '
        READ(LU1,*) FOUT,FAUX,FAUX,FAUX,FAUX,FAUX,FAUX
   13 CONTINUE
      CLOSE(LU1)
C     End of reading the 3-rd main input data file.
C     FOUT is the output index file.
C
C     Reading input SEP parameter files:
      CALL RSEP1(LU1,SEP1)
      CALL RSEP1(LU1,SEP2)
C     Numbers of gridpoints:
      CALL RSEP3I('N1',N1,1)
      CALL RSEP3I('N2',N2,1)
      CALL RSEP3I('N3',N3,1)
      CALL RSEP3I('L1',L1,1)
      CALL RSEP3I('L2',L2,1)
      CALL RSEP3I('L3',L3,1)
      IF(N1.LT.1.OR.N2.LT.1.OR.N3.LT.1.OR.
     *   L1.LT.1.OR.L2.LT.1.OR.L3.LT.1) THEN
C       NETIND-09
        PAUSE 'Error NETIND-09: Number of gridpoints is not positive'
        STOP
      END IF
C     Boundaries of the model volume:
      CALL RSEP3R('D1',D1,1.)
      CALL RSEP3R('D2',D2,1.)
      CALL RSEP3R('D3',D3,1.)
      CALL RSEP3R('O1',O1,0.)
      CALL RSEP3R('O2',O2,0.)
      CALL RSEP3R('O3',O3,0.)
      X1MIN=O1-0.5*D1
      X2MIN=O2-0.5*D2
      X3MIN=O3-0.5*D3
      X1MAX=X1MIN+FLOAT(N1)*D1
      X2MAX=X2MIN+FLOAT(N2)*D2
      X3MAX=X3MIN+FLOAT(N3)*D3
C     Input grid has N1*N2*N3 big bricks by L1*L2*L3 small bricks.
C     Boundaries of model volume: X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX.
C     Total number of big bricks
      NBIG=N1*N2*N3
C     Total number of small bricks (i.e. of gridpoints)
      NPOS=NBIG*L1*L2*L3
      IF(NPOS.GT.MPOS) THEN
C       NETIND-10
        PAUSE
     *  'Error NETIND-10: Too many gridpoints to be stored in array IND'
        STOP
      END IF
C
C     Maximum numbers of small bricks inside a big brick:
      CALL RSEP1(LU1,SEP3)
      CALL RSEP1(LU1,SEP4)
      CALL RSEP3I('L1MAX',L1MAX,0)
      CALL RSEP3I('L2MAX',L2MAX,0)
      CALL RSEP3I('L3MAX',L3MAX,0)
C     Output big brick may have at most L1MAX*L2MAX*L3MAX small bricks.
C     (0 means no limitation)
      IF(SEP4.EQ.' ') THEN
        SEP4=SEP3
      END IF
C     SEP4 is the name of the SEP parameter file to be updated.
C
C.......................................................................
C
C     Reading coordinates of the source point from 'src':
      WRITE(*,'(2A)') '+Reading source file:   ',FSRC(1:56)
      OPEN(LU1,FILE=FSRC)
      READ(LU1,*) (FAUX,I=1,20)
      TTMAX=0.
      AUX5=0.
      READ(LU1,*) FAUX,AUX1,AUX2,AUX3,TTMAX,AUX5
      TTMAX=TTMAX-AUX5
      CLOSE(LU1)
      AUX4=(X1MAX-X1MIN)/(1000.*FLOAT(N1*L1))
      AUX5=(X2MAX-X2MIN)/(1000.*FLOAT(N2*L2))
      AUX6=(X3MAX-X3MIN)/(1000.*FLOAT(N3*L3))
      CALL POSX(AUX1-AUX4,X1MIN,X1MAX,N1*L1,IN1)
      CALL POSX(AUX1+AUX4,X1MIN,X1MAX,N1*L1,IL1)
      CALL POSX(AUX2-AUX5,X2MIN,X2MAX,N2*L2,IN2)
      CALL POSX(AUX2+AUX5,X2MIN,X2MAX,N2*L2,IL2)
      CALL POSX(AUX3-AUX6,X3MIN,X3MAX,N3*L3,IN3)
      CALL POSX(AUX3+AUX6,X3MIN,X3MAX,N3*L3,IL3)
      ISRC=1+IN1+(IN2+IN3*N2*L2)*N1*L1
      ISRC1= IL1-IN1
      ISRC2=(IL2-IN2)*N1*L1
      ISRC3=(IL3-IN3)*N2*L2*N1*L1
C     Source point is situated in the ISRC-th small brick,
C     or in small bricks shifted by ISRC1 and/or ISRC2 and/or ISRC3.
C
C     Reading coordinates of the receiver point and time from 'END':
      WRITE(*,'(2A)') '+Reading endpoint file: ',FEND(1:56)
      OPEN(LU1,FILE=FEND)
      READ(LU1,*) (FAUX,I=1,20)
      AUX5=0.
      READ(LU1,*) FAUX,AUX1,AUX2,AUX3,AUX4,AUX5
      TTMAX=TTMAX+AUX4+AUX5
      CLOSE(LU1)
      AUX4=(X1MAX-X1MIN)/(1000.*FLOAT(N1*L1))
      AUX5=(X2MAX-X2MIN)/(1000.*FLOAT(N2*L2))
      AUX6=(X3MAX-X3MIN)/(1000.*FLOAT(N3*L3))
      CALL POSX(AUX1-AUX4,X1MIN,X1MAX,N1*L1,IN1)
      CALL POSX(AUX1+AUX4,X1MIN,X1MAX,N1*L1,IL1)
      CALL POSX(AUX2-AUX5,X2MIN,X2MAX,N2*L2,IN2)
      CALL POSX(AUX2+AUX5,X2MIN,X2MAX,N2*L2,IL2)
      CALL POSX(AUX3-AUX6,X3MIN,X3MAX,N3*L3,IN3)
      CALL POSX(AUX3+AUX6,X3MIN,X3MAX,N3*L3,IL3)
      IREC=1+IN1+(IN2+IN3*N2*L2)*N1*L1
      IREC1= IL1-IN1
      IREC2=(IL2-IN2)*N1*L1
      IREC3=(IL3-IN3)*N2*L2*N1*L1
C     Receiver point is situated in the IREC-th small brick,
C     or in small bricks shifted by irec1 and/or IREC2 and/or IREC3.
C     TTMAX is maximum sum of travel times, limiting the Fresnel volume.
C
C     READING INPUT INDEX FILE (DEFAULT: 1,2,3,4,...):
      WRITE(*,'(2A)') '+Reading input index file: ',FIND(1:53)
      DO 21 IBIG=1,NBIG
        IND(IBIG)=IBIG
   21 CONTINUE
      IF(FIND.NE.' ') THEN
        OPEN(LU1,FILE=FIND)
        READ(LU1,*) (IND(IBIG),IBIG=1,NBIG)
        CLOSE(LU1)
      END IF
C
C     Number of bricks covered by the network:
C     Big bricks:
      L4=0
      DO 22 IBIG=1,NBIG
        L4=MAX0(IND(IBIG),L4)
   22 CONTINUE
C     Small bricks (number of travel times to be read in):
      L1234=L1*L2*L3*L4
C
C     Upgrading small bricks to big bricks (updating 'index file'):
      WRITE(*,'(A)')
     *  '+Updating index file.                                         '
      IPOS=NPOS+1
      DO 36 IN3=N3-1,0,-1
        DO 35 IL3=L3-1,0,-1
          DO 34 IN2=N2-1,0,-1
            DO 33 IL2=L2-1,0,-1
              DO 32 IN1=N1,1,-1
                IADR=IND(IN1+N1*(IN2+N2*IN3))
                DO 31 IL1=L1-1,0,-1
                  IPOS=IPOS-1
                  IF(IADR.LE.0) THEN
                    IND(IPOS)=0
                  ELSE
                    IND(IPOS)=IL1+L1*(IL2+L2*(IL3+L3*IADR-L3))+1
                  END IF
   31           CONTINUE
   32         CONTINUE
   33       CONTINUE
   34     CONTINUE
   35   CONTINUE
   36 CONTINUE
C
C     Reading travel times:
      IF(L1234.GT.MTT) THEN
C       NETIND-11
        PAUSE
     *  'Error NETIND-11: Too many network nodes with given travel time'
        STOP
      END IF
      WRITE(*,'(2A)') '+Reading travel time field: ',FTT1(1:52)
      OPEN(LU1,FILE=FTT1)
      READ(LU1,*) (TT1(IADR),IADR=1,L1234)
      CLOSE(LU1)
      WRITE(*,'(2A)') '+Reading travel time field: ',FTT2(1:52)
      OPEN(LU1,FILE=FTT2)
      READ(LU1,*) (TT2(IADR),IADR=1,L1234)
      CLOSE(LU1)
C
C     Converting 'index file' into 'Fresnel volume index file':
      WRITE(*,'(A)')
     *  '+Labeling the Fresnel volume.                                 '
      L4=0
      DO 41 IPOS=1,NPOS
        IADR=IND(IPOS)
        IND(IPOS)=0
        IF(IADR.GT.0) THEN
          IF(TT1(IADR)+TT2(IADR).LE.TTMAX.OR.IPOS.EQ.ISRC
     *                                   .OR.IPOS.EQ.IREC) THEN
            L4=L4+1
            IND(IPOS)=L4
            IF(IPOS.EQ.ISRC) THEN
              IF(ISRC1.LE.0.AND.ISRC2.LE.0) THEN
                ISRC=ISRC+ISRC3
                ISRC3=-ISRC3
              END IF
              IF(ISRC1.LE.0) THEN
                ISRC=ISRC+ISRC2
                ISRC2=-ISRC2
              END IF
              ISRC=ISRC+ISRC1
              ISRC1=-ISRC1
            END IF
            IF(IPOS.EQ.IREC) THEN
              IF(IREC1.LE.0.AND.IREC2.LE.0) THEN
                IREC=IREC+IREC3
                IREC3=-IREC3
              END IF
              IF(IREC1.LE.0) THEN
                IREC=IREC+IREC2
                IREC2=-IREC2
              END IF
              IREC=IREC+IREC1
              IREC1=-IREC1
            END IF
          END IF
        END IF
   41 CONTINUE
C
C     Writing Fresnel volume index file:
      WRITE(*,'(2A)') '+Writing output index file: ',FOUT(1:52)
      OPEN(LU1,FILE=FOUT)
      WRITE(LU1,'(10I8)') (IND(IPOS),IPOS=1,NPOS)
      CLOSE(LU1)
C
C.......................................................................
C
      IREFL=IREFL+1
      IF(IREFL.LE.NREFL) GO TO 10
C     End of loop for reflections
C
C.......................................................................
C
C     New number of big bricks (N1*N2*N3):
      N1=N1*L1
      N2=N2*L2
      N3=N3*L3
C
C     New number of small bricks (L1*L2*L3):
      IF(MGRID.GT.0) THEN
        AUX1=FLOAT(MGRID/L4)
        IF(N1.EQ.1.OR.N2.EQ.1.OR.N3.EQ.1) THEN
          I=INT(SQRT(AUX1))
        ELSE
          I=INT(AUX1**0.333333)
        END IF
        L1=I
        L2=I
        L3=I
        IF(L1MAX.GT.1) THEN
          L1=MIN0(L1,L1MAX)
        END IF
        IF(L2MAX.GT.1) THEN
          L2=MIN0(L2,L2MAX)
        END IF
        IF(L3MAX.GT.1) THEN
          L3=MIN0(L3,L3MAX)
        END IF
      ELSE
        L1=MAX0(2,L1MAX)
        L2=MAX0(2,L2MAX)
        L3=MAX0(2,L3MAX)
        I=2
      END IF
      IF(N1.EQ.1) THEN
        L1=1
      END IF
      IF(N2.EQ.1) THEN
        L2=1
      END IF
      IF(N3.EQ.1) THEN
        L3=1
      END IF
      IF(MIND.GT.0) THEN
        AUX1=FLOAT(MIND/(N1*N2*N3))+0.5
        IF(N1.EQ.1.OR.N2.EQ.1.OR.N3.EQ.1) THEN
          J=INT(SQRT(AUX1))
        ELSE
          J=INT(AUX1**0.333333)
        END IF
        IF(J.GT.1) THEN
C         limiting the number of small bricks that are likely to become
C         big bricks in the next iteration:
          L1=MIN0(L1,J)
          L2=MIN0(L2,J)
          L3=MIN0(L3,J)
        END IF
      ELSE
        J=2
      END IF
C
C     Screen output:
      IF(MGRID.LE.0) THEN
        WRITE(*,'(A,I6,A,I7,A,3(I3,A),3(I7,A))')
     *     '+',L4,' of',N1*N2*N3,' big bricks,',L1,'*',L2,'*',L3,'*',L4,
     *                       '=',L1*L2*L3*L4,' small bricks'
      ELSE
        WRITE(*,'(A,I6,A,I7,A,3(I3,A),3(I7,A))')
     *     '+',L4,' of',N1*N2*N3,' big bricks,',L1,'*',L2,'*',L3,'*',L4,
     *                       '=',L1*L2*L3*L4,' of',MGRID,' small bricks'
      END IF
      WRITE(*,'(A)') ' in Fresnel volume.'
C
C     New grid dimensions:
      D1=(X1MAX-X1MIN)/FLOAT(N1)
      D2=(X2MAX-X2MIN)/FLOAT(N2)
      D3=(X3MAX-X3MIN)/FLOAT(N3)
      O1=X1MIN+0.5*D1
      O2=X2MIN+0.5*D2
      O3=X3MIN+0.5*D3
C
C     Updating 'SEP3':
      OPEN(LU1,FILE=SEP3)
C     Searching for the end of file
   90 CONTINUE
        READ(LU1,'(A)',END=91)
      GO TO 90
   91 CONTINUE
C     Appending new grid dimensions
      CALL WSEPI(LINE( 1:20),'N1',N1)
      CALL WSEPI(LINE(21:40),'N2',N2)
      CALL WSEPI(LINE(41:60),'N3',N3)
      WRITE(LU1,'(A)') LINE
      CALL WSEPI(LINE( 1:20),'L1',L1)
      CALL WSEPI(LINE(21:40),'L2',L2)
      CALL WSEPI(LINE(41:60),'L3',L3)
      WRITE(LU1,'(A)') LINE
      CALL WSEPR(LINE( 1:20),'D1',D1)
      CALL WSEPR(LINE(21:40),'D2',D2)
      CALL WSEPR(LINE(41:60),'D3',D3)
      WRITE(LU1,'(A)') LINE
      CALL WSEPR(LINE( 1:20),'O1',O1)
      CALL WSEPR(LINE(21:40),'O2',O2)
      CALL WSEPR(LINE(41:60),'O3',O3)
      WRITE(LU1,'(A)') LINE
      CLOSE(LU1)
C
C     End of computation:
      IF(J.GT.1) THEN
        WRITE(*,'(A)')
     *  '+End.                                                         '
      ELSE IF(N1*N2*N3.LE.MIND) THEN
        WRITE(*,'(2A)') '+*** One more iteration only -',
     *                  ' big bricks cannot be made smaller ***'
      ELSE
        PAUSE 'Warning: New big bricks are too small to fit in memory.'
      END IF
      IF(I.LE.1) THEN
        PAUSE 'Warning: Big bricks cannot be divided into small bricks.'
      END IF
      STOP
      END
C
C=======================================================================
C
      SUBROUTINE POSX(X,XMIN,XMAX,NLX,IX)
C
C Subroutine determining the grid interval along the axis.
C
C Input:
C     X...    A coordinate of a given point.
C     XMIN,XMAX... Limits of the grid line.
C     NLX...  The grid line is divided into n1*l1 grid intervals.
C
C Output:
C     IX...   The given point lies in the ix-th grid interval.
C
C Date: 1993, October 18
C coded by: Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     No auxiliary storage locations.
C
      IF(NLX.EQ.1) THEN
        IX=0
      ELSE
        IX=INT(FLOAT(NLX)*(X-XMIN)/(XMAX-XMIN))
        IF(IX.LT.0.OR.NLX.LT.IX) THEN
C         NETIND-12
          PAUSE
     *     'Error NETIND-12: Source or receiver point outside the model'
          STOP
        ELSE IF(IX.GE.NLX) THEN
          IX=NLX-1
        END IF
      END IF
      RETURN
      END
C
C=======================================================================
C
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'length.for'
C     length.for
*     INCLUDE 'forms.for'
C     forms.for
C
C=======================================================================
C