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