C
      SUBROUTINE INIT(PN,V)
C
C  ROUTINE FOR THE CALCULATION OF THE PHASE VELOCITIES FOR THE
C  SPECIFIED UNIT NORMAL TO THE PHASE FRONT
C  SOLUTION OF THE EIGENVALUE PROBLEM
C
      DIMENSION C(3,3),PN(3),V(3)
      DOUBLE PRECISION A1,A2,A3,R,S,T,Q,P,TH
      COMMON /AUXI/  IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT,
     1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOU,
     2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,MORI
C
      CALL CHRM1(C,PN,PN)
C
C     COMPUTATION OF THE ROOTS OF THE CHARACTERISTIC EQUATION WITH
C     CARDANS FORMULA
C
C     COMPUTATION OF THE COEFFICIENTS OF THE CHARACTERISTIC EQUATION
C     X**3+R*X**2+S*X+T=0
C
      C11=C(1,1)
      C12=C(1,2)
      C13=C(1,3)
      C22=C(2,2)
      C23=C(2,3)
      C33=C(3,3)
      A1=C22*C33-C23*C23
      A2=C11*C33-C13*C13
      A3=C11*C22-C12*C12
      R=-(C11+C22+C33)
      S=A1+A2+A3
      A2=C12*C33-C13*C23
      A3=C12*C23-C13*C22
      T=-C11*A1+C12*A2-C13*A3
C
C     SOLUTION OF CUBIC EQUATION FOLLOWING THE ALGORITHM
C     FROM NUMERICAL RECIPES
C
      Q=(R*R-3.*S)/9.
      P=(2.*R*R*R-9.*R*S+27.*T)/54.
      A1=Q*Q*Q-P*P
      PI=3.14159265
      IF(A1.LE.0.00001)THEN
        TH=PI
      END IF
      IF(A1.GT.0.00001)THEN
        TH=Q*Q*Q
        TH=P/DSQRT(TH)
        IF(TH.LE.(-1.))TH=-1.
        IF(TH.GE.(1.))TH=1.
        TH=DACOS(TH)
      END IF
      A1=-2.*DSQRT(Q)
      A2=-R/3.
      A3=TH/3.
      V(1)=A1*DCOS(A3)+A2
      A3=(TH+2.*PI)/3.
      V(2)=A1*DCOS(A3)+A2
      A3=(TH+4.*PI)/3.
      V(3)=A1*DCOS(A3)+A2
      X=V(1)
      AUX1=X**3+R*X**2+S*X+T
      X=V(2)
      AUX2=X**3+R*X**2+S*X+T
      X=V(3)
      AUX3=X**3+R*X**2+S*X+T
      IF(IPRINT.GT.2)WRITE(LOU,'(A/1X,3E15.6)')'INIT: PRECISION OF SOLUT
     *IONS OF CUBIC EQUATION ', AUX1,AUX2,AUX3
      RETURN
      END
C
C     *********************************************************
C
      SUBROUTINE INV(A,AINV,DET)
C
C     ROUTINE INV DETERMINESS AN ADJOINT MATRIX OF A REAL 3*3 MATRIX
C
      DIMENSION A(3,3),AINV(3,3)
      COMMON /AUXI/  IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT,
     1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOUT,
     2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,MORI
      COMMON /ZERO/ RNULL
C
      DET=A(1,1)*(A(2,2)*A(3,3)-A(2,3)*A(3,2))-A(2,1)*(A(1,2)*A(3,3)
     1-A(3,2)*A(1,3))+A(3,1)*(A(1,2)*A(2,3)-A(2,2)*A(1,3))
      INULL=0
      IF(ABS(DET).GT.RNULL) GOTO 100
      DO 10 K=1,3
      DO 10 J=1,3
      IF(ABS(A(J,K)).GT.RNULL) INULL=1
      AINV(J,K)=0.0
   10 CONTINUE
      IF(INULL.EQ.0) THEN
        IF(IPRINT.GT.2)WRITE(LOUT,1000) A
      ELSE
        IF(IPRINT.GT.2)WRITE(LOUT,1010)
      END IF
  100 AINV(1,1)=(A(2,2)*A(3,3)-A(3,2)*A(2,3))
      AINV(2,1)=-(A(2,1)*A(3,3)-A(2,3)*A(3,1))
      AINV(3,1)=(A(2,1)*A(3,2)-A(2,2)*A(3,1))
      AINV(1,2)=-(A(1,2)*A(3,3)-A(1,3)*A(3,2))
      AINV(2,2)=(A(1,1)*A(3,3)-A(1,3)*A(3,1))
      AINV(3,2)=-(A(1,1)*A(3,2)-A(1,2)*A(3,1))
      AINV(1,3)=(A(1,2)*A(2,3)-A(1,3)*A(2,2))
      AINV(2,3)=-(A(1,1)*A(2,3)-A(1,3)*A(2,1))
      AINV(3,3)=(A(1,1)*A(2,2)-A(1,2)*A(2,1))
      IF(IPRINT.LE.2)RETURN
      WRITE(LOUT,1030)A
      WRITE(LOUT,1031)AINV
      WRITE(LOUT,1032)DET
      DO 30  K=1,3
      E1=0.0
      E2=0.0
      E3=0.0
      DO 20 J=1,3
      E1=E1+A(J,K)*AINV(1,J)
      E2=E2+A(J,K)*AINV(2,J)
      E3=E3+A(J,K)*AINV(3,J)
  20  CONTINUE
      WRITE(LOUT,1033)E1,E2,E3
  30  CONTINUE
C
 1000 FORMAT(1x,' INV: ALL ELEMENTS OF MATRIX A ZERO',/,3(3F12.6,/))
 1010 FORMAT(1x,' INV :DETERMINANT OF MATRIX A ZERO')
 1030 FORMAT(1X,' INV: A'/3(3F16.8/))
 1031 FORMAT(1X,' INV: AINV'/3(3F16.8/))
 1032 FORMAT(1X,' INV: DET ',F16.8)
 1033 FORMAT(1X,' INV: ',3F16.8)
      RETURN
      END
