C SUBROUTINE FILE 'MEANS.FOR' CONTAINING SOME UTILITY PROGRAMS HELPFUL
C WHEN DEALING WITH THE MODEL.
C
C BY VLASTISLAV CERVENY, LUDEK KLIMES, IVAN PSENCIK
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 OF THE ALGORITHM.
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     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     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
C DATE: 1992, NOVEMBER 2
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 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 THE
C             BOUNDARIES X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX OF THE
C             COMPUTATIONAL VOLUME.
C     KBOUND..INDICES OF THE QUANTITIES CORRESPONDING TO THE ISOSURFACES.
C             THE COMPUTATIONAL VOLUME IS LIMITED BY THE 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... UNDEFINED.
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 ERROR MESSAGES:
C     581...  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.
C     582,583... EXCLUDED PROGRAM BRANCH:
C             THIS ERROR SHOULD NOT APPEAR.  CONTACT THE AUTHORS.
C
C DATE: 1992, DECEMBER 31
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     AUXILIARY STORAGE LOCATIONS FOR LOCAL MODEL PARAMETERS:
      REAL G(12),GAMMA(18),GSQRD,FAUX(10)
      REAL UP(10),US(10),RO,QP,QS,VP,VS,VD(10),QL
C     THESE AUXILIARY VARIABLES AND ARRAYS NEED NOT BE LOCATED IN A
C     COMMON BLOCK.  THERE IS NO REASON TO LOCATE THEM IN THE AUXILIARY
C     COMMON BLOCK /AUXMOD/ BUT TO SHARE THE MEMORY.
      COMMON/AUXMOD/ G,GAMMA,GSQRD,FAUX,UP,US,RO,QP,QS,VP,VS,VD,QL
C
C     G,GAMMA,GSQRD... SEE SUBROUTINE METRIC OF THE FILE 'METRIC.FOR'.
C     FAUX... AUXILIARY ARRAY TO STORE A FUNCTIONAL VALUE AND ITS
C             DERIVATIVES.
C     UP,US,RO,QP,QS... SEE SUBROUTINE PARM2 OF THE FILE 'PARM.FOR'.
C     UP,US,QP,QS,VP,VS,VD,QL... SEE SUBROUTINE VELOC OF THE FILE
C             'MODEL.FOR'.
C
C.......................................................................
C
C     OTHER AUXILIARY STORAGE LOCATIONS:
      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
C     (D) CHECK FOR CROSSING THE BOUNDARY OF THE COMPLEX BLOCK:
C     NOTE: IY(4)=ISB1, IY(5)=ICB1, IY(6)=ISRF.
      IF(X.GT.X0) THEN
        CALL BLOCK(Y(KDIM1),NOUGHT,IY(4),ISRF2,ISB2,ICB2,FAUX)
   40   CONTINUE
        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         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
                  CALL BLOCK(Y01(KDIM1),NOUGHT,ISB1,ISRF2,K1,K2,VD)
                  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)
            CALL BLOCK(Y(KDIM1),NOUGHT,IY(4),ISRF2,ISB2,ICB2,FAUX)
            IF(ISRF2.NE.0.AND.(ISRF2.NE.IY(6).OR.X.NE.X1).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)
              CALL BLOCK(YB(KDIM1),NOUGHT,IY(4),ISRF2,ISB2,ICB2,FAUX)
              IF(ISB2.EQ.IY(4)) THEN
                CALL SRFC2(IABS(IY(6)),YB(KDIM1),FAUX)
                IF(FAUX(1).NE.0.) THEN
                  PAUSE 'ERROR 582 IN CDE: EXCLUDED PROGRAM BRANCH'
                  STOP
                END IF
C               POINT XB IS SITUATED EXACTLY AT THE SURFACE IY(6)
                CALL BLOCK(YB(KDIM1),IY(6),IY(4),ISRF2,ISB2,ICB2,FAUX)
                IF(NOUGHT.NE.0.AND.NOUGHT.EQ.ISRF2) THEN
                  K1=ISB2
                  CALL BLOCK(YB(KDIM1),IY(6),K1,ISRF2,ISB2,ICB2,FAUX)
                END IF
              END IF
              IF(NOUGHT.EQ.0) THEN
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.
                CALL BLOCK(Y(KDIM1),IY(6),IY(4),ISRF2,K1,K2,FAUX)
                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:
                  CALL BLOCK(YB(KDIM1),IY(6),IY(4),ISRF2,K1,K2,FAUX)
                  IF(ISRF2.NE.0) THEN
