C 
C
C     *********************************************************
C
      SUBROUTINE SPLIN(X,FX,NMIN,NMAX)
C
      DIMENSION A(200),B(200),H(200),F(200),X(200),FX(200)
C
C     CUBIC SPLINE INTERPOLATION WITH ZERO CURVATURES AT
C     THE END POINTS
C
      IF((NMAX-NMIN).EQ.1)GO TO 4
      NMIN1=NMIN+1
      NMAX1=NMAX-1
      H(NMIN1)=X(NMIN1)-X(NMIN)
      D2=(FX(NMIN1)-FX(NMIN))/H(NMIN1)
      DO 1 I=NMIN1,NMAX1
      H(I+1)=X(I+1)-X(I)
      D1=D2
      D2=(FX(I+1)-FX(I))/H(I+1)
      B(I)=H(I)+H(I+1)
      FX(I)=3.*(D2-D1)/B(I)
      A(I)=H(I)/B(I)
    1 B(I)=H(I+1)/B(I)
    4 FX(NMIN)=0.
      FX(NMAX)=0.
      IF((NMAX-NMIN).EQ.1)RETURN
      H(NMIN)=0.
      F(NMIN)=0.
      DO 2 I=NMIN1,NMAX1
      XPOM=2.+A(I)*H(I-1)
      H(I)=-B(I)/XPOM
    2 F(I)=(FX(I)-A(I)*F(I-1))/XPOM
      DO 3 I=NMIN,NMAX1
      J=NMAX1-(I-NMIN)
    3 FX(J)=H(J)*FX(J+1)+F(J)
      RETURN
      END
C
C     *********************************************************
C
      SUBROUTINE STRESS(P,XN,U,TAU)
C
C  ROUTINE FOR THE COMPUTATION OF NORMAL STRESSES
C
      IMPLICIT COMPLEX (P,U,T,E)
      REAL N1,N2,N3
      DIMENSION XN(3),P(3),U(3),TAU(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,LOUT,
     2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori
      common/dens/rho(20)
C
      P1=P(1)
      P2=P(2)
      P3=P(3)
      U1=U(1)
      U2=U(2)
      U3=U(3)
      N1=XN(1)
      N2=XN(2)
      N3=XN(3)
      IF(IANI(LAY).EQ.0) GO TO 1000
C
C  ANISOTROPIC CASE
C
      RO=1.7+0.2*SQRT(A11)
      IF(IRHO.NE.0) RO=rho(lay)
C
C  DISPLACEMENT TENSOR
C
      E11=U1*P1
      E12=0.5*(U1*P2+U2*P1)
      E13=0.5*(U1*P3+U3*P1)
      E22=U2*P2
      E23=0.5*(U2*P3+U3*P2)
      E33=U3*P3
C
C  STRESS TENSOR
C
      P11=A11*E11+A12*E22+A13*E33+2.*A14*E23+2.*A15*E13+2.*A16*E12
      P22=A12*E11+A22*E22+A23*E33+2.*A24*E23+2.*A25*E13+2.*A26*E12
      P33=A13*E11+A23*E22+A33*E33+2.*A34*E23+2.*A35*E13+2.*A36*E12
      P23=A14*E11+A24*E22+A34*E33+2.*A44*E23+2.*A45*E13+2.*A46*E12
      P13=A15*E11+A25*E22+A35*E33+2.*A45*E23+2.*A55*E13+2.*A56*E12
      P12=A16*E11+A26*E22+A36*E33+2.*A46*E23+2.*A56*E13+2.*A66*E12
C
C  NORMAL STRESS
C
      TAU(1)=P11*N1+P12*N2+P13*N3
      TAU(2)=P12*N1+P22*N2+P23*N3
      TAU(3)=P13*N1+P23*N2+P33*N3
      TAU(1)=TAU(1)*RO
      TAU(2)=TAU(2)*RO
      TAU(3)=TAU(3)*RO
      RETURN
C
C  ISOTROPIC CASE
C
 1000 A12=A11-2.*A44
      RO=1.7+0.2*SQRT(A11)
      IF(IRHO.NE.0) RO=rho(lay)
      IF(ICOEF.EQ.0) GOTO 100
      WRITE(LOUT,'(A,4F12.7)') 'STRESS: A11,A44,A12,RO',A11,A44,A12,RO
      WRITE(LOUT,'(A,/,3(2F12.7,/))') 'STRESS: PSTR',P
      WRITE(LOUT,'(A,/,3(2F12.7,/))') 'STRESS: USTR',U
  100 TETA=U1*P1+U2*P2+U3*P3
      IF(ICOEF.GT.0)
     1WRITE(LOUT,'(A,/,5F12.7)') 'STRESS: UN,TETA',XN,TETA
      P11=A12*TETA+2.*A44*P1*U1
      P12=A44*(P1*U2+P2*U1)
      P13=A44*(P1*U3+P3*U1)
      P22=A12*TETA+2.*A44*P2*U2
      P23=A44*(P2*U3+P3*U2)
      P33=A12*TETA+2.*A44*P3*U3
      IF(ICOEF.GT.0)
     1WRITE(LOUT,'(A,/,6(2F12.7,/))')'STRESS: P11,P12,P13,P22,P23,P33'
     1,P11,P12,P13,P22,P23,P33
      TAU(1)=P11*N1+P12*N2+P13*N3
      TAU(2)=P12*N1+P22*N2+P23*N3
      TAU(3)=P13*N1+P23*N2+P33*N3
      TAU(1)=TAU(1)*RO
      TAU(2)=TAU(2)*RO
      TAU(3)=TAU(3)*RO
      RETURN
      END
C
C     *********************************************************
C
      subroutine surfpl(lin,lu3)
c
      dimension xx(300),yy(300),zz(300)
      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 /DIST/ XDST(200),NDSTX,XREPS,DST(2),NDST,REPS,LNDST,
     1xprf,yprf
      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)