C
C     *********************************************************
C
      SUBROUTINE OUT(X,Y,DERY,IHLF,NDIM,PRMT)
C
      DIMENSION Y(18),DERY(18),UN(3),POLD(3),PNEW(3),YA(3),YB(3),DEP(6),
     *PRMT(5),DN(3,3),DNG(3,2),XK(3),XKG(3),YK(3),YKG(3),TM(2),
     *XDAUX(3),YDAUX(3)
      COMMON /AUXI/  IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT,
     1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOU,
     2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,MORI
      INTEGER CODE
      COMMON /COD/  CODE(50,2),KREF,KC,ITYPE
      COMMON /DENS/ RHO(20)
      COMMON /INTRF/ Z(1000),SX(350),SY(350),NX(20),NY(20),BRD(6),NINT,
     1    XINTA
      COMPLEX PS
      COMMON /RAY/   AY(28,2000),DS(20,50),KINT(50),HHH(3,3),tmax,
     1               PS(3,7,50),IS(8,50),N,IREF,IND,IND1
      COMMON /RAY2/  DRY(3,2000)
      COMMON /ZERO/ RNULL
      COMMON /VSP/XVSP,YVSP,XNRM,YNRM,ICOD,IVSP
      COMMON /DYN/XDYN(3,3),ydyn(3,3)
      common/appr/ xold,xnew,yold(18),dold(18),ynew(18),dnew(18)
      COMMON /KM/KMAH,SPROLD,TSTOLD
C
      N=N+1
      NTR=20
      IF(IND.EQ.10)GO TO 25
      IF(N.GE.2000) GOTO 152
      IF(X.GE.TMAX)GO TO 155
C
C     TEST OF SATISFACTION OF THE EIKONAL EQUATION AND OF
C     THE DERIVATIVES OF THE EIKONAL EQUATION
C
      IF(IPREC.EQ.1.AND.NDER.GT.1)THEN
        TEST0=DERY(1)*Y(4)+DERY(2)*Y(5)+DERY(3)*Y(6)
        IF(ABS(TEST0).GT.RNULL)THEN
          DO 80 I=4,6
          Y(I)=Y(I)/TEST0
   80     CONTINUE
        END IF
        V2=1./(Y(4)*Y(4)+Y(5)*Y(5)+Y(6)*Y(6))
        TEST3=Y(7)*Y(4)+Y(8)*Y(5)+Y(9)*Y(6)
        IF(ABS(TEST3).GT.RNULL)THEN
          DO 81 I=4,6
          Y(I+3)=Y(I+3)-V2*TEST3*Y(I)
   81     CONTINUE
        END IF
        ALP1=DERY(4)*Y(7)+DERY(5)*Y(8)+DERY(6)*Y(9)
        ALP2=DERY(1)*Y(13)+DERY(2)*Y(14)+DERY(3)*Y(15)
        IF(ABS(ALP1-ALP2).GT.RNULL)THEN
          IF(ABS(ALP1).GT..000001)THEN
            ALP=ALP2/ALP1
            DO 82 I=7,9
            Y(I)=Y(I)*ALP
   82       CONTINUE
          ELSE
            DO 86 I=4,6
            DERY(I+3)=DERY(I+3)-V2*ALP2*Y(I)
   86       CONTINUE
          END IF
        END IF
        TEST4=Y(10)*Y(4)+Y(11)*Y(5)+Y(12)*Y(6)
        IF(ABS(TEST4).GT.RNULL)THEN
          DO 83 I=4,6
          Y(I+6)=Y(I+6)-V2*TEST4*Y(I)
   83     CONTINUE
        END IF
        ALP1=DERY(4)*Y(10)+DERY(5)*Y(11)+DERY(6)*Y(12)
        ALP2=DERY(1)*Y(16)+DERY(2)*Y(17)+DERY(3)*Y(18)
        IF(ABS(ALP1-ALP2).GT.RNULL)THEN
          IF(ABS(ALP1).GT..000001)THEN
            ALP=ALP2/ALP1
            DO 84 I=10,12
            Y(I)=Y(I)*ALP
   84       CONTINUE
          ELSE
            DO 87 I=4,6
            DERY(I+6)=DERY(I+6)-V2*ALP2*Y(I)
   87       CONTINUE
          END IF
        END IF
      END IF
      IF(IPREC.EQ.2.AND.NDER.GT.1)THEN
        TEST0=DERY(1)*Y(4)+DERY(2)*Y(5)+DERY(3)*Y(6)
        TEST1=DERY(1)*Y(13)+DERY(2)*Y(14)+DERY(3)*Y(15)
     1  -DERY(4)*Y(7)-DERY(5)*Y(8)-DERY(6)*Y(9)
        TEST2=DERY(1)*Y(16)+DERY(2)*Y(17)+DERY(3)*Y(18)
     1  -DERY(4)*Y(10)-DERY(5)*Y(11)-DERY(6)*Y(12)
        TEST3=Y(7)*Y(4)+Y(8)*Y(5)+Y(9)*Y(6)
        TEST4=Y(10)*Y(4)+Y(11)*Y(5)+Y(12)*Y(6)
        IF(ABS(TEST0-1.).GT.RNULL)WRITE(LOU,100)TEST0
        IF(ABS(TEST1).GT.RNULL.OR.ABS(TEST2).GT.RNULL)
     1  WRITE(LOU,101)TEST1,TEST2
        IF(ABS(TEST3).GT.RNULL.OR.ABS(TEST4).GT.RNULL)
     1  WRITE(LOU,102)TEST3,TEST4
  100   FORMAT(1X,'P_I*V_I=',E15.5)
  101   FORMAT(1X,'DER.EIKON.EQ.=',2E15.5)
  102   FORMAT(1X,'P_I*Q_{IJ}=',2E15.5)
      END IF
