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.20
C Date: 2000, May 29
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.......................................................................
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) 'SEP','NET1','NET2','NET3',/
C     'SEP'...String in apostrophes containing the name
C             the input file with the data specifying grid
C             File SEP, specifying the grid dimensions,
C             will be updated by appending the dimensions of the grid
C             for the next network ray tracing.
C             If the same file SEP is used in all iterations, it
C             accumulates the history of the grid dimensions.
C             dimensions and optionally some numerical parameters.
C             Description of file SEP
C             Additional parameters for this program
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 index
C             file IND from this file.
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     Defaults: 'SEP'='net.h', 'NET1'='net1.dat', 'NET2'='net2.dat',
C             'NET3'='net3.dat'.
C
C                                                     
C Data file SEP has the form of the SEP (Stanford Exploration Project)
C parameter file:
C Parameters common with program 'net.for'
C Additional parameters for this program:
C     L1MAX=integer, L2MAX=integer, L3MAX=integer... Output big brick
C             may have at most L1MAX*L2MAX*L3MAX small bricks,
C             0 means no limitation.
C             Do not specify different values of LiMAX in different
C             directions.
C             Defaults: L1MAX=2, L2MAX=2, L3MAX=2.
C     NL1MAX=integer, NL2MAX=integer, NL3MAX=integer... Maximum values
C             of N1*L1, N2*L2, N3*L3, respectively,
C             0 means no limitation.
C             These parameters enable to limit the density of the
C             gridpoints in order to stop calculations at a desired
C             accuracy level and prevent extremely expensive
C             calculation.  It is recommended to specify NL1MAX,
C             NL2MAX, NL3MAX as the same multiple of N1, N2, N3,
C             respectively.
C             Hint:
C               (a) Run 'net.for' on a reasonably dense grid.
C               (b) Look at the maximum errors of calculated travel
C                   times and estimate how many times they should be
C                   decreased.  Let us denote the value by TIMES here.
C               (c) For network ray tracing (NFSMAX.GE.0), which is
C                   a first-order method, select NLiMAX=TIMES*Ni,
C                   where N1,N2,N3 are the values used at (a).
C               (d) Leave default values of LiMAX.
C               (e) Approximate NLiMAX by new values
C                     NLiMAX=Ni*(LiMAX)**(ITER-1),
C                   i.e.,
C                     NLiMAX=Ni*2**(ITER-1) for default LiMAX=2,
C                   with reasonable values of N1, N2, N3.  If resulting
C                   values of N1, N2, N3, NL1MAX, NL2MAX, NL3MAX are
C                   compatible with the RAM allocated in
C                   ram.inc (see
C                   net.inc), they may be used
C                   to start ITER iterations of ray tracing within
C                   Fresnel volumes.
C             If the accuracy is low and you wish to maximize it at
C               given RAM, replace the above procedure by:
C               (a) Leave default values of LiMAX and NLiMAX.
C               (b) Choose initial values of N1,N2,N3 which satisfy
C                     N1*N2*N3=MRAM/4/LiMAX**(2*ITER-4)
C                   in 2-D, or
C                     N1*N2*N3=MRAM/5/LiMAX**(3*ITER-6)
C                   in 3-D, where MRAM is declared in
C                   ram.inc.
C                   Number ITER of iterations is now just an estimation.
C               (c) If ITER.LE.2, increase ITER by decreasing Ni, or
C                   choose maximum initial values of N1,N2,N3 which fit
C                   in the memory according to the description in
C                   net.inc).
C             For more details, refer to the considerations in
C               netfv.htm.
C             Defaults: NL1MAX=0, NL2MAX=0, NL3MAX=0.
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) '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 (2) NREFL,/
C     NREFL...Number of reflections.
C     Default: NREFL=0.
C     NET1, NET2 and NET3 must have the same NREFL.
C (3) Once (3.1), then NREFL-times (3.2) and (3.1):
C (3.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 SEP file.
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 (3.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
      INCLUDE 'ram.inc'
C     ram.inc
      INTEGER IRAM(MRAM)
      EQUIVALENCE (RAM,IRAM)
C
C.......................................................................
C
      CHARACTER*80 FSEP,FNET1,FNET2,FNET3,FSRC,FEND,FIND,FTT1,FTT2,FOUT
      CHARACTER*80 FRAYS,FEND3
      CHARACTER*80 FVEL1,FICB1,FIND2,FVEL2,FICB2,FERR,FPRED,FNFS
      CHARACTER*60 LINE
      CHARACTER*1  FAUX
      INTEGER LU1,LU2
      PARAMETER (LU1=1,LU2=2)
      INTEGER MSMALL,NREFL,IREFL
      INTEGER N1,N2,N3,L1,L2,L3,L4,L1234,NBIG,IBIG,NPOS,IPOS,IADR
      INTEGER L1MAX,L2MAX,L3MAX
      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,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     MSMALL..Maximum number of the small bricks within the Fresnel
C             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     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     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:
      FSEP= 'net.h'
      FNET1='net1.dat'
      FNET2='net2.dat'
      FNET3='net3.dat'
      WRITE(*,'(2A)')  '+NETIND: Enter 4 input filenames',
     *           ' [''net.h'' ''net1.dat'' ''net2.dat'' ''net3.dat'']: '
      READ(*,*) FSEP,FNET1,FNET2,FNET3
C     FSEP,FNET1,FNET2,FNET3 are input/output data files.
C
C.......................................................................
C
C     Loop over reflections
      IREFL=0
   10 CONTINUE
C
C.......................................................................
C
C     Reading input SEP parameter file:
      CALL RSEP1(LU1,FSEP)
C
      CALL RSEP3I('NFSMAX',NFSMAX,0)
      IF(NFSMAX.LT.0) THEN
C       NETIND-03
        CALL ERROR
     *           ('NETIND-03: Second-order method is not supported now')
C       Present coding of the second-order method (NFSMAX.LT.0) does not
C       trace rays and does not include error estimation, which prevents
C       determination of Fresnel volumes.
      END IF
C
C     Numbers of gridpoints:
      CALL RSEP3I('N1',N1,1)
      CALL RSEP3I('N2',N2,1)
      CALL RSEP3I('N3',N3,1)
      CALL RSEP3I('L1',L1,1)
      CALL RSEP3I('L2',L2,1)
      CALL RSEP3I('L3',L3,1)
      IF(N1.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
        CALL ERROR('NETIND-09: Number of gridpoints is not positive')
      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.MRAM) THEN
C       NETIND-10
        CALL ERROR
     *      ('NETIND-10: Too many gridpoints to be stored in array IND')
      END IF
C
C     Maximum numbers of small bricks inside a big brick:
      CALL RSEP3I('L1MAX',L1MAX,2)
      CALL RSEP3I('L2MAX',L2MAX,2)
      CALL RSEP3I('L3MAX',L3MAX,2)
C     Output big brick may have at most L1MAX*L2MAX*L3MAX small bricks.
C     (0 means no limitation)
C
C     Maximum density of the gridpoints:
      CALL RSEP3I('NL1MAX',NL1MAX,0)
      CALL RSEP3I('NL2MAX',NL2MAX,0)
      CALL RSEP3I('NL3MAX',NL3MAX,0)
C
C     Reading the 1-st input file NET1 for the NET program:
      OPEN(LU1,FILE=FNET1,STATUS='OLD')
C     (1) names of the files with source, receivers, rays, and errors:
      FEND=' '
      READ(LU1,*) FSRC,FAUX,FAUX,FEND
      IF(FEND.EQ.' ') THEN
C       NETIND-01
        CALL ERROR
     *       ('NETIND-01: Name of file with times at receivers missing')
      END IF
C     (2) number of reflections:
      NREFL=0
      READ(LU1,*) NREFL
C     (3) 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
          CALL ERROR('NETIND-02: No index file specified')
        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) Names of the files with source, receivers, rays, and errors:
      READ(LU1,*) FAUX,FAUX,FAUX,FAUX
C     (2) Number of reflections:
      I=0
      READ(LU1,*) I
      IF(I.NE.NREFL) THEN
C       NETIND-04
        CALL ERROR('NETIND-04: Different number of reflections in NET2')
      END IF
C     (3) 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
        CALL ERROR('NETIND-05: Different input index files')
      END IF
      IF(FVEL1.NE.FVEL2) THEN
C       NETIND-06
        CALL ERROR('NETIND-06: Different velocity files')
      END IF
      IF(FICB1.NE.FICB2) THEN
C       NETIND-07
        CALL ERROR('NETIND-07: Different block files')
      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) Names of the files with source, receivers, rays, and errors:
      FRAYS=' '
      FEND3=' '
      READ(LU1,*) FAUX,FAUX,FRAYS,FEND3
