C 
C Subroutine file 'means.for' containing some utility programs helpful
C when dealing with the model.
C
C Date: 1999, May 17
C Coded by Ludek Klimes
C
C.......................................................................
C
C This file consists of the following subroutines:
C     CDE...  Subroutine designed to search for the point of
C             intersection of the given curve with the boundaries of the
C             complex block.  The curve is interpolated from the two
C             given points.  The subroutine performs the steps 5.8.3(c),
C             (d) and (e) of the algorithm.
C             CDE
C     CROSS...Subroutine designed to find the point of intersection of a
C             curve with a surface (see C.R.T.5.8.4b).
C             CROSS
C     HIVD2...Subroutine performing the Hermite interpolation of a
C             vector and its derivatives using functional values and
C             derivatives at 2 given points (see C.R.T.5.8.4a).
C             HIVD2
C     SMVPRD..Subroutine designed to evaluate symmetric matrix by vector
C             product.  It may be, e.g., called after the invocation of
C             the METRIC subroutine to transform the covariant vectors
C             to the contravariant ones and vice versa.
C             SMVPRD
C
C=======================================================================
C
C     
C
      SUBROUTINE CDE(NOUGHT,NEND,KEND,NBOUND,KBOUND,BOUND,
     *      KDIM1,KDIM2,NDIM,IY,ERR,X0,X1,Y1,D1,X2,Y2,D2,X,Y,D,XB,YB,DB)
      INTEGER NOUGHT,NEND,KEND(*),NBOUND,KBOUND(*)
      INTEGER IY(8),KDIM1,KDIM2,NDIM,MY
      PARAMETER (MY=35)
      REAL BOUND(*),ERR,X0,X1,Y1(NDIM),D1(NDIM),X2,Y2(NDIM),D2(NDIM)
      REAL X,Y(NDIM),D(NDIM),XB,YB(NDIM),DB(NDIM)
C
C This subroutine determines the point of intersection of the given
C curve element with the boundary of the complex block.  The curve is
C interpolated from the two given points.  The subroutine performs the
C steps 5.8.3(c), (d) and (e) of the algorithm.
C
C General meaning of some arguments used:
C     X1,X2,X,XB... Independent variable along the curve.
C     Y1,Y2,Y,YB,D1,D2,D,DB... Arrays of the dimension NDIM.
C     Y1(KDIM1:KDIM2),Y2(KDIM1:KDIM2),Y(KDIM1:KDIM2),YB(KDIM1:KDIM2)...
C             Coordinates.
C     D1(1:NDIM),D2(1:NDIM),D(1:NDIM),DB(1:NDIM)... Derivatives of
C             arrays Y1,Y2,Y,YB with respect to the independent
C             variable.
C
C Input:
C     NOUGHT..Zero if the curve is not situated along the structural
C               interface,
C             otherwise the index of the surface on which the curve is
C               situated.
C     NEND... Number of end surfaces limiting the computational volume.
C             Usually NEND=0.
C     KEND... Contains the indices of end surfaces.
C             Array of dimension NEND, not used if NEND=0.
C     NBOUND..Number of isosurfaces of computed quantities, limiting the
C             computational volume.  These isosurfaces will likely be
C             the boundaries X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX of the
C             computational volume.
C     KBOUND..Indices of the quantities corresponding to the
C             isosurfaces.  The computational volume is limited by the
C             inequalities
C               Y(IABS(KBOUND(I))).LE.BOUND(I) for KBOUND(I).LT.0,
C               Y(IABS(KBOUND(I))).GE.BOUND(I) for KBOUND(I).GT.0.
C             Array of dimension NBOUND, not used if NBOUND=0.
C     BOUND...Values of the isosurfaces, likely coinciding with
C             boundaries X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX of the
C             computational volume.  In such a case, KBOUND would take
C             the values KDIM1,-KDIM1,KDIM1+1,-KDIM1-1,KDIM2,-KDIM2.
C             Array of dimension NBOUND, not used if NBOUND=0.
C     KDIM1,KDIM2... KDIM1-th to KDIM2-th elements of arrays Y1,Y2,Y,YB
C             contain coordinates.  KDIM2 is assumed to be KDIM1+2.
C     NDIM... Dimension of the arrays Y1,Y2,Y,YB,D1,D2,D,DB.  It must
C             not exceed the parameter MY.
C     IY...   Integer array of the dimension at least 8.
C     IY(4)=ISB1... Index of the simple block containing the point X1.
C     IY(5)=ICB1... Index of the complex block containing the point X1.
C             It may (but need not) be supplemented by a sign '+' for P
C             wave and sign '-' for S wave.
C     ERR...  Maximum error in independent variable for the
C             determination of the point of intersection.
C     X0...   For X.LE.X0 just the intersection with the boundary of
C             the computational volume is checked, not with the
C             structural interfaces.
C             For X1.LE.X0, the initial point X1 is checked for location
C             outside the complex block, but very close to its boundary.
C             In such a case, X0 should be close to X1 (within an
C             interval of ERR) and may be located inside the complex
C             block.
C     X1,Y1,D1... Quantities at the initial point of the curve element.
C     X2,Y2,D2... Quantities at another point of the curve.  The curve
C             is interpolated using the values and their derivatives at
C             the points X1 and X2.
C     X,Y,D...Quantities at the endpoint of the curve element.
C     XB,YB,DB... Values ignored.
C
C Output:
C     NOUGHT,KDIM1,KDIM2,NDIM,ERR,X0,X1,Y1,D1,X2,Y2,D2... Unchanged
C             input.
C     IY(1),IY(2),IY(3),IY(5)... Unchanged input.
C     IY(4)=ISB1... Index of the simple block containing the point X.
C Output if the endpoint of the curve element is situated in the complex
C block IY(5):
C     IY(6),IY(7),IY(8)... Unchanged input.
C     X,Y,D...Unchanged input.
C     XB,YB,DB... Copy of X,Y,D.
C Output if the endpoint of the curve element is situated in another
C complex block or outside the computational volume:
C     IY(6)=ISRF... Index of the surface at which the point of
C             intersection with the boundary of the complex block is
C             situated, supplemented by a sign '+' or '-' for the point
C             situated at the positive or negative side of the surface,
C             respectively.
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 point of intersection.
C             ISB2=0 for a free space on the other side of ISRF.
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 point of intersection.
C             ICB2=0 for a free space on the other side of ISRF.
C     X,Y,D...Values corresponding to the point of intersection of the
C             curve with the boundary of the complex block or the
C             computational volume.  If possible, this point is situated
C             inside the given complex block close to its boundary, or
C             directly on its boundary.
C     XB,YB,DB... Values corresponding to another approximation of the
C             point of intersection.  This point is situated outside the
C             the given complex block, close to its boundary or directly
C             on its boundary.
C             Note that XB-X should be less than ERR.
C
C Subroutines and external functions required:
      EXTERNAL NSRFC,BLOCK,SRFC2,CROSS
      INTEGER NSRFC