C
C  CHECK THE POSITION OF THE RAY WITH RESPECT TO BOUNDARIES OF THE MODEL
C
      NTR=21
      IF(N.EQ.1)GO TO 13
      IF(IREF.GT.1.AND.(N-1).EQ.KINT(IREF-1)) GOTO 13
      NTR=22
      IF(Y(1).LT.BRD(1).OR.Y(1).GT.BRD(2).OR.Y(2).LT.BRD(3).OR.Y(2).
     1GT.BRD(4)) GOTO 20
C
C  CHECK WETHER THE RAY CROSSED THE BOREHOLE
C
      IF(ICOD.EQ.0)GO TO 4
      IF(ICOD.GT.IREF)GO TO 4
      FXYZA=(AY(2,N-1)-XVSP)*XNRM+(AY(3,N-1)-YVSP)*YNRM
      FXYZB=(Y(1)-XVSP)*XNRM+(Y(2)-YVSP)*YNRM
      IF(FXYZA*FXYZB.GT.0.)GO TO 4
C
C  THE RAY CROSSED THE VERTICAL PLANE PERPENDICULAR TO THE LINE
C  SOURCE-BOREHOLE AND CONTAINING THE BOREHOLE
C  DETERMINATION OF THE POINT OF INTERSECTION
C
      XNEW=X
      KDIM=6
      IF(NDER.GT.1)KDIM=18
      DO 41 I=1,KDIM
        YNEW(I)=Y(I)
        DNEW(I)=DERY(I)
   41 CONTINUE
      XA=XOLD
      XB=XNEW
C
C  DETERMINATION OF THE POINT OF INTERSECTION WITH THE VERTICAL PLANE
C  PERPENDICULAR TO THE LINE SOURCE-BOREHOLE
C
      CALL ROOT(XA,FXYZA,XB,FXYZB,X,PRMT,-1)
      IF(IAC.GE.10)GO TO 153
      CALL APPROX(X,Y,DERY,KDIM)
      IND=43
    4 CONTINUE
C
C  CHECK WETHER THE RAY CROSSED AN INTERFACE
C
      NTR=23
      INTR=LAY+1
      CALL DISC(Y,DEP)
      ZPL=DEP(1)
      INTR=LAY
      CALL DISC(Y,DEP)
      ZPU=DEP(1)
      NTR=24
      IF(ZPL.LE.ZPU) THEN
C
C  TWO NEIGHBOURHOOD INTERFACES CROSS EACH OTHER
C
      WRITE(LOU,'(A)') ' MAUX, LAY, IND,IND1,X,Y, Z,ZU,ZL'
      WRITE(LOU,'(4I5,5F12.5,/)')MAUX,LAY,IND,IND1,Y(1),Y(2),Y(3),ZPU,
     1ZPL
      GOTO 150
      END IF
      NTR=25
      IS(3,IREF)=LAY
      IF(Y(3).LT.ZPU.OR.Y(3).GT.ZPL) GOTO 30
      NTR=27
C
C  THE RAY DID NOT CROSS AN INTERFACE
C
C
      IF(IND.EQ.1.OR.IND.EQ.2)GO TO 148
      IF(IND.EQ.30.OR.IND.EQ.31)GO TO 148
   13 CONTINUE
      AY(1,N)=X
      XOLD=X
      DO 15 I=1,6
        AY(I+1,N)=Y(I)
        YOLD(I)=Y(I)
        DOLD(I)=DERY(I)
   15 CONTINUE
      DRY(1,N)=DERY(1)
      DRY(2,N)=DERY(2)
      DRY(3,N)=DERY(3)
      IF(NDER.GT.1)THEN
C
C     DETERMINATION OF KMAH INDEX
C
        IF(N.GT.1)THEN
          SPR1=Y(8)*Y(12)-Y(9)*Y(11)
          SPR2=Y(9)*Y(10)-Y(7)*Y(12)
          SPR3=Y(7)*Y(11)-Y(8)*Y(10)
          SPR=SPR1*DERY(1)+SPR2*DERY(2)+SPR3*DERY(3)
          IF(SPR*SPROLD.LT.0..AND.(X-TSTOLD).GT.0.)KMAH=KMAH+1
          SPROLD=SPR
          TSTOLD=X
        END IF
        DO 14 I=1,12
          YOLD(i+6)=y(i+6)
          DOLD(I+6)=DERY(I+6)
   14   CONTINUE
      END IF
      NTR=99
      IF(IND.EQ.10)GO TO 25
      IF(IND.EQ.43)GO TO 25
C
C  WRITING IN AY FIELD
C
      CALL PARDIS (Y,1)
      RETURN
C
C  THE RAY CROSSED ONE OF THE VERTICAL BOUNDARIES OF THE MODEL
C
   20 IF(Y(1).LT.BRD(1)) IND=1
      IF(Y(1).GT.BRD(2)) IND=2
      IF(Y(2).LT.BRD(3)) IND=3
      IF(Y(2).GT.BRD(4)) IND=4
      IF(IND.EQ.3.OR.IND.EQ.4) GOTO 23
      AUX=(BRD(IND)-AY(2,N-1))/(Y(1)-AY(2,N-1))
      TR=X
      X=AY(1,N-1)+AUX*(X-AY(1,N-1))
      K=6
      IF(NDER.GT.1)K=18
      DO 21 I=1,K
   21 Y(I)=Y(I)+DERY(I)*(X-TR)
      Y(1)=BRD(IND)
      AY(1,N)=X
      DO 3 I=1,6
   3  AY(I+1,N)=Y(I)
      GO TO 4
  23  AUX=(BRD(IND)-AY(3,N-1))/(Y(2)-AY(3,N-1))
      TR=X
      X=AY(1,N-1)+AUX*(X-AY(1,N-1))
      K=6
      IF(NDER.GT.1)K=18
      DO 16 I=1,K
   16 Y(I)=Y(I)+DERY(I)*(X-TR)
      Y(2)=BRD(IND)
      AY(1,N)=X
      DO 17 I=1,6
  17  AY(I+1,N)=Y(I)
      IF(IND.EQ.3) IND=30
      IF(IND.EQ.4) IND=31
      GO TO 4