C     (2) Number of reflections:
      I=0
      READ(LU1,*) I
      IF(I.NE.NREFL) THEN
C       NETIND-08
        CALL ERROR('NETIND-08: Different number of reflections in NET3')
      END IF
C     (3) 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=' '
        FICB1=' '
        FERR=' '
        FPRED=' '
        FNFS=' '
        READ(LU1,*) FOUT,FAUX,FICB1,FAUX,FERR,FPRED,FNFS
   13 CONTINUE
      CLOSE(LU1)
C     End of reading the 3-rd main input data file.
C     FOUT is the output index file.
C
C     FSEP 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)') '+NETIND: Reading source file:   ',FSRC(1:48)
      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
C     TTERR=-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)') '+NETIND: Reading endpoint file: ',FEND(1:48)
      OPEN(LU1,FILE=FEND)
      READ(LU1,*) (FAUX,I=1,20)
      AUX5=0.
      READ(LU1,*) FAUX,AUX1,AUX2,AUX3,AUX4,AUX5
      TTMAX=TTMAX+AUX4+AUX5
C     TTERR=TTERR+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)') '+NETIND: Reading input index file: ',FIND(1:45)
      DO 21 IBIG=1,NBIG
        IRAM(IBIG)=IBIG
   21 CONTINUE
      IF(FIND.NE.' ') THEN
        OPEN(LU1,FILE=FIND)
        READ(LU1,*) (IRAM(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(IRAM(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)')
     *  '+NETIND: 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=IRAM(IN1+N1*(IN2+N2*IN3))
                DO 31 IL1=L1-1,0,-1
                  IPOS=IPOS-1
                  IF(IADR.LE.0) THEN
                    IRAM(IPOS)=0
                  ELSE
                    IRAM(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(NPOS+2*L1234.GT.MRAM) THEN
C       NETIND-11
        CALL ERROR
     *      ('NETIND-11: Too many network nodes with given travel time')
      END IF
      WRITE(*,'(2A)') '+NETIND: Reading travel time field: ',FTT1(1:44)
      OPEN(LU1,FILE=FTT1)
      READ(LU1,*) (RAM(IADR),IADR=NPOS+1,NPOS+L1234)
      CLOSE(LU1)
      WRITE(*,'(2A)') '+NETIND: Reading travel time field: ',FTT2(1:44)
      OPEN(LU1,FILE=FTT2)
      READ(LU1,*) (RAM(IADR),IADR=NPOS+L1234+1,NPOS+2*L1234)
      CLOSE(LU1)
C
C     Converting 'index file' into 'Fresnel volume index file':
      WRITE(*,'(A)')
     *  '+NETIND: Labeling the Fresnel volume...                       '
      L4=0
      DO 41 IPOS=1,NPOS
        IADR=IRAM(IPOS)
        IRAM(IPOS)=0
        IF(IADR.GT.0) THEN
          IF(RAM(NPOS+IADR)+RAM(NPOS+L1234+IADR).LE.TTMAX
     *                   .OR.IPOS.EQ.ISRC.OR.IPOS.EQ.IREC) THEN
            L4=L4+1
            IRAM(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)') '+NETIND: Writing output index file: ',FOUT(1:44)
      OPEN(LU1,FILE=FOUT)
      WRITE(LU1,'(10I8)') (IRAM(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
      IF(N1*N2*N3.GT.MRAM) THEN
C       NETIND-51
        CALL WARN
     *      ('NETIND-51: New big bricks are too small to fit in memory')
      END IF
C
C     Calculating next memory requirements of 'net.for':
      IF(NFSMAX.GE.0) THEN
C       Network ray tracing:
        IF(FICB1.EQ.' ') THEN
          M1ICB=0
        ELSE
          M1ICB=1
        END IF
        IF(FPRED.EQ.' '.AND.
     *     FERR.EQ.' '.AND.FRAYS.EQ.' '.AND.FEND3.EQ.' ') THEN
          M1PRED=0
        ELSE
          M1PRED=1
        END IF
        IF(FNFS.EQ.'*') THEN
          M1NFS=-1
        ELSE
          M1NFS=0
        END IF
      ELSE
C       Second-order grid travel-time tracing:
        M1ICB=0
        M1PRED=0
        M1NFS=0
      END IF
C
C     New number of small bricks (L1*L2*L3):
      MSMALL=(MRAM-N1*N2*N3)/(4+M1ICB+M1PRED+M1NFS)
      AUX1=FLOAT(MSMALL/L4)
      IF(N1.EQ.1.OR.N2.EQ.1.OR.N3.EQ.1) THEN
        L0=INT(SQRT(AUX1))
      ELSE
        L0=INT(AUX1**0.333333)
      END IF
      L1=L0
      L2=L0
      L3=L0
      IL1=2
      IL2=2
      IL3=2
      IF(N1.EQ.1) THEN
        L1=1
        IL1=1
      END IF
      IF(N2.EQ.1) THEN
        L2=1
        IL2=1
      END IF
      IF(N3.EQ.1) THEN
        L3=1
        IL3=1
      END IF
      IF((N1*N2*N3+(4+M1ICB+M1PRED+M1NFS)*L4)*IL1*IL2*IL3.LE.MRAM) THEN
C       L1MAX, L2MAX and L3MAX do not apply to the last iteration
        IF(L1MAX.GT.0) THEN
          L1=MIN0(L1,L1MAX)
        END IF
        IF(L2MAX.GT.0) THEN
          L2=MIN0(L2,L2MAX)
        END IF
        IF(L3MAX.GT.0) THEN
          L3=MIN0(L3,L3MAX)
        END IF
      END IF
      IF(NL1MAX.GT.0) THEN
        L1=MIN0(L1,NL1MAX/N1)
      END IF
      IF(NL2MAX.GT.0) THEN
        L2=MIN0(L2,NL2MAX/N2)
      END IF
      IF(NL3MAX.GT.0) THEN
        L3=MIN0(L3,NL3MAX/N3)
      END IF
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 SEP history file:
      OPEN(LU1,FILE=FSEP)
C     Searching for the end of file
   90 CONTINUE
        READ(LU1,'(A)',END=91)
      GO TO 90
   91 CONTINUE
C     Appending new grid dimensions
      WRITE(LU1,'(A)') '# netind:'
      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
      WRITE(LU1,'(A)') '# netind.'
      WRITE(LU1,'(A)')
      CLOSE(LU1)
C
C     Screen output:
      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',MSMALL,' small bricks'
      WRITE(*,'(A)') ' in Fresnel volume.'
C
C     End of computation:
      IF((N1.GT.1.AND.NL1MAX.GT.0.AND.2*N1*L1.GT.NL1MAX).OR.
     *   (N2.GT.1.AND.NL2MAX.GT.0.AND.2*N2*L2.GT.NL2MAX).OR.
     *   (N3.GT.1.AND.NL3MAX.GT.0.AND.2*N3*L3.GT.NL3MAX)) THEN
        WRITE(*,'(2A)') '+NETIND: *** One more iteration only -',
     *                  ' smaller bricks are not permitted ***'
      ELSE IF(N1*N2*N3*L1*L2*L3.GT.MRAM) THEN
        WRITE(*,'(2A)') '+NETIND: *** One more iteration only -',
     *                  ' more big bricks cannot fit in RAM ***'
      END IF
      IF(L0.LE.1) THEN
C       NETIND-52
        CALL WARN
     *     ('NETIND-52: Big bricks cannot be divided into small bricks')
      ELSE IF(L1.LE.1.AND.L2.LE.1.AND.L3.LE.1) THEN
C       NETIND-53
        CALL WARN
     *   ('NETIND-53: Big bricks will not be divided into small bricks')
      END IF
      WRITE(*,'(A)')
     *  ' NETIND: Done.                                                '
      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
          CALL ERROR
     *         ('NETIND-12: Source or receiver point outside the model')
        ELSE IF(IX.GE.NLX) THEN
          IX=NLX-1
        END IF
      END IF
      RETURN
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'length.for'
C     length.for
*     INCLUDE 'forms.for'
C     forms.for
C
C=======================================================================
C