C Subroutine file 'model.for' to specify a seismic model. C C Version: 7.10 C Date: 2014, February 14 C Coded by Ludek Klimes C C....................................................................... C C This file consists of the following external procedures: C MODEL1..Subroutine designed to read the input data for the C description of the model and to store them in the memory. C Subroutine MODEL1 reads the input data (1)-(8) listed C below itself, then calls subroutine SRFC1 (included in the C subroutine file 'srfc.for') to read the input data (9) for C smooth surfaces, and finally calls subroutine PARM1 C (included in the subroutine file 'parm.for') to read the C input data (10) for the parameters of the medium. C MODEL1 C NSRFC...Integer function returning the number of surfaces covering C structural interfaces. C NSRFC C BLOCK...Subroutine designed to determine the mutual position of a C point and a simple and a complex block. Just calls more C general subroutine BLOCKS which does the job. C BLOCK C BLOCKS..Subroutine designed to determine the mutual position of a C point and a simple and a complex block. A more general C version of subroutine BLOCK, designed especially to work C with curves situated at structural interfaces. C BLOCKS C ISIDE...Auxiliary function to subroutine BLOCKS. C ISIDE C INTERF..Auxiliary subroutine to subroutine BLOCKS. C INTERF C SEPAR...Subroutine determining the surface separating two given C simple blocks. C SEPAR C VELOC...Subroutine transforming the values of a medium parameters C into velocities and loss factors. C VELOC C FPOWER..Subroutine evaluating the value and, possibly, the three C first and six second partial derivatives of a function, C if the value and the three first and six second partial C derivatives of the POWER-th power of the function are C known. Particularly, this is an auxiliary subroutine C to the subroutine VELOC. C FPOWER C C Note: C The lines denoted by '*V' in the first two columns of file C 'model.for' in subroutines VELOC and FPOWER are designed to C calculate the model variations with respect to the model C parameters. C File 'modelv.for', intended for the model inversion, is created C from 'model.for' by replacing each '*V' in the first two columns C by spaces using program 'clean.for'. Subroutines VAR4 and VAR5 C of file 'var.for' may then be called to handle the variations. C C....................................................................... C C C Input data MODEL for the specification of the model: 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 TEXTM), the input parameter is of the C type REAL. C (1) TEXTM C The string containing the name of the model. Only the first 80 C characters of the string are significant. C (2) KOORS,NEXPV,NEXPQ,IVERT,/ C KOORS...Specifies the type of the coordinate system: C KOORS.LE.0: Cartesian coordinates (default). C KOORS.EQ.1: polar spherical coordinates in radians, C (X1,X2,X3)=(colatitude,longitude,radius). C KOORS.GE.2: geographic spherical coordinates in radians, C (X1,X2,X3)=(longitude,latitude,radius). C KOORS.EQ.3: polar spherical coordinates in specified C angular units, C (X1,X2,X3)=(colatitude,longitude,elevation). C The angular units are specified by parameter CIRCLE of C input data (2A). C The elevation is measured from the reference radius C specified by parameter RADIUS of input data (2A). C KOORS.GE.4: geographic spherical coordinates in specified C angular units, C (X1,X2,X3)=(longitude,latitude,elevation). C The angular units are specified by parameter CIRCLE of C input data (2A). C The elevation is measured from the reference radius C specified by parameter RADIUS of input data (2A). C If the coordinate system is right-handed (recommended), C all vectorial products are evaluated using the right-hand C rule, otherwise using the left-hand rule. C KOORS is passed to the subroutines of file 'metric.for' C by means of invocation of subroutine METR1 and presently C represents the only input data for the coordinate system. C Note that possible future additional data for the C coordinate system should be read by subroutine METR1 and C should be contained within input data (2A). C metric.for C NEXPV,NEXPQ... The default values are highly recommended! C Velocities powered to NEXPV and loss factors powered to C NEXPQ are reported by the subroutines evaluating isotropic C material parameters. C For example, unit values of NEXPV and NEXPQ indicate that C velocities and loss factors are the output parameters C of the subroutines evaluating isotropic material C parameters, indices equal -1 indicate reciprocal values of C these quantities, i.e. slownesses and quality factors. C When using the basic submitted version of the subroutine C file 'parm.for', the default values of NEXPV=1, NEXPQ=1 C are highly recommended. Other values make sense only if a C user is submitting his own subroutines to evaluate C isotropic material parameters which, e.g., output the C slowness instead of the velocity. In such a case, C switching NEXPV from 1 to -1 may avoid the modification C of the user's subroutines. C IVERT...Orientation of the vertical axis: C IVERT=0: unknown (default), C IVERT=+1: X1 vertical, pointing upwards, C IVERT=-1: X1 vertical, pointing downwards, C IVERT=+2: X2 vertical, pointing upwards, C IVERT=-2: X2 vertical, pointing downwards, C IVERT=+3: X3 vertical, pointing upwards, C IVERT=-3: X3 vertical, pointing downwards, C Has no influence on the calculations, and need not be C specified. If it is non-zero, it may be considered for C plotting purposes. C /... Obligatory slash at the end of line for future extensions. C Default: KOORS=0, NEXPV=1, NEXPQ=1, IVERT=0. C (2A) Additional data for the coordinate system. C No data (2A) are present if KOORS=0, KOORS=1 or KOORS=2. C The data are read by subroutine METR1. For their description C refer to file C metric.for. C (3) X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX C Boundaries of the model. C (4) NSRFC C Number of smooth surfaces in the model. The surfaces are indexed C sequentially in any order by positive integers ISRFC from 1 to C NSRFC. C It is recommended to define only surfaces covering structural C interfaces (model surfaces) in this data set. Auxiliary surfaces C related to particular source-receiver configurations, numerical C procedures, etc., should preferably be defined in other data sets. C (5) NSB C Number of simple blocks in the model, defined in this data set. C The defined blocks are indexed sequentially by positive integers C ISB from 1 to NSB in the same order as they are specified in data C (6). Intersecting simple blocks are allowed but not recommended. C All material simple blocks in the model must be defined. C Free-space blocks need not be defined in this data set but it is C recommended to define the free-space simple blocks in order to C speed up the determination of a free space during computations. C Defined free-space blocks also enhance the possibility to check C for the model consistency. If the free-space simple blocks are C not, for some reason, defined here, they may be defined in the C additional data file designed just for the consistency check by C program modchk.for C (6) NSB input operations (READ statements): C For each simple block with index ISB, the indices of the surfaces C forming the set F(+) and the indices of the surfaces forming the C set F(-). The indices of surfaces from F(+) must be positive, the C indices of surfaces from F(-) must be indicated by negative signs. C The indices may be specified in an arbitrary order and must be C terminated by a slash. C (7) NCB C Number of material complex blocks in the model. The material C complex blocks are indexed sequentially by positive integers ICB C from 1 to NCB. The free-space complex block is not indexed and C consists of all simple blocks not contained in the material C complex blocks. Space covered by no simple block is also deemed C to be a free space. C (8) NCB input operations (READ statements): C For each material complex block, the indices of material simple C blocks forming the complex block. The indices may be specified C in an arbitrary order and must be terminated by a slash. C Each material simple block must appear exactly once within these C data lines. Simple blocks defined by data (6) but not listed here C are deemed to be free-space simple blocks. C (9) The data specifying NSRFC functions F(X1,X2,X3) describing the C smooth surfaces in the model. The data are read by subroutine C SRFC1. For their description refer to subroutine SRFC1 (included C in the subroutine file 'srfc.for'). C srfc.for: Input data C (10) The data specifying the distribution of parameters of the model C in all NCB material complex blocks. The data are read by C subroutine PARM1. For their description refer to subroutine C PARM1 (included in the subroutine file 'parm.for'). C parm.for: Input data C For an example refer to the sample input data for the model. C Example of data set C MODEL. C C C Input data in the form of C SEP C parameter or history file: C NEGPAR=integer... Flag whether the negative values of material C parameters are allowed: C NEGPAR=0: Negative values of material parameters or zero C P-wave velocity are reported as errors. C NEGPAR=1: Negative values of material parameters or zero C P-wave velocity are not reported as errors. C Default: NEGPAR=0 C C....................................................................... C C Storage in the memory: C The input data (1) to (8) describing the structure of the model C are stored in common blocks /MODELT/ and /MODELC/ defined in the C include file 'model.inc'. C model.inc C C======================================================================= C C C SUBROUTINE MODEL1(LUN) INTEGER LUN C C Subroutine MODEL1 reads the input data (1)-(8) for the description C of the model and stores them in common blocks /MODELT/ and /MODELC/. C Then it calls subroutine SRFC1 (included in the subroutine file C 'srfc.for') to read the input data (9) for smooth surfaces and to C store them in the memory. Finally, it calls subroutine PARM1 C (included in the subroutine file 'parm.for') to read the input data C (10) for the parameters of the medium and to store them in the memory. C C Input: C LUN... Logical unit number of the external input device C containing the input data. C The input parameter is not altered. C C No output. C C Common blocks /MODELT/ and /MODELC/: INCLUDE 'model.inc' C model.inc C All the storage locations of the common blocks are defined in this C subroutine. C C Subroutines and external functions required: EXTERNAL RSEP3I,METR1,SRFC1,PARM1 C RSEP3I...File 'sep.for'. C METR1...File 'metric.for'. C SRFC1 and subsequent routines... File 'srfc.for' and subsequent C files (like 'val.for' and 'fit.for'). C PARM1 and subsequent routines... File 'parm.for' and subsequent C files (like 'val.for' and 'fit.for'). C C Date: 2014, February 14 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER I,J,L C READ(LUN,*) TEXTM I=0 NEXPV=1 NEXPQ=1 IVERT=0 READ(LUN,*) I,NEXPV,NEXPQ,IVERT CALL METR1(I,LUN) READ(LUN,*) BOUNDM DO 10, I=1,3 IF (BOUNDM(2*I).NE.BOUNDM(2*I-1)) THEN IF (ABS(BOUNDM(2*I)).GT.100.*ABS(BOUNDM(2*I)-BOUNDM(2*I-1))) * THEN C 316 CALL WARN('316 in MODEL: Large value of coordinates') C The absolute value of coordinates is too large compared with C the dimension of the model (the interval between the minimum C and the maximum coordinate). This may cause too large C rounding errors. The coordinates should be comparable C in order of magnitude with the dimension of the model. ENDIF ENDIF 10 CONTINUE READ(LUN,*) NSRFCS C C Simple blocks: C Number of simple blocks READ(LUN,*) NSB C Initializing memory for indices of surfaces bounding simple blocks L=NSB+1 DO 11 I=L,MSB KSB(I)=0 11 CONTINUE C Reading indices of surfaces bounding simple blocks: DO 14 J=1,NSB READ(LUN,*) (KSB(I),I=L,MSB) DO 12 I=L,MSB IF(IABS(KSB(I)).GT.NSRFCS) THEN C 311 CALL ERROR('311 in MODEL: Block bounded by wrong interface') C Index of the surface bounding the simple block is greater C than the specified number of surfaces. ELSE IF(KSB(I).EQ.0) THEN KSB(J)=I-1 L=I GO TO 13 END IF 12 CONTINUE GO TO 99 13 CONTINUE 14 CONTINUE C C Complex blocks: C Number of complex blocks READ(LUN,*) NCB C Initializing memory for indices of simple blocks forming c. blocks L=NCB+1 DO 21 I=L,MCB KCB(I)=0 21 CONTINUE C Reading indices of simple blocks forming complex blocks DO 24 J=1,NCB READ(LUN,*) (KCB(I),I=L,MCB) DO 22 I=L,MCB IF(KCB(I).LT.0.OR.KCB(I).GT.NSB) THEN C 312 CALL ERROR * ('312 in MODEL: C. block composed of wrong s. block.') C Complex block composed of wrong simple block: C Index of a simple block composing the complex block is C greater than the specified number of simple blocks. ELSE IF(KCB(I).EQ.0) THEN KCB(J)=I-1 L=I GO TO 23 END IF 22 CONTINUE GO TO 99 23 CONTINUE 24 CONTINUE C C Smooth surfaces: CALL SRFC1(LUN,NSRFCS) C C Material parameters: CALL PARM1(LUN,NCB) CALL RSEP3I('NEGPAR',NEGPAR,0) RETURN C 99 CONTINUE C 310 CALL ERROR('310 in MODEL1: Insufficient memory in /MODELC/') C Insufficient memory for the input data in common block /MODELC/. C The dimensions MSB or MCB of arrays KSB or KCB, respectively, C must be enlarged. Refer to include file C model.inc END C C======================================================================= C C C INTEGER FUNCTION NSRFC() C C Integer function NSRFC is designed to return the number of surfaces C covering structural interfaces. C C No input. C C Output: C NSRFC...Number of surfaces covering structural interfaces. C C Common block /MODELC/: INCLUDE 'model.inc' C model.inc C None of the storage locations of the common block are altered. C C No subroutines and external functions required. C C Date: 1989, December 19 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C No auxiliary storage locations. C NSRFC=NSRFCS RETURN END C C======================================================================= C C C SUBROUTINE BLOCK(COOR,ISRF1,ISB1,ISRF2,ISB2,ICB2) REAL COOR(3) INTEGER ISRF1,ISB1,ISB2,ICB2,ISRF2 C C This subroutine searches for the simple block and the complex block in C which a given point is situated. This routine may be also used to C determine the index of a block touching a specified block at a given C point situated on the boundary of the specified block (the situation C which may occur when a ray impinges on a boundary of a block). C Another function of the routine is to determine the index of the C surface bounding a given simple block and separating it from the given C point. C C Input: C COOR... Array containing coordinates X1, X2, X3 of a given point. C ISRF1...Controls the possible determination of the simple block C touching block ISB1 at the given point. C ISRF1.EQ.0: The simple block touching ISB1 at the given C point need not be determined. This is the most frequent C option, used, e.g., if the given point is not situated C at a structural interface or if the simple block C touching ISB1 is already known. C ISRF1.NE.0: Index of the surface at which the given point C is situated and which separates simple block ISB1 from C another block to be determined and returned in ISB2. C Since the given point may be situated on either side of C the surface, surface IABS(ISRF1) is skipped in the list C of surfaces limiting simple blocks when determining the C position of the given point with respect to simple C blocks. The sign of ISRF1 is ignored. C This option cannot be used together with ISB1.EQ.0. C ISB1... Index of the simple block which position with respect to C the given point will be checked. C ISB1.EQ.0: No simple block given. The simple block in C which the given point is situated will be determined C and returned in ISB2. ISRF2 will be set to zero. C The determination starts with initial estimate ISBOLD, C continuing with ISBOLD+1, ISBOLD-1, ISBOLD+2, ISBOLD-2, C etc. Here ISBOLD is the simple block returned in ISB2 C during the last invocation of this subroutine with C ISB1.EQ.0 on input and ISB2.NE.0 on output. ISBOLD=1 C before such an invocation. C ISB1.GT.0: The position of the given point with C respect to simple block ISB1 will be examined. C If the point is situated in ISB1, the output value C of ISRF2 is set to 0. C If the point is situated in ISB1 and ISRF1.EQ.0, C the output value of ISB2 is set to ISB1. C If the point is situated in ISB1 and ISRF1.NE.0, C the simple block containing the given point, C separated from block ISB1 by surface ISRF1 and not C separated from block ISB1 by another surface will be C determined and returned in ISB2. C If the point is not situated in ISB1, the surface C separating the point from ISB1 is output in ISRF2. C In addition, the simple block containing the given C point, separated from block ISB1 by surface ISRF2 and C not separated from block ISB1 by another surface will C be determined and returned in ISB2. C The determination starts with initial estimate ISB1, C continuing with ISB1+1, ISB1-1, ISB1+2, ISB1-2, etc. C None of the input parameters are altered. C C Output: C ISRF2...For the given point not situated inside block ISB1 or C at its boundary, ISRF2 has the meaning of the index of C one of the surfaces bounding simple block ISB1 and C separating the given point from simple block ISB1, C supplemented by a sign '+' or '-' for the simple block C ISB1 situated at the positive or negative side of the C surface, respectively. If the point is separated from C simple block ISB1 by several surfaces, the first one of C the surfaces, for which ISB2 (see description of ISB2 C below) will be positive, is preferred. Roughly speaking, C surface ISRF2 at which simple block ISB1 touches simple C block containing the given point is preferred. C ISRF2=0 if ISB1=0 or if the given point is situated inside C block ISB1 or at its boundary. C If ISRF1.NE.0, the position with respect to surface C IABS(ISRF1) is not checked. C ISB2... Index of the simple block containing the given point. C If ISB1.EQ.0: C ISB2 is the index of the simple block containing the C given point. C The first simple block of ISBOLD, ISBOLD+1, ISBOLD-1, C ISBOLD+2, ISBOLD-2, etc., containing the given point, C is returned in ISB2. Here ISBOLD is the simple block C returned in ISB2 during the previous invocation of this C subroutine with ISB1.EQ.0. ISBOLD=1 during such the C first invocation. C If ISB1.NE.0: C If ISRF1.EQ.0 and ISRF2.EQ.0, ISB2 must not be C separated from simple block ISB1 by any surface. C A surface "separates" two simple blocks if both blocks C are bounded by the surface and are situated at its C opposite sides. C If ISRF1.NE.0 and ISRF2.EQ.0, ISB2 must be separated C from simple block ISB1 by surface IABS(ISRF1) and C must not be separated by another surface. C If ISRF2.NE.0, ISB2 must be separated from simple block C ISB1 by surface ISRF2 and must not be separated by C another surface. C The first simple block of ISB1, ISB1+1, ISB1-1, C ISB1+2, ISB1-2, etc., having the above properties, is C returned in ISB2. C If there is no simple block of the above properties, C ISB2=0. C ICB2... Index of the complex block in which simple block ISB2 is C situated. ICB2=0 if ISB2=0 or if simple block ISB2 is not C situated in any material complex block. C C Examples: C For each point of a set of discrete points, we wish to determine C simple block ISB and complex block ICB in which the point C is situated. C (a) We thus use C CALL BLOCK(COOR,0,0,I,ISB,ICB) C When tracing a curve (e.g., a ray), we wish rather to find each C crossed interface than only to determine the simple blocks C at discrete points along the curve. C (b) If the previous point along the curve is situated in C simple block ISB, we use C CALL BLOCK(COOR,0,ISB,ISRF1,ISBNEW,ICBNEW) C and get ISRF1.EQ.0 if no surface limiting simple block ISB C has been crossed. C (c) If ISRF1.NE.0, IABS(ISRF1) is the surface which has been C crossed. We should determine the point of intersection of C the curve with surface IABS(ISRF1). Then we use C CALL BLOCK(COOR,ISRF1,ISB,ISRF2,ISB2,ICB2) C to find simple block ISB2 and complex block ICB2 at the C other side of surface ISRF1. If ISRF2.EQ.0, the blocks C are successfully determined. If ISRF2.NE.0, surface C IABS(ISRF2) has also been crossed and we have to repeat C step (c) with ISRF1=ISRF2. C When tracing a curve along surface ISRFA, we call subroutine C BLOCKS instead of BLOCK, with analogous arguments: C (b') ISRF(1)=0 C ISRF(2)=ISRFA C CALL BLOCKS(COOR,2,ISRF,ISB,ISRF1,ISBNEW,ICBNEW) C (c') ISRF(1)=ISRF1 C ISRF(2)=ISRFA C CALL BLOCKS(COOR,2,ISRF,ISB,ISRF2,ISB2,ICB2) C Analogously, when tracing the curve of intersection of surfaces C ISRFA and ISRFB, we call subroutine BLOCKS in the C following way: C (b") ISRF(1)=0 C ISRF(2)=ISRFA C ISRF(3)=ISRFB C CALL BLOCKS(COOR,3,ISRF,ISB,ISRF1,ISBNEW,ICBNEW) C (c") ISRF(1)=ISRF1 C ISRF(2)=ISRFA C ISRF(3)=ISRFB C CALL BLOCKS(COOR,3,ISRF,ISB,ISRF2,ISB2,ICB2) C C Subroutines and external functions required: EXTERNAL BLOCKS C BLOCKS... This file. C C Date: 1999, March 20 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C INTEGER ISRF(1) C ISRF(1)=ISRF1 CALL BLOCKS(COOR,1,ISRF,ISB1,ISRF2,ISB2,ICB2) RETURN END C C======================================================================= C C C SUBROUTINE BLOCKS(COOR,NSRF1,ISRF,ISB1,ISRF2,ISB2,ICB2) REAL COOR(3) INTEGER NSRF1,ISRF(NSRF1),ISB1,ISRF2,ISB2,ICB2 C C This subroutine is a generalization of subroutine BLOCK to follow C points along structural interfaces and edges. Subroutine BLOCKS C searches for the simple block and the complex block in which a given C point is situated. The given point may be situated in the vicinity C of one or more structural interfaces. In such case, blocks at the C proper side of the interface(s) will be searched for. This routine C may also be used to determine the index of a block touching a C specified block at a given point situated on the boundary of the C specified block (the situation which may occur when a ray impinges on C a boundary of a block). Another function of the routine is to C determine the index of the surface bounding a given simple block and C separating it from the given point. C C Input: C COOR... Array containing coordinates X1, X2, X3 of a given point. C NSRF1...Number of surface indices specified in array ISRF. C There must be NSRF1.GE.1. For ISB1.EQ.0., there must be C NSRF1.EQ.1 and ISRF(1).EQ.0. C ISRF... Array containing the indices of the surfaces at which the C given point is situated. Since the given point may be C situated on either side of the surfaces, the surfaces C listed in array ISRF are skipped in the list of surfaces C limiting simple blocks when determining the position of C the given point with respect to simple blocks. C The signs of the indices in array ISRF are ignored. C ISRF(1): Controls the possible determination of the C simple block touching block ISB1 at the given point. C ISRF(1).EQ.0: The simple block touching ISB1 at the C given point need not be determined. This is the most C frequent option, used, e.g., if the given point is not C situated at a structural interface or if the simple C block touching ISB1 is already known. C ISRF(1).NE.0: Index of the surface at which the given C point is situated and which separates simple block C ISB1 from another block to be determined. Refer to C the description of output parameter ISB2. C This option cannot be used together with ISB1.EQ.0. C ISRF(2) to ISRF(NSRF1): Indices of the surfaces at C which the given point is situated, other than ISRF(1). C These indices are designed to specify the surface along C which a curve is traced, or the surfaces which curve of C intersection is traced. C ISB1... Index of the simple block which position with respect to C the given point will be checked. C ISB1.EQ.0: No simple block given. The simple block in C which the given point is situated will be determined C and returned in ISB2. ISRF2 will be set to zero. C The determination starts with initial estimate ISBOLD, C continuing with ISBOLD+1, ISBOLD-1, ISBOLD+2, ISBOLD-2, C etc. Here ISBOLD is the simple block returned in ISB2 C during the last invocation of this subroutine with C ISB1.EQ.0 and ISB2.NE.0. ISBOLD=1 before such an C invocation. C ISB1.GT.0: The position of the given point with C respect to simple block ISB1 will be examined. C If the point is situated in ISB1, the output value C of ISRF2 is set to 0. C If the point is situated in ISB1 and ISRF(1).EQ.0, C the output value of ISB2 is set to ISB1. C If the point is situated in ISB1 and ISRF(1).NE.0, C the simple block containing the given point, C separated from block ISB1 by surface ISRF(1) and not C separated from block ISB1 by another surface will be C determined and returned in ISB2. C If the point is not situated in ISB1, the surface C separating the point from ISB1 is output in ISRF2. C In addition, the simple block containing the given C point, separated from block ISB1 by surface ISRF2 and C not separated from block ISB1 by another surface will C be determined and returned in ISB2. C The determination starts with initial estimate ISB1, C continuing with ISB1+1, ISB1-1, ISB1+2, ISB1-2, etc. C None of the input parameters are altered. C C Output: C ISRF2...For the given point not situated inside block ISB1 or C at its boundary, ISRF2 has the meaning of the index of C one of the surfaces bounding simple block ISB1 and C separating the given point from simple block ISB1, C supplemented by a sign '+' or '-' for the simple block C ISB1 situated at the positive or negative side of the C surface, respectively. If the point is separated from C simple block ISB1 by several surfaces, the first one with C positive ISB2 (see below) is preferred. C ISRF2=0 if ISB1=0 or if the given point is situated inside C block ISB1 or at its boundary. C Note that surfaces ISRF(1) to ISRF(NSRF1) are skipped C when determining the position of the given point with C respect to simple blocks. C ISB2... Index of the simple block containing the given point. C If ISB1.EQ.0: C ISB2 is the index of the simple block containing the C given point. C The first simple block of ISBOLD, ISBOLD+1, ISBOLD-1, C ISBOLD+2, ISBOLD-2, etc., containing the given point, C is returned in ISB2. Here ISBOLD is the simple block C returned in ISB2 during the previous invocation of this C subroutine with ISB1.EQ.0. ISBOLD=1 during such the C first invocation. C If ISB1.NE.0: C If ISRF(1).EQ.0 and ISRF2.EQ.0, ISB2 must not be C separated from simple block ISB1 by any surface. C A surface "separates" two simple blocks if both blocks C are bounded by the surface and are situated at its C opposite sides. C If ISRF(1).NE.0 and ISRF2.EQ.0, ISB2 must be separated C from simple block ISB1 by surface IABS(ISRF(1)) and C must not be separated by another surface. C If ISRF2.NE.0, ISB2 must be separated from simple block C ISB1 by surface ISRF2 and must not be separated by C another surface. C The first simple block of ISB1, ISB1+1, ISB1-1, C ISB1+2, ISB1-2, etc., having the above properties, is C returned in ISB2. C If there is no simple block of the above properties, C ISB2=0. C ICB2... Index of the complex block in which simple block ISB2 is C situated. ICB2=0 if ISB2=0 or if simple block ISB2 is not C situated in any material complex block. C C Common block /MODELC/: INCLUDE 'model.inc' C model.inc C None of the storage locations of the common block are altered. C C Subroutines and external functions required: EXTERNAL ISIDE,INTERF,SRFC2 INTEGER ISIDE C ISIDE,INTERF... This file. C SRFC2 and subsequent routines... File 'srfc.for' and subsequent C files (like 'val.for' and 'fit.for'). C C Date: 1999, March 20 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER MSEPAR PARAMETER (MSEPAR=20) INTEGER NSEPAR,KSEPAR(MSEPAR),JSEPAR,ISEPAR,ISIDE1,I,J,K,II(1) INTEGER ISBOLD,ISB0,ISRF1 SAVE ISBOLD DATA ISBOLD/1/ C C....................................................................... C ISRF1=ISRF(1) C C Checking input values: IF(ISRF1.LT.-NSRFCS.OR.ISRF1.GT.NSRFCS) THEN C 313 CALL ERROR('313 in BLOCK: Wrong index of surface') C Absolute value of the input parameter ISRFC1 (index of the C surface) is greater than the number NSRFC of the surfaces C covering structural interfaces. END IF IF(ISB1.LT.0.OR.ISB1.GT.NSB) THEN C 314 CALL ERROR('314 in BLOCK: Wrong index of simple block') C Parameter ISB1 (index of the simple block) is either C negative or greater than the number NSB of simple blocks. END IF IF(ISB1.EQ.0) THEN IF(NSRF1.NE.1.OR.ISRF1.NE.0) THEN C 315 CALL ERROR('315 in BLOCK: No simple block specified') C If no simple block ISB1 is specified, NSRF1 must be 1 and C ISRF1=ISRF(1) must be 0. END IF END IF C C Initial simple block ISB0: IF(ISB1.EQ.0) THEN ISB0=ISBOLD ELSE ISB0=ISB1 END IF C C Position of the given point with respect to simple block ISB0: ISRF2=0 CALL INTERF(COOR,NSRF1,ISRF,ISB0,MSEPAR,NSEPAR,KSEPAR) IF(NSEPAR.EQ.0) THEN C The point is inside simple block ISB0: IF(ISRF1.EQ.0) THEN ISB2=ISB0 GO TO 90 ELSE NSEPAR=1 KSEPAR(1)=ISRF1 END IF ELSE IF(ISB1.EQ.0) THEN NSEPAR=1 ELSE ISRF2=KSEPAR(1) END IF END IF C C Search for the simple block in which the given point is situated: DO 20 JSEPAR=1,NSEPAR IF(ISB1.NE.0) THEN ISEPAR=IABS(KSEPAR(JSEPAR)) ISIDE1=-ISIDE(ISEPAR,ISB1) END IF C Loop over simple blocks DO 19 J=1,MAX0(ISB0-1,NSB-ISB0) DO 18 ISB2=ISB0+J,ISB0-J,-2*J IF(ISB2.GT.0.AND.ISB2.LE.NSB) THEN C Selecting simple block ISB2 according to ISB1 IF(ISB1.NE.0) THEN C Loop for surfaces bounding block ISB1 DO 11 I=KSB(ISB1-1)+1,KSB(ISB1) C Skipping simple blocks separated from ISB1 by another C surface than ISEPAR K=KSB(I) IF(IABS(K).NE.ISEPAR) THEN IF(ISIDE(K,ISB1).EQ.-ISIDE(K,ISB2)) THEN GO TO 17 END IF END IF 11 CONTINUE C Skipping simple blocks not separated from ISB1 by ISEPAR IF(ISIDE(ISEPAR,ISB2).NE.ISIDE1) THEN GO TO 17 END IF END IF C Determining the position of the given point with respect C to the simple block CALL INTERF(COOR,NSRF1,ISRF,ISB2,1,I,II) IF(I.EQ.0) THEN IF(ISB1.EQ.0) THEN ISBOLD=ISB2 ELSE ISRF2=KSEPAR(JSEPAR) END IF GO TO 90 END IF 17 CONTINUE END IF 18 CONTINUE 19 CONTINUE 20 CONTINUE C No simple block has been found: ISB2=0 ICB2=0 RETURN C C Determination of the complex block: 90 CONTINUE DO 92 J=1,NCB DO 91 I=KCB(J-1)+1,KCB(J) IF(KCB(I).EQ.ISB2) THEN ICB2=J RETURN END IF 91 CONTINUE 92 CONTINUE C No complex block: ICB2=0 RETURN END C C======================================================================= C C C INTEGER FUNCTION ISIDE(ISRF,ISB) INTEGER ISRF,ISB C C This is an auxiliary function to the subroutine BLOCKS. C This function determines the mutual position of a surface and a simple C block. C C Input: C ISRF... Index of the surface. The sign is ignored. C ISB... Index of the simple block. C None of the input parameters are altered. C C Output: C ISIDE...ISIDE=-1: The simple block is bounded by the surface and C is situated on its negative side. C ISIDE= 1: The simple block is bounded by the surface and C is situated on its positive side. C ISIDE= 2: The simple block is not bounded by the surface. C C Common block /MODELC/: INCLUDE 'model.inc' C model.inc C None of the storage locations of the common block are altered. C C No subroutines and external functions required. C C Date: 1989, December 19 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER IS,LS,MS C LS=KSB(ISB-1)+1 MS=KSB(ISB) C C Loop for surfaces bounding simple block ISB: DO 1 IS=LS,MS IF(IABS(KSB(IS)).EQ.IABS(ISRF)) THEN ISIDE=ISIGN(1,KSB(IS)) RETURN END IF 1 CONTINUE C ISIDE=2 RETURN END C C======================================================================= C C C SUBROUTINE INTERF(COOR,NSRF1,ISRF1,ISB,MSRF2,NSRF2,ISRF2) REAL COOR(3) INTEGER NSRF1,ISRF1(NSRF1),ISB,MSRF2,NSRF2,ISRF2(MSRF2) C C This is an auxiliary subroutine to the subroutine BLOCKS. C This subroutine determines the position of a given point with respect C to a given simple block. C C Input: C COOR... Array containing coordinates X1, X2, X3 of a given point. C NSRF1...Number of surface indices specified in array ISRF1. C Must be NSRF1.GE.1. C ISRF1...Array containing the indices of the surfaces at which the C given point is situated. Since the given point may be C situated on either side of the surfaces, the surfaces C listed in array ISRF1 are skipped in the list of surfaces C limiting the simple block when determining the position of C the given point with respect to the simple block. C The signs of the indices in array ISRF1 are ignored. C Zero indices are allowed and have no neaning. C ISB... Index of the given simple block. C MSRF2...Dimension of array ISRF2. C None of the input parameters are altered. C C Output: C NSRF2...Number of surfaces separating the given point from simple C block ISB. C NSRF2=0 if the given point is situated inside simple C block ISB. C ISRF2...Indices of surfaces separating the given point from simple C block ISB, supplemented by sign '+' or '-' for simple C block ISB situated at the positive or negative side of the C surface, respectively. C The first MSRF2 such surfaces from the list of surfaces C limiting the simple block are reported. C C Common block /MODELC/: INCLUDE 'model.inc' C model.inc C None of the storage locations of the common block are altered. C C Subroutines and external functions required: EXTERNAL SRFC2 C SRFC2 and subsequent routines... File 'srfc.for' and subsequent C files (like 'val.for' and 'fit.for'). C C Date: 1999, March 20 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER IS,JS,KS,I REAL F(10) C F... Auxiliary array to contain the value and partial C derivatives F, F1, F3, F11, F12, F22, F13, F23, F33 of the C function describing surfaces at the given point. C C Speed-up storage locations: INTEGER MOLD PARAMETER (MOLD=20) INTEGER NOLD,IOLD(MOLD) REAL COLD(3) SAVE COLD,NOLD,IOLD DATA COLD/3*-999999./ C C....................................................................... C C Already examined surfaces: IF(COOR(1).NE.COLD(1).OR. * COOR(2).NE.COLD(2).OR.COOR(3).NE.COLD(3)) THEN NOLD=0 END IF C C Loop for surfaces bounding simple block ISB: NSRF2=0 DO 9 IS=KSB(ISB-1)+1,KSB(ISB) KS=KSB(IS) JS=IABS(KS) C Skipping surfaces of array ISRF1 DO 1 I=1,NSRF1 IF(JS.EQ.IABS(ISRF1(I))) THEN GO TO 8 END IF 1 CONTINUE C Surfaces already examined DO 2 I=1,NOLD IF(JS.EQ.IABS(IOLD(I))) THEN IF(IOLD(I)*KS.LT.0) THEN NSRF2=NSRF2+1 ISRF2(NSRF2)=KS IF(NSRF2.GE.MSRF2) THEN RETURN END IF END IF GO TO 8 END IF 2 CONTINUE C The surface has to be examined CALL SRFC2(JS,COOR,F) IF(NOLD.LT.MOLD) THEN IF(F(1).LT.0.) THEN NOLD=NOLD+1 IOLD(NOLD)=-JS ELSE IF(F(1).GT.0.) THEN NOLD=NOLD+1 IOLD(NOLD)= JS END IF END IF IF(F(1)*FLOAT(KS).LT.0.) THEN NSRF2=NSRF2+1 ISRF2(NSRF2)=KS IF(NSRF2.GE.MSRF2) THEN RETURN END IF END IF 8 CONTINUE 9 CONTINUE RETURN END C C======================================================================= C C C SUBROUTINE SEPAR(ISB1,ISB2,NSRF,ISRF) INTEGER ISB1,ISB2,NSRF,ISRF C C Subroutine determining the surface separating given simple blocks, C i.e., the surface limiting the simple blocks, which are situated one C at the positive side of the surface and the other at the negative side C of the surface. C C Input: C ISB1,ISB2... Indices of given simple blocks. C None of the input parameters are altered. C C Output: C NSRF... Number of surfaces separating the simple blocks. C ISRF... Index of a surface separating the simple blocks, C supplemented by sign minus if simple block ISB1 is C situated at its negative side. C If NSRF.EQ.0: ISRF contains unchanged input value. C If NSRF.GT.1: ISRF is one of the surfaces separating the C simple blocks. C C Common block /MODELC/: INCLUDE 'model.inc' C model.inc C None of the storage locations of the common block are altered. C C No subroutines and external functions required. C C Date: 1999, May 18 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER I1,I2,KSBI2 C NSRF=0 C C Loop for surfaces bounding simple block ISB2 DO 2 I2=KSB(ISB2-1)+1,KSB(ISB2) KSBI2=-KSB(I2) C Loop for surfaces bounding simple block ISB1 DO 1 I1=KSB(ISB1-1)+1,KSB(ISB1) IF(KSB(I1).EQ.KSBI2) THEN NSRF=NSRF+1 IF(NSRF.EQ.1) THEN ISRF=KSB(I1) END IF END IF 1 CONTINUE 2 CONTINUE RETURN END C C======================================================================= C C C SUBROUTINE VELOC(IWAVE,UP,US,QP,QS,VP,VS,VD,QL) INTEGER IWAVE REAL UP(10),US(10),QP,QS,VP,VS,VD(10),QL C C This subroutine transforms the values of parameters of the medium into C velocities and loss factors. C C Input: C IWAVE...Type of wave. C IWAVE.GE.0: P wave, C IWAVE.LT.0: S wave. C UP,US...Powers of P and S wave velocities and their first and C second partial derivatives (the exponent of the powers is C NEXPV, see 'Input data for the model'), in order U, U1, C U2, U3, U11, U12, U22, U13, U23, U33. C QP,QS...Powers of the loss factors of P and S waves (the exponent C of the powers is NEXPQ, see 'Input data for the model'). C None of the input parameters are altered. C C Output: C VP,VS...P and S wave velocities. C VD... Velocity and its first and second partial derivatives C corresponding to the wave specified by IWAVE, C in order V, V1, V2, V3, V11, V12, V22, V13, V23, V33. C QL... Loss factor corresponding to the wave specified by IWAVE. C C Common block /MODELC/: INCLUDE 'model.inc' C model.inc C None of the storage locations of the common block are altered. C C Subroutines and external functions required: EXTERNAL FPOWER C FPOWER...This file. C C Date: 2005, June 1 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage location: REAL POWER,AUX1(1),AUX2(1) C POWER=FLOAT(NEXPV) IF(IWAVE.GE.0) THEN CALL FPOWER(10,UP,POWER,VD) CALL VAR5(1,1) VP=VD(1) CALL FPOWER(1,US,POWER,AUX2) VS=AUX2(1) AUX1(1)=QP ELSE CALL FPOWER(1,UP,POWER,AUX2) VP=AUX2(1) CALL FPOWER(10,US,POWER,VD) CALL VAR5(2,2) VS=VD(1) AUX1(1)=QS END IF CALL FPOWER(1,AUX1,FLOAT(NEXPQ),AUX2) QL=AUX2(1) RETURN END C C======================================================================= C C C SUBROUTINE FPOWER(N,FINP,POWER,FOUT) INTEGER N REAL FINP(N),POWER,FOUT(N) C C This subroutine evaluates the value and, possibly, the three first and C six second partial derivatives of a function if the value and the C three first and six second partial derivatives of the POWER-th power C of the function are known. C C Input: C N... For N=1: only the function value is evaluated. The C derivatives are ignored. C For N=4: the value and the three first partial derivatives C are evaluated. C For N=10: the value and the three first and six second C partial derivatives are evaluated. C FINP... Array containing the value, the first and second partial C derivatives of the POWER-th power of the function to be C evaluated, in the order F, F1, F2, F3, F11, F12, F22, F13, C F23, F33. For N=1, only the function value is required. C POWER...The specified function is equal to the POWER-th power of C the corresponding physical quantity. C POWER=0: Zero output array FOUT is generated. C None of the input parameters are altered (except FINP if this C parameter and FOUT are identical in the calling sequence). C C Output: C FOUT... Array containing the value, the first and second partial C derivatives of the evaluated function, in the order F, F1, C F2, F3, F11, F12, F22, F13, F23, F33. This parameter may C coincide with FINP, in which case FINP is destroyed on C output. Note that this coincidence is an exception to C ANSI standard of FORTRAN 77. C C No subroutines and external functions required. C C Date: 1999, January 15 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER I REAL F,AUX1,AUX2 C IF(POWER.EQ.0.) THEN DO 1 I=1,N FOUT(I)=0. 1 CONTINUE ELSE IF(0.999.LT.POWER.AND.POWER.LT.1.001) THEN DO 2 I=1,N FOUT(I)=FINP(I) 2 CONTINUE CALL VAR4(0,1.) ELSE IF(FINP(1).LT.0.) THEN C 317 CALL ERROR('317 in FPOWER: Negative material parameter') C Nonunit power of a material parameter is not allowed to be C negative. The negative value may be caused by oscillatory C character of interpolated positive values. ELSE IF(FINP(1).EQ.0.) THEN IF(POWER.LT.0.) THEN C 318 CALL ERROR('318 in FPOWER: Zero inverse material parameter') C Negative power of a material parameter cannot be zero. ELSE FOUT(1)=0. DO 3 I=2,N IF(FINP(I).NE.0.) THEN C 319 CALL ERROR('319 in FPOWER: Zero material parameter') C Nonunit power of zero material parameter is not allowed to C have nonzero derivatives. END IF FOUT(I)=0. 3 CONTINUE END IF ELSE IF(-1.001.LT.POWER.AND.POWER.LT.-0.999) THEN F=1./FINP(1) ELSE F=FINP(1)**(1./POWER) END IF FOUT(1)=F IF(N.GT.1) THEN AUX1= F/(FINP(1)*POWER) AUX2= (POWER-1.)/F FOUT(2)=AUX1*FINP(2) FOUT(3)=AUX1*FINP(3) FOUT(4)=AUX1*FINP(4) IF(N.GT.4) THEN FOUT(5)=AUX1*FINP(5)-AUX2*FOUT(2)*FOUT(2) FOUT(6)=AUX1*FINP(6)-AUX2*FOUT(2)*FOUT(3) FOUT(7)=AUX1*FINP(7)-AUX2*FOUT(3)*FOUT(3) FOUT(8)=AUX1*FINP(8)-AUX2*FOUT(2)*FOUT(4) FOUT(9)=AUX1*FINP(9)-AUX2*FOUT(3)*FOUT(4) FOUT(10)=AUX1*FINP(10)-AUX2*FOUT(4)*FOUT(4) END IF CALL VAR4(0,AUX1) CALL VAR4(2,-AUX2*FOUT(2)) CALL VAR4(3,-AUX2*FOUT(3)) CALL VAR4(4,-AUX2*FOUT(4)) END IF END IF RETURN END C C======================================================================= C