C
C  THE RAY CROSSED AN INTERFACE
C
   30 CONTINUE
      IF(IREF.LE.1.OR.KC.GT.0) GOTO 39
      IF(N-KINT(IREF-1).EQ.2) THEN
      IND=9
      IND1=LAY
      GO TO 25
      END IF
C
C  WHICH INTERFACE WAS CROSSED
C
   39 INTR=LAY
      IF(Y(3).GE.ZPL) INTR=LAY+1
      XNEW=X
      KDIM=6
      IF(NDER.GT.1)KDIM=18
      DO 40 I=1,KDIM
        YNEW(I)=Y(I)
        DNEW(I)=DERY(I)
   40 CONTINUE
      DO 35 I=1,3
        YA(I)=YOLD(I)
        YB(I)=YNEW(I)
   35 CONTINUE
      CALL DISC(YA,DEP)
      FXYZA=DEP(1)-YA(3)
      CALL DISC(YB,DEP)
      FXYZB=DEP(1)-YB(3)
C
C  DETERMINATION OF THE POINT OF INTERSECTION WITH THE INTERFACE
C
      XA=XOLD
      XB=XNEW
      CALL ROOT(XA,FXYZA,XB,FXYZB,XINT,PRMT,1)
      IF(IAC.GE.10)GO TO 153
      NTR=10
C
C     THE QUANTITIES OF RAY TRACING AND DYNAMIC RAY TRACING
C     AT THE POINT OF INCIDENCE WILL BE DETERMINED
C
      IND1=INTR
      CALL APPROX(XINT,Y,DERY,KDIM)
      AY(1,N)=XINT
      DO 63 I=1,6
        AY(I+1,N)=y(I)
   63 CONTINUE
      CALL PARDIS(Y,1)
      CALL FCT(X,Y,DERY)
      DO 64 I=1,3
        DRY(I,N)=DERY(I)
   64 CONTINUE
      IF(NDER.GT.1)THEN
        IF(N.GT.1)THEN
          SPR1=Y(8)*Y(12)-Y(9)*Y(11)
          SPR2=Y(9)*Y(10)-Y(7)*Y(12)
          SPR3=Y(7)*Y(11)-Y(8)*Y(10)
          SPR=SPR1*DERY(1)+SPR2*DERY(2)+SPR3*DERY(3)
          IF(SPR*SPROLD.LT.0..AND.(X-TSTOLD).GT.0.)KMAH=KMAH+1
          SPROLD=SPR
          TSTOLD=X
        END IF
        DO 62 I=1,3
          XDYN(I,1)=Y(I+6)
          XDYN(I,2)=Y(I+9)
          XDYN(I,3)=DERY(I)
          YDYN(I,1)=Y(I+12)
          YDYN(I,2)=Y(I+15)
          YDYN(I,3)=DERY(I+3)
          XK(I)=DERY(I)
          YK(I)=DERY(I+3)
   62   CONTINUE
      END IF
C
C  DETERMINATION OF THE UNIT NORMAL TO THE INTERFACE AT THE POINT OF
C  INCIDENCE
C
      CALL DISC (Y,DEP)
      DFX=DEP(2)
      DFY=DEP(3)
      ROO=1+DFX*DFX+DFY*DFY
      UN3=ROO**(-0.5)
      UN(3)=-UN3
      UN(2)=UN3*DFY
      UN(1)=UN3*DFX
C
C  SCALAR PRODUCT OF GROUP VELOCITY WITH UNIT NORMAL VECTOR
C
      PN=DERY(1)*UN(1)+DERY(2)*UN(2)+DERY(3)*UN(3)
C
C  ORIENTATION OF THE UNIT NORMAL SO THAT IT POINTS TO THE MEDIUM IN
C  WHICH THE ENERGY OF INCIDENT WAVE PROPAGATES
C
      IF(PN.LT.0) GOTO 120
      UN(1)=-UN(1)
      UN(2)=-UN(2)
      UN(3)=-UN(3)
      UN3=-UN3
  120 CONTINUE
      IF(NDER.GT.1)THEN
        XIINC=Y(4)*UN(1)+Y(5)*UN(2)+Y(6)*UN(3)
        UN33=UN3*UN3*UN3
        DFXX=DEP(4)
        DFXY=DEP(5)
        DFYY=DEP(6)
        TM(1)=0.
        TM(2)=0.
        DO 126 J=1,2
          AUX=0.
          DO 125 I=1,3
            TM(J)=TM(J)+UN(I)*XDYN(I,J)
            AUX=AUX+UN(I)*XK(I)
  125     CONTINUE
          TM(J)=-TM(J)/AUX
  126   CONTINUE
C
C     DETERMINATION OF DERIVATIVES OF THE UNIT NORMAL TO INTERFACE
C
        AUX1=1.+DFX*DFX
        AUX2=1.+DFY*DFY
        AUX3=DFX*DFY
        DN(1,1)=(DFXX*AUX2-AUX3*DFXY)*UN33
        DN(1,2)=(DFXY*AUX2-AUX3*DFYY)*UN33
        DN(1,3)=0.
        DN(2,1)=(DFXY*AUX1-AUX3*DFXX)*UN33
        DN(2,2)=(DFYY*AUX1-AUX3*DFXY)*UN33
        DN(2,3)=0.
        DN(3,1)=(DFX*DFXX+DFY*DFXY)*UN33
        DN(3,2)=(DFX*DFXY+DFY*DFYY)*UN33
        DN(3,3)=0.
        DO 129 L=1,3
          DO 128 M=1,2
            DNG(L,M)=0.
            DO 127 K=1,3
              DNG(L,M)=DNG(L,M)+DN(L,K)*(XDYN(K,M)+XK(K)*TM(M))
  127       CONTINUE
  128     CONTINUE
  129   CONTINUE
      END IF
      NTR=20
      IF(KRE.EQ.1) GOTO 24
      NTR=30
      IF(KRE.EQ.0) GOTO 130
C
C  MULTIPLE REFLECTED WAVE
C
      IF((IREF+1).GT.KRE.AND.INTR.EQ.INT1) GOTO 151
