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)
      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.000001)then
        th=pi
      end if
      if(a1.gt.0.000001)then
        th=q*q*q
        th=p/sqrt(th)
        th=acos(th)
      end if
      a1=-2.*sqrt(q)
      a2=-r/3.
      a3=th/3.
      v(1)=a1*cos(a3)+a2
      a3=(th+2.*pi)/3.
      v(2)=a1*cos(a3)+a2
      a3=(th+4.*pi)/3.
      v(3)=a1*cos(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(4),NINT,
     1    XINTA
      COMPLEX PS
      COMMON /RAY/   AY(28,400),DS(20,20),KINT(20),HHH(3,3),tmax,
     1               PS(3,7,20),IS(8,20),N,IREF,IND,IND1
      COMMON /RAY2/  DRY(3,400)
      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.400) 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(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
      IF(MREG.EQ.0.AND.MDIM.NE.0) THEN
        POLD(1)=Y(4)
        POLD(2)=Y(5)
        POLD(3)=Y(6)
        IREF=IREF+1
        CALL TRANSL(Y,POLD,PNEW,UN,0,1)
        IREF=IREF-1
      END IF
      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,400),DS(20,20),KINT(20),HHH(3,3),tmax,
     1               PS(3,7,20),IS(8,20),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)
      ns=n2-n1
      if(ns.eq.0)ns=1
      if(ipol.eq.2)ns=1
      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,ns
      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.ne.0)write(lou,'(a,/,6f10.5)') ' x1,x2,x3,tx,ty,tz',
     1ay(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
      DVS1=AY(13,J)/v2
      DVS2=AY(14,J)/v2
 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
          DVS1=AY(14,J)/v2
          DVS2=AY(15,J)/v2
          GO TO 1030
        END IF
        IF(M1.EQ.5)THEN
          M1=7
          M2=5
          M3=6
          DVS1=AY(15,J)/v2
          DVS2=AY(13,J)/v2
          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
      IF(IPOL.EQ.0.AND.J.NE.N2) GO TO 1100
      if(ipol.eq.1.and.(j.ne.n1.and.j.ne.n2))go to 1100
      P1=AY(5,J)
      P2=AY(6,J)
      P3=AY(7,J)
      XL=SQRT(P1*P1+P2*P2)
      if(xl.lt..000001)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
      if(ipol.ne.0)write(lou,'(a,/,4f10.5/6f10.5)')
     1' 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,ns
      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
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.ne.0)write(lou,'(a,/,6f10.5)') ' x1,x2,x3,gx,gy,gz',
     1ay(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
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=.00500101
      YO=0.01000101
C
      IN=0
   50 X=XO
C
      XO=-10.0*YO
      YO=-10.0*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.0
      UY=0.0
      V =0.0
      YT=0.0
      XT=1.0
      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-05) 100,80,80
C
   80 ICT=ICT+1
      IF(ICT-500) 60,85,85
   85 IF(IFIT) 100,90,100
   90 IF(N-5) 50,95,95
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(x)-1.0D-5)121,121,122
  121 x=0.
  122 IF(DABS(Y)-1.0D-3*DABS(X))135,125,125
  125 ALPHA=X+X
      SUMSQ=X*X+Y*Y
      N=N-2
      GO TO 140
  130 X=0.0
      NX=NX-1
      NXX=NXX-1
  135 Y=0.0
      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.0
      GO TO 155
  165 IF(N) 20,20,45
      END
C