C
C Program GRDPTS to generate the file containing the coordinates of all
C gridpoints of the given grid.
C
C Version: 6.20
C Date: 2008, June 12
C
C Coded by: Ludek Klimes
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     http://sw3d.cz/staff/klimes.htm
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 Name of the optional input files:
C     GRD='string',GRD1='string',GRD2='string',...,GRD9='string'...
C             Strings with the names of the input data files containing
C             the grid values.
C             These files are used only if KOLUMN.NE.0, KOLUMN1.NE.0,
C             KOLUMN2.NE.0, ..., KOLUMN9.NE.0, respectively, see below.
C             The grids may also contain undefined values. In this case,
C             if KOLUMN*.NE.0, the gridpoints with undefined values are
C             not written to output file PTS.
C             If the filename is blank, all gridded values are assumed
C             undefined.
C             No default, GRD must be specified and cannot be blank if
C             KOLUMN*.NE.0.
C Names of the output files:
C     PTS='string'...Name of the output file with the coordinates
C             of the gridpoints written in the form
C             POINTS.
C             The file is not generated if PTS=' '.
C             Default: PTS='pts.out'
C     LIN='string'...Name of the optional output file with the
C             coordinates of the gridpoints written in the form
C             LINES. The N2*N3 lines
C             is generated, each line consisting of N1 points.
C             Default: LIN=' ' (no file generated).
C     PLGN='string'... Name of the optional output file specifying the
C             polygons coinciding with grid faces.
C             The polygons are described in terms of the indices of
C             their vertices.
C             The file is not generated if PLGN=' ' (default).
C             For 3-D virtual reality, the polygons should be divided
C             into triangles by program
C             trgl.for.
C             Since the vertices have no normals, the current version of
C             program trgl.for does not work with them and parameter
C             TRGL should be used instead of PLGN.
C             Description of file PLGN
C             Default: PLGN=' '
C     TRGL='string'... Name of the optional output file specifying the
C             triangles covering grid faces.
C             The triangles are described in terms of the indices of
C             their vertices.
C             The file is not generated if TRGL=' ' (default).
C             Description of file TRGL
C             Default: TRGL=' '
C Data specifying the form of the output file with points:
C     KOLUMN=integer,KOLUMN1=integer,KOLUMN2=integer,...,
C     KOLUMN9=integer... If non-zero, specifies the position in output
C             file PTS where to write the values at gridpoints.
C             KOLUMN*=0:  Grid values are not written to file PTS.
C             KOLUMN*=1, 2 or 3:  The corresponding coordinate is
C               replaced with the grid values.
C             KOLUMN*=4,5,...:  The grid values are written after the
C               three coordinates, into the specified numeric column.
C             KOLUMN*=-1, -2 or -3:  The grid values are added to
C               the corresponding coordinate.
C             Default: KOLUMN*=0
C Data specifying grid dimensions:
C     N1=positive integer... Number of gridpoints along the X1 axis.
C             Default: N1=1
C     N2=positive integer... Number of gridpoints along the X2 axis.
C             Default: N2=1
C     N3=positive integer... Number of gridpoints along the X3 axis.
C             Default: N3=1
C     O1=real... First coordinate of the grid origin (first point of the
C             grid).
C             Default: O1=0.
C     O2=real... Second coordinate of the grid origin.
C             Default: O2=0.
C     O3=real... Third coordinate of the grid origin.
C             Default: O3=0.
C     D1=real... Grid interval in the direction of the first coordinate
C             axis.
C             Default: D1=1.
C     D2=real... Grid interval in the direction of the second coordinate
C             axis.
C             Default: D2=1.
C     D3=real... Grid interval in the direction of the third coordinate
C             axis.
C             Default: D3=1.
C Optional parameters specifying the form of the real quantities
C written to the output formatted files:
C     MINDIG,MAXDIG=positive integers ... See the description in file
C             forms.for.
C
C                                                     
C Output file PTS with the gridpoints:
C (1) /
C (2) For each gridpoint data (2.1):
C (2.1) 'NNNNNN',X1,X2,X3,/ or
C       'NNNNNN',X1,X2,X3,X4,/ or
C       'NNNNNN',X1,X2,X3,X4,X5,/ etc.
C     'NNNNNN'...  Name of the point - six-digit integer index of the
C             gridpoint (larger grids than 999999 gridpoints are not
C             expected to be converted into this form suitable for
C             a reasonably small number of points).
C      X1,X2,X3... Coordinates of the gridpoint, or the value at
C             the gridpoint, see parameter KOLUMN.
C      X4,X5,...,X9...  Optional grid values.
C (3) /
C
C                                                     
C Optional output file LIN with the gridlines:
C (1) /
C (2) For each of the N2*N3 gridlines data (2.1) and (2.2):
C (2.1) 'NNNNNN',X1,X2,X3,/
C     'NNNNNN'...  Name of the line - six-digit integer index of the
C             gridline, the first three digits correspond to the index
C             along the second axis, and the second three digits
C             correspond to the third axis (larger grids than N2=999
C             and N3=999 are not expected to be converted into this form
C             suitable for a reasonably small number of points).
C      X1,X2,X3... Coordinates of the first gridpoint of the gridline,
C             or the value at the gridpoint, see parameter KOLUMN.
C (2.2) Points of the line - for each point of the line (2.2.1):
C (2.2.1) X1,X2,X3,/ or
C         X1,X2,X3,X4,/ or
C         X1,X2,X3,X4,X5,/ etc.
C      X1,X2,X3... Coordinates of the gridpoint of the gridline,
C             or the value at the gridpoint, see parameter KOLUMN.
C      X4,X5,...,X9... Optional grid values at the gridpoint of the
C             line.
C (2.3) / (a slash indicating end of the line).
C (3) / (a slash indicating end of file).
C
C                                                    
C Optional output file PLGN with the polygons:
C (1) For each grid face data (1.1):
C (1.1) I1,I2,I3,I4,/
C     I1,I2,I3,I4... Integer indices of 4 vertices of the rectangle
C             forming the grid face.
C             The vertices are stored in file PTS and are indexed by
C             positive integers according to their order.
C     /...    List of vertices is terminated by a slash.
C
C                                                    
C Optional output file TRGL with the triangles:
C (1) For each triangle data (1.1):
C (1.1) I1,I2,I3,/
C     I1,I2,I3... Integer indices of 3 vertices of one of two triangles
C             forming the grid face.
C             The vertices are stored in file PTS and are indexed by
C             positive integers according to their order.
C     /...    List of vertices is terminated by a slash.
C
C=======================================================================
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
      INTEGER IRAM(MRAM)
      EQUIVALENCE (IRAM,RAM)
