C 
C Program INTF to check the positions of given points with respect to
C interfaces in the model.
C
C Version: 5.10
C Date: 1997, September 26
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                                                    
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) 'MODEL','INTF','OUT',KSRFC,KOLUMN,/
C     'MODEL'... Input data file describing the model, it is described
C             in the subroutine file 'model.for'.
C             Description of file MODEL
C     'INTF'..String with the name of the input data file containing the
C             points situated at (or close to) interface(s).
C     'OUT'...String with the name of the output data file containing the
C             value F(X1,X2,X3) of the function describing the
C             corresponding surface, evaluated at the given points.
C     KSRFC...KSRFC=0: Points submitted in file 'INTF' correspond to
C               various surfaces in the model.  Each point thus must be
C               supplemented with the index of the surface.
C             KSRFC.NE.0: Points submitted in file 'INTF' correspond to
C               surface number IABS(KSRFC).
C             KSRFC.GT.0: Input file 'INTF' has format 'Points', see
C               'formsdoc.htm'.
C             KSRFC.LT.0: Input file 'INTF' has format 'Lines', see
C               'formsdoc.htm'.
C     KOLUMN..Specifies the position in output file 'OUT' where to write
C             value F(X1,X2,X3) of the function describing the
C             corresponding surface, evaluated at the given points.
C     Default: 'MODEL'='model.dat', 'INTF'='intf.dat', 'OUT'='intf.out',
C             KSRFC=0, KOLUMN=4.
C
C Input file 'INTF':
C For KSRFC.GT.0:
C     Input file 'INTF' has format POINTS, see 'formsdoc.htm'.
C     Description of form POINTS
C For KSRFC.LT.0:
C     Input file 'INTF' has format LINES, see 'formsdoc.htm'.
C     Description of form LINES
C For KSRFC.EQ.0:
C     Several lines terminated by a slash or EOF.  Each line corresponds
C     to a given point and contains 4 numbers:
C     ISRFC,X1,X2,X3
C     ISRFC...Index of the surface.
C     X1,X2,X3... Coordinates of the point.
C In all cases, the possible quantities following the coordinates are
C not considered.
C
C Output file 'OUT' has similar form as input file 'INTF', depending on
C KSRFC:
C For KOLUMN=1,2,3:
C     KOLUMN-th coordinate of the input point is replaced by the value
C     F(X1,X2,X3) of the function describing the corresponding surface.
C     This option may be used to project the points or lines given by 2
C     of 3 coordinates on the surfaces described in the simple form of
C     F(X1,X2,X3)=W(X1,X2)-X3 or F(X1,X2,X3)=W(X1,X3)-X2 or
C     F(X1,X2,X3)=W(X2,X3)-X1.
C For KSRFC=0 and KOLUMN=0:
C     All 3 coordinates are replaced by single value F(X1,X2,X3).
C     If the point is situated exactly at the surface, F(X1,X2,X3)=0.
C Otherwise:
C     The 3 coordinates are supplemented with the 4th number, the value
C     of F(X1,X2,X3) to check the position of the point with respect to
C     the surface.
C
C=======================================================================
C
C External procedures directly referred:
      EXTERNAL LENGTH,MODEL1,SRFC2
      INTEGER  LENGTH
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-----------------------------------------------------------------------
C
C     Filenames:
      CHARACTER*80 FILE0,FILE1,FILE2
C
C     Logical unit numbers:
      INTEGER LU1,LU2
      PARAMETER (LU1=11)
      PARAMETER (LU2=12)
C
C     Data:
      CHARACTER*80 TEXT
      LOGICAL NEWLIN
      INTEGER KSRFC,ISRFC,NPTS,I,KOLUMN
      REAL COOR(3),F(10),FMAX,FRMS,FABS,FAVE,GIANT
      PARAMETER (GIANT=9.9E9)
C
C.......................................................................
C
C     Opening data files and reading the input data:
C
C     Main input data file read from the interactive device (*):
      WRITE(*,'(A)')
     *      ' Enter model, point, output filenames, and surface index: '
      FILE0='model.dat'
      FILE1='intf.dat'
      FILE2='intf.out'
      KSRFC=0
      KOLUMN=4
      READ(*,*) FILE0,FILE1,FILE2,KSRFC,KOLUMN
C
      WRITE(*,'(A)')
     *  '+Reading input data for the model.                            '
      OPEN(LU1,FILE=FILE0,STATUS='OLD')
      CALL MODEL1(LU1)
      CLOSE(LU1)
C
      WRITE(*,'(A)')
     *  '+                                                             '
      OPEN(LU1,FILE=FILE1,STATUS='OLD')
      OPEN(LU2,FILE=FILE2)
C
      IF(KSRFC.NE.0) THEN
        ISRFC=IABS(KSRFC)
        READ(LU1,*,END=90) (TEXT,I=1,20)
        WRITE(LU2,'(A)') '/'
      END IF
C
C.......................................................................
C
C     Check of the positions
C
      NEWLIN=.TRUE.
      NPTS=0
      FMAX=0.
      FRMS=0.
      FABS=0.
      FAVE=0.
   10 CONTINUE
        IF(KSRFC.GT.0) THEN
C         Input format 'Points':
          TEXT='$'
          COOR(1)=0.
          COOR(2)=0.
          COOR(3)=0.
          READ(LU1,*,END=90) TEXT,COOR(1),COOR(2),COOR(3)
          IF(TEXT.EQ.'$') THEN
            GO TO 90
          END IF
        ELSE IF(KSRFC.LT.0) THEN