C
C  TERMINATION OF THE RAY AT AN INNER INTERFACE OR AT THE FREE SURFACE
C  OR AT THE BOTTOM OF THE MODEL
C
      IF((IREF+1).GT.KRE) IRR=IREF
      NTR=35
      IF((IREF+1).GT.KRE) GOTO 148
C
C  NCD : INDICATOR WHETHER A REFLECTION OR TRANSMISSION TAKES PLACE
C        WITH RESPECT TO THE CODE OF THE WAVE
C
      NCD=CODE(IREF+1,1)-CODE(IREF,1)
C
C  NCD1 : INDICATOR FOR THE TYPE OF THE WAVE IN THE NEW LAYER
C
      NCD1=CODE(IREF+1,2)-CODE(IREF,2)
      NTR=40
      IF(KC.GT.0.AND.IREF.EQ.1.AND.INTR.EQ.LAY) GOTO 151
      NTR=50
      IF(KC.LT.0.AND.IREF.EQ.1.AND.INTR.NE.LAY) GOTO 151
      NTR=60
      IF(IREF.EQ.1) GOTO 170
      NTR=70
C     IF(INTR.EQ.INT1.AND.NCD.NE.0.OR.NCD1.NE.0) GOTO 151
      IF(INTR.EQ.INT1.AND.NCD.NE.0)GOTO 151
      NTR=75
      IF(INTR.EQ.INT1.AND.NCD1.NE.0)GOTO 151
C
C  REFLECTION OR TRANSMISSION OF THE RAY
C
      NTR=80
      IF(INTR.NE.INT1.OR.NCD.NE.0.OR.NCD1.NE.0) GOTO 170
C
C  COMPOUND ELEMENT OF THE RAY
C
      IREFR=1
      KINT(IREF)=0
      IREF=IREF+1
      NCD=CODE(IREF+1,1)-CODE(IREF,1)
      NCD1=CODE(IREF+1,2)-CODE(IREF,2)
C
C  TERMINATION AT AN INNER INTERFACE
C
      IRR=IREF
      NTR=90
      IF((IREF+1).GT.KRE) GOTO 24
      NTR=100
      GOTO 170
C
C  REFRACTED WAVE
C
C  ORDINARY TERMINATION AT UPPER BOUNDARY
C
  130 NTR=110
      IF(INTR.EQ.LAY.AND.LAY.EQ.1) GOTO 148
      NCD=1
      NCD1=0
C
C  SPECIFICATION OF THE LAYER OF THE GENERATED WAVE
C
  170 INT1=INTR
      IF(IRHO.NE.0)DS(8,IREF)=RHO(LAY)
      IRR=IREF
      IREF=IREF+1
      NTR=120
      IF(NCD.EQ.0) GOTO 200
      NTR=130
      IF(INTR.EQ.LAY) GOTO 190
C
C  TRANSMISSION AT THE LOWER INTERFACE
C
      NTR=140
      IF(NCD.LT.0) GOTO 151
C
C  ORDINARY TERMINATION AT LOWER BOUNDARY
C
      NTR=150
      IF(KRE.LE.1.AND.INTR.EQ.NINT) GOTO 148
      NTR=160
      IF(INTR.EQ.NINT) GOTO 151
      LAY=LAY+1
      NTR=170
      GOTO 200
C
C  TRANSMISSION AT THE UPPER INTERFACE
C
  190 NTR=180
      IF(NCD.GT.0.AND.KRE.GT.1) GOTO 151
      NTR=190
      IF(KRE.LE.1.AND.LAY.EQ.1) GOTO 24
      NTR=200
      IF(LAY.EQ.1) GOTO 151
      LAY=LAY-1
C
C  TRANSFORMATION OF THE SLOWNESS VECTOR
C
  200 IF(INTR.EQ.NINT.AND.MDIM.NE.0) GOTO 154
      DO 210 I=1,3
      POLD(I)=Y(3+I)
  210 CONTINUE
      ITRANS=0
C
C  CHECK WHETHER A REFLECTION OR TRANSMISSION TAKES PLACE
C
      IF(KC.EQ.0) GOTO 217
      IF(CODE(IREF,1).EQ.CODE(IREF-1,1)) GOTO 218
  217 ITRANS=1
  218 CALL TRANSL (Y,POLD,PNEW,UN,ITRANS,1)
      IF(IND.EQ.9.OR.IND.EQ.10) GOTO 25
      DO 220 I=1,3
      Y(I+3)=PNEW(I)
  220 CONTINUE
      IF(NDER.GT.1)THEN
        CALL FCT(X,Y,DERY)
        DO 225 I=1,3
          XKG(I)=DERY(I)
          YKG(I)=DERY(I+3)
  225   CONTINUE
        XI=PNEW(1)*UN(1)+PNEW(2)*UN(2)+PNEW(3)*UN(3)
        DXTN=XKG(1)*UN(1)+XKG(2)*UN(2)+XKG(3)*UN(3)
        DO 228 M=1,2
          AUX=XKG(1)*DNG(1,M)+XKG(2)*DNG(2,M)+XKG(3)*DNG(3,M)
          DO 224 I=1,3
            XDAUX(I)=XDYN(I,M)
            YDAUX(I)=YDYN(I,M)
  224     CONTINUE
          DO 227 I=1,3
            AUXX=(DNG(I,M)-UN(I)*AUX/DXTN)*(XI-XIINC)
            DO 226 K=1,3
              AUXX=AUXX-UN(I)*(XKG(K)*(YDaux(K)+YK(K)*TM(M))-
     1                  YKG(K)*(XDaux(K)+XK(K)*TM(M)))/DXTN
  226       CONTINUE
            XDYN(I,M)=XDYN(I,M)+(XK(I)-XKG(I))*TM(M)
            YDYN(I,M)=YDYN(I,M)+(YK(I)-YKG(I))*TM(M)+AUXX
  227     CONTINUE
  228   CONTINUE
        DO 229 I=1,3
          Y(I+6)=XDYN(I,1)
          Y(I+9)=XDYN(I,2)
          Y(I+12)=YDYN(I,1)
          Y(I+15)=YDYN(I,2)
  229   CONTINUE
      END IF
      PRMT(5)=2.
      RETURN