C     NSRFC,BLOCK... File 'model.for'.
C     SRFC2 and subsequent routines... File 'srfc.for' and subsequent
C              files (like 'val.for' and 'fit.for').
C     CROSS,HIVD2... This file.
C
C Date: 1999, May 17
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations for local model parameters:  FAUX(10),
C       G(12),GAMMA(18),GSQRD,  UP(10),US(10),RO,QP,QS, VP,VS,VD(10),QL:
      INCLUDE 'auxmod.inc'
C     auxmod.inc
C
C.......................................................................
C
C     Other auxiliary storage locations:
      INTEGER NSRF1,ISRF(2)
      INTEGER ISRAUX,ISBAUX,ISB1,ISRF2,ISB2,ICB2,K0,K1,K2,K,I
      REAL XAUX,YAUX(MY),DAUX(MY),X01,Y01(MY),D01(MY)
C
C     ISB1,X01,Y01,D01,VD... Used only to check for the location of the
C             point X1.
C
C.......................................................................
C
C     (c) Check for crossing the coordinate boundaries of the
C     computational volume:
      DO 31 I=1,NBOUND
        K=IABS(KBOUND(I))
        IF((KBOUND(I).GT.0.AND.Y(K).LT.BOUND(I)).OR.
     *     (KBOUND(I).LT.0.AND.Y(K).GT.BOUND(I))) THEN
          IY(6)=100+I
          FAUX(1)=BOUND(I)
          CALL CROSS(SRFC2,IY(6),K,K,NDIM,ERR,X1,Y1,D1,X2,Y2,D2,
     *                                              X,Y,D,XB,YB,DB,FAUX)
        END IF
   31 CONTINUE
      IF(IY(6).GT.100) THEN
        IY(7)=IY(4)
        IY(8)=IABS(IY(5))
      END IF
      ISB1=IY(4)
C
      XB=X
      DO 39 I=1,NDIM
        YB(I)=Y(I)
        DB(I)=D(I)
   39 CONTINUE
C
C     (d) Check for crossing the boundary of the complex block:
C     Note: IY(4)=ISB1, IY(5)=ICB1, IY(6)=ISRF.
      IF(NOUGHT.EQ.0) THEN
        NSRF1=1
      ELSE
        NSRF1=2
        ISRF(2)=NOUGHT
      END IF
      IF(X.GT.X0) THEN
        ISRF(1)=0
        CALL BLOCKS(Y(KDIM1),NSRF1,ISRF,IY(4),ISRF2,ISB2,ICB2)
        IF(ISRF2.NE.0) THEN
