C 
C Program GRDPTS to generate the file containing the coordinates of all
C gridpoints of the given grid.
C
C Version: 5.50
C Date: 2001, January 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 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.
C             The file is not generated if PTS=' '.
C             Default: PTS='pts.out'
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 grid values.
C             KOLUMN=4:  The corresponding coordinate is written after
C               the three 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
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.
C      X4...  Optional grid value.
C (3) /
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
C     Filenames and parameters:
      CHARACTER*80 FILSEP,FGRD,FPTS,FPLGN,FTRGL
      INTEGER LU,LU2
      REAL UNDEF
      PARAMETER (LU=1,LU2=2,UNDEF=-999999999.)
C
C     Input data:
      INTEGER KOLUMN,N1,N2,N3
      REAL O1,O2,O3,D1,D2,D3
C
C     Other variables:
      CHARACTER*42 FORMAT
      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
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('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:
      IF(FPTS.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
        OPEN(LU,FILE=FPTS)
        WRITE(LU,'(A)') '/'
        FORMAT(1:10)='(A,I6.6,A,'
        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(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:
                CALL FORM2(NCOL,X,X,FORMAT(11:10+8*NCOL))
                WRITE(LU,FORMAT)
     *                        '''',I,''' ',X1,(' ',X(J4),J4=2,NCOL),' /'
                N=N+1
              END IF
   21       CONTINUE
   22     CONTINUE
   23   CONTINUE
        WRITE(LU,'(A)') '/'
        CLOSE(LU)
      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