C PROGRAM 'NETIND' TO GENERATE THE INDEX FILE MAPPING GRIDPOINTS ONTO
C THE NETWORK NODES SITUATED WITHIN THE FRESNEL VOLUME.
C
C VERSION: 2.00
C DATE: 1994, FEBRUARY 16
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: PSENCIK@CSEARN.BITNET, OR PSENCIK@EARN.CVUT.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 MAIN INPUT DATA READ FROM THE * DEVICE:
C     THE DATA ARE READ IN BY THE LIST DIRECTED INPUT (FREE FORMAT).
C     IN THE LIST OF INPUT DATA BELOW, EACH NUMBERED PARAGRAPH INDICATES
C     THE BEGINNING OF A NEW INPUT OPERATION (NEW READ STATEMENT).
C (1) 'NET1','NET2','NET3','NET4',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.
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.  THE FIRST LINE OF THE FILE WILL BE
C             UPDATED ACCORDING TO THE INPUT DATA AND THE SIZE OF THE
C             FRESNEL VOLUME.  OTHER LINES OF THE FILE WILL REMAIN
C             UNCHANGED.  THIS FILENAME MAY COINCIDE WITH 'NET1'.
C     'NET4'..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, OR A BLANK.  IF 'NET4' IS NOT BLANK
C             THE FILE WILL BE UPDATED LIKE 'NET3'.  THIS FILENAME MAY
C             COINCIDE WITH 'NET2'.
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',
C             'NET3'='NET3.OUT', 'NET4'='NET4.OUT', 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 STRUCTURE OF THE DATA FILES 'NET1', 'NET2', 'NET3', 'NET4:
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) 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 (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 OTHER DATA FILES REQUIRED:
C     'SRC'...COORDINATES OF THE SOURCE POINT.  THE NAME READ FROM FILE
C             'NET1', DATA (5), ITEM 'SRC'.  ONLY THE FIRST SOURCE
C             POINT IS BEING READ, THE FILE SHOULD NOT CONTAIN MORE
C             SOURCE POINTS.
C     'END'...COORDINATES OF THE RECEIVER POINT, TRAVEL TIME, AND ITS
C             ESTIMATED MAXIMUM ERROR.  THE NAME READ FROM FILE 'NET1',
C             DATA (5), ITEM 'END'. ONLY THE FIRST RECEIVER POINT IS
C             CONSIDERED.
C     'IND'...INPUT INDEX FILE.  THE NAME READ FROM FILE 'NET1',
C             DATA (7.1), ITEM 'IND(I)', AND CHECKED WITH FILE 'NET2',
C             DATA (7.1), ITEM 'IND(NREFL-I)'.  CONSIDERED ONLY IF THE
C             FILENAME IS NOT BLANK.
C     'TT1'...INPUT FILE CONTAINING TRAVEL-TIMES FROM THE SOURCE TO ALL
C             NETWORK NODES.  THE NAME READ FROM FILE 'NET1', DATA
C             (7.1), ITEM 'TT(I)'.
C     'TT2'...INPUT FILE CONTAINING TRAVEL-TIMES FROM THE RECEIVER TO
C             ALL NETWORK NODES.  THE NAME READ FROM FILE 'NET1', DATA
C             (7.1), ITEM 'TT(NREFL-I)'.
C     'OUT'...OUTPUT INDEX FILE.  THE NAME READ FROM FILE 'NET3',
C             DATA (7.1), ITEM 'IND(I)'.  THIS FILENAME MAY COINCIDE
C             WITH 'IND', IN SUCH A CASE THE OLD INDEX FILE 'IND' WILL
C             BE OVERWRITTEN BY THE NEW ONE.
C THE STRUCTURE OF THE ABOVE LISTED DATA FILES IS BRIEFLY DESCRIBED
C BELOW.  THE STRUCTURE OF THE FILES MENTIONED IN DATA FILES 'NET*' BUT
C NOT USED BY THIS PROGRAM IS NOT DESCRIBED HERE, REFER TO 'NET.FOR'.
C
C INPUT FILE 'SRC':
C (1) SEVERAL (0 TO 20) STRINGS TERMINATED BY / (A SLASH).
C (2) X1S,X2S,X3S,TTS,TTSERR/
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.  SEE DESCRIPTION OF 'END'.
C     TTSERR..INITIAL VALUE OF THE ARRIVAL TIME ERROR AT THE SOURCE
C             POINT.  IT IS LIKELY ZERO AT THE ACTUAL SOURCE.  SEE
C             DESCRIPTION OF 'END'.
C     DEFAULT: TTS=0, TTSERR=0.
C OTHER SOURCE POINTS ARE NOT READ AND SHOULD BE AVOIDED IN THESE DATA.
C
C INPUT FILE 'END':
C (1) SEVERAL (0 TO 20) STRINGS TERMINATED BY / (A SLASH).
C (2) X1R,X2R,X3R,TTMAX1,TTMAX2,/
C     X1R,X2R,X3R... COORDINATES OF THE RECEIVER.
C     TTMAX1+TTMAX2+TTS-TTSERR... REAL NUMBER: MAXIMUM SUM OF TRAVEL
C             TIMES, LIMITING THE FRESNEL VOLUME.  POSSIBLE PHYSICAL
C             MEANING:
C             TTMAX1=TT... FIRST-ARRIVAL TIME BETWEEN THE SOURCE AND THE
C               RECEIVER.
C             TTMAX1=TTERR... ESTIMATED MAXIMUM ERROR OF THE
C               FIRST-ARRIVAL TIME.
C             TTS,TTSERR... SEE DATA FILE 'SRC'.
C     DEFAULT: TTMAX2=0.
C DATA FOR OTHER RECEIVERS ARE IGNORED.
C
C INPUT INDEX FILE 'IND' AND OUTPUT INDEX FILE 'OUT':
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 'TT1' AND 'TT2' (TIME FIELDS):
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 (NUMBER OF SMALL BRICKS)
C             ON INPUT.
C
C-----------------------------------------------------------------------
C
      PROGRAM NETIND
C
C.......................................................................
C
      CHARACTER*80 FNET1,FNET2,FNET3,FNET4,FSRC,FEND,FIND,FTT1,FTT2,FOUT
      CHARACTER*80 FVEL1,FICB1,FIND2,FVEL2,FICB2
      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 X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,TT1(MTT),TT2(MTT),TTMAX
      REAL AUX1,AUX2,AUX3,AUX4,AUX5,AUX6
C
C     FNET1,FNET2,FNET3,FNET4... 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.OUT'
      FNET4='NET4.OUT'
      MIND=MPOS
      MGRID=MTT
      WRITE(*,'(2A)')  ' ENTER 4 NAMES OF INPUT FILES FOR ''NET'', AND '
     *                                     ,'MAX.NUMBER OF SMALL BRICKS'
      WRITE(*,'(2A,2(I8,A))')    ' (DEFAULT ''NET1.DAT'', ''NET2.DAT'','
     *              ,' ''NET3.OUT'', ''NET4.OUT'',',MIND,',',MGRID,'): '
      READ(*,*) FNET1,FNET2,FNET3,FNET4,MIND,MGRID
C     FNET1,FNET2,FNET3,FNET4 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 FOR THE NET PROGRAM:
      OPEN(LU1,FILE=FNET1,STATUS='OLD')
C     (1) NUMBERS OF GRIDPOINTS:
      N3=1
      L1=1
      L2=1
      L3=1
      READ(LU1,*) N1,N2,N3,L1,L2,L3
      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: NUMBER OF GRIDPOINTS IS NOT POSITIVE'
        STOP
      END IF
C     (2) BOUNDARIES OF MODEL VOLUME:
      READ(LU1,*) X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX
C     (3) SIZE OF THE FORWARD STAR:
      READ(LU1,*) I
C     (4) FACTORIZED ANISOTROPY:
      READ(LU1,*) (AUX1,I=1,21)
C     (5) NAMES OF THE FILES WITH SOURCE, RECEIVERS, RAYS, AND ERRORS:
      FEND=' '
      READ(LU1,*) FSRC,FAUX,FAUX,FEND
      IF(FEND.EQ.' ') THEN
        PAUSE'ERROR: NAME OF FILE WITH TIMES AT RECEIVERS NOT SPECIFIED'
        STOP
      END IF
C     (6) NUMBER OF REFLECTIONS:
      NREFL=0
      READ(LU1,*) NREFL
C     (7) 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
          PAUSE 'ERROR: NO INDEX FILE SPECIFIED'
          STOP
        END IF
      END IF
C
C     END OF READING THE 1-ST MAIN INPUT DATA FILE:
      CLOSE(LU1)
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     FSRC AND FEND ARE THE SOURCE AND ENDPOINT (REVCEIVER) 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     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
        PAUSE 'ERROR: TOO MANY GRIDPOINTS TO BE STORED IN ARRAY IND'
        STOP
      END IF
C
C     READING THE 2-ND INPUT FILE FOR THE NET PROGRAM:
      OPEN(LU1,FILE=FNET2,STATUS='OLD')
C     (1) NUMBERS OF GRIDPOINTS:
      IN3=1
      IL1=1
      IL2=1
      IL3=1
      READ(LU1,*) IN1,IN2,IN3,IL1,IL2,IL3
      IF(IN1.NE.N1.OR.IN2.NE.N2.OR.IN3.NE.N3.OR.
     *   IL1.NE.L1.OR.IL2.NE.L2.OR.IL3.NE.L3) THEN
        PAUSE 'ERROR: DIFFERENT NUMBERS OF GRIDPOINTS'
        STOP
      END IF
C     (2) BOUNDARIES OF MODEL VOLUME:
      READ(LU1,*) AUX1,AUX2,AUX3,AUX4,AUX5,AUX6
      IF(X1MIN.NE.AUX1.OR.X1MAX.NE.AUX2.OR.X2MIN.NE.AUX3.OR.
     *   X2MAX.NE.AUX4.OR.X3MIN.NE.AUX5.OR.X3MAX.NE.AUX6) THEN
        PAUSE 'ERROR: DIFFERENT SIZE OF MODEL VOLUME IN THE SECOND DATA'
        STOP
      END IF
C     (3) SIZE OF THE FORWARD STAR:
      READ(LU1,*) I
C     (4) FACTORIZED ANISOTROPY:
      READ(LU1,*) (AUX1,I=1,21)
C     (5) NAMES OF THE FILES WITH SOURCE, RECEIVERS, RAYS, AND ERRORS:
      READ(LU1,*) FAUX,FAUX,FAUX,FAUX
C     (6) NUMBER OF REFLECTIONS:
      I=0
      READ(LU1,*) I
      IF(I.NE.NREFL) THEN
        PAUSE'ERROR: DIFFERENT NUMBER OF REFLECTIONS IN THE SECOND DATA'
        STOP
      END IF
C     (7) 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
        PAUSE 'ERROR: DIFFERENT INPUT INDEX FILES'
        STOP
      END IF
      IF(FVEL1.NE.FVEL2) THEN
        PAUSE 'ERROR: DIFFERENT VELOCITY FILES'
        STOP
      END IF
      IF(FICB1.NE.FICB2) THEN
        PAUSE 'ERROR: DIFFERENT BLOCK FILES'
        STOP
      END IF
C
C     END OF READING THE 2-ND MAIN INPUT DATA FILE:
      CLOSE(LU1)
C     FTT2 IS THE INPUT TRAVEL-TIME FILE, WITH TIMES FROM THE RECEIVER.
C
C     READING THE 3-RD INPUT FILE FOR THE NET PROGRAM:
      OPEN(LU1,FILE=FNET3,STATUS='OLD')
C     (1) NUMBERS OF GRIDPOINTS:
      L1MAX=0
      L2MAX=0
      L3MAX=0
      READ(LU1,*) IN1,IN2,IN3,L1MAX,L2MAX,L3MAX
C     (2) BOUNDARIES OF MODEL VOLUME:
      READ(LU1,*) AUX1,AUX2,AUX3,AUX4,AUX5,AUX6
      IF(X1MIN.NE.AUX1.OR.X1MAX.NE.AUX2.OR.X2MIN.NE.AUX3.OR.
     *   X2MAX.NE.AUX4.OR.X3MIN.NE.AUX5.OR.X3MAX.NE.AUX6) THEN
        PAUSE 'ERROR: DIFFERENT SIZE OF MODEL VOLUME IN THE THIRD DATA'
        STOP
      END IF
C     (3) SIZE OF THE FORWARD STAR:
      READ(LU1,*) I
C     (4) FACTORIZED ANISOTROPY:
      READ(LU1,*) (AUX1,I=1,21)
C     (5) NAMES OF THE FILES WITH SOURCE, RECEIVERS, RAYS, AND ERRORS:
      READ(LU1,*) FAUX,FAUX,FAUX,FAUX
C     (6) NUMBER OF REFLECTIONS:
      I=0
      READ(LU1,*) I
      IF(I.NE.NREFL) THEN
        PAUSE 'ERROR: DIFFERENT NUMBER OF REFLECTIONS IN THE THIRD DATA'
        STOP
      END IF
C     (7) 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
C
C     END OF READING THE 3-RD MAIN INPUT DATA FILE:
      CLOSE(LU1)
C     OUTPUT BIG BRICK MAY HAVE AT MOST L1MAX*L2MAX*L3MAX SMALL BRICKS.
C     FOUT IS THE OUTPUT INDEX FILE.
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
        PAUSE 'ERROR: 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.                                 '
              write(*,*) irec,irec1,irec2,irec3
      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
              write(*,*) irec,irec1,irec2,irec3
              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     UPDATING 'NET3' AND 'NET4':
      OPEN(LU2,STATUS='SCRATCH')
      OPEN(LU1,FILE=FNET3,STATUS='OLD')
      READ(LU1,*) IN1,IN2,IN3,IL1,IL2,IL3
      CALL COPY(LU1,LU2)
      REWIND(LU2)
      REWIND(LU1)
      WRITE(LU1,'(6I5,A)') N1,N2,N3,L1,L2,L3,' /'
      CALL COPY(LU2,LU1)
      CLOSE(LU1)
      IF(FNET4.NE.' ') THEN
        REWIND(LU2)
        OPEN(LU1,FILE=FNET4)
        READ(LU1,*) IN1,IN2,IN3,IL1,IL2,IL3
        CALL COPY(LU1,LU2)
        REWIND(LU2)
        REWIND(LU1)
        WRITE(LU1,'(6I5,A)') N1,N2,N3,L1,L2,L3,' /'
        CALL COPY(LU2,LU1)
        CLOSE(LU1)
      END IF
      CLOSE(LU2)
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
          PAUSE 'ERROR: 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
      SUBROUTINE COPY(LU1,LU2)
C
C SUBROUTINE COPYING A SEQUENTIAL FORMATTED FILE TO ANOTHER SEQUENTIAL
C FORMATTED FILE.
C
C INPUT:
C     LU1...  LOGICAL UNIT NUMBER OF THE SOURCE FILE.  THE FILE IS
C             TREATED AS ALREADY OPEN.  THE RECORDS FROM THE CURRENT
C             POSITION TO THE END OF FILE ARE COPIED.
C     LU2...  LOGICAL UNIT NUMBER OF THE TARGET FILE.  THE FILE IS
C             TREATED AS ALREADY OPEN.  THE COPY IS APPENDED AT THE
C             CURRENT POSITION.
C
C NO OUTPUT.
C
C DATE: 1993, OCTOBER 18
C CODED BY: LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     AUXILIARY STORAGE LOCATIONS:
      INTEGER NLINE
      PARAMETER (NLINE=255)
      CHARACTER*(NLINE) LINE
      INTEGER I
C
   10 CONTINUE
        READ(LU1,'(A)',END=99) LINE
        DO 11 I=NLINE,1,-1
          IF(LINE(I:I).NE.' ') THEN
            WRITE(LU2,'(A)') LINE(1:I)
            GO TO 20
          END IF
   11   CONTINUE
        WRITE(LU2,'(A)')
   20   CONTINUE
      GO TO 10
C
   99 CONTINUE
      RETURN
      END
C
C=======================================================================
C