c
      dst(1)=0.
      mori=1
      ndst=1
      reps=.001
      tmax=.01
      ind=1
      isour=1
      read(lin,100)npar
      read(lin,102)x0,y0,z0,ddelta
      if(abs(x0).lt..0001.or.abs(y0).lt..0001.or.abs(z0).lt..0001)then
        x0=1.
        y0=1.
        z0=1.
      end if
      if(abs(ddelta).lt..000001)ddelta=.05
      write(lou,100)npar
      write(lou,102)x0,y0,z0,ddelta
      xprf=x0
      yprf=y0
      kref=1
      itype=0
      code(1,1)=1
    1 continue
      inum=0
      delta=0.
      itype=itype+1
      code(1,2)=itype
      if(npar.le.2)ndstx=0
      if(npar.eq.3)then
        ndstx=1
        nder=2
        mdim=2
      end if
    2 continue
      inum=inum+1
      if(npar.le.2)then
        call PROFIL(x0,y0,z0,0.,delta,PAZM,RANG,
     1  X,Y,Z,T,.2,.0001,0.,.6,.4,10,3,1,0,12,0)
        xx(inum)=ay(5,1)
        yy(inum)=ay(6,1)
        zz(inum)=ay(7,1)
      end if
      if(npar.eq.3)then
        call PROFIL(x0,y0,z0,0.,delta,PAZM,RANG,
     1  X,Y,Z,T,.004,.0001,-.5,1.,.6,10,3,1,0,12,0)
        if(ind.ne.0)then
          xx(inum)=dry(1,1)
          yy(inum)=dry(2,1)
          zz(inum)=dry(3,1)
        end if
      end if
      delta=delta+ddelta
      if(inum.lt.300.and.delta.lt.6.3)go to 2
      write(lu3,100)itype,npar,inum
      delta=0.
      do 3 i=1,inum
      write(lu3,101)delta,xx(i),yy(i),zz(i)
      delta=delta+ddelta
    3 continue
      if(itype.lt.3)go to 1
      itype=0
      write(lu3,100)itype
c
  100 format(16i5)
  101 format(4e15.5)
  102 format(4f10.5)
      return
      end
C
C     *********************************************************
C
      SUBROUTINE TRANSL(Y,P,PNEW,UN,ITRANS,IASW)