C         Boundary of the simple block is crossed
C         Note: in this routine, unlike in the paper on C.R.T., the
C         point of intersection with the boundary of the simple block is
C         found even if the boundary of the complex block is not
C         crossed.
C         (d1)
          ISRAUX=IY(6)
          ISBAUX=ISB2
          XAUX=X
          DO 41 I=1,NDIM
            YAUX(I)=Y(I)
            DAUX(I)=D(I)
   41     CONTINUE
C         5.30:
          CALL SRFC2(IABS(ISRF2),Y(KDIM1),FAUX)
C         Following loop is included to avoid infinite repeating of the
C         steps (d2) and (d3) of the algorithm
          DO 47 K=1,100
C           (d2)
            IY(6)=ISRF2
C           (d3)
C           Check for the location of the point X1 before calling cross
            IF(X1.LE.X0) THEN
C             X1 may be located outside the complex block
              CALL SRFC2(IABS(IY(6)),Y1(KDIM1),VD)
              IF(VD(1)*FAUX(1).GT.0.) THEN
C               Points X1 and X are located outside the complex block,
C               looking for the point X01 between X0 and X, situated
C               inside the complex block.
                X01=X0
   42           CONTINUE
                CALL HIVD2(KDIM2-KDIM1+1,X1,Y1(KDIM1),D1(KDIM1),X2,
     *                    Y2(KDIM1),D2(KDIM1),X01,Y01(KDIM1),D01(KDIM1))
                CALL SRFC2(IABS(IY(6)),Y01(KDIM1),VD)
                IF(VD(1)*FAUX(1).LE.0.) THEN
C                 Point X01 is likely located inside the complex block
                  ISRF(1)=0
                  CALL BLOCKS(Y01(KDIM1),NSRF1,ISRF,ISB1,ISRF2,K1,K2)
                  IF(ISRF2.EQ.0) THEN
C                   Point X01 is located inside the simple block ISB1,
C                   point X1 may be replaced by X01.
                    CALL HIVD2(NDIM,X1,Y1,D1,X2,Y2,D2,X01,Y01,D01)
                    CALL CROSS(SRFC2,IABS(IY(6)),KDIM1,KDIM2,NDIM,ERR,
     *                        X01,Y01,D01,X2,Y2,D2,X,Y,D,XB,YB,DB,FAUX)
                    GO TO 43
                  END IF
                ELSE
C                 Trying a new point X01
                  IF(X01.EQ.X0) THEN
                    X01=X01+ERR
                  ELSE
                    X01=X01+(X01-X0)
                  END IF
                  IF(X01.LT.X) THEN
                    GO TO 42
                  END IF
                END IF
              END IF
            END IF
            CALL CROSS(SRFC2,IABS(IY(6)),KDIM1,KDIM2,NDIM,ERR,
     *                            X1,Y1,D1,X2,Y2,D2,X,Y,D,XB,YB,DB,FAUX)
   43       CONTINUE
C           X and XB are the approximations of the point of intersection
C           with the surface IY(6)
            ISRF(1)=0
            CALL BLOCKS(Y(KDIM1),NSRF1,ISRF,IY(4),ISRF2,ISB2,ICB2)
            IF(ISRF2.EQ.IY(6).AND.X.NE.X1) THEN
C             587
              CALL ERROR('587 in CDE: Boundary point out of block')
C             This error should not appear.  Contact the authors.
            END IF
C           5.30:
            IF(ISRF2.NE.0) THEN
              CALL SRFC2(IABS(ISRF2),Y(KDIM1),FAUX)
            END IF
            IF(ISRF2.NE.0.AND.ISRF2.NE.IY(6).AND.
     *         FAUX(1)*(D(KDIM1)*FAUX(2)+D(KDIM1+1)*FAUX(3)
     *                                    +D(KDIM2)*FAUX(4)).GT.0.) THEN
C             (d3-i)
C             Point X is not situated at the boundary of the simple
C             block.  It is separated from the simple block IY(4) by
C             the surface ISRF2 situated before (not after) the point X.
C             Go to (d2)
            ELSE
C             X is situated at the boundary of the simple block IY(4)
              ISRF(1)=0
              CALL BLOCKS(YB(KDIM1),NSRF1,ISRF,IY(4),ISRF2,ISB2,ICB2)
              IF(ISB2.EQ.IY(4)) THEN
                CALL SRFC2(IABS(IY(6)),YB(KDIM1),FAUX)
                IF(FAUX(1).NE.0.) THEN
C                 582
                  CALL ERROR('582 in CDE: Excluded program branch')
