C Program 'grid' to generate the velocities in a rectangular grid C required for the full wave finite differences, the shortest path C calculation of seismic rays, eikonal equation 'finite differences', C or raster imaging of the model. Also an oblique vertical 2-D section C across the 3-D model may be gridded. C C The rectangular grid is specified in Cartesian coordinates, and then C transformed to the model coordinates. For the Cartesian coordinates C connected with a particular kind of curvilinear model coordinates see C subroutine CARTES of the file 'metric.for'. C The rectangular grid could also be specified in respect to the model C coordinates, and limited by coordinate planes specified in the model C coordinates. This option may be enabled by changing the first C executive statement ' LCART=.TRUE.' to ' LCART=.FALSE.'. C C Name of the main input data file read from the * external unit: C (1) 'GRID' C 'GRID'..Name of the main input data file described below. C Default: 'GRID'='grid.dat'. C C Main input data file 'GRID': C (1) 'MODEL','NET','IND','VEL','ICB','VEL1','VEL2','VEL3','VEL11', C 'VEL12','VEL22','VEL13','VEL23','VEL33',/ C 'MODEL'... Input data file describing the model, it is described C in the subroutine file 'model.for'. C 'NET'...Name of the input data file, described below, describing C the dimensions and density of a rectangular grid of C points. C 'IND'...If not blank, the name of the index file. C This option enables to specify other than rectangular C region covered by a rectangular grid: C The rectangular volume bounded by coordinate limits C X1MIN,X1MAX, X2MIN,X2MAX, AND X3MIN,X3MAX is divided into C N1*N2*N3 big bricks. Each element (index) of the index C file corresponds to one big brick. If it equals zero, C the big brick does not belong to the region in which the C velocity is discretized. C 'VEL'...Name of the output formatted file, C containing the velocities at the given gridpoints. C Velocity=0 is inserted in a free space. C If blank, the file is not created. C 'ICB'...Name of the output formatted file, C containing the indices of complex geological blocks at the C given gridpoints. C If blank, the file is not created. C 'VEL1','VEL2','VEL3','VEL11','VEL12','VEL22','VEL13','VEL23', C 'VEL33'... Names of the output formatted files containing C individual first or second partial velocity derivatives C at the given gridpoints. C If the filename is blank, the corresponding file is not C created. C Default: 'MODEL'='model.dat', 'NET'='net.dat', 'IND'=' ', C 'VEL'='vel.out', 'ICB'='icb.out', 'VEL1'=' ', 'VEL2'=' ', C 'VEL3'=' ', 'VEL11'=' ', 'VEL12'=' ', 'VEL22'=' ', C 'VEL13'=' ', 'VEL23'=' ', 'VEL33'=' '. C (2) (IPS(I),I=1,NPS),/ C IPS... List of indices of complex blocks terminated by a slash. C If a complex block is not listed, P-wave velocity is C evaluated at grid points. C IPS(I).LT.0: S-wave velocity is evaluated at grid points. C IPS(I).GT.0: free space (velocity=0) is assumed at grid C points. This option is likely be used to avoid C refraction in deep layers when computing reflected rays C by means of the shortest path network algorithm. C Default: P-wave velocity in all material complex blocks. C C Input file 'MODEL' is described within the FORTRAN77 source code file C 'model.for'. C C Input file 'NET': C The data are read in by the list directed input (free format). In C the list of input data below, each numbered paragraph indicates C the beginning of a new input operation (new READ statement). If C the first letter of the symbolic name of the input variable is C I-N, the corresponding value in input data must be of the type C INTEGER. Otherwise, the input parameter is of the type REAL. C (1) One of the following possibilities (1.1), (1.2), (1.3), or (1.4): C (1.1) N1,N2,N3,/ C a basic option: C N1,N2,N3... Numbers of gridpoints of the rectangular grid in the C directions of axes X1,X2,X3. C The rectangular volume bounded by coordinate limits C X1MIN,X1MAX, X2MIN,X2MAX, and X3MIN,X3MAX is divided into C N1*N2*N3 big bricks. C Gridpoints are defined as centres of the big bricks. C Default: N3=1. C (1.2) N1,N2,N3,L1,L2,L3,/ C A special option: C The rectangular volume bounded by coordinate limits C X1MIN,X1MAX, X2MIN,X2MAX, and X3MIN,X3MAX is divided into C N1*N2*N3 big bricks. C Then (if the numbers L1,L2,L3 are specified in addition to C N1,N2,N3) each big brick is subdivided into L1*L2*L3 small C bricks. C The output velocities are computed in the centres of small C bricks. C Outer loop is over big bricks, the discrete velocities C within each big brick being output consecutively. C N1,N2,N3... Numbers of big bricks in the X1,X2,X3-axis directions. C Must be positive. C N3 need not be specified for a 2-D model. C L1,L2,L3... Numbers of small bricks in one big brick in the C X1,X2,X3-axis directions. If specified, must be positive. C (If the numbers L1,L2,L3 are not specified, each big brick C contains just one small brick, as large as big one.) C Default: N3=1, L1=1, L2=1, L3=1. C (1.3) 0,N2,N3,/ C A special option: 2-D oblique vertical section in 3-D: C The rectangular vertical section bounded by the vertical C lines (X1,X2)=(X1MIN,X2MIN) and (X1,X2)=(X1MAX,X2MAX), C from X3=X3MIN to X3=X3MAX is divided into N2*N3 cells. C N2,N3.. Numbers of cells in the horizontal and vertical C directions. Must be positive. C Default: N3=1. C (1.4) 0,N2,N3,L1,L2,L3,/ C A special option: 2-D oblique vertical section in 3-D: C The rectangular vertical section bounded by the vertical C lines (X1,X2)=(X1MIN,X2MIN) and (X1,X2)=(X1MAX,X2MAX), C from X3=X3MIN to X3=X3MAX is divided into N2*N3 big cells, C each big cell being divided into L2*L3 small cells. C N2,N3.. Numbers of big cells in the horizontal and vertical C directions. Must be positive. C L1... Must equal 0 or 1. L1=0 is understood as L1=1. C L2,L3.. Numbers of small cells in one big cell in the horizontal C and vertical directions. If specified, must be positive. C Default: N3=1, L1=1, L2=1, L3=1. C (2) X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,/ C Boundaries of the volume covered by the rectangular homogeneous C grid. The X1 coordinates of the gridpoints are C X1MIN+0.5*DX1, X1MIN+1.5*DX1, ... ,X1MAX-0.5*DX1 , C where C DX1=(X1MAX-X1MIN)/N1, C or C DX1=(X1MAX-X1MIN)/N1*L1, C if L1 is specified, similarly for X2 and X3. C Default: X1MIN=0, X1MAX=0, X2MIN=0, X2MAX=0, X3MIN=0, X3MAX=0. C C Input file 'IND': C This option enables to specify other than rectangular C region covered by a rectangular grid: C The rectangular volume bounded by coordinate limits C X1MIN,X1MAX, X2MIN,X2MAX, and X3MIN,X3MAX is divided into C N1*N2*N3 big bricks. Each element (index) of the index C file corresponds to one big brick. If it equals zero, C the big brick does not belong to the region in which the C velocity is discretized. C Attention: The nonzero indices must be formed by the sequence C 1,2,3,... of positive integers. C (1) (IND(I),I=1,N1*N2*N3) C IND(I)..Zero if the I-th big brick does not belong to the region C in which the velocity is discretized. C Otherwise the index the big brick. The gridpoints within C the big brick are indexed by integers from C L1*L2*L3*(IND(I)-1)+1 to L1*L2*L3*IND(I). C Default: IND(I)=I. C C Output files 'VEL','VEL1','VEL2','VEL3','VEL11','VEL12','VEL22', C 'VEL13','VEL23': C (1) (V(I),I=1,L1234), where L1234 is the number of gridponts. C L1234=L1*L2*L3*L4. If the file 'IND' is not specified, C L4=N1*N2*N3 by the default. C V(I)... Velocity or its partial derivative at the I-th gridpoint. C Free space is indicated by a zero velocity or derivative C V(I)=0. C C Output file 'ICB': C (1) (ICB(I),I=1,L1234), where L1234 is the number of gridponts. C ICB(I)..Index of (geological) block in which the I-th gridpoint C is situated. C C Subroutines referenced: EXTERNAL KOOR,MODEL1,BLOCK,VELOC,PARM2,WARRAY INTEGER KOOR C KOOR... File 'metric.for'. C MODEL1,BLOCK,VELOC... File 'model.for'. C PARM2...File 'parm.for'. C WARRAY..File 'forms.for'. C Note that the above subroutines reference many other external C procedures from various source code files of the 'MODEL' subroutine C package. These indirectly referenced procedures are not named here, C but are listed in the particular subroutine source code files. C C Date: 1996, September 30 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Storage locations: C C Input data: CHARACTER*80 FGRID,FMODEL,FNET,FIND,FVEL,FICB,FVELD(9) INTEGER LU1,LU2,LUD(9),MPS,MIND PARAMETER (LU1=1,LU2=2,MPS=100,MIND=500000) INTEGER N1,N2,N3,L1,L2,L3,IPS(MPS),IND(MIND) REAL X1MIN,X2MIN,X3MIN,X1MAX,X2MAX,X3MAX C C LU1,LU2,LUD... Logical unit numbers. C C Others: LOGICAL LCART,LIND,LVEL0,LICB,LVELD,LVEL(9),LOBLIQ INTEGER MVEL PARAMETER (MVEL=10000) INTEGER N123,L1234,I1234,IN1,IN2,IN3,IL1,IL2,IL3,IBRICK,INDOLD INTEGER ICB(MVEL),NVEL,ISB1,ISRF2,ISB2,ICB2,I,II REAL VOUT(MVEL),VOUD(MVEL,9),DX1,DX2,DX3 REAL COOR(3),UP(10),US(10),VD(10),AUX0,AUX1,AUX2,AUX3,AUX4 REAL G(12),GAMMA(18),PDER(9),AUX11,AUX12,AUX22,AUX13,AUX23,AUX33 C C LCART.. True if the model specified in curvilinear coordinates is C gridded in Cartesian coordinates. C LIND... Indication of indexed grid to specify irregular subvolume C of the rectangular volume covered by the grid. C LVEL0,LICB,LVELD,LVEL... Indication of opening and generating C output files. C LOBLIQ..Indication of a special option enabling to grid an oblique C vertical profile. C ICB... Indices of complex blocks. C ISB1... Index of the simple block, see subroutine block. C ISRF2...Index of a surface, see subroutine block. C ISB2... Index of the simple block, see subroutine block. C ICB2... Index of the complex block, see subroutine block. C I... Index of a gridpoint, or loop variable. C DX1,DX2,DX3... Grid intervals. C VEL... Velocity. C COOR... Cordinates of a gridpoint. C UP,US,VD,AUX0,AUX1,AUX2,AUX3,AUX4... Auxiliary storage locations C for local model parameters or temporary variables. C G,GAMMA,PDER,AUX11,AUX12,AUX22,AUX13,AUX23,AUX33... Auxiliary C storage locations used in coordinate transformations. C DATA LUD/3,4,5,6,7,8,9,10,11/ C C....................................................................... C LCART=.TRUE. C C....................................................................... C C Opening files and reading input data: C C Name of main input data: WRITE(*,'(A)') ' Enter the name of the main input data file: ' FGRID='grid.dat' READ(*,*) FGRID WRITE(*,'(A)') '+ ' C C Reading main input data: OPEN(LU1,FILE=FGRID,STATUS='OLD') FMODEL='model.dat' FNET='net.dat' FIND=' ' FVEL='vel.out' FICB='icb.out' DO 10 I=1,9 FVELD(I)=' ' 10 CONTINUE READ(LU1,*) FMODEL,FNET,FIND,FVEL,FICB,(FVELD(I),I=1,9) DO 11 I=1,MPS IPS(I)=0 11 CONTINUE READ(LU1,*) IPS DO 15 I=1,MPS IF(IPS(I).NE.0) THEN IF(IABS(IPS(I)).LT.I) THEN IPS(IABS(IPS(I)))=IPS(I) ELSE IF(IABS(IPS(I)).GT.I) THEN DO 12 II=IABS(IPS(I)),MPS IF(IPS(II).EQ.0) THEN IPS(II)=IPS(I) IPS(I)=0 GO TO 13 END IF 12 CONTINUE PAUSE 'ERROR: TOO SMALL ARRAY IPS(MPS)' STOP 13 CONTINUE END IF END IF 15 CONTINUE DO 16 I=1,MPS IF(IPS(I).EQ.0) THEN IPS(I)=I ELSE IF(IPS(I).GT.0) THEN IPS(I)=0 END IF 16 CONTINUE CLOSE(LU1) C C Data for model: OPEN(LU1,FILE=FMODEL,STATUS='OLD') CALL MODEL1(LU1) CLOSE(LU1) IF(KOOR().EQ.0) THEN C No transformation between Cartesian and model coordinates LCART=.FALSE. END IF C C C Data for grid: OPEN(LU1,FILE=FNET,STATUS='OLD') N3=1 L1=1 L2=1 L3=1 READ(LU1,*) N1,N2,N3,L1,L2,L3 IF(N1.EQ.0) THEN C Vertical oblique profile IF(L1.NE.0.AND.L1.NE.1) THEN PAUSE 'ERROR: INCORRECT L1 FOR AN OBLIQUE PROFILE' STOP END IF LOBLIQ=.TRUE. N1=1 L1=1 ELSE LOBLIQ=.FALSE. END IF IF(N1.LT.1.OR.N2.LT.1.OR.N3.LT.1.OR.L1.LT.1.OR.L2.LT.1.OR.L3.LT.1) * THEN PAUSE 'ERROR: NON-POSITIVE NUMBER OF GRIDPOINTS' STOP END IF X1MIN=0. X1MAX=0. X2MIN=0. X2MAX=0. X3MIN=0. X3MAX=0. READ(LU1,*) X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX CLOSE(LU1) C C Reading the index array: IF(FIND.EQ.' ') THEN LIND=.FALSE. IND(1)=1 L1234=L1*L2*L3*N1*N2*N3 ELSE LIND=.TRUE. N123=N1*N2*N3 IF(N123.GT.MIND) THEN PAUSE 'ERROR: TOO MANY BIG BRICKS' STOP END IF DO 31 IBRICK=1,N123 IND(IBRICK)=IBRICK 31 CONTINUE OPEN(LU1,FILE=FIND,STATUS='OLD') READ(LU1,*) (IND(IBRICK),IBRICK=1,N123) CLOSE(LU1) L1234=0 DO 32 IBRICK=1,N123 IF(IND(IBRICK).GT.0) THEN L1234=L1234+1 END IF 32 CONTINUE L1234=L1*L2*L3*L1234 END IF C C Output file with velocities at gridpoints: IF(FVEL.EQ.' ') THEN LVEL0=.FALSE. ELSE LVEL0=.TRUE. OPEN(LU1,FILE=FVEL) END IF C C Output file with indices of complex blocks: IF(FICB.EQ.' ') THEN LICB=.FALSE. ELSE LICB=.TRUE. OPEN(LU2,FILE=FICB) END IF C C Output file with velocity derivatives at gridpoints: LVELD=.FALSE. DO 33 I=1,9 IF(FVELD(I).EQ.' ') THEN LVEL(I)=.FALSE. ELSE LVELD =.TRUE. LVEL(I)=.TRUE. OPEN(LUD(I),FILE=FVELD(I)) END IF 33 CONTINUE C C....................................................................... C C Loops over gridpoints: C IF(LOBLIQ) THEN DX1=(X1MAX-X1MIN)/FLOAT(N2*L2) ELSE DX1=(X1MAX-X1MIN)/FLOAT(N1*L1) END IF DX2=(X2MAX-X2MIN)/FLOAT(N2*L2) DX3=(X3MAX-X3MIN)/FLOAT(N3*L3) ISB1=0 NVEL=0 IBRICK=0 INDOLD=0 I1234=0 WRITE(*,'(''+'',I18,'' gridpoints of'',I8)') I1234,L1234 C C Loop over big bricks: DO 83 IN3=0,N3-1 DO 82 IN2=0,N2-1 DO 81 IN1=0,N1-1 C C Check for the computational volume: IF(LIND) THEN IBRICK=IBRICK+1 IF(IND(IBRICK).EQ.0) THEN GO TO 80 END IF IF(IND(IBRICK).NE.INDOLD+1) THEN PAUSE 'ERROR: INDICES NOT CONSECUTIVE' STOP END IF INDOLD=IND(IBRICK) END IF C C Loop over small bricks: DO 73 IL3=1,L3 COOR(3)=X3MIN+DX3*(FLOAT(IN3*L3+IL3)-0.5) DO 72 IL2=1,L2 COOR(2)=X2MIN+DX2*(FLOAT(IN2*L2+IL2)-0.5) DO 71 IL1=1,L1 IF(LOBLIQ) THEN COOR(1)=X1MIN+DX1*(FLOAT(IN2*L2+IL2)-0.5) ELSE COOR(1)=X1MIN+DX1*(FLOAT(IN1*L1+IL1)-0.5) END IF C C Transformation from Cartesian to model coordinates: IF(LCART) THEN G(1)=COOR(1) G(2)=COOR(2) G(3)=COOR(3) CALL CARTES(COOR,.FALSE.,G,PDER) END IF C C Velocity evaluation: CALL BLOCK(.TRUE.,COOR,0,ISB1,ISRF2,ISB2,ICB2,VD) ISB1=ISB2 IF(ICB2.EQ.0) THEN C free space: DO 41 I=1,10 VD(I)=0. 41 CONTINUE ELSE IF(IPS(ICB2).EQ.0) THEN C Deemed free space: ICB2=0 DO 42 I=1,10 VD(I)=0. 42 CONTINUE ELSE CALL PARM2(IABS(ICB2),COOR,UP,US,AUX0,AUX1,AUX2) CALL VELOC(IPS(ICB2),UP,US, * AUX1,AUX2,AUX3,AUX4,VD,AUX0) END IF C C Writing output files: IF(NVEL.EQ.MVEL) THEN WRITE(*,'(''+Writing: '',I8)') I1234 IF(LVEL0) THEN CALL WARRAY(LU1,' ','FORMATTED', * .FALSE.,0.,.FALSE.,0.,MVEL,VOUT) C For velocities up to 9.99999, the above statement C might be replaced, for instance, by: C WRITE(LU1,'(10F8.5)') VOUT END IF IF(LICB) THEN WRITE(LU2,'(20(1X,I2))') ICB END IF DO 51 I=1,9 IF(LVEL(I)) THEN CALL WARRAY(LUD(I),' ','FORMATTED', * .FALSE.,0.,.FALSE.,0.,MVEL,VOUD(1,I)) END IF 51 CONTINUE NVEL=0 WRITE(*,'(''+ '')') END IF NVEL=NVEL+1 VOUT(NVEL)=VD(1) ICB (NVEL)=ICB2 IF(LVELD) THEN IF(LCART) THEN C Transformation from model to Cartesian coordinates C covariant derivatives CALL METRIC(COOR,AUX1,G,GAMMA) AUX1=VD(2) AUX2=VD(3) AUX3=VD(4) AUX11=VD( 5)-GAMMA(1)*AUX1-GAMMA( 7)*AUX2 * -GAMMA(13)*AUX3 AUX12=VD( 6)-GAMMA(2)*AUX1-GAMMA( 8)*AUX2 * -GAMMA(14)*AUX3 AUX22=VD( 7)-GAMMA(3)*AUX1-GAMMA( 9)*AUX2 * -GAMMA(15)*AUX3 AUX13=VD( 8)-GAMMA(4)*AUX1-GAMMA(10)*AUX2 * -GAMMA(16)*AUX3 AUX23=VD( 9)-GAMMA(5)*AUX1-GAMMA(11)*AUX2 * -GAMMA(17)*AUX3 AUX33=VD(10)-GAMMA(6)*AUX1-GAMMA(12)*AUX2 * -GAMMA(18)*AUX3 C Transformation of derivatives VD(2)= AUX1*PDER(1)+ AUX2*PDER(2)+ AUX3*PDER(3) VD(3)= AUX1*PDER(4)+ AUX2*PDER(5)+ AUX3*PDER(6) VD(4)= AUX1*PDER(7)+ AUX2*PDER(8)+ AUX3*PDER(9) AUX1 =AUX11*PDER(1)+AUX12*PDER(2)+AUX13*PDER(3) AUX2 =AUX12*PDER(1)+AUX22*PDER(2)+AUX23*PDER(3) AUX3 =AUX13*PDER(1)+AUX23*PDER(2)+AUX33*PDER(3) VD(5)= AUX1*PDER(1)+ AUX2*PDER(2)+ AUX3*PDER(3) AUX1 =AUX11*PDER(4)+AUX12*PDER(5)+AUX13*PDER(6) AUX2 =AUX12*PDER(4)+AUX22*PDER(5)+AUX23*PDER(6) AUX3 =AUX13*PDER(4)+AUX23*PDER(5)+AUX33*PDER(6) VD(6)= AUX1*PDER(1)+ AUX2*PDER(2)+ AUX3*PDER(3) VD(7)= AUX1*PDER(4)+ AUX2*PDER(5)+ AUX3*PDER(6) AUX1 =AUX11*PDER(7)+AUX12*PDER(8)+AUX13*PDER(9) AUX2 =AUX12*PDER(7)+AUX22*PDER(8)+AUX23*PDER(9) AUX3 =AUX13*PDER(7)+AUX23*PDER(8)+AUX33*PDER(9) VD(8)= AUX1*PDER(1)+ AUX2*PDER(2)+ AUX3*PDER(3) VD(9)= AUX1*PDER(4)+ AUX2*PDER(5)+ AUX3*PDER(6) VD(10)=AUX1*PDER(7)+ AUX2*PDER(8)+ AUX3*PDER(9) END IF DO 61 I=1,9 IF(LVEL(I)) THEN VOUD(NVEL,I)=VD(I+1) END IF 61 CONTINUE END IF C C Screen output: I1234=I1234+1 IF(MOD(I1234,1000).EQ.0) THEN WRITE(*,'(''+'',I18)') I1234 END IF C 71 CONTINUE 72 CONTINUE 73 CONTINUE C 80 CONTINUE 81 CONTINUE 82 CONTINUE 83 CONTINUE C C Rest of output files: IF(NVEL.GT.0) THEN WRITE(*,'(''+Writing: '',I8)') I1234 IF(LVEL0) THEN CALL WARRAY(LU1,' ','FORMATTED', * .FALSE.,0.,.FALSE.,0.,NVEL,VOUT) C For velocities up to 9.99999, the above statement C might be replaced, for instance, by: C WRITE(LU1,'(10F8.5)') (VOUT(I),I=1,NVEL) END IF IF(LICB) THEN WRITE(LU2,'(20(1X,I2))') (ICB(I),I=1,NVEL) END IF DO 91 I=1,9 IF(LVEL(I)) THEN CALL WARRAY(LUD(I),' ','FORMATTED', * .FALSE.,0.,.FALSE.,0.,NVEL,VOUD(1,I)) END IF 91 CONTINUE NVEL=0 END IF C C Screen output: WRITE(*,'(''+'',I18)') I1234 STOP END C C======================================================================= C INCLUDE 'model.for' INCLUDE 'metric.for' INCLUDE 'srfc.for' INCLUDE 'parm.for' INCLUDE 'val.for' INCLUDE 'fit.for' INCLUDE 'forms.for' C C======================================================================= C