C
C Program GRDPTS to generate the file containing the coordinates of all
C gridpoints of the given grid.
C
C Version: 6.00
C Date: 2006, March 3
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 Name of the optional input file:
C     GRD='string'... String with the name of the input data file
C             containing the grid values.
C             This file is used only if KOLUMN.NE.0, see below.
C             The grid 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 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 ... 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:  The grid values are written after the three
C               coordinates, into the fourth numeric column.
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 in 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,/
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...  Optional grid value.
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,/ or
C       'NNNNNN',X1,X2,X3,X4,/
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      X4...  Optional grid value in the first gridpoint of the line.
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,/
C      X1,X2,X3... Coordinates of the gridpoint of the gridline,
C             or the value at the gridpoint, see parameter KOLUMN.
C      X4...  Optional grid value in the gridpoint of the 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,FLIN,FPTS,FPLGN,FTRGL
      INTEGER LU,LU2
      REAL UNDEF
      PARAMETER (LU=1,LU2=2)
C
C     Input data:
      INTEGER KOLUMN,N1,N2,N3
      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
      REAL X(4),X1,X2,X3,X4
      EQUIVALENCE (X(1),X1),(X(2),X2),(X(3),X3),(X(4),X4)
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 ,' ')
      CALL RSEP3T('PTS' ,FPTS ,'pts.out')
      CALL RSEP3T('LIN' ,FLIN,' ')
      CALL RSEP3T('PLGN',FPLGN,' ')
      CALL RSEP3T('TRGL',FTRGL,' ')
      CALL RSEP3I('KOLUMN',KOLUMN,0)
      IF (KOLUMN.LT.0.OR.4.LT.KOLUMN) THEN
C       GRDPTS-02
        CALL ERROR('GRDPTS-02: Wrong value of KOLUMN')
C       KOLUMN must be 0, 1, 2, 3 or 4.
      ENDIF
      IF (KOLUMN.NE.0.AND.FGRD.EQ.' ') THEN
C       GRDPTS-03
        CALL ERROR('GRDPTS-03: 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)
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.)
C
C     Writing output points or lines:
      IF((FPTS.NE.' ').OR.(FLIN.NE.' ')) THEN
        IF(KOLUMN.NE.0) THEN
          IF(N1*N2*N3.GT.MRAM) THEN
C           GRDPTS-04
            CALL ERROR('GRDPTS-04: 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.
C           You may wish to increase the dimension MRAM in file
C           ram.inc.
          END IF
C         Reading grid values:
          CALL RARRAY(LU,FGRD,'FORMATTED',.TRUE.,UNDEF,N1*N2*N3,RAM)
        END IF
        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
              IF(KOLUMN.EQ.0) THEN
                LWRITE=.TRUE.
              ELSE
                IF(RAM(I).NE.UNDEF) THEN
                  X(KOLUMN)=RAM(I)
                  IRAM(I)=N
                  LWRITE=.TRUE.
                ELSE
                  IRAM(I)=0
                  LWRITE=.FALSE.
                END IF
              END IF
              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(KOLUMN.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(KOLUMN.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(KOLUMN.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