      SUBROUTINE DISPL(P,U,I1)
C
C  ROUTINE FOR THE COMPUTATION OF DISPLACEMENT VECTOR OF GENERATED
C  WAVES AT AN INTERFACE. IT WORKS EVEN FOR A COMPLEX REFRACTION VECTOR.
C  ROUTINES CHRM AND POLAR WORK ONLY FOR REAL REFRACTION VECTORS
C
C  INPUT PARAMETERS:
C  P(1-3)... SLOWNESS VECTOR
C  I1... TYPE OF THE WAVE WITH THE SLOWNESS VECTOR P
C        I1=1... REFLECTED S WAVE (SV COMPONENT IN ISOTROPIC CASE)
C        I1=2... REFLECTED S WAVE (SH COMPONENT IN ISOTROPIC CASE)
C        I1=3... REFLECTED P WAVE
C        I1=4... TRANSMITTED S WAVE (SV COMPONENT IN ISOTROPIC CASE)
C        I1=5... TRANSMITTED S WAVE (SH COMPONENT IN ISOTROPIC CASE)
C        I1=6... TRANSMITTED P WAVE
C
C  OUTPUT PARAMETERS:
C  U(1-3)... CORRESPONDING DISPLACEMENT VECTOR
C
      complex uxn(3,3)
      COMPLEX U,U1,U2,U3,P,P1,P2,P3,P2P3,P1P2,P1P3,P1P1,P2P2,P3P3,C11,
     1  C12,C13,C22,C23,C33,C11N,C22N,C33N,C23SQ,C13SQ,C12SQ,CD11,CD12,
     2  CD13,CD22,CD23,CD33,CZ01,CZ02,CZ03,CDTR,cp1,cp2,cp3,cab
      DIMENSION U(3),P(3),VSQ(3),Y(6),UN(3),xx(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,IOUT,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
      COMPLEX PS
      COMMON /RAY/   AY(28,400),DS(20,20),KINT(20),HHH(3,3),tmax,
     1               PS(3,7,20),IS(8,20),DINC,DTOLD,N,IREF,IND,IND1
C
      IRF=IREF
      IF(I1.LE.3)THEN
        MC=IS(5,IRF)
      ELSE
        MC=IS(6,IRF)
      END IF
      ISW=0
C
C  THIS PART CORRESPONDS TO CHRM ,BUT FOR COMPLEX REFRACTION VECTORS
C
 900  P1=P(1)
      P2=P(2)
      P3=P(3)
 910  continue     
      IF(ICOEF.GT.0)WRITE(LOUT,'(A,3(/2F12.6))')'DISPL: START P',
     1p1,p2,p3 
      P2P3=P2*P3
      P1P2=P1*P2
      P1P3=P1*P3
      P1P1=P1*P1
      P2P2=P2*P2
      P3P3=P3*P3
      C11=P1P1*A11+P2P2*A66+P3P3*A55
     1+2.*P2P3*A56+2.*P1P3*A15+2.*P1P2*A16
      C22=P1P1*A66+P2P2*A22+P3P3*A44
     1+2.*P2P3*A24+2.*P1P3*A46+2.*P1P2*A26
      C33=P1P1*A55+P2P2*A44+P3P3*A33
     1+2.*P2P3*A34+2.*P1P3*A35+2.*P1P2*A45
      C23=P1P1*A56+P2P2*A24+P3P3*A34
     1   +P2P3*A2344+P1P3*A3645+P1P2*A2546
      C13=P1P1*A15+P2P2*A46+P3P3*A35
     1   +P2P3*A3645+P1P3*A1355+P1P2*A1456
      C12=P1P1*A16+P2P2*A26+P3P3*A45
     1   +P2P3*A2546+P1P3*A1456+P1P2*A1266
      C11N=C11-1.
      C22N=C22-1.
      C33N=C33-1.
      C23SQ=C23*C23
      C13SQ=C13*C13
      C12SQ=C12*C12
      CD11=C22N*C33N-C23SQ
      CD22=C11N*C33N-C13SQ
      CD33=C11N*C22N-C12SQ
      CD12=C13*C23-C12*C33N
      CD13=C12*C23-C13*C22N
      CD23=C12*C13-C23*C11N
      CDTR=CD11+CD22+CD33
      IF(ICOEF.GT.0)
     1WRITE(LOUT,'(A,4(/4F12.6))')'DISPL: CDIJ'
     1,CD11,CD12,CD13,CD22,CD23,CD33,CDTR
      if(abs(cdtr).lt..0000001)write(lout,'(A)')'DISPL: SHEAR WAVE
     1 SINGULARITY IN CALCULATION OF R/T COEFFICIENT' 
C
C  THIS PART CORRESPONDS TO ROUTINE POLAR BUT FOR COMPLEX POLARIZATION
C  VECTORS
C
      U1=CD11
      U2=CD12
      U3=CD13
      uxn(1,1)=u1
      uxn(1,2)=u2
      uxn(1,3)=u3
      XN=real(U1*CONJG(U1)+U2*CONJG(U2)+U3*CONJG(U3))
      xx(1)=xn
      IF(ICOEF.GT.0)WRITE(LOUT,'(A,F14.7)')'DISPL: XN',XN
      U1=CD12
      U2=CD22
      U3=CD23
      uxn(2,1)=u1
      uxn(2,2)=u2
      uxn(2,3)=u3
      XN=real(U1*CONJG(U1)+U2*CONJG(U2)+U3*CONJG(U3))
      xx(2)=xn
      IF(ICOEF.GT.0)WRITE(LOUT,'(A,F14.7)')'DISPL: XN',XN
      U1=CD13
      U2=CD23
      U3=CD33
      uxn(3,1)=u1
      uxn(3,2)=u2
      uxn(3,3)=u3
      XN=real(U1*CONJG(U1)+U2*CONJG(U2)+U3*CONJG(U3))
      xx(3)=xn
      IF(ICOEF.GT.0)WRITE(LOUT,'(A,F14.7)')'DISPL: XN',XN
      xn=0.
      do 920 l=1,3
        if(xn.ge.xx(l))go to 920
        xn=xx(l)
        nx=l
  920 continue
      xn=1./sqrt(xn)
      u1=uxn(nx,1)*xn      
      u2=uxn(nx,2)*xn      
      u3=uxn(nx,3)*xn      
      U(1)=U1
      U(2)=U2
      U(3)=U3
      IF(ICOEF.GT.0)WRITE(LOUT,'(A,3(/2F12.6))')'DISPL: U',U
C
C  CHECK OF PRECISION
C
      CZ01=C11N*u1+C12*u2+C13*u3
      CZ02=C12*u1+C22N*u2+C23*u3
      CZ03=C13*u1+C23*u2+C33N*u3
      IF(ICOEF.GT.0)
     1WRITE(LOUT,'(A,3(/2F14.7))')'DISPL: PRECISSION OF DISPL'
     1,CZ01,CZ02,CZ03
      IF(ISW.GT.0)THEN
        GOTO 1100
      END IF
      RETURN
C
C  COMPUTATION OF DISPLACEMENT VECTORS AT AN INTERFACE IF A SHEAR WAVE
C  SINGULARITY OCCURS.
C  FIRST: THE DISPLACEMENT VECTOR OF THE QP-WAVE IS COMPUTED FOR THE
C         NORMAL DIRECTION OF THE S-WAVE
C  SECOND: TWO DISPLACEMENT VECTORS ARE CONSTRUCTED SO THAT THEY ARE
C          PERPENDICULAR TO THE QP-WAVE DISPLACEMENT VECTOR
C          (THEY CORRESPOND TO SV AND SH COMPONENT)
C
 1000 ISW=1
      do 950 k=1,3
      Y(k)=AY(k+1,N)
      UN(k)=DS(k,IRF)
  950 continue
      Y(4)=REAL(P1)
      Y(5)=REAL(P2)
      Y(6)=REAL(P3)
      IF(ICOEF.GT.0)
     1WRITE(LOUT,'(A,3I5,6F12.6)')'DISPL: I1,IRF,MC,Y',I1,IRF,MC,Y
      UUN=Y(4)*Y(4)+Y(5)*Y(5)+Y(6)*Y(6)
      UUN=1./SQRT(UUN)
      Y(4)=Y(4)*UUN
      Y(5)=Y(5)*UUN
      Y(6)=Y(6)*UUN
      CALL INIT(Y,VSQ)
      V=AMIN1(VSQ(1),VSQ(2),VSQ(3))
      V=1./SQRT(V)
      IF(ICOEF.GT.0)WRITE(LOUT,'(A,3F12.4)')'DISPL: VSQ(1-3)=',VSQ
      IF(ICOEF.GT.0)WRITE(LOUT,'(A,F12.4)')'DISPL: V=',V
      p1=y(4)/v
      p2=y(5)/v
      p3=y(6)/v
      go to 910
 1100 if(i1.eq.1.or.i1.eq.4)then
        cp1=U1*U3
        cp2=U2*U3
        cp3=-(U1*U1+U2*U2)
        cab=sqrt(real(cp1*conjg(cp1)+cp2*conjg(cp2)+cp3*conjg(cp3)))
        U(1)=CP1/cab
        U(2)=CP2/cab
        U(3)=cp3/cab
        IF(ICOEF.GT.0)WRITE(LOUT,'(A,7(/2F12.6))')'DISPL: U,P',U,P
        RETURN
      end if
      if(i1.eq.2.or.i1.eq.5)then
        cp1=-U2
        cp2=U1
        cp3=cmplx(0.,0.)
        cab=sqrt(real(CP1*conjg(CP1)+CP2*conjg(CP2)+CP3*conjg(CP3)))
        U(1)=CP1/cab
        U(2)=CP2/cab
        U(3)=cmplx(0.,0.)
        IF(ICOEF.GT.0)WRITE(LOUT,'(A,7(/2F12.6))')' DISPL U,P',U,P
        RETURN
      end if
      END