C 
C File 'ray.for' for the complete tracing of a ray from a given point
C
C Date: 2004, June 11
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 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 data SRFC
C     This input is performed only if NSRFCA.GT.0, see (5).
C Example of 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.
C     YL(2)=VS... Velocity of S-waves at the point.
C     YL(3)=RO... Density at the point.
C     YL(4)=V1,YL(5)=V2,YL(6)=V3... Velocity derivatives in general
C             coordinates.
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)... Travel time.
C     Y(2)... Imaginary part of the complex-valued travel 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 slowness vector.
C     Y(9),Y(10),Y(11)... Covariant components of the polarization
C             vector E1 perpendicular to the ray.
C     ( Y(12),Y(16),Y(20),Y(24) )    Ray propagator matrix (i.e. the
C     ( Y(13),Y(17),Y(21),Y(25) ) ...matrix of fundamental solutions
C     ( Y(14),Y(18),Y(22),Y(26) )    of the dynamic ray tracing system.
C     ( Y(15),Y(19),Y(23),Y(27) )
C     Y(28) TO Y(NY), where NY=27+NAMPL... NAMPL real quantities
C             representing complex-valued vectorial reduced amplitudes.
C             The vectorial reduced amplitudes are specified in the
C             ray-centred coordinate system.
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                                                      
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     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
C Date: 1990, January 23
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(35),YY(5)
      INTEGER IY(12),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,32,33 are indicated in subroutine TRANS
C               (file 'trans.for').
C               The reasons IEND=39,40,50,60 are checked and indicated
C               by IY(6) in subroutine RAYCB (file 'raycb.for').
C               IEND=61 is indicated in this subroutine using IY(6).
C               IEND=71,72,73,74 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             30,32... Reflection or type conversion at the fictitious
C               part of the interface.
C             33... Zero reflection/transmission coefficient.
C             37... 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=37 is reported as an error
C               inside this subroutine RAY2 and should not occur as the
C               output value.
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             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: 2004, April 5
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.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.GE.37) THEN
C         570
          CALL ERROR('570 in RAY2: Wrong function of subroutine CODE')
C         Reason 37 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