C Subroutine file 'parm.for' for specification and interpolation of the C material parameters of the model in rectangular grids. C C By Vlastislav Cerveny, Ludek Klimes, Ivan Psencik C C This file consists of the following subroutines: C PARM1...Subroutine reading the input data for the material C parameters of the model. C PARM2...Subroutine evaluating the isotropic material parameters C including their first and second derivatives. C The functions may be embedded: the independent variable C of the function may be another material parameter of the C same complex block foregoing in the input data. C ADD10...Auxiliary subroutine to PARM2 summing 3 arrays of C dimension 10. C LIN10...Auxiliary subroutine to PARM2 evaluating the linear C combination of 3 arrays of dimension 10. C PARM3...Subroutine evaluating the anisotropic material parameters C including their first and second derivatives. C Subroutines PARM1, PARM2, and PARM3, supporting isotropic complete ray C tracing algorithm, anisotropic ray tracing and other seismic modelling C algorithms, only mediate the work of subroutines VAL1, VAL2 and FPOWER C which must be appended. In addition, subroutines CURVN1 (or its C alternative CURVB1), CURV2D (or its alternative CURVBD), SURFB1, C SURFBD, VAL3B1, VAL3BD, VGEN, TERMS, SNHCSH, TRIDEC, TRISOL, DSPLNZ, C INTRVL from the subroutine package 'FITPACK' by Alan Kaylor Cline, C Department of Computer Sciences, University of Texas at Austin, are C used. In the complete ray tracing, this software file 'parm.for' may C be replaced by any user-defined package containing subroutines PARM1 C and PARM2 with the same number, type and meaning of their parameters C as in this file. C C Attention: C The lines denoted by '*V' in the first 2 columns serve for the C evaluation of the variances of the propagation velocities with C respect to the model parameters for the purpose of inversion. C If the model variations with respect to the model parameters are C to be calculated, '*V' in the first 2 columns of these lines C should be replaced by 2 spaces. In such a case, subroutines VAR4 C and VAR5 of the file 'var.for' are called to handle the C variations. C C If model variations are taken into account: C Model variations are assumed to be stored while evaluating the C functions during the invocation of subroutine VAL of file C 'val.for' and subsequent routines of file 'fit.for'. C The variations of P-wave velocity are assumed to be stored in C register 1 of the system VAR*, the variations of S-wave velocity C are assumed to be stored in register 2 of the system VAR*. C Variations of the density and loss (or quality) factors are not C considered, although they may be stored in other registers. C Subroutines VAR4 and VAR5 are called within the subroutine PARM2 C in order to deal with the variations of P and S wave velocities. C C Input data (read in by subroutine PARM1): C These input data define the complex blocks. They are read in by C subroutine PARM1. The number NCB of the complex blocks to be C defined is an input argument of subroutine PARM1. The data are C read in by the list directed input (free format). C (1) NCB-times (i.e. once for each complex block) input data (1A)+(1B): C (1A) TEXTG,ICB C Identification of the complex block. C TEXTG...Any string. Its first 3 characters must differ from C 'VP ', 'VS ', 'DEN', 'QP ', 'QS ', C 'A11', 'A12', 'A22', 'A13', 'A23', 'A33', 'A14', 'A24', C 'A34', 'A44', 'A15', 'A25', 'A35', 'A45', 'A55', 'A16', C 'A26', 'A36', 'A46', 'A56', 'A66', C 'B11', 'B12', 'B22', 'B13', 'B23', 'B33', 'B14', 'B24', C 'B34', 'B44', 'B15', 'B25', 'B35', 'B45', 'B55', 'B16', C 'B26', 'B36', 'B46', 'B56', 'B66', C 'Q11', 'Q12', 'Q22', 'Q13', 'Q23', 'Q33', 'Q14', 'Q24', C 'Q34', 'Q44', 'Q15', 'Q25', 'Q35', 'Q45', 'Q55', 'Q16', C 'Q26', 'Q36', 'Q46', 'Q56', 'Q66', 'END'. C ICB... Index of the complex block. C (1B) Several times 'Input data for one material parameter', see below. C Isotropic complex block: C At least one of velocities 'VP ' and 'VS ' must be specified. C Unspecified isotropic elastic parameters ('VP ', 'VS ', 'QP ', C 'QS ') take their default values. Unspecified anisotropic C elastic parameters default to the isotropic medium. C Anisotropic complex block with given isotropic reference medium: C Isotropic complex block with one to all anisotropic elastic C parameters specified. C Anisotropic complex block: C At least 9 reduced anisotropic elastic parameters 'A11', 'A12', C 'A22', 'A13', 'A23', 'A33', 'A44', 'A55', and 'A66' must be C specified in the anisotropic complex block. Unspecified C anisotropic elastic parameters default to zeros. C (2) TEXTE,AUX C End of data. C TEXTE...String, the first 3 characters of which must be upper-case C 'END'. C AUX... Any number or a slash. C For an example refer to the sample input data for the model. C C Input data for one material parameter: C The data are read in by the list directed input (free format). In C the list of input data below, each numbered paragraph indicates C the beginning of a new input operation (new read statement). If C the first letter of the symbolic name of the input variable is C I-N, the corresponding value in input data must be of the type C INTEGER. Otherwise (except TEXTF), the input parameter is of the C type REAL. C (1) TEXTF,POWER C Physical meaning of the function. C TEXTF...String identifying which material parameter the function C describes. Only the first 3 characters are significant. C The first 3 characters of the string must be: C 'VP ' for P wave velocity, C 'VS ' for S wave velocity, C 'DEN' for density, C 'QP ' for P wave loss factor, C 'QS ' for S wave loss factor. C 'A11', 'A12', 'A22', 'A13', 'A23', 'A33', 'A14', 'A24', C 'A34', 'A44', 'A15', 'A25', 'A35', 'A45', 'A55', 'A16', C 'A26', 'A36', 'A46', 'A56', or 'A66' for reduced (i.e. C divided by the density) anisotropic elastic parameters C (components of the real part of the symmetric 6*6 C stiffness matrix divided by the density). C 'Q11', 'Q12', 'Q22', 'Q13', 'Q23', 'Q33', 'Q14', 'Q24', C 'Q34', 'Q44', 'Q15', 'Q25', 'Q35', 'Q45', 'Q55', 'Q16', C 'Q26', 'Q36', 'Q46', 'Q56', or 'Q66' for reduced (i.e. C divided by the density) imaginary anisotropic elastic C parameters (components of the imaginary part of the C symmetric 6*6 stiffness matrix divided by the density). C All strings must be entered in uppercase. C POWER...The specified function is equal to the POWER-th power of C the material parameter. C Examples: C For P wave velocity TEXTF='VP ' and POWER=1.0, C for P wave slowness TEXTF='VP ' and POWER=-1.0, C for S wave quadratic slowness TEXTF='VS ' and C POWER=-2.0, C for P wave loss factor TEXTF='QP ' and POWER=1.0, C for S wave quality factor TEXTF='QS ' and POWER=-1.0. C Default values of isotropic elastic parameters not specified in C the input data (subroutine PARM2): C Isotropic complex block or anisotropic complex block with C given isotropic reference medium (at least one of VP and C VS given): C P wave velocity: VP=1.73205*VS, C S wave velocity: VS=0.57735*VP, C Anisotropic complex block (neither VP nor VS given): C Default isotropic parameters are expressed in terms of C anisotropic parameters to allow for application of the C isotropic code to anisotropic models. C Average isotropic medium minimizing the squared norm of C the differrence between the christoffel matrices, C averaged over propagation directions: C AP=A11+A22+A33, C AL=A12+A13+A23, C AS=A44+A55+A66, C VP=SQRT( (3.*AP+2.*AL+4.*AS)/15 ), C VS=SQRT( ( AP - AL+3.*AS)/15 ), C This reference isotropic medium is designed especially C to serve as a basis for anisotropic perturbations. From C this point of view, the variations of the these default C isotropic material parameters with respect to the model C parameters, evaluated by means of subroutines VAR*, are C unnecessary. No variations of the above defaults are C thus set in PARM2. C In any case: C Density: RHO=1, C P wave loss factor (if S wave loss factor is specified): C QP=1.333333*QS*(US(1)/UP(1))**2, C S wave loss factor (if P wave loss factor is specified): C QS=0.750000*QP*(UP(1)/US(1))**2, C P and S wave loss factors (if none of them is given): C QP=0.0, QS=0.0. C Default values of anisotropic elastic parameters not specified in C the input data (subroutine PARM3): C Isotropic complex block or anisotropic complex block with C given isotropic reference medium (at least one of VP and C VS given): C Default anisotropic parameters are expressed in terms of C isotropic parameters for the sake of continuity with C isotropic models: C A11=A22=A33= VP*VP, C A12=A13=A23= VP*VP-2*VS*VS, C A44=A55=A66= VS*VS, C A14=A24=A34=A15=A25=A35=A45=A16=A26=A36=A46=A56= 0. C Anisotropic complex block (neither VP nor VS given): C Aij=0. C In any case: C RHO=1, C Qij=0 in version 5.00. C (2) IVAR1,IVAR2,IVAR3,SIGMA,POWERW,/ C The form of the function. C IVAR1,IVAR2,IVAR3... Denote the form of the function. The function C must be of the form C F(X1,X2,X3) = W(A1,A2,A3)-B1-B2-B3 . C X1, X2, X3 are the general coordinates. Each of A1, A2, C A3, B1, B2, B3 must be either: (a) one of general C coordinates X1, X2, X3, (b) another previously defined C function F(X1,X2,X3) of the same complex block, or (c) C must be left out. At most 3 of parameters A1-B3 may be of C kind (a) or (b). Note that IVAR1 controls the type of A1 C and B1, IVAR2 controls the type of A2 and B2, IVAR3 C controls the type of A3 and B3. C For IVAR1.EQ.0: A1, B1 are empty (left out). C For IVAR1.EQ.1: A1=X1, B1 is empty. C For IVAR1.EQ.2: A1=X2, B1 is empty. C For IVAR1.EQ.3: A1=X3, B1 is empty. C For IVAR1.GE.4: A1=F(X1,X2,X3), where F(X1,X2,X3) is C another function of the same complex block defined in C the input data as the (IVAR1-3)-th function of the C complex block. B1 is empty. C Example: C If density=1.7+0.2*VP then the interpolated function C is W(A1,A2,A3)=1.7+0.2*A1 with the independent C variable A1=VP(X1,X2,X3). This is specified by C IVAR1=4, IVAR2=0, IVAR3=0 if VP is the first read in C parameter, by IVAR1=5, IVAR2=0, IVAR3=0 if VP is the C second read in parameter, etc. C The possible alternatives are W(A1,A2,A3)=1.7+0.2*A2 C with A2=VP(X1,X2,X3) specified by IVAR1=0, C ivar2=(4 or 5 or the like), ivar3=0, and C W(A1,A2,A3)=1.7+0.2*A3 with A3=VP(X1,X2,X3) specified C by IVAR1=0, IVAR3=2, IVAR3=(4 or 5 or the like). C For IVAR1.EQ.-1: B1=X1, A1 is empty. C For IVAR1.EQ.-2: B1=X2, A1 is empty. C For IVAR1.EQ.-3: B1=X3, A1 is empty. C For IVAR1.GE.-4: B1=F(X1,X2,X3), where F(X1,X2,X3) is C another function of the same complex block defined in C the input data as the (-IVAR1-3)-th function of the C complex block. A1 is empty. C The meaning of the parameters IVAR2, IVAR3 is similar. C Examples: C IVAR1: IVAR2: IVAR3: the form of the function: C 1 2 3 F(X1,X2,X3)=W(X1,X2,X3) C 3 1 2 F(X1,X2,X3)=W(X3,X1,X2) C 1 2 0 F(X1,X2,X3)=W(X1,X2) C 5 0 0 F(X1,X2,X3)=W(F2(X1,X2,X3)), where C F2(X1,X2,X3) is the second material parameter of the C complex block defined in the input data. Function W is C interpolated by means of splines under tension. C SIGMA...Is the tension factor (its sign is ignored). This value C indicates the curviness desired. If ABS(SIGMA) is nearly C zero (e.g. 0.001), the resulting surface is approximately C the tensor product of cubic splines. If ABS(SIGMA) is C large (e.g. 50.), the resulting surface is approximately C tri-linear. If SIGMA equals zero, tensor products of C cubic splines result. A recommended value for SIGMA is C approximately 1. In absolute value. C POWERW..Given grid values (7) correspond to the POWERW-th power of C interpolated function W. The given grid values (7) are C thus raised to the (1/POWERW)-th power immediately after C reading and then interpolated. C /... Obligatory slash at the end of line for future extensions. C Default: IVAR1=0, IVAR2=0, IVAR3=0, SIGMA=0, POWERW=1. C (3) NX(1),...,NX(NVAR) C The numbers of grid coordinates for the interpolation. C This input is performed if at least one of IVAR1, IVAR2, IVAR3 is C positive. C Each of NX(1),...,NX(NVAR) corresponds to one positive value of C IVAR1, IVAR2, IVAR3 and specifies the number of grid coordinates C corresponding to that independent variable of function W, see (2). C The sign of NX(1),...,NX(NVAR) is ignored. NVAR (.LE.3) is the C number of positive values of the above quantities IVAR1, IVAR2, C IVAR3, i.e. the number of independent variables of function W, C see (1). C (4) X1(1),...,X1(NX(1)) C The grid coordinates corresponding to the first independent C variable of function W, see (2). C This input is performed if NX(1) is specified, see (3), and is not C zero. The grid coordinates may be specified in any order. C (5) X2(1),...,X2(NX(2)) C The grid coordinates corresponding to the second independent C variable of function W, see (2). C This input is performed if NX(2) is specified, see (3), and is not C zero. The grid coordinates may be specified in any order. C (6) X3(1),...,X3(NX(3)) C The grid coordinates corresponding to the third independent C variable of function W, see (2). C This input is performed if NX(3) is specified, see (3), and is not C zero. The grid coordinates may be specified in any order. C (7) (((W(I,J,K),I=1,MAX(NX(1),1)),J=1,MAX(NX(2),1)),K=1,MAX(NX(3),1)) C the values of function W at grid points. Function value W(i,j,k) C corresponds to point (X1(I),X2(J),X3(K)). C C Date: 1996, September 30 C Coded by Ludek Klimes C C======================================================================= C SUBROUTINE PARM1(LUN,NCB) INTEGER LUN,NCB C C This subroutine reads the input data for the distributions of the C material parameters: C In the isotropic complex blocks: C P and S wave velocities, density, P and S wave loss factors C in the anisotropic complex blocks: C 21 real parts of the reduced (divided by the density) elastic C parameters, 21 corresponding imaginary parts, density. C Reference isotropic P ans S velocities and loss factors may also be C specified together with anisotropic material parameters. C When reading the data, subroutine PARM1 also determines the parameters C necessary to compute an interpolatory function on a three dimensional C rectangular grid, and stores them in the memory. The function C determined can be represented as a tensor product of splines under C tension. The functions may be embedded. For actual interpolation of C material parameters it is necessary to call subroutine PARM2 for C isotropic model (or mean isotropic model corresponding to the C anisotropic model), or subroutine PARM3 for anisotropic material C parameters of the isotropic or anisotropic model. Subroutines PARM2 C and PARM3 also return the first and second partial derivatives of C propagation velocities or reduced elastic parameters. Subroutine C PARM1 may be called several times. The complex blocks are indexed C succesively, following the complex blocks defined during the previous C invocations. C C Input: C LUN... Logical unit number of the external input device C containing the input data. C NCB... Number of the material complex blocks for which the input C data are specified during the current invocation of PARM1. C None of the input parameters are altered. C C No output. C C Subroutines and external functions required: EXTERNAL VAL1 C VAL1, SORTV, READV... File 'val.for'. C CURVN1 or CURVB1 (alternatives), SURFB1, VAL3B1, SNHCSH, VGEN, C TERMS, TRIDEC, TRISOL... Subroutine package 'FITPACK' C (file 'fit.for'). C C Date: 1995, December 17 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: CHARACTER*3 TFUNCT(47) DATA TFUNCT/'VP ','VS ','DEN','QP ','QS ', *'A11','A12','A22','A13','A23','A33','A14','A24','A34','A44', *'A15','A25','A35','A45','A55','A16','A26','A36','A46','A56','A66', *'Q11','Q12','Q22','Q13','Q23','Q33','Q14','Q24','Q34','Q44', *'Q15','Q25','Q35','Q45','Q55','Q16','Q26','Q36','Q46','Q56','Q66'/ C CALL VAL1(LUN,2,NCB,47,TFUNCT) RETURN END C C======================================================================= C SUBROUTINE PARM2(ICB,COOR,UP,US,RHO,QP,QS) INTEGER ICB REAL COOR(3),UP(10),US(10),RHO,QP,QS C C This subroutine evaluates P and S wave velocities, density, and P and C S wave loss factors at a given point. The three first and six second C partial derivatives of the velocities are also evaluated. The C specified functions are represented as a tensor product of splines C under tension. The parameters may be dependent either on the general C coordinates or on the distribution of another parameter, e.g. C VS=0.577*VP or RHO=1.7+0.2*VP, where VP, VS and RHO are P and S C velocities and density. The coefficients of these functions are C prepared in subroutine PARM1, in which the input data concerning the C distribution of individual parameters within each complex block are C read in. The default values of parameters not specified in the input C data are: C P wave velocity: VP=1.73205*VS, C S wave velocity: VS=0.57735*VP, C Density: RHO=1.0, C P wave loss factor (if the S wave loss factor is specified): C QP=1.333333*QS*(US(1)/UP(1))**2, C S wave loss factor (if the P wave loss factor is specified): C QS=0.750000*QP*(UP(1)/US(1))**2, C P and S wave loss factors (if none of them is specified): C QP=0.0, QS=0.0. C Attention: The above default values are reasonable only if arguments C UP and US of this subroutine are velocities, and QP and QS loss C factors. In other words, these default settings are usefull only C if NEXPV=1 and NEXPQ=1 in the input data set model, line (2). C Note that at least one of the velocities must be specified in the C input data. P wave velocity must be positive, other material C parameters must be non-negative. C C Input: C ICB... Index of a complex block. C COOR... Array containing coordinates X1, X2, X3 of the given C point. C None of the input parameters are altered. C C Output: C UP,US...Powers of the P and S wave velocities (the exponent of the C power is NEXPV, see the input data for the model) and C their first and second partial derivatives in order U, U1, C U2, U3, U11, U12, U22, U13, U23, U33, at the given point. C RHO... Density at the given point. C QP,QS...Powers of the P and S wave loss factors (the exponent of C the power is NEXPQ, see the input data for the model) at C the given point. C C Subroutines and external functions required: EXTERNAL FPOWER,VAL2 C FPOWER...File 'model.for'. C VAL2... File 'val.for'. C CURV2D or CURVBD (alternatives), SURFBD, VAL3BD, SNHCSH, DSPLNZ, C INTRVL... Subroutine package 'FITPACK' (file 'fit.for'). C C Error messages: C 321... No velocity is defined: C Neither P nor S wave velocity in the current complex block C is defined in the input data. C 322... Prohibited material parameter: C P wave velocity must be positive, other material C parameters must be non-negative. C C Date: 1996, September 30 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER I REAL FAUX(10,47),POWER(47),POWER1,POWER2,POWER3,POWER4,POWER5 EQUIVALENCE (POWER(1),POWER1),(POWER(2),POWER2),(POWER(3),POWER3) EQUIVALENCE (POWER(4),POWER4),(POWER(5),POWER5) REAL AUX(50) C C Constants: REAL C1,C2,C3,C4 PARAMETER (C1=1./15.,C2=2./15.,C3=3./15.,C4=4./15.) C C....................................................................... C CALL VAL2(2,IABS(ICB),5,COOR,FAUX,POWER) C C Velocities: IF(POWER1.NE.0.) THEN IF(FAUX(1,1).LE.0.) GO TO 9 CALL FPOWER(10,FAUX(1,1),POWER1,UP) CALL VAR5(1,1) IF(POWER2.NE.0.) THEN IF(FAUX(1,2).LT.0.) GO TO 9 CALL FPOWER(10,FAUX(1,2),POWER2,US) CALL VAR5(2,2) ELSE DO 1 I=1,10 US(I)=0.57735*UP(I) 1 CONTINUE CALL VAR4(0,0.57735) CALL VAR5(2,1) END IF ELSE IF(POWER2.NE.0.) THEN IF(FAUX(1,2).LE.0.) GO TO 9 CALL FPOWER(10,FAUX(1,2),POWER2,US) CALL VAR5(2,2) DO 2 I=1,10 UP(I)=1.73205*US(I) 2 CONTINUE CALL VAR4(0,1.73205) CALL VAR5(1,2) ELSE IF(POWER(06).NE.0..AND.POWER(07).NE.0..AND. * POWER(08).NE.0..AND.POWER(09).NE.0..AND. * POWER(10).NE.0..AND.POWER(11).NE.0..AND. * POWER(15).NE.0..AND.POWER(20).NE.0..AND. * POWER(26).NE.0.) THEN C Isotropic reference medium to the anisotropic material: CALL FPOWER(10,FAUX(1,06),.5*POWER(06),AUX(01)) CALL FPOWER(10,FAUX(1,08),.5*POWER(08),AUX(31)) CALL FPOWER(10,FAUX(1,11),.5*POWER(11),AUX(41)) CALL ADD10(AUX(01),AUX(31),AUX(41)) CALL FPOWER(10,FAUX(1,07),.5*POWER(07),AUX(11)) CALL FPOWER(10,FAUX(1,09),.5*POWER(09),AUX(31)) CALL FPOWER(10,FAUX(1,10),.5*POWER(10),AUX(41)) CALL ADD10(AUX(11),AUX(31),AUX(41)) CALL FPOWER(10,FAUX(1,15),.5*POWER(15),AUX(21)) CALL FPOWER(10,FAUX(1,20),.5*POWER(20),AUX(31)) CALL FPOWER(10,FAUX(1,26),.5*POWER(26),AUX(41)) CALL ADD10(AUX(21),AUX(31),AUX(41)) CALL LIN10(C3,AUX(01), C2,AUX(11),C4,AUX(21),AUX(31)) CALL LIN10(C1,AUX(01),-C1,AUX(11),C3,AUX(21),AUX(41)) CALL FPOWER(10,AUX(31),2.,UP) CALL FPOWER(10,AUX(41),2.,US) ELSE PAUSE 'Error 321 in PARM2: No velocity is defined' STOP END IF END IF END IF C C Density: IF(POWER3.NE.0.) THEN IF(FAUX(1,3).LT.0.) GO TO 9 CALL FPOWER(1,FAUX(1,3),POWER3,AUX) RHO=AUX(1) ELSE RHO=1. END IF C C Loss factors: IF(POWER4.NE.0.) THEN IF(FAUX(1,4).LT.0.) GO TO 9 CALL FPOWER(1,FAUX(1,4),POWER4,AUX) QP=AUX(1) IF(POWER5.NE.0.) THEN IF(FAUX(1,5).LT.0.) GO TO 9 CALL FPOWER(1,FAUX(1,5),POWER5,AUX) QS=AUX(1) ELSE IF(US(1).GT.0.) THEN QS=0.750000*QP*(UP(1)/US(1))**2 ELSE QS=0. END IF END IF ELSE IF(POWER5.NE.0.) THEN IF(FAUX(1,5).LT.0.) GO TO 9 CALL FPOWER(1,FAUX(1,5),POWER5,AUX) QS=AUX(1) QP=1.333333*QS*(US(1)/UP(1))**2 ELSE QP=0. QS=0. END IF END IF RETURN C 9 CONTINUE WRITE(*,'('' X='',F9.3,'' Y='', * F9.3,'' Z='',F9.3/'' VP='',F7.3,'' VS='',F7.3,'' RO='', * F7.3,'' QP='',F7.3,'' QS='',F7.3)') COOR,(FAUX(1,I),I=1,5) PAUSE 'Error 322 in PARM2: Prohibited material parameter' STOP END C C----------------------------------------------------------------------- C SUBROUTINE ADD10(A,B,C) REAL A(10),B(10),C(10) C C Auxiliary subroutine to PARM2 summing 3 arrays of dimension 10. C C....................................................................... C INTEGER I C DO 10 I=1,10 A(I)=A(I)+B(I)+C(I) 10 CONTINUE RETURN END C C----------------------------------------------------------------------- C SUBROUTINE LIN10(C1,A1,C2,A2,C3,A3,A) REAL C1,A1(10),C2,A2(10),C3,A3(10),A(10) C C Auxiliary subroutine to PARM2 evaluating linear combination of C 3 arrays of dimension 10. C C....................................................................... C INTEGER I C DO 10 I=1,10 A(I)=C1*A1(I)+C2*A2(I)+C3*A3(I) 10 CONTINUE RETURN END C C======================================================================= C SUBROUTINE PARM3(ICB,COOR,A,RHO,Q) INTEGER ICB REAL COOR(3),A(10,21),RHO,Q(21) C C This subroutine is redundant for isotropic seismic models and codes. C C This subroutine evaluates 21 real parts of the reduced (divided by the C density) elastic parameters, 21 corresponding imaginary parts, and the C density at a given point. The three first and six second partial C derivatives of the 21 real parts of the reduced elastic parameters are C also evaluated. The specified functions are represented as a tensor C product of splines under tension. The parameters may be dependent C either on the general coordinates or on the distribution of another C parameter. The coefficients of these functions are prepared in C subroutine PARM1, in which the input data concerning the distribution C of individual parameters within each complex block are read in. C Variations of real parts of the reduced elastic parameters Aij with C respect to model parameters (evaluated by subroutines VAR*) are C stored in registers 6 to 26. C C Input: C ICB... Index of a complex block. C COOR... Array containing coordinates X1, X2, X3 of the given C point. C None of the input parameters are altered. C C Output: C A... Values, first and second partial derivatives of real C parts of 21 reduced (divided by the density) elastic C parameters. The order of the value, first and second C partial derivatives of each parameter Aij is: C Aij,Aij1,Aij2,Aij3,Aij11,Aij12,Aij22,Aij13,Aij23,Aij33. C The order of parameters (second array index) is: C A11,A12,A22,A13,A23,A33,A14,A24,A34,A44,A15,A25,A35,A45, C A55,A16,A26,A36,A46,A56,A66. C RHO... Density at the given point. C Q... Imaginary parts of 21 reduced elastic parameters at the C given point, ordered as C Q11,Q12,Q22,Q13,Q23,Q33,Q14,Q24,Q34,Q44,Q15,Q25,Q35,Q45, C Q55,Q16,Q26,Q36,Q46,Q56,Q66. C C Subroutines and external functions required: EXTERNAL FPOWER,VAL2 C FPOWER...File 'model.for'. C VAL2... File 'val.for'. C CURV2D or CURVBD (alternatives), SURFBD, VAL3BD, SNHCSH, DSPLNZ, C INTRVL... Subroutine package 'FITPACK' (file 'fit.for'). C C Error messages: C 323... Undefined elastic parameter: C If isotropic P-velocity VP is not specified in the input C data, A11, A22, and A33 must be specified in order to use C this subroutine. C 324... Undefined elastic parameter: C If isotropic s-velocity VS is not specified in the input C data, A44, A55, and A66 must be specified in order to use C this subroutine. C 325... Prohibited elastic parameter: C No component of both real and imaginary parts of the 6*6 C matrix of elastic parameters (stiffness matrix) may be C negative. C*********************************************************************** C n 326... Prohibited elastic parameter: C o c Both real and imaginary parts of the 6*6 matrix of elastic C t o parameters (stiffness matrix) must be positively C d semi-definite, and the 3*3 upper-left minor of its real C e part must be positively definite. The density must be C d positive. C*********************************************************************** C C Date: 1996, September 30 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER KA(21),IA,I,I1,I2 REAL FAUX(10,47),POWER(47),POWER3,AUX(10) EQUIVALENCE (POWER(3),POWER3) C Order of processing the reduced elastic parameters: DATA KA/1,3,6,21,15,10,2,4,5,20,19,14,16,17,18,11,12,13,7,8,9/ C C....................................................................... C CALL VAL2(2,IABS(ICB),47,COOR,FAUX,POWER) C C Density: IF(POWER3.NE.0.) THEN IF(FAUX(1,3).LT.0.) GO TO 99 CALL FPOWER(1,FAUX(1,3),POWER3,AUX) RHO=AUX(1) ELSE RHO=1. END IF C C Real-valued elastic parameters: DO 19 I2=1,21 IA=KA(I2) I=IA+5 IF(POWER(I).NE.0.) THEN IF(FAUX(1,I).LT.0.) GO TO 99 CALL FPOWER(10,FAUX(1,I),POWER(I),A(1,IA)) CALL VAR5(I,I) ELSE IF(I2.LE.3) THEN C Default diagonal stiffnesses A11,A22,A33: IF(POWER(1).NE.0.) THEN C Aij=VP*VP: IF(FAUX(1,1).LT.0.) GO TO 99 CALL FPOWER(10,FAUX(1,1),0.5*POWER(1),A(1,IA)) CALL VAR5(I,1) ELSE IF(POWER(2).NE.0.) THEN C Aij=3.*VS*VS: IF(FAUX(1,2).LT.0.) GO TO 99 CALL FPOWER(10,FAUX(1,2),0.5*POWER(2),A(1,IA)) CALL VAR5(I,2) DO 11 I1=1,10 A(I1,IA)=3.*A(I1,IA) 11 CONTINUE CALL VAR4(0,3.) CALL VAR5(I,I) ELSE PAUSE 'Error 323 in PARM3: Undefined elastic parameter' STOP END IF ELSE IF(I2.LE.6) THEN C Default diagonal stiffnesses A66,A55,A44: IF(POWER(2).NE.0.) THEN C Aij=VS*VS: IF(FAUX(1,2).LT.0.) GO TO 99 CALL FPOWER(10,FAUX(1,2),0.5*POWER(2),A(1,IA)) CALL VAR5(I,2) ELSE IF(POWER(1).NE.0.) THEN C Aij=.333333*VP*VP: IF(FAUX(1,1).LT.0.) GO TO 99 CALL FPOWER(10,FAUX(1,1),0.5*POWER(1),A(1,IA)) CALL VAR5(I,1) DO 12 I1=1,10 A(I1,IA)=.333333*A(I1,IA) 12 CONTINUE CALL VAR4(0,.333333) CALL VAR5(I,I) ELSE PAUSE 'Error 324 in PARM3: Undefined elastic parameter' STOP END IF ELSE IF(I2.LE.9) THEN C Default non-diagonal stiffnesses A12,A13,A23: IF(POWER(1).NE.0.) THEN IF(FAUX(1,1).LT.0.) GO TO 99 IF(POWER(2).NE.0.) THEN C Aij=VP*VP-2*VS*VS: IF(FAUX(1,2).LT.0.) GO TO 99 CALL FPOWER(10,FAUX(1,2),0.5*POWER(1),AUX) CALL VAR5(I,2) CALL VAR4(0,-2.) CALL VAR5(I,I) CALL FPOWER(10,FAUX(1,1),0.5*POWER(1),A(1,IA)) CALL VAR5(I,1) DO 13 I1=1,10 A(I1,IA)=A(I1,IA)-2.*AUX(I1) 13 CONTINUE ELSE C Aij=.333333*VP*VP: CALL FPOWER(10,FAUX(1,1),0.5*POWER(1),A(1,IA)) CALL VAR5(I,1) DO 14 I1=1,10 A(I1,IA)=.333333*A(I1,IA) 14 CONTINUE CALL VAR4(0,.333333) CALL VAR5(I,I) END IF ELSE IF(POWER(2).NE.0.) THEN C Aij=VS*VS: IF(FAUX(1,2).LT.0.) GO TO 99 CALL FPOWER(10,FAUX(1,2),0.5*POWER(2),A(1,IA)) CALL VAR5(I,2) ELSE PAUSE 'Error 324 in PARM3: Undefined elastic parameter' STOP END IF END IF ELSE C Default non-diagonal stiffnesses C A56,A46,A45,A16,A26,A36,A15,A25,A35,A14,A24,A34: DO 18 I1=1,10 A(I1,IA)=0. 18 CONTINUE END IF 19 CONTINUE C C Imaginary parts of the elastic parameters: DO 29 I2=1,21 I=I2+26 IF(POWER(I).NE.0.) THEN IF(FAUX(1,I).LT.0.) GO TO 99 CALL FPOWER(1,FAUX(1,I),POWER(I),Q(I2)) ELSE C Default: Q(I2)=0. END IF 29 CONTINUE C C Check for positive semidefinitness: C *** not coded *** RETURN C 99 CONTINUE WRITE(*,'('' X='',F9.3,'' Y='',F9.3,'' Z='',F9.3, * '' FAUX(1,'',I2,'')='',F7.3)') COOR,I,FAUX(1,I) PAUSE 'Error 325 in PARM3: Prohibited elastic parameter' STOP END C C======================================================================= C