C
C Program GRDISO for identification of points at isosurfaces C of 3-D gridded values. If one of two neigbouring gridpoints displays C the value greater than VALUE, and if the second one of the two C gridpoints displays the value less or eaqual to VALUE, the coordinates C of the point which lies between the two gridpoints and corresponds C to the value VALUE are computed and stored in the output file. C C Version: 5.50 C Date: 2001, April 5 C C Coded by Petr Bulant C Department of Geophysics, Charles University Prague, C Ke Karlovu 3, 121 16 Praha 2, Czech Republic, C E-mail: bulant@seis.karlov.mff.cuni.cz 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 input file: C GRD='string'... Name of the input ASCII file with the grid values. C Default: GRD='grd.out' C For general description of the files with gridded data refer C to file forms.htm. C Parameters defining the basic regular rectangular grid: C N1=positive integer... Number of gridpoints of the basic grid C along the X1 axis. Default: N1=1 C N2=positive integer... Number of gridpoints of the basic grid C along the X2 axis. Default: N2=1 C N3=positive integer... Number of gridpoints of the basic grid C along the X3 axis. Default: N3=1 C O1=real... X1 coordinate of the origin of the grid. Default: O1=0. C O2=real... X2 coordinate of the origin of the grid. Default: O2=0. C O3=real... X3 coordinate of the origin of the grid. Default: O3=0. C D1=real... Grid spacing along the X1 axis. Default: D1=1. C D2=real... Grid spacing along the X2 axis. Default: D2=1. C D3=real... Grid spacing along the X3 axis. Default: D3=1. C Data specifying dimensions of the grid, on which the calculation C of points at isosurfaces will be performed: C N1NEW=positive integer... Number of gridpoints along the X1 axis. C Default: N1NEW=N1 C N2NEW=positive integer... Number of gridpoints along the X2 axis. C Default: N2NEW=N2 C N3NEW=positive integer... Number of gridpoints along the X3 axis. C Default: N3NEW=N3 C NO1=positive integer... Index of the first gridpoint along C the X1 axis. C Default: NO1=1 C NO2=positive integer... Index of the first gridpoint along C the X2 axis. C Default: NO2=1 C NO3=positive integer... Index of the first gridpoint along C the X3 axis. C Default: NO3=1 C ND1=positive integer... Multiplication factor of the grid interval C along the X1 axis. C Default: ND1=1 C ND2=positive integer... Multiplication factor of the grid interval C along the X2 axis. C Default: ND2=1 C ND3=positive integer... Multiplication factor of the grid interval C along the X3 axis. C Default: ND3=1 C The grid for calculation should be always a subgrid of the original C grid. Two gridpoints located at the different sides of the C isosurface are searched for on the grid for calculation. C The coordinates of the point at isosurface are then calculated on C the original grid. C Value corresponding to the points to be calculated: C VALUE=real ... Value corresponding to the isosurface being C searched for. C Default: VALUE=0. C Mode of the calculation: C MODE=integer C MODE=0 ... The points corresponding to VALUE are C searched for at all gridlegs. C MODE=1,2,3 ... The points are searched for only C at the gridlegs parallel with X1, X2 or X3 axis. C All the points are written to the file PTS. C MODE=-1,-2,-3 ... The points are searched for only C at the gridlegs parallel with X1, X2 or X3 axis. C The points with values increasing with C the corresponding gridlegs are written to the file PTS, C other points to the file PTS1. C Default: MODE=0 C Names of the output files: C PTS='string'...Name of the output file with the coordinates C of the points. The points are the intersections of C the gridlegs with the isosurface of the value VALUE. C The file is not generated if PTS=' '. C Description of file PTS. C Default: PTS='pts.out' C PTS1='string'...Name of the second output file with the points C in case of negative MODE. C The file is not generated if PTS1=' '. C Description of file PTS1. C Default: PTS1='pts1.out' C C C Output file PTS or PTS1 with the points at isosurface: C (1) / C (2) For each point data (2.1): C (2.1) 'NNNNNN',X1,X2,X3,/ C 'NNNNNN'... Name of the point - six-digit integer index of the C point. C X1,X2,X3... Coordinates of the point. C (3) / C C----------------------------------------------------------------------- C Subroutines and external functions required: EXTERNAL GIGLEG,GICP, *ERROR,RSEP1,RSEP3T,RSEP3R,RSEP3I,FORM1 C GIGLEG,GICP ... This file. C ERROR ... File error.for. C RSEP1,RSEP3T,RSEP3R,RSEP3I ... C File sep.for. C FORM1 ... forms.for. C C Common block /GIC/: INCLUDE 'grdiso.inc' C grdiso.inc C C....................................................................... C Auxiliary storage locations: CHARACTER*80 FSEP,FGRD,FPTS,FPTS1 INTEGER LU,LU1,LENG,MODE PARAMETER (LU=1,LU1=2,LENG=6) REAL UNDEF,VALUE PARAMETER (UNDEF=-999999.) INTEGER I1,I2,I3,I4,J2,JCOOR,IP1,IP2,IP3,NPTS, * NLEG2,NLEG3 REAL XX(3),YY(3),ZZ(3),OUTMIN,OUTMAX,VAL1,VAL2,VAL3 CHARACTER*20 FORMAT,TEXTP C C....................................................................... C C Reading name of SEP file with input data: WRITE(*,'(A)') '+GRDISO: Enter input filename: ' FSEP=' ' READ(*,*) FSEP C C Reading all data from the SEP file into the memory: IF (FSEP.EQ.' ') THEN C GRDISO-01 CALL ERROR('GRDISO-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 C C Reading all the data from file FSEP to the memory C (SEP parameter file form): CALL RSEP1(LU,FSEP) 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',DD(1),1.) CALL RSEP3R('D2',DD(2),1.) CALL RSEP3R('D3',DD(3),1.) CALL RSEP3I('N1NEW',N1N,N1) CALL RSEP3I('N2NEW',N2N,N2) CALL RSEP3I('N3NEW',N3N,N3) IF (N1.EQ.N1N.AND.N2.EQ.N2N.AND.N3.EQ.N3N) THEN NO1=1 NO2=1 NO3=1 ND1=1 ND2=1 ND3=1 ELSE CALL RSEP3I('NO1',NO1,1) CALL RSEP3I('NO2',NO2,1) CALL RSEP3I('NO3',NO3,1) CALL RSEP3I('ND1',ND1,1) CALL RSEP3I('ND2',ND2,1) CALL RSEP3I('ND3',ND3,1) IF ((NO1+ND1*(N1N-1).GT.N1).OR. * (NO2+ND2*(N2N-1).GT.N2).OR. * (NO3+ND3*(N3N-1).GT.N3)) THEN C GRDISO-02 CALL ERROR('GRDISO-02: Wrong grid for calculation') C The grid for calculation must be a subgrid C of the original grid. ENDIF ENDIF C C Mode of the calculation and the isovalue: CALL RSEP3I('MODE',MODE,0) CALL RSEP3R('VALUE',VALUE,0.) C C Recalling the input and output filenames: CALL RSEP3T('GRD',FGRD,'grd.out') CALL RSEP3T('PTS',FPTS,'pts.out') IF (MODE.LT.0) CALL RSEP3T('PTS1',FPTS1,'pts1.out') C C C Reading input grid values: IF (N1*N2*N3.GT.MRAM) THEN C GRDISO-03 CALL ERROR('GRDISO-03: Small array RAM.') C Try to enlarge the dimension MRAM in the file C ram.inc. ENDIF CALL RARRAY(LU,FGRD,'FORMATTED',.TRUE.,UNDEF,N1*N2*N3,RAM) NPTS=0 C C Opening output files, initializing: WRITE(*,'(A)') '+GRDISO: Working ... ' IF (FPTS.NE.' ') THEN OPEN(LU,FILE=FPTS,FORM='FORMATTED') WRITE(LU,'(A)') '/' ENDIF IF ((MODE.LT.0).AND.(FPTS1.NE.' ')) THEN OPEN(LU1,FILE=FPTS1,FORM='FORMATTED') WRITE(LU1,'(A)') '/' ENDIF C N1N2= N1 * N2 C NLEG1=(N1N-1)*N2N*N3N NLEG2=N1N*(N2N-1)*N3N NLEG3=N1N*N2N*(N3N-1) NLEG12=NLEG1+NLEG2 NLEG=NLEG12+NLEG3 C NN11= N1N-1 NN1N2= N1N * N2N NN11N2= (N1N-1)* N2N NN1N21= N1N *(N2N-1) C DIP - number of points in RAM between neighboring points C in calculation grid DIP(1)=1 DIP(2)=N1 DIP(3)=N1*N2 C I1 - index of the gridleg being processed I2=1 I3=NLEG I4=IABS(MODE) IF (I4.EQ.1) THEN I3=NLEG1 ELSEIF (I4.EQ.2) THEN I2=NLEG1 I3=NLEG12 ELSEIF (I4.EQ.3) THEN I2=NLEG12 ENDIF C Loop along all gridlegs, recording points at isosurfaces: DO 29, I1=I2,I3 CALL GIGLEG(I1,IP1,IP2) VAL1=RAM(IP1) VAL2=RAM(IP2) IF (((VAL1.LE.VALUE).AND.(VAL2.GT.VALUE)).OR. * ((VAL1.GT.VALUE).AND.(VAL2.LE.VALUE))) THEN C The points of the gridleg are at different sides C of the isosurface - searching for two neighbouring C gridpoints of the original grid, which are also C at the different sides of the isosurface: C JCOOR - index of the coordinate axis JCOOR=3 IF (I1.LE.NLEG12) JCOOR=2 IF (I1.LE.NLEG1) JCOOR=1 C J2 - number of gridsteps of the original grid C between IP1 and IP2 10 CONTINUE CALL GICP(IP1,XX) CALL GICP(IP2,YY) J2=IABS(NINT((XX(JCOOR)-YY(JCOOR))/DD(JCOOR))) J2=J2/2 IF (J2.NE.0) THEN IP3=IP1+J2*DIP(JCOOR) VAL3=RAM(IP3) IF (((VAL1.LE.VALUE).AND.(VAL3.LE.VALUE)).OR. * ((VAL1.GT.VALUE).AND.(VAL3.GT.VALUE))) THEN IP1=IP3 VAL1=VAL3 ELSE IP2=IP3 VAL2=VAL3 ENDIF GOTO 10 ENDIF C End of the search for neighbouring gridpoints. C Computing the coordinates ZZ of the point at the isosurface: ZZ(1)=XX(1) ZZ(2)=XX(2) ZZ(3)=XX(3) IF ((VAL2-VAL1).NE.0.) ZZ(JCOOR)= * XX(JCOOR)+(VALUE-VAL1)*(YY(JCOOR)-XX(JCOOR))/(VAL2-VAL1) C Writing the point: C Name of the point: NPTS=NPTS+1 DO 204, I2=0,LENG-1 TEXTP(LENG-I2:LENG-I2)= * CHAR(ICHAR('0')+MOD(NPTS,10**(I2+1))/10**I2) 204 CONTINUE C Setting the output format: FORMAT='(3A,03(F00.0,1X),A)' OUTMIN=0. OUTMAX=0. DO 214, I2=1,3 IF (OUTMIN.GT.ZZ(I2)) OUTMIN=ZZ(I2) IF (OUTMAX.LT.ZZ(I2)) OUTMAX=ZZ(I2) 214 CONTINUE CALL FORM1(OUTMIN,OUTMAX,FORMAT(8:15)) FORMAT(14:17)= '1X),' C Writing: IF ((MODE.LT.0).AND.(VAL1.GT.VAL2).AND.(FPTS1.NE.' ')) THEN WRITE(LU1,FORMAT) '''',TEXTP(1:LENG),''' ',ZZ,'/' ELSEIF (FPTS.NE.' ') THEN WRITE(LU,FORMAT) '''',TEXTP(1:LENG),''' ',ZZ,'/' ENDIF ENDIF 29 CONTINUE IF (FPTS.NE.' ') THEN WRITE(LU,'(A)') '/' CLOSE(LU) ENDIF IF ((MODE.LT.0).AND.(FPTS1.NE.' ')) THEN WRITE(LU1,'(A)') '/' CLOSE(LU1) ENDIF C WRITE(*,'(A)') *'+GRDISO: Done. ' STOP END C C======================================================================= C SUBROUTINE GIGLEG(ILEG,IPOIN1,IPOIN2) C C----------------------------------------------------------------------- C INTEGER ILEG,IPOIN1,IPOIN2 C Input: C ILEG ... Index of gridleg. C Output: C IPOIN1,IPOIN2 ... Indices of gridpoints forming the gridleg. C....................................................................... C Common block /RAMC/ to store the information about points on C structural interfaces and common block /MSC/ to store auxiliary C quantities: INCLUDE 'grdiso.inc' C grdiso.inc. C ........................... C Auxiliary storage locations: INTEGER I31,I21,I1,I2,I3,J1,J2,J3,J21,J31,JCOOR,ILEG1 C....................................................................... C ILEG1=ILEG-1 IF (ILEG.LE.NLEG1) THEN I31=ILEG1 / NN11N2 I21=(ILEG1 - I31*NN11N2) / NN11 I1=ILEG - I31*NN11N2 - I21*NN11 JCOOR=1 ELSEIF (ILEG.LE.NLEG12) THEN I31=(ILEG1-NLEG1) / NN1N21 I21=((ILEG1-NLEG1) - I31*NN1N21) / N1N I1=(ILEG-NLEG1) - I31*NN1N21 - I21*N1N JCOOR=2 ELSEIF (ILEG.LE.NLEG) THEN I31=(ILEG1-NLEG12) / NN1N2 I21=((ILEG1-NLEG12) - I31*NN1N2) / N1N I1=(ILEG-NLEG12) - I31*NN1N2 - I21*N1N JCOOR=3 ELSE C GRDISO-04 CALL ERROR ('GRDISO-04: Wrong ILEG.') C This error should not appear. Contact the author. ENDIF I2=I21+1 I3=I31+1 J1=NO1+(I1-1)*ND1 J2=NO2+(I2-1)*ND2 J3=NO3+(I3-1)*ND3 J21=J2-1 J31=J3-1 IPOIN1=J31*N1N2+J21*N1+J1 IPOIN2=IPOIN1+DIP(JCOOR) RETURN END C C C======================================================================= C SUBROUTINE GICP(IPOIN,XX) C C----------------------------------------------------------------------- C INTEGER IPOIN REAL XX(3) C IPOIN ... Index of a point. C XX ... Coordinates of the point. C....................................................................... C Common block /RAMC/ to store the information about points on C structural interfaces and common block /MSC/ to store auxiliary C quantities: INCLUDE 'grdiso.inc' C grdiso.inc. C ........................... C Auxiliary storage locations: INTEGER I31,I21,I11,IPOIN1 C....................................................................... C IPOIN1=IPOIN-1 I31= IPOIN1 / N1N2 I21=(IPOIN1 - I31*N1N2) / N1 I11= IPOIN - I31*N1N2 - I21*N1 - 1 XX(1)=O1+FLOAT(I11)*DD(1) XX(2)=O2+FLOAT(I21)*DD(2) XX(3)=O3+FLOAT(I31)*DD(3) RETURN 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