C File 'ray.for' for the complete tracing of a ray from a given point C C Version: 7.10 C Date: 2014, May 20 C Coded by Ludek Klimes C C....................................................................... C C This file consists of: C Description of the quantities computed along a ray, see C C.R.T.5.2. C RAY1... Subroutine designed to read the input data for complete C ray tracing and to store them in the common blocks /RAYT/ C and /DCRT/, see C C.R.T.5.6. C RAY1 C RAY2... Subroutine designed to continue in the complete ray C tracing of a ray from the given point, see C C.R.T.5.7. C RAY2 C C....................................................................... C C C Input data DCRT for complete ray tracing: 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 TEXTR), the input parameter is of the C type REAL. C (1) TEXTR C String describing the data. Only the first 80 characters of the C string are significant. C (2) KSTORE,NEXPS,NHLF,MODCRT C Quantities controlling numerical integration. C The default values represented by a single slash are usually C satisfactory here. C KSTORE..Specifies whether the conversion coefficients are to be C considered, see C C.R.T.5.5.4 C and 5.6e: C KSTORE.LE.1... No amplitude conversion coefficients are C applied at the point of intersection with an interface. C KSTORE.GE.2... Amplitude conversion coefficients are C applied at the point of intersection with an interface. C Whenever KSTORE is changed, data (8) and (9) should be C adjusted properly. C NEXPS.. Specifies independent variable along rays, see C C.R.T.5.1. C Travel time given by NEXPS=0 (default) is highly C recommended. C NHLF... Maximum allowed number of halvings (bisections) of initial C increment of independent variable during numerical C integration. Default value is NHLF=5, and the authors C have never used different one. C MODCRT..Specifies the mode of complete ray tracing, and the ray C elements along which the computed quantities are stored. C 0...Initial mode, C storing the computed quantities along the whole rays, C no storing at the endpoints of the ray elements. C 1...Initial mode, C storing the quantities just along new ray elements, C no storing at the endpoints of the ray elements. C 2...Initial mode, C storing the quantities just along new ray elements, C storing at the endpoints of the new ray elements. C 3...Interpolation mode, C storing the quantities just along new ray elements, C storing at the endpoints of the new ray elements. C *** modcrt=3 is not enabled by the current version of C 'init.for' ***. C Default: KSTORE=0, NEXPS=0, NHLF=5, MODCRT=0. C C (3) STORE,STEP,UEB,UEBPP,UEBPH,UEBHH,UEBDRT C A single slash representing the default values is usually C satisfactory in place of UEBPP,UEBPH,UEBHH,UEBDRT here. C Quantities controlling numerical integration. C STORE...Step of independent variable for storing the computed C quantities along a ray, see C C.R.T.5.5.1. C For STORE=0 the quantities are not stored along rays. C STEP... Initial increment of independent variable for numerical C integration (expressed in travel-time units if NEXPS=0). C Recommended values: C STEP should never exceed the reciprocal value of the C size of the velocity gradient. C STEP should not exceed the grid interval (or preferably C half the grid interval) for B-spline or similar C interpolation divided by the corresponding velocity. C UEB... Upper error bound of travel time per one step of numerical C integration. Errors in the coordinates of points along C the ray are approximately transformed to units of travel C time and are also bounded by UEB. The error per step of C numerical integration is automatically kept within the C limit UEB if this does not require more than NHLF C bisections of the initial increment step. In the opposite C case the upper error bound is 2,4,8... times greater for C the computation of the rest of the ray. Thus the C computation of each ray is completed. C Recommended values: C UEB from TTERR/SQRT(N) for complex models with extensive C B-spline or similar grids to TTERR/N for simple models, C where N is the number of steps along a ray and TTERR is C the maximum desired travel-time error along the ray. C Here TTERR should be smaller (at least several times) C than the accuracy XERR of the two-point ray tracing C converted to the time units, see input data RPAR-(2). C Description of C RPAR-(2) C UEBPP,UEBPH,UEBHH... Maximum allowed accumulated deviations of the C two polarization vectors from the conditions of C orthonormality, see C C.R.T.5.8.2d C and 5.8.3i. C The accumulated deviations are expressed in time units. C If any of the accumulated deviations overflows its upper C bound, a warning is generated. C UEBDRT..Maximum allowed deviation of L.H.S. and R.H.S. in C.R.T., C eq.(5.11). C If the deviation of any component of the 4*4 matrix C exceeds UEBDRT, a warning is generated. C Default: UEBPP=0., UEBPH=0., UEBHH=0., UEBDRT=0.. C (4) X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,TMAX C The boundaries of the computational volume. C X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX... Coordinates specifying the C coordinate planes bounding the computational volume. The C coordinate planes are indexed 101, 102, 103, 104, 105 and C 106. Default values are the boundaries of the model, see C subroutine MODEL1 and the input data for the model. C TMAX... Maximum travel time for the computation of rays. The C corresponding isochrone is indexed by 107. Default value C is 999999. C (5) NSRFCA C Number of auxiliary surfaces, see C C.R.T.5.3. C The surfaces are indexed sequentially by positive integers C from NSRFC+1 to NSRFC+NSRFCA. Default value is 0. C (6) The indices of end surfaces. C Each index must have the sign corresponding to the side of the C surface at which the computational volume is situated. C The indices may be specified in an arbitrary order and must be C terminated by a slash. See C C.R.T.5.4f. C C (7) The indices of surfaces for storing computed quantities. C The indices of surfaces 1 to NSRFC covering structural interfaces C must include the proper sign. The indices may be specified in an C arbitrary order and must be terminated by a slash. Note that the C same surface may be specified several times, especially if the C quantities at different points of intersections with the surface C are stored in different files, see (8) and (9) below. For details C see C.R.T.5.5.2. C (8) For each above surface, the sequential number of the first point C of intersection of a ray with the surface in which the computed C quantities are stored into the corresponding file. A ray may have C several points of intersection with a specified surface. The C quantities need not be stored at all points of intersection. They C may be stored just from the point of intersection specified in C this input (8) to the one specified in the input (9) below. C Usually, user will wish to store the quantities at all points of C intersection, that corresponds to the values of 1 in this input C and is the default (indicated by a slash in the input data). C For a structural interface (not auxiliary surface), the positive C or negative side of the surface ISRFR is specified, then: C No point of intersection is accounted if the reflected ray is C situated at the other side of the surface ISRFR than the C specified one, and C one point of intersection is accounted if the reflected ray is C situated at the specified side of the surface ISRFR and C KSTORE.GE.2 in the data (2). C Two points of intersection are accounted if the reflected ray is C situated at the specified side of the surface ISRFR and C KSTORE.LE.1 in the data (2). C (9) For each above surface, the sequential number of the last point C of intersection of a ray with the surface in which the computed C quantities are stored into the corresponding file. Default values C of this input are 999999. C (10) The data specifying NSRFCA functions describing the auxiliary C surfaces. The data are read by subroutine SRFC1. For their C description refer to subroutine SRFC1. C Description of C data SRFC. C This input is performed only if NSRFCA.GT.0, see (5). C Example of C data set DCRT. C C....................................................................... C C Storage in the memory: C The input data (2) to (9) are stored in common block /DCRT/ C defined in the include file 'dcrt.inc'. C dcrt.inc C C======================================================================= C C Description of the quantities computed along a ray, see C C.R.T.5.2. C C C Local quantities: C (see C.R.T.5.5.4) C YL(1)=VP... Velocity of P waves at the point in an isotropic C medium. C Probably inconsistent use in an anisotropic medium: C Phase velocity of the P wave corresponding to the C direction of the reference slowness vector, or the P-wave C velocity in the reference isotropic medium. C YL(2)=VP... Velocity of S waves at the point in an isotropic C medium. C Probably inconsistent use in an anisotropic medium: C Common S-wave phase velocity corresponding to the C direction of the reference slowness vector, or the S-wave C velocity in the reference isotropic medium. C YL(3)=RO... Density at the point. C YL(4)=V1,YL(5)=V2,YL(6)=V3... Velocity derivatives in general C coordinates in an isotropic medium. C Probably inconsistent use in an anisotropic medium: C Spatial gradient of the P-wave or of the common S-wave C Hamiltonian function multiplied by the phase velocity C (both corresponding to the direction of the reference C slowness vector, or the velocity gradient in the reference C isotropic medium. C C C Quantities computed along a ray: C (see C.R.T.5.2) C Basic quantities computed along a ray: C (see C.R.T.5.2.1) C Y(1)... Reference travel time: C P-wave travel time or common S-wave travel time. C Y(2)... Imaginary part of the reference complex-valued travel C time. C Y(3),Y(4),Y(5)... Coordinates of points along the ray. C Y(6),Y(7),Y(8)... Covariant components of the reference slowness C vector. C Y(9),Y(10),Y(11)... Covariant components (with respect to the C Riemannian model metric tensor in model coordinates) of C the first contravariant basis vector of the reference C ray-centred coordinate system (Klimes, 2006, sec. 5.4). C The basis vector is normalized with respect to the C Riemannian model metric. C ( Y(12),Y(16),Y(20),Y(24) ) Reference propagator matrix of C ( Y(13),Y(17),Y(21),Y(25) )... geodetic deviation in the reference C ( Y(14),Y(18),Y(22),Y(26) ) ray-centred coordinates. C ( Y(15),Y(19),Y(23),Y(27) ) C Y(28) to Y(NY), where NY=27+NAMPL... NAMPL real quantities C representing the reference complex-valued vectorial C reduced amplitudes. C The reference vectorial reduced amplitudes are specified C with respect to the reference polarization vectors. C In isotropic media, the reference polarization vectors C coincide with the unit basis vectors of the ray-centred C coordinate system. C In anisotropic media, C the P-vave reference polarization vector is equal to the C P-wave eigenvector of the reference Christoffel matrix, C and the S-vave reference polarization vectors are equal C to the S-wave polarization vectors in the low-frequency C limit of the coupling ray theory for S waves, which lay C in the same plane as the S-wave eigenvectors of the C situated reference Christoffel matrix, but they do not C rotate about the common anisotropic reference ray. C Specification of Y(28) to Y(NY): C P wave at the initial point of the ray, C P wave at the point under consideration: C NAMPL=2, C Y(28)=REAL(A33), Y(29)=AIMAG(A33). C P wave at the initial point of the ray, C S wave at the point under consideration: C NAMPL=4, C Y(28)=REAL(A13), Y(29)=AIMAG(A13), C Y(30)=REAL(A23), Y(31)=AIMAG(A23). C S wave at the initial point of the ray, C P wave at the point under consideration: C NAMPL=4, C Y(28)=REAL(A31), Y(29)=AIMAG(A31), C Y(30)=REAL(A32), Y(31)=AIMAG(A32). C S wave at the initial point of the ray, C S wave at the point under consideration: C NAMPL=8, C Y(28)=REAL(A11), Y(29)=AIMAG(A11), C Y(30)=REAL(A21), Y(31)=AIMAG(A21), C Y(32)=REAL(A12), Y(33)=AIMAG(A12), C Y(34)=REAL(A22), Y(35)=AIMAG(A22). C Y(NY+1) to Y(35), where NY=27+NAMPL... Undefined. C Y(36) to Y(41)... Possible S-wave reference polarization vectors, C which are represented by the S-wave polarization vectors C in the low-frequency limit of the coupling ray theory C for S waves. C For an S wave in an anisotropic medium: C Components of the two S-wave polarization vectors. C The reference S-wave polarization vectors are situated C in the same plane as the S-wave eigenvectors of the C Christoffel matrix, but they do not rotate about the C common anisotropic reference ray. C For a P wave, or in an isotropic medium: C Undefined. C Y(42) to Y(44)... Possible P-wave reference polarization vector C coinciding with the P-wave eigenvector of the Christoffel C matrix. C In an anisotropic medium (both P and S waves): C Components of the P-wave eigenvector of the Christoffel C matrix. C In an isotropic medium: C Undefined. C C Auxiliary quantities computed along a ray: C (see C.R.T.5.2.2) C YY(1)...Independent variable along a ray. C YY(2)=UEBRAY... Upper error bound for ray tracing, which is equal C to the input value UEB at the initial point of the ray, C C.R.T.5.6j. C It is always doubled when the numerical integration C requires more than NHLF bisections of the initial C increment step, see C C.R.T.5.6g. C UEBRAY.GT.UEB at the endpoint of the ray indicates C a decreased accuracy of computation. C YY(3)=ERRPP,YY(4)=ERRPH,YY(5)=ERRHH... Deviations (in absolute C values) of the two computed basis vectors of the C ray-centred coordinate system from the conditions of C orthonormality accumulated along the ray. See C C.R.T.5.8.2d C and 5.8.3i. C Any of the deviations may be compared with the C corresponding specified limit UEBPP, UEBPH, UEBHH at the C endpoint of the ray. C C IY(1)=NY=27+NAMPL... Number of the basic quantities C describing the point of a ray, see above. C IY(2)=KODIND... Position in the code (index in array KODE) C corresponding to the considered element of a ray. Its C value is determined by subroutine CODE, see C C.R.T.4. C IY(3)=ICB0... Index of the complex block, from which the ray C entered the complex block in which the computed C 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 IY(4)=ISB1... Index of the simple block containing the C computed element of the ray. C IY(5)=ICB1... Index of the complex block containing the C computed element of the ray, supplemented by a sign '+' C for P wave and sign '-' for S wave. C IY(6)=ISRF... Index of the surface at which the endpoint of the C computed element of the ray is situated, supplemented by C a sign '+' or '-' for the endpoint situated at the C positive or negative side of the surface, respectively. C Zero inside the complex block. Note: the sign of this C quantity is the exception to the definition in the C original paper on C.R.T. C IABS(IY(6))=1,2,...,NSRFC: C One of the surfaces covering structural interfaces. C IABS(IY(6))=NSRFC+1,...,NSRFC+NSRFCA: C One of the auxiliary surfaces serving as the end C surfaces, see C.R.T.5.4f. C IY(6))=101,102,...,106: C The boundaries of the computational volume. C IY(6))=107: C The isochrone corresponding to the maximum travel time C for the computation of rays. C IY(6))=108: C Too small velocity. This usually occurs when a ray C is approaching the region of negative velocity. C Tracing of the ray should be terminated in order to C avoid too many steps of numerical integration (when C integrating with respect to travel time) and consecutive C program failure. C IY(6))=109: C The relative difference between the S-wave eigenvalues C of the Christoffel matrix is smaller than input SEP C parameter DSWAVE when tracing the anisotropic-ray-theory C S-wave rays (input SEP parameter KSWAVE=1 or KSWAVE=2). C IY(7)=ISB2... Index of the simple block touching the complex C block ICB1 from the other side of the surface ISRF at C the endpoint of the computed element of the ray. C ISB2=0 for a free space on the other side of ISRF. C Undefined inside the complex block, defined only at C the endpoint of the element of the ray. C IY(8)=ICB2... Index of the complex block touching the complex C block ICB1 from the other side of the surface ISRF at C the endpoint of the computed element of the ray. C ICB2=0 for a free space on the other side of ISRF. C Undefined inside the complex block, defined only at C the endpoint of the element of the ray. C IY(9)=IFCT... Number of invocations of subroutine FCT evaluating C the right-hand sides of the ordinary differential C equations, along the computed part of the ray. C IY(10)=IOUTP... Number of successful steps of the numerical C integration, along the ray (i.e. the number of invocations C of subroutine OUTP decreased by the number of invocations C of subroutine HPCG). C IY(11)=ITRANS... Number of transformations at an interface (i.e. C the number of invocations of subroutine TRANS). C IY(12)=KMAH... Number of caustic points along the ray (the index C of the ray trajectory). C IY(13)..Total number of the components of the polarization vectors C to be written to the output files, see Y(36) to Y(44). C IY(13)=0: Isotropic medium. C IY(13)=3: P wave in an anisotropic medium. C IY(13)=6: S wave in an anisotropic medium. C C Date: 2014, February 4 C Written by Vlastislav Cerveny, Ludek Klimes, Ivan Psencik C C======================================================================= C C C SUBROUTINE RAY1(LUN) INTEGER LUN C C Subroutine RAY1 reads the input data for complete ray tracing, see C C.R.T.5.6, and stores C them in common blocks /RAYT/ and /DCRT/. It is called once, before C the complete tracing of individual rays begins. Subroutine RAY1 calls C subroutine SRFC1 in order to read the input data for the auxiliary C surfaces, 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 block /MODELC/ (subroutine file 'model.for'): INCLUDE 'model.inc' C model.inc C None of the storage locations of the common block are altered. C C Common block /DCRT/: INCLUDE 'dcrt.inc' C dcrt.inc C All the storage locations of the common blocks are defined in this C subroutine. C C Subroutines and external functions required: INTEGER NSRFC EXTERNAL NSRFC,SRFC1 C NSRFC...File 'model.for' of the package 'MODEL'. C SRFC1 and subsequent routines... File 'srfc.for' and subsequent C files (like 'val.for' and 'fit.for') of the package C 'MODEL'. C C Date: 2004, May 5 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: CHARACTER*80 TEXTR INTEGER I C TEXTR...The name of the data. String of 80 characters. C C (1) String describing the data READ(LUN,*) TEXTR C (2,3) Quantities controlling numerical integration KSTORE=0 NEXPS=0 NHLF=5 MODCRT=0 READ(LUN,*) KSTORE,NEXPS,NHLF,MODCRT UEBPP =0. UEBPH =0. UEBHH =0. UEBDRT=0. READ(LUN,*) STORE,STEP,UEB,UEBPP,UEBPH,UEBHH,UEBDRT C (4) Coordinate planes and isochrone bounding the computational C volume DO 10 I=1,6 BOUNDR(I)=BOUNDM(I) 10 CONTINUE BOUNDR(7)=999999. READ(LUN,*) BOUNDR C (5) Number of auxiliary surfaces NSRFCA=0 READ(LUN,*) NSRFCA C C (6) End surfaces: C Initializing memory for indices of end surfaces DO 11 I=1,MEND KEND(I)=0 11 CONTINUE C Reading indices of end surfaces: READ(LUN,*) (KEND(I),I=1,MEND) DO 12 I=1,MEND IF(IABS(KEND(I)).GT.NSRFC()+NSRFCA) THEN C 563 CALL ERROR('563 in RAY1: Inadmissible end surface') C The index of an end surface in input data is greater than C the number of specified surfaces. ELSE IF(KEND(I).EQ.0) THEN NEND=I-1 GO TO 13 END IF 12 CONTINUE C 561 CALL ERROR('561 in RAY1: Insufficient memory in /DCRT/') C Insufficient memory for the input data in common block C /DCRT/. The dimension MEND of array KEND must be C enlarged. See the include file 'dcrt.inc'. C dcrt.inc 13 CONTINUE C C (7) Surfaces for storing computed quantities: C Initializing memory for indices of the surfaces DO 21 I=1,MSTOR KSTOR(I)=0 21 CONTINUE C Reading indices of the surfaces: READ(LUN,*) (KSTOR(I),I=1,MSTOR) DO 22 I=1,MSTOR/3+1 IF(KSTOR(I).LT.-NSRFC().OR.KSTOR(I).GT.NSRFC()+NSRFCA) THEN IF(KSTOR(I).LT.101.OR.KSTOR(I).GT.107) THEN C 564 CALL ERROR('564 in RAY1: Inadmissible store surface') C The index of a surface for storing the computed quantities C is greater than the number of specified surfaces and while C not specifying the computational volume boundary, or the C index of an auxiliary surface for storing is negative. END IF ELSE IF(KSTOR(I).EQ.0) THEN NSTOR=I-1 GO TO 23 END IF 22 CONTINUE C 562 CALL ERROR('562 in RAY1: Insufficient memory in /DCRT/') C Insufficient memory for the input data in common block /DCRT/. C The dimension MSTOR of array KSTOR must be enlarged. C See the include file 'dcrt.inc'. C dcrt.inc 23 CONTINUE C (8) Sequential numbers of the first points of intersection DO 28 I=NSTOR+1,2*NSTOR KSTOR(I)=1 28 CONTINUE READ(LUN,*) (KSTOR(I),I=NSTOR+1,2*NSTOR) C (9) Sequential numbers of the last points of intersection DO 29 I=2*NSTOR+1,3*NSTOR KSTOR(I)=999999 29 CONTINUE READ(LUN,*) (KSTOR(I),I=2*NSTOR+1,3*NSTOR) C C (10) Auxiliary surfaces: IF(NSRFCA.GT.0) THEN CALL SRFC1(LUN,NSRFCA) END IF C C SEP data: CALL RSEP3R('UEBMUL',UEBMUL,0.) RETURN END C C======================================================================= C C C SUBROUTINE RAY2(YL,Y,YY,IY,IEND) REAL YL(6),Y(44),YY(5) INTEGER IY(13),IEND C C This subroutine continues in the complete tracing of a ray from the C given point, see C C.R.T.5.7. C For MODCRT.LE.2 (see the input data (2)), C the whole ray is computed. For MODCRT.GE.3, just one ray element C is computed. C C Input: C YL... Array containing local quantities at the given point. C Y... Array containing basic quantities computed along the ray C at the given point. C YY... Array containing real auxiliary quantities computed along C the ray at the given point. C IY... Array containing integer auxiliary quantities computed C along the ray at the given point. C C Output: C YL... Array containing local quantities at the endpoint of the C ray. C Y... Array containing basic quantities computed along the ray C at the endpoint of the ray. C YY... Array containing real auxiliary quantities computed along C the ray at the endpoint of the ray. C IY... Array containing integer auxiliary quantities computed C along the ray at the endpoint of the ray. C IEND... Reason of the termination of the computation of a ray, see C C.R.T.5.4. C IEND=10,21,22,23,30 are indicated in subroutine CODE C (file 'code.for'). C IEND=24,25,26,27,32,33 are indicated in subroutine TRANS C (file 'trans.for'). C The reasons IEND=38,39,40,50,60 are checked and C indicated by IY(6) in subroutine RAYCB (file C 'raycb.for'). C IEND=61 is indicated in this subroutine using IY(6). C IEND=71,72,73,74,75 are indicated in subroutine INIT2 C (file 'init.for'). C IEND C 0... The computation of the ray may continue. C 10... Ray satisfies the whole code. C 21... The point of incidence is situated at a different C surface than that specified by the code. C 22... The next element of the ray is required by the code C to be situated in a complex block that does not touch C the point of incidence. C 23,24... Transmission is required by the code at a free C surface. C 25... Ray of the required reflected or transmitted wave C is not real-valued (overcritical reflection or C transmission). C 26... S wave in a liquid block is required by the code. C 27... The next element of the ray is required by the code C to be situated in a complex block that does not touch C the point of incidence. IEND=27 is indicated in C subroutine TRANS (file 'trans.for'), but is reported as C an error inside this subroutine RAY2 and should not C occur as the output value. C 30,32... Reflection or type conversion at the fictitious C part of the interface. C 33... Zero reflection/transmission coefficient. C 38... Relative difference between the S-wave eigenvalues C of the Christoffel matrix is smaller than input SEP C parameter DSWAVE when tracing the anisotropic-ray-theory C S-wave rays (input SEP parameter KSWAVE=1 or KSWAVE=2). C Technical details: C This reason of termination is checked and indicated by C IY(6)=109 on the output of subroutine RAYCB (file C 'raycb.for') called by subroutine RAY2 (this file). C When setting IEND=38 in subroutine RAY2, IY(6) is reset C to zero. C 39... Too small velocity. This usually occurs when a ray C is approaching the region of negative velocity. C Tracing of the ray should be terminated in order to C avoid too many steps of numerical integration (when C integrating with respect to travel time) and consecutive C program failure. C Technical details: C This reason of termination is checked in subroutine C OUTP (file 'raycb.for'), and indicated by IY(6)=108 on C the output of subroutine RAYCB (file 'raycb.for') called C by subroutine RAY2 (this file). When setting IEND=39 in C subroutine RAY2, IY(6) is reset to zero. C 40... Travel time greater than the specified limit. C 50... Ray has intersected a coordinate plane limiting the C computational volume for complete ray tracing. C 60,61... Ray has intersected a smooth surface limiting the C computational volume for complete ray tracing. C 71... There is no ray corresponding to the given ray C parameters. E.g., the given parameters do not belong to C the domain of the initial surface. C 72... Initial point of the ray is not situated in the C required complex block. C 73... Initial point of the ray is not situated in the C computational volume. C 74... Ray of the generated wave is not real-valued. C 75... Relative difference between the S-wave eigenvalues C of the Christoffel matrix at the initial point of the C anisotropic-ray-theory S-wave ray (input SEP parameter C KSWAVE=1 or KSWAVE=2) is smaller than input SEP C parameter DSWAVE. C This reason of termination is checked and indicated by C subroutine INIT2 (file 'init.for'). C Other values of IEND may indicate the program failure. C C Common block /DCRT/: INCLUDE 'dcrt.inc' C dcrt.inc C None of the storage locations of the common block are altered. C C Subroutines and external functions required: EXTERNAL NSRFC,KOOR,CODE,RAYCB,TRANS,CONV,RPAR31,RPAR32,RPAR33 EXTERNAL WRIT31,WRIT32,WRIT33 INTEGER NSRFC,KOOR C NSRFC,BLOCK,VELOC... File 'model.for' of the package 'MODEL'. C KOOR,METRIC... File 'metric.for' of the package 'MODEL'. C SRFC2 and subsequent routines... File 'srfc.for' and subsequent C files (like 'val.for' and 'fit.for') of the package C 'MODEL'. C PARM2 and subsequent routines... File 'parm.for' and subsequent C files (like 'val.for' and 'fit.for') of the package C 'MODEL'. C CDE,CROSS,HIVD2,SMVPRD... File 'means.for' of the package 'MODEL'. C HPCG... IBM Scientific Subroutine Package (file 'hpcg.for' of the C package 'MODEL'). C CODE... File 'code.for'. C RAYCB and subsequent routines... File 'raycb.for'. C TRANS,CONV and subsequent routines... File 'trans.for'. C RPAR31,RPAR32,RPAR33... File 'rpar.for'. C WRIT31,WRIT32,WRIT33 and subsequent routines... File 'writ.for'. C C Date: 2014, February 12 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER KODIND,ICBNEW,I,NAMPL REAL YC(12) C 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 I... Auxiliary loop variable. C NAMPL...Number of real quantities representing complex-valued C vectorial reduced amplitudes. See the subroutine CONV of C the subroutine file 'trans.for'. C YC... Array containing real quantities representing complex- C -valued vectorial reduced amplitudes. See C C.R.T.5.5.4 C and the subroutine CONV of the subroutine file C 'trans.for'. C C....................................................................... C 10 CONTINUE IF(IY(6).EQ.0) THEN C C (a1) Element of the ray inside a complex block CALL RAYCB(YL,Y,YY,IY) C C (b1) Storage of the computed quantities at the interface CALL RPAR31(YL,Y,YY,IY) IF(STORE.NE.0) THEN CALL WRIT31(YL,Y,YY,IY) END IF DO 21 I=1,NSTOR IF(KSTOR(I).EQ.IY(6)) THEN CALL RPAR32(I,YL,Y,YY,IY) CALL CONV(KSTORE,YL,Y,IY,NAMPL,YC) CALL WRIT32(I,YL,Y,YY,IY,NAMPL,YC) END IF IF(KSTOR(I).EQ.-IY(6).AND.KSTORE.GE.2) THEN IY(6)=-IY(6) CALL CONV(KSTORE,YL,Y,IY,NAMPL,YC) CALL WRIT32(I,YL,Y,YY,IY,NAMPL,YC) IY(6)=-IY(6) END IF 21 CONTINUE IF(IABS(IY(6)).LE.NSRFC()) THEN CALL RPAR33(YL,Y,YY,IY) CALL WRIT33(YL,Y,YY,IY) END IF C C (a2) Check for the boundary of the computational volume IF(IABS(IY(6)).GE.109) THEN IY(6)=0 IEND=38 GO TO 99 ELSE IF(IABS(IY(6)).GE.108) THEN IY(6)=0 IEND=39 GO TO 99 ELSE IF(IABS(IY(6)).GE.107) THEN IEND=40 GO TO 99 ELSE IF(IABS(IY(6)).GE.101) THEN IEND=50 GO TO 99 ELSE IF(IABS(IY(6)).GT.NSRFC()) THEN IEND=60 GO TO 99 ELSE IF(MODCRT.GE.3) THEN IEND=0 GO TO 99 END IF C ELSE C C (b2) Interpretation of the elementary waves code CALL CODE(IY,KODIND,ICBNEW,IEND) IF(IEND.NE.0) THEN GO TO 99 END IF C C (b3) Check for the end surfaces C Note: Reflection from an end surface is allowed. IF(IABS(IY(5)).NE.IABS(ICBNEW)) THEN DO 23 I=1,NEND IF(IABS(KEND(I)).EQ.IABS(IY(6))) THEN IEND=61 GO TO 99 END IF 23 CONTINUE END IF C C (b4) Transformation across the interface CALL TRANS(YL,Y,YY,IY,KODIND,ICBNEW,IEND) IF(IEND.EQ.27) THEN C 570 CALL ERROR('570 in RAY2: Wrong function of subroutine CODE') C Reason 27 of the termination of the ray computation C is reported by the subroutine TRANS, although it should be C reported previously by the subroutine CODE as the reason C 22 or 23, respectively. ELSE IF(IEND.NE.0) THEN GO TO 99 END IF C C (b5) Storage of the computed quantities at the interface CALL RPAR31(YL,Y,YY,IY) IF(STORE.NE.0) THEN CALL WRIT31(YL,Y,YY,IY) END IF DO 25 I=1,NSTOR IF(KSTOR(I).EQ.IY(6)) THEN CALL RPAR32(I,YL,Y,YY,IY) IF(KSTORE.LE.1) THEN CALL CONV(KSTORE,YL,Y,IY,NAMPL,YC) CALL WRIT32(I,YL,Y,YY,IY,NAMPL,YC) END IF END IF 25 CONTINUE C C (b6) The ray may continue in the next complex block IY(6)=0 C END IF GO TO 10 C 99 CONTINUE RETURN END C C======================================================================= C