C
C Subroutine file 'code.for' - codes for elementary waves.
C
C Date: 2000, May 10
C Coded by Ludek Klimes
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: 1996, September 30
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 PAUSE and STOP
C statements generating this error message may be disabled
C to save the computer memory. Then this subroutine will
C work, but considering only the first several elementary
C waves 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