C                   SURFACE ISRF2 SEPARATES THE POINT XB FROM THE SIMPLE
C                   BLOCK IY(4):
                    CALL BLOCK(Y(KDIM1),ISRF2,IY(4),K0,K1,K2,FAUX)
                    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
              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 46 I=1,NDIM
                  Y(I)=YAUX(I)
                  D(I)=DAUX(I)
   46           CONTINUE
                IF(ISB2.EQ.ISBAUX) THEN
C                 BOUNDARY OF THE COMPLEX BLOCK HAS NOT BEEN CROSSED
C                 DURING THE LAST STEP OF NUMERICAL INTEGRATION
                  IY(6)=ISRAUX
                  GO TO 49
                END IF
                CALL BLOCK(YAUX(KDIM1),NOUGHT,IY(4),ISRF2,ISB2,ICB2
     *                                                            ,FAUX)
                IF(ISRF2.EQ.0) THEN
C                 ISRF2 CAN BE ZERO ONLY IF ISB2.EQ.ISBAUX:
                  PAUSE 'ERROR 583 IN CDE: EXCLUDED PROGRAM BRANCH'
                  STOP
                END IF
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
          PAUSE 'ERROR 581 IN CDE: TOO MANY FICTITIOUS INTERFACES'
          STOP
   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
      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 ELSWHERE, 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, ACCIDENTALY, THE FUNCTION HAS THE SAME SIGN
C             AT THE POINTS X1 AND XA, THE VALUE AT X1 IS ASSUMED TO
C             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 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 ERROR MESSAGES:
C     584...  TOO MANY ITERATIONS:
C             MORE THAN 20 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.
C     585...  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.
C
C DATE: 1994, JANUARY 9
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     AUXILIARY STORAGE LOCATIONS:
      INTEGER I,ITER
      REAL X,FX,DFX,FB,DFB,FA,DFA,XC,XCB,XCA
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
      XB=X1
      X =X1
      DO 2 I=KDIM1,KDIM2
        YB(I)=Y1(I)
        DB(I)=D1(I)
    2 CONTINUE
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
        PAUSE 'ERROR 585 IN CROSS: ZERO-LENGTH INTERVAL'
        STOP
      END IF
C
C     ITERATIONS:
      DO 9 ITER=1,20
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
        ELSE IF(FA*FX.GE.0.) THEN
          IF(ITER.EQ.1) THEN
            IF(FA.EQ.0.) THEN
              X=XA
            ELSE
              FX=0.
            END IF
          ELSE
            XA=XB
            FA=FB
            DFA=DFB
          END IF
        END IF
        XB=X
        FB=FX
        DFB=DFX
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(ABS(XB-X1).LT.ABS(XA-X1)) THEN
C           POINT XA IS SITUATED AT THE OTHER SIDE OF THE SURFACE THAN
C           THE POINT X1
            XB=XA
            XA=X
          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)
          RETURN
        END IF
C
C       NEW APPROXIMATION:
        IF(MOD(ITER,2).EQ.1) THEN
C         REGULA FALSI:
          X=(FA*XB-FB*XA)/(FA-FB)
        ELSE
C         MODIFIED NEWTON-RAPHSON:
          XC=(XA+XB)/2.
          XCA=XA-FA/DFA+SIGN(ERR/50.,XA-XB)
          XCB=XB-FB/DFB+SIGN(ERR/50.,XA-XB)
          IF((XCA-XC)*(XCA-XB).LT.0.) THEN
            IF((XCB-XC)*(XCB-XB).LT.0.) THEN
              IF(ABS(XCA-XB).LT.ABS(XCB-XB)) THEN
                X=XCA
              ELSE
                X=XCB
              END IF
            ELSE
              X=XCA
            END IF
          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)
              END IF
            END IF
          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
      PAUSE 'ERROR 584 IN CROSS: TOO MANY ITERATIONS'
      STOP
      END
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
      SUBROUTINE SMVPRD(G,A1,A2,A3,B1,B2,B3)
      REAL G(6),A1,A2,A3,B1,B2,B3
C
C THIS SUBROUTINE IS AN AUXILIARY ROUTINE TO FCT. IT EVALUATES SYMMETRIC
C MATRIX BY VECTOR PRODUCT.
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
