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 http://sw3d.cz/staff/klimes.htm 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 C 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 C '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), C they may be used to start ITER iterations of ray C tracing within 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 C 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 C 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 C 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 C 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 C 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 C 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 C ICB(I) C 'TT(I)'... Name of the file containing travel-times at all C network nodes after I reflections. C Description of data C 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 C 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 C 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 C INTF(I) C NET1, NET2 and NET3: This line is skipped. C Examples of data sets C NET1, C NET2 and C 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