C
C.......................................................................
C
      EXTERNAL UARRAY
      REAL UARRAY
C     Filenames and parameters:
      CHARACTER*80 FILSEP,FGRD(0:9),FLIN,FPTS,FPLGN,FTRGL
      INTEGER LU,LU2
      REAL UNDEF
      PARAMETER (LU=1,LU2=2)
C
C     Input data:
      INTEGER KOLUMN(0:9),N1,N2,N3,N1N2N3
      REAL O1,O2,O3,D1,D2,D3
C
C     Other variables:
      CHARACTER*47 FORMA1,FORMA2,FORMA3
      LOGICAL LWRITE
      INTEGER NCOL,N,I,I1,I2,I3,J1,J2,J3,J4,K,KCOL
      REAL X(9),X1,X2,X3
      EQUIVALENCE (X(1),X1),(X(2),X2),(X(3),X3)
C
      UNDEF=UARRAY()
C
C-----------------------------------------------------------------------
C
C     Reading input SEP parameter file:
      WRITE(*,'(A)') '+GRDPTS: Enter input filename: '
      FILSEP=' '
      READ(*,*) FILSEP
      IF (FILSEP.EQ.' ') THEN
C       GRDPTS-01
        CALL ERROR('GRDPTS-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
      CALL RSEP1(LU,FILSEP)
      WRITE(*,'(A)') '+GRDPTS: Working ...           '
C
C     Reading input parameters from the SEP file:
      CALL RSEP3T('GRD' ,FGRD(0),' ')
      CALL RSEP3T('GRD1',FGRD(1),' ')
      CALL RSEP3T('GRD2',FGRD(2),' ')
      CALL RSEP3T('GRD3',FGRD(3),' ')
      CALL RSEP3T('GRD4',FGRD(4),' ')
      CALL RSEP3T('GRD5',FGRD(5),' ')
      CALL RSEP3T('GRD6',FGRD(6),' ')
      CALL RSEP3T('GRD7',FGRD(7),' ')
      CALL RSEP3T('GRD8',FGRD(8),' ')
      CALL RSEP3T('GRD9',FGRD(9),' ')
      CALL RSEP3T('PTS' ,FPTS ,'pts.out')
      CALL RSEP3T('LIN' ,FLIN,' ')
      CALL RSEP3T('PLGN',FPLGN,' ')
      CALL RSEP3T('TRGL',FTRGL,' ')
      CALL RSEP3I('KOLUMN' ,KOLUMN(0),0)
      CALL RSEP3I('KOLUMN1',KOLUMN(1),0)
      CALL RSEP3I('KOLUMN2',KOLUMN(2),0)
      CALL RSEP3I('KOLUMN3',KOLUMN(3),0)
      CALL RSEP3I('KOLUMN4',KOLUMN(4),0)
      CALL RSEP3I('KOLUMN5',KOLUMN(5),0)
      CALL RSEP3I('KOLUMN6',KOLUMN(6),0)
      CALL RSEP3I('KOLUMN7',KOLUMN(7),0)
      CALL RSEP3I('KOLUMN8',KOLUMN(8),0)
      CALL RSEP3I('KOLUMN9',KOLUMN(9),0)
      IF (KOLUMN(0).NE.0.AND.FGRD(0).EQ.' '.OR.
     *    KOLUMN(1).NE.0.AND.FGRD(1).EQ.' '.OR.
     *    KOLUMN(2).NE.0.AND.FGRD(2).EQ.' '.OR.
     *    KOLUMN(3).NE.0.AND.FGRD(3).EQ.' '.OR.
     *    KOLUMN(4).NE.0.AND.FGRD(4).EQ.' '.OR.
     *    KOLUMN(5).NE.0.AND.FGRD(5).EQ.' '.OR.
     *    KOLUMN(6).NE.0.AND.FGRD(6).EQ.' '.OR.
     *    KOLUMN(7).NE.0.AND.FGRD(7).EQ.' '.OR.
     *    KOLUMN(8).NE.0.AND.FGRD(8).EQ.' '.OR.
     *    KOLUMN(9).NE.0.AND.FGRD(9).EQ.' ') THEN
C       GRDPTS-02
        CALL ERROR('GRDPTS-02: File GRD not specified')
C       Input file GRD with gridded values must be specified if
C       KOLUMN*.NE.0.
C       There is no default filename.
      ENDIF
      NCOL=MAX0(3,KOLUMN(0),KOLUMN(1),KOLUMN(2),KOLUMN(3),KOLUMN(4),
     *            KOLUMN(5),KOLUMN(6),KOLUMN(7),KOLUMN(8),KOLUMN(9))
      K   =MIN0(0,KOLUMN(0),KOLUMN(1),KOLUMN(2),KOLUMN(3),KOLUMN(4),
     *            KOLUMN(5),KOLUMN(6),KOLUMN(7),KOLUMN(8),KOLUMN(9))
      KCOL=MAX0(0,KOLUMN(0),KOLUMN(1),KOLUMN(2),KOLUMN(3),KOLUMN(4),
     *            KOLUMN(5),KOLUMN(6),KOLUMN(7),KOLUMN(8),KOLUMN(9))-K
      IF(K.LT.-3.OR.NCOL.GT.9) THEN
C       GRDPTS-03
        CALL ERROR('GRDPTS-03: Wrong value of KOLUMN*')
C       KOLUMN* must be -3,-2,-1,0,1,2,3,...,9.
      ENDIF
      DO 11 K=1,NCOL
        X(K)=0.
   11 CONTINUE
      DO 12 K=0,9
        IF(KOLUMN(K).NE.0) THEN
          X(IABS(KOLUMN(K)))=X(IABS(KOLUMN(K)))+1.
        ENDIF
   12 CONTINUE
      DO 13 K=1,NCOL
        IF(X(K).GT.1.5) THEN
C         GRDPTS-04
          CALL ERROR('GRDPTS-04: Two equal values of KOLUMN*')
C         Two different grids cannot correspond to the same column.
        ENDIF
   13 CONTINUE
      DO 14 K=4,NCOL
        IF(X(K).LT.0.5) THEN
C         GRDPTS-05
          CALL ERROR('GRDPTS-05: Empty column in output files')
C         There is a column corresponding to no grid.
C         The values of KOLUMN* cannot form a gap.
        ENDIF
   14 CONTINUE
C     Data specifying grid dimensions
      CALL RSEP3I('N1',N1,1)
      CALL RSEP3I('N2',N2,1)
      CALL RSEP3I('N3',N3,1)
      CALL RSEP3R('O1',O1,0.)
      CALL RSEP3R('O2',O2,0.)
      CALL RSEP3R('O3',O3,0.)
      CALL RSEP3R('D1',D1,1.)
      CALL RSEP3R('D2',D2,1.)
      CALL RSEP3R('D3',D3,1.)
      N1N2N3=N1*N2*N3
C
C     Writing output points or lines:
      IF((FPTS.NE.' ').OR.(FLIN.NE.' ')) THEN
        DO 19 K=0,9
          IF(KOLUMN(K).NE.0) THEN
            IF((K+1)*N1N2N3.GT.MRAM) THEN
C             GRDPTS-06
              CALL ERROR('GRDPTS-06: Too small array RAM(MRAM)')
C             Array RAM(MRAM) allocated in include file 'ram.inc' is too
C             small to contain N1*N2*N3 grid values for each given grid.
C             You may wish to increase the dimension MRAM in file
C             ram.inc.
            END IF
C           Reading grid values:
            CALL RARRAY(LU,FGRD(K),'FORMATTED',.TRUE.,UNDEF,N1N2N3,
     *                  RAM(K*N1N2N3+1))
          END IF
   19   CONTINUE
        IF(FPTS.NE.' ') THEN
          OPEN(LU,FILE=FPTS)
          WRITE(LU,'(A)') ' /'
          FORMA1(1:10)='(A,I6.6,A,'
        ENDIF
        IF(FLIN.NE.' ') THEN
          OPEN(LU2,FILE=FLIN)
          WRITE(LU2,'(A)') ' /'
          FORMA2(1:15)='(A,I3.3,I3.3,A,'
          FORMA3(1:1)='('
        ENDIF
        I=0
        N=1
        DO 23 I3=0,N3-1
          X3=O3+D3*FLOAT(I3)
          DO 22 I2=0,N2-1
            X2=O2+D2*FLOAT(I2)
            DO 21 I1=0,N1-1
              X1=O1+D1*FLOAT(I1)
              I=I+1
              IF ((I1.EQ.0).AND.(FLIN.NE.' ')) THEN
                CALL FORM2(3,X,X,FORMA2(16:15+8*3))
                WRITE(LU2,FORMA2)
     *                   '''',I2+1,I3+1,''' ',X1,(' ',X(J4),J4=2,3),' /'
              ENDIF
              LWRITE=.TRUE.
              DO 20 K=0,9
                IF(KOLUMN(K).GT.0) THEN
                  IF(RAM(K*N1N2N3+I).NE.UNDEF) THEN
                    X(KOLUMN(K))=RAM(K*N1N2N3+I)
                    IF(K.EQ.0) THEN
                      IRAM(I)=N
                    END IF
                  ELSE
                    IRAM(I)=0
                    LWRITE=.FALSE.
                  END IF
                ELSE IF(KOLUMN(K).LT.0) THEN
                  IF(RAM(-K*N1N2N3+I).NE.UNDEF) THEN
                    X(-KOLUMN(K))=X(-KOLUMN(K))+RAM(-K*N1N2N3+I)
                    IF(K.EQ.0) THEN
                      IRAM(I)=N
                    END IF
                  ELSE
                    IRAM(I)=0
                    LWRITE=.FALSE.
                  END IF
                END IF
   20         CONTINUE
              IF(LWRITE) THEN
C               Writing:
                IF(FPTS.NE.' ') THEN
                  CALL FORM2(NCOL,X,X,FORMA1(11:10+8*NCOL))
                  WRITE(LU,FORMA1)
     *                        '''',I,''' ',X1,(' ',X(J4),J4=2,NCOL),' /'
                ENDIF
                IF(FLIN.NE.' ') THEN
                  CALL FORM2(NCOL,X,X,FORMA3(2:1+8*NCOL))
                  WRITE(LU2,FORMA3)
     *                                     X1,(' ',X(J4),J4=2,NCOL),' /'
                ENDIF
                N=N+1
              END IF
              IF ((I1.EQ.N1-1).AND.(FLIN.NE.' ')) THEN
                WRITE(LU2,'(A)') ' /'
              ENDIF
   21       CONTINUE
   22     CONTINUE
   23   CONTINUE
        IF(FPTS.NE.' ') THEN
          WRITE(LU,'(A)') ' /'
          CLOSE(LU)
        ENDIF
        IF(FLIN.NE.' ') THEN
          WRITE(LU2,'(A)') ' /'
          CLOSE(LU2)
        ENDIF
      END IF
C
C     Writing output polygons:
      IF(FPLGN.NE.' ') THEN
        OPEN(LU,FILE=FPLGN)
      END IF
      IF(FTRGL.NE.' ') THEN
        OPEN(LU2,FILE=FTRGL)
      END IF
      IF(FPLGN.NE.' '.OR.FTRGL.NE.' ') THEN
        DO 33 I3=0,N3-1
          DO 32 I2=0,N2-2
            DO 31 I1=0,N1-2
              I=1+I1+N1*(I2+N2*I3)
              J1=I
              J2=I+1
              J3=I+N1+1
              J4=I+N1
              IF(KCOL.NE.0) THEN
                J1=IRAM(J1)
                J2=IRAM(J2)
                J3=IRAM(J3)
                J4=IRAM(J4)
              END IF
              IF(FPLGN.NE.' ') THEN
                IF(J1.NE.0.AND.J2.NE.0.AND.J3.NE.0.AND.J4.NE.0) THEN
                  WRITE(LU,'(4(I6,A))') J1,' ',J2,' ',J3,' ',J4,' /'
                END IF
              END IF
              IF(FTRGL.NE.' ') THEN
                IF(J1.NE.0.AND.J2.NE.0.AND.J3.NE.0) THEN
                  WRITE(LU2,'(3(I6,A))') J1,' ',J2,' ',J3,' /'
                END IF
                IF(J1.NE.0.AND.J3.NE.0.AND.J4.NE.0) THEN
                  WRITE(LU2,'(3(I6,A))') J1,' ',J3,' ',J4,' /'
                END IF
              END IF
   31       CONTINUE
   32     CONTINUE
   33   CONTINUE
        DO 43 I3=0,N3-2
          DO 42 I2=0,N2-1
            DO 41 I1=0,N1-2
              I=1+I1+N1*(I2+N2*I3)
              J1=I
              J2=I+1
              J3=I+N1*N2+1
              J4=I+N1*N2
              IF(KCOL.NE.0) THEN
                J1=IRAM(J1)
                J2=IRAM(J2)
                J3=IRAM(J3)
                J4=IRAM(J4)
              END IF
              IF(FPLGN.NE.' ') THEN
                IF(J1.NE.0.AND.J2.NE.0.AND.J3.NE.0.AND.J4.NE.0) THEN
                  WRITE(LU,'(4(I6,A))') J1,' ',J2,' ',J3,' ',J4,' /'
                END IF
              END IF
              IF(FTRGL.NE.' ') THEN
                IF(J1.NE.0.AND.J2.NE.0.AND.J3.NE.0) THEN
                  WRITE(LU2,'(3(I6,A))') J1,' ',J2,' ',J3,' /'
                END IF
                IF(J1.NE.0.AND.J3.NE.0.AND.J4.NE.0) THEN
                  WRITE(LU2,'(3(I6,A))') J1,' ',J3,' ',J4,' /'
                END IF
              END IF
   41       CONTINUE
   42     CONTINUE
   43   CONTINUE
        DO 53 I3=0,N3-2
          DO 52 I2=0,N2-2
            DO 51 I1=0,N1-1
              I=1+I1+N1*(I2+N2*I3)
              J1=I
              J2=I+N1
              J3=I+N1*N2+N1
              J4=I+N1*N2
              IF(KCOL.NE.0) THEN
                J1=IRAM(J1)
                J2=IRAM(J2)
                J3=IRAM(J3)
                J4=IRAM(J4)
              END IF
              IF(FPLGN.NE.' ') THEN
                IF(J1.NE.0.AND.J2.NE.0.AND.J3.NE.0.AND.J4.NE.0) THEN
                  WRITE(LU,'(4(I6,A))') J1,' ',J2,' ',J3,' ',J4,' /'
                END IF
              END IF
              IF(FTRGL.NE.' ') THEN
                IF(J1.NE.0.AND.J2.NE.0.AND.J3.NE.0) THEN
                  WRITE(LU2,'(3(I6,A))') J1,' ',J2,' ',J3,' /'
                END IF
                IF(J1.NE.0.AND.J3.NE.0.AND.J4.NE.0) THEN
                  WRITE(LU2,'(3(I6,A))') J1,' ',J3,' ',J4,' /'
                END IF
              END IF
   51       CONTINUE
   52     CONTINUE
   53   CONTINUE
      END IF
      IF(FPLGN.NE.' ') THEN
        CLOSE(LU)
      END IF
      IF(FTRGL.NE.' ') THEN
        CLOSE(LU2)
      END IF
C
      WRITE(*,'(A)') '+GRDPTS: 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
C
C=======================================================================
C