C 
C Program INTF to check the positions of given points with respect to
C interfaces in the model.
C
C Version: 5.50
C Date: 2001, April 2
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.  Only the points within the model volume are
C             processed and written in this version.
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.  At most 98 columns are considered.
C             If the input file INTF has format LINES, the columns
C             are indexed from 1, in the other two cases the columns
C             are indexed from 0, the column 0 then contains either
C             the name of a point or the index of the surface and cannot
C             be modified.
C             If KOLUMN=0, input points are not modified.
C             For KSRFC=0, KOLUMN is allowed to be KOLUMN=-1, see the
C              description below.
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
C                                                 
C Output file INTFOUT has similar form as input file INTF, depending on
C KSRFC and KOLUMN:
C For KSRFC=0 and KOLUMN=-1:
C     Only the index of surface ISRFC and the value F(X1,X2,X3) of the
C     function describing the corresponding surface are written at
C     each line of the file.
C     If the point is situated exactly at the surface, F(X1,X2,X3)=0.
C For KOLUMN=0:
C     Output points situated inside the model volume coincide with the
C     input ones.  Points situated outside the model volume are not
C     written.
C Otherwise:
C     KOLUMN-th column of the input file is replaced by the value
C     F(X1,X2,X3) of the function describing the corresponding surface.
C     If the input file INTF contains less than (KOLUMN-1) values,
C     the missing values are replaced by the value of parameter ZERO,
C     see below.
C     For KOLUMN=1,2,3:
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 Note that only the points within the model volume are written to
C output file INTFOUT, which applies also for reference points of the
C lines if INTF is of format LINES.
C
C=======================================================================
C
C External procedures directly referred:
      EXTERNAL LENGTH,ERROR,RSEP1,RSEP3T,RSEP3I,MODEL1,SRFC2,FORM1
      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
      CHARACTER*20 FORMAT,FORMAR
      LOGICAL NEWLIN
      INTEGER KSRFC,ISRFC,NPTS,I,N,KOLUMN,J1
      REAL COOR(98),F(10),FMAX,FRMS,FABS,FAVE,UNDEF,ZERO,OUTMIN,OUTMAX
      PARAMETER (UNDEF=9.9E9, ZERO=0.)
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
      IF ((KOLUMN.LT.-1).OR.(KOLUMN.GT.98).OR.
     *    ((KOLUMN.EQ.-1).AND.(KSRFC.NE.0))) THEN
C       INTF-02
        CALL ERROR('INTF-02: Wrong value of KOLUMN')
C       Parameter KOLUMN may have values from 0 to 98 in case of
C       nonzero KSRFC, and from -1 to 98 if KSRFC=0.
      ENDIF
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
        COOR(1)=0.
        COOR(2)=0.
        COOR(3)=0.
        DO 12, I=4,98
          COOR(I)=UNDEF
  12    CONTINUE
        TEXT='$'
        IF(KSRFC.GT.0) THEN
C         Input format 'Points':
          READ(LU1,*,END=90) TEXT,(COOR(I),I=1,98)
          IF(TEXT.EQ.'$') THEN
            GO TO 90
          END IF
        ELSE IF(KSRFC.LT.0) THEN
C         Input format 'Lines':
          COOR(1)=UNDEF
          IF(NEWLIN) THEN
            READ(LU1,*,END=90) TEXT,(COOR(I),I=1,98)
            IF(TEXT.EQ.'$') THEN
              GO TO 90
            END IF
            NEWLIN=.FALSE.
            J1=LENGTH(TEXT)+1
            TEXT(J1:J1)=''''
            WRITE(LU2,'(2A)')  '''',TEXT(1:J1)
            IF(COOR(1).EQ.UNDEF) THEN
C             No reference point:
              WRITE(LU2,'(A)') '/'
              GO TO 10
            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
            READ(LU1,*,END=90) (COOR(I),I=1,98)
            IF(COOR(1).EQ.UNDEF) THEN
C             End of line:
              NEWLIN=.TRUE.
              WRITE(LU2,'(A)') '/'
              GO TO 10
            END IF
          END IF
        ELSE
