addrec.cal 0100666 0000766 0000766 00000000011 06724651754 012344 0 ustar klimes klimes @3=@1+$2 fd2d.for 0100666 0000766 0000766 00000254135 06701034720 011771 0 ustar klimes klimes C
C Program FD2D for 2-D P-SV elastic second-order finite differences C C Version: 5.30 C Date: 1999, April 2 C C Program written by Jiri Zahradnik C Department of Geophysics, Charles University Prague C V Holesovickach 2, 180 00 Praha 8, Czech Republic C Tel.: +420-2-21912546, Fax: +420-2-21912555 C E-mail: jz@karel.troja.mff.cuni.cz C C....................................................................... C C Contents: C Description of data files C General description C List of external procedures C C For the description of effective elastic parameters refer to file C modfd.for. C C....................................................................... C C Description of data files: C C Main input data read from the * input device: C The data are read in by the list directed input (free format), by C a single READ statement. The data consist of character strings C which have to be enclosed in apostrophes. Note that the C interactive * external unit may be piped or redirected to a file. C The following items are read: C (1) 'FDSEP' / C 'FDSEP'... String in apostrophes containing the name of the main C input data file with FD parameters in the SEP format. C The file specifies grid dimensions, FD time step, the form C and names of the files with effective elastic parameters C and snapshots, the type and position of the source, C the names of the files with receiver positions, snapshot C levels, the names of the files with synthetic seismograms C and snapshots, etc. C /... It is recommended to terminate these data by a slash. C Default: 'FDSEP'='fd.h'. C C Data file 'FDSEP' has the form of the SEP (Stanford Exploration C Project) 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 This program understands the following parameters common to nearly all C programs dealing with SEP-like gridded data formats: C N1=positive integer... Number of gridpoints of the basic grid C along the X1 axis. Note that 20 leftmost and 20 rightmost C gridlines are reserved for the dumpers and that an C explosive source requires 3 by 3 gridpoints. C Default: N1=1 C N2=positive integer... Number of gridpoints of the basic grid C along the X2 axis. In 2-D, leave the default value. C Default: N2=1 C N3=positive integer... Number of gridpoints of the basic grid C along the X3 axis. Note that the first (top) horizontal C gridline is assumed to be a free surface and that 20 C bottom gridlines are reserved for the dumpers. C Default: N3=1 C O1=real... X1 coordinate of the origin of the grid. C Default: O1=0. C O3=real... X3 coordinate of the origin of the grid. C Note that FD program 'fd2d.for' assumes the free Earth C surface at vertical coordinate X3=O3. C Default: O3=0. C D1=real... Grid spacing along the X1 axis. C The maximum relative error E in the wavefield is roughly C E = (2*pi*F)**3 * (D1/VSmin)**2 * Tmax / 24 , C where VSmin is the minimum S-wave velocity in the model, C Tmax=(NTFD-1)*DTFD is the maximum time and F is the mean C frequency of the elastic waves, often roughly equal to the C reference frequency SIGF of the source time function. C Equation C D1 = VSmin * SQRT( 3 * E / Tmax / (pi*F)**3 ) C may thus be used to estimate spatial steps D1 and D3. C Note that old-fashioned equation C D1 = VSmin / (10*Fmax) = VSmin / (30*SIGF) C corresponds to the maximum relative error C E = SIGF * Tmax * 1.15% C in the calculated wavefield. C Note that the grid intervals may also be negative. C Default: D1=1. C D3=real... Grid spacing along the X3 axis. C Absolute values of D1 and D3 must be equal. C Default: D3=1. C Number and size of FD time steps: C DTFD=real... Time step for finite differences. C Choose DTFD as to satisfy the stability condition C DTFD < D1 / (1.65*VPmax) , C close to the upper limit. C Here VPmax is the maximum P-wave velocity in the model. C NTFD=integer... Number of time levels. C Other numerical parameters: C LA1=integer... It is recommended not to specify this parameter. C The medium is assumed homogeneous (in the 'bedrock') for C the grid lines L > LA1-2. If there is no homogeneous C bedrock, choose the default of LA1=N3-1. C Default: LA1=N3-1 C Type and position of the source: C SRC='string'... Name of the input file with the position of the C point source. If SRC is blank, the default coordinates C (0,0,0) of the source are used. C Description of file SRC C Default: SRC=' '. C KSIG=integer... Type of the source time function: C KSIG=1: Gabor signal. C KSIG=3: Kuepper (Mueller) signal. C Default: KSIG=0 C SIGF=positive real... Reference frequency characterizing the C excitation wavelet. The Kuepper (Mueller) wavelet has C time duration 1/SIGF. For frequencies f > Fmax = 3*SIGF, C the spectral amplitude is negligibly small. C The excitation signal specifies the time dependence of the C density of the seismic force or the seismic moment along C a linear source. C To obtain the far-field displacement similar to that of C a point source, the excitation signal should specify the C time dependence of the halfth derivative of the point C seismic force or the point seismic moment. C Note that reference half-period TP=1/(2*SIGF) was used in C previous versions of this program instead of SIGF. C Default: SIGF=1. C SIGA=real... Reference amplitude of the excitation wavelet. C Maximum amplitude of the Kuepper wavelet for KSIG=3. C The amplitude describes the intensity of a linear seismic C force or moment. The 3-D amplitude corresponding to C a point source is asymptotically 1/SQRT(2*pi*SIGMA) C greater. Here SIGMA is the integral of the propagation C velocity along the raypath (in a homogeneous medium, C SIGMA=V*R, where V is the velocity and R the hypocentral C distance). Thus, if we wish to obtain the amplitude at C hypocentral distance R roughly corresponding to point C source of intensity SIGA=A, we may choose the line source C of intensity SIGA=A/SQRT(2*pi*SIGMA), where SIGMA C corresponds to the wave type and hypocentral distance. C Default: SIGA=1. C Files with receiver positions and synthetic seismograms: C REC='string'... Name of the input file with the positions of the C receivers. If REC is blank, no receivers are considered. C Description of file REC C Default: REC=' '. C SS='string'... Name of the output GSE file with the synthetic C seismograms. The file is created only if the receivers C are specified. C Default: SS='fd2d.gse'. C OUTU='string',OUTW='string'... Names of the output files with the C synthetic seismograms in the previously used format. Just C for the backward compatibility with program 'trans.for'. C It is recommended not to specify them. C Files given by OUTU and OUTW, formerly named 'outu.dat' C and 'outw.dat', correspond to the horizontal (U) and C vertical (W) component, respectively. These output files C are organized as 'time levels', not seismograms. It C means that the computed displacement values for all C receivers at a single time level are recorded, and this C process is repeated subsequently for the all computed C time levels. The output files can be processed by C program 'trans.for' which provides the output suitable C for plotting purposes (first column is time, second C column is the seismogram for the first receiver, etc.). C Attention: Only each 4th time level is on output. C The output time series have therefore the time step of C 4*DTFD. Files OUTU and OUTW may be tranformed by C 'trans.for' into 'a.dat' and 'b.dat', respectively. C These files may serve for plotting by GRAPHER. C Defaults: OUTU=' ', OUTW=' ' C Output files with snapshots: C O4=real... Time of the first snapshot. C Default: O4=0. C D4=real... Time interval to store the snapshots. C Default: D4=DTFD C N4=integer... Number of the snapshots. C Default: N4=NTFD C SNAP1='string'... Name of the output file with the N1*N3*N4 grid C values of the horizontal component of the displacement. C The file is not created if SNAP1=' '. C Default: SNAP1=' ' C SNAP3='string'... Name of the output file with the N1*N3*N4 grid C values of the vertical component of the displacement. C The file is not created if SNAP3=' '. C Default: SNAP3=' ' C Form of files with effective elastic parameters and snapshots: C FORM='string'... Form of the input and output files: C FORM='formatted': Formatted ASCII files, C FORM='unformatted': Unformatted files. C Input files with effective elastic parameters: C A11='string', B11='string', C11='string'... Strings in apostrophes C containing the names of the input files with (N1-1)*N2*N3 C values of the elastic parameters averaged in the direction C of the X1 axis, along the gridlines. C Defaults: A11=' ', B11=' ', C11=' '. C A33='string', B33='string', C33='string'... Strings in apostrophes C containing the names of the input files with N1*N2*(N3-1) C values of the elastic parameters averaged in the direction C of the X3 axis, along the gridlines. C Defaults: A33=' ', B33=' ', C33=' '. C A13='string', B13='string', C13='string'... Strings in apostrophes C containing the names of the input files with C (N1-1)*N2*(N3-1) values of the elastic parameters averaged C in the direction of the X1 axis, along the gridlines C shifted half a grid interval in the direction of the X3 C axis. C Defaults: A13=' ', B13=' ', C13=' '. C A31='string', B31='string', C31='string'... Strings in apostrophes C containing the names of the input files with C (N1-1)*N2*(N3-1) values of the elastic parameters averaged C in the direction of the X3 axis, along the gridlines C shifted half a grid interval in the direction of the X1 C axis. C Defaults: A31=' ', B31=' ', C31=' '. C Analogously A22, B22, C22, A12, B12, C12, A21, B21, C21, A23, B23, C C23, A32, B32 and C32. C DEN='string'... String in apostrophes containing the name of the C input file with N1*N2*N3 grid values of the density. C Default: DEN=' ' C Example of data set FDSEP C C C Input file SRC containing the coordinates of the point source: C (1) Several strings terminated by / (a slash). C The simplest way is to submit just the /. C (2) 'NAME',X1SRC,X2SRC,X3SRC,/ C 'NAME'... String reserved for the name of the initial point. C No meaning here, anything in apostrophes may be submitted. C X1SRC,X2SRC,X3SRC... Coordinates of a single initial point. C Default: X1SRC=0., X2SRC=0., X3SRC=0. C (3) / (a slash) or end of file. C Example of data set SRC C C C Input file REC containing receiver coordinates: C (1) Several strings terminated by / (a slash). C The simplest way is to submit just the /. C (2) Several times (2.1): C (2.1) 'NAMER(IR)',X1REC(IR),X2REC(IR),X3REC(IR),/ C 'NAMER(IR)'... String reserved for the name of the receiver. C No meaning here, anything in apostrophes may be submitted. C X1REC(IR),X2REC(IR),X3REC(IR)... Coordinates of the receiver. C Default: X1REC(IR)=0, X2REC(IR)=0, X3REC(IR)=0. C (3) / (a slash) or end of file. C Example of data file REC C C....................................................................... C C General description: C C Calculation of P-SV waves in 2-D geological structures by the C finite-difference method of the second order accuracy in both in space C and time. No absorption included. C C The explicit PS2 scheme described in Zahradnik and Priolo (1995) C is used. The uniform square grid: spatial step D1=D3. C C The computational region is a rectangle. C The top side is the flat free Earth's surface. C The left, right and bottom sides are nonreflecting after C Emerman and Stephen. Simultaneously, artificial dumpers, 20*D1 wide, C are applied at these three sides. C C The excitation is by a line source, described by Aboudi's body force C or seismic moment corresponding to the explosive source. C The excitation by a plane wave incident from below (well tested only C for the vertical incidence), realized by the algorithm of Alterman and C Karal, will be included later. C C Reference: C Zahradnik, J. and E. Priolo (1995): Heterogeneous C formulations of elastodynamic equations and finite-difference C schemes. Geophys. J. Int., 120, 663-676. C C If necessary, ask the author for more papers on the related subjects C (absorption, filtration, convolution with earthquake C ground motions, applications). C C Remarks: C Program FD2D should not be used to calculate elastic wavefields C in fluids! Spurious waves, other than P waves, and instabilities C may occur in fluids. C C The absorbing boundary conditions and dumpers do not work C perfectly! The reflections from the boundaries may be found in C the seismograms by calculating the seismograms with two different C positions of the boundaries and comparing the results. C C Mueller's time function C Program modification is possible for another time C function, see TIMFUN at the end of this code. C C The line source is of the explosive (compressional) type. An easy C program modification is possible for another source type, C e.g., the vertical force, etc. C C The nonreflecting boundaries after Emerman and Stephen are C applied. Program modification is possible for changing the right C side to Reynolds. C C The physical model is between the columns K=3 and K=NX-2. C When the right side is symmetrical, it is located at K=NX-2. C The columns K=1, 2, NX, NX-1 serve for the nonreflecting C (or symmetrical) sides. C C The left, right and bottom sides = nonreflecting + dumpers. C Dumpers are 20 gridpoints wide. Therefore, smaller models than C 41(horiz)x21(vert) are not possible. C A program modification is possible for making the left or right C side symmetrical and/or to disable the dumpers. C C It is useful to have some space, at least 5 gridlines, between C the investigated structure and the nonreflecting sides. C C The program has not been optimized. The main intention was to C permit many different models of the medium with no change in the C code. Therefore, the largest part of the code form complicated C subroutines manipulating with the parameters of the medium. C The code has been originally developed as a 4th-order version. C Although the present version is 2nd-order only, many formal C complications from the planned 4th order still remain in the code. C Ignore numerous warnings 'formal argument ... never used'. C C....................................................................... C C This file consists of the following external procedures: C FD2D... Main program allocating the memory for the FD code. C FD2D C FD... Subroutine containing the main finite-difference code C (originally the main program). C FD C CDUMP...Function computes spatial variation of artificial dumpers. C CDUMP C AUXL1...Subroutine designed to transfer the displacements from C arrays UX,UY,etc., into auxiliary variables UOLD,U1, C U2,etc. (Selecting values from vicinity of a grid point C at line L=1.) C AUXL1 C AUXL2...Subroutine designed to transfer the displacements from C arrays UX,UY,etc., into auxiliary variables UOLD,U1, C U2,etc. (Selecting values from vicinity of a grid point C at line L=2.) C AUXL2 C AUXLA...Subroutine designed to transfer the displacements from C arrays UX,UY,etc., into auxiliary variables UOLD,U1, C U2,etc. (Selecting values from vicinity of a grid point C at line L>2.) C AUXLA C INDEX...Subroutine designed to compute 1-D index ICEN of a grid C point K,L, indices II(1)-II(25) of the grid-point C neighbours, indices isa,isb,...ish of the neighbouring C small legs (h/2), computing indices ila,ilb,...ilh of C the neighbouring large legs (h). C INDEX C PARL1...Subroutine designed to transfer the effective parameters C from arrays AHL,AHS,etc. into auxiliary variables. C (Selecting values from vicinity of a grid point at C line L=1.) C PARL1 C PARL2...Subroutine designed to transfer the effective parameters C from arrays AHL,AHS,etc. into auxiliary variables. C (Selecting values from vicinity of a grid point at C line L=2.) C PARL2 C PARL... Subroutine designed to transfer the effective parameters C from arrays AH,BH,etc. into auxiliary variables. C (Selecting values from vicinity of a grid point at C line L>2.) C PARL C PARHOM..Subroutine analogous to PARL but for a grid point with C a completely homogeneous vicinity. C PARHOM C DIFL1...Subroutine devoted to finite-difference scheme for a grid C point at the free surface at line L=1 C DIFL1 C DIFL2...Subroutine devoted to finite-difference scheme for a grid C point in a heterogeneous medium at line L=2. C DIFL2 C DIFL... Subroutine devoted to finite-difference scheme for a grid C point in a heterogeneous medium at line L>2. C DIFL C DIFHOM..Subroutine devoted to finite-difference scheme for a grid C point in a homogeneous medium. C DIFHOM C FLERIB..Subroutine devoted to left and right nonreflecting C boundaries. C FLERIB C FBOB... Subroutine devoted to bottom nonreflecting boundaries. C FBOB C FCORNB..Subroutine devoted to left and right bottom corners of C the nonreflecting boundaries. C FCORNB C SAVER...Subroutine saving 'very old' values (from time level M-1) C at boundary lines for subsequent computations at C nonreflecting boundaries. C SAVER C TIMFUN..Function devoted to the excitation signal. C TIMFUN C C======================================================================= C C C PROGRAM FD2D C C Main program allocating the memory for the FD code. C C Subroutines and external functions required: EXTERNAL FD,RSEP1,RSEP3I C FD... This file. C RSEP1,RSEP3I... File 'sep.for'. C Note that the above subroutines reference many other external C procedures from various subroutine files. These indirectly C referenced procedures are not named here, but are listed in the C particular subroutine files. C C----------------------------------------------------------------------- C C Common block /RAMC/ to allocate large array RAM: INCLUDE 'ram.inc' C ram.inc INTEGER IRAM(MRAM) EQUIVALENCE (IRAM,RAM) C C----------------------------------------------------------------------- C C Input and output data filenames and logical unit numbers: INTEGER LU1,LU2,LU3,LU4 PARAMETER (LU1=1,LU2=2,LU3=3,LU4=4) CHARACTER*80 FSEP C LU1,LU2,LU3,LU4... Logical unit numbers to be used to read and C write data files. C FSEP... The name of the input data file with FD parameters in SEP C format. C INTEGER MREC,N1,N3,N,NTFD,NR,NS,NGRID,I1,I2,I3,I4 PARAMETER (MREC=1024) CHARACTER*6 REC(MREC) C C....................................................................... C C Main input data: C Default: FSEP='fd.h' C Reading main input data: WRITE(*,'(A)') ' FD2D: Enter name of input SEP file [''fd.h'']: ' READ (*,*) FSEP WRITE(*,'(A)') '+FD2D: Working... ' C C Reading all the data from input SEP parameter file: CALL RSEP1(LU1,FSEP) C C Data specifying numbers of gridpoints and FD time steps: CALL RSEP3I('N1' ,N1 ,1) CALL RSEP3I('N3' ,N3 ,1) CALL RSEP3I('NTFD',NTFD,1) C C Allocating the memory in array RAM: NGRID=19*N1*N3+10*N1+12*N3 NS=2*NTFD+5 IF(NGRID+NS.GT.MRAM) THEN C FD2D-01 CALL ERROR('FD2D-01: Small dimension MRAM in ram.inc') C Dimension MRAM must be at least C 19*N1*N3+10*N1+12*N3+(2*NTFD+5)*MAX(1,NREC), C where N1, N3 and NTFD are parameters in the main input file and C NREC is the number of receivers specified in the file which name C is given by parameter REC. END IF N =N1*N3 NR=MIN0((MRAM-NGRID)/NS,MREC) NS=NTFD*NR I1= 1+19*N I2=I1+10*N1 I3=I2+12*N3 I4=I3+ 5*NR C CALL FD(LU1,LU2,LU3,LU4,N1,N3,N,NR,NS, * RAM( 1),RAM( N + 1),RAM( 2*N + 1),RAM( 3*N + 1), * RAM( 4*N + 1),RAM( 5*N + 1),RAM( 6*N + 1),RAM( 7*N + 1), * RAM( 8*N + 1),RAM( 9*N + 1),RAM(10*N + 1),RAM(11*N + 1), * RAM(12*N + 1),RAM(13*N + 1),RAM(14*N + 1),RAM(15*N + 1), * RAM(16*N + 1),RAM(17*N + 1),RAM(18*N + 1), * RAM( I1),RAM( N1+I1),RAM( 2*N1+I1),RAM( 3*N1+I1), * RAM( 4*N1+I1),RAM( 5*N1+I1),RAM( 6*N1+I1), * RAM( I2),RAM( N3+I2),RAM( 2*N3+I2),RAM( 3*N3+I2), * RAM( 4*N3+I2),RAM( 5*N3+I2),RAM( 6*N3+I2),RAM( 7*N3+I2), * RAM( 8*N3+I2),RAM( 9*N3+I2),RAM(10*N3+I2),RAM(11*N3+I2), * IRAM( I3),IRAM( NR+I3),IRAM(2*NR+I3), * RAM( 3*NR+I3),RAM( 4*NR+I3), * RAM( I4),RAM( NS+I4),REC) C WRITE (*,'(A)') '+FD2D: Done. ' STOP END C C C======================================================================= C C C SUBROUTINE FD(LU1,LU2,LU3,LU4,KMAX,LMAX,KLMAX,MAXREC,MAXSS, * UX,UY,WX,WY,AHS,BHS,CHS,AVS,BVS,CVS, * AHL,BHL,CHL,AVL,BVL,CVL,DEN,DUMP,WORK, * UNZ1,UNZ2,UNZ,WNZ1,WNZ2,WNZ,EX, * UNX1,UNX2,UNX,WNX1,WNX2,WNX,UK1,UK2,UK3,WK1,WK2,WK3, * KRE,LRE,INREC,ULEV,WLEV,USS,WSS,REC) INTEGER LU1,LU2,LU3,LU4,KMAX,LMAX,KLMAX,MAXREC,MAXSS REAL UX(KLMAX),UY(KLMAX),WX(KLMAX),WY(KLMAX) REAL AHS(KLMAX),BHS(KLMAX),CHS(KLMAX) REAL AVS(KLMAX),BVS(KLMAX),CVS(KLMAX) REAL AHL(KLMAX),BHL(KLMAX),CHL(KLMAX) REAL AVL(KLMAX),BVL(KLMAX),CVL(KLMAX) REAL DEN(KLMAX),DUMP(KLMAX),WORK(KLMAX) REAL UNZ1(KMAX),UNZ2(KMAX),UNZ(KMAX) REAL WNZ1(KMAX),WNZ2(KMAX),WNZ(KMAX) REAL EX(4,KMAX) REAL UNX1(LMAX),UNX2(LMAX),UNX(LMAX) REAL WNX1(LMAX),WNX2(LMAX),WNX(LMAX) REAL UK1(LMAX),UK2(LMAX),UK3(LMAX) REAL WK1(LMAX),WK2(LMAX),WK3(LMAX) INTEGER KRE(MAXREC),LRE(MAXREC),INREC(MAXREC) REAL ULEV(MAXREC),WLEV(MAXREC) REAL USS(MAXSS),WSS(MAXSS) CHARACTER*(*) REC(MAXREC) C C Subroutine (originally main program) reading the input data and C computing P-SV waves in 2-D geological structures by the C finite-difference method. C C Input: C LU1,LU2,LU3,LU4... Logical unit numbers to be used to read and C write data files. C KMAX,LMAX,KLMAX,MAXREC,MAXSS... Array dimensions. C AHL,BHL,CHL, C AVL,BVL,CVL, C AHS,BHS,CHS, C AVS,BVS,CVS, C DEN,QFC...arrays with effective parameters C WORK... temporary working array C REC... C C No output. C C Subroutines and external functions required: EXTERNAL INDEX EXTERNAL RSEP1,RSEP3I,RSEP3R EXTERNAL FDGRID,FDOUT1,FDREC,FDIN C INDEX...This file. C RSEP1,RSEP3I,RSEP3R...File 'sep.for'. C FDGRID,FDOUT1,FDREC,FDIN..File 'fdaux.for'. C Note that the above subroutines reference many other external C procedures from various subroutine files. These indirectly C referenced procedures are not named here, but are listed in the C particular subroutine files. C C----------------------------------------------------------------------- C C Storage locations: C C Input and output data filenames: CHARACTER*80 FILNA1,FILNA2,OUTSS,FOUTU,FOUTW CHARACTER*11 FORM C DIMENSION II(25),IJ(25) DIMENSION INU(6),INW(6),UP(6),WP(6) C COMMON /DSPLC/ UOLD,WOLD, 1 U1,U2,U3,U4,U5,U6,U7,U8,U9,U10,U11,U12,U13,U14, 2 U15,U16,U17,U18,U19,U20,U21,U22,U23,U24,U25, 3 W1,W2,W3,W4,W5,W6,W7,W8,W9,W10,W11,W12,W13,W14, 4 W15,W16,W17,W18,W19,W20,W21,W22,W23,W24,W25 COMMON /PPAR/ pahs,pdes,rbcs,rfgs, 1 shs,sas,ses,sds,qgs,qfs,qbs,qcs, 2 pahl,pdel,rbcl,rfgl, 3 shl,sal,sel,sdl,qgl,qfl,qbl,qcl, 4 phom,qhom,rhom,shom,dnst COMMON /INDI/ ICEN,II,IJ, 1 isa,isb,isc,isd,ise,isf,isg,ish, 2 ila,ilb,ilc,ild,ile,ilf,ilg,ilh COMMON /BASE/ MPAR,NSIZE,NX,NZ,NZ5,NZ6,CN,CNST COMMON /DMP/ NZONE C C....................................................................... C C Data specifying FD grid dimensions CALL FDGRID(KMAX,LMAX,KLMAX,NX,NZ,DX) C Data required to write synthetic seismograms CALL RSEP3R('O1',O1,0.0) CALL RSEP3R('O3',O3,0.0) CALL RSEP3R('D1',D1,1.0) CALL RSEP3R('D3',D3,1.0) C C Number and size of FD time steps: CALL RSEP3I('NTFD',NT,1) CALL RSEP3R('DTFD',DT,1.) C C Type and position of the source: C (arguments: Name of parameter in input data, Variable, Default): C The excitation (plane waves P=0, S=1; line source=2): C (for future extension, do not specify in input data, use ISOURC=2) CALL RSEP3I('ISOURC',ISOURC,2) C Grid indices of the source: C (for backward compatibility, do not specify in input data) CALL RSEP3I('KSOUR',KSOUR,1) CALL RSEP3I('LSOUR',LSOUR,1) C Semiduration TP of the excitation wavelet: CALL RSEP3R('SIGF',SIGF,1.) TP=0.5/SIGF C Medium must be homogeneous for L.GE.LA1-2. Choose LA1.GE.4. C If a plane wave is specified, wavefront t=0 is at L=LA1+3 and C L=LA1 to LA1+3 are prescribed. CALL RSEP3I('LA1',LA1,NZ-1) C C Reading receiver (and source) positions: CALL FDREC(LU1,MAXREC,NREC,REC,KRE,LRE,KSOUR,LSOUR) IF(NT*NREC.GT.MAXSS) THEN C FD2D-02 CALL ERROR('FD2D-02: Small dimension MAXSS') END IF C C Reading data for snapshot levels: CALL RSEP3I('N4',N4,NT) CALL RSEP3R('D4',D4,DT) CALL RSEP3R('O4',O4,0.) IF(D4.LT.DT) THEN C FD2D-03 CALL ERROR('FD2D-03: Small time step D4 for snapshot levels') END IF IF(O4+FLOAT(N4-1)*D4.GE.(FLOAT(NT)-0.5)*DT) THEN C FD2D-04 CALL ERROR('FD2D-04: Large number N4 of snapshot levels') END IF C C Reading the gridded effective elastic parameters: CALL FDIN(LU1,NX,NZ,AHL,BHL,CHL,AVL,BVL,CVL, * AHS,BHS,CHS,AVS,BVS,CVS,DEN,UX,UY) C C The right, left and bottom sides of the comp. region are C nonreflecting when NSYM=0. The right and bottom are C nonreflecting but the left is symmetrical when NSYM=1. C The latter is available for the incident S-wave only. CALL RSEP3I('NSYM' ,NSYM ,0) C The incidence angle, positive clockwise (degrees). C For vertical incidence the angle is 0. CALL RSEP3R('ANGINC',ANGINC,0.) C Azimuth of the profile with respect to the incident ray (degrees). c For the ray in (x,z) plane choose 0. CALL RSEP3R('AZIPRO',AZIPRO,0.) C Incident-wave velocity (m/s) in bedrock: CALL RSEP3R('ROCKVE',ROCKVE,999999.) C C Opening output receiver files (old form): IF(NREC.GT.0) THEN CALL RSEP3T('OUTU',FOUTU,' ') IF(FOUTU.NE.' ') THEN OPEN(LU3,FILE=FOUTU) END IF CALL RSEP3T('OUTW',FOUTW,' ') IF(FOUTW.NE.' ') THEN OPEN(LU4,FILE=FOUTW) END IF END IF C C Opening output snapshot files: C Form of the output files: CALL RSEP3T('FORM',FORM,'FORMATTED') C Names of the output snapshot files: CALL RSEP3T('SNAP1',FILNA1,' ') CALL RSEP3T('SNAP3',FILNA2,' ') IF(FILNA1.NE.' ') THEN OPEN(LU1,FILE=FILNA1,FORM=FORM) END IF IF(FILNA2.NE.' ') THEN OPEN(LU2,FILE=FILNA2,FORM=FORM) END IF C C Auxiliary parameters C TSWI2=9999. C DO 350 I=1,NREC INREC(I)=LRE(I)+(KRE(I)-1)*NZ 350 CONTINUE C C Preparatory calculations for artificial dumpers C fixed with NZONE=20. C C goto 4001 NZONE=20 C NXV=NX-NZONE NZV=NZ-NZONE NXP=NXV+1 NZP=NZV+1 DO 3200 K=1,NX DO 3100 L=1,NZ IDU=L+(K-1)*NZ DUMP(IDU)=1. 3100 CONTINUE 3200 CONTINUE C C Change to K=1,NXV, and 'C' the loop 3700 if symmetrical left side. C DO 3400 K=NZONE+1,NXV DO 3300 L=NZP,NZ IDU=L+(K-1)*NZ DUMP(IDU)=CDUMP(NZ+1-L) 3300 CONTINUE 3400 CONTINUE DO 3700 K=1,NZONE DO 3500 L=1,NZV IDU=L+(K-1)*NZ DUMP(IDU)=CDUMP(K) 3500 CONTINUE DO 3600 L=NZP,NZ IDU=L+(K-1)*NZ DUMP(IDU)=CDUMP(K)*CDUMP(NZ+1-L) 3600 CONTINUE 3700 CONTINUE C DO 4000 K=NXP,NX DO 3800 L=1,NZV IDU=L+(K-1)*NZ DUMP(IDU)=CDUMP(NX+1-K) 3800 CONTINUE DO 3900 L=NZP,NZ IDU=L+(K-1)*NZ DUMP(IDU)=CDUMP(NX+1-K)*CDUMP(NZ+1-L) 3900 CONTINUE 4000 CONTINUE C 4001 CONTINUE NZ5=NZ-5 NZ6=2*NZ+1 CN=DT/DX CNST=CN*CN/12. NSIZE=NX*NZ LA2=LA1+1 LA3=LA1+2 LA4=LA1+3 C NSIZEA=NSIZE NXA=NX NZA=NZ C IF(ISOURC.NE.2) THEN PI=3.141593 ANG=PI*ANGINC/180. AZI=PI*AZIPRO/180. LZERO=LA3 IF(ANG.GE.0.) KZERO=1 IF(ANG.LT.0.) KZERO=NX SX=SIN(ANG)*COS(AZI)/ROCKVE SY=SIN(ANG)*SIN(AZI)/ROCKVE SZ=COS(ANG)/ROCKVE TIMEXC=SX*FLOAT(NX-1)*DX+SZ*3.*DX + 2.*TP C IF(ISOURC.EQ.1) GOTO 353 C The P-wave incidence UCOM=SIN(ANG) WCOM=-1.*COS(ANG) GOTO 354 C The S-wave incidence 353 UCOM=-1.*COS(ANG) WCOM=-1.*SIN(ANG) 354 CONTINUE ELSE C Explosive point source UCOM=0. WCOM=0. TIMEXC=0. IF(KSOUR.LT.2.OR.KSOUR.GT.NX-1.OR. * LSOUR.LT.2.OR.LSOUR.GT.NZ-1) THEN C FD2D-05 CALL ERROR('FD2D-05: Explosive source exceeds the grid') C Explosive source is composed of 3*3 gridpoints. It thus C cannot be centred at the boundary of the grid. It is also not C recommended to place the source into a dumper. The dumpers C are located along the bottom and sides of the grid and are C NZONE gridpoints thick (NZONE=20 in this version). END IF ICSR=LSOUR+(KSOUR-1)*NZ INU(1)=ICSR-NZ-1 INU(2)=ICSR-NZ INU(3)=ICSR-NZ+1 INU(4)=ICSR+NZ-1 INU(5)=ICSR+NZ INU(6)=ICSR+NZ+1 INW(1)=ICSR-NZ-1 INW(2)=ICSR-1 INW(3)=ICSR+NZ-1 INW(4)=ICSR-NZ+1 INW(5)=ICSR+1 INW(6)=ICSR+NZ+1 C Unit seismic moment CA=(DT/DX)**2/DEN(ICSR)/(4.*DX) CB=0.5*CA UP(1)=CB UP(2)=CA UP(3)=CB UP(4)=-CB UP(5)=-CA UP(6)=-CB WP(1)=CB WP(2)=CA WP(3)=CB WP(4)=-CB WP(5)=-CA WP(6)=-CB ENDIF C DO 360 I=1,NSIZE UX(I)=0. UY(I)=0. WX(I)=0. WY(I)=0. 360 CONTINUE C C....................................................................... C C Starting loop over time levels C DO 2000 M=1,NT WRITE(*,'(3(A,I6))') '+',M,' of',NT,' time levels' MPAR=1 IF((FLOAT(M)/2.-FLOAT(M/2)).LT.0.1) MPAR=2 C C Excitation signal at lines L=LA1 to LA4 C T=FLOAT(M)*DT DO 650 I=1,4 DO 640 K=1,NX EX(I,K)=0. 640 CONTINUE 650 CONTINUE IF(T.GT.TIMEXC) GOTO 900 TRET=T-DT C C We assume the incident wavefront striking the plane (x,z) at t=0 TAUY=SY*0. C DO 800 L=LA1,LA4 I=(L-LA1)+1 TAUZ=SZ*FLOAT(LZERO-L)*DX DO 700 K=1,NX TAUX=SX*FLOAT(K-KZERO)*DX TAU=TAUX+TAUY+TAUZ EX(I,K)=TIMFUN(TRET,TAU,TP) 700 CONTINUE 800 CONTINUE 900 CONTINUE C C Saving boundary-strip values from the 'very old' time level M-1 C IF(MPAR.EQ.1) THEN CALL SAVER(KMAX,LMAX,KLMAX,UY,WY, 1 UK1,UK2,UK3,UNX2,UNX1,UNX,UNZ2,UNZ1,UNZ, 2 WK1,WK2,WK3,WNX2,WNX1,WNX,WNZ2,WNZ1,WNZ) ELSE CALL SAVER(KMAX,LMAX,KLMAX,UX,WX, 1 UK1,UK2,UK3,UNX2,UNX1,UNX,UNZ2,UNZ1,UNZ, 2 WK1,WK2,WK3,WNX2,WNX1,WNX,WNZ2,WNZ1,WNZ) ENDIF C C Computing 'new' time level M+1 C C Regular grid points (i.e. all except nonreflecting boundaries) C C Line L=1 C L=1 DO 1100 K=3,NX-2 CALL INDEX(K,L) IF(MPAR.EQ.1) THEN CALL AUXL1(KLMAX,UX,UY,WX,WY) ELSE CALL AUXL1(KLMAX,UY,UX,WY,WX) ENDIF CALL PARL1(KLMAX,AHS,BVS,BVS,CHS,AHL,BVL,BVL,CHL,DEN) CALL DIFL1(UOLD,U3,U8,U13,U14,U15,U18,U23, 1 W3,W5,W8,W9,W13,W14,W15,W18,W19,W23,W25, 2 UNEW) CALL PARL1(KLMAX,BHS,CVS,AVS,BHS,BHL,CVL,AVL,BHL,DEN) CALL DIFL1(WOLD,W3,W8,W13,W14,W15,W18,W23, 1 U3,U5,U8,U9,U13,U14,U15,U18,U19,U23,U25, 2 WNEW) C IF(MPAR.EQ.1) THEN UY(ICEN)=UNEW WY(ICEN)=WNEW ELSE UX(ICEN)=UNEW WX(ICEN)=WNEW ENDIF 1100 CONTINUE C C Line L=2 C L=2 DO 1200 K=3,NX-2 CALL INDEX(K,L) IF(MPAR.EQ.1) THEN CALL AUXL2(KLMAX,UX,UY,WX,WY) ELSE CALL AUXL2(KLMAX,UY,UX,WY,WX) ENDIF CALL PARL2(KLMAX,AHS,BVS,BVS,CHS,AHL,BVL,BVL,CHL,DEN) CALL DIFL2(UOLD,U3,U8,U12,U13,U14,U15,U18, 1 U23, 2 W3,W5,W7,W8,W9,W12,W13,W14,W15,W17,W18,W19, 3 W23,W25, 4 UNEW) CALL PARL2(KLMAX,BHS,CVS,AVS,BHS,BHL,CVL,AVL,BHL,DEN) CALL DIFL2(WOLD,W3,W8,W12,W13,W14,W15,W18, 1 W23, 2 U3,U5,U7,U8,U9,U12,U13,U14,U15,U17,U18,U19, 3 U23,U25, 4 WNEW) C IF(MPAR.EQ.1) THEN UY(ICEN)=UNEW WY(ICEN)=WNEW ELSE UX(ICEN)=UNEW WX(ICEN)=WNEW ENDIF 1200 CONTINUE C C Lines L=3 to LA1-1 (heterogeneous, in general) C DO 1400 L=3,LA1-1 DO 1300 K=3,NX-2 c NOUT=0 c if(l.eq.3.or.l.eq.4) then c NOUT=1 c write(33,*) 'k,l',k,l c endif c CALL INDEX(K,L) IF(MPAR.EQ.1) THEN CALL AUXLA(KLMAX,KMAX,EX,K,L,LA1,LA2,LA3,LA4,UX,UY,WX,WY,T, 1 TIMEXC,UCOM,WCOM) ELSE CALL AUXLA(KLMAX,KMAX,EX,K,L,LA1,LA2,LA3,LA4,UY,UX,WY,WX,T, 1 TIMEXC,UCOM,WCOM) ENDIF CALL PARL(KLMAX,AHS,BVS,BVS,CHS,AHL,BVL,BVL,CHL,DEN) 1260 CALL DIFL(UOLD,U3,U8,U11,U12,U13,U14,U15,U18, 1 U23, 2 W1,W3,W5,W7,W8,W9,W11,W12,W13,W14,W15,W17,W18, 3 W19,W21,W23,W25, 4 UNEW,NOUT) CALL PARL(KLMAX,BHS,CVS,AVS,BHS,BHL,CVL,AVL,BHL,DEN) CALL DIFL(WOLD,W3,W8,W11,W12,W13,W14,W15,W18, 1 W23, 2 U1,U3,U5,U7,U8,U9,U11,U12,U13,U14,U15,U17,U18, 3 U19,U21,U23,U25, 4 WNEW,NOUT) C 1290 IF(MPAR.EQ.1) THEN UY(ICEN)=UNEW WY(ICEN)=WNEW ELSE UX(ICEN)=UNEW WX(ICEN)=WNEW ENDIF 1300 CONTINUE 1400 CONTINUE C C Lines LA1 to NZ-2 (homogeneous), optional excitation at LA1-LA4 C DO 1600 L=LA1,NZ-2 DO 1500 K=3,NX-2 CALL INDEX(K,L) IF(MPAR.EQ.1) THEN CALL AUXLA(KLMAX,KMAX,EX,K,L,LA1,LA2,LA3,LA4,UX,UY,WX,WY,T, 1 TIMEXC,UCOM,WCOM) ELSE CALL AUXLA(KLMAX,KMAX,EX,K,L,LA1,LA2,LA3,LA4,UY,UX,WY,WX,T, 1 TIMEXC,UCOM,WCOM) ENDIF CALL PARHOM(KLMAX,AHL,BVL,BVL,CHL,DEN) CALL DIFHOM(UOLD,U3,U8,U11,U12,U13,U14,U15,U18, 1 U23, 2 W1,W5,W7,W9,W17, 3 W19,W21,W25, 4 UNEW) CALL PARHOM(KLMAX,BHL,CVL,AVL,BHL,DEN) CALL DIFHOM(WOLD,W3,W8,W11,W12,W13,W14,W15,W18, 1 W23, 2 U1,U5,U7,U9,U17, 3 U19,U21,U25, 4 WNEW) C C IF(MPAR.EQ.1) THEN UY(ICEN)=UNEW WY(ICEN)=WNEW ELSE UX(ICEN)=UNEW WX(ICEN)=WNEW ENDIF 1500 CONTINUE 1600 CONTINUE C C C Nonreflecting boundaries - inner lines C C The right side E-S: CALL FLERIB(KLMAX,LMAX,1,NZA-3,NSIZEA-3*NZA,NSIZEA-2*NZA, 1 NSIZEA-3*NZA,UNX2,UNX1,WNX2,WNX1, 2 UX,UY,WX,WY,AHL,BHL,DEN) C The right side Reynolds C CALL FLERIB(KLMAX,LMAX,1,NZA-3,NSIZEA-3*NZA,NSIZEA-2*NZA, C 1 NSIZEA-3*NZA,UNX,UNX2,WNX,WNX2, C 2 UX,UY,WX,WY,AHL,BHL,DEN) C The left side E-S: CALL FLERIB(KLMAX,LMAX,1,NZA-3,2*NZA,NZA, 1 NZA,UK3,UK2,WK3,WK2, 2 UX,UY,WX,WY,AHL,BHL,DEN) C The bottom side E-S: CALL FBOB(KLMAX,KMAX,4,NXA-3,NZA-2,1,UNZ2,UNZ1,WNZ2,WNZ1, 1 UX,UY,WX,WY,AVL,BVL,DEN) C The right bottom corner E-S: CALL FCORNB(KLMAX,NSIZEA-NZA-1,NZA,NZA,0, 1 UX,UY,WX,WY,AHL,BHL,DEN) C The left bottom corner E-S: CALL FCORNB(KLMAX,2*NZA-1,0,-NZA,1, 1 UX,UY,WX,WY,AHL,BHL,DEN) C C Nonreflecting boundaries - outer lines C C The right side E-S: CALL FLERIB(KLMAX,LMAX,1,NZA-2,NSIZEA-2*NZA,NSIZEA-NZA, 1 NSIZEA-2*NZA,UNX1,UNX,WNX1,WNX, 2 UX,UY,WX,WY,AHL,BHL,DEN) C The right side Reynolds: C CALL FLERIB(KLMAX,LMAX,1,NZA-2,NSIZEA-2*NZA,NSIZEA-NZA, C 1 NSIZEA-2*NZA,UNX2,UNX1,WNX2,WNX1, C 2 UX,UY,WX,WY,AHL,BHL,DEN) C The left side E-S: CALL FLERIB(KLMAX,LMAX,1,NZA-2,NZA,0, 1 0,UK2,UK1,WK2,WK1, 2 UX,UY,WX,WY,AHL,BHL,DEN) C The bottom side E-S: CALL FBOB(KLMAX,KMAX,3,NXA-2,NZA-1,0,UNZ1,UNZ,WNZ1,WNZ, 1 UX,UY,WX,WY,AVL,BVL,DEN) C The right bottom corner E-S: CALL FCORNB(KLMAX,NSIZEA,NZA,NZA,0, 1 UX,UY,WX,WY,AHL,BHL,DEN) C The left bottom corner E-S: CALL FCORNB(KLMAX,NZA,0,-NZA,1, 1 UX,UY,WX,WY,AHL,BHL,DEN) C IF(NSYM.EQ.0) goto 4300 C C Symmetrical bottom c4200 DO 1624 L=NZ-1,NZ c DO 1623 K=3,NX-2 c IND=L+(K-1)*NZ c IF(L.EQ.NZ-1) GOTO 1622 c IF(MPAR.EQ.1) THEN c UY(IND)=UY(IND-4) c WY(IND)=-WY(IND-4) c ELSE c UX(IND)=UX(IND-4) c WX(IND)=-WX(IND-4) c ENDIF c GOTO 1623 c1622 IF(MPAR.EQ.1) THEN c UY(IND)=UY(IND-2) c WY(IND)=-WY(IND-2) c ELSE c UX(IND)=UX(IND-2) c WX(IND)=-WX(IND-2) c ENDIF c1623 CONTINUE c1624 CONTINUE C C Symmetrical left side ( GOOD FOR INCIDENT S WAVE ONLY !!!!!!!!) DO 1620 K=1,2 DO 1610 L=1,NZ IND=L+(K-1)*NZ IF(K.EQ.2) GOTO 1605 IF(MPAR.EQ.1) THEN UY(IND)=UY(IND+4*NZ) WY(IND)=-WY(IND+4*NZ) ELSE UX(IND)=UX(IND+4*NZ) WX(IND)=-WX(IND+4*NZ) ENDIF GOTO 1610 1605 IF(MPAR.EQ.1) THEN UY(IND)=UY(IND+2*NZ) WY(IND)=-WY(IND+2*NZ) ELSE UX(IND)=UX(IND+2*NZ) WX(IND)=-WX(IND+2*NZ) ENDIF 1610 CONTINUE 1620 CONTINUE c improved left side symmetry GOOD FOR INCIDENT S WAVE ONLY ! c For S incidence change u=0 for w=0 !!!!!!!!!! DO 1627 L=NZ-1,NZ K=3 IND=L+(K-1)*NZ IF(MPAR.EQ.1) THEN WY(IND)=0. ELSE WX(IND)=0. ENDIF 1627 CONTINUE C C Symmetrical right side (GOOD FOR INCIDENT S WAVE ONLY !!!!!!!!) c DO 1640 K=NX-1,NX c DO 1630 L=1,NZ c IND=L+(K-1)*NZ c IF(K.EQ.NX-1) GOTO 1625 c IF(MPAR.EQ.1) THEN c UY(IND)=UY(IND-4*NZ) c WY(IND)=-WY(IND-4*NZ) c ELSE c UX(IND)=UX(IND-4*NZ) c WX(IND)=-WX(IND-4*NZ) c ENDIF c GOTO 1630 c1625 IF(MPAR.EQ.1) THEN c UY(IND)=UY(IND-2*NZ) c WY(IND)=-WY(IND-2*NZ) c ELSE c UX(IND)=UX(IND-2*NZ) c WX(IND)=-WX(IND-2*NZ) c ENDIF c1630 CONTINUE c1640 CONTINUE c improved right side symmetry GOOD FOR INCIDENT S WAVE ONLY ! c For S incid change u=0 for w=0 !!!!!!!!!! c DO 1628 L=NZ-1,NZ c K=NX-2 c IND=L+(K-1)*NZ c IF(MPAR.EQ.1) THEN c WY(IND)=0. c ELSE c WX(IND)=0. c ENDIF c1628 CONTINUE C C 4300 CONTINUE C C Line source (a body force equivalent) C IF(ISOURC.NE.2) GOTO 1700 IF(T.GT.(2.*TP)) GOTO 1700 TIMF=TIMFUN(T,0.,TP) IF(MPAR.EQ.1) THEN C C Unit seismic force * CA=(DT/DX)**2/DEN(ICSR) C Single horizontal force * UY(ICSR)=UY(ICSR)-CA*TIMF C Single vertical force * WY(ICSR)=WY(ICSR)-CA*TIMF C DO 1650 I=1,6 CC Centre of compression (explosive line source) UY(INU(I))=UY(INU(I))-UP(I)*TIMF WY(INW(I))=WY(INW(I))-WP(I)*TIMF CC Centre of rotation CC UY(INW(I))=UY(INW(I))-WP(I)*TIMF CC WY(INU(I))=WY(INU(I))+UP(I)*TIMF CC Double couple CC UY(INW(I))=UY(INW(I))-WP(I)*TIMF CC WY(INU(I))=WY(INU(I))-UP(I)*TIMF CC Single horizontal force (extended) CC UY(INU(I))=UY(INU(I))-UP(I)*TIMF 1650 CONTINUE ELSE C C Single horizontal force * UX(ICSR)=UX(ICSR)-CA*TIMF C Single vertical force * WX(ICSR)=WX(ICSR)-CA*TIMF C DO 1660 I=1,6 CC Centre of compression UX(INU(I))=UX(INU(I))-UP(I)*TIMF WX(INW(I))=WX(INW(I))-WP(I)*TIMF CC Centre of rotation CC UX(INW(I))=UX(INW(I))-WP(I)*TIMF CC WX(INU(I))=WX(INU(I))+UP(I)*TIMF CC Double couple CC UX(INW(I))=UX(INW(I))-WP(I)*TIMF CC WX(INU(I))=WX(INU(I))-UP(I)*TIMF CC Single horizontal force (extended) CC UX(INU(I))=UX(INU(I))-UP(I)*TIMF 1660 CONTINUE ENDIF 1700 CONTINUE C C C Saving results for receivers: C IF(NREC.GT.0) THEN DO 1800 I=1,NREC J=(I-1)*NT+M IF(MPAR.EQ.1) THEN USS(J)=UY(INREC(I)) WSS(J)=WY(INREC(I)) ELSE USS(J)=UX(INREC(I)) WSS(J)=WX(INREC(I)) END IF IF(D1.LT.0.) THEN USS(J)=-USS(J) END IF IF(D3.LT.0.) THEN WSS(J)=-WSS(J) END IF 1800 CONTINUE if(mod(m,4).eq.0) then C Saving results for receivers at each 4th time level DO 1801 I=1,NREC IF(MPAR.EQ.1) THEN ULEV(I)=UY(INREC(I)) WLEV(I)=WY(INREC(I)) ELSE ULEV(I)=UX(INREC(I)) WLEV(I)=WX(INREC(I)) END IF 1801 CONTINUE IF(FOUTU.NE.' ') THEN WRITE(LU3,'(8E10.4)') (ULEV(I),I=1,NREC) ENDIF IF(FOUTW.NE.' ') THEN WRITE(LU4,'(8E10.4)') (WLEV(I),I=1,NREC) ENDIF endif ENDIF C C Writing a snapshot: I4=NINT((FLOAT(M-1)*DT-O4)/D4)+1 IT=NINT((FLOAT(I4-1)*D4+O4)/DT)+1 C I4 is the index of the snapshot nearest to the M-th time level, C IT is the time level nearest to the I4-th snapshot. IF(M.EQ.IT.AND.I4.LE.N4) THEN IF(FILNA1.NE.' ') THEN IF(MPAR.EQ.1) THEN CALL FDOUT1(LU1,NX,NZ,D1,' ',UY,WORK) ELSE CALL FDOUT1(LU1,NX,NZ,D1,' ',UX,WORK) END IF END IF IF(FILNA2.NE.' ') THEN IF(MPAR.EQ.1) THEN CALL FDOUT1(LU2,NX,NZ,D3,' ',WY,WORK) ELSE CALL FDOUT1(LU2,NX,NZ,D3,' ',WX,WORK) END IF END IF END IF C C Artificial dumpers at right and bottom C c IF(T.LT.TSWI2) GOTO 2000 DO 1950 IA=1,NSIZE UY(IA)=UY(IA)*DUMP(IA) UX(IA)=UX(IA)*DUMP(IA) WY(IA)=WY(IA)*DUMP(IA) WX(IA)=WX(IA)*DUMP(IA) 1950 CONTINUE C 2000 CONTINUE c call gettim(ihod2,imin2,isec2,idum) c write(*,*) ihod1,imin1,isec1 c write(*,*) ihod2,imin2,isec2 C IF(FILNA1.NE.' ') THEN CLOSE(LU1) END IF IF(FILNA2.NE.' ') THEN CLOSE(LU2) END IF C C Writing synthetic seismograms: IF(NREC.GT.0) THEN CALL RSEP3T('SS',OUTSS,'fd2d.gse') IF(OUTSS.NE.' ') THEN OPEN(LU1,FILE=OUTSS) CALL WGSE1(LU1,' ') DO 9000 I=1,NREC X=O1+D1*FLOAT(KRE(I)-1) Y=0. Z=O3+D3*FLOAT(LRE(I)-1) TL=0. CALL WGSE2(LU1,REC(I),' ',1,X,Y,Z,TL,DT,NT,USS(NT*(I-1)+1)) CALL WGSE2(LU1,REC(I),' ',3,X,Y,Z,TL,DT,NT,WSS(NT*(I-1)+1)) 9000 CONTINUE CALL WGSE3(LU1) CLOSE(LU1) END IF END IF C RETURN END C C======================================================================= C C C FUNCTION CDUMP(I) C C Purpose: C Spatial variation of artificial dumpers C COMMON /DMP/ NZONE PI=3.141593 A=PI*FLOAT(I)/FLOAT(NZONE) CDUMP=0.98+0.01*(1.-COS(A)) C CDUMP=0.94+0.03*(1.-COS(A)) RETURN END C DEBUG SUBCHK C END DEBUG C C======================================================================= C C C SUBROUTINE AUXL1(KLMAX,AX,AY,BX,BY) C C Purpose: C Transfer of the displacements from arrays UX,UY,etc. C into auxiliary variables UOLD,U1,U2,etc. C (Selecting values from vicinity of a grid point at line L=1.) C DIMENSION II(25),IJ(25) DIMENSION AX(KLMAX),AY(KLMAX),BX(KLMAX),BY(KLMAX) COMMON /INDI/ ICEN,II,IJ, 1 isa,isb,isc,isd,ise,isf,isg,ish, 2 ila,ilb,ilc,ild,ile,ilf,ilg,ilh COMMON /DSPLC/ UOLD,WOLD, 1 U1,U2,U3,U4,U5,U6,U7,U8,U9,U10,U11,U12,U13,U14, 2 U15,U16,U17,U18,U19,U20,U21,U22,U23,U24,U25, 3 W1,W2,W3,W4,W5,W6,W7,W8,W9,W10,W11,W12,W13,W14, 4 W15,W16,W17,W18,W19,W20,W21,W22,W23,W24,W25 C UOLD=AY(ICEN) U3=AX(II(3)) U5=AX(II(5)) U8=AX(II(8)) U9=AX(II(9)) U13=AX(II(13)) U14=AX(II(14)) U15=AX(II(15)) U18=AX(II(18)) U19=AX(II(19)) U23=AX(II(23)) U25=AX(II(25)) WOLD=BY(ICEN) W3=BX(II(3)) W5=BX(II(5)) W8=BX(II(8)) W9=BX(II(9)) W13=BX(II(13)) W14=BX(II(14)) W15=BX(II(15)) W18=BX(II(18)) W19=BX(II(19)) W23=BX(II(23)) W25=BX(II(25)) RETURN END C DEBUG SUBCHK C END DEBUG C C======================================================================= C C C SUBROUTINE AUXL2(KLMAX,AX,AY,BX,BY) C C Purpose: C Transfer of the displacements from arrays UX,UY,etc. C into auxiliary variables UOLD,U1,U2,etc. C (Selecting values from vicinity of a grid point at line L=2.) C DIMENSION II(25),IJ(25) DIMENSION AX(KLMAX),AY(KLMAX),BX(KLMAX),BY(KLMAX) COMMON /INDI/ ICEN,II,IJ, 1 isa,isb,isc,isd,ise,isf,isg,ish, 2 ila,ilb,ilc,ild,ile,ilf,ilg,ilh COMMON /DSPLC/ UOLD,WOLD, 1 U1,U2,U3,U4,U5,U6,U7,U8,U9,U10,U11,U12,U13,U14, 2 U15,U16,U17,U18,U19,U20,U21,U22,U23,U24,U25, 3 W1,W2,W3,W4,W5,W6,W7,W8,W9,W10,W11,W12,W13,W14, 4 W15,W16,W17,W18,W19,W20,W21,W22,W23,W24,W25 C UOLD=AY(ICEN) U3=AX(II(3)) U5=AX(II(5)) U7=AX(II(7)) U8=AX(II(8)) U9=AX(II(9)) U12=AX(II(12)) U13=AX(II(13)) U14=AX(II(14)) U15=AX(II(15)) U17=AX(II(17)) U18=AX(II(18)) U19=AX(II(19)) U23=AX(II(23)) U25=AX(II(25)) WOLD=BY(ICEN) W3=BX(II(3)) W5=BX(II(5)) W7=BX(II(7)) W8=BX(II(8)) W9=BX(II(9)) W12=BX(II(12)) W13=BX(II(13)) W14=BX(II(14)) W15=BX(II(15)) W17=BX(II(17)) W18=BX(II(18)) W19=BX(II(19)) W23=BX(II(23)) W25=BX(II(25)) RETURN END C DEBUG SUBCHK C END DEBUG C C======================================================================= C C C SUBROUTINE AUXLA(KLMAX,KMAX,EX,K,L,LA1,LA2,LA3,LA4,AX,AY,BX,BY, 1 T,TIMEXC,UCOM,WCOM) C C Purpose: C Transfer of the displacements from arrays UX,UY,etc. C into auxiliary variables UOLD,U1,U2,etc. C (Selecting values from vicinity of a grid point at line L>2.) C DIMENSION II(25),IJ(25) DIMENSION AX(KLMAX),AY(KLMAX),BX(KLMAX),BY(KLMAX) DIMENSION EX(4,KMAX) COMMON /INDI/ ICEN,II,IJ, 1 isa,isb,isc,isd,ise,isf,isg,ish, 2 ila,ilb,ilc,ild,ile,ilf,ilg,ilh COMMON /DSPLC/ UOLD,WOLD, 1 U1,U2,U3,U4,U5,U6,U7,U8,U9,U10,U11,U12,U13,U14, 2 U15,U16,U17,U18,U19,U20,U21,U22,U23,U24,U25, 3 W1,W2,W3,W4,W5,W6,W7,W8,W9,W10,W11,W12,W13,W14, 4 W15,W16,W17,W18,W19,W20,W21,W22,W23,W24,W25 C UOLD=AY(ICEN) U1=AX(II(1)) U3=AX(II(3)) U5=AX(II(5)) U7=AX(II(7)) U8=AX(II(8)) U9=AX(II(9)) U11=AX(II(11)) U12=AX(II(12)) U13=AX(II(13)) U14=AX(II(14)) U15=AX(II(15)) U17=AX(II(17)) U18=AX(II(18)) U19=AX(II(19)) U21=AX(II(21)) U23=AX(II(23)) U25=AX(II(25)) WOLD=BY(ICEN) W1=BX(II(1)) W3=BX(II(3)) W5=BX(II(5)) W7=BX(II(7)) W8=BX(II(8)) W9=BX(II(9)) W11=BX(II(11)) W12=BX(II(12)) W13=BX(II(13)) W14=BX(II(14)) W15=BX(II(15)) W17=BX(II(17)) W18=BX(II(18)) W19=BX(II(19)) W21=BX(II(21)) W23=BX(II(23)) W25=BX(II(25)) IF(L.LT.LA1.OR.L.GT.LA4) GOTO 500 IF(T.GT.TIMEXC) GOTO 500 IF(L.EQ.LA1) GOTO 100 IF(L.EQ.LA2) GOTO 200 IF(L.EQ.LA3) GOTO 300 IF(L.EQ.LA4) GOTO 400 100 DIF=EX(3,K-2) DIFU=DIF*UCOM DIFW=DIF*WCOM U5=U5+DIFU W5=W5+DIFW DIF=EX(3,K) DIFU=DIF*UCOM DIFW=DIF*WCOM U15=U15+DIFU W15=W15+DIFW DIF=EX(3,K+2) DIFU=DIF*UCOM DIFW=DIF*WCOM U25=U25+DIFU W25=W25+DIFW GOTO 500 200 DIF=EX(3,K-1) DIFU=DIF*UCOM DIFW=DIF*WCOM U9=U9+DIFU W9=W9+DIFW DIF=EX(3,K) DIFU=DIF*UCOM DIFW=DIF*WCOM U14=U14+DIFU W14=W14+DIFW DIF=EX(3,K+1) DIFU=DIF*UCOM DIFW=DIF*WCOM U19=U19+DIFU W19=W19+DIFW DIF=EX(4,K-2) DIFU=DIF*UCOM DIFW=DIF*WCOM U5=U5+DIFU W5=W5+DIFW DIF=EX(4,K) DIFU=DIF*UCOM DIFW=DIF*WCOM U15=U15+DIFU W15=W15+DIFW DIF=EX(4,K+2) DIFU=DIF*UCOM DIFW=DIF*WCOM U25=U25+DIFU W25=W25+DIFW GOTO 500 300 DIF=EX(1,K-2) DIFU=DIF*UCOM DIFW=DIF*WCOM U1=U1-DIFU W1=W1-DIFW DIF=EX(1,K) DIFU=DIF*UCOM DIFW=DIF*WCOM U11=U11-DIFU W11=W11-DIFW DIF=EX(1,K+2) DIFU=DIF*UCOM DIFW=DIF*WCOM U21=U21-DIFU W21=W21-DIFW DIF=EX(2,K-1) DIFU=DIF*UCOM DIFW=DIF*WCOM U7=U7-DIFU W7=W7-DIFW DIF=EX(2,K) DIFU=DIF*UCOM DIFW=DIF*WCOM U12=U12-DIFU W12=W12-DIFW DIF=EX(2,K+1) DIFU=DIF*UCOM DIFW=DIF*WCOM U17=U17-DIFU W17=W17-DIFW GOTO 500 400 DIF=EX(2,K-2) DIFU=DIF*UCOM DIFW=DIF*WCOM U1=U1-DIFU W1=W1-DIFW DIF=EX(2,K) DIFU=DIF*UCOM DIFW=DIF*WCOM U11=U11-DIFU W11=W11-DIFW DIF=EX(2,K+2) DIFU=DIF*UCOM DIFW=DIF*WCOM U21=U21-DIFU W21=W21-DIFW 500 CONTINUE RETURN END C DEBUG SUBCHK C END DEBUG C C======================================================================= C C C SUBROUTINE INDEX(K,L) C C Purpose: C Computing 1-D index ICEN of a grid point K,L, C computing indices II(1)-II(25) of the grid-point neighbours, C computing indices isa,isb,...ish of the neighbouring small legs c (h/2), C computing indices ila,ilb,...ilh of the neighbouring large legs C (h). C COMMON /INDI/ ICEN,II,IJ, 1 isa,isb,isc,isd,ise,isf,isg,ish, 2 ila,ilb,ilc,ild,ile,ilf,ilg,ilh COMMON /BASE/ MPAR,NSIZE,NX,NZ,NZ5,NZ6,CN,CNST DIMENSION II(25),IJ(25) c c ....f.......g.... c . f . g . c . f . g . c ........J.I C E N c . c C b . c . c E b . c ....c...N...b.... c c c ................. c . . . c eeeeeee . hhhhhhh c ........J.I C E N c . C . c ddddddd E aaaaaaa c ........N........ c c C ICEN=L+(K-1)*NZ ICEN2=ICEN-2 JCEN=ICEN-IFIX(FLOAT(ICEN)/FLOAT(NZ)) C DO 100 IPOM=1,25 IR=IFIX((FLOAT(IPOM)-0.00001)/5.) NZMUL=NZ5*IR-NZ6+IPOM II(IPOM)=ICEN2+NZMUL 100 CONTINUE C C Indices of the grid legs, see Fig. in PARL c ila, ilb, etc. = indices for legs at grid lines C ila=icen ilb=jcen ild=icen-nz ilf=jcen-1 c c isa, isb, etc. = indices for INTERlegs between grid lines c isa=icen-(k-1) isd=isa-(nz-1) ise=isd-1 ish=isa-1 C isb=jcen isc=jcen-(nz-1) isf=isc-1 isg=isb-1 c RETURN END C C======================================================================= C C C SUBROUTINE PARL1(KLMAX,P2,Q2,R2,S2,P4,Q4,R4,S4,D4) C C Purpose: C Transfer of the effective parameters from arrays AHL,AHS,etc. C into auxiliary variables ........ C (selecting values from vicinity of a grid point at line L=1) C COMMON /PPAR/ pahs,pdes,rbcs,rfgs, 1 shs,sas,ses,sds,qgs,qfs,qbs,qcs, 2 pahl,pdel,rbcl,rfgl, 3 shl,sal,sel,sdl,qgl,qfl,qbl,qcl, 4 phom,qhom,rhom,shom,dnst COMMON /INDI/ ICEN,II,IJ, 1 isa,isb,isc,isd,ise,isf,isg,ish, 2 ila,ilb,ilc,ild,ile,ilf,ilg,ilh DIMENSION P2(KLMAX),Q2(KLMAX),R2(KLMAX),S2(KLMAX) DIMENSION P4(KLMAX),Q4(KLMAX),R4(KLMAX),S4(KLMAX),D4(KLMAX) DIMENSION II(25),IJ(25) C c ##################################################### c Density for the central point dnst=d4(icen) c ##################################################### c c The parameters for the non-mixed derivatives Dx(p Dx) c horizontal legs de and ah pdes=p4(ild) pahs=p4(ila) c c The parameters for the non-mixed derivatives Dz(r Dz) c vertical legs fg and bc rfgs=0. rbcs=r4(ilb) c c The parameters for the mixed derivatives Dx(s Dz) c horizontal (inter)legs h,a,e,d shs=0. sas=s2(isa) ses=0. sds=s2(isd) c c equating legs h-a and d-e c shs=s4(ila) c sas=shs c sds=s4(ild) c ses=sds c c The parameters for the mixed derivatives Dz(q Dx) c vertical (inter)legs f,c,g,b qfs=0. qcs=q2(isc) qgs=0. qbs=q2(isb) c c equating legs f-g and b-c c qfs=0. c qgs=0. c qbs=q4(ilb) c qcs=qbs c c ###################################################### c The parameters for the outer formula (4th order) c not yet OK c ###################################################### c c RETURN END C C======================================================================= C C C SUBROUTINE PARL2(KLMAX,P2,Q2,R2,S2,P4,Q4,R4,S4,D4) C C Purpose: C Transfer of the effective parameters from arrays AHL,AHS,etc. C into auxiliary variables ........ C (selecting values from vicinity of a grid point at line L=2) C COMMON /PPAR/ pahs,pdes,rbcs,rfgs, 1 shs,sas,ses,sds,qgs,qfs,qbs,qcs, 2 pahl,pdel,rbcl,rfgl, 3 shl,sal,sel,sdl,qgl,qfl,qbl,qcl, 4 phom,qhom,rhom,shom,dnst COMMON /INDI/ ICEN,II,IJ, 1 isa,isb,isc,isd,ise,isf,isg,ish, 2 ila,ilb,ilc,ild,ile,ilf,ilg,ilh DIMENSION P2(KLMAX),Q2(KLMAX),R2(KLMAX),S2(KLMAX) DIMENSION P4(KLMAX),Q4(KLMAX),R4(KLMAX),S4(KLMAX),D4(KLMAX) DIMENSION II(25),IJ(25) C c ##################################################### c Density for the central point dnst=d4(icen) c ################################################## c c The parameters for the non-mixed derivatives Dx(p Dx) c horizontal legs de and ah pdes=p4(ild) pahs=p4(ila) c c The parameters for the non-mixed derivatives Dz(r Dz) c vertical legs fg and bc rfgs=r4(ilf) rbcs=r4(ilb) c c The parameters for the mixed derivatives Dx(s Dz) c horizontal (inter)legs h,a,e,d c shs=s2(ish) c sas=s2(isa) c ses=s2(ise) c sds=s2(isd) c equating legs h-a and d-e shs=s4(ila) c sas=shs sds=s4(ild) c ses=sds c c The parameters for the mixed derivatives Dz(q Dx) c vertical (inter)legs f,c,g,b c qfs=q2(isf) c qcs=q2(isc) c qgs=q2(isg) c qbs=q2(isb) c equating legs f-g and b-c qfs=q4(ilf) c qgs=qfs qbs=q4(ilb) c qcs=qbs c c ###################################################### c The parameters for the outer (large) square - 4th order c not yet OK c ###################################################### c RETURN END C C======================================================================= C C C SUBROUTINE PARL(KLMAX,P2,Q2,R2,S2,P4,Q4,R4,S4,D4) C C Purpose: C Transfer of the effective parameters from arrays AH,BH,etc. C into auxiliary variables ........ C (selecting values from vicinity of a grid point at line L>2) C COMMON /PPAR/ pahs,pdes,rbcs,rfgs, 1 shs,sas,ses,sds,qgs,qfs,qbs,qcs, 2 pahl,pdel,rbcl,rfgl, 3 shl,sal,sel,sdl,qgl,qfl,qbl,qcl, 4 phom,qhom,rhom,shom,dnst COMMON /INDI/ ICEN,II,IJ, 1 isa,isb,isc,isd,ise,isf,isg,ish, 2 ila,ilb,ilc,ild,ile,ilf,ilg,ilh DIMENSION P2(KLMAX),Q2(KLMAX),R2(KLMAX),S2(KLMAX) DIMENSION P4(KLMAX),Q4(KLMAX),R4(KLMAX),S4(KLMAX),D4(KLMAX) DIMENSION II(25),IJ(25) C c ##################################################### c Density for the central point, formally written c in the array indexed as 'large' horizontal parameters c and/or the array of the displacements dnst=d4(icen) c ##################################################### c c The parameters for the inner formula (2nd order) c - denotation of the used arrays: p4,q4,r4,s4 c - they correspond to the legs c coinciding with grid lines c - numbering of the legs by indices c ila, ilb, ild, ilf c - denotation of the used arrays: p2,q2,r2,s2 c - they correspond to the legs c half a way c between the grid lines c - numbering of the legs by indices c isa, isb, ... ish c ################################################## c c The parameters for the non-mixed derivatives Dx(p Dx) c horizontal legs de and ah pdes=p4(ild) pahs=p4(ila) c c The parameters for the non-mixed derivatives Dz(r Dz) c vertical legs fg and bc rfgs=r4(ilf) rbcs=r4(ilb) c c The parameters for the mixed derivatives Dx(s Dz) c horizontal (inter)legs h,a,e,d c shs=s2(ish) c sas=s2(isa) c ses=s2(ise) c sds=s2(isd) c equating legs h-a and d-e shs=s4(ila) c sas=shs sds=s4(ild) c ses=sds c c The parameters for the mixed derivatives Dz(q Dx) c vertical (inter)legs f,c,g,b c qfs=q2(isf) c qcs=q2(isc) c qgs=q2(isg) c qbs=q2(isb) c equating legs f-g and b-c qfs=q4(ilf) c qgs=qfs qbs=q4(ilb) c qcs=qbs c ###################################################### c The parameters for the outer formula (4th order) c for EQUATED legs only !!!! c ###################################################### c c The parameters for the non-mixed derivatives c double legs de and ah c pdel=1./(1./p4(ild)+1./p4(ild-nz)) c pahl=1./(1./p4(ila)+1./p4(ila+nz)) c c The parameters for the non-mixed derivatives c double legs fg and bc c rfgl=1./(1./r4(ilf)+1./r4(ilf-1)) c rbcl=1./(1./r4(ilb)+1./r4(ilb+1)) c c The parameters for the mixed derivatives c legs h,a,e,d c shl=1./(1./s4(ila)+1./s4(ila+nz)) cc sal=shl c sdl=1./(1./s4(ild)+1./s4(ild-nz)) cc sel=sdl c c The parameters for the mixed derivatives c legs f,c,g,b c qfl=1./(1./q4(ilf)+1./q4(ilf-1)) cc qgl=qfl c qbl=1./(1./q4(ilb)+1./q4(ilb+1)) cc qcl=qbl c RETURN END C DEBUG SUBCHK C END DEBUG C C======================================================================= C C C SUBROUTINE PARHOM(KLMAX,P4,Q4,R4,S4,D4) C C Purpose: C Analogous to PARL but for a grid point with a completely C homogeneous vicinity C COMMON /PPAR/ pahs,pdes,rbcs,rfgs, 1 shs,sas,ses,sds,qgs,qfs,qbs,qcs, 2 pahl,pdel,rbcl,rfgl, 3 shl,sal,sel,sdl,qgl,qfl,qbl,qcl, 4 phom,qhom,rhom,shom,dnst COMMON /INDI/ ICEN,II,IJ, 1 isa,isb,isc,isd,ise,isf,isg,ish, 2 ila,ilb,ilc,ild,ile,ilf,ilg,ilh DIMENSION P4(KLMAX),Q4(KLMAX),R4(KLMAX),S4(KLMAX) DIMENSION D4(KLMAX) DIMENSION II(25),IJ(25) C dnst=d4(icen) c phom=p4(ila) qhom=q4(ilb) rhom=r4(ilb) shom=s4(ila) RETURN END C C======================================================================= C C C SUBROUTINE DIFL1 (AOLD,A3,A8,A13,A14,A15, 1 A18,A23, 2 B3,B5,B8,B9,B13,B14,B15, 3 B18,B19,B23,B25, 4 DISNEW) C C Purpose: C The finite-difference scheme for a grid point at the free surface C at line L=1 c COMMON /PPAR/ pahs,pdes,rbcs,rfgs, 1 shs,sas,ses,sds,qgs,qfs,qbs,qcs, 2 pahl,pdel,rbcl,rfgl, 3 shl,sal,sel,sdl,qgl,qfl,qbl,qcl, 4 phom,qhom,rhom,shom,dnst COMMON /BASE/ MPAR,NSIZE,NX,NZ,NZ5,NZ6,CN,CNST C c inner part (step h in displ.): non mixed + 0.25*(mixed) c c Attention: Zero-valued parameters have been defined in PARL1. c They appear here only formally. c c pintp=pahs*(a18-a13)-pdes*(a13-a8) pintr=rbcs*(a14-a13)+rfgs*(0.-a13) pints=0.25*( shs*(b13+b18-0. -0. )+sas*(-b13-b18+b19+b14)- 3 ses*(b8+b13-0.-0. ) -sds*(-b8-b13+b14+b9) ) pintq=0.25*(-qgs*(0. +b18-0. -b13)-qfs*(-b8-0.+0. +b13)+ 5 qbs*(b18+b19-b13-b14)+qcs*(-b9-b8+b13+b14) ) c c outer part (step 2h in displ.): non mixed + 0.25*(mixed) c (FOR 4th order ONLY) c c poutp=pahl*(a23-a13)-pdel*(a13-a3) c poutr=rbcl*(a15-a13)+rfgl*(0.-a13) c pouts=0.25*( shl*(b13+b23-0. -0. )+sal*(-b13-b23+b25+b15)- c 3 sel*(b3+b13-0.-0. ) -sdl*(-b3-b13+b15+b5) ) c poutq=0.25*(-qgl*(0. +b23-0. -b13)-qfl*(-b3-0.+0. +b13)+ c 5 qbl*(b23+b25-b13-b15)+qcl*(-b5-b3+b13+b15) ) c c 4th order scheme c sleft=16.*pint-pout c c 2nd order (POUT shouldn't be computed in this case) c sleft=12.*pint c sleftp=12.*pintp sleftr=12.*pintr slefts=12.*pints sleftq=12.*pintq c sleft=sleftp+sleftr+slefts+sleftq c c DISNEW=2.*CNST*sleft+2.*A13-AOLD DISNEW=0. IF(DNST.NE.0.) DISNEW=(CNST/DNST)*sleft+2.*A13-AOLD RETURN END C C======================================================================= C C C SUBROUTINE DIFL2 (AOLD,A3,A8,A12,A13,A14,A15, 1 A18,A23, 2 B3,B5,B7,B8,B9,B12,B13,B14,B15,B17, 3 B18,B19,B23,B25, 4 DISNEW) C C Purpose: C The finite-difference scheme for a grid point in a heterogeneous C medium at line L=2 C COMMON /PPAR/ pahs,pdes,rbcs,rfgs, 1 shs,sas,ses,sds,qgs,qfs,qbs,qcs, 2 pahl,pdel,rbcl,rfgl, 3 shl,sal,sel,sdl,qgl,qfl,qbl,qcl, 4 phom,qhom,rhom,shom,dnst COMMON /BASE/ MPAR,NSIZE,NX,NZ,NZ5,NZ6,CN,CNST C c inner part (step h): non mixed + 0.25*(mixed) pintp=pahs*(a18-a13)-pdes*(a13-a8) pintr=rbcs*(a14-a13)+rfgs*(a12-a13) c pints=0.25*( shs*(b13+b18-b12-b17)+sas*(-b13-b18+b19+b14)- c 3 ses*(b8+b13-b7-b12) -sds*(-b8-b13+b14+b9) ) c pintq=0.25*(-qgs*(b17+b18-b12-b13)-qfs*(-b8-b7+b12+b13)+ c 5 qbs*(b18+b19-b13-b14)+qcs*(-b9-b8+b13+b14) ) c mixed derivative with equated legs pints=0.25*( shs*(-b12-b17+b19+b14)- 3 sds*(-b7-b12+b14+b9) ) pintq=0.25*(-qfs*(b17+b18-b8-b7)+ 5 qbs*(b18+b19-b9-b8) ) c outer part (step 2h): non mixed + 0.25*(mixed) c (FOR 4th order ONLY) c c !!!!!!!!!! Attention: Here rfgl, qgl, qfl are problematic !!!!!!!!!!!! c !!!!!!!!!! because of the coincidence of legs F and G with the surface c !!!!!!!!!! Possible actions: redefine them in PARL2, !!!!!!!!!!!!!!!!! c !!!!!!!!!! or use only the 2nd order formula here !!!!!!!!!!!!!!!!!!!! c c poutp=pahl*(a23-a13)-pdel*(a13-a3) c poutr=rbcl*(a15-a13)+rfgl*(0.-a13) c pouts=0.25*( shl*(b13+b23-0. -0. )+sal*(-b13-b23+b25+b15)- c 3 sel*(b3+b13-0.-0. ) -sdl*(-b3-b13+b15+b5) ) c poutq=0.25*(-qgl*(0. +b23-0. -b13)-qfl*(-b3-0.+0. +b13)+ c 5 qbl*(b23+b25-b13-b15)+qcl*(-b5-b3+b13+b15) ) c c 4th order scheme c sleft=16.*pint-pout c c 2nd order (POUT shouldn't be computed in this case) c sleft=12.*pint c c sleftp=12.*pintp sleftr=12.*pintr slefts=12.*pints sleftq=12.*pintq c sleft=sleftp+sleftr+slefts+sleftq c c DISNEW=CNST*sleft+2.*A13-AOLD DISNEW=0. IF(DNST.NE.0.) DISNEW=(CNST/DNST)*sleft+2.*A13-AOLD RETURN END C C======================================================================= C C C SUBROUTINE DIFL (AOLD,A3,A8,A11,A12,A13,A14,A15, 1 A18,A23, 2 B1,B3,B5,B7,B8,B9,B11,B12,B13,B14,B15,B17, 3 B18,B19,B21,B23,B25, 4 DISNEW,NOUT) C C Purpose: C The finite-difference scheme for a grid point in a heterogeneous C medium at line L>2 c c Displacement numbering when updating the central point (#13): c 1 6 11 16 21 c 2 7 12 17 22 c 3 8 13 18 23 c 4 9 14 19 24 A1, A2, etc. for non mixed derivatives, c 5 10 15 20 25 B1, B2, etc. for mixed derivatives. c c c c General denotation p q r s A B c means c when updating u: ah bv bv ch u w c when updating w: bh cv av bh w u c c here 1st char.: a=lambda+2mu, b=mu, c=lambda c 2nd char.: v=vertical legs, h=horizontal legs c c c c Denotation of the effective parameters for non-mixed derivatives: c c |---------|---------| 1st char.: the parameter (p or r) c | r | 2nd and 3rd char.: double legs (f and g, c | f | a and h, c | g | etc.) c ---pde----|---pah---- c | r | (p Ax)x, (r Az)z c | b | c | c | c |---------|---------| c c c c c Denotation of the effective parameters for mixed derivatives: c c |---------|---------| 1st char.: the parameter (q or s) c | qf | qg | 2nd char.: the leg (a,b,...g,h) c | qf | qg | c | qf | qg | c |---------|---------| c | qc | qb | (q Bx)z c | qc | qb | c | qc | qb | c |---------|---------| c c Denotation of the effective parameters for mixed derivatives: c c |---------|---------| 1st char.: the parameter (q or s) c | | | 2nd char.: the leg (a,b,...g,h) c | se se | sh sh | c | | | c |---------|---------| c | | | (s Bz)x c | sd sd | sa sa | c | | | c |---------|---------| c c c When using the parameters in the 'inner' formulas (step h) c they are denoted shs, sas, etc. c c When using the parameters in the 'outer' formulas (step 2h) c they are denoted shl, sal, etc. c c c c configuration of the grid points (for the displacements) c and the grid legs (for the parameters): c c c 1......6......11......16......21 c . f . g . c . f . g . c . f . g . c 2 eeee 7 eeee 12 hhhh 17 hhhh 22 Grid legs a,b,...g,h c . f f . g g . for the 'small' c . f ee.eee.hh..hh g . and 'large' squares, c . f f . g g . used in the inner and c 3......8......13......18......23 outer formulas, resp. c . d c . b b . c . d dd.ddd.aa..aa b . c . d c . b b . c 4 dddd 9 dddd 14 aaaa 19 aaaa 24 c . d . b . c . d . b . c . d . b . c 5.....10......15......20......25 c c c c CNST=dt*dt/dx*dx/12 c c c c C COMMON /PPAR/ pahs,pdes,rbcs,rfgs, 1 shs,sas,ses,sds,qgs,qfs,qbs,qcs, 2 pahl,pdel,rbcl,rfgl, 3 shl,sal,sel,sdl,qgl,qfl,qbl,qcl, 4 phom,qhom,rhom,shom,dnst COMMON /BASE/ MPAR,NSIZE,NX,NZ,NZ5,NZ6,CN,CNST C c inner part (step h): non mixed + 0.25*(mixed) pintp=pahs*(a18-a13)-pdes*(a13-a8) pintr=rbcs*(a14-a13)+rfgs*(a12-a13) c pints=0.25*( shs*(b13+b18-b12-b17)+sas*(-b13-b18+b19+b14)- c 3 ses*(b8+b13-b7-b12) -sds*(-b8-b13+b14+b9) ) c pintq=0.25*(-qgs*(b17+b18-b12-b13)-qfs*(-b8-b7+b12+b13)+ c 5 qbs*(b18+b19-b13-b14)+qcs*(-b9-b8+b13+b14) ) c mixed derivative with equated legs pints=0.25*( shs*(-b12-b17+b19+b14)- 3 sds*(-b7-b12+b14+b9) ) pintq=0.25*(-qfs*(b17+b18-b8-b7)+ 5 qbs*(b18+b19-b9-b8) ) c outer part (step 2h): non mixed + 0.25*(mixed) c (FOR 4th order ONLY) c c poutp=pahl*(a23-a13)-pdel*(a13-a3) c poutr=rbcl*(a15-a13)+rfgl*(a11-a13) cc pouts=0.25*( shl*(b13+b23-b11-b21)+sal*(-b13-b23+b25+b15)- cc 3 sel*(b3+b13-b1-b11) -sdl*(-b3-b13+b15+b5) ) cc poutq=0.25*(-qgl*(b21+b23-b11-b13)-qfl*(-b3-b1+b11+b13)+ cc 5 qbl*(b23+b25-b13-b15)+qcl*(-b5-b3+b13+b15) ) cc mixed derivative with equated legs c pouts=0.25*( shl*(-b12-b17+b19+b14)- c 3 sdl*(-b7-b12+b14+b9) ) c poutq=0.25*(-qfl*(b17+b18-b8-b7)+ c 5 qbl*(b18+b19-b9-b8) ) cc cc 4th order scheme cc sleft=16.*pint-pout cc cc 2nd order (POUT shouldn't be computed in this case) cc sleft=12.*pint cc c sleftp=12.*pintp sleftr=12.*pintr slefts=12.*pints sleftq=12.*pintq c sleft=sleftp+sleftr+slefts+sleftq c c DISNEW=CNST*sleft+2.*A13-AOLD DISNEW=0. IF(DNST.NE.0.) DISNEW=(CNST/DNST)*sleft+2.*A13-AOLD RETURN END C C======================================================================= C C C SUBROUTINE DIFHOM(AOLD,A3,A8,A11,A12,A13,A14,A15,A18,A23, 2 B1,B5,B7,B9,B17,B19,B21,B25, 4 DISNEW) C C Purpose: C The finite-difference scheme for a grid point in a homogeneous C medium C COMMON /PPAR/ pahs,pdes,rbcs,rfgs, 1 shs,sas,ses,sds,qgs,qfs,qbs,qcs, 2 pahl,pdel,rbcl,rfgl, 3 shl,sal,sel,sdl,qgl,qfl,qbl,qcl, 4 phom,qhom,rhom,shom,dnst COMMON /BASE/ MPAR,NSIZE,NX,NZ,NZ5,NZ6,CN,CNST C c inner part pint=phom*(A18-2.*A13+A8)+rhom*(A14-2.*A13+A12) 1 +0.25*(shom+qhom)*(b19-b17+b7-b9) c outer part (for 4th order only) c pout=phom*(A23-2.*A13+A3)+rhom*(A15-2.*A13+A11) c 1 +0.25*(shom+qhom)*(b25-b21+b1-b5) c c 4th order c sleft=16.*pint-pout c 2nd order (POUT shouldn't be computed in this case) sleft=12.*pint c c DISNEW=CNST*sleft+2.*A13-AOLD DISNEW=0. IF(DNST.NE.0.) DISNEW=(CNST/DNST)*sleft+2.*A13-AOLD RETURN END C DEBUG SUBCHK C END DEBUG C C======================================================================= C C C SUBROUTINE FLERIB(KLMAX,LMAX,LIM1,LIM2,IPT1,IPT2,IPT3, 1 ULINA,ULINB,WLINA,WLINB, 2 UX,UY,WX,WY,AH,BH,DENS) C C Purpose: C Left and right nonreflecting boundaries C DIMENSION AH(KLMAX),BH(KLMAX),DENS(KLMAX) DIMENSION UX(KLMAX),UY(KLMAX),WX(KLMAX),WY(KLMAX) DIMENSION ULINA(LMAX),ULINB(LMAX),WLINA(LMAX),WLINB(LMAX) COMMON /BASE/ MPAR,NSIZE,NX,NZ,NZ5,NZ6,CN,CNST C DO 500 L=LIM1,LIM2 C c ALFA=SQRT(AH(IPT3+L)) ALFA=0. IF(DENS(IPT3+L).NE.0.) ALFA=SQRT(AH(IPT3+L)/DENS(IPT3+L)) C2=CN*ALFA IF(L.EQ.1) C2=C2*1.414 C3=1./(1.+C2) C c BETA=SQRT(BH(IPT3+L)) BETA=0. IF(DENS(IPT3+L).NE.0.) BETA=SQRT(BH(IPT3+L)/DENS(IPT3+L)) C4=CN*BETA IF(L.EQ.1) C4=C4*1.414 C5=1./(1.+C4) C cc C6=(BETA-ALFA)/ALFA cc C7=(BETA-ALFA)/BETA C C IF(MPAR.EQ.1) THEN PAR1=UY(IPT1+L) PAR2=UX(IPT2+L) PAR3=UX(IPT1+L) PAR4=UX(IPT2+L+1) RAR1=WY(IPT1+L) RAR2=WX(IPT2+L) RAR3=WX(IPT1+L) RAR4=WX(IPT2+L+1) ELSE PAR1=UX(IPT1+L) PAR2=UY(IPT2+L) PAR3=UY(IPT1+L) PAR4=UY(IPT2+L+1) RAR1=WX(IPT1+L) RAR2=WY(IPT2+L) RAR3=WY(IPT1+L) RAR4=WY(IPT2+L+1) ENDIF C C Emerman, Stephen UAUX=C2*(PAR1-ULINA(L)+ULINB(L))- 1 (PAR1-2.*(PAR2+PAR3)+ULINB(L)+ULINA(L)) WAUX=C4*(RAR1-WLINA(L)+WLINB(L))- 1 (RAR1-2.*(RAR2+RAR3)+WLINB(L)+WLINA(L)) C C Stacey P3 C C3=1. C C5=1. C UAUX=-1.*C2*(PAR2-PAR3)+C6*C2*(RAR4-RAR2)+PAR2 C WAUX=-1.*C4*(RAR2-RAR3)+C7*C4*(PAR4-PAR2)+RAR2 C C Clayton, Engquist C C3=1. C C5=1. C UAUX=-1.*C2*(PAR2-PAR3)+PAR2 C WAUX=-1.*C4*(RAR2-RAR3)+RAR2 C C Reynolds C Attention: CHANGE ALSO CALL FLERIB AND IZ4 IN SAVER !!!!!!!!!!!! C C3=1. C C5=1. C UAUX=-1.*C2*(PAR2-PAR3-ULINB(L)+ULINA(L))+(PAR2+PAR3-ULINB(L)) C WAUX=-1.*C4*(RAR2-RAR3-WLINB(L)+WLINA(L))+(RAR2+RAR3-WLINB(L)) C IF(MPAR.EQ.1) THEN UY(IPT2+L)=C3*UAUX WY(IPT2+L)=C5*WAUX ELSE UX(IPT2+L)=C3*UAUX WX(IPT2+L)=C5*WAUX ENDIF 500 CONTINUE RETURN END C DEBUG SUBCHK C END DEBUG C C======================================================================= C C C SUBROUTINE FBOB(KLMAX,KMAX,KIM1,KIM2,IPT1,IPT2, 1 ULINA,ULINB,WLINA,WLINB, 2 UX,UY,WX,WY,AV,BV,DENS) C C Purpose: C Bottom nonreflecting boundaries C DIMENSION AV(KLMAX),BV(KLMAX),DENS(KLMAX) DIMENSION UX(KLMAX),UY(KLMAX),WX(KLMAX),WY(KLMAX) DIMENSION ULINA(KMAX),ULINB(KMAX),WLINA(KMAX),WLINB(KMAX) COMMON /BASE/ MPAR,NSIZE,NX,NZ,NZ5,NZ6,CN,CNST C DO 500 K=KIM1,KIM2 IPT=IPT1+(K-1)*NZ IPTV=K*(NZ-1)-IPT2 NZ99=K*NZ-IPT2 C c C2=CN*SQRT(BV(IPTV)) C2=0. IF(DENS(NZ99).NE.0.) C2=CN*SQRT(BV(IPTV)/DENS(NZ99)) C3=1./(1.+C2) c C4=CN*SQRT(AV(IPTV)) C4=0. IF(DENS(NZ99).NE.0.) C4=CN*SQRT(AV(IPTV)/DENS(NZ99)) C5=1./(1.+C4) C IF(MPAR.EQ.1) THEN PAR1=UY(IPT) PAR2=UX(IPT+1) PAR3=UX(IPT) RAR1=WY(IPT) RAR2=WX(IPT+1) RAR3=WX(IPT) ELSE PAR1=UX(IPT) PAR2=UY(IPT+1) PAR3=UY(IPT) RAR1=WX(IPT) RAR2=WY(IPT+1) RAR3=WY(IPT) ENDIF UAUX=C2*(PAR1-ULINA(K)+ULINB(K))- 1 (PAR1-2.*(PAR2+PAR3)+ULINB(K)+ULINA(K)) WAUX=C4*(RAR1-WLINA(K)+WLINB(K))- 1 (RAR1-2.*(RAR2+RAR3)+WLINB(K)+WLINA(K)) IF(MPAR.EQ.1) THEN UY(IPT+1)=C3*UAUX WY(IPT+1)=C5*WAUX ELSE UX(IPT+1)=C3*UAUX WX(IPT+1)=C5*WAUX ENDIF 500 CONTINUE RETURN END C DEBUG SUBCHK C END DEBUG C C======================================================================= C C C SUBROUTINE FCORNB(KLMAX,NIND,NDIF1,NDIF2,NLFT, 1 UX,UY,WX,WY,AH,BH,DENS) C C Purpose: C Left and right bottom corners of the nonreflecting boundaries C The 'corner' means 3 points: x , or x C x x x x C DIMENSION AH(KLMAX),BH(KLMAX),DENS(KLMAX) DIMENSION UX(KLMAX),UY(KLMAX),WX(KLMAX),WY(KLMAX) COMMON /BASE/ MPAR,NSIZE,NX,NZ,NZ5,NZ6,CN,CNST C c ALFA=SQRT(AH(NIND-NDIF1)) c BETA=SQRT(BH(NIND-NDIF1)) ALFA=0. NZ99=NIND-NDIF1 IF(DENS(NZ99).NE.0.) ALFA=SQRT(AH(NIND-NDIF1)/DENS(NZ99)) BETA=0. IF(DENS(NZ99).NE.0.) BETA=SQRT(BH(NIND-NDIF1)/DENS(NZ99)) C4=ALFA+BETA IF(NLFT.EQ.1) THEN C5=BETA-ALFA ELSE C5=ALFA-BETA ENDIF C6=CN/1.4142135/2. C************************** DO 500 I=1,3 IF(I.EQ.1) IND=NIND-1 IF(I.EQ.2) IND=NIND-NDIF2 IF(I.EQ.3) IND=NIND C************************** C IND=NIND C************************** IF(MPAR.EQ.1) THEN UAUX=C4*(UX(IND-1)-2.*UX(IND)+UX(IND-NDIF2))+ 1 C5*(WX(IND-1)-2.*WX(IND)+WX(IND-NDIF2)) WAUX=C4*(WX(IND-1)-2.*WX(IND)+WX(IND-NDIF2))+ 1 C5*(UX(IND-1)-2.*UX(IND)+UX(IND-NDIF2)) UY(IND)=C6*UAUX+UX(IND) WY(IND)=C6*WAUX+WX(IND) ELSE UAUX=C4*(UY(IND-1)-2.*UY(IND)+UY(IND-NDIF2))+ 1 C5*(WY(IND-1)-2.*WY(IND)+WY(IND-NDIF2)) WAUX=C4*(WY(IND-1)-2.*WY(IND)+WY(IND-NDIF2))+ 1 C5*(UY(IND-1)-2.*UY(IND)+UY(IND-NDIF2)) UX(IND)=C6*UAUX+UY(IND) WX(IND)=C6*WAUX+WY(IND) ENDIF C********************** 500 CONTINUE C*********************** RETURN END C DEBUG SUBCHK C END DEBUG C C======================================================================= C C C SUBROUTINE SAVER(KMAX,LMAX,KLMAX,A,B, 1 UK1,UK2,UK3,UNX2,UNX1,UNX,UNZ2,UNZ1,UNZ, 2 WK1,WK2,WK3,WNX2,WNX1,WNX,WNZ2,WNZ1,WNZ) C C Purpose: C Saving 'very old' values (from time level M-1) at boundary C lines for subsequent computations at nonreflecting boundaries C DIMENSION UK1(LMAX),UK2(LMAX),UK3(LMAX) DIMENSION UNX2(LMAX),UNX1(LMAX),UNX(LMAX) DIMENSION UNZ2(KMAX),UNZ1(KMAX),UNZ(KMAX) DIMENSION WK1(LMAX),WK2(LMAX),WK3(LMAX) DIMENSION WNX2(LMAX),WNX1(LMAX),WNX(LMAX) DIMENSION WNZ2(KMAX),WNZ1(KMAX),WNZ(KMAX) DIMENSION A(KLMAX),B(KLMAX) COMMON /BASE/ MPAR,NSIZE,NX,NZ,NZ5,NZ6,CN,CNST C DO 100 I=1,NZ UK1(I)=A(I) WK1(I)=B(I) UK2(I)=A(NZ+I) WK2(I)=B(NZ+I) IZ1=2*NZ+I UK3(I)=A(IZ1) WK3(I)=B(IZ1) IZ2=NSIZE-3*NZ+I UNX2(I)=A(IZ2) WNX2(I)=B(IZ2) IZ3=IZ2+NZ UNX1(I)=A(IZ3) WNX1(I)=B(IZ3) C For Emerman-Stephen: IZ4=IZ3+NZ C For Reynolds: C IZ4=NSIZE-4*NZ+I UNX(I)=A(IZ4) WNX(I)=B(IZ4) 100 CONTINUE DO 200 I=1,NX IZ5=NZ-2+(I-1)*NZ UNZ2(I)=A(IZ5) WNZ2(I)=B(IZ5) IZ6=IZ5+1 UNZ1(I)=A(IZ6) WNZ1(I)=B(IZ6) IZ7=IZ6+1 UNZ(I)=A(IZ7) WNZ(I)=B(IZ7) 200 CONTINUE RETURN END C DEBUG SUBCHK C END DEBUG C C======================================================================= C C C FUNCTION TIMFUN(T,TAU,TP) C C Purpose: C The excitation signal, which is the half-th integral of a seismic C force or the half-th derivative of a seismic moment C TIMFUN=0. IF(T.LT.TAU.OR.T.GT.TAU+2.*TP) RETURN PI=3.141593 C C Type of the signal CALL RSEP3I('KSIG',KSIG,3) C Reference frequency SIGF=1/(2*TP) C Reference time SIGT=TP C Reference amplitude CALL RSEP3R('SIGA',SIGA,1.) C Reference width in reference half-periods: CALL RSEP3R('SIGW',SIGW,1.) C Phase CALL RSEP3R('SIGPH',SIGPH,1.) C GO TO (1,1,10,1,30,1,1,1) KSIG+2 1 CONTINUE C FD2D-06 CALL ERROR('FD2D-06: Only signals KSIG=1,3 allowed') C C Gabor: C [Hron: frequency SIGF=1, width SIGW=4, phase SIGPH=-pi/2] 10 CONTINUE FNUL=1. ARA=2.*PI*FNUL*(T-TAU-TP) ARAG=ARA/SIGW ARAG=-1.*ARAG*ARAG TIMFUN=SIGA*SIN(ARA)*EXP(ARAG) RETURN C C [Jih: frequency SIGF=0.45, width SIGW=0.05, phase SIGPH=0] 11 CONTINUE FNUL=0.45 GAMA=0.05 ARA=2.*PI*FNUL*(T-TAU-TP) ARAG=ARA/GAMA ARAG=-1.*ARAG*ARAG TIMFUN=COS(ARA)*EXP(ARAG) RETURN C C [Bard, model B: C frequency SIGF=2.416666, width SIGW=7.5, phase SIGPH=-pi] 12 CONTINUE FNUL=2.416666 GAMA=7.5 ARA=2.*PI*FNUL*(T-TAU-TP) ARAG=ARA/GAMA ARAG=-1.*ARAG*ARAG TIMFUN=-1.*COS(ARA)*EXP(ARAG) RETURN C C Hermite-Gaussian (Ricker) [Gauss: order SIGW=1]: 21 CONTINUE GAMA=40. ARB=T-TAU-TP ARD=-2.*GAMA*ARB ARE=-1.*GAMA*ARB*ARB TIMFUN=ARD*EXP(ARE) RETURN C C Hermite-Gaussian (Ricker) [Ricker: order SIGW=2]: 22 CONTINUE PI2=PI*PI A=PI2*(T-TAU-TP)*(T-TAU-TP)/TP/TP XMAX=TP*(PI2-0.5)/SQRT(PI) TIMFUN=(PI2-A)*EXP(-1.*A)/XMAX RETURN C C Hermite-Gaussian (Ricker) [Bard, models E,F: SIGW=2] 23 CONTINUE c FNUL=0.862070 for test E c FNUL=0.714286 for test F ARG=PI*FNUL*(T-TAU-TP) ARG=ARG*ARG TIMFUN=-0.707107*(ARG-0.5)*EXP(-1.*ARG) RETURN C C Hermite-Gaussian (Ricker) [Kawase: SIGW=2] 24 CONTINUE FC=0.25 ARG=PI*PI*FC*FC*(T-TAU-TP)*(T-TAU-TP) TIMFUN=(2.*ARG-1.)*EXP(-1.*ARG) RETURN C C Kuepper (Mueller) (duration=2TP): 30 CONTINUE C Reference amplitude is the maximum amplitude AMPL=SIGA*4./SQRT(27.) DPI=PI+PI ARGMN1=DPI*(T-TAU)/(2.*TP) ARGMN2=2.*ARGMN1 TIMFUN=AMPL*(SIN(ARGMN1)-0.5*SIN(ARGMN2)) RETURN C c Priolo (max. frequency 22 Hz) fmax=22. eta=0.5 eps=1. s1=-fmax*fmax*(t-tp)*(t-tp)*eta s2=eps*pi*fmax*(t-tp) cc sign. changed in contrast to Priolo definition due to cc opposite sign. in my source (-exp changed to +exp) TIMFUN=+exp(s1)*fmax* 1 (2.*fmax*eta*cos(s2)*(t-tp)+eps*pi*sin(s2)) RETURN c C Aboudi TIMFUN=0. X=T-TAU IF(X.LT.0..OR.X.GT.2.*TP) RETURN A=ABS(X-TP) B=2./TP/TP IF(A.GE.0..AND.A.LE.TP/2.) GOTO 100 IF(A.GE.TP/2..AND.A.LE.TP) GOTO 200 100 TIMFUN=1.-B*A*A RETURN 200 C=A-TP/2. TIMFUN=1.-B*(A*A-2.*C*C) RETURN C RETURN END C C======================================================================= C INCLUDE 'error.for' C error.for INCLUDE 'fdaux.for' C fdaux.for INCLUDE 'gmtra.for' C gmtra.for INCLUDE 'sep.for' C sep.for INCLUDE 'gse.for' C gse.for INCLUDE 'forms.for' C forms.for INCLUDE 'length.for' C length.for C C======================================================================= Cfdaux.for 0100666 0000766 0000766 00000043372 06735567274 012306 0 ustar klimes klimes C
C Subroutine file 'fdaux.for' to mediate input and output of finite C diference programs C C Version: 5.40 C Date: 1999, June 28 C Coded by Ludek Klimes C C....................................................................... C C This Fortran77 file consists of the following external procedures: C FDIN... Subroutine to read the effective elastic parameters for C 2-D finite difference programs. C FDIN C FDIN1...Subroutine to read the values of a single effective C elastic parameter for 2-D finite difference programs. C FDIN1 C FDOUT1..Subroutine to write the discrete values of a single C quantity for 2-D finite difference programs. C FDOUT1 C FDGRID..Subroutine to read grid dimensions for 2-D finite C differences. C FDGRID C FDREC...Subroutine to read receiver and source positions for 2-D C finite differences. C FDREC C C======================================================================= C C C SUBROUTINE FDIN(LU,NX,NZ,AHL,BHL,CHL,AVL,BVL,CVL, * AHS,BHS,CHS,AVS,BVS,CVS,DEN,QFC,WORK) C INTEGER LU,NX,NZ REAL AHL(NZ ,NX-1),BHL(NZ ,NX-1),CHL(NZ ,NX-1) REAL AVL(NZ-1,NX ),BVL(NZ-1,NX ),CVL(NZ-1,NX ) REAL AHS(NZ-1,NX-1),BHS(NZ-1,NX-1),CHS(NZ-1,NX-1) REAL AVS(NZ-1,NX-1),BVS(NZ-1,NX-1),CVS(NZ-1,NX-1) REAL DEN(NZ ,NX ),QFC(NZ ,NX ),WORK(NZ ,NX ) C C Subroutine FDIN reads the effective elastic parameters for 2-D finite C difference programs. C C Input: C LU... Logical unit number to be used to read all files. C NX... Number of grid points in the horizontal direction. C NZ... Number of grid points in the vertical direction. C AHL,BHL,CHL,AVL,BVL,CVL,AHS,BHS,CHS,AVS,BVS,CVS,DEN,QFC,WORK... C Arrays of the dimensions as declared above or greater. C WORK is a temporary working array, its content will be C destroyed. C C Output: C AHL,BHL,CHL,AVL,BVL,CVL,AHS,BHS,CHS,AVS,BVS,CVS,DEN,QFC... Values C read from formatted or unformatted files given by SEP C parameters AHL,BHL,CHL,AVL,BVL,CVL,AHS,BHS,CHS,AVS,BVS, C CVS,DEN,QFC, respectively. C WORK... Undefined. C C Date: 1998, November 11 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C INTEGER NEWPAR C NEWPAR..NEWPAR=0: Average parameters along the legs are read. C NEWPAR=1: Average parameters over the grid facets are read C and the parameters corresponding to the legs are derived C from them. C C....................................................................... C CALL RSEP3I('NEWPAR',NEWPAR,0) IF(NEWPAR.EQ.0) THEN CALL FDIN1(LU,NX-1,NZ-1,'A13',AHS,WORK,0,AHS) CALL FDIN1(LU,NX-1,NZ-1,'B13',BHS,WORK,0,BHS) CALL FDIN1(LU,NX-1,NZ-1,'C13',CHS,WORK,0,CHS) CALL FDIN1(LU,NX-1,NZ-1,'A31',AVS,WORK,0,AHS) CALL FDIN1(LU,NX-1,NZ-1,'B31',BVS,WORK,0,BHS) CALL FDIN1(LU,NX-1,NZ-1,'C31',CVS,WORK,0,CHS) CALL FDIN1(LU,NX-1,NZ ,'A11',AHL,WORK,0,AHS) CALL FDIN1(LU,NX-1,NZ ,'B11',BHL,WORK,0,BHS) CALL FDIN1(LU,NX-1,NZ ,'C11',CHL,WORK,0,CHS) CALL FDIN1(LU,NX ,NZ-1,'A33',AVL,WORK,0,AVS) CALL FDIN1(LU,NX ,NZ-1,'B33',BVL,WORK,0,BVS) CALL FDIN1(LU,NX ,NZ-1,'C33',CVL,WORK,0,CVS) ELSE CALL FDIN1(LU,NX-1,NZ-1,'AA13',AHS,WORK,0,AHS) CALL FDIN1(LU,NX-1,NZ-1,'BB13',BHS,WORK,0,BHS) CALL FDIN1(LU,NX-1,NZ-1,'CC13',CHS,WORK,0,CHS) CALL FDIN1(LU,NX-1,NZ-1,' ',AVS,WORK,4,AHS) CALL FDIN1(LU,NX-1,NZ-1,' ',BVS,WORK,4,BHS) CALL FDIN1(LU,NX-1,NZ-1,' ',CVS,WORK,4,CHS) CALL FDIN1(LU,NX-1,NZ ,' ',AHL,WORK,2,AHS) CALL FDIN1(LU,NX-1,NZ ,' ',BHL,WORK,2,BHS) CALL FDIN1(LU,NX-1,NZ ,' ',CHL,WORK,2,CHS) CALL FDIN1(LU,NX ,NZ-1,' ',AVL,WORK,1,AVS) CALL FDIN1(LU,NX ,NZ-1,' ',BVL,WORK,1,BVS) CALL FDIN1(LU,NX ,NZ-1,' ',CVL,WORK,1,CVS) END IF CALL FDIN1(LU,NX ,NZ ,'DEN',DEN,WORK,0,DEN) CALL FDIN1(LU,NX ,NZ ,'QFC',QFC,WORK,0,QFC) RETURN END C C======================================================================= C C C SUBROUTINE FDIN1(LU,N1,N2,NAME,ARRAY,WORK,KDEF,DEF) C INTEGER LU,N1,N2,KDEF CHARACTER*(*) NAME REAL ARRAY(N2*N1),WORK(N1*N2),DEF(*) C C Subroutine FDIN1 reads the values of a single effective elastic C parameter for 2-D finite difference programs. C C Input: C LU... Logical unit number to be used for input. C N1... Number of values in the horizontal direction. C N2... Number of values in the vertical direction. C NAME... String containing the name of the parameter specifying the C name of the file containing the grid values to be read. C Except for its case, it should match the parameter name in C a previously read input SEP parameter file. If parameter C NAME is not defined in the input SEP parameter file, no C input is performed by this subroutine and the output C arrays are filled with zeros. C KDEF... Kind of default values: C KDEF=0: Zeros. Array DEF is not used. C KDEF=1: Array DEF should have dimensions DEF(N2,N1-1). C Then the default is ARRAY(I2,1)=DEF(I2,1), C ARRAY(I2,N1)=DEF(I2,N1-1), and for 1.LT.I1.LT.N1 C ARRAY(I2,I1)=(DEF(I2,I1-1)+DEF(I2,I1))/2. C KDEF=2: Array DEF should have dimensions DEF(N2-1,N1). C Then the default is ARRAY(1,I1)=0.5*DEF(1,I1), C ARRAY(N2,I1)=DEF(N2-1,I1), and for 1.LT.I2.LT.N2 C ARRAY(I2,I1)=(DEF(I2-1,I1)+DEF(I2,I1))/2. C KDEF=4: Array DEF should have dimensions DEF(N2,N1). C Then the default is ARRAY(I2,I1)=DEF(I2,I1). C DEF... Array to specify default values. C C Output: C ARRAY...Array of N2*N1 values. The values are read from the file C which name is given by parameter NAME in SEP input data C read previously to the invocation of this subroutine. C Since the inner loop in the input file is over the C horizontal grid lines (N1 values), whereas the inner loop C in 2-D FD programs is over the vertical grid lines C (N2 values), N1*N2 matrix WORK read from the file is C transposed to yield N2*N1 output array ARRAY. C WORK... Temporary storage array containing the grid values in the C same order as in the input file. C C Date: 1998, November 11 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: CHARACTER*80 FILE CHARACTER*11 FORM INTEGER I1,I2,I,J C C....................................................................... C C Name and form of the input data file: IF(NAME.EQ.' ') THEN FILE=' ' ELSE CALL RSEP3T(NAME,FILE,' ') CALL RSEP3T('FORM',FORM,'FORMATTED') END IF C C Check whether the filename is specified in the input SEP data: IF(FILE.EQ.' ') THEN C C No input file specified: I=1 J=1 IF(KDEF.EQ.1) THEN C Averaging along the X1 axis: DO 12 I2=1,N2 ARRAY(I)=DEF(J) DO 11 I1=2,N1-1 I=I+N2 J=J+N2 ARRAY(I)=0.5*(DEF(J-N2)+DEF(J)) 11 CONTINUE I=I+N2 ARRAY(I)=DEF(J) I=I-N2*(N1-1)+1 J=J-N2*(N1-2)+1 12 CONTINUE ELSE IF(KDEF.EQ.2) THEN C Averaging along the X2 axis: DO 22 I1=1,N1 ARRAY(I)=0.5*DEF(J) DO 21 I2=2,N2-1 I=I+1 J=J+1 ARRAY(I)=0.5*(DEF(J-1)+DEF(J)) 21 CONTINUE I=I+1 ARRAY(I)=DEF(J) I=I+1 J=J+1 22 CONTINUE ELSE IF(KDEF.EQ.4) THEN C Copying: DO 40 I=1,N2*N1 ARRAY(I)=DEF(I) 40 CONTINUE ELSE C Zeros: DO 90 I=1,N1*N2 ARRAY(I)=0. 90 CONTINUE END IF C C Transposing the array: * CALL GMTRA(ARRAY,WORK,N2,N1) C ELSE C C Reading the array (undefined values are replaced by zeros): CALL RARRAY(LU,FILE,FORM,.TRUE.,0.,N1*N2,WORK) C C Transposing the array: CALL GMTRA(WORK,ARRAY,N1,N2) C END IF RETURN END C C======================================================================= C C C SUBROUTINE FDOUT1(LU,N1,N2,SIGNUM,NAME,ARRAY,WORK) C INTEGER LU,N1,N2 CHARACTER*(*) NAME REAL SIGNUM,ARRAY(N2*N1),WORK(N1*N2) C C Subroutine FDOUT1 to write the discrete values of a single quantity C for 2-D finite difference programs. C C Input: C LU... Logical unit number to be used for output. C N1... Number of values in the horizontal direction. C N2... Number of values in the vertical direction. C SIGNUM..If SIGNUM is negative, the array values will be written C with opposite signs. C NAME... String containing the name of the parameter specifying the C name of the output file to be written. C Except for its case, it should match the parameter name in C a previously read input SEP parameter file. If parameter C NAME is not defined in the input SEP parameter file, no C output is performed. C ARRAY...Array of N2*N1 values. The values will be written to the C file which name is given by parameter NAME in SEP input C data read previously to the invocation of this subroutine. C Since the inner loop in the input file is over the C horizontal grid lines (N1 values), whereas the inner loop C in 2-D FD programs is over the vertical grid lines C (N2 values), N2*N1 matrix ARRAY is transposed to yield C N1*N2 array WORK written to the output file. C WORK... Temporary working array, its content will be destroyed. C C Output: C WORK... Undefined. C C Date: 1999, April 6 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: CHARACTER*80 FILE CHARACTER*11 FORM INTEGER I C C....................................................................... C C Name of the output data file: IF(NAME.EQ.' ') THEN C The output file is already open: FILE=' ' ELSE CALL RSEP3T(NAME,FILE,' ') C Check whether the filename is specified in the input SEP data: IF(FILE.EQ.' ') THEN RETURN END IF END IF C C Form of the output data file: CALL RSEP3T('FORM',FORM,'FORMATTED') C C Transposing the array: CALL GMTRA(ARRAY,WORK,N2,N1) C IF(SIGNUM.LT.0) THEN DO 10 I=1,N1*N2 WORK(I)=-WORK(I) 10 CONTINUE END IF C C Writing the array: CALL WARRAY(LU,FILE,FORM,.FALSE.,0.,.FALSE.,0.,N1*N2,WORK) C RETURN END C C======================================================================= C C C SUBROUTINE FDGRID(KMAX,LMAX,KLMAX,NX,NZ,DX) C INTEGER KMAX,LMAX,KLMAX,NX,NZ REAL DX C C Subroutine to read grid dimensions for 2-D finite differences. C C Input: C KMAX... Maximum number of vertical gridlines. C LMAX... Maximum number of horizontal gridlines. C KLMAX...Maximum number of gridpoints. C C Output: C NX... Number of vertical gridlines. C NZ... Number of horizontal gridlines. C DX... Absolute value of the grid interval. The absolute values C of the horizontal and vertical grid intervals must be C equal. C C Date: 1998, September 11 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER N1,N2,N3 REAL D1,D2,D3,DZ C C....................................................................... C C Recalling the data specifying grid dimensions C (arguments: Name of parameter in input data, Variable, Default): CALL RSEP3I('N1',N1,1) CALL RSEP3I('N2',N2,1) CALL RSEP3I('N3',N3,1) CALL RSEP3R('D1',D1,1.) CALL RSEP3R('D2',D2,1.) CALL RSEP3R('D3',D3,1.) IF (N1.LE.1) THEN NX=N2 NZ=N3 DX=ABS(D2) DZ=ABS(D3) ELSE IF (N2.LE.1) THEN NX=N1 NZ=N3 DX=ABS(D1) DZ=ABS(D3) ELSE IF (N3.LE.1) THEN NX=N1 NZ=N2 DX=ABS(D1) DZ=ABS(D2) ELSE C FDAUX-01 CALL ERROR('FDAUX-01: 3-D grid specified for 2-D FD') C Exactly one of input parameters N1, N2 ,N3 must equal 1 for 2-D C finite differences. END IF IF (NX.LE.1.OR.NZ.LE.1) THEN C FDAUX-02 CALL ERROR('FDAUX-02: 1-D grid specified for 2-D FD') C Exactly one of input parameters N1, N2 ,N3 must equal 1 for 2-D C finite differences. END IF IF (NX.GT.KMAX) THEN C FDAUX-03 CALL ERROR('FDAUX-03: Small dimension KMAX') C Number NX of vertical gridlines exceeds maximum number KMAX of C vertical gridlines. END IF IF (NZ.GT.LMAX) THEN C FDAUX-04 CALL ERROR('FDAUX-04: Small dimension LMAX') C Number NZ of horizontal gridlines exceeds maximum number LMAX of C horizontal gridlines. END IF IF (NX*NZ.GT.KLMAX) THEN C FDAUX-05 CALL ERROR('FDAUX-05: Small dimension KLMAX') C Number N1*N2*N3 of gridpoints exceeds maximum number KLMAX of C gridpoints. END IF IF (DX.NE.DZ) THEN C FDAUX-06 CALL ERROR('FDAUX-06: Different grid intervals') C Horizontal and vertical grid intervals must be equal for C 2-D finite-difference programs. END IF C RETURN END C C======================================================================= C C C SUBROUTINE FDREC(LU,MAXREC,NREC,REC,KRE,LRE,KSOUR,LSOUR) C INTEGER LU,MAXREC,NREC,KRE(MAXREC),LRE(MAXREC),KSOUR,LSOUR CHARACTER*(*) REC(MAXREC) C C Subroutine to read receiver and source positions for 2-D finite C differences. C C Input: C LU... Logical unit number to be used for input. C MAXREC..Maximum number of receivers. C C Output: C NREC... Number of receivers. C KREC,LREC... Horizontal and vertical indices of receiver C gridpoints. C KSOUR,LSOUR... Horizontal and vertical indices of the source C gridpoint. C C Date: 1998, October 7 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C REAL UNDEF PARAMETER (UNDEF=-999999.) C C Auxiliary storage locations: CHARACTER*80 FREC CHARACTER*6 TEXT INTEGER N1,N3 REAL D1,D3,O1,O3,X1,X2,X3 C C....................................................................... C CALL RSEP3I('N1' ,N1 ,1 ) CALL RSEP3I('N3' ,N3 ,1 ) CALL RSEP3R('D1' ,D1 ,1. ) CALL RSEP3R('D3' ,D3 ,1. ) CALL RSEP3R('O1' ,O1 ,0. ) CALL RSEP3R('O3' ,O3 ,0. ) C NREC=0 CALL RSEP3T('REC',FREC,' ') IF(FREC.NE.' ') THEN OPEN(LU,FILE=FREC,STATUS='OLD') READ(LU,*) (TEXT,I=1,20) 11 CONTINUE X1=UNDEF X2=UNDEF X3=UNDEF READ(LU,*,END=12) TEXT,X1,X2,X3 IF(X1.EQ.UNDEF.AND.X2.EQ.UNDEF.AND.X3.EQ.UNDEF) THEN GO TO 12 END IF IF(NREC.GE.MAXREC) THEN C FDAUX-07 CALL ERROR('FDAUX-07: Small dimension MAXREC') C Dimension MRAM must be at least C 19*N1*N3+10*N1+12*N3+(2*NTFD+5)*MAX(1,NREC), C where N1, N3 and NTFD are parameters in the main input file C and NREC is the number of receivers specified in the file C which name is given by parameter REC. END IF NREC=NREC+1 REC(NREC)=TEXT IF(X1.EQ.UNDEF) X1=0. IF(X2.EQ.UNDEF) X2=0. IF(X3.EQ.UNDEF) X3=0. KRE(NREC)=NINT((X1-O1)/D1)+1 LRE(NREC)=NINT((X3-O3)/D3)+1 IF(KRE(NREC).LT.1.OR.KRE(NREC).GT.N1.OR. * LRE(NREC).LT.1.OR.LRE(NREC).GT.N3) THEN C FDAUX-08 CALL ERROR('FDAUX-08: Receiver point outside the FD grid') END IF GO TO 11 12 CONTINUE CLOSE(LU) END IF C CALL RSEP3T('SRC',FREC,' ') IF(FREC.NE.' ') THEN OPEN(LU,FILE=FREC,STATUS='OLD') READ(LU,*) (TEXT,I=1,20) X1=0. X2=0. X3=0. READ(LU,*) TEXT,X1,X2,X3 CLOSE(LU) KSOUR=NINT((X1-O1)/D1)+1 LSOUR=NINT((X3-O3)/D3)+1 IF(KSOUR.LT.1.OR.KSOUR.GT.N1.OR. * LSOUR.LT.1.OR.LSOUR.GT.N3) THEN C FDAUX-09 CALL ERROR('FDAUX-09: Source point outside the FD grid') END IF END IF RETURN END C C======================================================================= Cfdcolor.bat 0100666 0000766 0000766 00000000133 06507336746 012566 0 ustar klimes klimes GS386 -sDEVICE=pcx24b -r72 -g480x480 -dNOPAUSE -sOutputFile=fdcolor.pcx fdcolor.ps -c quit fdcolor.ps 0100666 0000766 0000766 00000003664 06621450556 012447 0 ustar klimes klimes %! %%BoundingBox: 14 14 468 468 % For bitmap 480*480 %----------------------------------------------------------- /NSAT 1 def /NHUE 250 def % First version of the figures /NSAT 60 def /NHUE 250 def % Second version of the figures /NSAT 7 def /NHUE 36 def %----------------------------------------------------------- /cm {2.54 div 72 mul} def /X0 8.5 cm def /Y0 8.5 cm def % For bitmap 480*480 %X0 10.5 cm def /Y0 17.5 cm def % For hardcopy /R 6.0 cm def % Radius of the colour circle /A 8.0 cm def % Half-axes /H 1 NHUE div def % Hue increment /C0 -5 24 div def % Colour corresponding to the divergence /DOTRADIUS 0.150 cm def /FONT {27 mul} def /CSIZE {30 mul} def %----------------------------------------------------------- /DOT {DOTRADIUS 0 360 arc fill} def /WRITE {1.0 FONT selectfont show} def /COLOR {360 div dup 1 1 sethsbcolor C0 sub 360 mul dup cos R 0.7 CSIZE add mul X0 add exch sin R 18 add mul Y0 add moveto -0.33 CSIZE -0.37 CSIZE rmoveto /Helvetica 1.0 CSIZE selectfont show} def % Plotting colour circle: 1 1 NSAT div neg 0 {dup /SAT exch def R mul /RR exch def % H 2 div neg dup C0 sub 360 mul exch H neg dup C0 sub 360 mul exch H add H 1 { dup SAT 1 sethsbcolor C0 sub 360 mul dup X0 Y0 moveto X0 Y0 RR 6 -2 roll arc closepath fill } for } for % Plotting dots: 0 0 0 setrgbcolor 15 15 360 { dup cos R mul X0 add exch sin R mul Y0 add DOT } for % Coordinate axes: 0 0 0 setrgbcolor X0 A sub Y0 moveto X0 A add Y0 lineto currentpoint stroke moveto -1.9 FONT -0.9 FONT rmoveto (div) /Helvetica WRITE (u) /Helvetica-BoldOblique WRITE X0 Y0 A sub moveto X0 Y0 A add lineto currentpoint stroke moveto -4.5 FONT -0.9 FONT rmoveto (SAMP ) /Courier 0.8 FONT selectfont show (rot) /Helvetica WRITE (u) /Helvetica-BoldOblique WRITE % Abbreviations of colour names: (R) 0 COLOR (Y) 60 COLOR (G) 120 COLOR (C) 180 COLOR (B) 240 COLOR (M) 300 COLOR showpage %%EOF fderr.htm 0100666 0000766 0000766 00000002564 06613251360 012256 0 ustar klimes klimes
# # Perl script to convert PCX snapshots into GIF (just for MS-DOS!) # # Version: 5.20 # Date: 1998, October 8 # # Usage: # perl fdgifdos.pl SEP # Parameters: # SEP... Name of the history file. # Parameters read from the history file: # SNAPS=string... Mask for the names of PostScript snapshot files. # Default: SNAPS='fd00*.ps' # MOVIE=string... Names of the animated GIF file. # Default: MOVIE='fd.gif' # Include files required: # sep.pl # Programs required: # display # gifmerge # # ====================================================================== # Main program 'fdgifdos.pl': # ~~~~~~~~~~~~~~~~~~~~~~~~~~~ require 'sep.pl'; $SEP="$ARGV[0]"; # &RSEP1($SEP); &RSEP3('SNAPS',$SNAPS,'fd00*.ps'); &RSEP3('MOVIE',$MOVIE,'fd.gif'); # unlink(<${SNAPS}>); # $SNAPS=~s/.ps//; open(LU,"|display -b gif $SNAPS.pcx ."); close(LU) || die "Error"; # # Merging GIF files into one animated GIF file: print "creating '$MOVIE'..."; open(LU,"|gifmerge -20 -l0 $SNAPS.gif>$MOVIE"); close(LU) || die "Error"; print "\nFDGIFDOS: Done.\n"; # ====================================================================== 1; #fdgiflux.pl 0100666 0000766 0000766 00000002472 06617213412 012605 0 ustar klimes klimes #!perl #
# # Perl script to convert PCX snapshots into GIF (just for Linux!) # # Version: 5.20 # Date: 1998, October 30 # # Usage: # perl fdgiflux.pl SEP # Parameters: # SEP... Name of the history file. # Parameters read from the history file: # SNAPS=string... Mask for the names of PostScript snapshot files. # Default: SNAPS='fd00*.ps' # MOVIE=string... Names of the animated GIF file. # Default: MOVIE='fd.gif' # Include files required: # sep.pl # Programs required: # convert # gifmerge # # ====================================================================== # Main program 'fdgiflux.pl': # ~~~~~~~~~~~~~~~~~~~~~~~~~~~ require 'sep.pl'; $SEP="$ARGV[0]"; # &RSEP1($SEP); &RSEP3('SNAPS',$SNAPS,'fd00*.ps'); &RSEP3('MOVIE',$MOVIE,'fd.gif'); # foreach $SNAP (<${SNAPS}>) { $SNAP=~s/.ps//; open(LU,"|convert -verbose $SNAP.pcx $SNAP.gif"); close(LU) || die "Error"; } # # Merging GIF files into one animated GIF file: $SNAPS=~s/.ps//; open(LU,"|gifmerge -20 -l0 $SNAPS.gif>$MOVIE"); close(LU) || die "Error"; print "\nFDGIFLUX: Done.\n"; # ====================================================================== 1; #fd.h 0100666 0000766 0000766 00000010263 07114363234 011200 0 ustar klimes klimes # SEP history file with demo data for 2-D FD # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # Input files required #chk.pl: "fd/" "modfd.dat" #chk.pl: "fd/" "fdmod.dat" #chk.pl: "fd/" "fdsrc.dat" #chk.pl: "fd/" "fdrec.dat" #chk.pl: "fd/" "snapdiv.cal" #chk.pl: "fd/" "snaprot.cal" #chk.pl: "fd/" "snap.cal" #chk.pl: "fd/" "addrec.cal" #chk.pl: "fd/" "fdmax.pl" #chk.pl: "forms/" "sep.pl" # Data specifying the model by means of the MODEL package MODEL='modfd.dat' # Alternative (old) model specification for program 'fdmod.for' FDMOD='fdmod.dat' # (this backward compatibility is included just for test purposes) # Grid dimensions, FD time steps, snapshots N1=111 N3=41 NTFD=1601 N4=321 D1=1.000 D3=1.000 DTFD=0.000860 D4=0.004300 O1=0.000 O3=0.000 O4=0.000 # Type and position of the point source SRC='fdsrc.dat' # file with the position of the point source KSIG=3 # Kuepper signal (with 2 local maxima) SIGF=3.333333 # ref. frequency of the Kuepper signal SIGA=1.1241E10 # SIGA=SQRT(27)*DEN*DX**3/DT**2, DEN=1600 # Files with receiver positions and synthetic seismograms REC='fdrec.dat' SS='fd.gse' # Form of files with effective elastic parameters and snapshots FORM='formatted' # Files with snapshots SNAP1='fd1.out' SNAP3='fd3.out' # Files with effective elastic parameters A11='a11.out' B11='b11.out' C11='c11.out' A13='a13.out' B13='b13.out' C13='c13.out' A31='a31.out' B31='b31.out' C31='c31.out' A33='a33.out' B33='b33.out' C33='c33.out' DEN='den.out' # NEWPAR=1 AA13='aa13.out' BB13='bb13.out' CC13='cc13.out' # Velocities for snapshot processing VP='vp.out' VS='vs.out' # Data to control seismogram plotting (program SP) SP1='fd1.ps' SP2=' ' SP3='fd3.ps' KODESP=1 SPCHRH=0.25 SPTMIN=0.000 SPTMAX=1.375 SPTLEN=13.75 SPTDIV=13.75 SPTSUB=2 SPXMIN=-10 SPXMAX=120 SPXLEN=13.0 SPXDIV=13 SPXSUB=1 NORMSP=1 SPAMP=0.08 # Programs to execute modfd: 'fd.h' / # gridded elastic parameters #fdmod: 'fd.h' / # just for test purposes fd2d: 'fd.h' / # FD synthetic seismograms sp: 'fd.h' / # seismogram plotting #....................................................................... # Programs to execute to prepare snapshots for plotting GRD='fd1.out' GRD1='fd11.out' GRD3='fd13.out' grdfd: 'fd.h' / GRD='fd3.out' GRD1='fd31.out' GRD3='fd33.out' grdfd: 'fd.h' / CAL='snapdiv.cal' GRD1='fd11.out' GRD2='fd33.out' GRD3='vp.out' GRD4='den.out' GRD5='fdp.out' grdcal: 'fd.h' / CAL='snaprot.cal' GRD1='fd13.out' GRD2='fd31.out' GRD3='vs.out' GRD4='den.out' GRD5='fds.out' grdcal: 'fd.h' / # Universal parameters to display snapshots UNIT='pt' YSIGN=-1 VREF=0.000 CREF=-.208333 # Colour of positive divergence is -5/24 VCIRC=6.283185 VWHITE=0.000 R=0 G=0 B=0 # Undefined values will be black # Parameters to display these snapshots SAMP=0.500 # S-wave amplification factor # 2-D equipartition law: SAMP=VS/VP # Seismic force: SAMP=(VS/VP)**1.5 VSAT=193178 VSAT4=2.000 # Colour saturation (amplitude) control # VSAT=SQRT(1600*(100/0.35)**3) N1NEW=1 N2NEW=1 N3NEW=1 GNORM=999. fdmax.pl: fd.h fdp.out fds.out fdmax.ps # Now display 'fdmax.ps' to check the values of SAMP, VSAT and VSAT4! # Denoting receivers by undefined values and making them black PTS='fdrec.dat' PTSGRD='ptsgrd.out' #ptsgrd: 'fd.h' / CAL='addrec.cal' GRD1='fdp.out' GRD2='ptsgrd.out' GRD3='fdprec.out' #grdcal: 'fd.h' / # Programs to execute to plot snapshots CAL='snap.cal' GRD1='fdprec.out' GRD2='fds.out' GRD3='fda.out' GRD4='fdn.out' #grdcal: 'fd.h' / GRD='fda.out' PS='fd000001.ps' GRDPS2='fdn.out' #grdps: 'fd.h' / # Mask for snapshot PS files and the name of the animated GIF file SNAPS='fd00*.ps' MOVIE='fd.gif' # Conversion of PS files to animated GIF (special software needed) #fdpcx.pl: fd.h # conversion to PCX #fdgiflux.pl: fd.h # conversion to GIF (Linux) #fdgifdos.pl: fd.h # conversion to GIF (MS-DOS) fd.htm 0100666 0000766 0000766 00000030437 07310045616 011545 0 ustar klimes klimes
By Jiri Zahradnik
Package FD contains programs for 2-D P-SV elastic second-order finite differences.
See file fdver.htm for the list of released versions and changes made in this version.
Package FD employs package FORMS for unified memory management and unified compilation.
All Fortran 77 source code and include files of the FD package are assumed to be located in a single directory together with all source code and include files of the FORMS and MODEL packages when being compiled and linked. The files with main programs contain, at their ends, Fortran 90 INCLUDE command for all subroutine files required. In this way, each program may simply be compiled and linked as a single file.
All filenames are assumed to be expressed in lowercase.
Fortran 77 source code files have extension '.for'. The corresponding files with specifications of the COMMON blocks have extension '.inc' and are included in the Fortran 77 source code by means of Fortran 90 INCLUDE command.
Package FD consists of packages FORMS and MODEL and of the following Fortran 77 and Perl source code and demo files:
The development of this package has been partially supported by:
# # Perl script to display the colour saturation of the snapshots # # Version: 5.50 # Date: 2001, March 27 # # Usage: # perl fdmax.pl SEP DIV ROT PS # Parameters: # SEP... Name of the history file. # DIV... Filename of the 4-D data cube of the "normalized" wavefield # divergence (P waves). # ROT... Filename of the 4-D data cube of the X2 component of the # "normalized" wavefield rotation (S waves). # PS... Output PostScript figure to be displayed and edited. # Refer to 'print' commands at the end of this file. # # ====================================================================== # @ARGVsave=@ARGV; @ARGV=(); require 'go.pl'; &CHK("forms/","sep.pl"); require 'sep.pl'; @ARGV=@ARGVsave; # &CHK("fd/","fdmax.pro"); &CHK("" ,"$ARGV[0]" ); &CHK("" ,"$ARGV[1]" ); &CHK("" ,"$ARGV[2]" ); # &RSEP1($ARGV[0]); &RSEP3("N4" ,$N4 ,1); &RSEP3("D4" ,$D4 ,1); &RSEP3("SAMP" ,$SAMP ,1); &RSEP3("VSAT" ,$VSAT ,1); &RSEP3("VSAT4",$VSAT4,0); # &ECHO( ">$ARGV[3]","%!"); &ECHO(">>$ARGV[3]","/N4 $N4 def"); &ECHO(">>$ARGV[3]","/D4 $D4 def"); &ECHO(">>$ARGV[3]","/SAMP $SAMP def"); &ECHO(">>$ARGV[3]","/VSAT $VSAT def"); &ECHO(">>$ARGV[3]","/VSAT4 $VSAT4 def"); # open(LU,">fdmax.h"); print LU "N1=111 N3=41 N4=321 \n"; print LU "D1=1.000 D3=1.000 D4=0.004300 \n"; print LU "O1=0.000 O3=0.000 O4=0.000 \n"; print LU "N1NEW=1 N2NEW=1 N3NEW=1 GNORM=999. \n"; print LU "GRD='$ARGV[1]' GRDNEW='fdmax.tmp' \n"; close (LU) || die "Error in writing history file fdmax.h"; &RUN("grdnorm","'fdmax.h' /"); &APPEND("$ARGV[3]","fdmax.tmp"); # open(LU,">>fdmax.h"); print LU "GRD='$ARGV[2]' GRDNEW='fdmax.tmp' \n"; close (LU) || die "Error in writing history file fdmax.h"; &RUN("grdnorm","'fdmax.h' /"); &APPEND("$ARGV[3]","fdmax.tmp"); &APPEND("$ARGV[3]","fdmax.pro"); # print "-----------------------------------------------------------------\n"; print "Display and edit file '$ARGV[3]' to adjust the colour saturation!\n\n"; print "Black line in figure '$ARGV[3]' corresponds to the maximum colour\n"; print "saturation, green dots to P waves and blue dots to S waves.\n"; print "Change SAMP to amplify S waves with respect to P waves,\n"; print "VSAT to adjust the initial level of the saturation line and\n"; print "VSAT4 to adjust its exponential decay.\n\n"; print "When figure '$ARGV[3]' is OK, modify SAMP, VSAT and VSAT4 in the\n"; print "history file accordingly!\n"; # ====================================================================== 1; #fdmax.pro 0100666 0000766 0000766 00000001712 06606556772 012276 0 ustar klimes klimes %! % *** Insert definitions of N4, D4, SAMP, VSAT, VSAT4 here *** % % *** Insert P and S norms here *** % %----------------------------------------------------------- /cm {2.54 div 72 mul} def /HSIZE 16 cm def /VSIZE 10 cm def /HOFFSET 2.5 cm def /VOFFSET 2.5 cm def /DOTRADIUS 0.050 cm def %----------------------------------------------------------- /EXP 2.718281828 D4 VSAT4 mul neg exp def /DOT {DOTRADIUS 0 360 arc fill} def % S wave norms: 0 0 1 setrgbcolor N4 1 sub -1 0 { N4 1 sub div HSIZE mul HOFFSET add exch SAMP mul VSAT div VSIZE mul VOFFSET add DOT } for % P wave norms: 0 1 0 setrgbcolor N4 1 sub -1 0 { N4 1 sub div HSIZE mul HOFFSET add exch VSAT div VSIZE mul VOFFSET add DOT } for % Plot of the maximum colour saturation: 0 0 0 setrgbcolor HOFFSET HSIZE add VOFFSET moveto HOFFSET VOFFSET lineto 0 1 N4 1 sub {dup N4 1 sub div HSIZE mul HOFFSET add exch EXP exch exp VSIZE mul VOFFSET add lineto } for closepath stroke showpage %%EOF fdmod.dat 0100666 0000766 0000766 00000003427 06477135164 012237 0 ustar klimes klimes Test example NOBL= 6 NOED,VS,PSRAT= 4 100.00 2. 1600. XED,ZED,NOVIC= 20.00 0.00 2 XED,ZED,NOVIC= 60.00 5.00 3 XED,ZED,NOVIC= 110.00 5.00 0 XED,ZED,NOVIC= 110.00 0.00 0 NOED,VS,PSRAT= 6 150.00 2. 1700. XED,ZED,NOVIC= 10.00 0.00 4 XED,ZED,NOVIC= 30.00 15.00 5 XED,ZED,NOVIC= 50.00 20.00 6 XED,ZED,NOVIC= 60.00 20.00 3 XED,ZED,NOVIC= 60.00 5.00 1 XED,ZED,NOVIC= 20.00 0.00 0 NOED,VS,PSRAT= 4 150.00 2. 1700. XED,ZED,NOVIC= 60.00 5.00 2 XED,ZED,NOVIC= 60.00 20.00 6 XED,ZED,NOVIC= 110.00 20.00 0 XED,ZED,NOVIC= 110.00 5.00 1 NOED,VS,PSRAT= 5 350.00 2. 1800. XED,ZED,NOVIC= 0.00 0.00 0 XED,ZED,NOVIC= 0.00 40.00 0 XED,ZED,NOVIC= 30.00 40.00 5 XED,ZED,NOVIC= 30.00 15.00 2 XED,ZED,NOVIC= 10.00 0.00 0 NOED,VS,PSRAT= 4 350.00 2. 1800. XED,ZED,NOVIC= 30.00 15.00 4 XED,ZED,NOVIC= 30.00 40.00 0 XED,ZED,NOVIC= 50.00 40.00 6 XED,ZED,NOVIC= 50.00 20.00 2 NOED,VS,PSRAT= 5 350.00 2. 1800. XED,ZED,NOVIC= 50.00 20.00 5 XED,ZED,NOVIC= 50.00 40.00 0 XED,ZED,NOVIC= 110.00 40.00 0 XED,ZED,NOVIC= 110.00 20.00 3 XED,ZED,NOVIC= 60.00 20.00 2 fdmod.for 0100666 0000766 0000766 00000104024 06701034672 012240 0 ustar klimes klimes C
C Program FDMOD to generate finite-difference effective parameters C (harmonic averages of elastic parameters) in 2-D models specified in C terms of homogeneous convex polygons C C Version: 5.30 C Date: 1999, April 2 C C Program written by Jiri Zahradnik C Department of Geophysics, Charles University Prague C V Holesovickach 2, 180 00 Praha 8, Czech Republic C Tel.: +420-2-21912546, Fax: +420-2-21912555 C E-mail: jz@karel.troja.mff.cuni.cz C Data input and output modified by Ludek Klimes C C....................................................................... C C Contents: C Description of the polygonal model C Description of data files C List of external procedures C C For the description of effective elastic parameters refer to file C modfd.for. C C....................................................................... C C Description of the polygonal model: C C The medium is composed of homogeneous polygonal blocks. The C blocks do not overlap each other and cover the computational region. C Convex blocks only are allowed (internal angles.le.180 degrees). C For this reason even some physically homogeneous parts C of the model must be subdivided in several (formal) blocks. C Numbering of the blocks is arbitrary. C Each block is characterized by the physical coordinates of its C nodes (in meters), S-wave velocity (m/s), P/S velocity ratio, C and density (kg/m**3). The nodes are numbered counterclockwise. C Also specified for each block segment is the serial number of the C neighbouring block. C See the accompanying example, or ask the author of the program C for more details. C C Description and remarks on FD method C C....................................................................... C C Description of data files: C C Main input data read from the * input device: C The data are read in by the list directed input (free format), by C a single READ statement. The data consist of character strings C which have to be enclosed in apostrophes. Note that the C interactive * external unit may be piped or redirected to a file. C The following items are read: C (1) 'FDSEP' / C 'FDSEP'... String in apostrophes containing the name of the main C input data file with FD parameters in the SEP format. C The file specifies the filename of the polygonal model, C dimensions of the grid of N1*N2*N3 gridpoints, the form C and names of the files with effective elastic parameters, C etc. C /... It is recommended to terminate these data by a slash. C Default: 'FDSEP'='fd.h'. C C Data file 'FDSEP' has the form of the SEP (Stanford Exploration C Project) 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 the model by means of the MODEL package: C FDMOD='string'... Name of the input data file describing a 2-D C model composed of homogeneous convex polygons. C The data file contains the name of the model, the number C of blocks, and for each block: C number of edges of the block, shear-wave velocity C (e.g., in m/s) inside the block, P/S velocity ratio inside C the block, density (e.g., in kg/m**3) inside the block, C x and z coordinates (e.g., in m) of the initial point of C segment, sequential number of the neighbouring block C (outside put 0). C Parameter FDMOD is required by this program (but no other C one). C Example of data set FDMOD C Default: FDMOD=' ' C This program understands the following parameters common to nearly all C programs dealing with SEP-like gridded data formats: C N1=positive integer... Number of gridpoints of the basic grid C along the X1 axis. Default: N1=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 O3=real... X3 coordinate of the origin of the grid. C Note that FD program 'fd2d.for' assumes the free Earth C surface at vertical coordinate X3=O3. C Default: O3=0. C D1=real... Grid spacing along the X1 axis. Default: D1=1. C D3=real... Grid spacing along the X3 axis. Default: D3=1. C The grid intervals may also be negative. C And the following parameters specific to effective-parameter FD codes: C FORM='string'... Form of the output files: C FORM='formatted': Formatted ASCII output files, C FORM='unformatted': Unformatted output files. C A11='string', B11='string', C11='string'... Strings in apostrophes C containing the names of the output files with (N1-1)*N2*N3 C values of the elastic parameters averaged in the direction C the X1 axis, along the gridlines. C Defaults: A11=' ', B11=' ', C11=' '. C A33='string', B33='string', C33='string'... Strings in apostrophes C containing the names of the output files with N1*N2*(N3-1) C values of the elastic parameters averaged in the direction C the X3 axis, along the gridlines. C Defaults: A33=' ', B33=' ', C33=' '. C A13='string', B13='string', C13='string'... Strings in apostrophes C containing the names of the output files with C (N1-1)*N2*(N3-1) values of the elastic parameters averaged C in the direction the X1 axis, along the gridlines shifted C half a grid interval in the direction of the X3 axis. C Defaults: A13=' ', B13=' ', C13=' '. C A31='string', B31='string', C31='string'... Strings in apostrophes C containing the names of the output files with C (N1-1)*N2*(N3-1) values of the elastic parameters averaged C in the direction the X3 axis, along the gridlines shifted C half a grid interval in the direction of the X1 axis. C Defaults: A31=' ', B31=' ', C31=' '. C DEN='string'... String in apostrophes containing the name of the C output ASCII file with N1*N2*N3 grid values of the C density. C Default: DEN=' '. C QFC='string'... String in apostrophes containing the name of the C output ASCII file with N1*N2*N3 grid values of the P-wave C quality factor. C Default: QFC=' ' C Example of data set FDSEP C C....................................................................... C C This file consists of the following external procedures: C FDMOD...Main program allocating the memory for FDPAR. C FDMOD C FDPAR...Subroutine to generate the parameters (originally the main C program). C FDPAR C XLEG... Subroutine for computation of effective (mean) squared C velocities for a grid leg. C XLEG C GRPT... Subroutine for location of grid point (XPT,ZPT) with C respect to block IBL. C GRPT C TWOPT...Subroutine for location of a grid leg with respect to C boundary of block IB. C TWOPT C C======================================================================= C C C PROGRAM FDMOD C C Main program allocating the memory. C C Subroutines and external functions required: EXTERNAL FDPAR,RSEP1,RSEP3I C FDPAR...This file. C RSEP1,RSEP3I... File 'sep.for'. C Note that the above subroutines reference many other external C procedures from various subroutine files. These indirectly C referenced procedures are not named here, but are listed in the C particular subroutine files. C C----------------------------------------------------------------------- C C Common block /RAMC/ to allocate large array RAM: INCLUDE 'ram.inc' C ram.inc C C----------------------------------------------------------------------- C C Input data filenames and logical unit numbers: INTEGER LU PARAMETER (LU=1) CHARACTER*80 FSEP C LU... Logical unit number to be used to read and write all data C files. C FSEP... The name of the input data file with FD parameters in SEP C format. C INTEGER N1,N3,N,NR,NS,NGRID,I1,I2,I3,I4 C C....................................................................... C C Main input data: C Default: FSEP='fd.h' C Reading main input data: WRITE(*,'(A)') ' FDMOD: Enter name of input SEP file [''fd.h'']: ' READ (*,*) FSEP WRITE(*,'(A)') '+FDMOD: Working... ' C C Reading all the data from input SEP parameter file: CALL RSEP1(LU1,FSEP) C C Data specifying numbers of gridpoints: CALL RSEP3I('N1' ,N1 ,1) CALL RSEP3I('N2' ,N3 ,1) C C Allocating the memory in array RAM: NGRID=15*N1*N3 IF(NGRID.GT.MRAM) THEN C FDMOD-01 CALL ERROR('FDMOD-01: Small dimension MRAM in ram.inc') C Dimension MRAM must be at least C 15*N1*N3 , C where N1 and N3 are parameters in the main input file. END IF N =N1*N3 C CALL FDPAR(LU,N, * RAM( 1),RAM( N + 1),RAM( 2*N + 1),RAM( 3*N + 1), * RAM( 4*N + 1),RAM( 5*N + 1),RAM( 6*N + 1),RAM( 7*N + 1), * RAM( 8*N + 1),RAM( 9*N + 1),RAM(10*N + 1),RAM(11*N + 1), * RAM(12*N + 1),RAM(13*N + 1),RAM(14*N + 1)) C WRITE(*,'(A)') '+FDMOD: Done. ' STOP END C C======================================================================= C C C SUBROUTINE FDPAR(LU,KLMAX, * AHS,BHS,CHS,AVS,BVS,CVS, * AHL,BHL,CHL,AVL,BVL,CVL,DEN,QFC,WORK) INTEGER LU,KLMAX REAL AHS(KLMAX),BHS(KLMAX),CHS(KLMAX) REAL AVS(KLMAX),BVS(KLMAX),CVS(KLMAX) REAL AHL(KLMAX),BHL(KLMAX),CHL(KLMAX) REAL AVL(KLMAX),BVL(KLMAX),CVL(KLMAX) REAL DEN(KLMAX),QFC(KLMAX),WORK(KLMAX) C C Subroutine (originally main program) to read the input data and C generate finite-difference effective parameters. C C Input: C LU... Logical unit number to be used to read and write data C files. C KLMAX...Array dimension. C AHL,BHL,CHL, C AVL,BVL,CVL, C AHS,BHS,CHS, C AVS,BVS,CVS, C DEN,QFC...arrays with effective parameters C WORK... temporary working array C C No output. C C Subroutines and external functions required: EXTERNAL XLEG,GRPT,TWOPT EXTERNAL RSEP3T,FDGRID,FDOUT1 C XLEG... This file. C GRPT... This file. C TWOPT.. This file. C RSEP3T..File 'sep.for'. C FDGRID..File 'fdaux.for'. C FDOUT1..File 'fdaux.for'. C Note that the above subroutines reference many other external C procedures from various subroutine files. These indirectly C referenced procedures are not named here, but are listed in the C particular subroutine files. C C C----------------------------------------------------------------------- C C Storage locations: C C Maximum array dimensions: PARAMETER (MAXBL=50,MAXEDG=10) C MAXBL.. Max number of polygonal blocks. C MAXEDG..Max number of edges for any block. C C Arrays: DIMENSION NOED(0:MAXBL),VS(0:MAXBL),PSRAT(0:MAXBL),HU(0:MAXBL) DIMENSION XED(MAXEDG,MAXBL),ZED(MAXEDG,MAXBL),NOVIC(MAXEDG,MAXBL) C COMMON /BLO1/ NOBL,NOED,VS,PSRAT,DX,HU COMMON /BLO2/ XED,ZED,NOVIC C NOBL... Number of blocks. C NOED... Number of edges of the block. C VS... Shear-wave velocity (m/s) inside the block. C PSRAT.. P/S velocity ratio inside the block. C HU... Density (kg/m**3) inside the block. C XED,ZED. X and Z coordinates (m) of the initial point of segment. C NOVIC... Serial number of the neigh. block (outside put 0). C CHARACTER*80 FMOD CHARACTER *5 NAMDUM C FMOD... The name of the input data file with model parameters in C list directed (free) format. C NAMDUM... Auxiliary storage location. C C....................................................................... C C Data specifying grid dimensions CALL FDGRID(KLMAX,KLMAX,KLMAX,NX,NZ,DX) C C Reading data describing a model composed of homogeneous polygons: CALL RSEP3T('FDMOD',FMOD,' ') OPEN(LU,FILE=FMOD,STATUS='OLD') C READ(1,'(1X,A5)') NAMDUM C READ(1,'(25X,I5)') NSYM NSYM=0 C READ(1,'(25X,I5)') NOBL DO 150 I=1,NOBL READ(1,'(25X,I5,3F10.0)') NOED(I),VS(I),PSRAT(I),HU(I) DO 100 J=1,NOED(I) READ(1,'(25X,2F10.0,I5)') XED(J,I),ZED(J,I),NOVIC(J,I) 100 CONTINUE 150 CONTINUE C WRITE(*,'(A)') '+FDMOD: Computing effective parameters ' C C Auxiliary parameters. LACC=0 C DO 310 I=1,NOBL C Attention: Now VS no more contains shear velocity but rigidity C and PSRAT is squared Vp/Vs ratio. The corresponding C change was made in FLERIB,FBOB,FCORNB. C VS(I)=VS(I)*VS(I) VS(I)=VS(I)*VS(I)*HU(I) PSRAT(I)=PSRAT(I)*PSRAT(I) 310 CONTINUE C NOED(0)=0 VS(0)=0. PSRAT(0)=1. C DXA=DX C C Vertical and horizontal parameters of medium for INTERlegs C (i.e. for legs between grid lines) NPUP=(NX-1)*(NZ-1) DO 430 IPU=1,NPUP NVER=NZ-1 NP=IFIX((IPU-0.01)/NVER) XSOUR=DX/2.+FLOAT(NP)*DX NPOR=IPU-NP*NVER ZSOUR=FLOAT(NPOR-1)*DX CALL XLEG(0,XSOUR,ZSOUR,DXA,AL,BE,RH) AVS(IPU)=AL BVS(IPU)=BE CVS(IPU)=AL-2.*BE NHOR=NZ-1 NP=IFIX((IPU-0.01)/NHOR) XSOUR=FLOAT(NP)*DX NPOR=IPU-NP*NHOR ZSOUR=DX/2.+FLOAT(NPOR-1)*DX CALL XLEG(1,XSOUR,ZSOUR,DXA,AL,BE,RH) AHS(IPU)=AL BHS(IPU)=BE CHS(IPU)=AL-2.*BE 430 CONTINUE C C Vertical parameters of medium for LARGE squares (= grid legs) C (for ignoring detailed computations, choose LACC=0) DO 450 K=2,NX-1 XOLD=FLOAT(K-1)*DX LAU=(K-1)*NZ XNZ=FLOAT(NZ) DO 400 L=1,NZ-1 ICEN=L+LAU JCEN=ICEN-IFIX(FLOAT(ICEN)/XNZ) ZOLD=FLOAT(L-1)*DX IF(L.LT.LACC) THEN CALL XLEG(0,XOLD,ZOLD,DXA/2.,AL1,BE1,RH1) ZNEW=ZOLD+DX/2. CALL XLEG(0,XOLD,ZNEW,DXA/2.,AL2,BE2,RH2) AL=(AL1+AL2)/2. BE=(BE1+BE2)/2. ELSE CALL XLEG(0,XOLD,ZOLD,DXA,AL,BE,RH) ENDIF AVL(JCEN)=AL BVL(JCEN)=BE CVL(JCEN)=AL-2.*BE 400 CONTINUE 450 CONTINUE C C Horizontal parameters of medium for LARGE squares (=grid legs) C (incl. density) DO 550 K=1,NX-1 XOLD=FLOAT(K-1)*DX LAU=(K-1)*NZ DO 500 L=1,NZ-1 ICEN=L+LAU ZOLD=FLOAT(L-1)*DX IF(L.LE.LACC) THEN CALL XLEG(1,XOLD,ZOLD,DXA/2.,AL1,BE1,RH1) XNEW=XOLD+DX/2. CALL XLEG(1,XNEW,ZOLD,DXA/2.,AL2,BE2,RH2) AL=(AL1+AL2)/2. BE=(BE1+BE2)/2. ELSE CALL XLEG(1,XOLD,ZOLD,DXA,AL,BE,RH) ENDIF AHL(ICEN)=AL BHL(ICEN)=BE CHL(ICEN)=AL-2.*BE DEN(ICEN)=RH 500 CONTINUE 550 CONTINUE C C Vertical and horizontal parameters for legs in the strip between C the inner and outer nonreflecting boundaries (a continuation) C (incl. density for horizontal legs) DO 600 L=1,NZ-1 ICEN=L ICEN2=ICEN+NZ JCEN=L JCEN2=JCEN+NZ-1 AVL(JCEN)=AVL(JCEN2) BVL(JCEN)=BVL(JCEN2) CVL(JCEN)=CVL(JCEN2) 600 CONTINUE LAU=(NX-1)*NZ DO 610 L=1,NZ-1 ICEN=L+LAU JCEN=ICEN-IFIX(FLOAT(ICEN)/FLOAT(NZ)) JCEN2=JCEN-NZ+1 AVL(JCEN)=AVL(JCEN2) BVL(JCEN)=BVL(JCEN2) CVL(JCEN)=CVL(JCEN2) 610 CONTINUE DO 620 K=1,NX-1 ICEN=K*NZ ICEN2=ICEN-1 AHL(ICEN)=AHL(ICEN2) BHL(ICEN)=BHL(ICEN2) CHL(ICEN)=CHL(ICEN2) DEN(ICEN)=DEN(ICEN2) 620 CONTINUE C C Halving the density at the free surface IZ=1: DO 70 I=1,NX*NZ,NZ DEN(I)=0.5*DEN(I) 70 CONTINUE C C No attenuation is considered in this version DO 80 I=1,NX*NZ QFC(I)=0. 80 CONTINUE C C Writing the gridded effective elastic parameters: CALL FDOUT1(LU,NX-1,NZ ,1.,'A11',AHL,WORK) CALL FDOUT1(LU,NX-1,NZ ,1.,'B11',BHL,WORK) CALL FDOUT1(LU,NX-1,NZ ,1.,'C11',CHL,WORK) CALL FDOUT1(LU,NX ,NZ-1,1.,'A33',AVL,WORK) CALL FDOUT1(LU,NX ,NZ-1,1.,'B33',BVL,WORK) CALL FDOUT1(LU,NX ,NZ-1,1.,'C33',CVL,WORK) CALL FDOUT1(LU,NX-1,NZ-1,1.,'A13',AHS,WORK) CALL FDOUT1(LU,NX-1,NZ-1,1.,'B13',BHS,WORK) CALL FDOUT1(LU,NX-1,NZ-1,1.,'C13',CHS,WORK) CALL FDOUT1(LU,NX-1,NZ-1,1.,'A31',AVS,WORK) CALL FDOUT1(LU,NX-1,NZ-1,1.,'B31',BVS,WORK) CALL FDOUT1(LU,NX-1,NZ-1,1.,'C31',CVS,WORK) CALL FDOUT1(LU,NX ,NZ ,1.,'DEN',DEN,WORK) CALL FDOUT1(LU,NX ,NZ ,1.,'QFC',QFC,WORK) C RETURN END C C======================================================================= C C C SUBROUTINE XLEG(KEY,XPT,ZPT,DEL,ALF,BET,RHO) C C Purpose: C Effective (mean) squared velocities for a grid leg C C Input: C KEY=0.....vertical leg C =1.....horizontal leg C XPT,ZPT...x,z coordinates of the grid point C DEL.......length of the leg C C Output: C ALF,BET...P and S squared velocities for the leg C C DIMENSION NOED(0:50),VS(0:50),PSRAT(0:50),HU(0:50) DIMENSION XED(10,50),ZED(10,50),NOVIC(10,50) COMMON /BLO1/ NOBL,NOED,VS,PSRAT,DX,HU COMMON /BLO2/ XED,ZED,NOVIC C RHO=0. C C An artificial shortening of the leg to avoid numerical problems C !!!!!!Some models possible with 0.01 instead of this 0.03!!!!!! C DELSA=DEL-0.03*DEL cc DELSA=DEL cc DELSA=DEL-0.0001*DEL C C Location of the grid point with respect to a block C DO 100 J=1,NOBL DIST=0. CALL GRPT(XPT,ZPT,J,IPRES,IVIC,DIST,ISEG,IBLVC,IORT) DIST=DIST/DELSA IF(IPRES.EQ.0) GOTO 100 IBL=J GOTO 200 100 CONTINUE WRITE(*,'(A,2F10.3/A)') * '+FDMOD: Gridpoint not located successfully at x,z=',XPT,ZPT,' ' 200 CONTINUE IF(IPRES.EQ.2) GOTO 1000 C C The grid point inside block IBL C RHO=HU(IBL) CALL TWOPT(XPT,ZPT,IBL,DELSA,KEY,0,ISG,IVC,PART) IF(ISG.NE.0) GOTO 400 C C No intersection of the leg with any boundary segment C 300 BET=VS(IBL) ALF=BET*PSRAT(IBL) RETURN C C The leg is intersecting segment ISG C Locating the end point of the leg with respect to block IVC C 400 IF(KEY.EQ.0) THEN XEND=XPT ZEND=ZPT+DELSA ELSE XEND=XPT+DELSA ZEND=ZPT ENDIF DIS2=0. CALL GRPT(XEND,ZEND,IVC,IPR2,IVI2,DIS2,ISE2,IBL2,IOR2) DIS2=DIS2/DELSA IF(IPR2.EQ.0) GOTO 600 IF(IPR2.EQ.1) GOTO 450 IF(IPR2.EQ.2.AND.KEY.EQ.0.AND.IOR2.EQ.0) GOTO 500 IF(IPR2.EQ.2.AND.KEY.EQ.1.AND.IOR2.EQ.1) GOTO 500 C C The end point located inside block IVC C (or on its boundary segment not coinciding with the leg) C 450 VSIN=VS(IBL) VSOU=VS(IVC) VPIN=VSIN*PSRAT(IBL) VPOU=VSOU*PSRAT(IVC) c BET=PART*VSIN+(1.-PART)*VSOU c ALF=PART*VPIN+(1.-PART)*VPOU BET=1./(PART/VSIN+(1.-PART)/VSOU) ALF=1./(PART/VPIN+(1.-PART)/VPOU) RETURN C C The end point located on a vertical/horizontal boundary C between blocks IVC and IBL2 C 500 VSIN=VS(IBL) VSOU1=VS(IVC) VSOU2=VS(IBL2) VSOU=(VSOU1+VSOU2)/2. VPIN=VSIN*PSRAT(IBL) VPOU=(VSOU1*PSRAT(IVC)+VSOU2*PSRAT(IBL2))/2. c BET=PART*VSIN+(1.-PART)*VSOU c ALF=PART*VPIN+(1.-PART)*VPOU BET=1./(PART/VSIN+(1.-PART)/VSOU) ALF=1./(PART/VPIN+(1.-PART)/VPOU) RETURN C C The end point located outside block IVC,inside block IBLNEW C 600 DO 700 J=1,NOBL CALL GRPT(XEND,ZEND,J,IPR3,IVI3,DIS3,ISE3,IBL3,IOR3) IF(IPR3.EQ.0) GOTO 700 IBLNEW=J GOTO 800 700 CONTINUE WRITE(*,'(A,2F10.3/A)') *'+FDMOD: Gridpoint not located successfully at x,z=',XEND,ZEND,' ' C C The end point in block IBLNEW, the leg in blocks IBL, IVC, IBLNEW C 800 CALL TWOPT(XEND,ZEND,IBLNEW,DELSA,KEY,1,ISG2,IVC2, PART2) VSIN=VS(IBL) VSOU1=VS(IVC) VSOU2=VS(IBLNEW) VPIN=VSIN*PSRAT(IBL) VPOU1=VSOU1*PSRAT(IVC) VPOU2=VSOU2*PSRAT(IBLNEW) c BET=PART*VSIN+(1.-PART-PART2)*VSOU1+PART2*VSOU2 c ALF=PART*VPIN+(1.-PART-PART2)*VPOU1+PART2*VPOU2 BET=1./(PART/VSIN+(1.-PART-PART2)/VSOU1+PART2/VSOU2) ALF=1./(PART/VPIN+(1.-PART-PART2)/VPOU1+PART2/VPOU2) RETURN C C The grid point located on the boundary of block IBL C (on segment ISEG, neighbouring to block IBLVC) C Locating the end point with respect to block IBL C 1000 CONTINUE if(ibl.eq.0) rho=hu(iblvc) if(iblvc.eq.0) rho=hu(ibl) if(ibl.ne.0.and.iblvc.ne.0) RHO=(HU(IBL)+HU(IBLVC))/2. IF(KEY.EQ.0) THEN XEND=XPT ZEND=ZPT+DELSA ELSE XEND=XPT+DELSA ZEND=ZPT ENDIF DIS2=0. CALL GRPT(XEND,ZEND,IBL,IPR2,IVI2,DIS2,ISE2,IBL2,IOR2) DIS2=DIS2/DELSA IF(IPR2.EQ.0) GOTO 1800 IF(IPR2.EQ.1) GOTO 1600 IF(IPR2.EQ.2.AND.KEY.EQ.0.AND.IOR2.EQ.0) GOTO 1700 IF(IPR2.EQ.2.AND.KEY.EQ.1.AND.IOR2.EQ.1) GOTO 1700 C C The end point located inside block IBL C (or on its boundary segment not coinciding with the leg) C 1600 BET=VS(IBL) ALF=BET*PSRAT(IBL) RETURN C C The leg is coinciding with a segment C 1700 VSIN=VS(IBL) C VSOU=VS(IBLVC) VSOU=VS(IBL2) VPIN=VSIN*PSRAT(IBL) C VPOU=VSOU*PSRAT(IBLVC) VPOU=VSOU*PSRAT(IBL2) BET=(VSIN+VSOU)/2. ALF=(VPIN+VPOU)/2. RETURN C C The end point located outside block IBL C Locating the end point with respect to block IBLVC C 1800 DIS3=0. CALL GRPT(XEND,ZEND,IBLVC,IPR3,IVI3,DIS3,ISE3,IBL3,IOR3) DIS3=DIS3/DELSA IF(IPR3.EQ.0) GOTO 2100 IF(IPR3.EQ.1) GOTO 1850 IF(IPR3.EQ.2) GOTO 1900 C C The end point located inside block IBLVC C 1850 BET=VS(IBLVC) ALF=BET*PSRAT(IBLVC) RETURN C C The end point located on segment ISE3 of block IBLVC C 1900 IF(KEY.EQ.0.AND.IOR3.EQ.0) GOTO 2000 IF(KEY.EQ.1.AND.IOR3.EQ.1) GOTO 2000 C 1950 BET=VS(IBLVC) ALF=BET*PSRAT(IBLVC) RETURN C C The end point located on a vertical/horizontal boundary C between blocks IBLVC and IBL3 C 2000 VSIN=VS(IBL) VSOU=(VS(IBLVC)+VS(IBL3))/2. VPIN=VS(IBL)*PSRAT(IBL) VPOU=(VS(IBLVC)*PSRAT(IBLVC)+VS(IBL3)*PSRAT(IBL3))/2. c BET=DIST*VSIN+(1.-DIST)*VSOU c ALF=DIST*VPIN+(1.-DIST)*VPOU BET=1./(DIST/VSIN+(1.-DIST)/VSOU) ALF=1./(DIST/VPIN+(1.-DIST)/VPOU) RETURN C C The end point located outside block IBLVC, inside block IBLNEW C 2100 DO 2150 J=1,NOBL DIS5=0. CALL GRPT(XEND,ZEND,J,IPR5,IVI5,DIS5,ISE5,IBL5,IOR5) DIS5=DIS5/DELSA IF(IPR5.EQ.0) GOTO 2150 IBLNEW=J GOTO 2170 2150 CONTINUE 2170 CONTINUE IF(KEY.EQ.0.AND.IORT.EQ.0) GOTO 2200 IF(KEY.EQ.1.AND.IORT.EQ.1) GOTO 2200 GOTO 2400 2200 IF(IPR5.EQ.1) GOTO 2300 IF(IPR5.EQ.2) GOTO 2350 C C The grid point located on a vertical/horizontal boundary C between blocks IBL and IBLVC C The end point inside block IBLNEW C 2300 VSIN=(VS(IBL)+VS(IBLVC))/2. VSOU=VS(IBLNEW) VPIN=(VS(IBL)*PSRAT(IBL)+VS(IBLVC)*PSRAT(IBLVC))/2. VPOU=VSOU*PSRAT(IBLNEW) c BET=DIST*VSIN+(1.-DIST)*VSOU c ALF=DIST*VPIN+(1.-DIST)*VPOU BET=1./(DIST/VSIN+(1.-DIST)/VSOU) ALF=1./(DIST/VPIN+(1.-DIST)/VPOU) RETURN C C The grid point located on a vertical/horizontal boundary C between blocks IBL and IBLVC C The end point located on a vertical/horizontal boundary C between blocks IBLNEW and IBL5 C 2350 VSIN=(VS(IBL)+VS(IBLVC))/2. VSOU=(VS(IBLNEW)+VS(IBL5))/2. VPIN=(VS(IBL)*PSRAT(IBL)+VS(IBLVC)*PSRAT(IBLVC))/2. VPOU=(VS(IBLNEW)*PSRAT(IBLNEW)+VS(IBL5)*PSRAT(IBL5))/2. c BET=DIST*VSIN+(1.-DIST)*VSOU c ALF=DIST*VPIN+(1.-DIST)*VPOU BET=1./(DIST/VSIN+(1.-DIST)/VSOU) ALF=1./(DIST/VPIN+(1.-DIST)/VPOU) RETURN C C The end point located inside block IBLNEW, C the leg in blocks IBLNEW and IVC4 (where IVC4=IBL or IBLVC) C 2400 IF(IPR5.EQ.1) GOTO 2600 IF(IPR5.EQ.2.AND.KEY.EQ.0.AND.IOR5.EQ.0) GOTO 2800 IF(IPR5.EQ.2.AND.KEY.EQ.1.AND.IOR5.EQ.1) GOTO 2800 2600 CALL TWOPT(XEND,ZEND,IBLNEW,DELSA,KEY,1,ISG4,IVC4,PART4) VSIN=VS(IVC4) VSOU=VS(IBLNEW) VPIN=VSIN*PSRAT(IVC4) VPOU=VSOU*PSRAT(IBLNEW) c BET=(1.-PART4)*VSIN+PART4*VSOU c ALF=(1.-PART4)*VPIN+PART4*VPOU c some models with inclined boundaries ending at left or right sides c require the following modification: c if(vsin.eq.0.) then c BET=1./(PART4/VSOU) c ALF=1./(PART4/VPOU) c RETURN c endif BET=1./((1.-PART4)/VSIN+PART4/VSOU) ALF=1./((1.-PART4)/VPIN+PART4/VPOU) RETURN 2800 CALL TWOPT(XEND,ZEND,IBL,DELSA,KEY,1,ISG5,IVC5,PART5) IF(ISG5.NE.0) THEN IVXX=IBL PARTX=PART5 GOTO 2900 ELSE CALL TWOPT(XEND,ZEND,IBLVC,DELSA,KEY,1,ISG6,IVC6,PART6) IF(ISG6.NE.0) THEN IVXX=IBLVC PARTX=PART6 GOTO 2900 ELSE c 2850 BET=(VS(IBLNEW)+VS(IBL5))/2. ALF=(VS(IBLNEW)*PSRAT(IBLNEW)+VS(IBL5)*PSRAT(IBL5))/2. RETURN ENDIF ENDIF 2900 VSIN=VS(IVXX) VSOU=(VS(IBLNEW)+VS(IBL5))/2. VPIN=VSIN*PSRAT(IVXX) VPOU=(VS(IBLNEW)*PSRAT(IBLNEW)+VS(IBL5)*PSRAT(IBL5))/2. c BET=(1.-PARTX)*VSIN+PARTX*VSOU c ALF=(1.-PARTX)*VPIN+PARTX*VPOU BET=1./((1.-PARTX)/VSIN+PARTX/VSOU) ALF=1./((1.-PARTX)/VPIN+PARTX/VPOU) RETURN END C DEBUG SUBCHK C END DEBUG C C======================================================================= C C C SUBROUTINE GRPT(XPT,ZPT,IBL,IPRES,IVIC,DIST,ISEG,IBLVC,IORT) C C Purpose: C Location of grid point (XPT,ZPT) with respect to block IBL C C Input: C XPT,ZPT......x,z coordinates of the grid point C IBL..........the block to be tested C C Output: C IPRES=0......the grid point located outside the block C =1......the grid point located inside the block C =2......the grid point located on the block boundary C C IVIC=0.......the grid point not in vicinity (.LT.DX) of any edge C .ne.0.......the grid point in vicinity of edge IVIC C C DIST.........distance from the edge (for IVIC.ne.0) C C ISEG.........the grid point on segment ISEG (for IPRES=2) C C IBLVC........the block neighbouring to ISEG (for IPRES=2) C C IORT=0.......a vertical segment (for IPRES=2) C =1.......a horizontal segment (for IPRES=2) C =2.......an oblique segment (for IPRES=2) C C DIMENSION NOED(0:50),VS(0:50),PSRAT(0:50),HU(0:50) DIMENSION XED(10,50),ZED(10,50),NOVIC(10,50) COMMON /BLO1/ NOBL,NOED,VS,PSRAT,DX,HU COMMON /BLO2/ XED,ZED,NOVIC C I=IBL IF(I.EQ.0) THEN IPRES=0 RETURN ENDIF NEDG=NOED(I) DO 200 J=1,NEDG C C Unit position vector (PX,PZ) of grid point (XPT,ZPT) C with respect to the origin of J-th segment of I-th block C PX=XPT-XED(J,I) PZ=ZPT-ZED(J,I) SIZ=SQRT(PX*PX+PZ*PZ) IF(ABS(SIZ).LT.0.00001) GOTO 300 PX=PX/SIZ PZ=PZ/SIZ C C Unit vector (BX,BZ) of J-th segment of I-th block C IF(J.LT.NEDG) THEN BX=XED(J+1,I)-XED(J,I) BZ=ZED(J+1,I)-ZED(J,I) ELSE BX=XED(1,I)-XED(J,I) BZ=ZED(1,I)-ZED(J,I) ENDIF SIZ=SQRT(BX*BX+BZ*BZ) BX=BX/SIZ BZ=BZ/SIZ C C Y-component (=PROD) of the vector product (PX,PZ)x(BX,BZ) C PROD=PZ*BX-PX*BZ IF(PROD.LT.-0.00001) GOTO 200 IF(ABS(PROD).LT.0.00001) GOTO 300 C C The grid point located outside the block C 250 IPRES=0 GOTO 600 C C The grid point possibly located on J-th segment C 300 IPRES=2 ISEG=J IBLVC=NOVIC(J,I) C C End points of J-th segment, orientation of the segment C X1=XED(J,I) Z1=ZED(J,I) IF(J.LT.NEDG) THEN X2=XED(J+1,I) Z2=ZED(J+1,I) ELSE X2=XED(1,I) Z2=ZED(1,I) ENDIF SEGLEN=SQRT((X2-X1)**2+(Z2-Z1)**2) IORT=2 IF(ABS(X1-X2).LT.0.00001) IORT=0 IF(ABS(Z1-Z2).LT.0.00001) IORT=1 C C Distances of the grid point from the end points C DIST1=SQRT((XPT-X1)**2+(ZPT-Z1)**2) DIST2=SQRT((XPT-X2)**2+(ZPT-Z2)**2) ROZ1=DIST1-SEGLEN ROZ2=DIST2-SEGLEN C IF(DIST1.GT.SEGLEN.OR.DIST2.GT.SEGLEN) GOTO 250 IF(ROZ1.GT.0.0001.OR.ROZ2.GT.0.0001) GOTO 250 C C The grid point certainly located on J-th segment C IVIC=0 IF(DIST1.LT.DX) GOTO 400 IF(DIST2.LT.DX) GOTO 500 GOTO 600 400 IVIC=J DIST=DIST1 GOTO 600 500 IF(J.LT.NEDG) THEN IVIC=J+1 ELSE IVIC=1 ENDIF DIST=DIST2 GOTO 600 200 CONTINUE C C The grid point 'on the left' with respect to all segments, C i.e., the grid point located inside I-th block C IPRES=1 C C Distances from all edges of the block C (The program cannot detect proximity of two edges simultan.) C IVIC=0 DO 550 J=1,NEDG DISTX=SQRT((XPT-XED(J,I))**2+(ZPT-ZED(J,I))**2) IF(DISTX.GT.DX) GOTO 550 IVIC=J DIST=DISTX GOTO 600 550 CONTINUE 600 RETURN END C DEBUG SUBCHK C END DEBUG C C======================================================================= C C C SUBROUTINE TWOPT(X,Z,IB,DELSA,KY,KYOR,ISG,IVC,PART) C C Purpose: C Location of a grid leg with respect to boundary of block IB C C Input: C X,Z........x,z coordinates of the grid point C IB.........the block to which the grid point is belonging C DELSA........length of the leg C KY=0.......vertical leg C =1.......horizontal leg C KYOR=0.....leg oriented down, or right, of the grid point C =1.....up or left C C Output: C ISG=0......no intersection of the leg with the block boundary C .ne.0....intersection with ISG-th segment C IVC........the neighbouring block (for ISG.ne.0) C PART.......a relative distance of the grid point from C the intersection (for ISG.ne.0) C C DIMENSION NOED(0:50),VS(0:50),PSRAT(0:50),HU(0:50) DIMENSION XED(10,50),ZED(10,50),NOVIC(10,50) COMMON /BLO1/ NOBL,NOED,VS,PSRAT,DX,HU COMMON /BLO2/ XED,ZED,NOVIC C IF(IB.EQ.0) THEN ISG=0 RETURN ENDIF C C End point (XX,ZZ) of the leg (the leg from (X,Z) to (XX,ZZ)) C IF(KY.EQ.0) GOTO 100 IF(KY.EQ.1) GOTO 200 100 XX=X IF(KYOR.EQ.0) THEN ZZ=Z+DELSA ELSE ZZ=Z-DELSA ENDIF GOTO 300 200 ZZ=Z IF(KYOR.EQ.0) THEN XX=X+DELSA ELSE XX=X-DELSA ENDIF C Intersection of the leg with the block boundary C (legs coinciding with boundaries are not allowed) C 300 NEDG=NOED(IB) DO 500 J=1,NEDG C C The end points (X1,Z1) and (X2,Z2) of J-th segment C X1=XED(J,IB) Z1=ZED(J,IB) IF(J.LT.NEDG) THEN X2=XED(J+1,IB) Z2=ZED(J+1,IB) ELSE X2=XED(1,IB) Z2=ZED(1,IB) ENDIF IF(KY.EQ.0) GOTO 600 IF(KY.EQ.1) GOTO 700 600 IF(X1.LT.X2) THEN XMIN=X1 XMAX=X2 ELSE XMIN=X2 XMAX=X1 ENDIF IF(X.LT.XMIN.OR.X.GT.XMAX) GOTO 500 ZSEC=Z1+(X-X1)*(Z2-Z1)/(X2-X1) IF(KYOR.EQ.0) THEN ZBOT=Z ZTOP=ZZ ELSE ZBOT=ZZ ZTOP=Z ENDIF IF(ZSEC.LT.ZBOT.OR.ZSEC.GT.ZTOP) GOTO 500 ISG=J IVC=NOVIC(J,IB) PART=ABS(ZSEC-Z)/DELSA GOTO 1000 700 IF(Z1.LT.Z2) THEN ZMIN=Z1 ZMAX=Z2 ELSE ZMIN=Z2 ZMAX=Z1 ENDIF IF(Z.LT.ZMIN.OR.Z.GT.ZMAX) GOTO 500 XSEC=X1+(Z-Z1)*(X2-X1)/(Z2-Z1) IF(KYOR.EQ.0) THEN XBOT=X XTOP=XX ELSE XBOT=XX XTOP=X ENDIF IF(XSEC.LT.XBOT.OR.XSEC.GT.XTOP) GOTO 500 ISG=J IVC=NOVIC(J,IB) PART=ABS(XSEC-X)/DELSA GOTO 1000 500 CONTINUE ISG=0 1000 RETURN END C DEBUG SUBCHK C END DEBUG C C======================================================================= C INCLUDE 'error.for' C error.for INCLUDE 'fdaux.for' C fdaux.for INCLUDE 'gmtra.for' C gmtra.for INCLUDE 'sep.for' C sep.for INCLUDE 'forms.for' C forms.for INCLUDE 'length.for' C length.for C C======================================================================= Cfdpcx.bat 0100666 0000766 0000766 00000000761 06606557000 012235 0 ustar klimes klimes @ECHO OFF IF "%2"=="" ECHO Error: No RAM disk given IF "%2"=="" GOTO :END REM ........................................................ COPY %1 %2 COPY fdpcx.pl %2 COPY sep.pl %2 COPY C:\PERL\perl.exe %2 COPY C:\PERL\perlglob.exe %2 COPY C:\GS\*.* %2 COPY fdgifdos.pl %2 %2 COPY gs386.exe gs.exe perl fdpcx.pl %1 ECHO To create GIF files type: perl fdgifdos.pl %1 REM ........................................................ :END fdpcx.pl 0100666 0000766 0000766 00000002554 06616221142 012100 0 ustar klimes klimes #!perl #
# # Perl script to convert PostScript snapshots into PCX # # Version: 5.20 # Date: 1998, October 30 # # Usage: # perl fdpcx.pl SEP # Parameters: # SEP... Name of the history file. # Parameters read from the history file: # N1=integer... Number of gridpoints (pixels) in horizotal direction. # Default: N1=1 # N3=integer... Number of gridpoints (pixels) in vertical direction. # Default: N3=1 # SNAPS=string... Mask for the names of PostScript snapshot files. # Default: SNAPS='fd00*.ps' # Include files required: # sep.pl # Programs required: # gs # # ====================================================================== # Main program 'fdpcx.pl': # ~~~~~~~~~~~~~~~~~~~~~~~~ require 'sep.pl'; $SEP="$ARGV[0]"; # &RSEP1($SEP); &RSEP3('N1',$N1,1); &RSEP3('N3',$N3,1); &RSEP3('SNAPS',$SNAPS,'fd00*.ps'); # foreach $SNAP (<${SNAPS}>) { $SNAP=~s/.ps//; print "FDPCX: converting $SNAP\n"; open(LU,"|gs -sDEVICE=pcx24b -r72 -g${N1}x${N3} -dNOPAUSE". " -sOutputFile=$SNAP.pcx $SNAP.ps -c quit"); close(LU) || die "Error"; } print "FDPCX: Done.\n"; # ====================================================================== 1; #fdrec.dat 0100666 0000766 0000766 00000000304 06573134542 012214 0 ustar klimes klimes / '000' 0 0 0 / '010' 10 0 0 / '020' 20 0 0 / '030' 30 0 0 / '040' 40 0 0 / '050' 50 0 0 / '060' 60 0 0 / '070' 70 0 0 / '080' 80 0 0 / '090' 90 0 0 / '100' 100 0 0 / '110' 110 0 0 / / fdsnaps.htm 0100666 0000766 0000766 00000001363 07170240150 012600 0 ustar klimes klimes
Demo model for 2-D finite differences
Finite-difference snapshots
The frame rate should be 5 snaps per second.
All 320 snaps should thus be displayed in 64s.
However, the frame rate may be decreased by a slow browser.
Colour scale for 2-D elastic snapshots
ffd.pl 0100666 0000766 0000766 00000000604 06607027626 011540 0 ustar klimes klimes #!perl #Released versions of package FD
5.20 (1998, October): *** First release consisting of program and subroutines files 'modfd.for','fdmod.for','fd2d.for','fdaux.for','gmtra.for', 'newpar.for','newpar.pl','fdmax.pl', demo data, scripts and history files. *** 5.30 (1999, June): 'fd2d.for', 'fdaux.for': Erroneous sign in output snapshots fixed. Comments in 'fdaux.for' improved. 'fdmod.for','newpar.for': Updated according to changes in 'fdaux.for'. 'modfd.for': Updated according to changes in 'model.for'. 'fdmax.pl: Updated according to changes in 'grdnorm.for'. 'snap.cal': Confused div and rot corrected. 'fd.h': Snapshots supplemented with the receivers. 'modfdps.h': Colours of complex blocks changed. *** new *** 'addrec.cal': Script for 'grdcal.for' to overlay snapshots with receiver positions. 5.40 (2000, May): 'modfd.for': Problems with gridpoints at interfaces fixed. 'fdaux.for': Erroneous averaging along the X1 axis fixed. 5.50 (2001, June): 'fdmax.pl': Missing apostrophes at data filenames fixed. 'fdsnaps.htm': Mistype in a comment line corrected.
# # Perl script to compile package FD by means of Perl script 'f.pl' # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ require 'f.pl'; # &COMPILE('modfd'); &COMPILE('fdmod'); &COMPILE('fd2d'); &COMPILE('newpar'); # 1; #gmtra.for 0100666 0000766 0000766 00000002241 06477135146 012267 0 ustar klimes klimes C SUBROUTINE 'GMTRA' OF THE IBM SCIENTIFIC SUBROUTINE PACKAGE. C C NOTE: TO CONFORM WITH THE FORTRAN77 STANDARD, DUMMY ARRAY DIMENSIONS C (1) HAVE BEEN CHANGED TO (*). C C .................................................................. C C SUBROUTINE GMTRA C C PURPOSE C TRANSPOSE A GENERAL MATRIX C C USAGE C CALL GMTRA(A,R,N,M) C C DESCRIPTION OF PARAMETERS C A - NAME OF MATRIX TO BE TRANSPOSED C R - NAME OF RESULTANT MATRIX C N - NUMBER OF ROWS OF A AND COLUMNS OF R C M - NUMBER OF COLUMNS OF A AND ROWS OF R C C REMARKS C MATRIX R CANNOT BE IN THE SAME LOCATION AS MATRIX A C MATRICES A AND R MUST BE STORED AS GENERAL MATRICES C C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED C NONE C C METHOD C TRANSPOSE N BY M MATRIX A TO FORM M BY N MATRIX R C C .................................................................. C SUBROUTINE GMTRA(A,R,N,M) DIMENSION A(*),R(*) C IR=0 DO 10 I=1,N IJ=I-N DO 10 J=1,M IJ=IJ+N IR=IR+1 10 R(IR)=A(IJ) RETURN END modfd.dat 0100666 0000766 0000766 00000004542 06532432366 012232 0 ustar klimes klimes 'Demo model for 2-D elastic finite differences' / 0 110 0 40 0 1 (boundaries of the model) 6 SURFACES 6 SIMPLE BLOCKS: -1 2 3 / -1 -2 3 4 / -3 4 5 6 / -1 -4 5 / -5 6 / -6 / 3 COMPLEX BLOCKS: 1 / 2 3 / 4 5 6 / 'EARTH SURFACE' 1 -3 0 0 0 / (i.e. W-X3=0, tension=0) 0 (W) 'SURFACE' 2 1 -3 0 0 / (i.e. W(X1)-X3=0, tension=0) 2 (numbers of grid points) 20 60 (X1 grid coordinates) 0 5 (X3 grid coordinates) 'SURFACE' 3 -3 0 0 0 / (i.e. W-X3=0, tension=0) 5 (W) 'SURFACE' 4 1 -3 0 0 / (i.e. W(X1)-X3=0, tension=0) 2 (numbers of grid points) 10 30 (X1 grid coordinates) 0 15 (X3 grid coordinates) 'SURFACE' 5 1 -3 0 0 / (i.e. W(X1)-X3=0, tension=0) 2 (numbers of grid points) 30 50 (X1 grid coordinates) 15 20 (X3 grid coordinates) 'SURFACE' 6 -3 0 0 0 / (i.e. W-X3=0, tension=0) 20 (W) 'END OF SURFACES' / 'COMPLEX BLOCK' 1 'VP ' 1 (P wave velocity) 0 0 0 0 / (i.e. VP=constant, tension=0) 200 (the value of velocity) 'VS ' 1 (S wave velocity) 0 0 0 0 / (i.e. VS=constant, tension=0) 100 (the value of velocity) 'DEN' 1 (density) 0 0 0 0 / (i.e. density=constant, tension=0) 1600 (the value of density) 'COMPLEX BLOCK' 2 'VP ' 1 (P wave velocity) 0 0 0 0 / (i.e. VP=constant, tension=0) 300 (the value of velocity) 'VS ' 1 (S wave velocity) 0 0 0 0 / (i.e. VS=constant, tension=0) 150 (the value of velocity) 'DEN' 1 (density) 0 0 0 0 / (i.e. density=constant, tension=0) 1700 (the value of density) 'COMPLEX BLOCK' 3 'VP ' 1 (P wave velocity) 0 0 0 0 / (i.e. VP=constant, tension=0) 700 (the value of velocity) 'VS ' 1 (S wave velocity) 0 0 0 0 / (i.e. VS=constant, tension=0) 350 (the value of velocity) 'DEN' 1 (density) 0 0 0 0 / (i.e. density=constant, tension=0) 1800 (the value of density) 'END OF COMPLEX BLOCKS, END OF THE INPUT DATA FOR THE MODEL' / modfd.for 0100666 0000766 0000766 00000075665 07073775642 012276 0 ustar klimes klimes C
C Program MODFD to calculate effective elastic parameters (harmonic C averages of elastic parameters) for 2-D elastic finite differences C in models specified by means of the MODEL package C C Version: 5.40 C Date: 2000, April 9 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 Contents: C Description of data files C Description of effective elast. parameters C List of external procedures C C....................................................................... C C Description of data files: C C Main input data read from the * input device: C The data are read in by the list directed input (free format), by C a single READ statement. The data consist of character strings C which have to be enclosed in apostrophes. Note that the C interactive * external unit may be piped or redirected to a file. C The following items are read: C (1) 'FDSEP' / C 'FDSEP'... String in apostrophes containing the name of the main C input data file with FD parameters in the SEP format. C The file specifies the filename of the data for the model, C dimensions of the grid of N1*N2*N3 gridpoints, the form C and names of the files with effective elastic parameters, C etc. C /... It is recommended to terminate these data by a slash. C Default: 'FDSEP'='fd.h'. C C Data file 'FDSEP' has the form of the SEP (Stanford Exploration C Project) 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 the model by means of the MODEL package: C MODEL='string'... Name of the input file with the data specifying C the model. Must be specified. C Description of MODEL C Example of data set MODEL C Default: MODEL=' ' C This program understands the following parameters common to nearly all C programs dealing with SEP-like gridded data formats: 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. C Note that FD program 'fd2d.for' assumes the free Earth C surface at vertical coordinate X3=O3. C 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 The grid intervals may also be negative. C And the following parameters specific to effective-parameter FD codes: C FORM='string'... Form of the output files: C FORM='formatted': Formatted ASCII output files, C FORM='unformatted': Unformatted output files. C A11='string', B11='string', C11='string'... Strings in apostrophes C containing the names of the output files with (N1-1)*N2*N3 C values of the elastic parameters averaged in the direction C the X1 axis, along the gridlines. C Defaults: A11=' ', B11=' ', C11=' '. C A33='string', B33='string', C33='string'... Strings in apostrophes C containing the names of the output files with N1*N2*(N3-1) C values of the elastic parameters averaged in the direction C the X3 axis, along the gridlines. C Defaults: A33=' ', B33=' ', C33=' '. C A13='string', B13='string', C13='string'... Strings in apostrophes C containing the names of the output files with C (N1-1)*N2*(N3-1) values of the elastic parameters averaged C in the direction the X1 axis, along the gridlines shifted C half a grid interval in the direction of the X3 axis. C Defaults: A13=' ', B13=' ', C13=' '. C A31='string', B31='string', C31='string'... Strings in apostrophes C containing the names of the output files with C (N1-1)*N2*(N3-1) values of the elastic parameters averaged C in the direction the X3 axis, along the gridlines shifted C half a grid interval in the direction of the X1 axis. C Defaults: A31=' ', B31=' ', C31=' '. C Analogously A22, B22, C22, A12, B12, C12, A21, B21, C21, A23, B23, C C23, A32, B32 and C32. C DEN='string'... String in apostrophes containing the name of the C output ASCII file with N1*N2*N3 grid values of the C density. C Default: DEN=' '. C QFC='string'... String in apostrophes containing the name of the C output ASCII file with N1*N2*N3 grid values of the P-wave C quality factor. C Default: QFC=' '. C VP='string'... String in apostrophes containing the name of the C output ASCII file with N1*N2*N3 grid values of the P wave C velocity. C Default: VP=' '. C VS='string'... String in apostrophes containing the name of the C output ASCII file with N1*N2*N3 grid values of the S wave C velocity. C Default: VS=' '. C Example of data set FDSEP C C....................................................................... C C Description of effective elastic parameters: C C Harmonic averages of the following elastic parameters are evaluated: C A = 1 / average( 1/(lambda+2*mu) ) C B = 1 / average( 1/mu ) C C = 1 / average( 1/lambda ) C Elastic parameters A,B,C are averaged between gridpoints along legs C in direction X1 (yielding parameters A1j,B1j,C1j), direction X2 C (yielding parameters A2j,B2j,C2j), and direction X3 (yielding C parameters A3j,B3j,C3j). Parameters Aij,Bij,Cij correspond to the C grid legs between gridpoints if i=j. For j different from i, the C interlegs are considered. The interlegs Aij,Bij,Cij are centred C between the grid legs Aii,Bii,Cii with respect to the Xj axis. C For example, A11,B11,C11 correspond to the grid legs in the direction C of the X1 axis, whereas A13,B13,C13 correspond to the interlegs in the C direction of the X1 axis, shifted vertically with respect to the C former grid legs if the X3 axis is vertical. C C The density is discretized at gridpoints of the basic grid N1*N2*N3: C DEN = density C C The P-wave quality factor is discretized at gridpoints of the basic C grid N1*N2*N3 to be substituted for the finite-difference quantity C QFC = Q/f, C where Q is the quality factor and f is frequency. C C For the legs or gridpoints situated at the structural interfaces, the C parameters are calculated independently at both the sides of the C interface. The results are then averaged arithmetically. C Please, note that this is the difference with respect to program FDMOD C by Jiri Zahradnik. The arithmetic averages should be calculated first C and harmonically averaged along grid legs. Moreover, when calculating C the density at the edges, the slopes of interfaces should be taken C into account instead of averaging the values from different blocks C with the same weight. C C The number of output values of each average elastic parameter is thus: C for A11,B11,C11: (N1-1)* N2 * N3 values, C for A12,B12,C12,A21,B21,C21: (N1-1)*(N2-1)* N3 values, C for A22,B22,C22: N1 *(N2-1)* N3 values, C for A13,B13,C13,A31,B31,C31: (N1-1)* N2 *(N3-1) values, C for A23,B23,C23,A32,B32,C32: N1 *(N2-1)*(N3-1) values, C for A33,B33,C33: N1 * N2 *(N3-1) values, C for DEN,QFC: N1 * N2 * N3 values, C where the dimensions of the basic grid are N1*N2*N3. C C....................................................................... C C This file consists of the following external procedures: C MODFD...Main program. C MODFD C LEGS... Subroutine to evaluate reciprocal values of elastic C parameters A=lambda+2*mu, B=mu and C=lambda along C the legs of a given grid. C LEGS C ELPAR...Subroutine to evaluate reciprocal values of elastic C parameters A=lambda+2*mu, B=mu and C=lambda at a C given point. C ELPAR C C======================================================================= C C C PROGRAM MODFD C C....................................................................... C C Input and output data files: CHARACTER*80 FSEP,FMOD CHARACTER*80 FA11,FB11,FC11,FA12,FB12,FC12,FA13,FB13,FC13 CHARACTER*80 FA21,FB21,FC21,FA22,FB22,FC22,FA23,FB23,FC23 CHARACTER*80 FA31,FB31,FC31,FA32,FB32,FC32,FA33,FB33,FC33 CHARACTER*80 FDEN,FQFC,FVP,FVS INTEGER LU PARAMETER (LU=1) C C Input data and parameters: INTEGER N1,N2,N3 REAL O1,O2,O3,D1,D2,D3,S1,S2,S3 C C Auxiliary storage locations: INTEGER IOUT,NOUT,NALL C C....................................................................... C C Reading main input data: FSEP='fd.h' WRITE(*,'(A)') ' MODFD: Enter names of SEP file [''fd.h'']: ' READ (*,*) FSEP WRITE(*,'(A)') '+MODFD: Working... ' 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',D1,1.) CALL RSEP3R('D2',D2,1.) CALL RSEP3R('D3',D3,1.) C Shifting grid origin by half the grid intervals for interlegs S1=O1+0.5*D1 S2=O2+0.5*D2 S3=O3+0.5*D3 C C Reading the data for the model: CALL RSEP3T('MODEL',FMOD,' ') OPEN(LU,FILE=FMOD,STATUS='OLD') CALL MODEL1(LU) CLOSE(LU) C C Recalling the output filenames: CALL RSEP3T('A11',FA11,' ') CALL RSEP3T('B11',FB11,' ') CALL RSEP3T('C11',FC11,' ') CALL RSEP3T('A12',FA12,' ') CALL RSEP3T('B12',FB12,' ') CALL RSEP3T('C12',FC12,' ') CALL RSEP3T('A13',FA13,' ') CALL RSEP3T('B13',FB13,' ') CALL RSEP3T('C13',FC13,' ') CALL RSEP3T('A21',FA21,' ') CALL RSEP3T('B21',FB21,' ') CALL RSEP3T('C21',FC21,' ') CALL RSEP3T('A22',FA22,' ') CALL RSEP3T('B22',FB22,' ') CALL RSEP3T('C22',FC22,' ') CALL RSEP3T('A23',FA23,' ') CALL RSEP3T('B23',FB23,' ') CALL RSEP3T('C23',FC23,' ') CALL RSEP3T('A31',FA31,' ') CALL RSEP3T('B31',FB31,' ') CALL RSEP3T('C31',FC31,' ') CALL RSEP3T('A32',FA32,' ') CALL RSEP3T('B32',FB32,' ') CALL RSEP3T('C32',FC32,' ') CALL RSEP3T('A33',FA33,' ') CALL RSEP3T('B33',FB33,' ') CALL RSEP3T('C33',FC33,' ') CALL RSEP3T('DEN',FDEN,' ') CALL RSEP3T('QFC',FQFC,' ') CALL RSEP3T('VP' ,FVP ,' ') CALL RSEP3T('VS' ,FVS ,' ') C C....................................................................... C C Calculating elastic parameters: C C Preparing screen output with the progress report: C Percents of calculated values IOUT=0 C Number of calculated values NOUT=0 C Number of values to be calculated NALL=0 IF(FA11.NE.' '.OR.FB11.NE.' '.OR.FC11.NE.' ' * .OR.FDEN.NE.' '.OR.FQFC.NE.' ' * .OR.FVP .NE.' '.OR.FVS .NE.' ') THEN NALL=NALL+(N1-1)* N2 * N3 END IF IF(FA12.NE.' '.OR.FB12.NE.' '.OR.FC12.NE.' ') THEN NALL=NALL+(N1-1)*(N2-1)* N3 END IF IF(FA13.NE.' '.OR.FB13.NE.' '.OR.FC13.NE.' ') THEN NALL=NALL+(N1-1)* N2 *(N3-1) END IF IF(FA21.NE.' '.OR.FB21.NE.' '.OR.FC21.NE.' ') THEN NALL=NALL+(N1-1)*(N2-1)* N3 END IF IF(FA22.NE.' '.OR.FB22.NE.' '.OR.FC22.NE.' ') THEN NALL=NALL+ N1 *(N2-1)* N3 END IF IF(FA23.NE.' '.OR.FB23.NE.' '.OR.FC23.NE.' ') THEN NALL=NALL+ N1 *(N2-1)*(N3-1) END IF IF(FA31.NE.' '.OR.FB31.NE.' '.OR.FC31.NE.' ') THEN NALL=NALL+(N1-1)* N2 *(N3-1) END IF IF(FA32.NE.' '.OR.FB32.NE.' '.OR.FC32.NE.' ') THEN NALL=NALL+ N1 *(N2-1)*(N3-1) END IF IF(FA33.NE.' '.OR.FB33.NE.' '.OR.FC33.NE.' ') THEN NALL=NALL+ N1 * N2 *(N3-1) END IF NALL=NALL*4 C C Initial screen output: WRITE(*,'(A,I3,A)') '+',IOUT, * ' per cent evaluated ' C C Calculating elastic parameters along grid legs in direction X1: CALL LEGS(LU,FA11,FB11,FC11,FDEN,FQFC,FVP,FVS, * 1,N1,O1,D1,2,N2 ,O2,D2,3,N3 ,O3,D3,IOUT,NOUT,NALL) CALL LEGS(LU,FA12,FB12,FC12,' ',' ',' ',' ', * 1,N1,O1,D1,2,N2-1,S2,D2,3,N3 ,O3,D3,IOUT,NOUT,NALL) CALL LEGS(LU,FA13,FB13,FC13,' ',' ',' ',' ', * 1,N1,O1,D1,2,N2 ,O2,D2,3,N3-1,S3,D3,IOUT,NOUT,NALL) C C Calculating elastic parameters along grid legs in direction X2: CALL LEGS(LU,FA21,FB21,FC21,' ',' ',' ',' ', * 2,N2,O2,D2,1,N1-1,S1,D1,3,N3 ,O3,D3,IOUT,NOUT,NALL) CALL LEGS(LU,FA22,FB22,FC22,' ',' ',' ',' ', * 2,N2,O2,D2,1,N1 ,O1,D1,3,N3 ,O3,D3,IOUT,NOUT,NALL) CALL LEGS(LU,FA23,FB23,FC23,' ',' ',' ',' ', * 2,N2,O2,D2,1,N1 ,O1,D1,3,N3-1,S3,D3,IOUT,NOUT,NALL) C C Calculating elastic parameters along grid legs in direction X3: CALL LEGS(LU,FA31,FB31,FC31,' ',' ',' ',' ', * 3,N3,O3,D3,1,N1-1,S1,D1,2,N2 ,O2,D2,IOUT,NOUT,NALL) CALL LEGS(LU,FA32,FB12,FC32,' ',' ',' ',' ', * 3,N3,O3,D3,1,N1 ,O1,D1,2,N2-1,S2,D2,IOUT,NOUT,NALL) CALL LEGS(LU,FA33,FB33,FC33,' ',' ',' ',' ', * 3,N3,O3,D3,1,N1 ,O1,D1,2,N2 ,O2,D2,IOUT,NOUT,NALL) C WRITE(*,'(A)') '+MODFD: Done. ' STOP END C C======================================================================= C C C SUBROUTINE LEGS(LU,FA,FB,FC,FDEN,FQFC,FVP,FVS, * K1,N1,O1,D1,K2,N2,O2,D2,K3,N3,O3,D3,IOUT,NOUT,NALL) C CHARACTER*(*) FA,FB,FC,FDEN,FQFC,FVP,FVS INTEGER LU,K1,N1,K2,N2,K3,N3,IOUT,NOUT,NALL REAL O1,D1,O2,D2,O3,D3 C C Subroutine to evaluate reciprocal values of elastic parameters C A=lambda+2*mu, B=mu and C=lambda along the legs of a given grid. C C Input: C LU... Logical unit number to be used for the output. C FA,FB,FC,FDEN,FQFC... Output file names. File of blank names will C not be generated. FDEN and FQFC may be non-blank only if C K1=1, K2=2, K3=3. C K1... Index of the coordinate axis parallel with the grid legs. C K2,K3...Indices of the other two axes in an ascending order. C N1... Number of grid legs along a gridline increased by 1 (i.e., C the number of grid points in the corresponding direction). C N2,N3...Numbers of gridlines. C O1,O2,O3... K1,K2,K3-th coordinates of the first gridpoint of the C first gridline. C D1,D2,D3... Grid intervals along the K1,K2,K3-th coordinate axes. C IOUT,NOUT,NALL... Integers to control the screen output with the C progress report. C IOUT: Percents of calculated values. C NOUT: Number of calculated values. C NALL: Number of values to be calculated. C C Output: C IOUT,NOUT,NALL... Updated input values. C C Date: 2000, April 9 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Common block /RAMC/ to allocate large array RAM: INCLUDE 'ram.inc' C ram.inc C C----------------------------------------------------------------------- C C Input data and parameters: CHARACTER*80 FORM REAL ERR,W1,W2,W3 C C Averaged parameters: LOGICAL LA,LB,LC REAL AOLD,BOLD,COLD,ANEW,BNEW,CNEW,A,B,C,DEN,QFC,VP,VS C C Auxiliary storage locations: INTEGER NN1,N1N2N3,J1,J2,J3,I1,I2,I3,IW2,IW3,IA,IB,IC,ID,IE,IF,IG INTEGER IH,ITER C Indices of blocks, surfaces, etc.: INTEGER IY(8),ISBOLD,ICBOLD,ISBNEW,ICBNEW,ISRF,IDUMMY(1) C Coordinates: REAL XOLD(3),XNEW(3),XTMP(3),X(3),XOLD1,XNEW1,XTMP1,X1,X2,X3 C Directional vectors etc.: REAL DOLD(3),DNEW(3),DTMP(3),D(3),DIST,H,DUMMY(1) C C....................................................................... C C Upper error bound for determination of the interface position: ERR=0.001000*ABS(D1) W1=AMAX1(ABS(O1),ABS(O1+FLOAT(N1-1)*D1)) W2=AMAX1(ABS(O2),ABS(O2+FLOAT(N2-1)*D2)) W3=AMAX1(ABS(O3),ABS(O3+FLOAT(N3-1)*D3)) W1=AMAX1(0.0000001*W1,0.000001*ABS(D1)) W2=AMAX1(0.0000001*W2,0.000001*ABS(D2)) W3=AMAX1(0.0000001*W3,0.000001*ABS(D3)) W1=SIGN(W1,D1) C C Check whether the output files should be generated: IF(FA.EQ.' '.AND.FB .EQ.' '.AND.FC .EQ.' ' * .AND.FDEN.EQ.' '.AND.FQFC.EQ.' ' * .AND.FVP .EQ.' '.AND.FVS .EQ.' ') THEN RETURN END IF NN1=N1-1 N1N2N3=NN1*N2*N3 C C Increments of indices in array RAM: IF(K1.EQ.1.AND.K2.EQ.2.AND.K3.EQ.3) THEN J1=1 J2=NN1 J3=0 ELSE IF(K1.EQ.2.AND.K2.EQ.1.AND.K3.EQ.3) THEN J1=N2 J2=1 J3=NN1*N2-N2 ELSE IF(K1.EQ.3.AND.K2.EQ.1.AND.K3.EQ.2) THEN J1=N2*N3 J2=1 J3=0 ELSE RETURN END IF C C....................................................................... C C Calculating elastic parameters along grid legs in direction K1: IA=1 IB=IA+N1N2N3 IC=IB+N1N2N3 ID=IC+N1N2N3 IE=ID+N1*N2*N3 IF=IE+N1*N2*N3 IF(FVP.EQ.' ') THEN IG=IF+1 ELSE IG=IF+N1*N2*N3 END IF IF(FVS.EQ.' ') THEN IH=IG+1 ELSE IH=IG+N1*N2*N3 END IF IF(IH.GT.MRAM) THEN C MODFD-01 CALL ERROR('MODFD-01: Insufficient memory') C Increase the dimension MRAM in include file RAM.INC END IF ISBOLD=0 DOLD(1)=0. DOLD(2)=0. DOLD(3)=0. DOLD(K1)=1. DNEW(1)=0. DNEW(2)=0. DNEW(3)=0. DNEW(K1)=0. DTMP(1)=0. DTMP(2)=0. DTMP(3)=0. DTMP(K1)=0. C DO 33 I3=0,N3-1 X3=O3+D3*FLOAT(I3) DO 32 I2=0,N2-1 X2=O2+D2*FLOAT(I2) DO 23 IW3=-1,1,2 XNEW(K3)=X3+W3*FLOAT(IW3) XOLD(K3)=XNEW(K3) XTMP(K3)=XNEW(K3) DO 22 IW2=-1,1,2 XNEW(K2)=X2+W2*FLOAT(IW2) XOLD(K2)=XNEW(K2) XTMP(K2)=XNEW(K2) C DO 21 I1=0,NN1 XNEW(K1)=O1+D1*FLOAT(I1)-W1 XOLD(K1)=O1+D1*FLOAT(I1-1)+W1 IF(IW2.EQ.-1.AND.IW3.EQ.-1) THEN RAM(IA)=0. RAM(IB)=0. RAM(IC)=0. RAM(ID)=0. RAM(IE)=0. IF(FVP.NE.' ') THEN RAM(IF)=0. END IF IF(FVS.NE.' ') THEN RAM(IG)=0. END IF END IF C C Elastic parameters A,B,C are harmonically averaged along C four grid legs close to the leg between gridpoints C XOLD and XNEW. The resulting elastic parameters are the C arithmetic averages of the respective four harmonic C averages. C C Initial values for averaging: A=0. B=0. C=0. H=0. C C Evaluation of block indices and material parameters: CALL BLOCK(XNEW,0,0,ISRF,ISBNEW,ICBNEW) CALL ELPAR * (ICBNEW,XNEW,LA,LB,LC,ANEW,BNEW,CNEW,DEN,QFC,VP,VS) C C Density and P-wave quality factor at the gridpoint: RAM(ID)=RAM(ID)+0.125*DEN RAM(IE)=RAM(IE)+0.125*QFC IF(FVP.NE.' ') THEN RAM(IF)=RAM(IF)+0.125*VP END IF IF(FVS.NE.' ') THEN RAM(IG)=RAM(IG)+0.125*VS END IF C C Harmonic average of elastic parameters: IF(I1.GT.0) THEN ITER=0 10 CONTINUE IF(ISBNEW.NE.ISBOLD) THEN ITER=ITER+1 IF(ITER.GT.100) THEN C MODFD-02 CALL ERROR * ('MODFD-02: More than 100 intersections') C More than 100 points of intersection of the C gridline element with the structural interfaces. C Check the input data and then contact the author. END IF C C Determining the interface between points XOLD, XNEW: IY(4)=ISBOLD IY(5)=ICBOLD IY(6)=0 XTMP(K1)=XNEW(K1) XNEW1=XNEW(K1) XOLD1=XOLD(K1) XTMP1=XTMP(K1) CALL CDE(0,0,IDUMMY,0,IDUMMY,DUMMY,1,3,3,IY,ERR, * XOLD1,XOLD1,XOLD,DOLD,XNEW1,XNEW,DNEW, * XTMP1,XTMP,DTMP,X1,X,D) IF(IY(6).EQ.0) THEN C No more structural interface: GO TO 20 END IF C X is the point of intersection with the interface. C C Contributions to the average: DIST=X(K1)-XOLD(K1) A=A+DIST*AOLD B=B+DIST*BOLD C=C+DIST*COLD XOLD(K1)=X(K1) CALL ELPAR * (ICBOLD,XOLD,LA,LB,LC,AOLD,BOLD,COLD,DEN,QFC,VP,VS) A=A+DIST*AOLD B=B+DIST*BOLD C=C+DIST*COLD H=H+DIST ISBOLD=IY(7) ICBOLD=IY(8) CALL ELPAR * (ICBOLD,XOLD,LA,LB,LC,AOLD,BOLD,COLD,DEN,QFC,VP,VS) C Point XOLD is now shifted to the point of C intersection. C GO TO 10 END IF 20 CONTINUE C C Last contributions to the average: DIST=XNEW(K1)-XOLD(K1) A=A+DIST*(AOLD+ANEW) B=B+DIST*(BOLD+BNEW) C=C+DIST*(COLD+CNEW) H=H+DIST C C Storing average elastic parameters in the memory: IF(LA) THEN RAM(IA)=RAM(IA)+0.5*H/A END IF IF(LB) THEN RAM(IB)=RAM(IB)+0.5*H/B END IF IF(LC) THEN RAM(IC)=RAM(IC)+0.5*H/C END IF C C Screen output: NOUT=NOUT+1 IF(100*NOUT.GT.IOUT*NALL) THEN IOUT=IOUT+1 WRITE(*,'(A,I3)') '+',IOUT END IF IA=IA+J1 IB=IB+J1 IC=IC+J1 END IF C C Proceeding to the initial point of the next interval: XOLD(K1)=O1+D1*FLOAT(I1)+W1 LA=.TRUE. LB=.TRUE. LC=.TRUE. C C Evaluation of block indices and material parameters: CALL BLOCK(XOLD,0,0,ISRF,ISBOLD,ICBOLD) CALL ELPAR * (ICBOLD,XOLD,LA,LB,LC,AOLD,BOLD,COLD,DEN,QFC,VP,VS) C C Density and P-wave quality factor at the gridpoint: RAM(ID)=RAM(ID)+0.125*DEN RAM(IE)=RAM(IE)+0.125*QFC IF(FVP.NE.' ') THEN RAM(IF)=RAM(IF)+0.125*VP END IF IF(FVS.NE.' ') THEN RAM(IG)=RAM(IG)+0.125*VS END IF C ID=ID+1 IE=IE+1 IF=IF+1 IG=IG+1 21 CONTINUE IA=IA-J1*NN1 IB=IB-J1*NN1 IC=IC-J1*NN1 ID=ID-N1 IE=IE-N1 IF=IF-N1 IG=IG-N1 22 CONTINUE 23 CONTINUE IA=IA+J2 IB=IB+J2 IC=IC+J2 ID=ID+N1 IE=IE+N1 IF=IF+N1 IG=IG+N1 32 CONTINUE IA=IA+J3 IB=IB+J3 IC=IC+J3 33 CONTINUE C C Form of the output files: CALL RSEP3T('FORM',FORM,'formatted') C C Writing the output: IA=1 IB=IA+N1N2N3 IC=IB+N1N2N3 ID=IC+N1N2N3 IE=ID+N1*N2*N3 IF=IE+N1*N2*N3 IG=IF+N1*N2*N3 IF(FDEN.NE.' ') THEN CALL WARRAY(LU,FDEN,FORM,.TRUE.,0.,.FALSE.,0., * N1*N2*N3,RAM(ID)) END IF IF(FQFC.NE.' ') THEN CALL WARRAY(LU,FQFC,FORM,.TRUE.,0.,.FALSE.,0., * N1*N2*N3,RAM(IE)) END IF IF(FVP .NE.' ') THEN CALL WARRAY(LU,FVP ,FORM,.TRUE.,0.,.FALSE.,0., * N1*N2*N3,RAM(IF)) END IF IF(FVS .NE.' ') THEN CALL WARRAY(LU,FVS ,FORM,.TRUE.,0.,.FALSE.,0., * N1*N2*N3,RAM(IG)) END IF IF(NN1.GT.1) THEN IF(FA.NE.' ') THEN CALL WARRAY(LU,FA,FORM,.TRUE.,0.,.FALSE.,0.,N1N2N3,RAM(IA)) END IF IF(FB.NE.' ') THEN CALL WARRAY(LU,FB,FORM,.TRUE.,0.,.FALSE.,0.,N1N2N3,RAM(IB)) END IF IF(FC.NE.' ') THEN CALL WARRAY(LU,FC,FORM,.TRUE.,0.,.FALSE.,0.,N1N2N3,RAM(IC)) END IF END IF WRITE(*,'(A,I3,A)') '+',IOUT, * ' per cent evaluated ' C RETURN END C C======================================================================= C C C SUBROUTINE ELPAR(ICB,X,LA,LB,LC,AINV,BINV,CINV,DEN,QFC,VP,VS) INTEGER ICB REAL X(3),AINV,BINV,CINV,DEN,QFC,VP,VS LOGICAL LA,LB,LC C C Subroutine to evaluate reciprocal values of elastic parameters C A=lambda+2*mu, B=mu and C=lambda at a given point. C C Input: C ICB... Index of the complex block. C X... Coordinates of a given point. C LA,LB,LC... Defined logical variables. C C Output: C LA... Unchanged input value if A is positive at the given point, C otherwise switched to .FALSE. C LB... Unchanged input value if B is positive at the given point, C otherwise switched to .FALSE. C LC... Unchanged input value if C is positive at the given point, C otherwise switched to .FALSE. C AINV... AINV=1/A if A=lambda+2*mu is positive, otherwise AINV=0. C BINV... BINV=1/B if B=mu is positive, otherwise BINV=0. C CINV... CINV=1/C if C=lambda is positive, otherwise CINV=0. C DEN... Density at the given point. C QFC... P-wave quality factor at the given point if it is finite. C Otherwise, QFC=0. C VP... P-wave velocity at the given point. C VS... S-wave velocity at the given point. C C Date: 1998, September 14 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations for local model parameters: REAL UP(10),US(10),QP,QS,VD(10),QL C C....................................................................... C IF(ICB.EQ.0) THEN C Free space: AINV=0. BINV=0. CINV=0. QL =0. DEN=0. VP =0. VS =0. ELSE C Material complex block: CALL PARM2(ICB,X,UP,US,DEN,QP,QS) CALL VELOC(ICB,UP,US,QP,QS,VP,VS,VD,QL) AINV=VP*VP*DEN BINV=VS*VS*DEN CINV=AINV-2.*BINV END IF IF(AINV.GT.0.) THEN AINV=1./AINV ELSE LA=.FALSE. END IF IF(BINV.GT.0.) THEN BINV=1./BINV ELSE LB=.FALSE. END IF IF(CINV.GT.0.) THEN CINV=1./CINV ELSE LC=.FALSE. END IF IF(QL.GT.0.) THEN QFC=1./QL ELSE QFC=0. END IF C RETURN END C C======================================================================= C INCLUDE 'error.for' C error.for INCLUDE 'fdaux.for' C fdaux.for INCLUDE 'gmtra.for' C gmtra.for INCLUDE 'model.for' C model.for INCLUDE 'metric.for' C metric.for INCLUDE 'srfc.for' C srfc.for INCLUDE 'parm.for' C parm.for INCLUDE 'val.for' C val.for INCLUDE 'fit.for' C fit.for INCLUDE 'means.for' C means.for INCLUDE 'sep.for' C sep.for INCLUDE 'forms.for' C forms.for INCLUDE 'length.for' C length.for C C======================================================================= Cmodfdps.h 0100666 0000766 0000766 00000005052 07114363302 012237 0 ustar klimes klimes # History file to generate plots of elastic parameters for 'fd.h' # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # Input files required: #chk.pl: "forms/" "echo.pl" #chk.pl: "forms/" "append.pl" # den.out,a11.out,b11.out,c11.out,a13.out,b13.out,c13.out,a31.out, # b31.out,c31.out,a33.out,b33.out,c33.out # Output files: # den.ps,m11.ps,m13.ps,m31.ps,m33.ps # Temporary files: # m11.out, m13.out, m31.out, m33.out, grdcal.tmp # Data specifying the model MODEL='modfd.dat' # Indices of complex blocks (just to display the model): YSIGN=-1 VCIRC=-6 VREF=1 CREF=0.333333 # block 1 is green N1=440 N3=160 # dense grid for smooth D1=0.250 D3=0.250 # interfaces O1=0.125 O3=0.125 HSIZE=11.0 # just for paper hardcopies #UNIT='pt' # just for bitmaps ICB='icb.out' grid: 'modfdps.h' / # Gridding block indices GRD='icb.out' PS='fdicb.ps' grdps: 'modfdps.h' / # Plotting block indices SNAPS='fdicb.ps' fdpcx.pl: modfdps.h # conversion to PCX (GhostScript required) # Plotting density discretized for 2-D FD program VREF=1600 CREF=0.666667 VCIRC=-300 VMIN=0 N1=111 N2=41 N3=1 HSIZE=11.1 # Just for paper hardcopies GRD='den.out' PS='den.ps' grdps: 'modfdps.h' / # Plotting effective elastic parameters discretized for 2-D FD program VREF=0 CREF=0.666667 VCIRC=-540000000 VMIN=0 NH=1 echo.pl: "@2=@1/3" ">grdcal.tmp" N1=110 N2=41 N3=1 HSIZE=11.0 # Just for paper hardcopies CAL='grdcal.tmp' GRD1='a11.out' GRD2='m11.out' grdcal: 'modfdps.h' / append.pl: m11.out b11.out append.pl: m11.out c11.out N3=3 GRD='m11.out' PS='m11.ps' grdps: 'modfdps.h' / N1=111 N2=40 N3=1 HSIZE=11.1 # Just for paper hardcopies GRD1='a33.out' GRD2='m33.out' grdcal: 'modfdps.h' / append.pl: m33.out b33.out append.pl: m33.out c33.out N3=3 GRD='m33.out' PS='m33.ps' grdps: 'modfdps.h' / N1=110 N2=40 N3=1 HSIZE=11.0 # Just for paper hardcopies GRD1='a13.out' GRD2='m13.out' grdcal: 'modfdps.h' / GRD1='a31.out' GRD2='m31.out' grdcal: 'modfdps.h' / append.pl: m13.out b13.out append.pl: m13.out c13.out append.pl: m31.out b31.out append.pl: m31.out c31.out N3=3 GRD='m13.out' PS='m13.ps' grdps: 'modfdps.h' / GRD='m31.out' PS='m31.ps' grdps: 'modfdps.h' / newpar.for 0100666 0000766 0000766 00000015441 06723406714 012453 0 ustar klimes klimes C
C Program NEWPAR to modify the effective elastic parameters (harmonic C averages of elastic parameters) for 2-D elastic finite differences C C Date: 1999, May 28 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 Main input data read from the * input device: C The data are read in by the list directed input (free format), by C a single READ statement. The data consist of character strings C which have to be enclosed in apostrophes. Note that the C interactive * external unit may be piped or redirected to a file. C The following items are read: C 'FSEP'..String in apostrophes containing the name of the input C file with the data specifying the dimensions of the basic C grid of N1*N2*N3 gridpoints. C /... Input data line should be terminated by a slash. C Defaults: 'FSEP'='fd.h'. C C Data file 'FSEP' specifying grid dimensions have the form of the SEP C (Stanford Exploration Project) 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 This program understands the following parameters common to nearly all C programs dealing with SEP-like gridded data formats: 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 And the following parameters specific to effective-parameter FD codes: C FORM='string'... Form of the input and output files: C FORM='formatted': Formatted ASCII files, C FORM='unformatted': Unformatted files. C A13='string', B13='string', C13='string'... Strings in apostrophes C containing the names of the input files with C (N1-1)*N2*(N3-1) values of the elastic parameters averaged C in the direction the X1 axis, along the gridlines shifted C half a grid interval in the direction of the X3 axis. C Defaults: A13=' ', B13=' ', C13=' '. C A31='string', B31='string', C31='string'... Strings in apostrophes C containing the names of the input files with C (N1-1)*N2*(N3-1) values of the elastic parameters averaged C in the direction the X3 axis, along the gridlines shifted C half a grid interval in the direction of the X1 axis. C Defaults: A31=' ', B31=' ', C31=' '. C AA13='string', BB13='string', CC13='string'... Strings in c apostrophes containing the names of the output files with C (N1-1)*N2*(N3-1) values of the average values of elastic C parameters A13 and A31, B13 and B31, C13 and C31, C respectively. C Defaults: AA13=' ', BB13=' ', CC13=' '. C Analogously A12, B12, C12, A21, B21, C21, AA12, BB12, CC12, A23, C B23, C23, A32, B32, C32, AA23, BB23 and CC23. C C======================================================================= C C Common block /RAMC/ to allocate large array RAM: INCLUDE 'ram.inc' C ram.inc C C----------------------------------------------------------------------- C C Input and output data filenames and logical unit numbers: CHARACTER*80 FSEP INTEGER LU PARAMETER (LU=1) C FSEP... The name of the input data file with FD parameters in SEP C format. C LU... Logical unit number to be used to read and write data C files. C INTEGER NX,NZ,N,I REAL DX C C....................................................................... C C Main input data: C Default: FSEP='fd.h' C Reading main input data: WRITE (*,'(2A)') * ' Enter name of input SEP file [''fd.h'']: ' READ (*,*) FSEP C C Reading all the data from input SEP parameter file: CALL RSEP1(LU,FSEP) C C Data specifying grid dimensions CALL FDGRID(MRAM/3,MRAM/3,MRAM/3,NX,NZ,DX) N=(NX-1)*(NZ-1) C C Modifying the gridded effective elastic parameters A: CALL FDIN1(LU,NX-1,NZ-1,'A13',RAM ,RAM(2*N+1),0,RAM) CALL FDIN1(LU,NX-1,NZ-1,'A31',RAM(N+1),RAM(2*N+1),0,RAM) DO 10 I=1,N RAM(I)=0.5*(RAM(I)+RAM(N+I)) 10 CONTINUE CALL FDOUT1(LU,NX-1,NZ-1,1.,'AA13',RAM,RAM(2*N+1)) C C Modifying the gridded effective elastic parameters B: CALL FDIN1(LU,NX-1,NZ-1,'B13',RAM ,RAM(2*N+1),0,RAM) CALL FDIN1(LU,NX-1,NZ-1,'B31',RAM(N+1),RAM(2*N+1),0,RAM) DO 20 I=1,N RAM(I)=0.5*(RAM(I)+RAM(N+I)) 20 CONTINUE CALL FDOUT1(LU,NX-1,NZ-1,1.,'BB13',RAM,RAM(2*N+1)) C C Modifying the gridded effective elastic parameters C: CALL FDIN1(LU,NX-1,NZ-1,'C13',RAM ,RAM(2*N+1),0,RAM) CALL FDIN1(LU,NX-1,NZ-1,'C31',RAM(N+1),RAM(2*N+1),0,RAM) DO 30 I=1,N RAM(I)=0.5*(RAM(I)+RAM(N+I)) 30 CONTINUE CALL FDOUT1(LU,NX-1,NZ-1,1.,'CC13',RAM,RAM(2*N+1)) C WRITE(*,'(A)') '+NEWPAR: done. ' 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 INCLUDE 'fdaux.for' C fdaux.for INCLUDE 'gmtra.for' C gmtra.for C C======================================================================= Cnewpar.pl 0100666 0000766 0000766 00000003043 06607027620 012267 0 ustar klimes klimes #!perl #
# # Perl script to modify the effective elastic parameters (harmonic # averages of elastic parameters) for 2-D elastic finite differences # # Version: 5.20 # Date: 1998, October 8 # #======================================================================= # require 'sep.pl'; # # Reading main input data from the command line: $SEP=$ARGV[0]; # # Reading input SEP header file: &RSEP1($SEP); &RSEP3('N1',$N1,1); &RSEP3('N3',$N3,1); &RSEP3('A13' ,$A13 ,' '); &RSEP3('A31' ,$A31 ,' '); &RSEP3('AA13',$AA13,' '); &RSEP3('B13' ,$B13 ,' '); &RSEP3('B31' ,$B31 ,' '); &RSEP3('BB13',$BB13,' '); &RSEP3('C13' ,$C13 ,' '); &RSEP3('C31' ,$C31 ,' '); &RSEP3('CC13',$CC13,' '); # Creating temporary SEP header file for grdcal.for: $n1=$N1-1; $n3=$N3-1; open (LU,'>sep.tmp'); print LU "N1=$n1 N3=$n3\n"; close(LU) || die "Error"; # Creating temporary command file for grdcal.for: open (LU,'>newpar.tmp'); print LU 'A=$1+$2 ',"\n"; print LU '$3=0.5*A',"\n"; print LU 'A=$4+$5 ',"\n"; print LU '$6=0.5*A',"\n"; print LU 'A=$7+$8 ',"\n"; print LU '$9=0.5*A',"\n"; close(LU) || die "Error"; # Running grdcal.for: open (LU,'|grdcal'); print LU "'sep.tmp' 'newpar.tmp' ", "'$A13' '$A31' '$AA13' ", "'$B13' '$B31' '$BB13' ", "'$C13' '$C31' '$CC13' /"; close(LU) || die "Error"; #======================================================================= 1; #snap.cal 0100666 0000766 0000766 00000001312 07115104652 012051 0 ustar klimes klimes Transformation of normalized divergence and rotation into norm and angle ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Input: @1,@2...Divergence and rotation normalized so as to locally preserve maximum amplitudes during the propagation. SAMP... SEP parameter to amplify the rotation before further processing. Output: @3... Angle in radians, measured from divergence towards the rotation. @4... Norm, i.e. SQRT(div**2+(rot*SAMP)**2) Commands for grdcal.for ~~~~~~~~~~~~~~~~~~~~~~~ div=@1 rot=@2*SAMP div2=div*div rot2=rot*rot norm2=div2+rot2 absdiv=ABS(div) absdiv=AMAX1(absdiv,0.000001) div=SIGN(absdiv,div) @3=ATAN2(rot,div) @4=SQRT(norm2) snapdiv.cal 0100666 0000766 0000766 00000001064 07115104644 012561 0 ustar klimes klimes Transformation of 2-D displacement into "normalized" divergence ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Input: @1... Displacement derivative dU1/dX1 @2... Displacement derivative dU3/dX3 $3... P wave velocity $4... Density Output: @5... Divergence of the input wavefield, multiplied by SQRT(density*VP**3) in order to locally preserve maximum amplitudes during the propagation across interfaces. Commands for grdcal.for ~~~~~~~~~~~~~~~~~~~~~~~ div=@1+@2 a=$3*$4 a=SQRT(a) a=$3*a @5=div*a snaprot.cal 0100666 0000766 0000766 00000001114 07115104670 012576 0 ustar klimes klimes Transformation of 2-D displacement into "normalized" rotation ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Input: @1... Displacement derivative dU1/dX3 @2... Displacement derivative dU3/dX1 $3... S wave velocity $4... Density Output: @5... X2 component of the rotation of the input wavefield, multiplied by SQRT(density*VS**3) in order to locally preserve maximum amplitudes during the propagation across interfaces. Commands for grdcal.for ~~~~~~~~~~~~~~~~~~~~~~~ rot=@1-@2 a=$3*$4 a=SQRT(a) a=$3*a @5=rot*a