C
C  ROUTINE FOR THE COMPUTATION OF THE TRANSFORMATION OF THE SLOWNESS
C  VECTOR AT AN INTERFACE
C
      SAVE A,AI,B,BI,C,CI,CD,CDI,XCOE,AT,BT,DETB
      SAVE XKO,PN,Y1,Z,NDER1,LAY1,ISG,IMAG
      DOUBLE PRECISION XCOE,COE
      DIMENSION A(3,3),AI(3,3),B(3,3),BI(3,3),C(3,3),CI(3,3),CD(3,3),
     *         CDI(3,3),AC(3,3),ACI(3,3),XH1(3),XH2(3),
     *         P(3),PNEW(3),UN(3),Y(18),PN(3),Y1(18),DERY(18),XCOE(7),
     *         COE(7),XSI(3),ISG(6),IR(3),IT(3),ICR(3),ICT(3),XKO(7)
      COMPLEX Z(6),Z0,CSQ,CSI(3),XPR,IMAG
      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,LOUT,
     2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori
      INTEGER CODE
      COMMON /COD/  CODE(50,2),KREF,KC,ITYPE
      COMMON /DENS/ RHO(20)
      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
C  IASW :      SWITCH INDICATING FROM WHICH PROGRAM TRANSL IS CALLED
C              IASW=0   CALLED FROM AMPL
C              IASW=1   CALLED FROM OUT
C
      IMAG=(0.,1.)
      IRF=IREF-1
C
C     FILLING DS FIELD
C
      IF(IASW.GT.0) THEN
        DS(1,IRF)=UN(1)
        DS(2,IRF)=UN(2)
        DS(3,IRF)=UN(3)
        IF(IRHO.EQ.0)DS(8,IRF)=0.2*SQRT(A11)+1.7
        KINT(IRF)=N
        IS(1,IRF)=LAY
        IS(2,IRF)=ITRANS
      END IF
      CALL PARDIS(Y,0)
      IF(ITRANS.EQ.0)IS(7,IRF)=LAY
      IF(ITRANS.NE.0)IS(8,IRF)=LAY
      DO 1 K=1,3
      IT(K)=0
      IR(K)=0
      ICT(K)=0
      ICR(K)=0
      ISG(K)=0
      ISG(K+3)=0
   1  CONTINUE
      KDIM=6
      IF(NDER.GT.1)KDIM=18
      DO 2 K=1,KDIM
   2  Y1(K)=Y(K)
      PUN=P(1)*UN(1)+P(2)*UN(2)+P(3)*UN(3)
      PU1=PUN*UN(1)
      PU2=PUN*UN(2)
      PU3=PUN*UN(3)
      PN(1)=P(1)-PU1
      PN(2)=P(2)-PU2
      PN(3)=P(3)-PU3
      IF(IASW.EQ.0) GO TO 9
      NDER1=NDER
      NDER=1
      lay1=lay
      lay=is(3,irf)
      CALL FCT(X,Y1,DERY)
      if(ind.eq.10)return
      DUN=DERY(1)*UN(1)+DERY(2)*UN(2)+DERY(3)*UN(3)
      lay=lay1
      CALL FCT(X,Y1,DERY)
      if(ind.eq.10)return
      NDER=NDER1
      DS(4,IRF)=PN(1)
      DS(5,IRF)=PN(2)
      DS(6,IRF)=PN(3)
      DS(7,IRF)=abs(dun)
   9  IF(KC.NE.0) ITYPE=CODE(IREF,2)
      IF(IREF.GT.KREF) ITYPE=CODE(KREF,2)
      IF(IANI(LAY).EQ.0) GOTO 230
