C Program MODCHK to perform the consistency check of the data describing C the structure of the model. C C Version: 5.50 C Date: 2001, April 5 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 Logical test for possible cavities in the model: C C The MODEL package does not require the free-space block to be C explicitly specified. All simple blocks, not defined in the data C file describing the model, are deemed to be free-space blocks. C In such a case, it is hard to check for the possible cavities C in the model. Although the 'modchk' program reports all simple C blocks that could, in principle, form such cavities, most of the C reported simple blocks are parts of the desired free space and C often mutually overlap. The number of reported undefined simple C blocks is thus often too large to be comfortably checked by a C user. C C It is thus recommended to specify free-space simple blocks in C addition to the material simple blocks. It may be done directly C in the data file specifying the model. All given simple blocks C which do not form material complex blocks are deemed to be C free-space simple blocks. Since the free-space simple blocks need C not be explicitly specified in the model data, they may also be C specified in a separate file submitted to the 'modchk' program. C The table of simple blocks in the separate file is just a C continuation of the corresponding table in the model data file. C Then the list of undefined free-space simple blocks reported by C the 'modchk' program should be much more useful. C C Although the undefined free-space blocks need not necessarily C physically exist and form cavities in the particular model, they C should be removed. They may often be unified with a neighbouring C simple block from which they are separated by an interface. C Especially if a user knows that an undefined simple block does C not physically exist, there is often no reason to separate it C from neighbouring defined simple blocks. If an undefined simple C block is allowed to form free space, the list of free-space blocks C should be updated. C C It is recommended to fix all undefined free-space simple blocks C before performing detailed numerical test for overlapping simple C and complex blocks in the model. C C Example: Assume three mutually intersecting interfaces 1, 2, 3, C limiting four simple blocks in the following way: C +3- C / C / C / | C / | C (+1,+3) / | C / | C / | C + / (+1,-2,-3)| C 1-----------------|(+2,-3) C - (-1,-2) | C | C -2+ C Then simple block (-1,+2,+3) is an undefined free-space simple C block and is reported by program 'modchk'. Although simple block C (-1,+2,+3) probably does not exist, it is recommended to make it C defined. The best way would be to unify it with an existing C simple block. Unfortunately, block (-1,+2,+3) cannot by simply C unified with any of the above simple blocks (+1,+3), (-1,-2), C (+2,-3), (+1,-2,-3). There are still at least two possibilities C how to fix this problem: C (a) C A safer but slower way is to define also simple block (-1,+2,+3) C and to attach it to a material of free-space complex block. C (b) C If we are sure that simple block (-1,+2,+3) does not exist and C thus cannot form an undesired cavity in the model, we may leave it C undefined free-space simple block, and run the 'modchk' program C with parameter LFREE=0 and dense test grid. C C Numerical test for overlapping simple or complex blocks in the model: C C The model is covered by a regular rectangular grid of points. C The position of each gridpoint with respect to all simple blocks C is then determined. Normally, each gridpoint situated inside C two or more simple blocks is reported. Here the word 'inside', C naturally, does not include the boundary of a simple block. C C Although overlapping simple blocks in the model are allowed if C they form the same complex blocks, the are not recommended. C If the overlapping simple blocks are present in the model, the C test for overlapping simple blocks may be disabled. In such a C case, only overlapping complex blocks are reported. C C It is reasonable to start the numerical test for overlapping C blocks with a small numbers of gridpoints (e.g., 10*10*10), and C increase the number of gridpoints if the model passes successfully C the first test, and so on. The number of gridpoints for the final C test is limited especially by the computational time. C C....................................................................... C C Description of data files: C C Input data read from the standard input device (*): C The data are read by the list directed input (free format) and C consist of a single string 'SEP': C 'SEP'...String in apostrophes containing the name of the input C SEP parameter or history file with the input data. C No default, 'SEP' must be specified and cannot be blank. C C C Input data file 'SEP': C File 'SEP' has the form of the SEP C parameter file. The parameters, which do not differ from their C defaults, need not be specified in file 'SEP'. C Names of the file with model and of the output file: C MODEL='string' ... Name of the input data file specifying C the model to be checked for consistency. C Default: MODEL='model.dat' C Description of MODEL C Example of MODEL C MODCHK='string' ... Name of the input data file containing C the additional free-space blocks. C Description of input file MODCHK C Default: MODCHK=' ' (no additional free-space blocks) C MODLOG='string' ... Name of the output data file with the C report on the consistency check. C Default: MODLOG='modchk.out' C Data for the model consistency check: C N1=integer, N2=integer, N3=integer ... Specification of the C regular rectangular grid of points for the numerical test C for overlapping simple or complex blocks in the model. C The model volume is divided into N1*N2*N3 rectangular C cells. The (N1+1)*(N2+1)*(N3+1) gridpoints are the corner C points of the cells. C If one of N1,N2,N3 is zero, the model is assumed 2-D. C The corresponding 2-D test grid then passes through the C centre of the 3-D model volume. C If all N1,N2,N3 are zero, the numerical test for C overlapping simple or complex blocks is not performed. C Defaults: N1=0, N2=0, N3=0 C LOVER=integer ... Value indicating whether the overlapping simple C blocks are allowed in the model. C LOVER=0: Overlapping simple blocks are reported. C LOVER=1: Overlapping complex blocks are reported. C Default: LOVER=0 C LFREE=integer ... Value indicating whether undefined free-space C simple blocks are to be reported together with overlapping C blocks. C LFREE=0: Undefined free-space simple blocks are not C reported. C LFREE=1: Undefined free-space simple blocks are reported. C Default: LFREE=0 C C C Input file 'MODCHK' specifying additional free-space simple blocks: C (1) NSBADD C Number of additional free-space simple blocks defined within this C data file. See also description of data (5) in 'model.for'. C (2) NSBADD input operations (READ statements): C For each simple block with index ISB, the indices of the surfaces C forming the set F(+) and the indices of the surfaces forming the C set F(-). The indices of surfaces from F(+) must be positive, the C indices of surfaces from F(-) must be indicated by negative signs. C The indices may be specified in an arbitrary order and must be C terminated by a slash. These data lines form the continuation of C of data (6) described in 'model.for'. C C======================================================================= C C Common blocks /MODELT/ and /MODELC/: INCLUDE 'model.inc' C model.inc C C----------------------------------------------------------------------- C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc C C Allocation of the working array: INTEGER MUB PARAMETER (MUB=MRAM) INTEGER KUB(MUB) EQUIVALENCE (KUB,RAM) C C----------------------------------------------------------------------- C CHARACTER*80 FILSEP INTEGER LU0 PARAMETER (LU0=1) C Input data: CHARACTER*80 FMODEL,FCHK,FOUT INTEGER ILOVER,ILFREE LOGICAL LOVER,LFREE INTEGER LU1,N1,N2,N3 PARAMETER (LU1=1) C C Temporary storage locations: INTEGER NSBOLD,I,J,K,L C C Logical test for undefined free-space simple blocks: INTEGER NUB,KUBNUB,NUNDEF,ISRF,ISB INTEGER MININT,MINSUB,MINKUB,NUMINT,NUMSUB C NUB... Starting index of the last candidate for undefined block C in array KUB. Each candidate occupies NSB+1 locations. C KUB(NUB)... Number of interfaces limiting the last candidate. C KUB(NUB+1:NUB+KUB(NUB))... Indices of interfaces limiting the last C candidate, including signs. C KUB(NUB+KUB(NUB)+1,NUB+NSB)... Indices of simple blocks with which C complements the candidate should be intersected. C Similarly for preceding candidates, if any. C NUNDEF... Number of undefined free-space simple blocks. C C Numerical test for overlapping blocks: INTEGER I1,I2,I3,I123,N123,NB,IB(100) REAL COOR(3) C C----------------------------------------------------------------------- C C Reading name of SEP file with input data: WRITE(*,'(A)') '+MODCHK: Enter input filename: ' FILSEP=' ' READ(*,*) FILSEP WRITE(*,'(A)') '+MODCHK: Working ... ' C C Reading all data from the SEP file into the memory: IF (FILSEP.NE.' ') THEN CALL RSEP1(LU0,FILSEP) ELSE C MODCHK-05 CALL ERROR('MODCHK-05: SEP file not given') C Input file in the form of the SEP (Stanford Exploration Project) C parameter or history file must be specified. C There is no default filename. ENDIF C C Reading input parameters from the SEP file: CALL RSEP3T('MODEL',FMODEL,'model.dat') CALL RSEP3T('MODCHK',FCHK,' ') CALL RSEP3T('MODLOG',FOUT,'modchk.out') CALL RSEP3I('N1',N1,0) CALL RSEP3I('N2',N2,0) CALL RSEP3I('N3',N3,0) CALL RSEP3I('LOVER',ILOVER,0) CALL RSEP3I('LFREE',ILFREE,0) IF (ILOVER.EQ.0) THEN LOVER=.FALSE. ELSE LOVER=.TRUE. ENDIF IF (ILFREE.EQ.0) THEN LFREE=.FALSE. ELSE LFREE=.TRUE. ENDIF C C Reading data for model: OPEN(LU1,FILE=FMODEL,STATUS='OLD') CALL MODEL1(LU1) CLOSE(LU1) C C Reading additional free-space simple blocks: NSBOLD=NSB IF(FCHK.NE.' ') THEN OPEN(LU1,FILE=FCHK,STATUS='OLD') C Number of additional simple blocks READ(LU1,*) K IF(KSB(NSB)+K.GT.MSB) THEN C MODCHK-01 CALL ERROR('MODCHK-01: Insufficient memory in /MODELC/') END IF DO 10 J=KSB(NSB),NSB+1,-1 KSB(J+K)=KSB(J) 10 CONTINUE DO 11 J=NSB,1,-1 KSB(J)=KSB(J)+K 11 CONTINUE NSB=NSB+K C Reading indices of surfaces bounding additional simple blocks: L=KSB(NSBOLD)+1 DO 14 J=NSBOLD+1,NSB READ(LU1,*) (KSB(I),I=L,MSB) DO 12 I=L,MSB IF(IABS(KSB(I)).GT.NSRFCS) THEN C MODCHK-02 CALL ERROR('MODCHK-02: Block bounded by wrong interface') ELSE IF(KSB(I).EQ.0) THEN KSB(J)=I-1 L=I GO TO 13 END IF 12 CONTINUE C MODCHK-03 CALL ERROR('MODCHK-03: Insufficient memory in /MODELC/') 13 CONTINUE 14 CONTINUE END IF CLOSE(LU1) C C Opening output report file: OPEN(LU1,FILE=FOUT) WRITE(LU1,'(A)') 'Report on the consistency check of the model' WRITE(LU1,'(A)') WRITE(LU1,'(I3,A)') NSBOLD, * ' simple blocks defined in model data file' WRITE(LU1,'(A)') FMODEL IF(FCHK.EQ.' ') THEN WRITE(LU1,'(2A)') 'None additional free-space simple block', * ' defined for the consistency check.' ELSE WRITE(LU1,'(I3,A)') NSB-NSBOLD, * ' additional free-space simple blocks defined in file' WRITE(LU1,'(A)') FCHK END IF WRITE(LU1,'(A)') C C....................................................................... C C Defining free-space complex block: C C Adding the free-space complex block: DO 20 J=KCB(NCB),NCB+1,-1 KCB(J+1)=KCB(J) 20 CONTINUE DO 21 J=NCB,1,-1 KCB(J)=KCB(J)+1 21 CONTINUE L=KCB(NCB) NCB=NCB+1 KCB(NCB)=L C C Creating the list of free-space simple blocks: C List of all simple blocks: DO 22 J=1,NSB KCB(L+J)=J 22 CONTINUE C Removing material simple blocks from the list: DO 23 J=NCB+1,L IF(KCB(L+KCB(J)).EQ.0) THEN WRITE(LU1,'(A,I3,A)') 'Simple block',KCB(J), * ' is situated in two different complex blocks!' ELSE KCB(L+KCB(J))=0 END IF 23 CONTINUE C List of remaining (i.e. free-space) simple blocks: DO 24 J=1,NSB IF(KCB(L+J).NE.0) THEN KCB(NCB)=KCB(NCB)+1 KCB(KCB(NCB))=KCB(L+J) END IF 24 CONTINUE C IF(KCB(NCB).EQ.L) THEN WRITE(LU1,'(A)') * 'There is no explicitly defined free space in the model.' ELSE WRITE(LU1,'(A)') 'Free space is composed of simple blocks:' WRITE(LU1,'(20I4)') (KCB(J),J=L+1,KCB(NCB)) END IF WRITE(LU1,'(A)') C C....................................................................... C WRITE(*,'(A)') '+Logical test for undefined free-space blocks. ' WRITE(LU1,'(2A)') 'Undefined free-space simple blocks ', * '(indices of surfaces limiting each one):' NUNDEF=0 C NUB=1 KUB(NUB)=0 DO 30 I=1,NSB KUB(NUB+I)=I 30 CONTINUE C 40 CONTINUE IF(NUB.LT.0) THEN C No candidates for undefined free-space simple block left: GO TO 70 END IF C KUBNUB=KUB(NUB) IF(KUBNUB.GE.NSB) THEN C Printing undefined free-space simple block: I=NUB DO 41 J=NUB+1,NUB+NSB IF(KUB(J).NE.0) THEN I=I+1 KUB(I)=KUB(J) END IF 41 CONTINUE NUNDEF=NUNDEF+1 WRITE(LU1,'(20I4/99(4X,20I4))') (KUB(J),J=NUB+1,I) NUB=NUB-NSB-1 GO TO 40 END IF C C Selecting simple block with minimum number MININT of C intersections: MININT=999999 C Loop over simple blocks DO 55 K=NUB+KUBNUB+1,NUB+NSB ISB=KUB(K) NUMINT=0 C Loop over surfaces limiting the simple block C (each surface separates the simple block from its complement) DO 53 J=KSB(ISB-1)+1,KSB(ISB) ISRF=KSB(J) C Loop over surfaces limiting the candidate DO 51 I=NUB+1,NUB+KUBNUB IF(KUB(I).EQ.ISRF) THEN C Empty intersection of the candidate and the complement GO TO 52 ELSE IF(KUB(I).EQ.-ISRF) THEN C The candidate is a subset of the complement NUMINT=1 NUMSUB=1 GO TO 54 END IF 51 CONTINUE NUMINT=NUMINT+1 52 CONTINUE 53 CONTINUE NUMSUB=0 54 CONTINUE IF(NUMINT.LT.MININT) THEN MININT=NUMINT MINSUB=NUMSUB MINKUB=K END IF 55 CONTINUE C IF(MININT.EQ.0) THEN C Candidate discarded: NUB=NUB-NSB-1 ELSE IF(MINSUB.GT.0) THEN C The candidate is a subset of one of the complements: KUB(MINKUB)=KUB(NUB+KUBNUB+1) KUB(NUB)=KUBNUB+1 KUB(NUB+KUBNUB+1)=0 ELSE C Candidate is intersected with each complement: ISB=KUB(MINKUB) KUB(MINKUB)=KUB(NUB+KUBNUB+1) KUB(NUB)=KUBNUB+1 C Candidate is MININT times reproduced IF(NUB+(ISB+1)*MININT-1.GT.MUB) THEN C MODCHK-04 CALL ERROR('MODCHK-04: Array dimension MUB too small') END IF DO 60 I=NUB+NSB+1,NUB+(NSB+1)*MININT-1 KUB(I)=KUB(I-NSB-1) 60 CONTINUE C Loop over surfaces limiting the simple block C (each surface separates the simple block from its complement) K=NUB DO 63 J=KSB(ISB-1)+1,KSB(ISB) ISRF=KSB(J) C Loop over surfaces limiting the candidate DO 61 I=K+1,K+KUBNUB IF(KUB(I).EQ.ISRF) THEN C Empty intersection of the candidate and the complement GO TO 62 END IF 61 CONTINUE KUB(NUB+KUBNUB+1)=-ISRF NUB=NUB+NSB+1 62 CONTINUE 63 CONTINUE NUB=NUB-NSB-1 END IF GO TO 40 C 70 CONTINUE C IF(NUNDEF.LE.0) THEN WRITE(LU1,'(A)') * 'O.K. There is no undefined free space in the model.' ELSE WRITE(LU1,'(2A)') 'Please, check carefully the above list ', * 'and edit the data for simple blocks!' END IF WRITE(LU1,'(A)') C C....................................................................... C C Test for overlapping simple blocks or complex blocks: C IF(N1.LE.0.AND.N2.LE.0.AND.N3.LE.0) THEN WRITE(LU1,'(A)') * 'Numerical test for overlapping blocks not performed.' ELSE WRITE(*,'(A)') '+ 0% Numerical test for overlapping blocks. ' C IF(LOVER) THEN WRITE(LU1,'(A)') * 'Overlapping simple blocks are allowed in the model,' WRITE(LU1,'(A)') * 'test for overlapping simple blocks is not performed!' WRITE(LU1,'(A)') WRITE(LU1,'(2A)') * 'List of points situated in more than one complex block', * ' (0=free space):' ELSE WRITE(LU1,'(A)') * 'List of points situated in more than one simple block:' END IF C I123=0 N123=(N1+1)*(N2+1)*(N3+1) DO 83 I3=0,N3 IF(N3.LE.0) THEN COOR(3)=(BOUNDM(5)+BOUNDM(6))/2. ELSE COOR(3)=BOUNDM(5)+(BOUNDM(6)-BOUNDM(5))*FLOAT(I3)/FLOAT(N3) END IF DO 82 I2=0,N2 IF(N2.LE.0) THEN COOR(2)=(BOUNDM(3)+BOUNDM(4))/2. ELSE COOR(2)=BOUNDM(3) * +(BOUNDM(4)-BOUNDM(3))*FLOAT(I2)/FLOAT(N2) END IF DO 81 I1=0,N1 IF(N1.LE.0) THEN COOR(1)=(BOUNDM(1)+BOUNDM(2))/2. ELSE COOR(1)=BOUNDM(1) * +(BOUNDM(2)-BOUNDM(1))*FLOAT(I1)/FLOAT(N1) END IF CALL CHK(COOR,NB,IB,LOVER) C Subroutine CHK IF(NB.GT.1.OR.(LFREE.AND.NB.LT.1)) THEN WRITE(LU1,'(3F9.3,3X,15I3/(30X,15I3))') * (COOR(I),I=1,3),(IB(I),I=1,NB) END IF C Displaying progress on the console: I123=I123+100 IF(MOD(I123,N123).LT.100) THEN WRITE(*,'(A,I3)') '+',I123/N123 END IF 81 CONTINUE 82 CONTINUE 83 CONTINUE C WRITE(LU1,'(A)') '---------' END IF CLOSE(LU1) WRITE(*,'(A)') '+Done. ' STOP END C C======================================================================= C C C SUBROUTINE CHK(COOR,NB,IB,LOVER) REAL COOR(3) INTEGER NB,IB(*) LOGICAL LOVER C C This is an auxiliary subroutine to program MODCHK. It determines the C position of a given point with respect to all simple (LOVER=.FALSE.) C or complex (LOVER=.TRUE.) blocks. C C Input: C COOR... Array containing coordinates X1, X2, X3 of a given point. C LOVER...Logical value indicating whether the overlapping simple C blocks are allowed in the model. C LOVER=.FALSE.: Simple blocks containing the given point C are reported. C LOVER=.TRUE.: Complex blocks containing the given point C are reported. C None of the input parameters are altered. C C Output: C NB... Number of blocks inside which the given point is situated. C IB(1:NB)... Indices of blocks inside which the given point is C situated. C C Common block /MODELC/: INCLUDE 'model.inc' C model.inc C None of the storage locations of the common block are altered. C C Subroutines and external functions required: EXTERNAL NSRFC,SRFC2 INTEGER NSRFC C SRFC2 and subsequent routines... File 'srfc.for' and subsequent C files (like 'val.for' and 'fit.for'). C C Date: 1996, September 30 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER ISB,ICB,I REAL F(10),F1(100) C C....................................................................... C C Values of functions describing the surfaces in the model: DO 11 I=1,NSRFC() CALL SRFC2(I,COOR,F) F1(I)=F(1) 11 CONTINUE C C Loop for simple blocks: NB=0 DO 29 ISB=1,NSB C Loop for surfaces bounding simple block ISB: DO 21 I=KSB(ISB-1)+1,KSB(ISB) IF(F1(IABS(KSB(I)))*FLOAT(KSB(I)).LE.0.) THEN C The point is not inside the simple block. GO TO 25 END IF 21 CONTINUE C The point is inside the simple block. NB=NB+1 IB(NB)=ISB 25 CONTINUE 29 CONTINUE C IF (LOVER.AND.NB.GT.0) THEN C Loop for simple blocks inside which the point is situated: DO 39 ISB=1,NB C Loop for complex blocks: DO 32 ICB=1,NCB C Loop for simple blocks composing complex block ICB: DO 31 I=KCB(ICB-1)+1,KCB(ICB) IF(KCB(I).EQ.IB(ISB)) THEN C The point is inside the complex block. IB(ISB)=ICB GO TO 35 END IF 31 CONTINUE 32 CONTINUE 35 CONTINUE 39 CONTINUE C Removing identical complex blocks: ICB=1 DO 49 ISB=2,NB DO 42 I=1,ICB IF(IB(ISB).EQ.IB(I)) THEN C Repeated index of complex block. GO TO 45 END IF 42 CONTINUE ICB=ICB+1 IB(ICB)=IB(ISB) 45 CONTINUE 49 CONTINUE NB=ICB C Renaming free-space complex block from NCB to 0: DO 51 ICB=1,NB IF(IB(ICB).EQ.NCB) THEN IB(ICB)=0 END IF 51 CONTINUE END IF C RETURN END C C======================================================================= C INCLUDE 'error.for' C error.for INCLUDE 'sep.for' C sep.for INCLUDE 'length.for' C length.for INCLUDE 'model.for' C model.for INCLUDE 'metric.for' C metric.for INCLUDE 'srfc.for' C srfc.for INCLUDE 'parm.for' C parm.for INCLUDE 'val.for' C val.for INCLUDE 'fit.for' C fit.for C C======================================================================= C