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 http://sw3d.cz/staff/klimes.htm 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 C MODEL C Example of data set C 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======================================================================= C