C 
C Program INTF to check the positions of given points with respect to
C interfaces in the model.
C
C Version: 5.40
C Date: 2000, May 29
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 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 input and output files:
C     MODEL='string'... Name of the input file with the data specifying
C             the model.
C             Description of MODEL
C             Example of MODEL
C             Default: MODEL='model.dat'
C     INTF='string'... Name of the input data file containing the
C             points situated at (or close to) interface(s).
C             Description of file INTF
C             Default: INTF='intf.dat'
C     INTFOUT='string'... 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             Description of file INTFOUT
C             Default: INTFOUT='intf.out'
C Data specific to this program:
C     KSRFC=integer ... Parameter describing the form of file INTF:
C             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
C               POINTS.
C             KSRFC.LT.0: Input file INTF has format
C               LINES.
C             Default: KSRFC=0
C     KOLUMN=integer ... Specifies the position in output file INTFOUT
C             where to write value F(X1,X2,X3) of the function
C             describing the corresponding surface, evaluated at the
C             given points.
C             Default: KOLUMN=4
C
C                                                    
C Input file INTF:
C For KSRFC.GT.0:
C     Input file INTF has format POINTS, see 'formsdat.htm'.
C     Description of form POINTS
C For KSRFC.LT.0:
C     Input file INTF has format LINES, see 'formsdat.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                                                 
C Output file INTFOUT has similar form as input file INTF, depending on
C KSRFC and KOLUMN:
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
      CHARACTER*80 FILSEP
      INTEGER LU0
      PARAMETER (LU0=1)
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     Reading name of SEP file with input data:
      WRITE(*,'(A)') '+INTF: Enter input filename: '
      FILSEP=' '
      READ(*,*) FILSEP
      WRITE(*,'(A)') '+INTF: Working ...           '
C
C     Reading all data from the SEP file into the memory:
      IF (FILSEP.NE.' ') THEN
        CALL RSEP1(LU0,FILSEP)
      ELSE
C       INTF-01
        CALL ERROR('INTF-01: 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',FILE0,'model.dat')
      CALL RSEP3T('INTF',FILE1,'intf.dat')
      CALL RSEP3T('INTFOUT',FILE2,'intf.out')
      CALL RSEP3I('KSRFC',KSRFC,0)
      CALL RSEP3I('KOLUMN',KOLUMN,4)
C
      OPEN(LU1,FILE=FILE0,STATUS='OLD')
      CALL MODEL1(LU1)
      CLOSE(LU1)
C
      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)')
     *                           ISRFC,COOR(1),COOR(2),COOR(3)
                ELSE IF(KOLUMN.EQ.0) THEN
                  WRITE(LU2,'(I3,2F12.6)')               ISRFC,F(1)
                ELSE
                  WRITE(LU2,'(I3,4F12.6)')
     *                           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
      WRITE(*,'(A)') ' INTF: Done.'
      STOP
      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