C
C  ORDINARY TERMINATION OF THE RAY
C
 148  IF(IND.EQ.1.OR.IND.EQ.2.OR.IND.EQ.30.OR.IND.EQ.31) GOTO 25
  24  IND=INTR+100
      IF((MREG.EQ.0.OR.MREG.GE.1).AND.MDIM.NE.0) THEN
        POLD(1)=Y(4)
        POLD(2)=Y(5)
        POLD(3)=Y(6)
        IS(3,IREF)=LAY
        IREF=IREF+1
        CALL TRANSL(Y,POLD,PNEW,UN,0,1)
        IREF=IREF-1
      END IF
      IF(KRE.LE.1) IRR=IREF
      IF(INTR.EQ.1) IND=3
      IF(INTR.EQ.NINT) IND=4
      KINT(IRR)=N
      IF(IND.NE.3) GOTO 25
C
C  COMPUTATION OF REFLECTED WAVES AT THE EARTH SURFACE FOR COEFFICIENTS
C  OF CONVERSION
C
      GOTO 25
  150 IND=20
      GOTO 25
  151 IND=8
      GOTO 25
  152 IND=13
      GOTO 25
  153 IND=39
      GOTO 25
  154 IND=15
      GOTO 25
  155 IND=12
   25 PRMT(5)=3.0
      IF(IND.NE.3.AND.IND.NE.43) RETURN
      ITYPE=CODE(IREF,2)
      CALL PARDIS(Y,2)
      RETURN
      END
C
C     *********************************************************
C
      SUBROUTINE POLAR (N1,N2,NN,I)