C                 This error should not appear.  Contact the authors.
                END IF
C               Point XB is situated exactly at the surface IY(6)
                ISRF(1)=IY(6)
                CALL BLOCKS(YB(KDIM1),NSRF1,ISRF,IY(4),ISRF2,ISB2,ICB2)
              END IF
C             Near the edge of a simple block, two different surfaces
C             bounding simple block IY(4) may be situated between
C             points X and XB.  In such a case, only one of them
C             separates simple block IY(4) from ISB2 and IY(6) may
C             be the second.
              ISRF(1)=IY(6)
              CALL BLOCKS(Y(KDIM1),NSRF1,ISRF,IY(4),ISRF2,K1,K2)
              IF(ISRF2.EQ.0.AND.ISB2.NE.K1) THEN
C               Surface IY(6) does not form the interface between
C               simple blocks IY(4) and ISB2:
                ISRF(1)=IY(6)
                CALL BLOCKS(YB(KDIM1),NSRF1,ISRF,IY(4),ISRF2,K1,K2)
                IF(ISRF2.NE.0) THEN
C                 Surface ISRF2 separates the point XB from the simple
C                 block IY(4):
                  ISRF(1)=ISRF2
                  CALL BLOCKS(Y(KDIM1),NSRF1,ISRF,IY(4),K0,K1,K2)
                  IF(K0.EQ.0.AND.ISB2.EQ.K1) THEN
C                   Surface ISRF2 forms the interface between simple
C                   blocks IY(4) and ISB2:
                    IY(6)=ISRF2
                  END IF
                END IF
              END IF
              IF(ICB2.EQ.IABS(IY(5))) THEN
C               (d3-ii)
C               X,Y is situated at the boundary of the simple block
C               but not situated at the boundary of the complex block
                IY(4)=ISB2
                X=XAUX
                DO 45 I=1,NDIM
                  Y(I)=YAUX(I)
                  D(I)=DAUX(I)
   45           CONTINUE
                IF(ISB2.EQ.ISBAUX) THEN
C                 Boundary of the complex block has not been crossed
C                 during the last step of numerical integration
                  XB=X
                  DO 46 I=1,NDIM
                    YB(I)=Y(I)
                    DB(I)=D(I)
   46             CONTINUE
                  IY(6)=ISRAUX
                  GO TO 49
                END IF
                ISRF(1)=0
                CALL BLOCKS
     *                    (YAUX(KDIM1),NSRF1,ISRF,IY(4),ISRF2,ISB2,ICB2)
                IF(ISRF2.EQ.0) THEN
C                 ISRF2 can be zero only if ISBAUX.EQ.0:
                  IF(ISBAUX.EQ.0) THEN
C                   5.30:
                    CALL BLOCKS
     *                    (Y(KDIM1),NSRF1,ISRF,IY(4),ISRF2,ISB2,ICB2)
                    IF(ISRF2.EQ.0) THEN
                      GO TO 49
                    END IF
                  ELSE
                    WRITE(*,'(20(A,I3))') ' ISRFC1=',ISRAUX,
     *                                    ' ISB1=',ISBAUX,
     *                                    ' ICB1=',IY(5),
     *                                    ' ISRFC2=',ISRF2,
     *                                    ' ISB2=',IY(4),
     *                                    ' ICB2=',ICB2
C                   583
                    CALL ERROR('583 in CDE: Excluded program branch')
C                   This error should not appear.  Contact the authors.
                  END IF
                END IF
C               5.30:
                CALL SRFC2(IABS(ISRF2),YAUX(KDIM1),FAUX)
C               Go to (d2)
              ELSE
C               (d3-iii)
C               X is situated at the boundary of the simple block and
C               X is situated at the boundary of the complex block
                GO TO 48
              END IF
            END IF
   47     CONTINUE
C           581
            CALL ERROR('581 in CDE: Too many fictitious interfaces')
C           More than 100 fictitious interfaces crossed during one
C           step of the numerical integration.
C           This error should not appear.  Contact the authors.
   48     CONTINUE
          IY(7)=ISB2
          IY(8)=ICB2
        END IF
   49   CONTINUE
      END IF
C
C     (e) Check for crossing the end surfaces
      DO 51 I=1,NEND
        IF(IABS(KEND(I)).GT.NSRFC()) THEN
          CALL SRFC2(IABS(KEND(I)),Y(KDIM1),FAUX)
          IF(FAUX(1)*FLOAT(KEND(I)).LE.0.)  THEN
            IY(6)=IABS(KEND(I))
            CALL CROSS(SRFC2,IY(6),KDIM1,KDIM2,NDIM,ERR,
     *                            X1,Y1,D1,X2,Y2,D2,X,Y,D,XB,YB,DB,FAUX)
          END IF
        END IF
   51 CONTINUE
