C 
C Program to generate the file containing the coordinates of all
C gridpoints of the given grid.
C
C Version: 5.20
C Date: 1998, June 24
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                                                    
C Description of the data files:
C
C The data are read in by the list directed input (free format).
C In the description of data files, each numbered paragraph indicates
C the beginning of a new input operation (new READ statement).
C If the symbolic name of the input variable is enclosed in apostrophes,
C the corresponding value in input data is of the type CHARACTER, i.e.
C it should be a character string enclosed in apostrophes.  If the first
C letter of the symbolic name is I-N, the corresponding value is of the
C type INTEGER.  Otherwise, the input parameter is of the type REAL and
C may or may not contain a decimal point.
C
C Input data read from the * external unit:
C     The interactive * external unit may also be redirected to the file
C     containing the relevant data.
C (1) 'SEP','PTS',/
C     'SEP'...String in apostrophes containing the name of the input
C             file with the data specifying grid dimensions.
C             Description of file SEP
C     'PTS'...String in apostrophes containing the name of the output
C             file with the coordinates of the gridpoints.
C     /...    Input data line must be terminated by a slash.
C     Default: 'SEP'='grd.h', 'PTS'='pts.out'.
C
C                                                     
C Data file SEP has the form of the SEP (Stanford Exploration Project)
C parameter file:
C     All the data are specified in the form of PARAMETER=VALUE, e.g.
C     N1=50, with PARAMETER directly preceding = without intervening
C     spaces and with VALUE directly following = without intervening
C     spaces.  The PARAMETER=VALUE couple must be delimited by a space
C     or comma from both sides.
C     The PARAMETER string is not case-sensitive.
C     PARAMETER= followed by a space resets the default parameter value.
C     All other text in the input files is ignored.  The file thus may
C     contain unused data or comments without leading comment character.
C     Everything between comment character # and the end of the
C     respective line is ignored, too.
C     The PARAMETER=VALUE couples may be specified in any order.
C     The last appearance takes precedence.
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,/
C     'NNNNNN'...  Name of the point - six-digit integer index of the
C                  gridpoint (larger grids than 999999 gridpoints are
C                  not expected to be converted into this form suitable
C                  for a reasonably small number of points).
C      X1,X2,X3... Coordinates of the shifted point,
C                  X1=X10+X11*C1+X12*C2,
C                  X2=X20+X21*C1+X22*C2,
C                  X3=X30+X31*C1+X32*C2.
C (3) /
C
C-----------------------------------------------------------------------
C
      CHARACTER*80 FGRD,FPTS
      INTEGER LU
      PARAMETER (LU=1)
C
      CHARACTER*34 FORMAT
      INTEGER I1,I2,I3,I
      REAL X(3),X1,X2,X3,O1,O2,O3,D1,D2,D3
      EQUIVALENCE (X(1),X1),(X(2),X2),(X(3),X3)
C
C.......................................................................
C
C     Reading main input data:
      WRITE (*,'(2A)')
     *  ' Enter filenames of Grid header + Output points, /: '
      FGRD='grd.h'
      FPTS='pts.out'
      READ (*,*) FGRD,FPTS
C
C     Reading all the data from file FGRD to the memory
C     (SEP parameter file form):
      CALL RSEP1(LU,FGRD)
C
C     Recalling the data specifying grid dimensions
C     (arguments: Name of value in input data, Variable, Default):
      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:
      OPEN(LU,FILE=FPTS)
      WRITE(LU,'(A)') '/'
      FORMAT(1:10)='(A,I6.6,A,'
      I=0
      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
C           Writing:
            CALL FORM2(3,X,X,FORMAT(11:34))
            WRITE(LU,FORMAT) '''',I,''' ',X1,' ',X2,' ',X3,' /'
   21     CONTINUE
   22   CONTINUE
   23 CONTINUE
      WRITE(LU,'(A)') '/'
      CLOSE(LU)
C
      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