C
C   ROUTINE COMPUTING POLARIZATION VECTORS
C
C   IF(IPOL.EQ.0) POLARIZATION VECTORS ARE NOT PRINTED
C   IF(IPOL.EQ.1) POLARIZATION VECTORS ARE PRINTED ONLY
C   FOR THE POINTS OF INCIDENCE OF A RAY AT INTERFACES
C   IF(IPOL.EQ.2) POLARIZATION VECTORS ARE PRINTED FOR ALL
C   COMPUTED POINTS ALONG THE RAY
C   N1,N2 - SUCCESSIVE NUMBERS OF THE FIRST AND THE LAST POINT
C   OF AN ELEMENT OF A RAY (IN CASE IPOL.EQ.0, N1=N2 AND THIS
C   NUMBER CORRESPONDS TO THE POINT OF INCIDENCE OF A RAY AT
C   AN INTERFACE; NEGATIVE N2 - INDICATION OF COMPOSED ELEMENT
C   OF THE RAY
C   NN - TOTAL NUMBER OF POINTS ALONG THE RAY
C   I - NUMBER OF THE ELEMENT OF THE RAY
C
C   CALLED FROM ROUTINES: AMPL
C   ROUTINES CALLED:      CHRM
C
      DIMENSION PXN(3,3),Y(6),B1(3),B2(3),XX(3),E(21),E1(3),E2(3),T(3)
      COMMON /APROX/ A11,A12,A13,A14,A15,A16,A22,A23,A24,A25,A26,A33,
     1               A34,A35,A36,A44,A45,A46,A55,A56,A66,
     1               DXA11,DXA12,DXA13,DXA14,DXA15,DXA16,DXA22,DXA23,
     1               DXA24,DXA25,DXA26,DXA33,DXA34,DXA35,DXA36,DXA44,
     1               DXA45,DXA46,DXA55,DXA56,DXA66,
     1               DYA11,DYA12,DYA13,DYA14,DYA15,DYA16,DYA22,DYA23,
     1               DYA24,DYA25,DYA26,DYA33,DYA34,DYA35,DYA36,DYA44,
     1               DYA45,DYA46,DYA55,DYA56,DYA66,
     1               DZA11,DZA12,DZA13,DZA14,DZA15,DZA16,DZA22,DZA23,
     1               DZA24,DZA25,DZA26,DZA33,DZA34,DZA35,DZA36,DZA44,
     1               DZA45,DZA46,DZA55,DZA56,DZA66,
     1               A2546,A1266,A1355,A1456,A3645,A2344
      COMMON /AUXI/  IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT,
     1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOU,
     2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,MORI
      INTEGER CODE
      COMMON /COD/  CODE(50,2),KREF,KC,ITYPE
      COMMON /DJK/   D11,D12,D13,D22,D23,D33,DTR
      COMMON /GAM/   C11,C12,C13,C22,C23,C33
      COMPLEX PS
      COMMON /RAY/   AY(28,2000),DS(20,50),KINT(50),HHH(3,3),TMAX,
     1               PS(3,7,50),IS(8,50),N,IREF,IND,IND1
      COMMON /ZERO/ RNULL
C
      EQUIVALENCE(E(1),A11),(E(2),A12),(E(3),A13),(E(4),A14),(E(5),A15)
     1           ,(E(6),A16),(E(7),A22),(E(8),A23),(E(9),A24)
     2     ,(E(10),A25),(E(11),A26),(E(12),A33),(E(13),A34),(E(14),A35)
     2     ,(E(15),A36),(E(16),A44),(E(17),A45),(E(18),A46),(E(19),A55)
     2     ,(E(20),A56),(E(21),A66)
C
      IF(N2.LT.0)THEN
        LAY=IS(3,I-1)
        N2=-N2
      ELSE
        LAY=IS(3,I)
      END IF
      IF(IPOL.GT.0)
     1WRITE(LOU,'(A,6I5)')' LAY,I,N1,N2,ITYPE,IANI',LAY,I,N1,N2,ITYPE,
     2IANI(LAY)
      IF(IANI(LAY).NE.0) GO TO 2000
C
C  ISOTROPIC CASE
C
      IF(NN.EQ.N2) GO TO 1005
      IF(ITYPE.NE.3) GO TO 1020
C
C  P-WAVE POLARIZATION VECTOR
C
 1005 DO 1010 J=N1,N2
      V=1./SQRT(AY(5,J)*AY(5,J)+AY(6,J)*AY(6,J)+AY(7,J)*AY(7,J))
      T(1)=V*AY(5,J)
      T(2)=V*AY(6,J)
      T(3)=V*AY(7,J)
      IF((IPOL.EQ.1.AND.J.EQ.N2).OR.IPOL.EQ.2)
     1WRITE(LOU,'(A,/,6F10.5)') ' X1,X2,X3,TX,TY,TZ',
     2AY(2,J),AY(3,J),AY(4,J),T
 1010 CONTINUE
      DO 1015 K=1,3
      HHH(1,K)=0.
      HHH(2,K)=0.
      HHH(3,K)=T(K)
 1015 CONTINUE
      IF(N2.EQ.NN) GO TO 1020
      RETURN
C
C  S-WAVE POLARIZATION VECTORS
C
 1020 CONTINUE
      DO 1100 J=N1,N2
      V=SQRT(AY(12,J))
      V2=2.*V
      M1=5
      M2=6
      M3=7
      IF(ABS(V2).GT..000001)THEN
        DVS1=AY(13,J)/V2
        DVS2=AY(14,J)/V2
      ELSE
        DVS1=0.
        DVS2=0.
      END IF
 1030 P1=AY(M1,J)
      P2=AY(M2,J)
      P3=AY(M3,J)
      XL=SQRT(P1*P1+P2*P2)
      IF(XL.LT.RNULL)THEN
        IF(M1.EQ.4)THEN
          M1=6
          M2=7
          M3=5
          IF(ABS(V2).GT..000001)THEN
            DVS1=AY(14,J)/V2
            DVS2=AY(15,J)/V2
          ELSE
            DVS1=0.
            DVS2=0.
          END IF
          GO TO 1030
        END IF
        IF(M1.EQ.5)THEN
          M1=7
          M2=5
          M3=6
          IF(ABS(V2).GT..000001)THEN
            DVS1=AY(15,J)/V2
            DVS2=AY(13,J)/V2
          ELSE
            DVS1=0.
            DVS2=0.
          END IF
          GO TO 1030
        END IF
      END IF
C
C  COMPUTATION OF ANGLE XI SPECIFYING POLARIZATION VECTORS
C
      IF(IPOL.EQ.2)WRITE(LOU,'(A,/,6F12.6)')
     1' P1,P2,P3,XL,D1,D2',P1,P2,P3,XL,DVS1,DVS2
      IF(J.EQ.N1)THEN
        X1=0.
        TIME=0.
        XO=(P3*(DVS1*P2-DVS2*P1))/(XL*XL)
        XI=0.
        GO TO 1040
      END IF
      TIME=AY(1,J)-AY(1,J-1)
      X1=(P3*(DVS1*P2-DVS2*P1))/(XL*XL)
      XI=.5*(X1+XO)*TIME+XI
 1040 IF(IPOL.EQ.2)WRITE(LOU,'(A,/,6F12.6)')
     1' XI0,XI1,TSTEP,XISTEP',XO,X1,TIME,XI
      XO=X1
      P1=AY(5,J)
      P2=AY(6,J)
      P3=AY(7,J)
      XL=SQRT(P1*P1+P2*P2)
      IF(XL.LT..0001)THEN
        B1(1)=1.
        B1(2)=0.
        B1(3)=0.
        B2(1)=0.
        B2(2)=1.
        IF(P3.LT.0.)B2(2)=-1.
        B2(3)=0.
      ELSE
        U=SQRT(P1*P1+P2*P2+P3*P3)
        P13=P1*P3/XL
        P23=P2*P3/XL
        B1(1)=P13/U
        B1(2)=P23/U
        B1(3)=-xl/U
        B2(1)=-P2/XL
        B2(2)=P1/XL
        B2(3)=0.
      END IF
      SN=SIN(XI)
      CS=COS(XI)
      DO 1050 K=1,3
      E1(K)=CS*B1(K)-SN*B2(K)
      E2(K)=SN*B1(K)+CS*B2(K)
 1050 CONTINUE
      AY(17,J)=E1(1)
      AY(18,J)=E1(2)
      AY(19,J)=E1(3)
      AY(20,J)=E2(1)
      AY(21,J)=E2(2)
      AY(22,J)=E2(3)
      IF((IPOL.EQ.1.AND.J.EQ.N2).OR.IPOL.EQ.2)
     1WRITE(LOU,'(A,/,4F10.5/6F10.5)')
     2' TIME,X1,X2,X3/  E1X,E1Y,E1Z,E2X,E2Y,E2Z',(AY(K,J),K=1,4),E1,E2
 1100 CONTINUE
      DO 1110 K=1,3
      HHH(1,K)=E1(K)
      HHH(2,K)=E2(K)
      IF(N2.EQ.NN) GO TO 1110
      HHH(3,K)=0.
 1110 CONTINUE
      RETURN
C
C  ANISOTROPIC CASE
C
C
C  COMPUTATION OF POLARIZATION VECTORS
C
 2000 DO 2300 J=N1,N2
      DO 2110 K=1,21
      E(K)=AY(K+7,J)
 2110 CONTINUE
      DO 2120 K=1,6
      Y(K)=AY(K+1,J)
 2120 CONTINUE
      A2546=A25+A46
      A1266=A12+A66
      A1355=A13+A55
      A1456=A14+A56
      A3645=A36+A45
      A2344=A23+A44
      CALL CHRM(Y)
 2125 IF(ABS(DTR).LT..000001)GO TO 2135
      PXN(1,1)=D11
      PXN(1,2)=D12
      PXN(1,3)=D13
      XX(1)=D11*D11+D12*D12+D13*D13
      PXN(2,1)=D12
      PXN(2,2)=D22
      PXN(2,3)=D23
      XX(2)=D12*D12+D22*D22+D23*D23
      PXN(3,1)=D13
      PXN(3,2)=D23
      PXN(3,3)=D33
      XX(3)=D13*D13+D23*D23+D33*D33
      XN=0.
      DO 2130 L=1,3
      IF(XN.GE.XX(L))GO TO 2130
      XN=XX(L)
      NX=L
 2130 CONTINUE
      XN=SQRT(XN)
      IF(XN.GT.RNULL)GO TO 2140
 2135 CONTINUE
      WRITE(LOU,'(A,A,5I5)')
     1' XN.LT.RNULL IN POLAR  -  SHEAR-WAVE SINGULARITY ???',
     2'  LAY,N,N1,N2',LAY,J,J,N1,N2
      IND=10
      RETURN
 2140 P1=PXN(NX,1)/XN
      P2=PXN(NX,2)/XN
      P3=PXN(NX,3)/XN
      IF(J.NE.N1)THEN
        XN=P1OLD*P1+P2OLD*P2+P3OLD*P3
        IF(XN.LT.0.)THEN
          P1=-P1
          P2=-P2
          P3=-P3
        ENDIF
      ENDIF
      P1OLD=P1
      P2OLD=P2
      P3OLD=P3
C
C  CHECK OF PRECISION
C
      Z01=(C11-1.)*P1+C12*P2+C13*P3
      Z02=C12*P1+(C22-1.)*P2+C23*P3
      Z03=C13*P1+C23*P2+(C33-1.)*P3
      Z01=ABS(Z01)
      Z02=ABS(Z02)
      Z03=ABS(Z03)
      IF(Z01.GT.0.01.OR.Z02.GT.0.01.OR.Z03.GT.0.01) THEN
        WRITE(LOU,'(A,4I5)') ' CHRISTOFFEL EQUATION IS SATISFIED
     1  WITH PRECISION LESS THAN 0.01: LAY,ITYPE,RAY ELEMENT,NPOINT ',
     2  LAY,ITYPE,I,J
      END IF
      IF((IPOL.EQ.1.AND.J.EQ.N2).OR.IPOL.EQ.2)
     1WRITE(LOU,'(A,/,6F10.5)')
     2' X1,X2,X3,GX,GY,GZ',AY(2,J),AY(3,J),AY(4,J),P1,P2,P3
      IF(J.EQ.N2)THEN
        HHH(ITYPE,1)=P1
        HHH(ITYPE,2)=P2
        HHH(ITYPE,3)=P3
      END IF
 2300 CONTINUE
      RETURN
      END
C
C     *********************************************************
C
      SUBROUTINE POLRT(XCOF,COF,M,ZERO,IER)
      DIMENSION XCOF(7),COF(7)
      COMPLEX ZERO(6)
      DOUBLE PRECISION XO,YO,X,Y,XPR,YPR,UX,UY,V,YT,XT,U,XT2,YT2,SUMSQ,
     /                 DX,DY,TEMP,ALPHA,XCOF,COF,DABS
C
      IFIT=0
      N=M
      IER=0
      IF(XCOF(N+1)) 10,25,10
   10 IF(N) 15,15,32
C
   15 IER=1
   20 RETURN
C
   25 IER=4
      GO TO 20
C
   30 IER=2
      GO TO 20
   32 IF(N-36) 35,35,30
   35 NX=N
      NXX=N+1
      N2=1
      KJ1=N+1
      DO 40 L=1,KJ1
      MT=KJ1-L+1
   40 COF(MT)=XCOF(L)
C
   45 XO=-.00500101D0
      YO=0.01000101D0
C
      IN=0
   50 X=XO
C
      XO=-10.D0*YO
      YO=-10.D0*X
C
      X=XO
      Y=YO
      IN=IN+1
      GO TO 59
   55 IFIT=1
      XPR=X
      YPR=Y
C
   59 ICT=0
   60 UX=0.D0
      UY=0.D0
      V =0.D0
      YT=0.D0
      XT=1.D0
      U=COF(N+1)
      IF(U) 65,130,65
   65 DO 70 I=1,N
      L=N-I+1
      TEMP=COF(L)
      XT2=X*XT-Y*YT
      YT2=X*YT+Y*XT
      U=U+TEMP*XT2
      V=V+TEMP*YT2
      FI=I
      UX=UX+FI*XT*TEMP
      UY=UY-FI*YT*TEMP
      XT=XT2
   70 YT=YT2
      SUMSQ=UX*UX+UY*UY
      IF(SUMSQ) 75,110,75
   75 DX=(V*UY-U*UX)/SUMSQ
      X=X+DX
      DY=-(U*UY+V*UX)/SUMSQ
      Y=Y+DY
   78 IF(DABS(DY)+DABS(DX)-1.0D-12) 100,80,80
C
   80 ICT=ICT+1
      IF(ICT-500) 60,85,85
   85 IF(IFIT) 100,90,100
   90 IF(IN-5) 50,95,50
C
   95 IER=3
      GO TO 20
  100 DO 105 L=1,NXX
      MT=KJ1-L+1
      TEMP=XCOF(MT)
      XCOF(MT)=COF(L)
  105 COF(L)=TEMP
      ITEMP=N
      N=NX
      NX=ITEMP
      IF(IFIT) 120,55,120
  110 IF(IFIT) 115,50,115
  115 X=XPR
      Y=YPR
  120 IFIT=0
      IF(DABS(Y)-1.0D-10*DABS(X))135,125,125
  125 ALPHA=X+X
      SUMSQ=X*X+Y*Y
      N=N-2
      GO TO 140
  130 X=0.D0
      NX=NX-1
      NXX=NXX-1
  135 Y=0.D0
      SUMSQ=0.0
      ALPHA=X
      N=N-1
  140 COF(2)=COF(2)+ALPHA*COF(1)
  145 DO 150 L=2,N
  150 COF(L+1)=COF(L+1)+ALPHA*COF(L)-SUMSQ*COF(L-1)
  155 ZERO(N2)=CMPLX(X,Y)
      N2=N2+1
      IF(SUMSQ) 160,165,160
  160 Y=-Y
      SUMSQ=0.D0
      GO TO 155
  165 IF(N) 20,20,45
      END
C