C
C Subroutine file 'code.for' - codes for elementary waves. C C Version: 6.50 C Date: 2011, May 19 C C Coded by Ludek Klimes C Department of Geophysics, Charles University in Prague, C Ke Karlovu 3, 121 16 Praha 2, Czech Republic, C http://sw3d.cz/staff/klimes.htm C C....................................................................... C C This file consists of: C Codes for elementary waves - general description C Codes for elementary waves C CODE1...Subroutine designed to read the input data for the codes C of elementary waves and to store them in the common block C /COD/. C CODE1 C LCODE...Integer function returning the length of the code of the C current elementary wave. C LCODE C NELEM...Integer function returning the number of ray elements C corresponding to the given position in the code of the C elementary wave. C NELEM C CODE... Subroutine designed to transform the used numerical code C of the elementary wave under consideration into C instructions specifying the behaviour of the wave at the C initial points of rays and at all points of incidence of C the rays at interfaces (boundaries of complex blocks). C CODE C C....................................................................... C C C Input data CODE for the codes of elementary waves: 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 TEXTC), the input parameter is of the C type REAL. C (1) TEXTC C String describing the data. Only the first 80 characters of the C string are significant. C (2) NSKIP C The number of elementary waves to be skipped over. The C (NSKIP+1)-th elementary wave in these input data will be the first C computed wave. Default value: 0. C For each elementary wave IWAVE the following data (3): C (3) KODTYP,KODE C KODTYP..Determines the type of the code C 1... Block-surface code 1 C 2... Basic block code 2 C 3... Compound-element block code 3 C 4... Generalized layer code 4 C 5... Reflection-transmission code 5 C KODE... The code of the elementary wave in the form of a finite C sequence of integers terminated by a slash. C Attention: C Generally, different types of the code may be combined C within the same input data, but the user has then be very C careful. If combining the different types of the code, C output parameters IWAVE0 and IKODE of subroutine CODE1 may C have no sense. Then: C If the parameter MODCRT of the data set 'DCRT' (see source C code file 'ray.for') is zero, the output parameters C IWAVE0 and IKODE of subroutine CODE1 are hopefully C unused. C If the parameter MODCRT of the data set 'DCRT' is not C zero, mixing the types of the wave code should be C avoided unless the user finds a way how to do it C reasonably. C (4) a slash. C Example of data set CODE C C....................................................................... C C Storage in the memory: C The input data (2) for the computed elementary wave are stored in C the common block /COD/ defined in the include file 'code.inc'. C code.inc C C======================================================================= C C C Codes for elementary waves - general description C C We consider ordinary seismic body waves, such as refracted, C primarily or multiply reflected, possibly converted waves. In C general, incidence of a wave at an interface (boundary of a C complex block) produces four waves, reflected P and S, and C transmitted P and S waves. When performing complete ray tracing, C we must know a priori which of the four generated waves to follow. C This decision must be made at each interface. The alphanumeric C string specifying the behaviour of a ray from its initial point to C its endpoint is a 'code'. In the following, we shall consider the C code to be a sequence of nonzero integers. C C The term 'elementary wave' does not have unique meaning in the C literature. Here we apply the term elementary wave to that part C of the wavefield that is described by one specific code. Since C there may be various types of codes, there is also a variety of C divisions of the wave field into elementary waves. C C We introduce the term 'element of a ray', which has an important C meaning in the construction of codes. By an element of a ray, we C denote that part of the ray that is situated in one complex block C between two successive points of reflection/transmission, or C between the initial point or endpoint of the ray and the closest C point of reflection/transmission, or between the initial point and C the endpoint of the ray, if the ray is entirely situated in one C complex block. C C In the following, we present several possible types of numerical C codes of elementary waves. C C Examples of code types C C (1) Block-surface code C Two integers are used for any element of the ray. The C absolute value of the first integer specifies the index of C the complex block, in which the element of the ray is C situated, the sign specifies the type of the wave (plus C sign... P wave, minus sign... S wave). The second C number is the index of the smooth surface, at which the C element terminates. The code of the ray is a chain of the C above doublets describing the ray from its initial point C to its endpoint. The index of a surface, at which the C last element of a ray terminates, need not be specified. C In this case the ray is allowed to terminate on any C surface bounding the complex block in which the last C element of the ray is situated. Application of this code C is straightforward. C C If any complex block is bounded by several smooth C surfaces, the block-surface code may divide artificially C the wave field into several elementary waves. This would C result in repeated application of the complete ray tracing C to all of the elementary waves, what may be time C consuming. In such situations, it is therefore desirable C to use other codes, for which the elementary wave has a C more general meaning. C C C (2) Basic block code C This code may be obtained from the block-surface code if C the indices of the surfaces, at which individual elements C of a ray terminate, are omitted. Then, any element of the C ray is described by a single integer specifying the C complex block and wave type, and it may terminate at any C surface bounding the complex block. The code of the ray C is again a chain of the numbers describing successively C its individual elements. C C Elementary waves specified by the block-surface code, C passing through the same complex blocks, and reflected C from different surfaces bounding a complex block, are C united into one elementary wave in the basic block code. C Similarly, elementary waves specified by the block-surface C code and transmitted from a block into the neighbouring C block across different surfaces separating these blocks, C are united into one elementary wave by the basic block C code. C C C (3) Compound-element block code C Let us introduce the terms simple and compound element of C a ray. For this purpose, we call the 'lower-index C boundary' of a complex block that part of its boundary, C which separates the complex block either from complex C blocks with lower indices or from a free space. The C remaining part of the boundary is called the C 'higher-index boundary' of the complex block. The initial C point of a ray is treated in the same way as the points C situated at a lower-index boundary. C C An element of a ray is called 'simple element' if its C initial point (e.g. the point where the wave enters the C corresponding complex block) and its endpoint (e.g. the C point at which the wave leaves the complex block) are C situated one on the lower-index and the other on the C higher-index boundary of the complex block. The C 'compound element' is such an element for which its C initial point and its endpoint are both situated either C on the lower-index or on the higher-index boundary of the C complex block. The compound element is formally C considered as two simple elements. C C Often it is convenient to work with waves, the rays of C which have either a compound element in a complex block or C two simple elements of the same wave type (unconverted C reflection) in the same block, as with one elementary C wave. The compound-element block code makes a division of C the wave field into such elementary waves possible. C C We introduce the compound-element block code as follows. C For any simple element of a ray one number is used. The C absolute value of this number specifies the index of the C complex block, in which the simple element is situated. C The sign of this number specifies the type of the wave C (plus sign... P wave, minus sign... S wave). For any C compound element of a ray, two identical numbers are used. C A compound element and a doublet of simple elements, C described by the same numbers, are not distinguished. The C code of a ray is a chain of these numbers describing the C ray from its initial point to its endpoint. The use of C the compound-element block code leads to more general C elementary waves. C C C (4) Generalized layer code C Another code may be obtained from the compound-element C block code if the modified interpretation of the model C structure is used. The modified interpretation consists C in assuming an existence of fictitious parts of blocks of C a zero thickness situated between every neighbouring C complex blocks, the indices of which are not sequentially C ordered. The number of fictitious blocks is considered to C be just the number, which is necessary to fill the gap C between indices of the two neighbouring complex blocks. C The ray can only pass through the fictitious blocks, no C reflection in fictitious blocks being allowed. Similarly, C no conversion is allowed neither at interfaces between C fictitious blocks nor at the interface at which the ray C leaves a fictitious block (conversion may take place only C at the interface at which the ray enters the fictitious C block). If any of these prohibited situations is C specified by the code, it leads to the termination of C computations of the corresponding ray. C C This code is a generalization of the code used in the 2-D C program package 'SEIS83' and is described in Cerveny, C Molotkov and Psencik: Ray Method in Seismology, Charles C University, Prague 1977. It enables any block structure C to be interpreted as locally layered structure with C fictitious layers. Thus, the routines for automatic or C semiautomatic generation of ray codes for layered C structures may be used for this code even in general block C structures. Its effective use is conditioned by proper C indexing of complex blocks, which minimizes the number of C fictitious layers. For example, in a layered structure, C the blocks (layers) should be indexed sequentially from C the top to the bottom. C C C (5) Reflection-transmission code C For any element of a ray one integer is used. The first C element of the ray, i.e. the element containing the C initial point of the ray, is denoted by 1 for P and by -1 C for S wave. Any other element is denoted by 1 or -1 (P or C S wave) if the ray is transmitted at the initial point of C the element, and by 2 or -2 (P or S wave) if the ray is C reflected at the initial point of the element. The code C of a ray is a chain of the above numbers corresponding to C individual elements of the ray from its initial point to C its endpoint. In a layered model, this code is very C convenient for the description of refracted and primary C reflected waves. The application of this code is C straightforward. C C Specification of an elementary wave C an elementary wave is specified by the following data: C (1) Integer KODTYP, which determines the type of the code C 1... Block-surface code, C 2... Basic block code, C 3... Compound-element block code, C 4... Generalized layer code, C 5... Reflection-transmission code. C (2) Integer array KODE containing the code of the elementary C wave. The code is a finite sequence of nonzero integers. C A zero indicates the end of a code in the input data. C The data (1) and (2) should be stored in common block C ------------------------------- C COMMON /COD/ KODTYP,KODE(MKODE), see 'code.inc'. C code.inc C ------------------------------- C The dimension MKODE of the array KODE may be adjusted by the user. C C Date: 1988, June 3 C Written by Vlastislav Cerveny, Ludek Klimes, Ivan Psencik C C======================================================================= C C C SUBROUTINE CODE1(LUN,IWAVE,IWAVE0,IKODE) INTEGER LUN,IWAVE,IWAVE0,IKODE C C This subroutine is called when starting the computation of a new C elementary wave. It stores the code of the elementary wave into the C common block /COD/. C C Input: C LUN... Logical unit number of the external input device C containing the input data. C IWAVE...Zero when starting the complete ray tracing program, C otherwise the index of the last elementary wave which has C been computed (i.e. the output from the previous C invocation of the subroutine CODE1). C C Output: C IWAVE...Zero if all required elementary waves are computed and the C complete ray tracing program will be terminated. C Otherwise, the index of the elementary wave which will be C computed (i.e. NSKIP+1 for the first computed elementary C wave, where NSKIP is the number of elementary waves having C been skipped over, otherwise the input value increased by C one). C IWAVE0..Index of the already computed elementary wave having the C most numerous common elements with the current elementary C wave. In the case of several possibilities, the first C computed wave of them is taken. C IKODE...The length of the common part of the codes of the IWAVE-th C and IWAVE0-th elementary waves. C C Common block /COD/: INCLUDE 'code.inc' C code.inc C All the storage locations of the common block are defined in this C subroutine. C C Subroutines and external functions required: EXTERNAL LCODE INTEGER LCODE C LCODE.. This file. C C Date: 2011, February 4 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Storage for previous codes: INTEGER MKODES PARAMETER (MKODES=1024) INTEGER KODES(0:MKODES) SAVE KODES C C Auxiliary storage locations: CHARACTER*80 TEXTC INTEGER NSKIP,I,J,K,L SAVE NSKIP C C TEXTC...The name of the data. String of 80 characters. C NSKIP...The number of elementary waves to be skipped over. C I,J... Auxiliary loop variables. C K... Auxiliary variable - index in array KODES. C L... Auxiliary variable - length of the current code. C C....................................................................... C IF(IWAVE.EQ.0) THEN C Reading the name of the input data READ(LUN,*) TEXTC C C Reading the number of elementary waves to be skipped over NSKIP=0 READ(LUN,*) NSKIP C No codes of elementary waves are read and stored now KODES(0)=0 END IF C C Reading the code of the elementary wave 10 CONTINUE KODTYP=0 DO 11 I=1,MKODE KODE(I)=0 11 CONTINUE READ(LUN,*) KODTYP,KODE IF(KODTYP.EQ.0) THEN IWAVE=0 RETURN ELSE K=KODES(IWAVE)+1 L=LCODE() IWAVE=IWAVE+1 IF(K+L.GT.MKODES) THEN C 401 CALL ERROR('401 in CODE1: Insufficient memory for codes') C The dimension MKODES of the array KODES(0:MKODES) in C this subroutine should be increased. The dimension MKODES C should at least equal to the number of elementary waves C plus the total number of all indices forming the codes of C the waves. C Note: C If the parameter MODCRT of the data set 'DCRT' (see source C code file 'ray.for') is zero, the output parameters IWAVE0 C and IKODE of this subroutine are likely unused. Then, C instead of increasing MKODES, the above statements C generating this error message may be disabled to save C the computer memory. Then this subroutine will work, C but considering only the first several elementary waves C when calculating output values of IWAVE0 and IKODE. ELSE C Storing the codes DO 12 I=L,1,-1 KODES(K+I)=KODE(I) 12 CONTINUE DO 13 I=K-1,IWAVE,-1 KODES(I+1)=KODES(I) 13 CONTINUE KODES(IWAVE)=K+L DO 14 I=IWAVE-1,0,-1 KODES(I)=KODES(I)+1 14 CONTINUE END IF END IF IF(IWAVE.LE.NSKIP) GO TO 10 C C Determining elementary wave having the most numerous common C elements with the current elementary wave IKODE=0 IWAVE0=0 L=LCODE() DO 39 J=1,MIN0(KODES(0),IWAVE-1) K=KODES(J-1) DO 31 I=1,MIN0(KODES(J)-K,L) IF(KODE(I).NE.KODES(K+I)) THEN GO TO 32 END IF 31 CONTINUE 32 CONTINUE I=I-1 IF(I.GT.IKODE) THEN IKODE=I IWAVE0=J END IF 39 CONTINUE RETURN END C C======================================================================= C C C INTEGER FUNCTION LCODE() C C Integer function LCODE is designed to return the length of the code of C the current elementary wave. C C No input. C C Output: C LCODE...Length of the code of the current elementary wave, i.e. C the count of the numeric items in the array KODE which C describes the behaviour of a ray at interfaces. The array C is called here the code of elementary wave. C C Common block /COD/: INCLUDE 'code.inc' C code.inc C None of the storage locations of the common block are altered. C C No subroutines and external functions required. C C Date: 1996, June 12 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER I C DO 1 I=1,MKODE IF(KODE(I).EQ.0) THEN LCODE=I-1 GO TO 2 END IF 1 CONTINUE LCODE=MKODE 2 CONTINUE RETURN END C C======================================================================= C C C INTEGER FUNCTION NELEM(KODIND) INTEGER KODIND C C This function is designed to return the maximum possible number of ray C elements corresponding to the given position in the code of the C elementary wave. C C Input: C KODIND..Position in the code (index in array KODE) C None of the input parameters are altered. C C Output: C NELEM...Number of ray elements between the initial point of the C ray and the end of the element corresponding to the given C position in the code. Note that, for a particular ray, C some of these possible ray elements may be left out or C coupled into one virtual element. C C Common block /COD/: INCLUDE 'code.inc' C code.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 IF(KODTYP.EQ.1) THEN NELEM=(KODIND+1)/2 ELSE NELEM=KODIND END IF RETURN END C C======================================================================= C C C SUBROUTINE CODE(IY,KODIND,ICBNEW,IEND) INTEGER IY(12),KODIND,ICBNEW,IEND C C This subroutine transforms the used numerical code of the elementary C wave under consideration into instructions specifying the behaviour C of the wave at the initial points of rays and at all points of C incidence of the rays at interfaces (boundaries of complex blocks). C Thus, it is first called at the initial point of a ray and then C successively at all points of incidence. C C Input: C IY... Integer array. Its numerical storage units 2, 3, 5, 6, 8 C must be defined as follows: C IY(2)... Position in the code (index in array KODE) C corresponding to the last computed element of a ray. C IY(2)=0 at the initial point of the ray, C IY(2)=KODIND from the last invocation of subroutine CODE C at other points of the ray. C IY(3)=ICB0... Index of the neighbouring complex block, C from which the ray entered the complex block in which C the last computed element of the ray is situated. C IY(3)=0 before leaving the complex block, in which the C initial point of the ray is situated. C At the initial point of the ray (IY(2)=0), IY(3) is C ignored. C IY(5)=ICB1... Index of the complex block containing the C computed element of the ray, supplemented by the sign C '+' for P wave and sign '-' for S wave. C ICB1 is ignored at the initial point of the ray C (IY(2)=0). C IY(6)=ISRF... Index of the surface at which the endpoint C of the last computed element of the ray is situated. C The sign of IY(6) is ignored. C IY(6)=0 at the initial point of the ray. C IY(8)=ICB2... Index of the complex block touching the C complex block ICB1 from the other side of the surface C ISRF at the endpoint of the last computed element of the C ray. IY(8)=0 for a free space. C At the initial point of the ray, ICB2 is the index of C the complex block containing the initial point. C The input parameter is not altered. C C Output: C KODIND..Position in the code (index in the array KODE) C corresponding to the next element of the ray. C ICBNEW..The index of the complex block in which the next C element of the ray is to be situated, supplemented C by the sign "+" for P wave or "-" for S wave. C IEND... Information on the process of the interpretation C of the code: C IEND.EQ.0... Computation of the ray continues. C IEND.NE.0... The computation of the ray terminates. C Different values of IEND specify the reason for the C termination: C 10... Ray satisfies the whole code (regular termination C of the ray). C 21... The point of incidence is situated at a different C surface than that required by the code. C 22... The next element of the ray is required by the C code to be situated in a complex block that does not C touch the point of incidence. C 23... Transmission is required by the code at a free C surface, or the initial point of the ray is situated C in free space. C 30... Reflection or wave conversion at the fictitious C part of the interface. C C Common block /COD/: C ------------------------------------------------------------------ INCLUDE 'code.inc' C code.inc C ------------------------------------------------------------------ C KODTYP..Determines the type of the code: C 1... Block-surface code, C 2... Basic block code, C 3... Compound-element block code, C 4... Generalized layer code, C 5... Reflection-transmission code. C KODE... Array containing the code of the elementary wave. C The code is a finite sequence of nonzero integers. C Zero indicates the end of a code. C None of the storage locations of the common block are altered. C C No subroutines and external functions required. C C Date: 1990, November 18 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER I,J C KODIND=IY(2) GO TO (10,20,30,30,50),KODTYP C C Block-surface code: 10 CONTINUE IF(KODIND.GT.0) THEN C Check on the surface: KODIND=KODIND+1 IF(KODE(KODIND).EQ.0) THEN C End of code IEND=10 RETURN ELSE IF(KODE(KODIND).NE.IABS(IY(6))) THEN C Wrong surface IEND=21 RETURN END IF END IF C The rest of interpretation is the same as in the basic block code C C Basic block code: 20 CONTINUE KODIND=KODIND+1 ICBNEW=KODE(KODIND) IF(ICBNEW.EQ.0) THEN C End of code IEND=10 ELSE IF(IABS(ICBNEW).EQ.IY(8)) THEN C Transmission IEND=0 ELSE IF(KODIND.GT.1.AND.IABS(ICBNEW).EQ.IABS(IY(5))) THEN C Reflection IEND=0 ELSE C Required complex block is not attainable IEND=22 END IF RETURN C C Compound-element block code: 30 CONTINUE IF(KODIND.GT.0) THEN I=IABS(IY(5)) J=(I-IABS(IY(3)))*(I-IY(8)) DO 31 I=KODIND-1,1,-1 IF(KODE(I).NE.KODE(KODIND)) GO TO 32 J=-J 31 CONTINUE 32 CONTINUE IF(J.GT.0) THEN C Compound element: KODIND=KODIND+1 IF(KODE(KODIND).NE.IY(5)) THEN C Wrong surface IEND=21 RETURN END IF END IF END IF C Rest of interpretation is the same as in Basic block code IF(KODTYP.EQ.3) GO TO 20 C C Generalized layer code: C Interpretation of compound elements is the same as in C Compound-element block code, then: KODIND=KODIND+1 ICBNEW=KODE(KODIND) IF(ICBNEW.EQ.0) THEN C End of code IEND=10 ELSE IF(KODIND.EQ.1) THEN C Initial point of a ray: IF(IABS(ICBNEW).EQ.IY(8)) THEN C Initial point of a ray is situated in the required c.b. IEND=0 ELSE C Initial point of a ray is not situated in the required c.b. IEND=22 END IF ELSE IF(IABS(ICBNEW).EQ.IABS(IY(5))) THEN C Reflection IEND=0 ELSE C Possible transmission: J=ISIGN(1,IY(8)-IABS(IY(5))) IF(IABS(ICBNEW).EQ.IABS(IY(5))+J) THEN C Loop for fictitious parts of blocks: DO 41 I=IABS(IY(5))+J+J,IY(8),J KODIND=KODIND+1 ICBNEW=KODE(KODIND) IF(ICBNEW.EQ.0) THEN C End of code in the fictitious part of block IEND=10 RETURN ELSE IF(ICBNEW.EQ.ISIGN(I,KODE(KODIND-1))) THEN C Transmission from the fictitious part of block IEND=0 ELSE C Termination of the ray computation in the fictitious part C of the block IEND=30 RETURN END IF 41 CONTINUE ELSE C Required complex block is not attainable IEND=22 END IF END IF RETURN C C Reflection-transmission code: 50 CONTINUE KODIND=KODIND+1 I=KODE(KODIND) IF(IABS(I).EQ.1) THEN C Transmission: IF(IY(8).GT.0) THEN C Transmission into material block ICBNEW=ISIGN(IY(8),I) IEND=0 ELSE C Transmission into free space IEND=23 END IF ELSE IF(IABS(I).EQ.2.AND.KODIND.GT.1) THEN C Reflection ICBNEW=ISIGN(IABS(IY(5)),I) IEND=0 ELSE C End of code IEND=10 END IF RETURN END C C======================================================================= C