C
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE CROSS(SRFC2,ISRF,KDIM1,KDIM2,NDIM,
     *                        ERR,X1,Y1,D1,X2,Y2,D2,XA,YA,DA,XB,YB,DB,F)
      EXTERNAL SRFC2
      INTEGER ISRF,KDIM1,KDIM2,NDIM
      REAL ERR,X1,Y1(NDIM),D1(NDIM),X2,Y2(NDIM),D2(NDIM)
      REAL XA,YA(NDIM),DA(NDIM),XB,YB(NDIM),DB(NDIM),F(10)
C
C This subroutine finds the point of intersection of a curve with
C a surface (see C.R.T.5.8.4b).  The curve is parametrized by an
C independent variable X and evaluated by the Hermite interpolation from
C the two given points.  The surface is specified in an implicit way by
C subroutine SRFC2 which is described elsewhere, or may coincide with an
C isosurface of a computed quantity (e.g. with a coordinate plane).
C
C Input:
C     SRFC2...Name of the external procedure evaluating the function
C             describing the surface ISRF, for ISRF.LE.100.
C     ISRF... Index of the surface.
C             For ISRF.LE.100:
C               The surface coincides with the zero isosurface of the
C               function no. ISRF evaluated by the subroutine SRFC2.
C             For ISRF.GT.100:
C               The surface coincides with an isosurface Y(KDIM1)=F(1).
C     KDIM1,KDIM2... Indices of quantities on which the function
C             describing the surface is dependent:
C             For ISRF.LE.100:
C               KDIM1-th to KDIM2-th elements of arrays Y1, Y2, YA, and
C               YB contain coordinates.  KDIM2 is assumed to be KDIM1+1
C               (2 coordinates) or KDIM1+2 (3 coordinates).
C               The function describing the surface ISRF is determined
C               by the subroutine SRFC2.
C             For ISRF.GT.100:
C               The surface coincides with an isosurface Y(KDIM1)=F(1).
C               KDIM2 is assumed to equal KDIM1.
C             When searching for the point of intersection, only
C             quantities Y(KDIM1:KDIM2) and their derivatives are
C             interpolated along the curve.
C     NDIM... Dimension of arrays Y1,D1,Y2,D2,YA,DA,YB,DB.
C             At the point of intersection, quantities Y(1:NDIM) and
C             their derivatives are interpolated.
C     ERR...  Maximum error in independent variable for the
C             determination of the point of intersection.
C     X1...   Independent variable corresponding to the first point
C             given for the interpolation of the curve.
C     Y1...   Array containing dependent variables at point X1.
C             Y1(KDIM1) to Y1(KDIM2) must contain the coordinates of
C             point X1.
C     D1...   Array containing the derivatives of the dependent
C             variables at point X1.
C     X2...   Independent variable corresponding to the second point
C             given for the interpolation of the curve.
C     Y2...   Array containing dependent variables at point X2.
C             Y2(KDIM1) to Y2(KDIM2) must contain the coordinates of
C             point X2.
C     D2...   Array containing the derivatives of the dependent
C             variables at point X2.
C     XA...   Independent variable corresponding to the point of the
C             curve at which the function specifying the surface has the
C             opposite sign than at X1. Then the point of intersection
C             is being found between the points X1 and XA.  The found
C             approximations XA and XB of the point of intersection are
C             situated close to the surface, XA at the same side as the
C             given point X1, XB at the oposite side than the given
C             point X1.  If, accidentally, the function has the same
C             sign at the points X1 and XA, the value at X1 is assumed
C             to equal zero.  Then XA=X1 and XB=X1 on output.
C     YA...   Array of the dimension at least NDIM.  YA(KDIM1) to
C             YA(KDIM2) must contain the coordinates of the point.
C             Other storage locations may be undefined.
C     DA...   Array of the dimension at least NDIM, containing the
C             derivatives of the dependent variables at point X.
C             DA(KDIM1) to DA(KDIM2) must contain the derivatives of
C             coordinates with respect to X, at the point XA.  Other
C             storage locations may be undefined.
C     YB...   Array of the dimension at least NDIM.
C     DB...   Array of the dimension at least NDIM.
C     F...    For ISRF.LE.100:
C               Array containing the value, and at least first
C               derivatives of the function ISRF specifying the surface,
C               at point XA.
C             For ISRF.GT.100:
C               The surface coincides with an isosurface Y(KDIM1)=F(1).
C               F(2) to F(10) need not be defined.
C ISRF, KDIM1, KDIM2, NDIM, ERR, X1, Y1, D1, X2, Y2, D2 are unaltered.
C
C Output:
C     XA...   Independent variable corresponding to the approximation of
C             the point of intersection, situated close to the surface
C             at the same side as the given point X1.  XA=X1 if the
C             function had at point X1 the same sign as F(1) on input.
C     YA...   Array containing dependent variables at the point XA.
C     DA...   Array containing the derivatives of the dependent
C             variables at the point XA.
C     XB...   Independent variable corresponding to the approximation of
C             the point of intersection, situated close to the surface
C             at the oposite side than the given point X1.  XB=X1 if the
C             function had at point X1 the same sign as F(1) on input.
C     YB...   Array containing dependent variables at the point XB.
C     DB...   Array containing the derivatives of the dependent
C             variables at the point XB.
C     F...    Undefined.
C
C Subroutines and external functions required:
      EXTERNAL HIVD2
