C                                                            <HTML><PRE>
C Subroutine file 'rpar.for' to control the take-off parameters of rays.
C
C Date: 2006, October 8
C Coded by Ludek Klimes and Petr Bulant
C
C.......................................................................
C
C This file consists of the following subroutines:
C     RPAR1...Subroutine designed to read the input data for the
C             take-off parameters of rays and to store them in common
C             block /RPARD/, to prepare the parameters for the
C             subroutine RPAR2 and to store them in the common block
C             /RPARC/.  It is called when starting the computation of a
C             new elementary wave.
C             <A HREF="#RPAR1">RPAR1</A>
C     RPAR2...Subroutine designed to specify the take-off parameters of
C             individual rays.
C             <A HREF="#RPAR2">RPAR2</A>
C     RPAR30..Subroutine performing the numerical quadrature of the
C             squared matrix of geometrical spreading, useful for the
C             estimination of the L2 norm of the distance between
C             neighbouring rays in order to control the fatness of ray
C             tubes.  This subroutine is designed to be called by
C             subroutine PSHIFT of file 'raycb.for'.  Subroutine PSHIFT
C             is called both with steps STEP and STORE along the ray.
C             <A HREF="#RPAR30">RPAR30</A>
C     RPAR31..Subroutine called with constant step STORE of the
C             independent variable along the ray, and at the points of
C             intersection with interfaces either before and after the
C             transformation, i.e. called simultaneously with the
C             subroutine WRIT31, see also
C             <A HREF="../../../papers/b88vc200.htm#28">C.R.T.5.5.1</A>.
C             Subroutime RPAR31 accumulates the ray histories and
C             collects other quantities useful for the determination of
C             the take-off parameters in the memory.
C             <A HREF="#RPAR31">RPAR31</A>
C     RPAR32..Subroutine called at the points of intersection of the ray
C             with specified interfaces, i.e. called simultaneously with
C             the subroutine WRIT32, see also
C             <A HREF="../../../papers/b88vc200.htm#29">C.R.T.5.5.2</A>.
C             It may be used to keep the quantities useful for the
C             determination of the take-off parameters in the memory.
C             <A HREF="#RPAR32">RPAR32</A>
C     RPAR33..Subroutine called at the endpoints of the elements of
C             rays, situated at structural interfaces, i.e. called
C             simultaneously with the subroutine WRIT33, see also
C             <A HREF="../../../papers/b88vc200.htm#29">C.R.T.5.5.3</A>.
C             It may be used to keep the quantities useful for the
C             determination of the take-off parameters in the memory.
C             <A HREF="#RPAR33">RPAR33</A>
C     RPAR4...Subroutine called after finishing the computation of each
C             ray.  It defines four storage locations of the common
C             block /INITC/ introduced in the file 'init.for'.
C             <A HREF="#RPAR4">RPAR4</A>
C     XFUN... Subroutine returning the value and first derivatives of
C             the Xi functions describing the profiles.  It is called by
C             the subroutine RPAR32, and may be user-defined.  In this
C             version, it calls the subroutine SRFC2 from the subroutine
C             file 'srfc.for', or directly returns one of the model
C             coordinates.
C             <A HREF="#XFUN">XFUN</A>
C     IHIST...Auxiliary integer function to the subroutine RPAR4 to
C             determine the ray-history index.  The ray-history index is
C             used to distinguish the rays with different histories for
C             the purposes of the subroutine RPAR2.
C             <A HREF="#IHIST">IHIST</A>
C     SRPARH..Auxiliary subroutine to the subroutine RPAR31 to shift the
C             storage locations of the common block /RPARH/.
C             <A HREF="#SRPARH">SRPARH</A>
C
C Note:
C     The lines denoted by '*' in the first column contain an undebugged
C     code intended for future extensions.
C
C.......................................................................
C
C Description of data files:
C                                                    <A NAME="data"></A>
C Input data RPAR for the take-off parameters of rays:
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 TEXTP), the input parameter is of the
C     type REAL.
C (1) TEXTP
C     String describing the data.  Only the first 80 characters of the
C     string are significant.
C                                                   <A NAME="data2"></A>
C (2) ISRFR,ISRFX1,ISRFX2,XERR,AERR,PRM0(I),I=1,MPRM0,/
C     ISRFR...Index of the reference surface at which the computed
C             quantities are stored and the given profiles or receivers
C             for two-point ray tracing are situated.  The profiles are
C             the crosssections of the surface ISRFR and the isosurfaces
C             of the X1 function corresponding to given values (profile
C             coordinates).  The receivers are the crosssections of the
C             surface ISRFR and the isosurfaces of the X1 and X2
C             functions corresponding to given pairs of values (receiver
C             coordinates).
C             ISRFR should be one of the indices listed in the input
C             data DCRT-(7) (the indices for storing computed
C             quantities), including the sign, see the subroutine file
C             'ray.for':
C             <A HREF="ray.for#data7">Description of data DCRT-(7)</A>
C             For ISRFR=0 there are no profiles for boundary-value ray
C               tracing, ISRFR not listed in the input data DCRT-(7) has
C               the same effect.
C             ISRFR=101,102,103,104,105, or 106 specify the boundary of
C               the computational volume.
C             Example: if the profile for the boundary-value ray tracing
C               is situated at the model (earth's) surface, ISRFR is the
C               index of the model surface.  Its sign must correspond to
C               the side of the surface ISRFR facing the material block.
C     ISRFX1..Index of the function that becomes to be the X1 function.
C             The X1 function is defined here by the subroutine XFUN.
C             <A HREF="#XFUN">XFUN</A>
C             ISRFX1=0:  No boundary-value (two-point) ray tracing and
C               no determination of the boundaries between the ray
C               histories.
C               The X1 function is zero and has no meaning.
C               This option may be used only if the boundary-value rays
C               are not searched for, i.e. if ISRFR.EQ.0 or REC.EQ.' '
C               or IPOINT=0 at all input data lines (4.1).
C             For ISRFX1=-1, -2, or -3,  the X1 function coincides with
C               the -ISRFX1-th model coordinate.
C             For ISRFX1 positive,  the subroutine XFUN is an interface
C               subroutine to the subroutine SRFC2, called with value
C               ISRFX1 assigned to the dummy argument ISRFC, see the
C               subroutine file 'srfc.for'.
C               <A HREF="../model/srfc.for">srfc.for</A>
C               ISRFX1 must not exceed the number of the defined
C               surfaces more than by 1.
C             Hint: Even if there is no two-point ray tracing, specify,
C               e.g., ISRFX1=-1 (default).
C     ISRFX2..Index of the function that becomes to be the X2 function.
C             ISRFX2=0: Initial-value ray tracing (for ISRFX1=0) or
C               one-parametric shooting to profiles corresponding to
C               isosurfaces of the X1 function.  The one-parametric
C               shooting to profiles reduces to the 2-D two-point ray
C               tracing in 2-D models in which the receivers at the
C               reference surface are determined by a single coordinate
C               (X1 function).
C             Otherwise:  Two-parametric 3-D two-point ray tracing to
C               the receivers defined as the crosssections of the
C               surface ISRFR and the isosurfaces of the X1 and X2
C               functions corresponding to given pairs of values
C               (receiver coordinates).
C               Meaning of ISRFX2 is then analogous to ISRFX1:
C               For ISRFX2=-1, -2, or -3,  the X2 function coincides
C                 with the -ISRFX2-th model coordinate.
C               For ISRFX2 positive, the subroutine XFUN is an interface
C                 subroutine to the subroutine SRFC2, called with value
C                 ISRFX2 assigned to the dummy argument ISRFC, see the
C                 subroutine file 'srfc.for'.
C                 <A HREF="../model/srfc.for">srfc.for</A>
C                 ISRFX2 must not exceed the number of the defined
C                 surfaces (including ISRFX1) more than by 1.
C             Hint: Even if there is no two-point ray tracing, specify,
C               e.g., ISRFX2=-2 (default) for two-parametric shooting.
C     XERR... Specifies the accuracy of the boundary-value ray tracing.
C             If the distance of the ray from the given profile or
C             receiver does not exceed XERR, the ray is considered to be
C             a boundary-value ray.  The distance of the ray from the
C             given profile is measured by means of the difference of
C             the X1 function.  Similarly, the distance of the ray from
C             the given receiver is measured in units of the X1 and X2
C             functions (e.g., length units for negative both ISRFX1 and
C             ISRFX2).
C             XERR should be reasonably greater than the maximum error
C             along any ray, converted to the units of XERR.
C             For the maximum error UEB per one step of numerical
C             integration refer to input data DCRT-(3) in 'ray.for'.
C             <A HREF="ray.for#data3">Description of data DCRT-(3)</A>
C             For two-parametric control of the take-off ray parameters
C             also rays distant more than XERR may be considered, see
C             parameter PRM0(5) below.
C     AERR... Specifies the accuracy of the determination of the
C             boundary rays between the different ray histories.
C             The distance between rays in the ray-parameter domain
C             is measured with respect to the metric tensor determined
C             by the input data (4) described below.  The distance
C             between two consecutive rays of different histories
C             (boundary rays) is bisectioned until it does not exceed
C             AERR.  Roughly speaking, the width of a demarcation belt
C             between different ray histories should not exceed AERR
C             times the size of the basic triangles in the ray-parameter
C             domain.
C             No boundary rays are determined for AERR.GT.1 (e.g.
C             AERR=999999.).  Note that IPOINT=0 has the same effect as
C             AERR.GT.1.
C     PRM0(1)... Maximum allowed distance of the boundary ray from the
C             shadow zone (measured on the reference surface).
C             The distance is not measured if PRM0(1)=0 (default).
C             Used only by two-parametric control of the take-off ray
C             parameters.
C     PRM0(2)... Maximum allowed length of sides of homogeneous
C             triangles (measured on the reference surface).
C             The triangles are not measured if PRM0(2)=0 (default).
C             Used only by two-parametric control of the take-off ray
C             parameters.
C     PRM0(3)... Controls the determination of the ray histories.
C             PRM0(3).EQ.0:
C               Rays terminating at different boundaries of the
C               computational volume have different histories.
C               This option is suitable for two-point ray tracing.
C               PRM0(3).EQ.0 corresponds to CRT version 5.00 and is the
C               default in version 5.10.
C             PRM0(3).NE.0 (e.g., PRM0(3)=1):
C               Only structural interfaces and the reference surface
C               (if specified) influence ray histories of unsuccessful
C               rays.  The boundaries of the computational volume and
C               end surfaces are not taken into account.  This option
C               limits the number of thin spaces, uncovered by ray
C               tubes, situated between ray tubes corresponding to
C               different ray histories.
C               This option is thus suitable for controlled
C               initial-value ray tracing followed by the travel-time
C               interpolation within ray tubes.  In such a case, it is
C               also recommended to set IPOINT=999999 (default) in data
C               (4.1) below in order to disable successful rays.
C               PRM0(3).EQ.1 corresponds to CRT version 4.10,
C             PRM0(3).GE.1.5 (e.g., PRM0(3)=2):
C               Index of the surface, at which the ray terminates, is
C               removed from ray histories of unsuccessful rays.
C               This option further limits the number of uncovered thin
C               spaces between ray tubes corresponding to different ray
C               histories and is thus suitable for controlled
C               initial-value ray tracing followed by the travel-time
C               interpolation within ray tubes.
C               It is recommended to set IPOINT=999999 (default) in data
C               (4.1) below in order to disable successful rays.
C             PRM0(3).GE.2.5 (e.g., PRM0(3)=4):
C               Reason of the termination of the computation of a ray is
C               removed from ray histories.  This option further limits
C               the number of uncovered thin spaces between ray tubes
C               corresponding to different ray histories and is thus
C               suitable for controlled initial-value ray tracing
C               followed by the travel-time interpolation within ray
C               tubes.
C               There are 2 suboptions:
C               PRM0(3).LT.3.5 (e.g., PRM0(3)=3):
C                 Distinguishing overcritical incidence (IEND=25) from
C                 other reasons of the termination of the computation of
C                 a ray (recommended for interpolations within ray
C                 cells).
C               PRM0(3).GE.3.5 (e.g., PRM0(3)=4):
C                 Reason of the termination of the computation of a ray
C                 is removed from ray histories at all.
C     PRM0(4)... Controls the maximum length of the sides of homogeneous
C             triangles in the ray-tube metric. The ray-tube
C             metric contains an information about the geometrical
C             spreading corresponding to the ray. Maximum length
C             of the homogeneous triangle's sides in the
C             ray-tube metric thus controls a width of the
C             corresponding ray tube.
C             PRM0(4) is thus intended to be used for controlled
C             initial-value ray tracing followed by the interpolation
C             within ray tubes.
C             There are 3 options:
C             PRM0(4)=0: Reccommended for two-point ray tracing.
C               Measuring of triangles in the ray-tube metric is
C               disabled.
C             PRM0(4).GT.0: Reccommended for controlled initial-value
C               ray tracing. Maximum length of the sides of triangles
C               in the ray-tube metric is PRM0(4), the metric
C               describes an average squared thickness of ray tubes.
C               This means, that in a homogeneous model, e.g., the
C               maximum thickness of the ray tubes will be about
C               PRM0(4)*SQRT(3).
C             PRM0(4).LT.0: Maximum length of the sides of triangles
C               in the ray-tube metric is ABS(PRM0(4)), the metric
C               describes a thickness of ray tubes at their ends. This
C               may be useful, e.g., for controlled initial-value
C               ray tracing of a weakly refracted wave.
C             Used only by two-parametric control of the take-off ray
C             parameters.
C     PRM0(5)... Controls the maximum distance of the boundary-value
C             rays from the given receivers.  If the program fails
C             to find a ray distant less then XERR, the nearest ray
C             is found.  If the nearest ray is distant less than
C             PRM0(5), it is considered to be a two-point ray and
C             a warning is generated.
C             PRM0(5), if specified, should be positive and greater
C             than XERR.  The default PRM0(5)=0. is understood as
C             infinite PRM0(5), and the nearest ray is always
C             considered to be a two-point ray.
C             Used only by two-parametric two-point ray tracing.
C     PRM0(I),I=6,MPRM0... Reserved for future extension of the
C             numerical parameters for two-point ray tracing.
C Default values:
C     ISRFR=0, ISRFX1=-1, ISRFX2=-2, XERR=0.000, AERR=999999,
C     PRM0(I)=0,I=1,MPRM0.
C             The boundary rays are not determined if AERR.GT.1.
C               Neither are they if IPOINT=0, see input data (4.1).
C             The boundary-value rays are not determined if ISRFR.EQ.0,
C               if ISRFR is not listed in the input data DCRT-(7),
C               or if REC.EQ.' '.
C               Neither are they if IPOINT=0, see input data (4.1).
C             Thus, just the initial-value ray tracing is performed if
C               (ISRFR.EQ.0 or ISRFR not in DCRT-(7) or REC.EQ.' ')
C               and AERR.GT.1.
C               Just the initial-value ray tracing is also performed if
C               IPOINT=0, see input data (4.1).
C               During the initial-value ray tracing, only the basic net
C               of rays with the given constant steps is generated.
C (3) The data specifying the functions No. ISRFX1 and ISRFX2.  This
C     auxiliary function is specified by the subroutine SRFC2 and the
C     data are read by subroutine SRFC1.  For the description of the
C     data refer to the subroutine SRFC1 of the file 'srfc.for'.  This
C     input is performed only if ISRFX1 or ISRFX2 is positive and the
C     function ISRFX1 or ISRFX2 is not defined in the data sets
C     <A HREF="crtin.for#data">MODEL or DCRT</A>.
C     In such a case, the index of the first
C     undefined function of ISRFX1 or ISRFX2 must equal the number NSRFC
C     of the surfaces defined in the data set MODEL + the number
C     ISRFCA of the auxiliary surfaces defined in the data set
C     DCRT + 1, and so on.
C     <A HREF="../model/srfc.for#data">Data for 'srfc.for'</A>
C     Example: Imagine that the profile will be constructed as the
C             intersection of the surface ISRFR with an auxiliary
C             vertical plane.  Using the basic version of the file
C             'srfc.for', the auxiliary vertical plane passing through
C             the points (X1,X2,X3)=(X1A,X2A,X3A) and
C             (X1,X2,X3)=(X1B,X2B,X3B) may be defined using subroutine
C             file 'srfc.for' by the following input data:
C               For X1A.NE.X1B:
C                   1 -2 0 0
C                   2
C                   X1A X1B
C                   X2A X2B
C                 The interpolated function then returns the oriented
C                 distance from the vertical plane, measured in the
C                 direction of the X2-axis but with opposite sign.
C               For X2A.NE.X2B:
C                   2 -1 0 0
C                   2
C                   X2A X2B
C                   X1A X1B
C                 The interpolated function then returns the oriented
C                 distance from the vertical plane, measured in the
C                 direction of the X1-axis but with opposite sign.
C                                                   <A NAME="data4"></A>
C (4) For each elementary wave IWAVE the following data (4.1) and (4.2):
C (4.1) Any times:
C     PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM,IPOINT,
C     PRM1(I),I=1,MPRM1,/
C     PAR1L,PAR2L... Ray parameters corresponding to the origin (0,0) of
C             the domain of the normalized ray parameters.
C             The kind of ray parameters used is given by index INIPAR
C             read by <A HREF="init.for">init.for</A> and described in
C             <A HREF="crtin.for#INIT">crtin.for</A>.
C             Imagine the normalized ray parameters as two coordinates
C             ranging from 0 to 1.  The normalized ray-parameter domain
C             is then a unit square.  Normalized ray parameters (0,0)
C             are mapped to ray parameters (PAR1L,PAR2L), normalized ray
C             parameters (1,0) are mapped to ray parameters
C             (PAR1A,PAR2A) and normalized ray parameters (0,1) are
C             mapped to ray parameters (PAR1B,PAR2B).  The mapping is
C             linear.  Any couple of normalized ray parameters can thus
C             be simply mapped to the ray parameters (e.g., to the
C             colatitude and longitude).
C     PAR1A,PAR2A... Ray parameters corresponding to the right bottom
C             corner (1,0) of the domain of the normalized ray
C             parameters.
C     PAR1B,PAR2B... Ray parameters corresponding to the left top corner
C             (0,1) of the domain of the normalized ray parameters.
C     ANUM,BNUM... Specify the 2*2 metric tensor G in the normalized ray
C             parameter domain.
C             The unit square of the normalized ray parameters is
C             divided into ANUM times BNUM rectangles.
C             For initial-value ray tracing or one-parametric shooting,
C             the normalized ray parameters corresponding to the
C             vertices of the rectangles are mapped to the ray
C             parameters and the corresponding basic rays are traced.
C             For two-parametric shooting, the size of the rectangles
C             define the metric for the coverage of the ray-parameter
C             domain by a triangularized network of basic rays.
C     ANUM,BNUM in RP2D (Initial-value ray tracing or one-parametric
C             shooting):
C             The components of the metric tensor are
C             G11=ANUM*ANUM, G12=0., G22=BNUM*BNUM.  The system of basic
C             rays forms a regular rectangular grid in the normalized
C             ray parameter domain.  The grid cells are unit squares
C             with respect to the above metric:
C             For ANUM.GT.0 and BNUM.GT.0: the two-parametric basic net
C               of INT(ANUM)*INT(ANUM) rays is given by the take-off ray
C               parameters
C                 PAR1=PAR1L+A*(PAR1A-PAR1L)/ANUM+B*(PAR1B-PAR1L)/BNUM,
C                 PAR2=PAR2L+A*(PAR2A-PAR2L)/ANUM+B*(PAR2B-PAR2L)/BNUM,
C               where
C                 A=0,1,2,...,INT(ANUM) and B=0,1,2,...,INT(BNUM).
C               The basic net of rays corresponds to the initial-value
C               ray tracing.
C               Besides the basic net of rays, also additional rays are
C               computed in order to find the boundary rays between the
C               bundles of rays with different histories, and to find
C               the boundary-value rays.  All additional rays are given
C               by
C                 B=0,1,2,...,INT(BNUM)
C               and real A from the interval 0 to INT(ANUM), forming
C               the cuts of the ray-parameter surface given by fixed B
C               and parametrized by A.  The boundary rays and the
C               boundary-value rays are searched for along the mentioned
C               cuts.  Under the boundary-value rays we understand here
C               the rays from the above mentioned cuts passing through
C               the XERR-vicinity of the given profiles.
C             For ANUM.GT.0 and BNUM.EQ.0: The one-parametric basic net
C               of INT(ANUM) rays is given by the take-off ray
C               parameters
C                 PAR1=PAR1L+A*(PAR1A-PAR1L)/ANUM,
C                 PAR2=PAR2L+A*(PAR2A-PAR2L)/ANUM,
C               where
C                 A=0,1,2,...,INT(ANUM).
C               Additional rays are computed as for ANUM.GT.0 and
C               BNUM.GT.0.
C             For ANUM.EQ.0 and BNUM.GT.0: Similarly for ANUM.GT.0 and
C               BNUM.EQ.0, but no additional rays are computed.
C             For ANUM.EQ.0 and BNUM.EQ.0: Just a single ray given by
C               PAR1=PAR1L and PAR2=PAR2L is computed.
C     ANUM,BNUM in RP3D (two-parametric shooting):
C             ANUM,BNUM must be positive for two-parametric shooting.
C             The components of the metric tensor are
C             G11=ANUM*ANUM*E11, G12=ANUM*BNUM*E12, G22=BNUM*BNUM*E22,
C             where E is some basic metric tensor, e.g. the unit matrix
C             or the Cartesian metric tensor on the unit sphere around
C             the source.  The system of basic rays forms basic
C             triangles covering the normalized ray parameter domain.
C             The basic triangles are constructed in a way to be roughly
C             equilateral with unit sides with respect to metric G.
C             Recommendations for the choice of ANUM and BNUM:
C               For strong lateral heterogeneities, the ratio of ANUM
C               and BNUM should roughly correspond to the ratio of the
C               lengths of sides of the ray-parameter domain.
C               For moderate lateral heterogeneities, the number
C               corresponding to the longitude in spherical coordinates
C               may be smaller.
C               The initial value of ANUM*BNUM may be chosen, e.g.,
C               between 20*20 to 40*40, or roughly corresponding to
C               ANUM*BNUM=(1/AERR)**2 for reasonable values of AERR.
C               ANUM and BNUM may then be adjusted according to the
C               complexity of the ray tracing.
C               If ANUM*BNUM is great, especially if it is approaching
C               30 or even 50 per cent of the number of all rays traced
C               for the elementary wave, ANUM and BNUM should be
C               decreased.
C               If ANUM*BNUM is small with respect to the number of the
C               rays traced for the elementary wave, ANUM and BNUM may
C               be increased.
C     IPOINT..IPOINT=0:  No ray history is recorded, all rays have the
C               same index 0 and no boundary-value ray tracing is
C               performed.
C             Otherwise:  The IABS(IPOINT)-th point of intersection of
C               a ray with the surface ISRFR is considered for the
C               boundary-value ray tracing.
C             IPOINT negative:  Only those rays satisfying the whole
C               code (IEND=10) are considered to be successful.  All
C               other rays are rejected (IMPORTANT OPTION).
C               Another effect of this option:
C                 If a ray terminates at a coordinate plane limiting the
C                 computational volume (IEND=50), subroutine CODE is
C                 called once again to check whether the ray satisfies
C                 the whole code.
C                 If yes, IEND=50 is replaced with IEND=10.
C             IPOINT=large (e.g. IPOINT=999999):
C               Ray history is recorded, boundary rays may
C               be determined but no boundary-value ray tracing is
C               performed because all rays are unsuccessful.
C               This is the recommended option for controlled
C               initial-value ray tracing followed by the travel-time
C               interpolation within ray tubes.
C             In most cases, user will wish the value of IPOINT=1.
C             For the structural interface (not auxiliary surface), the
C             positive or negative side of the surface ISRFR is
C             specified, then:
C               No point of intersection is accounted if the reflected
C                 ray is situated at the other side of the surface ISRFR
C                 than the specified one, and
C               two points of intersection are accounted if the
C                 reflected ray is situated at the specified side of the
C                 surface ISRFR.
C             Just the initial-value ray tracing is performed if
C             IPOINT=0.
C     PRM1(I),I=1,MPRM1... Reserved for future extension of the
C             numerical parameters for two-point ray tracing, specific
C             to the calculated elementary wave.
C Default values:
C     PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM
C             Zeroes when starting the complete ray tracing program,
C             otherwise the previous values.
C     IPOINT=999999 when starting the complete ray tracing program,
C             otherwise the previous value.
C     PRM1(I)=0, I=1,MPRM1.
C (4.2) /(a slash)
C <A HREF="../data/len/len-rpar.dat">Example of data RPAR</A> for
C     initial-value ray tracing.
C <A HREF="../data/len/len-rpa.dat" >Example of data RPAR</A> for 3-D
C     two-point ray tracing.
C
C                                                     <A NAME="REC"></A>
C Input file 'REC' containing profile (or receiver) coordinates:
C (1) Several strings terminated by / (a slash).
C             The simplest way is to submit just the /.
C (2) Several times (2.1):
C (2.1) 'NAMER(IR)',XR(1,IR),XR(2,IR),XR(3,IR),XR(4,IR),XR(5,IR),/
C     'NAMER(IR)'... String reserved for the name of the receiver.
C             No meaning here, anything in apostrophes may be submitted.
C     XR(1,IR),XR(2,IR),XR(3,IR)... Model coordinates of the receiver.
C             For ISRFX1.LT.0:  XREC1(IR)=XR(-ISRFX1,IR) is the first
C               coordinate of the point along the reference surface.
C             For ISRFX2.LT.0:  XREC2(IR)=XR(-ISRFX2,IR) is the second
C               coordinate of the point along the reference surface.
C             Values of other model coordinates have no meaning here.
C     XR(4,IR),XR(5,IR)... Coordinates of the receiver along the
C             reference surface other than the model coordinates.
C             For ISRFX1.GT.0:  XREC1(IR)=XR(4,IR) is the first
C               coordinate of the point along the reference surface.
C             For ISRFX2.GT.0:  XREC2(IR)=XR(5,IR) is the second
C               coordinate of the point along the reference surface.
C             Other values have no meaning here and need not be
C               specified.
C     Default: X1R(IR)=0, X2R(IR)=0, X3R(IR)=0, X4R(IR)=0, X5R(IR)=0.
C (3) / (a slash) or end of file.
C <A HREF="../data/len/len-rec.dat">Example of data file REC</A>
C
C.......................................................................
C
C Storage in the memory:
C     The input data (2) are stored in common block /RPARD/ defined
C     in the include file 'rpard.inc'.
C     <A HREF="rpard.inc">rpard.inc</A>
C     Common block /RPARC/ to collect the quantities needed for
C     boundary-value ray tracing and common block /RPARH/ describing the
C     histories of rays are defined in the include file 'rparc.inc'.
C     <A HREF="rparc.inc">rparc.inc</A>
C     The histories of rays are used by function IHIST to determine the
C     ray-history index. Different values of ray-history index are
C     attached to rays with different histories.
C
C=======================================================================
C
C     <A NAME="RPAR1"></A>
C
      SUBROUTINE RPAR1(LUN,IWAVE)
      INTEGER LUN,IWAVE
C
C This subroutine is called when starting the computation of a new
C elementary wave.  It prepares the parameters for the subroutine RPAR2
C to specify the take-off parameters of the individual rays.
C
C Input:
C     LUN...  Logical unit number of the external input device
C             containing the input data.
C     IWAVE...Index of the elementary wave that is going to be computed,
C             i.e. the output from the last invocation of the subroutine
C             CODE1.
C
C No output.
C
C Common block /DCRT/:
      INCLUDE 'dcrt.inc'
C     <A HREF="dcrt.inc">dcrt.inc</A>
C None of the storage locations of the common block are altered.
C
C Common block /RPARD/:
      INCLUDE 'rpard.inc'
C     <A HREF="rpard.inc">rpard.inc</A>
C Nearly all the storage locations of the common block are defined in
C this subroutine.
C
C Common block /RPARC/:
      INCLUDE 'rparc.inc'
C     <A HREF="rparc.inc">rparc.inc</A>
C None of the storage locations of the common block, except LURPAR, are
C altered.  LURPAR is set to LUN.
C
C Subroutines and external functions required:
      EXTERNAL NSRFC,RSEP3T
      INTEGER NSRFC
C     NSRFC... File 'model.for' of the package 'MODEL'.
C
C Date: 1999, September 23
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER LUREC
      PARAMETER (LUREC=9)
      CHARACTER*80 TEXTP
      INTEGER I,J
      REAL AUX(5)
C
C     LUREC...Logical unit number of the receiver file.  The file is
C             opened and closed here.
C     TEXTP...The name of the data. String of 80 characters.
C     I,J...  Auxiliary loop variables, or temporary variables.
C     AUX...  Temporary storage locations to read the receiver
C             coordinates.
C
C.......................................................................
C
      IF(IWAVE.EQ.0) THEN
C
C       (1) Reading the name of the input data:
        READ(LUN,*) TEXTP
C       (2):
        ISRFR=0
        ISRFX(1)=-1
        ISRFX(2)=-2
        XERR=0.000
        AERR=999999.
        DO 10 I=1,MPRM0
          PRM0(I)=0.
   10   CONTINUE
        READ(LUN,*) ISRFR,ISRFX(1),ISRFX(2),XERR,AERR,PRM0
C
C       (3) Auxiliary functions no. ISRFX(1) and ISRFX(2):
        IF(ISRFX(1).EQ.0.AND.ISRFX(2).NE.0) THEN
C         <A NAME="error642"></A><FONT COLOR="RED">642</FONT>
          CALL ERROR('642 in RPAR1: Zero ISRFX1, nonzero ISRFX2')
C         For initial-value ray tracing given by ISRFX1=0, there must be
C         ISRFX2=0 in the input data set RPAR.
        END IF
        I=NSRFC()+NSRFCA+1
        IF(ISRFX(1).EQ.I) THEN
          CALL SRFC1(LUN,1)
          I=I+1
        ELSE IF(ISRFX(1).LT.-3.OR.ISRFX(1).GT.I) THEN
C         <A NAME="error643"></A><FONT COLOR="RED">643</FONT>
          CALL ERROR('643 in RPAR1: Wrong value of ISRFX1')
C         The function no. ISRFX1 specified by the routines of the file
C         'srfc.for' is defined neither in the input data set
C         <A HREF="crtin.for#data">MODEL nor DCRT</A> and cannot be
C         defined in this data set.
C         The index ISRFX1 is too large or is less than -3.
        END IF
        IF(ISRFX(2).EQ.I) THEN
          CALL SRFC1(LUN,1)
        ELSE IF(ISRFX(2).LT.-3.OR.ISRFX(2).GT.I) THEN
C         <A NAME="error644"></A><FONT COLOR="RED">644</FONT>
          CALL ERROR('644 in RPAR1: Wrong value of ISRFX2')
C         The function no. ISRFX2 specified by the routines of the file
C         'srfc.for' is defined neither in the input data set
C         <A HREF="crtin.for#data">MODEL nor DCRT</A> and cannot be
C         defined in this data set.
C         The index ISRFX2 is too large or is less than -3.
        END IF
C
C       Default values of data (5):
        PAR1L=0.
        PAR2L=0.
        PAR1A=0.
        PAR2A=0.
        PAR1B=0.
        PAR2B=0.
        ANUM=0.
        BNUM=0.
        IPOINT=999999
C
C       Receivers:
        CALL RSEP3T('REC',TEXTP,' ')
        IF(TEXTP.EQ.' ') THEN
          NREC=0
        ELSE
          OPEN(LUREC,FILE=TEXTP,STATUS='OLD')
          READ(LUREC,*) (TEXTP,I=1,20)
          DO 11 I=1,MREC+1
            TEXTP='$'
            AUX(1)=0.
            AUX(2)=0.
            AUX(3)=0.
            AUX(4)=0.
            AUX(5)=0.
            READ(LUREC,*,END=12) TEXTP,(AUX(J),J=1,5)
            IF(TEXTP.EQ.'$') THEN
              GO TO 12
            END IF
            IF(I.GT.MREC) THEN
C             <A NAME="error641"></A><FONT COLOR="RED">641</FONT>
              CALL ERROR
     *                 ('641 in RPAR1: Too small array XREC in /RPARD/')
C             The number of given receivers or profiles is exceeding the
C             dimension MREC of the array XREC in the common block
C             /RPARD/.  MREC should thus be adjusted accordingly.
C             <A HREF="rpard.inc">/RPARD/</A>
            END IF
            IF(ISRFX(1).LT.0) THEN
              XREC(1,I)=AUX(-ISRFX(1))
            ELSE IF(ISRFX(1).GT.0) THEN
              XREC(1,I)=AUX(4)
            END IF
            IF(ISRFX(2).LT.0) THEN
              XREC(2,I)=AUX(-ISRFX(2))
            ELSE IF(ISRFX(2).GT.0) THEN
              XREC(2,I)=AUX(5)
            END IF
   11     CONTINUE
   12     CONTINUE
          NREC=I-1
          CLOSE(LUREC)
        END IF
C
      END IF
      LURPAR=LUN
      RETURN
      END
C
C=======================================================================
C
C     <A NAME="RPAR2"></A>
C
      SUBROUTINE RPAR2(IRAY,PAR1,PAR2)
      REAL PAR1,PAR2
      INTEGER IRAY
C
C This subroutine determines the take-off parameters of the ray.
C
C Input:
C     IRAY... Number of the already computed rays.  IRAY=0 when
C             beginning the computation of a new elementary wave.
C             Otherwise, the output from the previous invocation of
C             RPAR2.
C
C Output:
C     IRAY... IRAY=0 if all rays are computed and the computation of the
C             elementary wave will be terminated.
C             Otherwise, input value increased by 1.
C     PAR1,PAR2... Take-off parameters of the ray.
C
C Common block /RPARD/:
      INCLUDE 'rpard.inc'
C     <A HREF="rpard.inc">rpard.inc</A>
C Storage locations PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM,
C IPOINT,PRM1 of the common block may be redefined in this subroutine.
C
C Common blocks /RPARC/ and /RPARH/:
      INCLUDE 'rparc.inc'
C     <A HREF="rparc.inc">rparc.inc</A>
C Storage locations IRAY0,JRAY,JTYPE,G1,G2 of common block /RPARC/ are
C altered.
C Storage locations JPOINT and KMAH of common block /RPARH/ are
C annulled.
C Storage locations NRPARH and IRPARH(0) of common block /RPARH/ may be
C altered.
C
C Subroutines and external functions required:
      EXTERNAL RP2D
C     RP2D... One-parametric boundary-value ray tracing.  This file.
C
C Date: 2006, April 19
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER I,IAUX
      REAL AUX1,AUX2,AUX3,AUX4,AUX5,AUX6,AUX7,AUX8
C
C.......................................................................
C
      JPOINT=0
      KMAH=0
C
      IF(IRAY.EQ.0) THEN
C       New elementary wave
        IRAY0=0
        IF(ISRFX(2).EQ.0) THEN
          G2=0.
        END IF
        GO TO 90
      END IF
C
C     New ray - determination of shooting parameters G1,G2:
   10 CONTINUE
        JRAY=IRAY-IRAY0
        IF(ISRFX(2).EQ.0) THEN
C         One-parametric shooting (or initial-value ray tracing):
          CALL RP2D(IRAY0,JRAY,G1)
          IF(JRAY.EQ.0) THEN
C           End of shooting - new one-parametric slice or new data
            G2=G2+1.
            IF(G2.GT.BNUM+0.001) THEN
C             New data or new elementary wave
              GO TO 90
            ELSE
C             New one-parametric slice
              IRAY0=IRAY
              GO TO 10
            END IF
          END IF
        ELSE
C         Two-parametric shooting:
          CALL RP3D(JRAY,JTYPE,G1,G2)
          IF(JRAY.EQ.0) THEN
C           End of shooting - new data or new elementary wave
            GO TO 90
          END IF
        END IF
        IRAY=IRAY0+JRAY
C
C       New ray - determination of take-off parameters PAR1,PAR2:
        NRPARH=IRPARH(IRPARH(0))
        IF(ISRFX(2).EQ.0) THEN
C         One-parametric shooting:
C         Projection of the shooting domain (0,ANUM)*(0,BNUM)
          PAR1=PAR1L
          PAR2=PAR2L
          IF(ANUM.NE.0.) THEN
            PAR1=PAR1+(PAR1A-PAR1L)*G1/ANUM
            PAR2=PAR2+(PAR2A-PAR2L)*G1/ANUM
          END IF
          IF(BNUM.NE.0.) THEN
            PAR1=PAR1+(PAR1B-PAR1L)*G2/BNUM
            PAR2=PAR2+(PAR2B-PAR2L)*G2/BNUM
          END IF
        ELSE
C         Two-parametric shooting:
C         Projection of the shooting domain (0,1)*(0,1)
          PAR1=PAR1L+(PAR1A-PAR1L)*G1+(PAR1B-PAR1L)*G2
          PAR2=PAR2L+(PAR2A-PAR2L)*G1+(PAR2B-PAR2L)*G2
        END IF
        RETURN
C
C       Reading input data
   90   CONTINUE
        AUX1=PAR1L
        AUX2=PAR2L
        AUX3=PAR1A
        AUX4=PAR2A
        AUX5=PAR1B
        AUX6=PAR2B
        AUX7=ANUM
        AUX8=BNUM
        IAUX=IPOINT
        DO 91 I=1,MPRM1
          PRM1(I)=0.
   91   CONTINUE
        READ(LURPAR,*)
     *         PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM,IPOINT,PRM1
        IF(IRAY.NE.0) THEN
          IF(PAR1L.EQ.AUX1.AND.PAR2L.EQ.AUX2.AND.PAR1A.EQ.AUX3.AND.
     *       PAR2A.EQ.AUX4.AND.PAR1B.EQ.AUX5.AND.PAR2B.EQ.AUX6.AND.
     *       ANUM.EQ.AUX7.AND.BNUM.EQ.AUX8.AND.IAUX.EQ.IPOINT) THEN
C           New elementary wave
            IRAY=0
C           Writing the table of ray histories:
            WRITE(LUWARN(0),'(A)') 'Table of ray histories:'
            DO 92 I=1,IRPARH(0)
              WRITE(LUWARN(0),'(I4,A,I3,19I4,(4X,19I4))')
     *          INDMIN+I,':',(IRPARH(IAUX),IAUX=IRPARH(I-1)+1,IRPARH(I))
   92       CONTINUE
            RETURN
          END IF
        END IF
        IF(ISRFX(1).NE.0) THEN
          IF(PAR1L.EQ.PAR1A.AND.PAR2L.EQ.PAR2A) THEN
C           <A NAME="error645"></A><FONT COLOR="RED">645</FONT>
            CALL ERROR('645 in RPAR2: Zero ray-parameter domain')
C           (PAR1A,PAR2A) cannot equal (PAR1L,PAR2L) for one or two
C           parametric boundary-value ray tracing indicated by non-zero
C           ISRFX1.
          END IF
        END IF
        IF(ISRFX(2).NE.0) THEN
          IF(PAR1L.EQ.PAR1B.AND.PAR2L.EQ.PAR2B) THEN
C           <A NAME="error646"></A><FONT COLOR="RED">646</FONT>
            CALL ERROR('646 in RPAR2: Zero ray-parameter domain')
C           (PAR1B,PAR2B) cannot equal (PAR1L,PAR2L) for two-parametric
C           boundary-value ray tracing indicated by non-zero ISRFX2.
          END IF
          IF(ANUM.LE.0..OR.BNUM.LE.0.) THEN
C           <A NAME="error653"></A><FONT COLOR="RED">653</FONT>
            CALL ERROR('653 in RPAR2: Non-positive ANUM or BNUM')
C           Both ANUM and BNUM must be positive for two-parametric
C           boundary-value ray tracing indicated by non-zero ISRFX2.
          END IF
        END IF
        IF(ISRFX(1).EQ.0.AND.ISRFR.NE.0.AND.NREC.NE.0.AND.IPOINT.NE.0)
     *                                                              THEN
C         <A NAME="error647"></A><FONT COLOR="RED">647</FONT>
          CALL ERROR('647 in RPAR2: Zero value of ISRFX1')
C         ISRFX1 must not be zero if the boundary-value ray tracing
C         is to be performed, i.e. if ISRFR.NE.0 and NREC.NE.0 and
C         IPOINT.NE.0 in the input data.
        END IF
C       Initializing ray histories
        INDMIN=0
        IRPARH(0)=0
C       New initial ray
        IRAY0=IRAY
        G2=0.
      GO TO 10
C
      END
C
C=======================================================================
C
C     <A NAME="RPAR30"></A>
C
      SUBROUTINE RPAR30(TTINI,TT,Q11,Q21,Q12,Q22)
      REAL TTINI,TT,Q11,Q21,Q12,Q22
C
C This subroutine performs the numerical quadrature of the squared
C matrix of geometrical spreading, useful for the estimination of the
C L2 norm of the distance between neighbouring rays in order to control
C the fatness of ray tubes.  This subroutine is designed to be called
C by subroutine PSHIFT of file 'raycb.for'.  Subroutine PSHIFT is called
C both with steps STEP and STORE along the ray.
C
C Input:
C     TTINI...Travel time at the initial point of the ray, used to
C             identify the inititial point of a ray.
C     TT...   Travel time (independent variable for the numerical
C             quadrature along the ray).
C     Q11,Q21,Q12,Q22... Matrix of geometrical spreading.
C None of the input parameters are altered.
C
C No output.
C
C Common block /RPARD/:
      INCLUDE 'rpard.inc'
C     <A HREF="rpard.inc">rpard.inc</A>
C None of the storage locations of the common block are altered.
C
C Common block /RPARC/:
      INCLUDE 'rparc.inc'
C     <A HREF="rparc.inc">rparc.inc</A>
C Storage locations of the common block may be altered.
C
C No subroutines and external functions required.
C
C Date: 1999, May 25
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
      REAL TT1,GG111,GG112,GG122,AUX
      REAL TT2,GG211,GG212,GG222
      SAVE TT2,GG211,GG212,GG222
C
      IF(TT.EQ.TTINI) THEN
        GG11=0.
        GG12=0.
        GG22=0.
      ELSE
        TT1  =TT2
        GG111=GG211
        GG112=GG212
        GG122=GG222
      END IF
      TT2  =TT
      GG211=Q11*Q11+Q21*Q21
      GG212=Q11*Q12+Q21*Q22
      GG222=Q12*Q12+Q22*Q22
      IF(TT.NE.TTINI) THEN
        IF(PRM0(4).LT.0) THEN
C         Metric at the last point
          GG11=GG211
          GG12=GG212
          GG22=GG222
        ELSE
C         Average squared metric
          AUX=(TT2-TT1)/2.
          GG11=GG11+(GG111+GG211)*AUX
          GG12=GG12+(GG112+GG212)*AUX
          GG22=GG22+(GG122+GG222)*AUX
        END IF
      END IF
      RETURN
      END
C
C=======================================================================
C
C     <A NAME="RPAR31"></A>
C
      SUBROUTINE RPAR31(YL,Y,YY,IY)
      REAL YL(6),Y(35),YY(5)
      INTEGER IY(12)
C
C This subroutine may be used to store the quantities useful for the
C determination of the take-off parameters in the memory.  It is called
C with constant step STORE (see the input data in the subroutine file
C 'ray.for') of the independent variable along the ray (if STORE.NE.0),
C and at the points of intersection with interfaces either before and
C after the transformation (in any case).
C
C Input:
C     YL...   Array containing local quantities at the point of the ray.
C     Y...    Array containing basic quantities computed along the ray.
C     YY...   Array containing real auxiliary quantities computed along
C             the ray.
C     IY...   Array containing integer auxiliary quantities computed
C             along the ray.
C None of the input parameters are altered.
C
C No output.
C
C Common block /RPARD/:
      INCLUDE 'rpard.inc'
C     <A HREF="rpard.inc">rpard.inc</A>
C None of the storage locations of the common block are altered.
C
C Common block /RPARH/:
      INCLUDE 'rparc.inc'
C     <A HREF="rparc.inc">rparc.inc</A>
C Storage locations of the common block may be altered.
C
C Subroutines and external functions required:
      EXTERNAL NSRFC
      INTEGER NSRFC
C     NSRFC... File 'model.for' of the package 'MODEL'.
C
C Date: 1997, October 15
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     No auxiliary storage locations.
C
      IF(IPOINT.NE.0) THEN
        IF(NRPARH.EQ.IRPARH(IRPARH(0))) THEN
C         Initial element of the ray
          NRPARH=NRPARH+1
          CALL SRPARH
          IRPARH(NRPARH)=IY(5)
        END IF
        IF(JPOINT.LT.IABS(IPOINT)) THEN
          IF(IY(6).NE.0.AND.(IABS(IY(6)).LE.NSRFC().OR.
     *                       PRM0(3).EQ.0.0.OR.
     *                       PRM0(3).GE.1.5            )) THEN
C           Structural interface (or model boundary)
            NRPARH=NRPARH+1
            CALL SRPARH
            IF(MOD(NRPARH-IRPARH(IRPARH(0)),2).EQ.0) THEN
C             Incidence at a structural interface
              IRPARH(NRPARH)=IY(6)
            ELSE
C             Leaving a structural interface
              IRPARH(NRPARH)=IY(5)
            END IF
          END IF
        END IF
      END IF
      RETURN
      END
C
C=======================================================================
C
C     <A NAME="RPAR32"></A>
C
      SUBROUTINE RPAR32(ISTOR,YL,Y,YY,IY)
      INTEGER ISTOR,IY(12)
      REAL YL(6),Y(27),YY(5)
C
C This subroutine may be used to store the quantities useful for the
C determination of the take-off parameters in the memory.  It is called
C at the points of intersection of the ray with interfaces specified in
C the input data in the subroutine file 'ray.for' for storing the
C computed quantities, see also
C <A HREF="../../../papers/b88vc200.htm#29">C.R.T.5.5.2</A>.
C
C Input:
C     ISTOR...The sequential number in the input data of the specified
C             surface.
C     YL...   Array containing local quantities at the point of the ray.
C     Y...    Array containing basic quantities computed along the ray.
C     YY...   Array containing real auxiliary quantities computed along
C             the ray.
C     IY...   Array containing integer auxiliary quantities computed
C             along the ray.
C None of the input parameters are altered.
C
C No output.
C
C Common block /DCRT/ (subroutine file 'ray.for'):
      INCLUDE 'dcrt.inc'
C     <A HREF="dcrt.inc">dcrt.inc</A>
C None of the storage locations of the common block are altered.
C
C Common block /INITC/ (see subroutine file 'init.for'):
      INCLUDE 'initc.inc'
C     <A HREF="initc.inc">initc.inc</A>
C None of the storage locations of the common block are altered.
C
C Common block /RPARD/:
      INCLUDE 'rpard.inc'
C     <A HREF="rpard.inc">rpard.inc</A>
C None of the storage locations of the common block are altered.
C
C Common block /RPARC/ and /RPARH/:
      INCLUDE 'rparc.inc'
C     <A HREF="rparc.inc">rparc.inc</A>
C Storage locations X1,X2,X1G1,X2G1,X1G2,X2G2 of common block /RPARC/
C may be altered.
C Storage locations JPOINT and KMAH of common block /RPARH/ may be
C altered.
C
C Subroutines and external functions required:
      EXTERNAL METRIC,SMVPRD,XFUN
C     METRIC... File 'metric.for' of the package 'MODEL'.
C     SMVPRD... File 'means.for' of the package 'MODEL'.
C     XFUN... This file.
C
C Date: 1997, July 6
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Substitution in the common block /RPARC/:
      REAL X(2),XG(4)
      EQUIVALENCE (X(1),X1),(XG(1),X1G1)
C
C     Auxiliary storage locations:
      INTEGER NX,IX
      REAL GSQRD,G(12),AUX(18),AUX1,AUX2,AUX3
      EQUIVALENCE (AUX(11),AUX1),(AUX(12),AUX2),(AUX(13),AUX3)
      REAL H42,H52,H62,RH1,RH2,RR,XH1,XH2,XR,SIDE1,SIDE2
C
C     NX...   Number of defined Xi functions.
C     IX...   Auxiliary loop variable
C     GSQRD,G... See subroutine METRIC of the file 'metric.for'.
C     AUX,AUX1,AUX2,AUX3... Auxiliary storage locations:
C             AUX(1:18)... Christoffel symbols,
C             AUX(1:4)... Value and gradient of the function specifying
C               the reference surface ISRFR,
C             AUX(5:10)... Second derivatives of the function specifying
C               the reference surface ISRFR,
C             AUX(5:7)... Gradient of the Xi function.
C             AUX(11:13)... AUX1,AUX2,AUX3,
C             AUX1,AUX2,AUX3... Contravariant derivatives the functions
C               etc.
C     H42,H52,H62... Contravariant components of the 2-nd basis vector
C             of R.C.C.S.
C     RH1,RH2... Derivatives of the function describing reference the
C             surface ISRFR in R.C.C.S.
C     RR...   Norm of the gradient of the function describing reference
C             the surface ISRFR in R.C.C.S.
C     XH1,XH2... Derivatives of the Xi functions in R.C.C.S.
C     XR...   Projection of the gradient of the Xi functions onto the
C             normal to the reference surface.
C     SIDE1,SIDE2... Vectorial side of the shooting domain, in the ray
C             parameters PAR1,PAR2.
C
C.......................................................................
C
      IF(KSTOR(ISTOR).EQ.ISRFR) THEN
        DO 1 IX=1,ISTOR-1
          IF(KSTOR(IX).EQ.ISRFR) THEN
C           The surface is already encountered
            GO TO 9
          END IF
    1   CONTINUE
        JPOINT=JPOINT+1
        IF(JPOINT.EQ.IABS(IPOINT)) THEN
C         Successful ray:
C
C         KMAH index:
          KMAH=IY(12)
C
C         Metric tenzor
          CALL METRIC(Y(3),GSQRD,G,AUX)
C
C         Contravariant components of the 2-nd basis vector of R.C.C.S.:
          CALL SMVPRD(G(7),Y(6),Y(7),Y(8),AUX1,AUX2,AUX3)
          AUX1=GSQRD*SQRT(AUX1*Y(6)+AUX2*Y(7)+AUX3*Y(8))
          H42=(Y(7)*Y(11)-Y(8)*Y(10))/AUX1
          H52=(Y(8)*Y( 9)-Y(6)*Y(11))/AUX1
          H62=(Y(6)*Y(10)-Y(7)*Y( 9))/AUX1
C
C         Gradient of the function describing reference surface ISRFR:
C         Covariant components
          IF(IABS(ISRFR).LE.100) THEN
            CALL SRFC2(IABS(ISRFR),Y(3),AUX)
          ELSE
C           Reference surface limits the computational volume
            AUX(2)=0.
            AUX(3)=0.
            AUX(4)=0.
            AUX((IABS(ISRFR)-97)/2)=1.
          END IF
C         Contravariant components
          CALL SMVPRD(G(7),AUX(2),AUX(3),AUX(4),AUX1,AUX2,AUX3)
C         Gradient in R.C.C.S.
          RH1=AUX1*Y(9)+AUX2*Y(10)+AUX3*Y(11)
          RH2=AUX(2)*H42+AUX(3)*H52+AUX(4)*H62
C         Projection onto the normal to the reference surface
          RR =AUX1*AUX(2)+AUX2*AUX(3)+AUX3*AUX(4)
C
C         Number of Xi functions:
          IF(ISRFX(2).EQ.0) THEN
            IF(ISRFX(1).EQ.0) THEN
              NX=0
            ELSE
              NX=1
            END IF
          ELSE
            NX=2
          END IF
C
C         Loop for Xi functions:
          DO 5 IX=1,NX
C
C           Value and gradient of the Xi functions:
C           Covariant components
            CALL XFUN(IX,Y(3),X(IX),AUX(5))
C           Contravariant components
            CALL SMVPRD(G(7),AUX(5),AUX(6),AUX(7),AUX1,AUX2,AUX3)
C           Gradient in R.C.C.S.
            XH1=AUX1*Y(9)+AUX2*Y(10)+AUX3*Y(11)
            XH2=AUX(5)*H42+AUX(6)*H52+AUX(7)*H62
C           Projection onto the normal to the reference surface
            XR =AUX1*AUX(2)+AUX2*AUX(3)+AUX3*AUX(4)
C
C           Derivatives of the Xi function along the reference surface
C           in R.C.C.S.:
            AUX1=XH1*RH2-XH2*RH1
            AUX2=1.-(RH1*RH1+RH2*RH2)/RR
            IF(AUX2.LE.0.) THEN
C             AUX2 is the square of cosine of the incidence angle at the
C             reference surface.
C             <A NAME="error648"></A><FONT COLOR="RED">648</FONT>
              CALL ERROR
     *               ('648 in RPAR32: Ray tangent to reference surface')
C             If the source is situated at the reference surface, the
C             rays tangent to the reference surface should be avoided
C             by adjusting the ray-parameter domain properly.
C             The angular deviation of rays from the reference surface
C             should be at least 0.001 radians.
C             <A HREF="#data4">Data for the ray-parameter domain</A>
C             If this is not the case, please, contact the author.
            END IF
            XH1=(XH1-(XR*RH1+AUX1*RH2)/RR)/AUX2
            XH2=(XH2-(XR*RH2-AUX1*RH1)/RR)/AUX2
C
C           Derivative of the Xi function with respect to the shooting
C           parameters:
            SIDE1=PAR1A-PAR1L
            SIDE2=PAR2A-PAR2L
            XG(IX)=  ((XH1*Y(12)+XH2*Y(13))*(YI(12)*SIDE1+YI(16)*SIDE2)
     *               +(XH1*Y(16)+XH2*Y(17))*(YI(13)*SIDE1+YI(17)*SIDE2)
     *               +(XH1*Y(20)+XH2*Y(21))*(YI(14)*SIDE1+YI(18)*SIDE2)
     *               +(XH1*Y(24)+XH2*Y(25))*(YI(15)*SIDE1+YI(19)*SIDE2))
            SIDE1=PAR1B-PAR1L
            SIDE2=PAR2B-PAR2L
            XG(IX+2)=((XH1*Y(12)+XH2*Y(13))*(YI(12)*SIDE1+YI(16)*SIDE2)
     *               +(XH1*Y(16)+XH2*Y(17))*(YI(13)*SIDE1+YI(17)*SIDE2)
     *               +(XH1*Y(20)+XH2*Y(21))*(YI(14)*SIDE1+YI(18)*SIDE2)
     *               +(XH1*Y(24)+XH2*Y(25))*(YI(15)*SIDE1+YI(19)*SIDE2))
C
    5     CONTINUE
C
C         Undefined functional values:
          DO 6 IX=NX+1,2
            X(IX)=0.
            XG(IX)=0.
            XG(IX+2)=0.
    6     CONTINUE
C
        END IF
      END IF
    9 CONTINUE
      RETURN
      END
C
C=======================================================================
C
C     <A NAME="RPAR33"></A>
C
      SUBROUTINE RPAR33(YL,Y,YY,IY)
      REAL YL(6),Y(35),YY(5)
      INTEGER IY(12)
C
C This subroutine may be used to store the quantities useful for the
C determination of the take-off parameters in the memory.  It is called
C at the endpoints of the elements of rays, situated at structural
C interfaces, see also
C <A HREF="../../../papers/b88vc200.htm#29">C.R.T.5.5.3</A>.
C
C Input:
C     YL...   Array containing local quantities at the point of the ray.
C     Y...    Array containing basic quantities computed along the ray.
C     YY...   Array containing real auxiliary quantities computed along
C             the ray.
C     IY...   Array containing integer auxiliary quantities computed
C             along the ray.  For the definition of arrays YL, Y, YY and
C             IY see the file 'ray.for'.
C None of the input parameters are altered.
C
C No output.
C
C No subroutines and external functions required.
C
C Date: 1990, January 22
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
      RETURN
      END
C
C=======================================================================
C
C     <A NAME="RPAR4"></A>
C
      SUBROUTINE RPAR4(IRAY,PAR1,PAR2,YL,Y,YY,IY,IEND,ISHEET,IREC)
C
      REAL PAR1,PAR2,YL(6),Y(35),YY(5)
      INTEGER IRAY,IY(12),IEND,ISHEET,IREC
C
C This subroutine is called after finishing the computation of a ray.
C It defines the storage locations YI(22) to YI(25) of the common block
C /INITC/.
C
C Input:
C     IRAY... Index of the ray, i.e. the number of the already computed
C             rays, i.e. the output from the last invocation of RPAR2.
C     PAR1,PAR2... Take-off parameters of the ray.
C     YL...   Array containing local quantities at the point of the ray.
C     Y...    Array containing basic quantities computed along the ray.
C     YY...   Array containing real auxiliary quantities computed along
C             the ray.
C     IY...   Array containing integer auxiliary quantities computed
C             along the ray.  For the definition of arrays YL, Y, YY and
C             IY see the file 'ray.for'.
C     IEND... Reason of the termination of the computation of the ray,
C             see
C             <A HREF="../../../papers/b88vc200.htm#27">C.R.T.5.4</A>.
C             See subroutine RAY in the subroutine file 'ray.for' for
C             <A HREF="ray.for#IEND">detailed description of IEND</A>.
C
C Output:
C     ISHEET..Ray-history index.  The different ray histories are
c             consecutively indexed by positive integers 1,2,3,...
C             according to their appearance during ray tracing.
C             The ray histories are indexed independently within each
C             elementary wave.
C             The ray-history indices are complemented with sign:
C             Positive - successful ray (crossing reference surface),
C             negative - unsuccessful ray (terminating before crossing
C             reference surface).
C                                                    <A NAME="IREC"></A>
C     IREC... Index of the receiver for a two-point ray,
C             0 for a basic ray (and other rays in 2-D),
C             -1 for other rays.
C
C Common block /INITC/ (see subroutine file 'init.for'):
      INCLUDE 'initc.inc'
C     <A HREF="initc.inc">initc.inc</A>
C     ------------------------------------------------------------------
C     YI(22)..Area of the element of the ray-parameter surface,
C             corresponding to the ray,
C     YI(23),YI(24),YI(25)... Components 11, 12, 22 of the symmetric
C             matrix inverse to the specific moment of the element of
C             the ray-parameter surface corresponding to the ray, see
C             <A HREF="../../../papers/b88vc200.htm#52">C.R.T.6.1</A>.
C For one-parametric shooting, locations YI(22) to YI(25) of the common
C block are defined in this subroutine.
C For two-parametric shooting, locations YI(22) to YI(25) of the common
C block are set to zeros in this subroutine.
C Other storage locations remain unchanged.
C
C Common block /RPARD/:
      INCLUDE 'rpard.inc'
C     <A HREF="rpard.inc">rpard.inc</A>
C None of the storage locations of the common block are altered.
C
C Common block /RPARC/:
      INCLUDE 'rparc.inc'
C     <A HREF="rparc.inc">rparc.inc</A>
C Storage locations of the common block are not altered.
C
C Subroutines and external functions required:
      EXTERNAL IHIST,KOOR,METRIC,SMVPRD,RPMEM
      INTEGER  IHIST,KOOR
C     IHIST...This file.
C     KOOR,METRIC... File 'metric.for' of the package 'MODEL'.
C     SMVPRD..File 'means.for' of the package 'MODEL'.
C     RPMEM...File 'rp3d.for' of the package 'CRT'.
C
C Date: 2006, February 23
C Coded by Ludek Klimes and Petr Bulant
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      REAL GSQRD,G(12),AUX(18)
      REAL X1A
      REAL RG11,RG12,RG22,SG11,SG12,SG22,SIDE1,SIDE2,AUX1,AUX2,AUX3,AUX4
      REAL XG11,XG12,XG22,G1X1,G2X1,G1X2,G2X2
C
C     GSQRD,G,AUX... For subroutine METRIC.
C     X1A...  Derivative of the X1 function with respect to the
C             shooting parameter A.
C     RG11,RG12,RG22... Ray take-off parameter metric tensor.
C     SG11,SG12,SG22... Shooting parameter metric tensor.
C     SIDE1,SIDE2,AUX1,AUX2,AUX3... Temporary variables.
C     G1X1,G1X2,G2X1,G2X2... Derivatives of shooting parameters with
C             respect to the Xi functions.
C
C.......................................................................
C
C     Ray-history index: positive - successful ray,
C     otherwise - unsuccessful ray
      IF(JPOINT.LT.IABS(IPOINT)) THEN
C     KMAH index of an unsuccessful ray:
        KMAH=IY(12)
      END IF
      ISHEET=IHIST(IY,IEND)
C
C     Calculating shooting-parameter metric tensor, storing the results:
      IF(ISRFX(2).EQ.0) THEN
C       Storing the quantities for one-parametric shooting:
        IF(ANUM.GT.0.) THEN
          X1A=X1G1/ANUM
        ELSE
          X1A=0.
        END IF
        CALL RP2DM(JRAY,ISHEET,X1,X1A,IREC)
C       X1,X1G1... Determined in RPAR32 (if defined),
C       ISHEET... Determined in this subroutine.
C       IREC... Determined in RP2DM.
      ELSE
C       (a) metric tensor per ray take-off parameters:
C       Check for non-existing rays:
        IF(IEND.LT.70) THEN
C         Cartesian metric tensor (zero for a point source)
          RG11=      YI(12)*YI(12)+YI(13)*YI(13)
          RG12=      YI(12)*YI(16)+YI(13)*YI(17)
          RG22=      YI(16)*YI(16)+YI(17)*YI(17)
          IF(ABS(YI(12)*YI(17)-YI(13)*YI(14)).LT.0.000001) THEN
C           Plus unit-sphere metric tensor for point or line source
            IF(KOOR().EQ.0) THEN
              AUX1=YI(6)**2+YI(7)**2+YI(8)**2
            ELSE
              CALL METRIC(YI(3),GSQRD,G,AUX)
              CALL SMVPRD(G(7),YI(6),YI(7),YI(8),AUX1,AUX2,AUX3)
              AUX1=AUX1*YI(6)+AUX2*YI(7)+AUX3*YI(8)
            END IF
            RG11=RG11+(YI(14)*YI(14)+YI(15)*YI(15))/AUX1
            RG12=RG12+(YI(14)*YI(18)+YI(15)*YI(19))/AUX1
            RG22=RG22+(YI(18)*YI(18)+YI(19)*YI(19))/AUX1
          END IF
C         The above two metric tensors (equal along an exploding unit
C         sphere source) should be combined in a more sophisticated way.
C         **************************************************************
C         The metric tensor should be determined in 'init.for' !
C         **************************************************************
        ELSE
          RG11=1.
          RG12=0.
          RG22=1.
        END IF
        IF(IEND.GE.70) THEN
C         The ray has not been traced, GGij have not been initialized.
          GG11=0.
          GG12=0.
          GG22=0.
        END IF
        IF(PRM0(4).GE.0.) THEN
          IF(Y(1)-YI(1).GT.0.) THEN
            GG11=GG11/(Y(1)-YI(1))
            GG12=GG12/(Y(1)-YI(1))
            GG22=GG22/(Y(1)-YI(1))
          ELSE
            GG11=0.
            GG12=0.
            GG22=0.
          END IF
        END IF
C       (b) Metric tensor per shooting parameters:
C       SCALE renormalizes the unit-sphere metric in radians to the
C       metric in some scaled units, proportional to radians by the
C       factors of scale in the directions of both the shooting
C       parameters.
        SIDE1=PAR1A-PAR1L
        SIDE2=PAR2A-PAR2L
        AUX3=GG11*SIDE1+GG12*SIDE2
        AUX4=GG12*SIDE1+GG22*SIDE2
        XG11=SIDE1*AUX3+SIDE2*AUX4
        SCALE=ANUM/SQRT(SIDE1*SIDE1+SIDE2*SIDE2)
        SIDE1=SIDE1*SCALE
        SIDE2=SIDE2*SCALE
        AUX1=RG11*SIDE1+RG12*SIDE2
        AUX2=RG12*SIDE1+RG22*SIDE2
        SG11=SIDE1*AUX1+SIDE2*AUX2
        SIDE1=PAR1B-PAR1L
        SIDE2=PAR2B-PAR2L
        XG12=SIDE1*AUX3+SIDE2*AUX4
        AUX3=GG11*SIDE1+GG12*SIDE2
        AUX4=GG12*SIDE1+GG22*SIDE2
        XG22=SIDE1*AUX3+SIDE2*AUX4
        SCALE=BNUM/SQRT(SIDE1*SIDE1+SIDE2*SIDE2)
        SIDE1=SIDE1*SCALE
        SIDE2=SIDE2*SCALE
        SG12=SIDE1*AUX1+SIDE2*AUX2
        AUX1=RG11*SIDE1+RG12*SIDE2
        AUX2=RG12*SIDE1+RG22*SIDE2
        SG22=SIDE1*AUX1+SIDE2*AUX2
C       Another extremely simple metric, conformal with the metric used
C       at one-parametric shooting, and coinciding with the above scaled
C       metric for two-parametric shooting in the case of a small
C       rectangular shooting domain situated at the equator:
C       SG11=ANUM*ANUM
C       SG12=0.
C       SG22=BNUM*BNUM
C       Regularization near the poles:
        SG11=SG11+ANUM*ANUM/100.
        SG22=SG22+BNUM*BNUM/100.
C       (c) Calculating the derivatives for two-parametric shooting:
        IF(ISHEET.GT.0) THEN
          AUX1=X1G1*X2G2-X2G1*X1G2
          IF(AUX1.NE.0.) THEN
            G1X1= X2G2/AUX1
            G2X1=-X2G1/AUX1
            G1X2=-X1G2/AUX1
            G2X2= X1G1/AUX1
          ELSE
C           Caustic:
            AUX1=X1G1*X1G1+X2G1*X2G1+X1G2*X1G2+X2G2*X2G2
            IF(AUX1.NE.0.) THEN
C             Line caustic - generalized inverse matrix:
              G1X1=X1G1/AUX1
              G2X1=X1G2/AUX1
              G1X2=X2G1/AUX1
              G2X2=X2G2/AUX1
            ELSE
C             Point caustic:
              G1X1=0.
              G2X1=0.
              G1X2=0.
              G2X2=0.
            END IF
          END IF
        ELSE
          X1=0.
          X2=0.
          G1X1=0.
          G2X1=0.
          G1X2=0.
          G2X2=0.
        END IF
C       (c) Storing the quantities for two-parametric shooting:
        CALL RPMEM(JRAY,JTYPE,ISHEET,G1,G2,SG11,SG12,SG22,
     *                         XG11,XG12,XG22,X1,X2,G1X1,G2X1,G1X2,G2X2)
C       JRAY,JTYPE,G1,G2... Determined in RPAR2 (output of RP3D),
C       X1,X2,X1G1,X2G1,X1G2,X2G2... Determined in RPAR32 (if defined),
C       ISHEET,SG11,SG12,SG22,G1X1,G2X1,G1X2,G2X2... Determined in this
C       subroutine.
C       (d) Identifying basic and two-point rays:
        IF (JTYPE.LT.-1000) THEN
C         Two-point ray:
          IREC=-JTYPE-1000
        ELSEIF (JTYPE.EQ.0) THEN
C         Basic ray:
          IREC=0
        ELSE
C         Auxiliary or boundary ray:
          IREC=-1
        END IF
      END IF
C
C     Determination of the area YI(22) and inverse specific moment
C     YI(23:25) of element of the ray-parameter surface, corresponding
C     to the ray:
      IF(ISRFX(2).EQ.0) THEN
        CALL RP2DG(YI(22),YI(23),YI(24),YI(25))
      ELSE
        YI(22)=0.
        YI(23)=0.
        YI(24)=0.
        YI(25)=0.
C       CALL RP3DG(YI(22),YI(23),YI(24),YI(25)) *** for future extension
      END IF
C
C     Storing additional quantities related to the shooting algorithm:
      YI(26)=G1
      YI(27)=G2
      YI(28)=X1
      YI(29)=X2
C
      RETURN
      END
C
C=======================================================================
C
C     <A NAME="XFUN"></A>
C
      SUBROUTINE XFUN(IFUN,COOR,X,XDER)
      INTEGER IFUN
      REAL COOR(3),X,XDER(3)
C
C Interface subroutine designed to return the value and first
C derivatives of the Xi functions describing the profiles.  It is called
C by the subroutine RPAR32, and may be user-defined.
C
C Input:
C     IFUN... 1 for the first (X1 function), and 2 for the second
C             receiver parameter (X2 function).
C     COOR... Array containing coordinates X1, X2, X3 of the given
C             point.
C None of the input parameters are altered.
C
C Output:
C     X...    Value of the corresponding Xi function at the given point.
C     XDER... Array containing the derivatives of the Xi function with
C             respect to the general coordinates X1, X2, X3.
C
C Common block /RPARD/:
      INCLUDE 'rpard.inc'
C     <A HREF="rpard.inc">rpard.inc</A>
C None of the storage locations of the common block are altered.
C
C Subroutines and external functions required:
      EXTERNAL SRFC2
C     SRFC2... File 'srfc.for' of the package 'MODEL'.
C
C Date: 1993, December 20
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      REAL FAUX(10)
C
      IF(ISRFX(IFUN).GT.0) THEN
        CALL SRFC2(ISRFX(IFUN),COOR,FAUX)
        X=FAUX(1)
        XDER(1)=FAUX(2)
        XDER(2)=FAUX(3)
        XDER(3)=FAUX(4)
      ELSE
        X=0.
        XDER(1)=0.
        XDER(2)=0.
        XDER(3)=0.
        IF(ISRFX(IFUN).LT.0) THEN
          X=COOR(-ISRFX(IFUN))
          XDER(-ISRFX(IFUN))=1.
        END IF
      END IF
      RETURN
      END
C
C=======================================================================
C
C     <A NAME="IHIST"></A>
C
      INTEGER FUNCTION IHIST(IY,IEND)
      INTEGER IY(*),IEND
C
C Auxiliary integer function to the subroutine RPAR4 designed to
C determine the ray-history index.  The ray-history index is used to
C distinguish between the rays with different histories for the purposes
C of the subroutine RPAR2.
C
C The ray history is recorded as a sequence of integers.  In this
C version, the history consists of:
C (A) Indices of complex blocks along the ray (including the sign),
C (B) Indices of structural interfaces intersected by the ray (including
C     the sign), and index of the end surface limiting the computational
C     volume, if intersected.
C (C) Reason IEND (see subroutine RAY2 of 'ray.for' and INIT2 of
C     'init.for') if the ray is not successful (does not reach the
C     reference surface).
C (D) KMAH index.
C Note - differences in versions 4.10 and older:
C (B) Index of the end surface not considered in the history.
C (C) IEND not considered if IEND.GE.40.
C Note - differences in versions 5.00 and older:
C (D) KMAH index is not included in the history of an unsuccessful ray.
C     It saves CPU time but may cause missing two-point rays, especially
C     close to the boundary of the computational volume.  This option
C     have been discarded in ver. 5.10 because of introducing controlled
C     initial-value ray tracing followed by interpolation within ray
C     cells.
C
C Input:
C     IY...   Array containing integer auxiliary quantities computed
C             along the ray.  For the definition of array IY see file
C             'ray.for'.
C     IEND... Reason of the termination of the computation of the ray,
C             see
C             <A HREF="../../../papers/b88vc200.htm#27">C.R.T.5.4</A>.
C             See subroutine RAY in the subroutine file 'ray.for' for
C             <A HREF="ray.for#IEND">detailed description of IEND</A>.
C None of the input parameters are altered.
C
C Output:
C     IHIST...Positive (1,2,3,...) for successful ray,
C             otherwise (0 or -1,-2,-3,...) for unsuccessful ray.
C             Successful ray is understood here to be the ray having at
C             least IABS(IPOINT) points of intersection with the
C             specified surface ISRFR, see the input data (2).
C             Different values of IHIST correspond to rays with
C             different histories.
C
C Common block /RPARD/:
      INCLUDE 'rpard.inc'
C     <A HREF="rpard.inc">rpard.inc</A>
C None of the storage locations of the common block are altered.
C
C Common block /RPARH/:
      INCLUDE 'rparc.inc'
C     <A HREF="rparc.inc">rparc.inc</A>
C Storage locations of the common block may be altered.
C
C Subroutines and external functions required:
      EXTERNAL SRPARH,CODE
C     SRPARH... This file.
C     CODE...   File 'code.for' of package 'CRT'.
C
C Date: 2000, May 14
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER I,K
C
      IF(IPOINT.EQ.0) THEN
        IHIST=0
      ELSE
        IF(JPOINT.LT.IABS(IPOINT)) THEN
C         Unsuccessful ray:
          IF(PRM0(3).GE.1.5) THEN
C           Removing the index of the surface, at which the ray
C           terminated, from the ray history:
            IF(IEND.LT.70) THEN
              IF(NRPARH-IRPARH(IRPARH(0)).EQ.0) THEN
C               <A NAME="error651"></A><FONT COLOR="RED">651</FONT>
                CALL ERROR('651 in IHIST: Empty ray history')
C               The history of a ray, which has defined initial
C               conditions (i.e., IEND.LT.70), should contain at
C               least two items.
C               This error should not appear.  Contact the authors.
              END IF
              IF(MOD(NRPARH-IRPARH(IRPARH(0)),2).EQ.1) THEN
C               <A NAME="error652"></A><FONT COLOR="RED">652</FONT>
                CALL ERROR('652 in IHIST: Odd length of ray history')
C               At this moment, the history of the ray should contain
C               an even number of items.
C               This error should not appear.  Contact the authors.
              END IF
              NRPARH=NRPARH-1
            END IF
          END IF
          IF(PRM0(3).LT.2.5) THEN
C           Including reason IEND of the termination of the computation
C           of a ray into the ray history:
            NRPARH=NRPARH+1
            CALL SRPARH
            IRPARH(NRPARH)=IEND
          ELSE IF(PRM0(3).LT.3.5) THEN
C           Distinguishing overcritical incidence (IEND=25) from other
C           reasons of the termination of the computation of a ray:
            IF(IEND.EQ.25) THEN
              NRPARH=NRPARH+1
              CALL SRPARH
              IRPARH(NRPARH)=IEND
            END IF
          END IF
C         Following 3 lines include KMAH index also in the history of an
C         unsuccessful ray:
          IF(IEND.LT.70) THEN
            NRPARH=NRPARH+1
            CALL SRPARH
            IRPARH(NRPARH)=KMAH
          END IF
        ELSE
C         Successful ray:
          NRPARH=NRPARH+1
          CALL SRPARH
          IRPARH(NRPARH)=KMAH
        END IF
C
        DO 13 K=1,IRPARH(0)
          IF(IRPARH(K)-IRPARH(K-1).EQ.NRPARH-IRPARH(IRPARH(0))) THEN
            DO 11 I=IRPARH(K-1)+1,IRPARH(K)
              IF(IRPARH(I).NE.IRPARH(I-IRPARH(K)+NRPARH)) THEN
                GO TO 12
              END IF
   11       CONTINUE
C           The history of the last ray is equal to the K-th ray history
            IHIST=INDMIN+K
            GO TO 20
          END IF
   12     CONTINUE
   13   CONTINUE
C       New ray history:
        NRPARH=NRPARH+1
        CALL SRPARH
        DO 18 I=NRPARH,IRPARH(0)+2,-1
          IRPARH(I)=IRPARH(I-1)
   18   CONTINUE
        DO 19 I=1,IRPARH(0)
          IRPARH(I)=IRPARH(I)+1
   19   CONTINUE
        IRPARH(0)=IRPARH(0)+1
        IRPARH(IRPARH(0))=NRPARH
        IHIST=INDMIN+IRPARH(0)
CH      WRITE(LUWARN(0),'(A/I4,A,I3,19I4,(4X,19I4))')
CH   *              'New ray history:',
CH   *              IHIST,':',(IRPARH(I),I=IRPARH(IRPARH(0)-1)+1,NRPARH)
C
   20   CONTINUE
        IF(JPOINT.LT.IABS(IPOINT)) THEN
          IHIST=-IHIST
        END IF
        IF(IPOINT.LT.0) THEN
          IF(IEND.EQ.50) THEN
            CALL CODE(IY,K,I,IEND)
            IF(IEND.NE.10) THEN
              END=50
            END IF
          END IF
          IF(IEND.NE.10) THEN
            IHIST=-IABS(IHIST)
          END IF
        END IF
      END IF
      RETURN
      END
C
C=======================================================================
C
C     <A NAME="SRPARH"></A>
C
      SUBROUTINE SRPARH()
C
C Auxiliary subroutine to the subroutine RPAR31 designed to shift the
C storage locations of the common block /RPARH/.
C
C No input.
C
C No output.
C
C Common block /RPARH/:
      INCLUDE 'rparc.inc'
C     <A HREF="rparc.inc">rparc.inc</A>
C Storage locations of the common block may be altered.
C
C No subroutines and external functions required.
C
C Date: 2000, May 14
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER ISHIFT,I
C
C     ISHIFT..Shift in the array IRPARH.
C     I...    Auxiliary loop variable.
C
C.......................................................................
C
    1 CONTINUE
      IF(NRPARH.GT.MRPARH) THEN
        IF(IRPARH(0).LE.2) THEN
C         <A NAME="error649"></A><FONT COLOR="RED">649</FONT>
          CALL ERROR('649 in SRPARH: Too small array IRPARH in /RPARH/')
C         Array IRPARH of the common block /RPARH/ designed to store
C         different ray histories is too small to contain the last three
C         different ray histories.  The dimension MRPARH of array
C         IRPARH should be adjusted in all declarations of /RPARH/.
        END IF
        INDMIN=INDMIN+1
        IRPARH(0)=IRPARH(0)-1
        ISHIFT=IRPARH(1)-IRPARH(0)
        NRPARH=NRPARH-ISHIFT
        DO 2 I=1,IRPARH(0)
          IRPARH(I)=IRPARH(I+1)-ISHIFT
    2   CONTINUE
        DO 3 I=IRPARH(0)+1,NRPARH-1
          IRPARH(I)=IRPARH(I+ISHIFT)
    3   CONTINUE
C       <A NAME="error650"></A><FONT COLOR="RED">650</FONT>
        WRITE(LUWARN(0),'(A/I4,A,I3,19I4,(4X,19I4))')
     *                  'Warning 650 in SRPARH: Discarded ray history:',
     *                  INDMIN,':',(IRPARH(I),I=IRPARH(0)+1,IRPARH(1))
C         Array IRPARH of the common block /RPARH/ designed to store
C         different ray histories is too small for all ray histories.
C         Some ray histories have thus been discarded.  Dimension MRPARH
C         of array IRPARH may be adjusted to prevent problems due to
C         discarded ray histories.
        GO TO 1
      END IF
      RETURN
      END
C
C=======================================================================
C                                                                 </PRE>