C         Input format 'Lines':
   20     CONTINUE
          IF(NEWLIN) THEN
            COOR(1)=GIANT
            COOR(2)=0.
            COOR(3)=0.
            TEXT='$'
            READ(LU1,*,END=90) TEXT,COOR(1),COOR(2),COOR(3)
            IF(TEXT.EQ.'$') THEN
              GO TO 90
            END IF
            NEWLIN=.FALSE.
            I=LENGTH(TEXT)+1
            TEXT(I:I)=''''
            WRITE(LU2,'(2A)')  '''',TEXT(1:I)
            IF(COOR(1).EQ.GIANT) THEN
C             No reference point:
              WRITE(LU2,'(A)') '/'
              GO TO 20
            ELSE
C             Check for the position of the reference point
              IF(BOUNDM(1).GT.COOR(1).OR.COOR(1).GT.BOUNDM(2).OR.
     *           BOUNDM(3).GT.COOR(2).OR.COOR(2).GT.BOUNDM(4).OR.
     *           BOUNDM(5).GT.COOR(3).OR.COOR(3).GT.BOUNDM(6)) THEN
                WRITE(LU2,'(A)') '/'
              END IF
            END IF
          ELSE
            COOR(1)=GIANT
            COOR(2)=0.
            COOR(3)=0.
            READ(LU1,*,END=90) COOR(1),COOR(2),COOR(3)
            IF(COOR(1).EQ.GIANT) THEN
C             End of line:
              NEWLIN=.TRUE.
              WRITE(LU2,'(A)') '/'
              GO TO 20
            END IF
          END IF
        ELSE
C         Input format for mixed surfaces:
          ISRFC=0
          READ(LU1,*,END=90) ISRFC,COOR(1),COOR(2),COOR(3)
          IF(ISRFC.EQ.0) THEN
            GO TO 90
          END IF
        END IF
        IF(BOUNDM(1).LE.COOR(1).AND.COOR(1).LE.BOUNDM(2)) THEN
          IF(BOUNDM(3).LE.COOR(2).AND.COOR(2).LE.BOUNDM(4)) THEN
            IF(BOUNDM(5).LE.COOR(3).AND.COOR(3).LE.BOUNDM(6)) THEN
              NPTS=NPTS+1
              CALL SRFC2(ISRFC,COOR,F)
              FMAX=AMAX1(ABS(F(1)),FMAX)
              FRMS=FRMS+F(1)*F(1)
              FABS=FABS+ ABS(F(1))
              FAVE=FAVE+     F(1)
              IF(KSRFC.GT.0) THEN
C               Output format 'Points':
                I=LENGTH(TEXT)+1
                TEXT(I:I)=''''
                I=MAX0(I,9)
                IF(1.LE.KOLUMN.AND.KOLUMN.LE.3) THEN
                  COOR(KOLUMN)=F(1)
                  WRITE(LU2,'(2A,3F12.6,A)')
     *                  '''',TEXT(1:I),COOR(1),COOR(2),COOR(3),     ' /'
                ELSE
                  WRITE(LU2,'(2A,4F12.6,A)')
     *                  '''',TEXT(1:I),COOR(1),COOR(2),COOR(3),F(1),' /'
                END IF
              ELSE IF(KSRFC.LT.0) THEN
C               Output format 'lines':
                IF(1.LE.KOLUMN.AND.KOLUMN.LE.3) THEN
                  COOR(KOLUMN)=F(1)
                  WRITE(LU2,'(   3F12.6,A)')
     *                                 COOR(1),COOR(2),COOR(3),     ' /'
                ELSE
                  WRITE(LU2,'(   4F12.6,A)')
     *                                 COOR(1),COOR(2),COOR(3),F(1),' /'
                END IF
              ELSE
C               Output format for mixed surfaces:
                IF(1.LE.KOLUMN.AND.KOLUMN.LE.3) THEN
                  COOR(KOLUMN)=F(1)
                  WRITE(LU2,'(I3,3F12.6,A)')
     *                           ISRFC,COOR(1),COOR(2),COOR(3)
                ELSE IF(KOLUMN.EQ.0) THEN
                  WRITE(LU2,'(I3,2F12.6)')               ISRFC,F(1)
                ELSE
                  WRITE(LU2,'(   4F12.6,A)')
     *                           ISRFC,COOR(1),COOR(2),COOR(3),F(1)
                END IF
              END IF
C             WRITE(*,'(A,I5,I3,2F12.6)') '+',NPTS,ISRFC,F(1),FMAX
            END IF
          END IF
        END IF
      GO TO 10
C
   90 CONTINUE
      IF(NPTS.GT.0) THEN
        FRMS=SQRT(FRMS/FLOAT(NPTS))
        FABS=     FABS/FLOAT(NPTS)
        FAVE=     FAVE/FLOAT(NPTS)
      END IF
      WRITE(LU2,'(A,I4,9(A,F10.6))') '/',NPTS,' POINTS, MAX=',FMAX,
     *                     ', RMS=',FRMS,', ABS=',FABS,', AVERAGE=',FAVE
      WRITE( * ,'(A,I4,9(A,F10.6))') '+',NPTS,' POINTS, MAX=',FMAX,
     *                     ', RMS=',FRMS,', ABS=',FABS,', AVERAGE=',FAVE
      STOP
      END
C
C=======================================================================
C
      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
      INCLUDE 'length.for'
C     length.for
C
C=======================================================================
C