C     SRFC2 and subsequent routines... File 'srfc.for' and subsequent
C              files (like 'val.for' and 'fit.for') - dummy argument.
C     HIVD2... This file.
C
C Date: 1996, July 10
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER I,ITER
      REAL FAOLD,X,FX,DFX,FB,DFB,FA,DFA,XC,XCB
*old  REAL                                    ,XCA
*!!!  REAL XD,XD2,XE,AUX0,AUX1,AUX2,AUX3
*out  real xam(50),fam(50),dfam(50),xbm(50),fbm(50),dfbm(50)
*out  integer icross,jcross,mcross
*out  data    icross,jcross,mcross/0,0,0/
C
C.......................................................................
C
C     Initial values:
      IF(ISRF.LE.100) THEN
        FA=F(1)
        DFA=F(2)*DA(KDIM1)
        DO 1 I=1,KDIM2-KDIM1
          DFA=DFA+F(2+I)*DA(KDIM1+I)
    1   CONTINUE
      ELSE
        FA=YA(KDIM1)-F(1)
        DFA=DA(KDIM1)
      END IF
      FAOLD=FA
      XB=X1
      X =X1
      DO 2 I=KDIM1,KDIM2
        YB(I)=Y1(I)
        DB(I)=D1(I)
    2 CONTINUE
C     In the beginning, both points XB and X are identical with point X1
C
C     Check for zero intervals:
      IF(XA.EQ.X1) THEN
        DO 3 I=1,NDIM
          YA(I)=Y1(I)
          DA(I)=D1(I)
          YB(I)=Y1(I)
          DB(I)=D1(I)
    3   CONTINUE
        RETURN
      END IF
      IF(X2.EQ.X1) THEN
C       585
        CALL ERROR('585 in CROSS: Zero-length interval')
C       Two points X1 and X2 used for interpolation of the given
C       line (e.g. the velocity isoline, or the ray), when
C       searching for its intersection with the interface, are
C       identical and no interpolation is possible.  Likely a bug
C       in the procedure calling this subroutine.
      END IF
C
C     Iterations:
      DO 9 ITER=1,50
C
C       Functional value and derivative at X:
        IF(ISRF.LE.100) THEN
          CALL SRFC2(ISRF,YB(KDIM1),F)
          FX=F(1)
          DFX=F(2)*DB(KDIM1)
          DO 5 I=1,KDIM2-KDIM1
            DFX=DFX+F(2+I)*DB(KDIM1+I)
    5     CONTINUE
        ELSE
          FX=YB(KDIM1)-F(1)
          DFX=DB(KDIM1)
        END IF
C
C       Selection of points:
        IF(FX.EQ.0.) THEN
          XA=X
          FA=FX
        ELSE IF(FA*FX.GE.0.) THEN
          IF(ITER.EQ.1) THEN
C           Input points X1 and XA are situated at the same side
            IF(FA.EQ.0.) THEN
              X=XA
              FX=FA
            ELSE
              XA=X
              FA=FX
            END IF
          ELSE
C           Here FA, FX and FB should be non-zero due to previous checks
            XA=XB
            FA=FB
            DFA=DFB
          END IF
        END IF
        XB=X
        FB=FX
        DFB=DFX
C
*out    xam(iter)=xa
*out    fam(iter)=fa
*out    dfam(iter)=dfa
*out    xbm(iter)=xb
*out    fbm(iter)=fb
*out    dfbm(iter)=dfb
c
C       New point or end of iterations:
        IF(ABS(XB-XA).LE.ERR) THEN
C         Point of intersection is found within the specified error err
C         *** end of iterations ***
          IF(FB*FAOLD.LT.0.) THEN
C           Point XA is situated at the other side of the surface than
C           the point X1 - changing XA and XB
            X =XA
            XA=XB
            XB=X
          END IF
          IF((XB-XA)*(XB-X1).LT.0.) THEN