C
C  TRANSMISSION OR REFLECTION FOR THE ANISOTROPIC CASE
C
      IF(MTRNS.GT.0) THEN
        WRITE(LOUT,999)UN,PN,(Y(K),K=1,KDIM),P
 999    FORMAT(' UN,PN,Y,P',/,6(E12.5,2X))
        WRITE(LOUT,'(A)')' ITYPE,IREF,LAY,IANI(LAY),N,KINT(IREF-1)'
        WRITE(LOUT,'(6I5)')ITYPE,IREF,LAY,IANI(LAY),N,KINT(IREF-1)
        WRITE(LOUT,'(A,/,10F10.5,/,11F10.5)')' A(I,J,K,L)=',A11,A12,A13,
     1  A14,A15,A16,A22,A23,A24,A25,A26,A33,A34,A35,A36,A44,A45,A46,A55,
     2  A56,A66
      END IF
      CALL CHRM1(A,UN,UN)
      AT=A(1,1)+A(2,2)+A(3,3)
      CALL INV(A,AI,DET)
      XCOE(7)=DET
      CALL CHRM1(C,PN,PN)
      CALL INV(C,CI,DET)
      DO 10 J=1,3
      DO 10 K=1,3
      CD(J,K)=C(J,K)
      IF(J.EQ.K)CD(J,K)=CD(J,K)-1
   10 CONTINUE
      CALL INV(CD,CDI,DET)
      XCOE(1)=DET
      CALL CHRM1(B,PN,UN)
      CALL INV(B,BI,DET)
      BT=B(1,1)+B(2,2)+B(3,3)
      DETB=DET
      DO 20 J=1,3
      DO 20 K=1,3
      AC(J,K)=A(J,K)+C(J,K)
  20  CONTINUE
      CALL INV(AC,ACI,DET)
      DO 22 J=1,3
      XH1(J)=0.0
      DO 21 M=1,3
  21  XH1(J)=AI(M,J)*B(J,M)+XH1(J)
  22  CONTINUE
      XCOE(6)=XH1(1)+XH1(2)+XH1(3)
      XCOE(6)=2.*XCOE(6)
      DO 30 J=1,3
      XH1(J)=0.0
      DO 29 M=1,3
  29  XH1(J)=AI(M,J)*CD(J,M)+XH1(J)
  30  CONTINUE
      DO 40 J=1,3
      XH2(J)=0.0
      DO 39 M=1,3
  39  XH2(J)=BI(M,J)*A(J,M)+XH2(J)
  40  CONTINUE
      XCOE(5)=0.0
      DO 60 K=1,3
      XCOE(5)=XH1(K)+4.*XH2(K)+XCOE(5)
  60  CONTINUE
      DO 70 K=1,3
      DO 70 J=1,3
      ACI(J,K)=ACI(J,K)-AI(J,K)-CI(J,K)
  70  CONTINUE
      DO 80 J=1,3
      XH1(J)=0.0
      DO 79 M=1,3
  79  XH1(J)=ACI(M,J)*B(J,M)+XH1(J)
  80  CONTINUE
      DO 90 J=1,3
      XH2(J)=0.0
      DO 89 M=1,3
  89  XH2(J)=A(M,J)*B(J,M)+XH2(J)
  90  CONTINUE
      XCOE(4)=0.0
      DO 100 J=1,3
      XCOE(4)=XH1(J)+XH2(J)+XCOE(4)
 100  CONTINUE
      XCOE(4)=4.*DETB+XCOE(4)-AT*BT
      XCOE(4)=2.*XCOE(4)
      DO 110 J=1,3
      XH1(J)=0.0
      DO 109 M=1,3
 109  XH1(J)=CDI(M,J)*A(J,M)+XH1(J)
 110  CONTINUE
      DO 120 J=1,3
      XH2(J)=0.0
      DO 119 M=1,3
 119  XH2(J)=BI(M,J)*CD(J,M)+XH2(J)
 120  CONTINUE
      XCOE(3)=0.0
      DO 130 J=1,3
      XCOE(3)=XH1(J)+4.*XH2(J)+XCOE(3)
 130  CONTINUE
      DO 140 J=1,3
      XH1(J)=0.0
      DO 139 M=1,3
 139  XH1(J)=CDI(M,J)*B(J,M)+XH1(J)
 140  CONTINUE
      XCOE(2)=XH1(1)+XH1(2)+XH1(3)
      XCOE(2)=2.*XCOE(2)
      DO 141 K=1,7
      M=8-K
  141 XKO(M)=XCOE(K)
      CALL POLRT(XCOE,COE,6,Z,IER)
      iww=0
  142 IF(MTRNS.GT.0)then
        WRITE(LOUT,1000)XKO
        WRITE(LOUT,1010) Z
 1000   FORMAT(' COEF. OF POLYNOMIAL',/,4(E12.5,4X))
 1010   FORMAT(' ROOTS OF POLYNOMIAL',/,4(E12.5,4X))
      end if
      IF(IER.EQ.131) THEN
        WRITE(LOUT,1030)
 1030   FORMAT(/' NOT ALL ZEROS FOUND , PROGRAM TERMINATES'//)
        STOP
      END IF
      IF(IER.EQ.130) THEN
        WRITE(LOUT,1040)
 1040   FORMAT(/' DEGREE OF POLYNOMAL IS REDUCED %%%%%%%%%%%%%%%%'//)
      END IF
      DO 150 J=1,6
      ZI=AIMAG(Z(J))
      IF(ABS(ZI).GT..01)THEN
        ISG(J)=0
        GOTO 165
      else
        z(j)=z(j)-zi*IMAG
      END IF
C
C CHECK OF REAL-VALUED SLOWNESS VECTORS OF GENERATED WAVES
C
      Z0=Z(J)
      xpr=xko(7)+z0*(xko(6)+z0*(xko(5)+z0*(xko(4)+z0*(xko(3)+z0*(xko(2)+
     *z0*xko(1))))))
      IF(MTRNS.ne.0.or.(iww.eq.1.and.cabs(xpr).gt..00001))then
        WRITE(LOUT,1041) XPR,Z0
 1041   FORMAT(' PRECISSION OF ROOTS',/,4(E12.5,4X))
      end if
      DO 160 K=1,3
      Y1(K+3)=PN(K)+REAL(Z(J))*UN(K)
  160 CONTINUE
      NDER1=NDER
      NDER=1
      CALL FCT(X,Y1,DERY)
      IF(IND.EQ.10)RETURN
      NDER=NDER1
      AAA=Y1(4)*DERY(1)+Y1(5)*DERY(2)+Y1(6)*DERY(3)
      IF(MTRNS.GT.0)WRITE(LOUT,2000)AAA,Y1(4),Y1(5),Y1(6)
 2000 FORMAT(1X,'GROUP*SLOWNESS, SLOWNESS'/1X,4E15.5)
      IF(ABS(AAA-1.).GT.1.0E-02)THEN
        IND=10
        RETURN
      END IF
      SG=UN(1)*DERY(1)+UN(2)*DERY(2)+UN(3)*DERY(3)
      IF(MTRNS.GT.0)WRITE(LOUT,1042)SG
 1042 FORMAT(' DIRECTION OF RAY VECTOR (UN(I)*DERY(I),I=1,3',F10.5)
      IF(MTRNS.GT.0) WRITE(LOUT,1043)(DERY(i),i=1,3)
 1043 FORMAT(' RAY VECTOR (DERY)'/3F10.5)
      IF(SG.GE.0.)ISG(J)=1
      IF(SG.LT.0.)ISG(J)=-1
  165 CONTINUE
  150 CONTINUE
      MCT=0
      MCR=0
C
C CHECK OF COMPLEX-VALUED SLOWNESS VECTORS OF GENERATED WAVES
C
      DO 179 K=1,6
      IF(ISG(K).NE.0) GOTO 179
      DIR=AIMAG(Z(K))*UN(3)
      IF(UN(3).GT.0.0) GOTO 175
C
C  INCIDENT WAVE STRIKES THE INTERFACE FROM ABOVE
C
      IF(DIR.LT.0.0)THEN
        MCR=MCR+1
        ICR(MCR)=K
      END IF
      IF(DIR.GT.0.0)THEN
        MCT=MCT+1
        ICT(MCT)=K
      END IF
      GOTO 179
C
C  INCIDENT WAVE STRIKES THE INTERFACE FROM BELOW
C
  175 IF(DIR.GT.0.0)THEN
        MCR=MCR+1
        ICR(MCR)=K
      END IF
      IF(DIR.LT.0.0)THEN
        MCT=MCT+1
        ICT(MCT)=K
      END IF
  179 CONTINUE
      MT=0
      DO 180 K=1,6
      IF(ISG(K).EQ.(-1))THEN
        MT=MT+1
        IT(MT)=K
      END IF
  180 CONTINUE
      MR=0
      DO 190 K=1,6
      IF(ISG(K).EQ.1)THEN
        MR=MR+1
        IR(MR)=K
      END IF
  190 CONTINUE
C
      IF(MTRNS.EQ.0) GOTO 191
      WRITE(LOUT,1021)ISG
 1021 FORMAT(' INDICATIONS,ISG(6),MCR ICR(3),MCT ICT(3),MR IR(3)
     1MT IT(3) ',/,6I5)
      WRITE(LOUT,1020)MCR,ICR
      WRITE(LOUT,1020)MCT,ICT
      WRITE(LOUT,1020)MR,IR
      WRITE(LOUT,1020)MT,IT
1020  FORMAT(14I5)
 191  IF(ITRANS.EQ.0)IS(5,IRF)=MCR
      IF(ITRANS.GT.0)IS(6,IRF)=MCT
      IF(ITRANS.EQ.0)GO TO 200
C
C     TRANSMISSION
C
      IF(MT.EQ.3) THEN
      Z1=REAL(Z(IT(1)))
      Z2=REAL(Z(IT(2)))
      Z3=REAL(Z(IT(3)))
      XSI(3)=AMAX1(Z1,Z2,Z3)
      XSI(1)=AMIN1(Z1,Z2,Z3)
      XSI(2)=Z1+Z2+Z3-XSI(3)-XSI(1)
      END IF
      IF(MT.EQ.2) THEN
      Z1=REAL(Z(IT(1)))
      Z2=REAL(Z(IT(2)))
      XSI(2)=AMAX1(Z1,Z2)
      XSI(1)=Z1+Z2-XSI(2)
      CSI(1)=Z(ICT(1))
      END IF
      IF(MT.EQ.1) THEN
      XSI(1)=REAL(Z(IT(1)))
      IF(CABS(Z(ICT(1))).GT.CABS(Z(ICT(2)))) THEN
      CSI(1)=Z(ICT(2))
      CSI(2)=Z(ICT(1))
      ELSE
      CSI(1)=Z(ICT(1))
      CSI(2)=Z(ICT(2))
      END IF
      END IF
      IF(MT.EQ.0) THEN
      CSI(1)=CMPLX(999.,999.)
      CSI(3)=CMPLX(0.,0.)
      DO 189 K=1,3
      IF(CABS(Z(ICT(K))).LT.CABS(CSI(1))) CSI(1)=Z(ICT(K))
      IF(CABS(Z(ICT(K))).GT.CABS(CSI(3))) CSI(3)=Z(ICT(K))
  189 CONTINUE
      CSI(2)=Z(ICT(1))+Z(ICT(2))+Z(ICT(3))-CSI(1)-CSI(3)
      END IF
      IF(MDIM.LT.1) GOTO 196
C
C     COMPUTATION OF TRANSMITTED SLOWNESS VECTORS FOR EVALUATING
C     AMPLITUDES IN ROUTINE AMPL
C
      DO 192 K=1,MT
      DO 193 L=1,3
      PS(L,K+3,IRF)=PN(L)+XSI(K)*UN(L)
  193 CONTINUE
  192 CONTINUE
      DO 194 K=1,MCT
      I=ICT(K)
      M=MT+K
      DO 195 L=1,3
      PS(L,M+3,IRF)=PN(L)+CSI(K)*UN(L)
  195 CONTINUE
  194 CONTINUE
      IF(IASW.EQ.0)RETURN
C
C     OVERCRITICAL INCIDENCE
C
  196 IF(ITYPE.GT.MT) IND=9
      IF(IND.EQ.9) RETURN
C
      GOTO 210
C
C     REFLECTION
C
  200 IF(MR.EQ.3) THEN
      Z1=REAL(Z(IR(1)))
      Z2=REAL(Z(IR(2)))
      Z3=REAL(Z(IR(3)))
      XSI(3)=AMIN1(Z1,Z2,Z3)
      XSI(1)=AMAX1(Z1,Z2,Z3)
      XSI(2)=Z1+Z2+Z3-XSI(3)-XSI(1)
      END IF
      IF(MR.EQ.2) THEN
      Z1=REAL(Z(IR(1)))
      Z2=REAL(Z(IR(2)))
      XSI(2)=AMIN1(Z1,Z2)
      XSI(1)=Z1+Z2-XSI(2)
      CSI(1)=Z(ICR(1))
      END IF
      IF(MR.EQ.1) THEN
      XSI(1)=REAL(Z(IR(1)))
      IF(CABS(Z(ICR(1))).GT.CABS(Z(ICR(2)))) THEN
      CSI(1)=Z(ICR(2))
      CSI(2)=Z(ICR(1))
      ELSE
      CSI(1)=Z(ICR(1))
      CSI(2)=Z(ICR(2))
      END IF
      END IF
      IF(MR.EQ.0) THEN
      CSI(1)=CMPLX(0.0)
      CSI(3)=CMPLX(999.,999.)
      DO 211 K=1,3
      IF(CABS(Z(ICR(K))).GT.CABS(CSI(1))) CSI(3)=Z(ICR(K))
      IF(CABS(Z(ICR(K))).LT.CABS(CSI(3))) CSI(1)=Z(ICR(K))
  211 CONTINUE
      CSI(2)=Z(ICR(1))+Z(ICR(2))+Z(ICR(3))-CSI(1)-CSI(3)
      END IF
C
C     COMPUTATION OF REFLECTED SLOWNESS VECTORS FOR EVALUATING
C     AMPLITUDES IN ROUTINE AMPL
C
      IF(MDIM.LT.1) GOTO 209
      DO 201 K=1,MR
      DO 202 L=1,3
      PS(L,K,IRF)=PN(L)+XSI(K)*UN(L)
  202 CONTINUE
  201 CONTINUE
      DO 203 K=1,MCR
      I=ICR(K)
      M=MR+K
      DO 204 L=1,3
      PS(L,M,IRF)=PN(L)+CSI(K)*UN(L)
  204 CONTINUE
  203 CONTINUE
      IF(IASW.EQ.0) RETURN
C
C     OVERCRITICAL INCIDENCE
C
 209  IF(ITYPE.GT.MR) IND=9
      IF(IND.EQ.9) RETURN
C
  210 IF(MTRNS.GT.0) WRITE(LOUT,1111) XSI,Z1,Z2,Z3,ITYPE
1111  FORMAT(' XSI,Z1,Z2,Z3,ITYPE',/,6E14.6,I10)
      XSSI=XSI(ITYPE)
      PU1=XSSI*UN(1)
      PU2=XSSI*UN(2)
      PU3=XSSI*UN(3)
      PNEW(1)=PN(1)+PU1
      PNEW(2)=PN(2)+PU2
      PNEW(3)=PN(3)+PU3
      Y(4)=PNEW(1)
      Y(5)=PNEW(2)
      Y(6)=PNEW(3)
      NDER1=NDER
      NDER=1
      CALL FCT(X,Y,DERY)
      if(ind.eq.10)return
      NDER=NDER1
      DUN=DERY(1)*UN(1)+DERY(2)*UN(2)+DERY(3)*UN(3)
      RHO2=0.2*SQRT(A11)+1.7
      IF(IRHO.NE.0) RHO2=RHO(LAY)
      DS(10,IRF)=abs(dun)
      DS(11,IRF)=RHO2
      IF(MTRNS.GT.0)WRITE(LOUT,1100)PNEW
 1100 FORMAT(' PNEW',/,3F12.8)
      RETURN
C
C  COMPUTATION OF THE TRANSFORMATION OF THE SLOWNESS VECTOR
C  IN THE ISOTROPIC CASE
C
  230 IF(ITYPE.EQ.3)CI2=1./A11
      IF(ITYPE.NE.3)CI2=1./A44
      PP=0.0
      DO 250 K=1,3
      PP=PP+P(K)*P(K)
  250 CONTINUE
      PUN2=PUN*PUN
      ROO=CI2-PP+PUN2
      ROT1=1./A11-PP+PUN2
      ROT2=1./A44-PP+PUN2
      IF(MTRNS.GT.0) WRITE(LOUT,1065) ROO,ROT1,ROT2
 1065 FORMAT('ROO,ROT1,ROT2',3F10.5)
C  ROT1.LT.ROT2 FOR PHYSICAL SUFFCIENT ELAST. COEF.
      IF(ITRANS.EQ.0.AND.MDIM.GE.1) THEN
      IF(ROT1.GT.0.AND.ROT2.GT.0) MR=3
      IF(ROT1.GT.0.AND.ROT2.GT.0) MCR=0
      IF(ROT1.LE.0.AND.ROT2.GT.0) MR=2
      IF(ROT1.LE.0.AND.ROT2.GT.0) MCR=1
      IF(ROT1.LE.0.AND.ROT2.LE.0) MCR=3
      IF(ROT1.LE.0.AND.ROT2.LE.0) MR=0
      IS(5,IRF)=MCR
      IF(MTRNS.GT.0) WRITE(LOUT,1066)MR,MCR
 1066 FORMAT(' MR,MCR',/,2I5)
      DO 251 K=1,MR
      IF(K+MCR.EQ.1) RO=ROT1
      IF(K+MCR.GT.1) RO=ROT2
      SQU=SQRT(RO)
      IF(MTRNS.GT.0) WRITE(LOUT,1062) SQU
 1062 FORMAT(' SQU',/,2F12.6)
      IS(5,IRF)=MCR
      J=MR+1-K
      DO 252 L=1,3
      PS(L,J,IRF)=P(L)-(PUN-SQU)*UN(L)
  252 CONTINUE
  251 CONTINUE
      DO 253 K=1,MCR
      M=MR+K
      IF(M.EQ.3) CSQ=CSQRT(CMPLX(ROT1,0.0))
      IF(M.NE.3) CSQ=CSQRT(CMPLX(ROT2,0.0))
      IF(MTRNS.GT.0) WRITE(LOUT,1064)CSQ
 1064 FORMAT(' CSQ',2F10.5)
      DO 254 L=1,3
      PS(L,M,IRF)=P(L)-(PUN-CSQ)*UN(L)
  254 CONTINUE
  253 CONTINUE
      IF(IASW.EQ.0) RETURN
      END IF
      IF(ITRANS.EQ.1.AND.MDIM.GE.1) THEN
      IF(ROT1.GT.0.AND.ROT2.GT.0) MT=3
      IF(ROT1.GT.0.AND.ROT2.GT.0) MCT=0
      IF(ROT1.LE.0.AND.ROT2.GT.0) MT=2
      IF(ROT1.LE.0.AND.ROT2.GT.0) MCT=1
      IF(ROT1.LE.0.AND.ROT2.LE.0) MCT=3
      IF(ROT1.LE.0.AND.ROT2.LE.0) MT=0
      IS(6,IRF)=MCT
      IF(MTRNS.GT.0) WRITE(LOUT,1061)MT,MCT
 1061 FORMAT(' MT,MCT',/,2I5)
      DO 255 K=1,MT
      IF(K+MCT.EQ.1) RO=ROT1
      IF(K+MCT.NE.1) RO=ROT2
      SQU=SQRT(RO)
      IF(MTRNS.GT.0) WRITE(LOUT,1062)SQU
      J=MT+1-K
      DO 256 L=1,3
      PS(L,J+3,IRF)=P(L)-(PUN+SQU)*UN(L)
  256 CONTINUE
  255 CONTINUE
      DO 257 K=1,MCT
      M=MT+K
      IF(M.EQ.3)CSQ=CSQRT(CMPLX(ROT1,0.0))
      IF(M.NE.3)CSQ=CSQRT(CMPLX(ROT2,0.0))
      IF(MTRNS.GT.0) WRITE(LOUT,1064)CSQ
      DO 258 L=1,3
      PS(L,M+3,IRF)=P(L)-(PUN+CSQ)*UN(L)
  258 CONTINUE
  257 CONTINUE
      IF(IASW.EQ.0) RETURN
      END IF
      IF(ROO.LE.0.0) GOTO 300
      SQU=SQRT(ROO)
      IF(ITRANS.EQ.1) GOTO 280
C
C  REFLECTION
C
      PU1=PUN-SQU
      DO 260 K=1,3
      PNEW(K)=P(K)-PU1*UN(K)
  260 CONTINUE
      GOTO 275
C
C  TRANSMISSION
C
  280 PU1=PUN+SQU
      DO 270 K=1,3
      PNEW(K)=P(K)-PU1*UN(K)
  270 CONTINUE
  275 IF(MTRNS.GT.0) WRITE(LOUT,1063)PNEW
 1063 FORMAT(' PNEW',/,3F12.6)
      PP=0.0
      DO 290 K=1,3
      PP=PP+Pnew(K)*Pnew(K)
  290 CONTINUE
      RHO2=0.2*SQRT(A11)+1.7
      IF(IRHO.NE.0) RHO2=RHO(LAY)
      PUNG=PNEW(1)*UN(1)+PNEW(2)*UN(2)+PNEW(3)*UN(3)
      DS(10,IRF)=abs(pung)/pp
      DS(11,IRF)=RHO2
      RETURN
  300 IND=9
      NTR=102
      RETURN
      END
C