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