C           Point XB is closer to X1 than point XA
            IF(FA*FB.LT.0.) THEN
C             Points XA and XB cannot be changed
C             586
              CALL ERROR('586 in CROSS: Reverse order of points')
C             A pair of close points XA and XB situated at different
C             sides of the surface has been found, but point XA situated
C             at the same side as point X1 is far from X1 than point XB.
C             This error should not appear.  Contact the authors.
            ELSE
C             Points XA and XB are situated at the same side of the
C             surface (and, hopefully, FA and FB should be zero)
C             - changing XA and XB
              X =XA
              XA=XB
              XB=X
            END IF
          END IF
          CALL HIVD2(NDIM,X1,Y1,D1,X2,Y2,D2,XA,YA,DA)
          CALL HIVD2(NDIM,X1,Y1,D1,X2,Y2,D2,XB,YB,DB)
*out      if(icross.eq.0) then
*out        open(57,file='cross.out')
*out      end if
*out      icross=icross+1
*out      jcross=jcross+iter-1
*out      mcross=max0(mcross,iter-1)
*out      if(mod(icross,100).eq.0) then
*out        write(57,*) icross,jcross,mcross,float(jcross)/float(icross)
*out      end if
*out      if(iter.gt.20) then
*out        open(58,file='error.out')
*out        do 7 I=1,iter
*out          write(58,*) xam(i),xbm(i),fam(i),fbm(i),dfam(i),dfbm(i)
*out7       continue
*out        endfile  (58)
*out        backspace(58)
*out        pause 'Error new in CROSS: More than 20 iterations'
*out      end if
          RETURN
        END IF
C
cccc    (a) Odstranit *old
cccc    (b) Odstranit *old, odstranit *???
cccc    (c) Odstranit *old, odstranit *???, odstranit *new (nedulezite)
cccc    (d) Vratit *old, vratit *???, vratit *new
cccc    (e) Odstranit *!!!
cccc    (f) Odstranit *!!!, odstranit *old
cccc
C       New approximation:
        IF(MOD(ITER,2).EQ.1) THEN
C         Regula falsi:
          X=(FA*XB-FB*XA)/(FA-FB)
C!!!      CHECKING THE ACCURACY OF THE REGULA FALSI METHOD:
*!!!      IF(ITER.EQ.1) THEN
C!!!        COEFFICIENTS OF THE CUBIC TAYLOR EXPANSION AT XC=(XA+XB)/2.:
*!!!        XD=XB-XA
*!!!        XD2=XD*XD/2.
*!!!        AUX3=(DFB+DFA)/2.
*!!!        AUX2=(DFB-DFA)/XD
*!!!        AUX1=(FB-FA)/XD
*!!!        AUX0=(FB+FA)/2.-AUX2*XD2
*!!!        AUX1=1.5*AUX1-0.5*AUX3
*!!!        AUX3=(AUX3-AUX1)/XD2
C!!!        COEFFICIENTS OF THE QUADRATIC TAYLOR EXPANSION AT X:
*!!!        XD=X-XC
*!!!        XD=X-XC
*!!!        XD2=XD*XD/2.
*!!!        XE=SIGN(1.,AUX1)
*!!!        AUX0=AUX0+AUX1*XD+AUX2*XD2+AUX3*XD*XD2/3.
*!!!        AUX1=AUX1+AUX2*XD+AUX3*XD2
*!!!        IF(ABS(AUX0).GT.0.5*ERR*ABS(AUX1)) THEN
C!!!          THE ACCURACY SHOULD BE IMPROVED:
*!!!          AUX2=AUX2+AUX3*XD
*!!!          AUX3=AUX1**2-2.*AUX0*AUX2
C!!!          HERE 2**(24/2)=4096 ASSUMES 24 BIT FLOATING-POINT ACCURACY
*!!!          IF(4096.*ABS(AUX0*AUX2).GT.AUX1**2.AND.AUX3.GE.0.) THEN
*!!!            XE=X+(XE*SQRT(AUX3)-AUX1)/AUX2
*!!!          ELSE
*!!!            XE=X-AUX0/AUX1
*!!!          END IF
*!!!          IF((XE-XA)*(XE-XB).LE.0.) THEN
*!!!            X=XE
*!!!          END IF
*!!!        END IF
*!!!      END IF
        ELSE
C         Modified Newton-Raphson:
          XC=(XA+XB)/2.
*old      XCA=XA-FA/DFA+SIGN(ERR/50.,XA-XB)
          XCB=XB-FB/DFB+SIGN(ERR/50.,XA-XB)
