C
C Program GRDPTS to generate the file containing the coordinates of all C gridpoints of the given grid. C C Version: 5.90 C Date: 2005, April 25 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 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 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 C Filenames and parameters: CHARACTER*80 FILSEP,FGRD,FLIN,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*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 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