C
C *** SYNTHETIC WAVES FROM THE HASKELL,S UNILATERAL DISLOCATION SOURCE
C     IN THE 3D LAYERED HALF-SPACE USING THE EXACT GREEN FUNCTION 
C     BASED ON THE RECIPROCITY RELATION ***
C
      SUBROUTINE HKUXYZ(DOM,MOM,NOM,ML,DNS,VP,VS,THK
     *                 ,LS,SDEP,LO,ODEP,MOB,NOB,ROB,RAD
     *                 ,FLL,FLW,STR,DIP,FLD,RST,VR,NVR
     *                 ,UX,UY,UZ)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VP(ML),VS(ML),CPX,SM,SM0,STL,STT,CCC,AAA
     *          ,UX(MOM,MOB),UY(MOM,MOB),UZ(MOM,MOB)
     *          ,QX1(0:3),QY1(0:3),QZ1(0:3),VX1(0:3),VY1(0:3)
     *          ,WX1(0:3),WY1(0:3),WZ1(0:3)
     *          ,QX2(0:3),QY2(0:3),QZ2(0:3),VX2(0:3),VY2(0:3)
     *          ,WX2(0:3),WY2(0:3),WZ2(0:3)
     *          ,XX(3),XY(3),XZ(3),YX(3),YY(3),YZ(3)
     *          ,ZX(3),ZY(3),ZZ(3)
      DIMENSION DNS(ML),THK(ML),ROB(MOB)
      COMMON / CN / C1,N1,C2,N2 /FDAT/ OM,NL,RR
C
      CPX=(0.D0,1.D0)
C
      CR=DCOS(RAD)
      SR=DSIN(RAD)
C
      CST= DCOS(STR)
      SST= DSIN(STR)
C
      EX= CST
      EY= SST
C
      CDP= DCOS(DIP)
      SDP= DSIN(DIP)
C
      XN=-SDP*SST
      YN= SDP*CST
      ZN=-CDP
C
      SM = DNS(LS)*VS(LS)*VS(LS)
      SM0= FLD*SM*FLL*FLW
C
      IF(NVR.EQ.1)THEN
        TL=0.D0
        WRITE(*,*)' TL=0 !!'
      ELSE
        TL=FLL/VR
      END IF
C
      DO 100 IOB=1,NOB
C
        RR=ROB(IOB)
        CRX=-SR*SR/RR
        CRY= SR*CR/RR
        SRX= CRY     
        SRY=-CR*CR/RR
C
        DO 110 IOM=1,NOM
C
          OM=DOM*DFLOAT(IOM)
          OM2=OM*OM
          WRITE(*,*)' OM =',OM
C
          TTK= OM*RST/2.D0
          IF(TTK.EQ.0.D0)THEN
            STT= (1.D0,0.D0)
          ELSE
            STT= DSIN(TTK)*CDEXP(-CPX*TTK)/TTK
          END IF
C
          TTL= OM*TL/2.D0
          IF(TTL.EQ.0.D0)THEN
            STL= (1.D0,0.D0)
          ELSE
            STL= DSIN(TTL)*CDEXP(-CPX*TTL)/TTL
          END IF
C
          CCC=-CPX*SM0*STL*STT/OM
C
          WNB=OM/C2
          WNA=WNB/(1000.D0*N1)
          WNC=OM/C1
C
          IF((LO.LT.LS).OR.((LO.EQ.LS).AND.(ODEP.LE.SDEP)))THEN
C
            CALL QVWXYZ(WNA,WNB,N1,DNS,VP,VS,THK,ML,IER
     *           ,LS,SDEP,LO,ODEP,1,CR,SR,CRX,SRX,CRY,SRY
     *           ,QX1,QY1,QZ1,VX1,VY1,WX1,WY1,WZ1)
            IF(IER.NE.0)WRITE(*,*)' 1:IOB,IOM,IER',IOB,IOM,IER
            CALL QVWXYZ(WNB,WNC,N2,DNS,VP,VS,THK,ML,IER
     *           ,LS,SDEP,LO,ODEP,1,CR,SR,CRX,SRX,CRY,SRY
     *           ,QX2,QY2,QZ2,VX2,VY2,WX2,WY2,WZ2)
            IF(IER.NE.0)WRITE(*,*)' 2:IOB,IOM,IER',IOB,IOM,IER
C
          ELSE
C
            CALL QVWXYZ(WNA,WNB,N1,DNS,VP,VS,THK,ML,IER
     *           ,LO,ODEP,LS,SDEP,2,CR,SR,CRX,SRX,CRY,SRY
     *           ,QX1,QY1,QZ1,VX1,VY1,WX1,WY1,WZ1)
            IF(IER.NE.0)WRITE(*,*)' 1:IOB,IOM,IER',IOB,IOM,IER
            CALL QVWXYZ(WNB,WNC,N2,DNS,VP,VS,THK,ML,IER
     *           ,LO,ODEP,LS,SDEP,2,CR,SR,CRX,SRX,CRY,SRY
     *           ,QX2,QY2,QZ2,VX2,VY2,WX2,WY2,WZ2)
            IF(IER.NE.0)WRITE(*,*)' 2:IOB,IOM,IER',IOB,IOM,IER
