C Program 'modchk' to perform the consistency check of the data
C describing the structure of the model.
C
C Version: 5.00
C Date: 1996, September 30
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 demmed 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 neigbouring
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 usually 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 ovelaping simple
C and complex blocks in the model.
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 succesfully
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 the data files:
C
C The data are read in by the list directed input (free format).
C In the description of data files, each numbered paragraph indicates
C the beginning of a new input operation (new READ statement).
C If the symbolic name of the input variable is enclosed in apostrophes,
C the corresponding value in input data is of the type CHARACTER, i.e.
C it should be a character string enclosed in apostrophes. If the first
C letter of the symbolic name is I-N, the corresponding value is of the
C type INTEGER. Otherwise, the input parameter is of the type REAL and
C may or may not contain a decimal point.
C
C Input data read from the * external unit:
C The interactive * external unit may also be redirected to the file
C containing the relevant data.
C (1) 'FMODEL','FCHK','FOUT',N1,N2,N3,LOVER,/
C 'FMODEL'... String with the name of the input data file specifying
C the model to be checked for consistency.
C 'FCHK'..String with the name of the input data file containing the
C additional free-space blocks.
C 'FCHK'=' ': No additional free-space blocks defined.
C 'FOUT'..String with the name of the output data file witth the
C report on the consistency check.
C N1,N2,N3... Specification of the regular rectangular grid of
C points for the numerical test for overlapping simple or
C complex blocks in the model. The model volume is divided
C into N1*N2*N3 rectangular cells. The (N1+1)*(N2+1)*(N3+1)
C gridpoints are the corner 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 LOVER...Logical value indicating whether the overlapping simple
C blocks are allowed in the model.
C LOVER=.FALSE.: Overlapping simple blocks are reported.
C LOVER=.TRUE.: Overlapping complex blocks are reported.
C Default: 'FMODEL'='model.dat', 'FCHK'=' ', 'FOUT'='modchk.out',
C N1=0, N2=0, N3=0, LOVER=.FALSE.
C
C Input data file 'FMODEL' specifying the model:
C Refer to the description in 'model.for'.
C
C Input file 'FCHK' specifying the 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
C-----------------------------------------------------------------------
C
C Input data:
CHARACTER*80 FMODEL,FCHK,FOUT
LOGICAL LOVER
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 MUB
PARAMETER (MUB=10240)
INTEGER KUB(MUB),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 input data:
C
C Main input data:
WRITE(*,'(A)') ' Enter FMODEL,FCHK,FOUT,N1,N2,N3,LOVER: '
FMODEL='model.dat'
FCHK=' '
FOUT='modchk.out'
N1=0
N2=0
N3=0
LOVER=.FALSE.
READ(*,*) FMODEL,FCHK,FOUT,N1,N2,N3,LOVER
WRITE(*,'(A)') '+Reading input data. '
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
PAUSE 'Error 310 in MODCHK: Insufficient memory in /MODELC/'
STOP
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
PAUSE
* 'Error 311 in MODCHK: Block bounded by wrong interface'
STOP
ELSE IF(KSB(I).EQ.0) THEN
KSB(J)=I-1
L=I
GO TO 13
END IF
12 CONTINUE
PAUSE 'Error 310 in MODCHK: Insufficient memory in /MODELC/'
STOP
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
PAUSE 'Error in MODCHK: Array dimension MUB too small'
STOP
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)') '+ % 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)
IF(NB.GT.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
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 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) 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 'model.for'
INCLUDE 'metric.for'
INCLUDE 'srfc.for'
INCLUDE 'parm.for'
INCLUDE 'val.for'
INCLUDE 'fit.for'
C
C=======================================================================
C