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