C
          END IF
C
          XX(1)=(QX1(1)+QX2(1))*CR -(VX1(1)+VX2(1))*SR
     *          +(QX1(0)+QX2(0))*CRX-(VX1(0)+VX2(0))*SRX
          YX(1)=(QX1(1)+QX2(1))*SR +(VX1(1)+VX2(1))*CR
     *          +(QX1(0)+QX2(0))*SRX+(VX1(0)+VX2(0))*CRX
          ZX(1)= WX1(1)+WX2(1)
          XY(1)=(QY1(1)+QY2(1))*CR -(VY1(1)+VY2(1))*SR
     *          +(QY1(0)+QY2(0))*CRX-(VY1(0)+VY2(0))*SRX
          YY(1)=(QY1(1)+QY2(1))*SR +(VY1(1)+VY2(1))*CR
     *          +(QY1(0)+QY2(0))*SRX+(VY1(0)+VY2(0))*CRX
          ZY(1)= WY1(1)+WY2(1)
          XZ(1)=(QZ1(1)+QZ2(1))*CR +(QZ1(0)+QZ2(0))*CRX
          YZ(1)=(QZ1(1)+QZ2(1))*SR +(QZ1(0)+QZ2(0))*SRX
          ZZ(1)= WZ1(1)+WZ2(1)
C
          XX(2)=(QX1(2)+QX2(2))*CR -(VX1(2)+VX2(2))*SR
     *          +(QX1(0)+QX2(0))*CRY-(VX1(0)+VX2(0))*SRY
          YX(2)=(QX1(2)+QX2(2))*SR +(VX1(2)+VX2(2))*CR
     *          +(QX1(0)+QX2(0))*SRY+(VX1(0)+VX2(0))*CRY
          ZX(2)= WX1(2)+WX2(2)
          XY(2)=(QY1(2)+QY2(2))*CR -(VY1(2)+VY2(2))*SR
     *          +(QY1(0)+QY2(0))*CRY-(VY1(0)+VY2(0))*SRY
          YY(2)=(QY1(2)+QY2(2))*SR +(VY1(2)+VY2(2))*CR
     *          +(QY1(0)+QY2(0))*SRY+(VY1(0)+VY2(0))*CRY
          ZY(2)= WY1(2)+WY2(2)
          XZ(2)=(QZ1(2)+QZ2(2))*CR +(QZ1(0)+QZ2(0))*CRY
          YZ(2)=(QZ1(2)+QZ2(2))*SR +(QZ1(0)+QZ2(0))*SRY
          ZZ(2)= WZ1(2)+WZ2(2)
C
          XX(3)=(QX1(3)+QX2(3))*CR-(VX1(3)+VX2(3))*SR
          YX(3)=(QX1(3)+QX2(3))*SR+(VX1(3)+VX2(3))*CR
          ZX(3)= WX1(3)+WX2(3)
          XY(3)=(QY1(3)+QY2(3))*CR-(VY1(3)+VY2(3))*SR
          YY(3)=(QY1(3)+QY2(3))*SR+(VY1(3)+VY2(3))*CR
          ZY(3)= WY1(3)+WY2(3)
C          XZ(3)=(QZ1(3)+QZ2(3))*CR
C          YZ(3)=(QZ1(3)+QZ2(3))*SR
C          ZZ(3)= WZ1(3)+WZ2(3)
C
          AAA= 2.D0*EX*XN*XX(1)+(EX*YN+EY*XN)*XX(2)+EX*ZN*XX(3)
     *        +(EY*XN+EX*YN)*XY(1)+2.D0*EY*YN*XY(2)+EY*ZN*XY(3)
     *        +(EZ*XN+EX*ZN)*XZ(1)+(EZ*YN+EY*ZN)*XZ(2)
          UX(IOM,IOB)= CCC*AAA
C
          AAA= 2.D0*EX*XN*YX(1)+(EX*YN+EY*XN)*YX(2)+EX*ZN*YX(3)
     *        +(EY*XN+EX*YN)*YY(1)+2.D0*EY*YN*YY(2)+EY*ZN*YY(3)
     *        +(EZ*XN+EX*ZN)*YZ(1)+(EZ*YN+EY*ZN)*YZ(2)
          UY(IOM,IOB)= CCC*AAA
C
          AAA= 2.D0*EX*XN*ZX(1)+(EX*YN+EY*XN)*ZX(2)+EX*ZN*ZX(3)
     *        +(EY*XN+EX*YN)*ZY(1)+2.D0*EY*YN*ZY(2)+EY*ZN*ZY(3)
     *        +(EZ*XN+EX*ZN)*ZZ(1)+(EZ*YN+EY*ZN)*ZZ(2)
          UZ(IOM,IOB)= CCC*AAA
C
  110   CONTINUE
C
  100 CONTINUE
C
      RETURN
      END