C         (KSRFC=0) Input format for mixed surfaces:
          ISRFC=0
          READ(LU1,*,END=90) ISRFC,(COOR(I),I=1,98)
          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)
              N=MAX0(KOLUMN,0)
              DO 14, I=KOLUMN+1,98
                IF (COOR(I).NE.UNDEF) N=I
  14          CONTINUE
              OUTMIN=F(1)
              OUTMAX=F(1)
              DO 16, I=1,N
                IF (COOR(I).EQ.UNDEF) COOR(I)=ZERO
                IF((COOR(I).LT.OUTMIN).AND.(I.NE.KOLUMN)) OUTMIN=COOR(I)
                IF((COOR(I).GT.OUTMAX).AND.(I.NE.KOLUMN)) OUTMAX=COOR(I)
  16          CONTINUE
              FORMAR='F00.0,1X'
              CALL FORM1(OUTMIN,OUTMAX,FORMAR(1:8))
              FORMAR(7:8)='1X'
              IF(KSRFC.GT.0) THEN
C               Output format 'Points':
                J1=LENGTH(TEXT)+1
                TEXT(J1:J1)=''''
                J1=MAX0(J1,9)
                IF (KOLUMN.NE.0) THEN
                  COOR(KOLUMN)=F(1)
                ENDIF
                FORMAT='(2A,00(F00.0,1X),A)'
                FORMAT(6:6)=CHAR(ICHAR('0')+MOD(N/1,10))
                FORMAT(5:5)=CHAR(ICHAR('0')+MOD(N/10,10))
                FORMAT(8:15)=FORMAR(1:8)
                WRITE(LU2,FORMAT) '''',TEXT(1:J1),(COOR(I),I=1,N),' /'
              ELSE IF(KSRFC.LT.0) THEN
C               Output format 'lines':
                IF (KOLUMN.NE.0) THEN
                  COOR(KOLUMN)=F(1)
                ENDIF
                FORMAT='(00(F00.0,1X),A)'
                FORMAT(3:3)=CHAR(ICHAR('0')+MOD(N/1,10))
                FORMAT(2:2)=CHAR(ICHAR('0')+MOD(N/10,10))
                FORMAT(5:12)=FORMAR(1:8)
                WRITE(LU2,FORMAT) (COOR(I),I=1,N),' /'
              ELSE
C               (KSRFC=0) Output format for mixed surfaces:
                IF (KOLUMN.EQ.-1) THEN
                  FORMAT='(I3,1F00.0)'
                  FORMAT(6:10)=FORMAR(1:5)
                  WRITE(LU2,FORMAT) ISRFC,F(1)
                ELSE
                  IF (KOLUMN.NE.0) THEN
                    COOR(KOLUMN)=F(1)
                  ENDIF
                  FORMAT='(I3,00(F00.0,1X))'
                  FORMAT(6:6)=CHAR(ICHAR('0')+MOD(N/1,10))
                  FORMAT(5:5)=CHAR(ICHAR('0')+MOD(N/10,10))
                  FORMAT(8:15)=FORMAR(1:8)
                  WRITE(LU2,FORMAT) ISRFC,(COOR(I),I=1,N)
                END IF
              END IF
            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
      OUTMIN=AMIN1(FMAX,FRMS,FABS,FAVE)
      OUTMAX=AMAX1(FMAX,FRMS,FABS,FAVE)
      CALL FORM1(OUTMIN,OUTMAX,FORMAR(1:8))
      FORMAT='(A,I0,4(A,F00.0))'
      I=INT(ALOG10(FLOAT(NPTS)))+1
      IF (I.GT.9) THEN
C       INTF-03
        CALL ERROR('INTF-03: Too many points in file INTF')
C       This format specification allows for maximum of 100 000 000
C       of points in file INTF.
      ENDIF
      FORMAT(5:5)=CHAR(ICHAR('0')+I)
      FORMAT(11:15)=FORMAR(1:5)
      WRITE(LU2,FORMAT) '/',NPTS,' POINTS, MAX=',FMAX,
     *                     ', RMS=',FRMS,', ABS=',FABS,', AVERAGE=',FAVE
      WRITE( * ,FORMAT) '+',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 'forms.for'
C     forms.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