*old      IF((XCA-XC)*(XCA-XB).LT.0.) THEN
*old        IF((XCB-XC)*(XCB-XB).LT.0.) THEN
*old          IF(ABS(XCA-XB).LT.ABS(XCB-XB)) THEN
*old            X=XCA
*old          ELSE
*old            X=XCB
*old          END IF
*old        ELSE
*old          X=XCA
*old        END IF
*old      ELSE
            IF((XCB-XC)*(XCB-XB).LT.0.) THEN
              X=XCB
            ELSE
              X=XC
*???          IF(ABS(XCB-XB).LT.SQRT(ABS(XC-XB)*ERR)) THEN
C               Attempt to halve the number of iterations
*???            X=XB+SIGN(SQRT(ABS(XC-XB)*ERR),XA-XB)
*new            IF((X-XC)*(X-XB).GE.0.) THEN
*new              X=XC
*new              pause 'New warning in CROSS'
*new            END IF
*???          END IF
            END IF
*old      END IF
        END IF
C
C       Interpolation of the ray:
        CALL HIVD2(KDIM2-KDIM1+1,X1,Y1(KDIM1),D1(KDIM1),
     *                     X2,Y2(KDIM1),D2(KDIM1),X,YB(KDIM1),DB(KDIM1))
C
    9 CONTINUE
C     End of loop for iterations
*out  open(58,file='error.out')
*out  do 8 iter=1,50
*out    write(58,*) xam(iter),xbm(iter),fam(iter),fbm(iter),
*out *              dfam(iter),dfbm(iter)
*out8 continue
C       584
        CALL ERROR('584 in CROSS: Too many iterations')
C       More than 50 iterations when determining the point of
C       intersection of the ray with a surface.  This may be
C       caused by too small upper error bound UEB in the input
C       data 'dcrt.dat' - less than rounding errors.  Usually,
C       this error should not appear.  Contact the authors.
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE HIVD2(NDIM,X1,Y1,D1,X2,Y2,D2,X,Y,D)
      INTEGER NDIM
      REAL X1,Y1(NDIM),D1(NDIM),X2,Y2(NDIM),D2(NDIM)
      REAL X,Y(NDIM),D(NDIM)
C
C This subroutine performs Hermite interpolation of a vector and its
C derivatives using functional values and derivatives at 2 given points.
C
C Input:
C     NDIM... Dimension of arrays Y1,D1,Y2,D2,Y,D (the number of
C             independent variables).
C     X1...   Independent variable corresponding to the first given
C             point.
C     Y1...   Array containing functional values at the first given
C             point.
C     D1...   Array containing the derivatives at the first given point.
C     X2...   Independent variable corresponding to the second given
C             point.
C     Y2...   Array containing functional values at the second given
C             point.
C     D2...   Array containing the derivatives at the second given point
C     X...    Independent variable of the point at which the
C             interpolated vector is to be evaluated.
C None of the input parameters are altered.
C
C Output:
C     Y...    Array containing interpolated functional values at X.
C     D...    Array containing the derivatives of the interpolated
C             functional values at X.
C
C No subroutines and external functions required.
C
C Date: 1989, October 20
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER I
      REAL A,B,A1,A2,B1,B2,DA1,DB1,DB2
C
C     Substitutions:
      A=(X-X2)/(X1-X2)
      B=(A-1.)*A
C
C     Basic functions:
      A1=(A-B-B)*A
      A2=1.-A1
      B1=B*(X-X2)
      B2=B*(X-X1)
C     Derivatives of basic functions:
      DA1=6.*B/(X2-X1)
      DB1=3.*B+A
      DB2=3.*B+1.-A
C
C     Interpolation:
      DO 1 I=1,NDIM
        Y(I)=A1*Y1(I)+A2*Y2(I)+B1*D1(I)+B2*D2(I)
        D(I)=DA1*(Y1(I)-Y2(I))+DB1*D1(I)+DB2*D2(I)
    1 CONTINUE
C
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE SMVPRD(G,A1,A2,A3,B1,B2,B3)
      REAL G(6),A1,A2,A3,B1,B2,B3
C
C This subroutine is designed to evaluate symmetric matrix by vector
C product.  It may be, e.g., called after the invocation of the METRIC
C subroutine to transform the covariant vectors to the contravariant
C ones and vice versa.
C
C Input:
C   G...      Array containing components G11, G12, G22, G13, G23, G33
C             of the 3*3 symmetric matrix.
C   A1,A2,A3... Components of the 3-vector.
C
C Output:
C   B1,B2,B3... Components of vector B=G*A.
C
C No subroutines and external functions required.
C
C Date: 1989, September 5
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
      B1=G(1)*A1+G(2)*A2+G(4)*A3
      B2=G(2)*A1+G(3)*A2+G(5)*A3
      B3=G(4)*A1+G(5)*A2+G(6)*A3
      RETURN
      END
C
C=======================================================================
C