C     
C *** SYNTHESIZED WAVES FROM MULTI-DIPOLE SOURCES IN A LAYERED HALF-SPACE 
C     USING THE R/T COEFFICIENT METHOD AND THE ASYMPTOTIC METHOD 
C     (HISADA, Bull. Seismo. Soc. Am., V.84, Oct. 1994, for Part 1)
C     (HISADA, Bull. Seismo. Soc. Am., V.85, Aug. 1995, for Part 2)
C     (HISADA, Bull. Seismo. Soc. Am., V.85, Oct. 1995, for comments)
C
C *** FOR BUGS OR COMMENTS, PLEASE E-MAIL IN "hisada@cc.kogakuin.ac.jp",
C     OR MAIL TO "Yoshiaki Hisada, Dept. of Architecture,
C     Kogakuin Univeristy, Nishi-Shinjuku 1-24-2, Tokyo 163-91, Japan"
C
C                      (made in Jan, 1995; Modififed in Aug, 1995)
C
      PROGRAM grfault
C
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(MOM=300,MOB=20)
      COMPLEX*16 GX(MOM,MOB),GY(MOM,MOB),GZ(MOM,MOB)
C
      WRITE(*,1000)
 1000 FORMAT(//' GRFAULT: TO COMPUTE WAVES FROM MULTI-HASKELL-SOURCES'/
     *'   IN A LAYERED HALF-SPACE USING THE R/T MATRIX METHOD'/
     *'   AND THE ASYMPTOTIC TECHNIQUE BY HISADA (1995).'/
     *'   FOR QUESTIONS OR BUGS, E-MAIL TO "hisada@cc.kogakuin.ac.jp".')
C
      CALL INDAT(MOM,MOB)
      CALL HKQVW(MOM,MOB,GX,GY,GZ)
      CALL OTDAT(MOM,MOB,GX,GY,GZ)
C
      WRITE(*,*)
      WRITE(*,*)
      WRITE(*,*)
     *'Output data are stored in "grfault.dxyz" and "grfault.dqvw".'
      WRITE(*,*)
      WRITE(*,*)
     *'Please use "grfft.f" to obtain seismograms in the time domain.'
      WRITE(*,*)
C
      STOP
      END
C
C
C *** DATA INPUT ***
C
      SUBROUTINE INDAT(MOM,MOB)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VP,VS
      COMMON /FREQ/DOM,IOM,NOM /OBSP/ODEP,XO(20),YO(20),NOB,LO
     * /STRC/DNS(20),VP(20),VS(20),THK(20),NL
     *    /INTG/C11,C12,N11,N12 /INTG2/WNB2,WNC2,I12,N21,N22
     *    /FLT/XF(30),YF(30),LF(30),DP(30),STR(30),DIP(30),RAK(30)
     *      ,FL(30),FW(30),FD(30),RT(30),VR(30),NVR(30),TD(30),NF,IS 
     *    /SIGN/ISIGN /OTST/IOT,IST
C
C *** PARAMETERS IN COMMOM STATEMENT ***
C
      NML =20
      NMOB=20
      NFLT=30
      IF(MOB.LT.NMOB)NMOB=MOB
C
C *** OMEGA DATA ***
C
      OPEN (UNIT=7,FILE='grfault.in',STATUS='OLD')
C
      READ(7,1100)IOM,NOM,DOM
      WRITE(10,*)'FROM grfault.in:'
      WRITE(10,*)
      WRITE(10,*)
      WRITE(10,*)'IOM,NOM,DOM=',IOM,NOM,DOM
 1100 FORMAT(/2I10,F10.0)
C
      IF(NOM.GT.MOM)THEN
        WRITE(10,*)' NOM>MOM !!!!!!!!!!!!!'
        STOP
      END IF
C
C *** MEDIUM DATA ***
C
      READ(7,1200)NL
      DO 100 IL=1,NL
        READ(7,1210)DNS(IL),VPP,QP,VSS,QS,THK(IL)
        VP(IL)=VPP*DCMPLX(1.D0,-0.5D0/QP)
        VS(IL)=VSS*DCMPLX(1.D0,-0.5D0/QS)
  100 CONTINUE
 1200 FORMAT(/I10)
 1210 FORMAT(6E8.0)
C
      WRITE(10,2200)NL
      DO 200 IL=1,NL
        WRITE(10,2210)IL,DNS(IL),VP(IL),VS(IL),THK(IL)
  200 CONTINUE
 2200 FORMAT(//' Total Layers Number            =',I5/
     *       /' Layer Num., Density , Vp(real), Vp(imag)'
     *        ', Vs(real), Vs(imag),Thickness')
 2210 FORMAT(I5,5X,6F10.3)
C
      IF(NL.GT.NML)THEN
        WRITE(10,*)' NL>NML !!!!!!!!!!!!!'
        STOP
      END IF
C
C *** FAULT DATA ***
C
      READ(7,1300)NF
C
      DO 250 I=1,NF
C
        READ(7,1305)TD(I)
        READ(7,1310)XF(I),YF(I)
        READ(7,1320)DP(I),STR(I)
        READ(7,1320)DIP(I),RAK(I)
        STR(I)=STR(I)*0.01745329252d0
        DIP(I)=DIP(I)*0.01745329252d0
        RAK(I)=RAK(I)*0.01745329252d0
        READ(7,1320)FL(I),FW(I),FD(I)
        READ(7,1330)RT(I),VR(I)
        READ(7,1340)NVR(I)
C
        DRU=0.0
        DRL=THK(1)
        DO 260 L=1,NL
          IF((DP(I).GE.DRU).AND.(DP(I).LT.DRL))LF(I)=L
          IF((DP(I).GE.DRU).AND.(L.EQ.NL))     LF(I)=L
          DRU=DRL
          DRL=DRL+THK(L+1)
  260   CONTINUE
C
  250 CONTINUE
C
 1300 FORMAT(/I10)
 1305 FORMAT(F10.0)
 1310 FORMAT(2F10.0)
 1320 FORMAT(3F10.0)
 1330 FORMAT(2F10.0)
 1340 FORMAT(I10)
C
      WRITE(10,2300)NF
      DO 270 I=1,NF
      WRITE(10,2305)I
      WRITE(10,2310)XF(I),YF(I)
      WRITE(10,2320)DP(I),LF(I)
      WRITE(10,2330)STR(I),DIP(I)
      WRITE(10,2335)RAK(I)
      WRITE(10,2340)FL(I),FW(I)
      WRITE(10,2350)FD(I),RT(I)
      IF(NVR(I).NE.1)THEN
        WRITE(10,2360)VR(I)
      ELSE
        WRITE(10,2370)
      END IF
  270 CONTINUE
C
 2300 FORMAT(//' Source Data (MULTI-HASKELL''S SOURCES):'
     *        '  : Number of Source =',I3)
 2305 FORMAT(/'   Source Number                            =',I5)
 2310 FORMAT(/'    Source X-Axis(m)=',F12.3,', Y-Axis (m)   =',F12.3)
 2320 FORMAT( '    Source Depth (m)=',F12.3,', Source Layer =',I5)
 2330 FORMAT( '    Strike Angle(rd)=',F12.8,', Dip Angle(rd)=',F12.8)
 2335 FORMAT( '    Rake Angle  (rd)=',F12.8)
 2340 FORMAT( '    Fault Length (m)=',F12.3,', Width(m)     =',F12.3)
 2350 FORMAT( '    Dislocation  (m)=',F12.5,', Rise Time(s) =',F12.8)
 2360 FORMAT( '    Rupture Velocity (m/s) =',F12.3)
 2370 FORMAT( '    Rupture Velocity (m/s) = Infinite ')
C
      IF(NF.GT.NFLT)THEN
        WRITE(10,*)' NF.GT.NFLT !!!!!!'
        STOP
      END IF
C
C *** NUMERICAL INTEGRAL DATA ***
C
      READ(7,1400)I12
      WRITE(10,2400)I12
 1400 FORMAT(/I10)
 2400 FORMAT(//' INTEGRAL DATA: PATTERN=',I10)
C
      IF(I12.EQ.1)THEN
C
        READ(7,1410)C11,N11,C12,N12
        READ(7,*)
        WRITE(10,2410)C11,N11,C12,N12
        WRITE(10,2420)N11*2+1,N12*2+1
C
      ELSE
C
        READ(7,*)
        READ(7,1410)WNB2,N21,WNC2,N22
        WRITE(10,2430)WNB2,N21,WNC2,N22
        WRITE(10,2420)N11*2+1,N12*2+1
C
      END IF
C
 1410 FORMAT(F10.0,I10,F10.0,I10)
 2410 FORMAT(/'  C1,N1,C2,N2=',F10.2,I5,F10.2,I5)
 2420 FORMAT( '  ACTUAL INTEGRAL POINTS FOR N1 AND N2=',2I5)
 2430 FORMAT(/'  WNB,N1,WNC,N2=',F10.7,I5,F10.7,I5)
C
C *** OP DATA ***
C
      READ(7,1500)LO,ODEP
      WRITE(10,2500)LO,ODEP
 1500 FORMAT(/I10,F10.0)
 2500 FORMAT(//' OP LAYER, OP DEPTH FROM UPPER BOUNDARY=',I5,F10.3)
C
      IF((ODEP.EQ.THK(LO)).AND.(LO.NE.NL))THEN
        ODEP=0.D0
        LO=LO+1
        WRITE(10,2510)LO,ODEP
      END IF
 2510 FORMAT(/' NEW OP LAYER, OP DEPTH FROM UPPER BOUNDARY=',I5,F10.3)
C
      READ(7,1510)NOB
C
      IF(NOB.GT.NMOB)THEN
        WRITE(10,*)' NOB>NMOB !!!!!!!!!!!!!'
        STOP
      END IF
C
      DO 300 I=1,NOB,2
        READ(7,1520)XO(I),YO(I),XO(I+1),YO(I+1)
  300 CONTINUE
C
 1510 FORMAT(I10)
 1520 FORMAT(4F10.0)
C
C *** DATA FOR THE ASYMPTOTIC ***
C
      READ(7,1600)IST
      IF(IST.EQ.1)THEN
        WRITE(10,2600)
      ELSE
        WRITE(10,2610)
      END IF
 1600 FORMAT(/I10)
 2600 FORMAT(//' THE ASYMPTOTIC METHOD IS USED')
 2610 FORMAT(//' THE ASYMPTOTIC METHOD IS NOT USED')
C
C *** CHANGE OF SIGNS FOR IMAGINARY PARTS ***
C
      READ(7,1700)ISIGN
      IF(ISIGN.EQ.1)WRITE(10,2700)
      IF(ISIGN.EQ.0)WRITE(10,2710)
 1700 FORMAT(/I10)
 2700 FORMAT(//' SIGNS OF IMAGINARY PARTS FOR RESULTS ARE CHANGED')
 2710 FORMAT(//' SIGNS OF IMAGINARY PARTS FOR RESULTS ARE NOT CHANGED')
C
      CLOSE(UNIT=7)
C
C *** OUTPUT OF INTEGRAND OF U,V,AND W ***
C
   50 WRITE(*,7000)
      READ(*,*)IOT
 7000 FORMAT(/' Output Integrands?: Enter No=0 or Yes=1: ')
C
      IF((IOT.EQ.1).AND.(NOB.NE.1))THEN
        WRITE(*,*)' NOB.NE.1 !!!!'
        GOTO 50
      END IF
      IF((IOT.EQ.1).AND.(NOM.NE.1))THEN
        WRITE(*,*)' NOM.NE.1 !!!!'
        GOTO 50
      END IF
C
      RETURN
      END
C
C
C
C *** DATA OUTPUT ***
C
      SUBROUTINE OTDAT(MOM,MOB,GX,GY,GZ)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VP,VS,GX(MOM,MOB),GY(MOM,MOB),GZ(MOM,MOB),AA,CC
      COMMON /FREQ/DOM,IOM,NOM /OBSP/ODEP,XO(20),YO(20),NOB,LO
     *    /STRC/DNS(20),VP(20),VS(20),THK(20),NL
     *    /FLT/XF(30),YF(30),LF(30),DP(30),STR(30),DIP(30),RAK(30)
     *      ,FL(30),FW(30),FD(30),RT(30),VR(30),NVR(30),TD(30),NF,IS 
     *    /SIGN/ISIGN
C
C *** CHECK ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(10,*)'FROM OTDAT:NL>NC',NL,NC
         STOP
      END IF
C
      NC=20
      IF(NOB.GT.NC)THEN
         WRITE(10,*)'FROM OTDAT:NOB>NC',NOB,NC
         STOP
      END IF
C
      NC=30
      IF(NF.GT.NC)THEN
         WRITE(10,*)'FROM OTDAT:NF>NC',NF,NC
         STOP
      END IF
C
C *** OUTPUT OF DATA (GX, GY AND GZ)***
C
      IF(ISIGN.EQ.1)THEN
        DO 100 IOB=1,NOB
        DO 110 K=IOM,NOM
          GX(K,IOB)=DCMPLX(DREAL(GX(K,IOB)),-DIMAG(GX(K,IOB)))
          GY(K,IOB)=DCMPLX(DREAL(GY(K,IOB)),-DIMAG(GY(K,IOB)))
          GZ(K,IOB)=DCMPLX(DREAL(GZ(K,IOB)),-DIMAG(GZ(K,IOB)))
  110   CONTINUE
  100   CONTINUE
      END IF
C      
      OPEN(UNIT=8,FILE='grfault.dxyz',STATUS='UNKNOWN')
        WRITE(8,2500)DOM,IOM,NOM
        WRITE(8,2510)NOB
        DO 500 IOB=1,NOB
          DO 510 K=IOM,NOM
            WRITE(8,2520)GX(K,IOB),GY(K,IOB),GZ(K,IOB)
  510     CONTINUE
  500   CONTINUE
      CLOSE(8)
C
 2500 FORMAT(1X,F10.6,2I10)
 2510 FORMAT(1X,I10,///)
 2520 FORMAT(1X,6E12.5)
C
C *** OUTPUT OF DATA (GQ, GV AND GW)***
C
        DO 600 IOB=1,NOB
          RX=XO(IOB)-XF(1)
          RY=YO(IOB)-YF(1)
          RR=SQRT(RX*RX+RY*RY)
          CR=RX/RR
          SR=RY/RR
          DO 610 K=IOM,NOM
            AA=GX(K,IOB)
            CC=GY(K,IOB)
            GX(K,IOB)= AA*CR+CC*SR
            GY(K,IOB)=-AA*SR+CC*CR
  610     CONTINUE
  600   CONTINUE
C      
      OPEN(8,FILE='grfault.dqvw',STATUS='UNKNOWN')
        WRITE(8,3500)DOM,IOM,NOM
        WRITE(8,3510)NOB
        DO 800 IOB=1,NOB
          DO 810 K=IOM,NOM
            WRITE(8,3520)GX(K,IOB),GY(K,IOB),GZ(K,IOB)
  810     CONTINUE
  800   CONTINUE
      CLOSE(8)
C
 3500 FORMAT(1X,F10.6,2I10)
 3510 FORMAT(1X,I10,///)
 3520 FORMAT(1X,6E12.5)
C
      RETURN
      END
C
C
C
C *** COMPUTATION OF GREEN FUNCTIONS USING THE R/T METHOD ***
C
      SUBROUTINE HKQVW(MOM,MOB,GGX,GGY,GGZ)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 GGX(MOM,MOB),GGY(MOM,MOB),GGZ(MOM,MOB),VP,VS
     *      ,QX1(0:2),QZ1(0:2),VX1(0:2),WX1(0:2),WZ1(0:2)
     *      ,QX2(0:2),QZ2(0:2),VX2(0:2),WX2(0:2),WZ2(0:2)
     *      ,QXS(0:2),QZS(0:2),VXS(0:2),WXS(0:2),WZS(0:2)
     *      ,SM,RAM,SMO,EDT,GX,GY,GZ,QA,QB,VA,VB,WA,WB
     *      ,QX(0:3),QY(0:3),QZ(0:3),VX(0:3),VY(0:3)
     *      ,WX(0:3),WY(0:3),WZ(0:3),UX(3,3),UY(3,3),UZ(3,3)
      DIMENSION EN(3,3)
      COMMON /FREQ/DOM,IOM,NOM /OBSP/ODEP,XO(20),YO(20),NOB,LO
     *    /STRC/DNS(20),VP(20),VS(20),THK(20),NL
     *    /INTG/C11,C12,N11,N12 /INTG2/WNB2,WNC2,I12,N21,N22
     *    /FLT/XF(30),YF(30),LF(30),DP(30),STR(30),DIP(30),RAK(30)
     *      ,FL(30),FW(30),FD(30),RT(30),VR(30),NVR(30),TD(30),NF,IS 
     *    /OTST/IOT,IST  /FCQVW/SM,RAM,RR  /SOLY/SDEP,FDEP,LS,LFL
C
C *** CHECK OF DIMENSION ***
C
      NC=20
      IF(NL.GE.NC)THEN
        WRITE(*,*)' FROM HKQVW:NL>NC !!!!!!!!!!'
        WRITE(*,*)' NL,NC=',NL,NC
        STOP
      END IF
C
      NC=20
      IF(NOB.GE.NC)THEN
        WRITE(*,*)' FROM HKQVW:NOB>NC !!!!!!!!!!'
        WRITE(*,*)' NOB,NC=',NOB,NC
        STOP
      END IF
C
      NC=30
      IF(NF.GE.NC)THEN
        WRITE(*,*)' FROM HKQVW:NF>NC !!!!!!!!!!'
        WRITE(*,*)' NF,NC=',NF,NC
        STOP
      END IF
C
C *** GREEN'S FUNCTION ***
C
C   * output of integrands for check *
C
      IF(IOT.eq.1)THEN
      OPEN (10,FILE='QX.OUT',STATUS='NEW')
      OPEN (11,FILE='VX.OUT',STATUS='NEW')
      OPEN (12,FILE='WX.OUT',STATUS='NEW')
      OPEN (13,FILE='QZ.OUT',STATUS='NEW')
      OPEN (14,FILE='WZ.OUT',STATUS='NEW')
      OPEN (15,FILE='QXR.OUT',STATUS='NEW')
      OPEN (16,FILE='VXR.OUT',STATUS='NEW')
      OPEN (17,FILE='WXR.OUT',STATUS='NEW')
      OPEN (18,FILE='QZR.OUT',STATUS='NEW')
      OPEN (19,FILE='WZR.OUT',STATUS='NEW')
      OPEN (20,FILE='QXZ.OUT',STATUS='NEW')
      OPEN (21,FILE='VXZ.OUT',STATUS='NEW')
      OPEN (22,FILE='WXZ.OUT',STATUS='NEW')
      OPEN (23,FILE='QZZ.OUT',STATUS='NEW')
      OPEN (24,FILE='WZZ.OUT',STATUS='NEW')
      END IF
C
c   * the layer and depth of the observation points *
c
      LS=LO
      SDEP=ODEP
C
      WRITE(*,*)
      DO 100 IOB=1,NOB
c
        WRITE(*,*)' Now, Receiver Number =',IOB
C
        DO 110 IO=IOM,NOM
C
          OM=DOM*DFLOAT(IO)
          WRITE(*,*)'   OM=',OM
C
          IF(I12.EQ.1)THEN
            N1=N11
            N2=N12
            WNB=OM/C12
            WNA=WNB/(1000.D0*N1)
c            WRITE(*,*)'HKQVW:WNA=WNB/(1000*N1)'
            WNC=OM/C11
          ELSE
            N1=N21
            N2=N22
            WNB=WNB2
            WNA=WNB/(1000.D0*N1)
c            WRITE(*,*)'HKQVW:WNA=WNB/(1000*N1)'
            WNC=WNC2
          END IF
C
          GGX(IO,IOB)=(0.D0,0.D0)
          GGY(IO,IOB)=(0.D0,0.D0)
          GGZ(IO,IOB)=(0.D0,0.D0)
C
          DO 120 IS=1,NF
C
            RX=XF(IS)-XO(IOB)
            RY=YF(IS)-YO(IOB)
            RR=SQRT(RX*RX+RY*RY)
C
            LFL=LF(IS)
            FDEP=DP(IS)
            DO 122 L=2,LF(IS)
              FDEP=FDEP-THK(L-1)
  122       CONTINUE
C
            CALL ENMO(OM,EN,SMO)
            EDT=CDEXP(DCMPLX(0.D0,OM*TD(IS)))
            SMO=SMO*EDT
C
c         * the wavenumber integration for range 1 where the poles are included *
            CALL QVWRZ(OM,WNA,WNB,N1,QX1,QZ1,VX1,WX1,WZ1)
C
c         * the wavenumber integration for range 2 where the poles are not included*
            CALL QVWRZ(OM,WNB,WNC,N2,QX2,QZ2,VX2,WX2,WZ2)
C
            DO 126 J=0,2
              QXS(J)=(0.D0,0.D0)
              VXS(J)=(0.D0,0.D0)
              WXS(J)=(0.D0,0.D0)
              QZS(J)=(0.D0,0.D0)
              WZS(J)=(0.D0,0.D0)
  126       CONTINUE   
C
c         * the wavenumber integration of the asymptotic solutions *
            IF((IST.EQ.1).AND.(IABS(LS-LFL).LE.1))THEN
              CALL AINTRZ(LS,NL,QXS,QZS,VXS,WXS,WZS)
            END IF
C
c         * see Hisada(1994) eq.(26) *
            CR=RX/RR
            SR=RY/RR
            C2=CR*CR
            S2=SR*SR
            CS=CR*SR
            CX=S2/RR
            SY=C2/RR
            CY=-CS/RR
            SX=CY
C
c         * Urx, Ury, Urz: see Hisada(1994) eq.(11) *
            QA=QX1(0)+QX2(0)+QXS(0)
              QX(0)=QA*CR
              QY(0)=QA*SR
              QZ(0)=-(QZ1(0)+QZ2(0)+QZS(0))
C
c         * Utx, Uty     : see Hisada(1994) eq.(11) *
            VA=VX1(0)+VX2(0)+VXS(0)
              VX(0)=-VA*SR
              VY(0)= VA*CR
C
c         * Uzx, Uzy, Uzz: see Hisada(1994) eq.(11) *
            WA=WX1(0)+WX2(0)+WXS(0)
              WX(0)=-WA*CR
              WY(0)=-WA*SR
              WZ(0)=-(WZ1(0)+WZ2(0)+WZS(0))
c
c         * Urx,x, Ury,x, Urx,y, Ury,y, Urz,x Urz,y: see Hisada(1994) eq.(27) *
            QB=QX1(1)+QX2(1)+QXS(1)
              QX(1)=QB*C2+QA*CX
              QY(1)=QB*CS+QA*SX
              QX(2)=QB*CS+QA*CY
              QY(2)=QB*S2+QA*SY
            QB=QZ1(1)+QZ2(1)+QZS(1)
              QZ(1)=-QB*CR
              QZ(2)=-QB*SR
c
c         * Utx,x, Uty,x, Utx,y, Uty,y: see Hisada(1994) eq.(27) *
            VB=VX1(1)+VX2(1)+VXS(1)
              VX(1)=-(VB*CS+VA*SX)
              VY(1)=  VB*C2+VA*CX
              VX(2)=-(VB*S2+VA*SY)
              VY(2)=  VB*CS+VA*CY
c
c         * Uzx,x, Uzy,x, Uzx,y, Uzy,y, Uzz,x Uzz,y: see Hisada(1994) eq.(27) *
            WB=WX1(1)+WX2(1)+WXS(1)
              WX(1)=-(WB*C2+WA*CX)
              WY(1)=-(WB*CS+WA*SX)
              WX(2)=-(WB*CS+WA*CY)
              WY(2)=-(WB*S2+WA*SY)
            WB=WZ1(1)+WZ2(1)+WZS(1)
              WZ(1)=-WB*CR
              WZ(2)=-WB*SR
c
c         * Urx,z, Ury,z, Urz,z: see Hisada(1994) eq.(29) *
            QB=QX1(2)+QX2(2)+QXS(2)
              QX(3)= QB*CR
              QY(3)= QB*SR
              QZ(3)=-(QZ1(2)+QZ2(2)+QZS(2))
c
c         * Utx,z, Uty,z: see Hisada(1994) eq.(29) *
            VB=VX1(2)+VX2(2)+VXS(2)
              VX(3)=-VB*SR
              VY(3)= VB*CR
c
c         * Uzx,z, Uzy,z, Uzz,z: see Hisada(1994) eq.(29) *
            WB=WX1(2)+WX2(2)+WXS(2)
              WX(3)=-WB*CR
              WY(3)=-WB*SR
              WZ(3)=-(WZ1(2)+WZ2(2)+WZS(2))
C
c         * Uxx,x, Uxx,y, Uxy,x, Uxy,y, Uxz,x, Uxz,y: see Hisada(1994) eq.(24) *
            UX(1,1)=QX(1)*CR+QX(0)*CX-VX(1)*SR-VX(0)*SX
            UX(1,2)=QX(2)*CR+QX(0)*CY-VX(2)*SR-VX(0)*SY
            UY(1,1)=QY(1)*CR+QY(0)*CX-VY(1)*SR-VY(0)*SX
            UY(1,2)=QY(2)*CR+QY(0)*CY-VY(2)*SR-VY(0)*SY
            UZ(1,1)=QZ(1)*CR+QZ(0)*CX
            UZ(1,2)=QZ(2)*CR+QZ(0)*CY
C
c         * Uyx,x, Uyx,y, Uyy,x, Uyy,y, Uyz,x, Uyz,y: see Hisada(1994) eq.(24) *
            UX(2,1)=QX(1)*SR+QX(0)*SX+VX(1)*CR+VX(0)*CX
            UX(2,2)=QX(2)*SR+QX(0)*SY+VX(2)*CR+VX(0)*CY
            UY(2,1)=QY(1)*SR+QY(0)*SX+VY(1)*CR+VY(0)*CX
            UY(2,2)=QY(2)*SR+QY(0)*SY+VY(2)*CR+VY(0)*CY
            UZ(2,1)=QZ(1)*SR+QZ(0)*SX
            UZ(2,2)=QZ(2)*SR+QZ(0)*SY
C
c         * Uzx,x, Uzx,y, Uzy,x, Uzy,y, Uzz,x, Uzz,y: see Hisada(1994) eq.(27) *
            UX(3,1)=WX(1)
            UX(3,2)=WX(2)
            UY(3,1)=WY(1)
            UY(3,2)=WY(2)
            UZ(3,1)=WZ(1)
            UZ(3,2)=WZ(2)
C
c         * Uxx,z, Uxy,z, Uxz,z, Uyx,z, Uyy,z, Uyz,z: see Hisada(1994) eq.(25) *
            UX(1,3)=QX(3)*CR-VX(3)*SR
            UY(1,3)=QY(3)*CR-VY(3)*SR
            UZ(1,3)=QZ(3)*CR
            UX(2,3)=QX(3)*SR+VX(3)*CR
            UY(2,3)=QY(3)*SR+VY(3)*CR
            UZ(2,3)=QZ(3)*SR
C
c         * Uzx,z, Uzy,z, Uzz,z: see Hisada(1994) eq.(29) *
            UX(3,3)=WX(3)
            UY(3,3)=WY(3)
            UZ(3,3)=WZ(3)
C
c         * Gk=Mij*Ujk,i: see Hisada(1994) eq.(22) *
c
            GX=(0.D0,0.D0)
            GY=(0.D0,0.D0)
            GZ=(0.D0,0.D0)
            DO 128 I=1,3
            DO 129 J=1,3
              GX=GX+EN(I,J)*UX(J,I)
              GY=GY+EN(I,J)*UY(J,I)
              GZ=GZ+EN(I,J)*UZ(J,I)
  129       CONTINUE
  128       CONTINUE
            GX=GX*SMO
            GY=GY*SMO
            GZ=GZ*SMO
C
            GGX(IO,IOB)=GGX(IO,IOB)+GX
            GGY(IO,IOB)=GGY(IO,IOB)+GY
            GGZ(IO,IOB)=GGZ(IO,IOB)+GZ
C
  120     CONTINUE
C
  110   CONTINUE
C
  100 CONTINUE
C
      IF(IOT.eq.1)THEN
      CLOSE(10)
      CLOSE(11)
      CLOSE(12)
      CLOSE(13)
      CLOSE(14)
      CLOSE(15)
      CLOSE(16)
      CLOSE(17)
      CLOSE(18)
      CLOSE(19)
      CLOSE(20)
      CLOSE(21)
      CLOSE(22)
      CLOSE(23)
      CLOSE(24)
      END IF
C
      RETURN
      END
C
C
C
C *** EN AND SMO (MOMENT TENSOR FOR THE SOURCE) ***
C
      SUBROUTINE ENMO(OM,EN,SMO)
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION EN(3,3)
      COMPLEX*16 VS,VP,SMO,CPX,SM,RAM,STT,STL
      COMMON/STRC/DNS(20),VP(20),VS(20),THK(20),NL
     *    /FLT/XF(30),YF(30),LF(30),DP(30),STR(30),DIP(30),RAK(30)
     *      ,FL(30),FW(30),FD(30),RT(30),VR(30),NVR(30),TD(30),NF,IS 
     *    /FCQVW/SM,RAM,RR
C
C *** CHECK OF DIMENSION ***
C
      NC=20
      IF(NL.GE.NC)THEN
        WRITE(*,*)' FROM ENMO:NL>NC !!!!!!!!!!'
        WRITE(*,*)' NL,NC=',NL,NC
        STOP
      END IF
C
      NC=30
      IF(NF.GE.NC)THEN
        WRITE(*,*)' FROM ENMO:NF>NC !!!!!!!!!!'
        WRITE(*,*)' NF,NC=',NF,NC
        STOP
      END IF
C
      CPX=(0.D0,1.D0)
C
C *** EN ***
C
      CST= DCOS(STR(IS))
      SST= DSIN(STR(IS))
C
      CDP= DCOS(DIP(IS))
      SDP= DSIN(DIP(IS))
C
      CRK= DCOS(RAK(IS))
      SRK= DSIN(RAK(IS))
C
      XN=-SDP*SST
      YN= SDP*CST
      ZN=-CDP
C
      EX= CRK*CST+SRK*CDP*SST
      EY= CRK*SST-SRK*CDP*CST
      EZ=-SRK*SDP
C
      EN(1,1)= 2.D0*EX*XN
      EN(1,2)= EX*YN +EY*XN
      EN(1,3)= EX*ZN +EZ*XN
      EN(2,1)= EY*XN +EX*YN
      EN(2,2)= 2.D0*EY*YN
      EN(2,3)= EY*ZN +EZ*YN
      EN(3,1)= EZ*XN +EX*ZN
      EN(3,2)= EZ*YN +EY*ZN
      EN(3,3)= 2.D0*EZ*ZN
C
C *** SMO ***
C
      L=LF(IS)
      SM = DNS(L)*VS(L)*VS(L)
      RAM= DNS(L)*VP(L)*VP(L)-2.D0*SM
C
      SMO= FD(IS)*SM*FL(IS)*FW(IS)
C
      IF(NVR(IS).EQ.1)THEN
        STL= (1.D0,0.D0)
      ELSE
        TL=FL(IS)/VR(IS)
        TTL= OM*TL/2.D0
        STL= DSIN(TTL)*CDEXP(CPX*TTL)/TTL
      END IF
C
      TTK= OM*RT(IS)/2.D0
      IF(TTK.EQ.0.D0)THEN
        STT= (1.D0,0.D0)
      ELSE
        STT= DSIN(TTK)*CDEXP(CPX*TTK)/TTK
      END IF
C
      SMO=CPX*SMO*STL*STT/OM
C
      RETURN
      END
C
C
C *** THE COMPUTATION OF Q, V, AND W AND THEIR DERIVATIVES 
C     USING SIMPSON'S METHOD ***
C
      SUBROUTINE QVWRZ(OM,A,B,N,GQX,GQZ,GVX,GWX,GWZ)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 GQX(0:2),GQZ(0:2),GVX(0:2),GWX(0:2),GWZ(0:2)
     *          ,QXA(0:2),QZA(0:2),VXA(0:2),WXA(0:2),WZA(0:2)
     *          ,QXB(0:2),QZB(0:2),VXB(0:2),WXB(0:2),WZB(0:2)
     *          ,QXC(0:2),QZC(0:2),VXC(0:2),WXC(0:2),WZC(0:2)
     *          ,VP,VS,VP2,VS2,SM1,X1,SM2,X2,SM3,X3,RM1,RM2,RM3
      COMMON/OTST/IOT,IST /SOLY/SDEP,ODEP,LS,LO
     *      /STRC/DNS(20),VP(20),VS(20),THK(20),NL
     *      /ASYMP/SM1,X1,SM2,X2,SM3,X3 /RAM/RM1,RM2,RM3
     *      /DPSO /DEP1,DEP2,SDP,ODP
C
C *** CHECK OF DIMENSION ***
C
      NC=20
      IF(NL.GE.NC)THEN
        WRITE(*,*)' FROM QVWRZ:NL>NC !!!!!!!!!!'
        WRITE(*,*)' NL,NC=',NL,NC
        STOP
      END IF
C
      IF(N.LT.1)THEN
        WRITE(*,*)' FROM QVWRZ: N.LT.1 !!!!!!!'
        RETURN
      END IF
      IF(A.EQ.B)THEN
        WRITE(*,*)' FROM QVWRZ: A.EQ.B !!!!!!!'
        RETURN
      END IF
C
      DO 95 J=0,2
      GQX(J)=(0.D0,0.D0)
      GQZ(J)=(0.D0,0.D0)
      GVX(J)=(0.D0,0.D0)
      GWX(J)=(0.D0,0.D0)
      GWZ(J)=(0.D0,0.D0)
   95 CONTINUE
C
      IASYM=0
      IF((IST.EQ.1).AND.(IABS(LS-LO).LE.1))THEN
C
        IASYM=1
C
        SDP=0.D0
        DEP1=0.D0
        DEP2=0.D0
        DO 200 I=1,LS-1
        DEP1=DEP1+THK(I)
  200   CONTINUE
        SDP=DEP1+SDEP
        DEP2=DEP1+THK(LS)
C
        ODP=0.D0
        DO 210 I=1,LO-1
        ODP=ODP+THK(I)
  210   CONTINUE
        ODP=ODP+ODEP
C         
        VP2=VP(LS)*VP(LS)
        VS2=VS(LS)*VS(LS)
        SM2=DNS(LS)*VS2
        RM2=DNS(LS)*VP(LS)*VP(LS)-2.D0*SM2
        X2=(1.D0+VS2/VP2)/(1.D0-VS2/VP2)
C
        IF(LS.GE.2)THEN
          VP2=VP(LS-1)*VP(LS-1)
          VS2=VS(LS-1)*VS(LS-1)
          SM1=DNS(LS-1)*VS2
          RM1=DNS(LS-1)*VP(LS-1)*VP(LS-1)-2.D0*SM1
          X1=(1.D0+VS2/VP2)/(1.D0-VS2/VP2)
        END IF
C
        IF(LS.NE.NL)THEN
          VP2=VP(LS+1)*VP(LS+1)
          VS2=VS(LS+1)*VS(LS+1)
          SM3=DNS(LS+1)*VS2
          RM3=DNS(LS+1)*VP(LS+1)*VP(LS+1)-2.D0*SM3
          X3=(1.D0+VS2/VP2)/(1.D0-VS2/VP2)
        END IF
C
        CALL CASYM(LS,NL)
C
      END IF
C
      H=(B-A)/N
      CALL FQVWRZ(IASYM,NL,OM,A,QXA,QZA,VXA,WXA,WZA)
      DO  300 I=1,N
        X=A+H*(FLOAT(I)-0.5D0)
        CALL FQVWRZ(IASYM,NL,OM,X,QXB,QZB,VXB,WXB,WZB)
        X=A+H*FLOAT(I)
        CALL FQVWRZ(IASYM,NL,OM,X,QXC,QZC,VXC,WXC,WZC)
        DO 310 J=0,2
         GQX(J)= GQX(J)+H*(QXA(J)+4.D0*QXB(J)+QXC(J))/6.D0
         GVX(J)= GVX(J)+H*(VXA(J)+4.D0*VXB(J)+VXC(J))/6.D0
         GWX(J)= GWX(J)+H*(WXA(J)+4.D0*WXB(J)+WXC(J))/6.D0
         GQZ(J)= GQZ(J)+H*(QZA(J)+4.D0*QZB(J)+QZC(J))/6.D0
         GWZ(J)= GWZ(J)+H*(WZA(J)+4.D0*WZB(J)+WZC(J))/6.D0
         QXA(J)= QXC(J)
         VXA(J)= VXC(J)
         WXA(J)= WXC(J)
         QZA(J)= QZC(J)
         WZA(J)= WZC(J)
  310   CONTINUE
  300 CONTINUE
C
      RETURN
      END
C
C
C
C *** THE INTEGRANDS OF Q, V, AND W AND THEIR DERIVATIVES ***
C
      SUBROUTINE FQVWRZ(IASYM,NL,OM,WN,QX,QZ,VX,WX,WZ)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 QX(0:2),QZ(0:2),VX(0:2),WX(0:2),WZ(0:2)
     *   ,VLX(2),URX(2,2),URZ(2,2),VLXZ,URXZ(2),URZZ(2),SM,RAM
     *   ,H1X,H2X,V1X,V2X,V3X,V4X,V1Z,V2Z,V3Z,V4Z
      DIMENSION BJ0(0:1),BJ1(0:1),BJK(0:1),BJD(0:1)
      COMMON /OTST/IOT,IST /FCQVW/SM,RAM,RR /SOLY/SDEP,ODEP,LS,LO
C
      CALL RTSH(OM,WN,VLX)
      CALL RTPSV(OM,WN,URX,URZ)
C
      IF(IASYM.EQ.1)THEN
        CALL ASYM(LS,NL,WN,H1X,H2X,V1X,V2X,V3X,V4X,V1Z,V2Z,V3Z,V4Z)
C
        VLX(1)=VLX(1)-H1X
        URX(1,1)=URX(1,1)-V1X
        URX(2,1)=URX(2,1)-V2X
        URZ(1,1)=URZ(1,1)-V1Z
        URZ(2,1)=URZ(2,1)-V2Z
        VLX(2)=VLX(2)-H2X
        URX(1,2)=URX(1,2)-V3X
        URX(2,2)=URX(2,2)-V4X
        URZ(1,2)=URZ(1,2)-V3Z
        URZ(2,2)=URZ(2,2)-V4Z

      END IF
C
      VLXZ=VLX(2)/SM
      URXZ(1)=URX(1,2)/SM+WN*URX(2,1)
      URZZ(1)=URZ(1,2)/SM+WN*URZ(2,1)
      URXZ(2)=(URX(2,2)-WN*RAM*URX(1,1))/(RAM+2.D0*SM)
      URZZ(2)=(URZ(2,2)-WN*RAM*URZ(1,1))/(RAM+2.D0*SM)
C
      RK=WN*RR
      BJ0(0)=BESJ0(RK)
      BJ1(0)=BESJ1(RK)
      BJK(0)=BJ1(0)/RK
      BJD(0)=BJ0(0)-BJK(0)
C
      BJ0(1)=-BJ1(0)*WN
      BJ1(1)= BJD(0)*WN
      BJK(1)=(BJ1(1)-BJ1(0)/RR)/RK
      BJD(1)= BJ0(1)-BJK(1)
C
      QX(0)= URX(1,1)*BJD(0)+VLX(1)*BJK(0)
      VX(0)= URX(1,1)*BJK(0)+VLX(1)*BJD(0)
      WX(0)= URX(2,1)*BJ1(0)
      QZ(0)= URZ(1,1)*BJ1(0)
      WZ(0)= URZ(2,1)*BJ0(0)
C
      QX(1)= URX(1,1)*BJD(1)+VLX(1)*BJK(1)
      VX(1)= URX(1,1)*BJK(1)+VLX(1)*BJD(1)
      WX(1)= URX(2,1)*BJ1(1)
      QZ(1)= URZ(1,1)*BJ1(1)
      WZ(1)= URZ(2,1)*BJ0(1)
C
      QX(2)= URXZ(1)*BJD(0)+VLXZ*BJK(0)
      VX(2)= URXZ(1)*BJK(0)+VLXZ*BJD(0)
      WX(2)= URXZ(2)*BJ1(0)
      QZ(2)= URZZ(1)*BJ1(0)
      WZ(2)= URZZ(2)*BJ0(0)
C
      IF(IOT.eq.1)THEN
        WRITE(10,1100)WN,QX(0)
        WRITE(11,1100)WN,VX(0)
        WRITE(12,1100)WN,WX(0)
        WRITE(13,1100)WN,QZ(0)
        WRITE(14,1100)WN,WZ(0)
        WRITE(15,1100)WN,QX(1)
        WRITE(16,1100)WN,VX(1)
        WRITE(17,1100)WN,WX(1)
        WRITE(18,1100)WN,QZ(1)
        WRITE(19,1100)WN,WZ(1)
        WRITE(20,1100)WN,QX(2)
        WRITE(21,1100)WN,VX(2)
        WRITE(22,1100)WN,WX(2)
        WRITE(23,1100)WN,QZ(2)
        WRITE(24,1100)WN,WZ(2)
      END IF
 1100 FORMAT(1X,3D14.5)
C
      RETURN
      END
C
C
C *** COEFFICIENTS FOR ASYMPTOTIC SOLUTIONS ***
C
      SUBROUTINE CASYM(LS,NL)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 SM1,X1,SM2,X2,SM3,X3,BB(4,4),BB1(4,4),BB2(4,4)
     *   ,H1T,H2T,H1D,H2D,RH,H1R1,H2R1,H1R2,H2R2       
     *   ,CUPX,CUSX,CUP1X,CUS1X,CUP2X,CUS2X
     *   ,CDPX,CDSX,CDP1X,CDS1X,CDP2X,CDS2X
     *   ,CUPZ,CUSZ,CUP1Z,CUS1Z,CUP2Z,CUS2Z
     *   ,CDPZ,CDSZ,CDP1Z,CDS1Z,CDP2Z,CDS2Z
      COMMON/ASYMP/SM1,X1,SM2,X2,SM3,X3
     *      /DPSO /DEP1,DEP2,SDP,ODP
     *      /HX  /H1T,H2T,H1D,H2D,H1R1,H2R1,H1R2,H2R2       
     *      /VUX /CUPX,CUSX,CUP1X,CUS1X,CUP2X,CUS2X
     *      /VDX /CDPX,CDSX,CDP1X,CDS1X,CDP2X,CDS2X
     *      /VUZ /CUPZ,CUSZ,CUP1Z,CUS1Z,CUP2Z,CUS2Z
     *      /VDZ /CDPZ,CDSZ,CDP1Z,CDS1Z,CDP2Z,CDS2Z
C
C *** ODP < DEP1 ***
C
      IF(ODP.GE.DEP1) GOTO 10
      IF(LS.EQ.1) GOTO 10
C
C  ** QX **
C
C   * SH (TRANSMITTED) *
C
      H1T=0.15915494309D0/(SM1+SM2)
      H2T=SM1*H1T
C
C   * P-SV (TRANSMITTED) *
C
      CUPX=0.25d0/(3.141592654d0*SM2*(X2+1.0d0))
      CUSX=CUPX*(X2-1.0)
C
      CALL EINVE(SM1,SM2,X1,X2,BB)
C
      CUP1X=BB(3,3)*CUPX+BB(3,4)*CUSX
      CUS1X=BB(4,3)*CUPX+BB(4,4)*CUSX
      CUP2X=BB(3,4)*CUPX
      CUS2X=BB(4,4)*CUPX
C
C  ** QZ **
C
C   * P-SV (TRANSMITTED) *
C
      CUSZ=0.25/(3.141592654d0*SM2*(X2+1.0d0))
      CUPZ=-CUSZ
C
      CUP1Z=BB(3,3)*CUPZ+BB(3,4)*CUSZ
      CUS1Z=BB(4,3)*CUPZ+BB(4,4)*CUSZ
      CUP2Z=BB(3,4)*CUPZ
      CUS2Z=BB(4,4)*CUPZ
C
      GOTO 50
C
C *** DEP1 < (OR =) ODP < SDP ***
C
   10 IF(ODP.GE.SDP) GOTO 20
C
C  ** QX **
C
C   * SH *
C
      H2D=0.07957747155D0
      H1D=H2D/SM2
      RH=(SM2-SM1)/(SM2+SM1)
      IF(LS.EQ.1)RH=(1.D0,0.D0)
      H1R1= RH*H1D
      H2R1=-RH*H2D
      RH=(SM3-SM2)/(SM3+SM2)
      IF(LS.EQ.NL)RH=(0.D0,0.D0)
      H1R2=-RH*H1D
      H2R2=-RH*H2D
C
C   * P-SV (UP-GOING DIRECT) *
C
      CUPX=0.25d0/(3.141592654d0*SM2*(X2+1.0d0))
      CUSX=CUPX*(X2-1.0)
C
C   * P-SV (DOWN-GOING REFLECTED FROM UPPER BOUNDARY) *
C
      IF(LS.EQ.1)THEN
        BB1(1,3)=-(X2-2.0)
        BB1(1,4)=2.D0
        BB1(2,3)=-0.5D0*(X2-1.D0)*(X2-3.D0)
        BB1(2,4)=X2-2.D0
      ELSE 
        CALL EINVE(SM1,SM2,X1,X2,BB1)
      END IF
C
      CDP1X=BB1(1,3)*CUPX+BB1(1,4)*CUSX
      CDS1X=BB1(2,3)*CUPX+BB1(2,4)*CUSX
      CDP2X=BB1(1,4)*CUPX
      CDS2X=BB1(2,4)*CUPX
C
C   * P-SV (UP-GOING REFLECTED FROM LOWER BOUNDARY) *
C
      IF(LS.EQ.NL)GOTO 15
C
      CDPX=CUPX
      CDSX=CUSX
C
      CALL EINVE(SM2,SM3,X2,X3,BB2)
C
      CUP1X=BB2(3,1)*CDPX+BB2(3,2)*CDSX
      CUS1X=BB2(4,1)*CDPX+BB2(4,2)*CDSX
      CUP2X=BB2(3,2)*CDPX
      CUS2X=BB2(4,2)*CDPX
C
C  ** QZ **
C
   15 CUSZ=0.25d0/(3.141592654d0*SM2*(X2+1.0d0))
      CUPZ=-CUSZ
C
      CDP1Z=BB1(1,3)*CUPZ+BB1(1,4)*CUSZ
      CDS1Z=BB1(2,3)*CUPZ+BB1(2,4)*CUSZ
      CDP2Z=BB1(1,4)*CUPZ
      CDS2Z=BB1(2,4)*CUPZ
C
      IF(LS.EQ.NL)GOTO 50
C
      CDPZ=CUSZ
      CDSZ=-CDPZ
C
      CUP1Z=BB2(3,1)*CDPZ+BB2(3,2)*CDSZ
      CUS1Z=BB2(4,1)*CDPZ+BB2(4,2)*CDSZ
      CUP2Z=BB2(3,2)*CDPZ
      CUS2Z=BB2(4,2)*CDPZ
C
      GOTO 50
C
C *** SDP < (OR =) ODP < DEP2 ***
C
   20 IF((LS.NE.NL).AND.(ODP.GE.DEP2)) GOTO 30
C
C  ** QX **
C
C   * SH *
C
      H2D=0.07957747155D0
      H1D=H2D/SM2
      RH=(SM2-SM1)/(SM2+SM1)
      IF(LS.EQ.1)RH=(1.D0,0.D0)
      H1R1= RH*H1D
      H2R1=-RH*H2D
      RH=(SM3-SM2)/(SM3+SM2)
      IF(LS.EQ.NL)RH=(0.D0,0.D0)
      H1R2=-RH*H1D
      H2R2=-RH*H2D
C
C   * P-SV (DOWN-GOING DIRECT) *
C
      CDPX=0.25d0/(3.141592654d0*SM2*(X2+1.0d0))
      CDSX=CDPX*(X2-1.0d0)
C
C   * P-SV (UP-GOING REFLECTED FROM LOWER BOUNDARY) *
C
      IF(LS.EQ.NL)GOTO 25
C
      CALL EINVE(SM2,SM3,X2,X3,BB2)
C
      CUP1X=BB2(3,1)*CDPX+BB2(3,2)*CDSX
      CUS1X=BB2(4,1)*CDPX+BB2(4,2)*CDSX
      CUP2X=BB2(3,2)*CDPX
      CUS2X=BB2(4,2)*CDPX
C
C   * P-SV (DOWN-GOING REFLECTED FROM UPPER BOUNDARY) *
C
   25 CUPX=CDPX
      CUSX=CDSX
C
      IF(LS.EQ.1)THEN
        BB1(1,3)=-(X2-2.0)
        BB1(1,4)=2.D0
        BB1(2,3)=-0.5D0*(X2-1.D0)*(X2-3.D0)
        BB1(2,4)=X2-2.D0
      ELSE 
        CALL EINVE(SM1,SM2,X1,X2,BB1)
      END IF
C
      CDP1X=BB1(1,3)*CUPX+BB1(1,4)*CUSX
      CDS1X=BB1(2,3)*CUPX+BB1(2,4)*CUSX
      CDP2X=BB1(1,4)*CUPX
      CDS2X=BB1(2,4)*CUPX
C
C  ** QZ **
C
      CDPZ=0.25d0/(3.141592654d0*SM2*(X2+1.0d0))
      CDSZ=-CDPZ
C
      IF(LS.EQ.NL)GOTO 27
C
      CUP1Z=BB2(3,1)*CDPZ+BB2(3,2)*CDSZ
      CUS1Z=BB2(4,1)*CDPZ+BB2(4,2)*CDSZ
      CUP2Z=BB2(3,2)*CDPZ
      CUS2Z=BB2(4,2)*CDPZ
C
   27 CUPZ=-CDPZ
      CUSZ=-CDSZ
C
      CDP1Z=BB1(1,3)*CUPZ+BB1(1,4)*CUSZ
      CDS1Z=BB1(2,3)*CUPZ+BB1(2,4)*CUSZ
      CDP2Z=BB1(1,4)*CUPZ
      CDS2Z=BB1(2,4)*CUPZ
C
      GOTO 50
C
C *** SDP < DEP2 < (OR =) ODP ***
C
C  ** QX **
C
C   * SH (TRANSMITTED) *
C
   30 H1T=0.15915494309D0/(SM3+SM2)
      H2T=-SM3*H1T
C
C   * P-SV (TRANSMITTED) *
C
      CDPX=0.25d0/(3.141592654d0*SM2*(X2+1.0d0))
      CDSX=CDPX*(X2-1.0d0)
C
      CALL EINVE(SM2,SM3,X2,X3,BB)
C
      CDP1X=BB(1,1)*CDPX+BB(1,2)*CDSX
      CDS1X=BB(2,1)*CDPX+BB(2,2)*CDSX
      CDP2X=BB(1,2)*CDPX
      CDS2X=BB(2,2)*CDPX
C
C  ** QZ **
C
      CDPZ=0.25d0/(3.141592654d0*SM2*(X2+1.0d0))
      CDSZ=-CDPZ
C
      CALL EINVE(SM2,SM3,X2,X3,BB)
C
      CDP1Z=BB(1,1)*CDPZ+BB(1,2)*CDSZ
      CDS1Z=BB(2,1)*CDPZ+BB(2,2)*CDSZ
      CDP2Z=BB(1,2)*CDPZ
      CDS2Z=BB(2,2)*CDPZ
C
   50 RETURN
C
      END
C
C
C *** PRODUCT OF INV E(L+1) * E(L) (SEE APPENDIX B IN HISADA, 1995)***
C
      SUBROUTINE EINVE(SM1,SM2,X1,X2,BB)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 BB(4,4),SM1,SM2,X1,X2,A,A2,B
C
      A=SM1/SM2
      A2=A*A
      B=0.5/(X1+A2*X2+A*(X1*X2+1.0))
C
      BB(1,1)=2.0*(X1+1.0)*A*(A+X1)
      BB(2,1)=(X1+1.0)*A*(X1*X2-3.0*X1+2.0-A*(X1*X2-3.0*X2+2.0))
      BB(3,1)=2.0*(X1-2.0)*(1.0-A)*(1.0+A*X2)
      BB(4,1)=A*(X2-1.0)*(X1*X1-3.0*X1+4.0)
     *         -(X1-1.0)*(4.0+A2*(X1-3.0)*X2)
C
      BB(1,2)=(0.0,0.0)
      BB(2,2)=2.0*(X1+1.0)*A*(1.0+A*X2)
      BB(3,2)=4.0*(A-1.0)*(1.0+A*X2)
      BB(4,2)=-BB(3,1)
C
      BB(1,3)=2.0*(X2-2.0)*(A-1.0)*(A+X1)
      BB(2,3)=A*(X1-1.0)*(X2*X2-3.0*X2+4.0)
     *         -(X2-1.0)*(4.0*A2+(X2-3.0)*X1)
      BB(3,3)=2.0*(X2+1.0)*(1.0+A*X2)
      BB(4,3)=(X2+1.0)*(A*(X1*X2-3.0*X2+2.0)-X1*X2+3.0*X1-2.0)
C
      BB(1,4)=4.0*(1.0-A)*(A+X1)
      BB(2,4)=-BB(1,3)
      BB(3,4)=(0.0,0.0)
      BB(4,4)=2.0*(X2+1.0)*(A+X1)
C
      DO 100 J=1,4
      DO 100 I=1,4
      BB(I,J)=BB(I,J)*B
  100 CONTINUE
C
      RETURN
      END
C
C
C *** ASTMPTOTIC SOLUTION ***
C
      SUBROUTINE ASYM(LS,NL,WN,H1X,H2X
     *               ,V1X,V2X,V3X,V4X,V1Z,V2Z,V3Z,V4Z)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 SM1,SM2,SM3,X1,X2,X3,A,B
     *   ,H1T,H2T,H1D,H2D,H1R1,H2R1,H1R2,H2R2       
     *   ,CUPX,CUSX,CUP1X,CUS1X,CUP2X,CUS2X
     *   ,CDPX,CDSX,CDP1X,CDS1X,CDP2X,CDS2X
     *   ,CUPZ,CUSZ,CUP1Z,CUS1Z,CUP2Z,CUS2Z
     *   ,CDPZ,CDSZ,CDP1Z,CDS1Z,CDP2Z,CDS2Z
     *   ,CSD,CPD,CST,CPT,CSR1,CPR1,CSR2,CPR2
     *   ,V1X,V2X,V3X,V4X,V1Z,V2Z,V3Z,V4Z,H1X,H2X
      COMMON /ASYMP/SM1,X1,SM2,X2,SM3,X3
     *       /DPSO /DEP1,DEP2,SDP,ODP
     *      /HX  /H1T,H2T,H1D,H2D,H1R1,H2R1,H1R2,H2R2       
     *      /VUX /CUPX,CUSX,CUP1X,CUS1X,CUP2X,CUS2X
     *      /VDX /CDPX,CDSX,CDP1X,CDS1X,CDP2X,CDS2X
     *      /VUZ /CUPZ,CUSZ,CUP1Z,CUS1Z,CUP2Z,CUS2Z
     *      /VDZ /CDPZ,CDSZ,CDP1Z,CDS1Z,CDP2Z,CDS2Z
C
C *** DISP-STRESS VECTOR DUE TO Qx FOR ODP < DEP1 ***
C
      IF(ODP.GE.DEP1)GOTO 10
      IF(LS.EQ.1)GOTO 10
C
C   * UP-GOING WAVE DUE TO TRANSMITTED WAVE *
C
      WHZ=WN*(SDP-ODP)
      WHZ1=WN*(SDP-DEP1)
      WZ1Z=WN*(DEP1-ODP)
      EX=DEXP(-WHZ)
C
      CPT=EX*(CUP1X-WHZ1*CUP2X)
      CST=EX*(CUS1X-WZ1Z*CUP1X-WHZ1*CUS2X+WZ1Z*WHZ1*CUP2X)
C
      H1X=EX*H1T
      H2X=WN*EX*H2T
C
      A=X1-1.0
      B=X1-3.0
      V1X=CPT+CST
      V2X=A*CPT-CST
      V3X=WN*SM1*(-B*CPT+2.0*CST)
      V4X=WN*SM1*( A*CPT-2.0*CST)
C
C   * DISP-STRESS VECTOR DUE TO Qz *
C
      CPT=EX*(CUP1Z-WHZ1*CUP2Z)
      CST=EX*(CUS1Z-WZ1Z*CUP1Z-WHZ1*CUS2Z+WZ1Z*WHZ1*CUP2Z)
C
      V1Z=CPT+CST
      V2Z=A*CPT-CST
      V3Z=WN*SM1*(-B*CPT+2.0*CST)
      V4Z=WN*SM1*( A*CPT-2.0*CST)
C
      GOTO 50
C
C *** DISP-STRESS VECTOR DUE TO Qx FOR DEP1 < (OR =) ODP < SDP ***
C
   10 IF(ODP.GE.SDP)GOTO 20
C
C   * UP-GOING WAVE DUE TO DIRECT WAVE *
C
      WHZ=WN*(SDP-ODP)
      EXD=DEXP(-WHZ)
C
      CPD=EXD*CUPX
      CSD=EXD*(CUSX-WHZ*CUPX)
C
C   * DOWN-GOING WAVE DUE TO REFLECT WAVE FROM UPPER BOUNDARY *
C
      WHZ1=WN*(SDP-DEP1)
      WZZ1=WN*(ODP-DEP1)
      EXR1=DEXP(-WHZ1-WZZ1)
C
      CPR1=EXR1*(CDP1X-WHZ1*CDP2X)
      CSR1=EXR1*(CDS1X-WHZ1*CDS2X-WZZ1*CDP1X+WZZ1*WHZ1*CDP2X)
C
C   * UP-GOING WAVE DUE TO REFLECT WAVE FROM LOWER BOUNDARY *
C
      IF(LS.EQ.NL)THEN
        CPR2=(0.D0,0.D0)
        CSR2=(0.D0,0.D0)
      ELSE
        WZ2H=WN*(DEP2-SDP)
        WZ2Z=WN*(DEP2-ODP)
        EXR2=DEXP(-WZ2H-WZ2Z)
        CPR2=EXR2*(CUP1X-WZ2H*CUP2X)
        CSR2=EXR2*(CUS1X-WZ2H*CUS2X-WZ2Z*CUP1X+WZ2Z*WZ2H*CUP2X)
      END IF
C
C   * DISP-STRESS VECTORS *
C
      H1X=H1D*EXD+H1R1*EXR1+H1R2*EXR2
      H2X=WN*(H2D*EXD+H2R1*EXR1+H2R2*EXR2)
C
      A=X2-1.0
      B=X2-3.0
      V1X=CPR1+CSR1+CPD+CSD+CPR2+CSR2
      V2X=-A*CPR1+CSR1+A*(CPD+CPR2)-CSD-CSR2
      V3X=WN*SM2*(B*CPR1-2.0*CSR1-B*(CPD+CPR2)+2.0*(CSD+CSR2))
      V4X=WN*SM2*(A*CPR1-2.0*CSR1+A*(CPD+CPR2)-2.0*(CSD+CSR2))
C
C   * DISP-STRESS VECTOR DUE TO Qz *
C
      CPD=EXD*CUPZ
      CSD=EXD*(CUSZ-WHZ*CUPZ)
C
      CPR1=EXR1*(CDP1Z-WHZ1*CDP2Z)
      CSR1=EXR1*(CDS1Z-WHZ1*CDS2Z-WZZ1*CDP1Z+WZZ1*WHZ1*CDP2Z)
C
      IF(LS.EQ.NL)THEN
        CPR2=(0.D0,0.D0)
        CSR2=(0.D0,0.D0)
      ELSE
        CPR2=EXR2*(CUP1Z-WZ2H*CUP2Z)
        CSR2=EXR2*(CUS1Z-WZ2H*CUS2Z-WZ2Z*CUP1Z+WZ2Z*WZ2H*CUP2Z)
      END IF
C
      V1Z=CPR1+CSR1+CPD+CSD+CPR2+CSR2
      V2Z=-A*CPR1+CSR1+A*(CPD+CPR2)-CSD-CSR2
      V3Z=WN*SM2*(B*CPR1-2.0*CSR1-B*(CPD+CPR2)+2.0*(CSD+CSR2))
      V4Z=WN*SM2*(A*CPR1-2.0*CSR1+A*(CPD+CPR2)-2.0*(CSD+CSR2))
C
      GOTO 50
C
C *** DISP-STRESS VECTOR DUE TO Qx FOR SDP < (OR =) ODP < DEP2 ***
C
   20 IF((LS.NE.NL).AND.(ODP.GE.DEP2))GOTO 30
C
C   * DOWN-GOING WAVE DUE TO DIRECT WAVE *
C
      WZH=WN*(ODP-SDP)
      EXD=DEXP(-WZH)
C
      CPD=EXD*CDPX
      CSD=EXD*(CDSX-WZH*CDPX)
C
C   * UP-GOING WAVE DUE TO REFLECT WAVE FROM LOWER BOUNDARY *
C
      IF(LS.EQ.NL)THEN
        CPR2=(0.D0,0.D0)
        CSR2=(0.D0,0.D0)
      ELSE
        WZ2H=WN*(DEP2-SDP)
        WZ2Z=WN*(DEP2-ODP)
        EXR2=DEXP(-WZ2H-WZ2Z)
        CPR2=EXR2*(CUP1X-WZ2H*CUP2X)
        CSR2=EXR2*(CUS1X-WZ2H*CUS2X-WZ2Z*CUP1X+WZ2Z*WZ2H*CUP2X)
      END IF
C
C   * DOWN-GOING WAVE DUE TO REFLECT WAVE FROM UPPER BOUNDARY *
C
      WHZ1=WN*(SDP-DEP1)
      WZZ1=WN*(ODP-DEP1)
      EXR1=DEXP(-WHZ1-WZZ1)
C
      CPR1=EXR1*(CDP1X-WHZ1*CDP2X)
      CSR1=EXR1*(CDS1X-WHZ1*CDS2X-WZZ1*CDP1X+WZZ1*WHZ1*CDP2X)
C
C   * DISP-STRESS VECTORS *
C
      H1X=H1D*EXD+H1R1*EXR1+H1R2*EXR2
      H2X=WN*(-H2D*EXD+H2R1*EXR1+H2R2*EXR2)
C
      A=X2-1.0
      B=X2-3.0
      V1X=CPD+CSD+CPR1+CSR1+CPR2+CSR2
      V2X=-A*(CPD+CPR1)+CSD+CSR1+A*CPR2-CSR2
      V3X=WN*SM2*(B*(CPD+CPR1)-2.0*(CSD+CSR1)-B*CPR2+2.0*CSR2)
      V4X=WN*SM2*(A*(CPD+CPR1)-2.0*(CSD+CSR1)+A*CPR2-2.0*CSR2)
C
C   * DISP-STRESS VECTOR DUE TO Qz *
C
      CPD=EXD*CDPZ
      CSD=EXD*(CDSZ-WZH*CDPZ)
C
      IF(LS.EQ.NL)THEN
        CPR2=(0.D0,0.D0)
        CSR2=(0.D0,0.D0)
      ELSE
        CPD=EXD*CDPZ
        CSD=EXD*(CDSZ-WZH*CDPZ)
        CPR2=EXR2*(CUP1Z-WZ2H*CUP2Z)
        CSR2=EXR2*(CUS1Z-WZ2H*CUS2Z-WZ2Z*CUP1Z+WZ2Z*WZ2H*CUP2Z)
      END IF
C
      CPR1=EXR1*(CDP1Z-WHZ1*CDP2Z)
      CSR1=EXR1*(CDS1Z-WHZ1*CDS2Z-WZZ1*CDP1Z+WZZ1*WHZ1*CDP2Z)
C
      V1Z=CPD+CSD+CPR1+CSR1+CPR2+CSR2
      V2Z=-A*(CPD+CPR1)+CSD+CSR1+A*CPR2-CSR2
      V3Z=WN*SM2*(B*(CPD+CPR1)-2.0*(CSD+CSR1)-B*CPR2+2.0*CSR2)
      V4Z=WN*SM2*(A*(CPD+CPR1)-2.0*(CSD+CSR1)+A*CPR2-2.0*CSR2)
C
      GOTO 50
C
C *** DISP-STRESS VECTOR DUE TO Qx FOR DEP2 < (OR =) ODP ***
C
   30 WZH=WN*(ODP-SDP)
      WZ2H=WN*(DEP2-SDP)
      WZZ2=WN*(ODP-DEP2)
      EX=DEXP(-WZH)
C
      CPT=EX*(CDP1X-WZ2H*CDP2X)
      CST=EX*(CDS1X-WZZ2*CDP1X-WZ2H*CDS2X+WZZ2*WZ2H*CDP2X)
C
      H1X=EX*H1T
      H2X=WN*EX*H2T
C
      A=X3-1.0
      B=X3-3.0
      V1X=CPT+CST
      V2X=-A*CPT+CST
      V3X=WN*SM3*(B*CPT-2.0*CST)
      V4X=WN*SM3*(A*CPT-2.0*CST)
C
C   * DISP-STRESS VECTOR DUE TO Qz *
C
      CPT=EX*(CDP1Z-WZ2H*CDP2Z)
      CST=EX*(CDS1Z-WZZ2*CDP1Z-WZ2H*CDS2Z+WZZ2*WZ2H*CDP2Z)
C
      V1Z=CPT+CST
      V2Z=-A*CPT+CST
      V3Z=WN*SM3*(B*CPT-2.0*CST)
      V4Z=WN*SM3*(A*CPT-2.0*CST)
C
   50 RETURN
      END
C
C
C *** ANALITICAL INTEGRATIONS OF ASYMPTOTIC SOLUTIONS ***
C
      SUBROUTINE AINTRZ(LS,NL,QXS,QZS,VXS,WXS,WZS)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 QXS(0:2),QZS(0:2),VXS(0:2),WXS(0:2),WZS(0:2)
     *   ,SM1,X1,SM2,X2,SM3,X3,RM1,RM2,RM3,SM,RM,DC1,DC2,DC3
     *   ,H1T,H2T,H1D,H2D,H1R1,H2R1,H1R2,H2R2       
     *   ,CUPX,CUSX,CUP1X,CUS1X,CUP2X,CUS2X
     *   ,CDPX,CDSX,CDP1X,CDS1X,CDP2X,CDS2X
     *   ,CUPZ,CUSZ,CUP1Z,CUS1Z,CUP2Z,CUS2Z
     *   ,CDPZ,CDSZ,CDP1Z,CDS1Z,CDP2Z,CDS2Z
     *   ,DCUPT,DCUST,DCUPD,DCUSD,DCUPR,DCUSR
     *   ,DCDPT,DCDST,DCDPD,DCDSD,DCDPR,DCDSR
     *   ,D1MH1,DD0H1,D11H1,D1MH2,DD0H2,D1MH1R,DD0H1R,D1MH1Z,DD0H1Z
     *   ,D1MV1,DD0V1,D11V1,D01V2,D10V2,D1MV3,DD0V3,D10V4
     *   ,D1MV1R,DD0V1R,D10V2R,D1MV1Z,DD0V1Z,D10V2Z
     *   ,D01V1,D10V1,D00V2,D11V2,D10V3,D00V4
     *   ,D10V1R,D00V2R,D10V1Z,D00V2Z
      COMMON/ASYMP/SM1,X1,SM2,X2,SM3,X3 /RAM/RM1,RM2,RM3
     *      /DPSO /DEP1,DEP2,SDP,ODP /FCQVW/SM,RM,R
     *      /HX  /H1T,H2T,H1D,H2D,H1R1,H2R1,H1R2,H2R2       
     *      /VUX /CUPX,CUSX,CUP1X,CUS1X,CUP2X,CUS2X
     *      /VDX /CDPX,CDSX,CDP1X,CDS1X,CDP2X,CDS2X
     *      /VUZ /CUPZ,CUSZ,CUP1Z,CUS1Z,CUP2Z,CUS2Z
     *      /VDZ /CDPZ,CDSZ,CDP1Z,CDS1Z,CDP2Z,CDS2Z
C
C *** ANALITICAL INTEGRATIONS ***
C
C   * INTEGRATIONS FOR DIRECT WAVE *
C
      ZH=DABS(SDP-ODP)
      ZH2=ZH*ZH
      RR=DSQRT(R*R+ZH2)
      R2=RR*RR
      R3=R2*RR
      R5=R3*R2
C
      S00=1.D0/RR
      S01=ZH/R3
      S02=(3.D0*ZH2/R2-1.D0)/R3
      S03=3.D0*ZH*(5.D0*ZH2/R2-3.D0)/R5
C
      S1M=(RR-ZH)/R
      S10=S1M/RR
      S11=R/R3
      S12=3.D0*R*ZH/R5
      S13=3.D0*R*(5.D0*ZH2/R2-1.D0)/R5
C
C   * INTEGRATIONS FOR REFLECTED WAVE FROM UPPER BOUNDARY *
C
      ZH=SDP+ODP-2.D0*DEP1
      ZH2=ZH*ZH
      RR=DSQRT(R*R+ZH2)
      R2=RR*RR
      R3=R2*RR
      R5=R3*R2
C
      S001=1.D0/RR
      S011=ZH/R3
      S021=(3.D0*ZH2/R2-1.D0)/R3
      S031=3.D0*ZH*(5.D0*ZH2/R2-3.D0)/R5
C
      S1M1=(RR-ZH)/R
      S101=S1M1/RR
      S111=R/R3
      S121=3.D0*R*ZH/R5
      S131=3.D0*R*(5.D0*ZH2/R2-1.D0)/R5
C
C   * INTEGRATIONS FOR REFLECTED WAVE FROM LOWER BOUNDARY *
C
      ZH=2.D0*DEP2-SDP-ODP
      ZH2=ZH*ZH
      RR=DSQRT(R*R+ZH2)
      R2=RR*RR
      R3=R2*RR
      R5=R3*R2
C
      S002=1.D0/RR
      S012=ZH/R3
      S022=(3.D0*ZH2/R2-1.D0)/R3
      S032=3.D0*ZH*(5.D0*ZH2/R2-3.D0)/R5
C
      S1M2=(RR-ZH)/R
      S102=S1M2/RR
      S112=R/R3
      S122=3.D0*R*ZH/R5
      S132=3.D0*R*(5.D0*ZH2/R2-1.D0)/R5
C
C *** ODP < DEP1 ***
C
      IF(ODP.GE.DEP1) GOTO 10
      IF(LS.EQ.1) GOTO 10
C
C  ** QX **
C
C   * SH (TRANSMITTED) *
C
      D1MH1=H1T*S1M
      DD0H1=H1T*(S00-S1M/R)
      QXS(0)=D1MH1/R
      VXS(0)=DD0H1
C
      D11H1=H1T*S11
      D1MH1R=DD0H1
      DD0H1R=(-DD0H1+D1MH1/R)/R-D11H1
      QXS(1)=(-D1MH1/R+D1MH1R)/R
      VXS(1)=DD0H1R
C
      D1MH2=H2T*S10
      DD0H2=H2T*(S01-S10/R)
      D1MH1Z=D1MH2/SM1
      DD0H1Z=DD0H2/SM1
      QXS(2)=D1MH1Z/R
      VXS(2)=DD0H1Z
C
C   * P-SV (TRANSMITTED) *
C
      HD1=SDP-DEP1
      D1Z=DEP1-ODP
      DC1=HD1*CUP2X
      DC2=HD1*CUS2X+D1Z*CUP1X
      DC3=D1Z*DC1
      DCUPT=CUP1X*S1M-DC1*S10
      DCUST=CUS1X*S1M-DC2*S10+DC3*S11
      D1MV1=DCUPT+DCUST
C
      DCUPT=CUP1X*S11-DC1*S12
      DCUST=CUS1X*S11-DC2*S12+DC3*S13
      D11V1=DCUPT+DCUST
      D10V4=SM1*((X1-1.0)*DCUPT-2.0*DCUST)
C
      S1=S00-S1M/R
      S2=S01-S10/R
      S3=S02-S11/R
      DCUPT=CUP1X*S1-DC1*S2
      DCUST=CUS1X*S1-DC2*S2+DC3*S3
      DD0V1=DCUPT+DCUST
C
      DCUPT=CUP1X*S01-DC1*S02
      DCUST=CUS1X*S01-DC2*S02+DC3*S03
      D01V2=(X1-1.0)*DCUPT-DCUST
C
      DCUPT=CUP1X*S10-DC1*S11
      DCUST=CUS1X*S10-DC2*S11+DC3*S12
      D10V2=(X1-1.0)*DCUPT-DCUST
      D1MV3=SM1*(-(X1-3.0)*DCUPT+2.0*DCUST)
C
      S1=S01-S10/R
      S2=S02-S11/R
      S3=S03-S12/R
      DCUPT=CUP1X*S1-DC1*S2
      DCUST=CUS1X*S1-DC2*S2+DC3*S3
      DD0V3=SM1*(-(X1-3.0)*DCUPT+2.0*DCUST)
C
      QXS(0)=QXS(0)+DD0V1
      VXS(0)=VXS(0)+D1MV1/R
      WXS(0)=D10V2
C
      D1MV1R=DD0V1
      DD0V1R=(-DD0V1+D1MV1/R)/R-D11V1
      D10V2R=D01V2-D10V2/R
C
      QXS(1)=QXS(1)+DD0V1R
      VXS(1)=VXS(1)+(-D1MV1/R+D1MV1R)/R
      WXS(1)=D10V2R
C
      D1MV1Z=D1MV3/SM1+D10V2
      DD0V1Z=DD0V3/SM1+D01V2-D10V2/R
      D10V2Z=(D10V4-RM1*D11V1)/(RM1+2.D0*SM1) 
C
      QXS(2)=QXS(2)+DD0V1Z
      VXS(2)=VXS(2)+D1MV1Z/R
      WXS(2)=D10V2Z
C
C  ** QZ **
C
C   * P-SV (TRANSMITTED) *
C
      DC1=HD1*CUP2Z
      DC2=HD1*CUS2Z+D1Z*CUP1Z
      DC3=D1Z*DC1
      DCUPT=CUP1Z*S01-DC1*S02
      DCUST=CUS1Z*S01-DC2*S02+DC3*S03
      D01V1=DCUPT+DCUST
      D00V4=SM1*((X1-1.0)*DCUPT-2.0*DCUST)
C
      DCUPT=CUP1Z*S10-DC1*S11
      DCUST=CUS1Z*S10-DC2*S11+DC3*S12
      D10V1=DCUPT+DCUST
C
      DCUPT=CUP1Z*S00-DC1*S01
      DCUST=CUS1Z*S00-DC2*S01+DC3*S02
      D00V2=(X1-1.0)*DCUPT-DCUST
C
      DCUPT=CUP1Z*S11-DC1*S12
      DCUST=CUS1Z*S11-DC2*S12+DC3*S13
      D11V2=(X1-1.0)*DCUPT-DCUST
      D10V3=SM1*(-(X1-3.0)*DCUPT+2.0*DCUST)
C
      QZS(0)=D10V1
      WZS(0)=D00V2
C
      D10V1R=D01V1-D10V1/R
      D00V2R=-D11V2
C
      QZS(1)=D10V1R
      WZS(1)=D00V2R
C
      D10V1Z=D10V3/SM1+D11V2
      D00V2Z=(D00V4-RM1*D01V1)/(RM1+2.D0*SM1)
C
      QZS(2)=D10V1Z
      WZS(2)=D00V2Z
C
      GOTO 50
C
C *** DEP1 < (OR =) ODP < SDP ***
C
   10 IF(ODP.GE.SDP) GOTO 20
C
C  ** QX **
C
C   * SH *
C
      D1MH1=H1D*S1M+H1R1*S1M1+H1R2*S1M2
      DD0H1=H1D*(S00-S1M/R)+H1R1*(S001-S1M1/R)+H1R2*(S002-S1M2/R)
      QXS(0)=D1MH1/R
      VXS(0)=DD0H1
C
      D11H1=H1D*S11+H1R1*S111+H1R2*S112
      D1MH1R=DD0H1
      DD0H1R=(-DD0H1+D1MH1/R)/R-D11H1
      QXS(1)=(-D1MH1/R+D1MH1R)/R
      VXS(1)=DD0H1R
C
      D1MH2=H2D*S10+H2R1*S101+H2R2*S102
      DD0H2=H2D*(S01-S10/R)+H2R1*(S011-S101/R)+H2R2*(S012-S102/R)
      D1MH1Z=D1MH2/SM2
      DD0H1Z=DD0H2/SM2
      QXS(2)=D1MH1Z/R
      VXS(2)=DD0H1Z
C
C   * P-SV (UP-GOING DIRECT) *
C
      DC1=(SDP-ODP)*CUPX
      DCUPD=CUPX*S1M
      DCUSD=CUSX*S1M-DC1*S10
      D1MV1=DCUPD+DCUSD
C
      DCUPD=CUPX*S11
      DCUSD=CUSX*S11-DC1*S12
      D11V1=DCUPD+DCUSD
      D10V4=SM2*((X2-1.0)*DCUPD-2.0*DCUSD)
C
      S1=S00-S1M/R
      S2=S01-S10/R
      DCUPD=CUPX*S1
      DCUSD=CUSX*S1-DC1*S2
      DD0V1=DCUPD+DCUSD
C
      DCUPD=CUPX*S01
      DCUSD=CUSX*S01-DC1*S02
      D01V2=(X2-1.0)*DCUPD-DCUSD
C
      DCUPD=CUPX*S10
      DCUSD=CUSX*S10-DC1*S11
      D10V2=(X2-1.0)*DCUPD-DCUSD
      D1MV3=SM2*(-(X2-3.0)*DCUPD+2.0*DCUSD)
C
      S1=S01-S10/R
      S2=S02-S11/R
      DCUPD=CUPX*S1
      DCUSD=CUSX*S1-DC1*S2
      DD0V3=SM2*(-(X2-3.0)*DCUPD+2.0*DCUSD)
C
      QXS(0)=QXS(0)+DD0V1
      VXS(0)=VXS(0)+D1MV1/R
      WXS(0)=D10V2
C
      D1MV1R=DD0V1
      DD0V1R=(-DD0V1+D1MV1/R)/R-D11V1
      D10V2R=D01V2-D10V2/R
C
      QXS(1)=QXS(1)+DD0V1R
      VXS(1)=VXS(1)+(-D1MV1/R+D1MV1R)/R
      WXS(1)=D10V2R
C
      D1MV1Z=D1MV3/SM2+D10V2
      DD0V1Z=DD0V3/SM2+D01V2-D10V2/R
      D10V2Z=(D10V4-RM2*D11V1)/(RM2+2.D0*SM2) 
C
      QXS(2)=QXS(2)+DD0V1Z
      VXS(2)=VXS(2)+D1MV1Z/R
      WXS(2)=D10V2Z
C
C   * P-SV (DOWN-GOING REFLECTED FROM UPPER BOUNDARY) *
C
      HD1=SDP-DEP1
      ZD1=ODP-DEP1
      DC1=HD1*CDP2X
      DC2=HD1*CDS2X+ZD1*CDP1X
      DC3=ZD1*DC1
      DCDPR=CDP1X*S1M1-DC1*S101
      DCDSR=CDS1X*S1M1-DC2*S101+DC3*S111
      D1MV1=DCDPR+DCDSR
C
      DCDPR=CDP1X*S111-DC1*S121
      DCDSR=CDS1X*S111-DC2*S121+DC3*S131
      D11V1=DCDPR+DCDSR
      D10V4=SM2*((X2-1.0)*DCDPR-2.0*DCDSR)
C
      S1=S001-S1M1/R
      S2=S011-S101/R
      S3=S021-S111/R
      DCDPR=CDP1X*S1-DC1*S2
      DCDSR=CDS1X*S1-DC2*S2+DC3*S3
      DD0V1=DCDPR+DCDSR
C
      DCDPR=CDP1X*S011-DC1*S021
      DCDSR=CDS1X*S011-DC2*S021+DC3*S031
      D01V2=-(X2-1.0)*DCDPR+DCDSR
C
      DCDPR=CDP1X*S101-DC1*S111
      DCDSR=CDS1X*S101-DC2*S111+DC3*S121
      D10V2=-(X2-1.0)*DCDPR+DCDSR
      D1MV3=SM2*((X2-3.0)*DCDPR-2.0*DCDSR)
C
      S1=S011-S101/R
      S2=S021-S111/R
      S3=S031-S121/R
      DCDPR=CDP1X*S1-DC1*S2
      DCDSR=CDS1X*S1-DC2*S2+DC3*S3
      DD0V3=SM2*((X2-3.0)*DCDPR-2.0*DCDSR)
C
      QXS(0)=QXS(0)+DD0V1
      VXS(0)=VXS(0)+D1MV1/R
      WXS(0)=WXS(0)+D10V2
C
      D1MV1R=DD0V1
      DD0V1R=(-DD0V1+D1MV1/R)/R-D11V1
      D10V2R=D01V2-D10V2/R
C
      QXS(1)=QXS(1)+DD0V1R
      VXS(1)=VXS(1)+(-D1MV1/R+D1MV1R)/R
      WXS(1)=WXS(1)+D10V2R
C
      D1MV1Z=D1MV3/SM2+D10V2
      DD0V1Z=DD0V3/SM2+D01V2-D10V2/R
      D10V2Z=(D10V4-RM2*D11V1)/(RM2+2.D0*SM2) 
C
      QXS(2)=QXS(2)+DD0V1Z
      VXS(2)=VXS(2)+D1MV1Z/R
      WXS(2)=WXS(2)+D10V2Z
C
C   * P-SV (UP-GOING REFLECTED FROM LOWER BOUNDARY) *
C
      IF(LS.EQ.NL)GOTO 15
C
      D2H=DEP2-SDP
      D2Z=DEP2-ODP
      DC1=D2H*CUP2X
      DC2=D2H*CUS2X+D2Z*CUP1X
      DC3=D2Z*DC1
      DCUPR=CUP1X*S1M2-DC1*S102
      DCUSR=CUS1X*S1M2-DC2*S102+DC3*S112
      D1MV1=DCUPR+DCUSR
C
      DCUPR=CUP1X*S112-DC1*S122
      DCUSR=CUS1X*S112-DC2*S122+DC3*S132
      D11V1=DCUPR+DCUSR
      D10V4=SM2*((X2-1.0)*DCUPR-2.0*DCUSR)
C
      S1=S002-S1M2/R
      S2=S012-S102/R
      S3=S022-S112/R
      DCUPR=CUP1X*S1-DC1*S2
      DCUSR=CUS1X*S1-DC2*S2+DC3*S3
      DD0V1=DCUPR+DCUSR
C
      DCUPR=CUP1X*S012-DC1*S022
      DCUSR=CUS1X*S012-DC2*S022+DC3*S032
      D01V2=(X2-1.0)*DCUPR-DCUSR
C
      DCUPR=CUP1X*S102-DC1*S112
      DCUSR=CUS1X*S102-DC2*S112+DC3*S122
      D10V2=(X2-1.0)*DCUPR-DCUSR
      D1MV3=SM2*(-(X2-3.0)*DCUPR+2.0*DCUSR)
C
      S1=S012-S102/R
      S2=S022-S112/R
      S3=S032-S122/R
      DCUPR=CUP1X*S1-DC1*S2
      DCUSR=CUS1X*S1-DC2*S2+DC3*S3
      DD0V3=SM2*(-(X2-3.0)*DCUPR+2.0*DCUSR)
C
      QXS(0)=QXS(0)+DD0V1
      VXS(0)=VXS(0)+D1MV1/R
      WXS(0)=WXS(0)+D10V2
C
      D1MV1R=DD0V1
      DD0V1R=(-DD0V1+D1MV1/R)/R-D11V1
      D10V2R=D01V2-D10V2/R
C
      QXS(1)=QXS(1)+DD0V1R
      VXS(1)=VXS(1)+(-D1MV1/R+D1MV1R)/R
      WXS(1)=WXS(1)+D10V2R
C
      D1MV1Z=D1MV3/SM2+D10V2
      DD0V1Z=DD0V3/SM2+D01V2-D10V2/R
      D10V2Z=(D10V4-RM2*D11V1)/(RM2+2.D0*SM2) 
C
      QXS(2)=QXS(2)+DD0V1Z
      VXS(2)=VXS(2)+D1MV1Z/R
      WXS(2)=WXS(2)+D10V2Z
C
C  ** QZ **
C
   15 DC1=(SDP-ODP)*CUPZ
      DCUPD=CUPZ*S01
      DCUSD=CUSZ*S01-DC1*S02
      D01V1=DCUPD+DCUSD
      D00V4=SM2*((X2-1.0)*DCUPD-2.0*DCUSD)
C
      DCUPD=CUPZ*S10
      DCUSD=CUSZ*S10-DC1*S11
      D10V1=DCUPD+DCUSD
C
      DCUPD=CUPZ*S00
      DCUSD=CUSZ*S00-DC1*S01
      D00V2=(X2-1.0)*DCUPD-DCUSD
C
      DCUPD=CUPZ*S11
      DCUSD=CUSZ*S11-DC1*S12
      D11V2=(X2-1.0)*DCUPD-DCUSD
      D10V3=SM2*(-(X2-3.0)*DCUPD+2.0*DCUSD)
C
      QZS(0)=D10V1
      WZS(0)=D00V2
C
      D10V1R=D01V1-D10V1/R
      D00V2R=-D11V2
C
      QZS(1)=D10V1R
      WZS(1)=D00V2R
C
      D10V1Z=D10V3/SM2+D11V2
      D00V2Z=(D00V4-RM2*D01V1)/(RM2+2.D0*SM2)
C
      QZS(2)=D10V1Z
      WZS(2)=D00V2Z
C
      DC1=HD1*CDP2Z
      DC2=HD1*CDS2Z+ZD1*CDP1Z
      DC3=ZD1*DC1
      DCDPR=CDP1Z*S011-DC1*S021
      DCDSR=CDS1Z*S011-DC2*S021+DC3*S031
      D01V1=DCDPR+DCDSR
      D00V4=SM2*((X2-1.0)*DCDPR-2.0*DCDSR)
C
      DCDPR=CDP1Z*S101-DC1*S111
      DCDSR=CDS1Z*S101-DC2*S111+DC3*S121
      D10V1=DCDPR+DCDSR
C
      DCDPR=CDP1Z*S001-DC1*S011
      DCDSR=CDS1Z*S001-DC2*S011+DC3*S021
      D00V2=-(X2-1.0)*DCDPR+DCDSR
C
      DCDPR=CDP1Z*S111-DC1*S121
      DCDSR=CDS1Z*S111-DC2*S121+DC3*S131
      D11V2=-(X2-1.0)*DCDPR+DCDSR
      D10V3=SM2*((X2-3.0)*DCDPR-2.0*DCDSR)
C
      QZS(0)=QZS(0)+D10V1
      WZS(0)=WZS(0)+D00V2
C
      D10V1R=D01V1-D10V1/R
      D00V2R=-D11V2
C
      QZS(1)=QZS(1)+D10V1R
      WZS(1)=WZS(1)+D00V2R
C
      D10V1Z=D10V3/SM2+D11V2
      D00V2Z=(D00V4-RM2*D01V1)/(RM2+2.D0*SM2)
C
      QZS(2)=QZS(2)+D10V1Z
      WZS(2)=WZS(2)+D00V2Z
C
      IF(LS.EQ.NL)GOTO 50
C
      DC1=D2H*CUP2Z
      DC2=D2H*CUS2Z+D2Z*CUP1Z
      DC3=D2Z*DC1
      DCUPR=CUP1Z*S012-DC1*S022
      DCUSR=CUS1Z*S012-DC2*S022+DC3*S032
      D01V1=DCUPR+DCUSR
      D00V4=SM2*((X2-1.0)*DCUPR-2.0*DCUSR)
C
      DCUPR=CUP1Z*S102-DC1*S112
      DCUSR=CUS1Z*S102-DC2*S112+DC3*S122
      D10V1=DCUPR+DCUSR
C
      DCUPR=CUP1Z*S002-DC1*S012
      DCUSR=CUS1Z*S002-DC2*S012+DC3*S022
      D00V2=(X2-1.0)*DCUPR-DCUSR
C
      DCUPR=CUP1Z*S112-DC1*S122
      DCUSR=CUS1Z*S112-DC2*S122+DC3*S132
      D11V2=(X2-1.0)*DCUPR-DCUSR
      D10V3=SM2*(-(X2-3.0)*DCUPR+2.0*DCUSR)
C
      QZS(0)=QZS(0)+D10V1
      WZS(0)=WZS(0)+D00V2
C
      D10V1R=D01V1-D10V1/R
      D00V2R=-D11V2
C
      QZS(1)=QZS(1)+D10V1R
      WZS(1)=WZS(1)+D00V2R
C
      D10V1Z=D10V3/SM2+D11V2
      D00V2Z=(D00V4-RM2*D01V1)/(RM2+2.D0*SM2)
C
      QZS(2)=QZS(2)+D10V1Z
      WZS(2)=WZS(2)+D00V2Z
C
      GOTO 50
C
C *** SDP < (OR =) ODP < DEP2 ***
C
   20 IF((LS.NE.NL).AND.(ODP.GE.DEP2)) GOTO 30
C
C  ** QX **
C
C   * SH *
C
      D1MH1=H1D*S1M+H1R1*S1M1+H1R2*S1M2
      DD0H1=H1D*(S00-S1M/R)+H1R1*(S001-S1M1/R)+H1R2*(S002-S1M2/R)
      QXS(0)=D1MH1/R
      VXS(0)=DD0H1
C
      D11H1=H1D*S11+H1R1*S111+H1R2*S112
      D1MH1R=DD0H1
      DD0H1R=(-DD0H1+D1MH1/R)/R-D11H1
      QXS(1)=(-D1MH1/R+D1MH1R)/R
      VXS(1)=DD0H1R
C
      D1MH2=-H2D*S10+H2R1*S101+H2R2*S102
      DD0H2=-H2D*(S01-S10/R)+H2R1*(S011-S101/R)+H2R2*(S012-S102/R)
      D1MH1Z=D1MH2/SM2
      DD0H1Z=DD0H2/SM2
      QXS(2)=D1MH1Z/R
      VXS(2)=DD0H1Z
C
C   * P-SV (DOWN-GOING DIRECT) *
C
      DC1=(ODP-SDP)*CDPX
      DCDPD=CDPX*S1M
      DCDSD=CDSX*S1M-DC1*S10
      D1MV1=DCDPD+DCDSD
C
      DCDPD=CDPX*S11
      DCDSD=CDSX*S11-DC1*S12
      D11V1=DCDPD+DCDSD
      D10V4=SM2*((X2-1.0)*DCDPD-2.0*DCDSD)
C
      S1=S00-S1M/R
      S2=S01-S10/R
      DCDPD=CDPX*S1
      DCDSD=CDSX*S1-DC1*S2
      DD0V1=DCDPD+DCDSD
C
      DCDPD=CDPX*S01
      DCDSD=CDSX*S01-DC1*S02
      D01V2=-(X2-1.0)*DCDPD+DCDSD
C
      DCDPD=CDPX*S10
      DCDSD=CDSX*S10-DC1*S11
      D10V2=-(X2-1.0)*DCDPD+DCDSD
      D1MV3=SM2*((X2-3.0)*DCDPD-2.0*DCDSD)
C
      S1=S01-S10/R
      S2=S02-S11/R
      DCDPD=CDPX*S1
      DCDSD=CDSX*S1-DC1*S2
      DD0V3=SM2*((X2-3.0)*DCDPD-2.0*DCDSD)
C
      QXS(0)=QXS(0)+DD0V1
      VXS(0)=VXS(0)+D1MV1/R
      WXS(0)=D10V2
C
      D1MV1R=DD0V1
      DD0V1R=(-DD0V1+D1MV1/R)/R-D11V1
      D10V2R=D01V2-D10V2/R
C
      QXS(1)=QXS(1)+DD0V1R
      VXS(1)=VXS(1)+(-D1MV1/R+D1MV1R)/R
      WXS(1)=D10V2R
C
      D1MV1Z=D1MV3/SM2+D10V2
      DD0V1Z=DD0V3/SM2+D01V2-D10V2/R
      D10V2Z=(D10V4-RM2*D11V1)/(RM2+2.D0*SM2) 
C
      QXS(2)=QXS(2)+DD0V1Z
      VXS(2)=VXS(2)+D1MV1Z/R
      WXS(2)=D10V2Z
C
C   * P-SV (UP-GOING REFLECTED FROM LOWER BOUNDARY) *
C
      IF(LS.EQ.NL)GOTO 25
C
      D2H=DEP2-SDP
      D2Z=DEP2-ODP
      DC1=D2H*CUP2X
      DC2=D2H*CUS2X+D2Z*CUP1X
      DC3=D2Z*DC1
      DCUPR=CUP1X*S1M2-DC1*S102
      DCUSR=CUS1X*S1M2-DC2*S102+DC3*S112
      D1MV1=DCUPR+DCUSR
C
      DCUPR=CUP1X*S112-DC1*S122
      DCUSR=CUS1X*S112-DC2*S122+DC3*S132
      D11V1=DCUPR+DCUSR
      D10V4=SM2*((X2-1.0)*DCUPR-2.0*DCUSR)
C
      S1=S002-S1M2/R
      S2=S012-S102/R
      S3=S022-S112/R
      DCUPR=CUP1X*S1-DC1*S2
      DCUSR=CUS1X*S1-DC2*S2+DC3*S3
      DD0V1=DCUPR+DCUSR
C
      DCUPR=CUP1X*S012-DC1*S022
      DCUSR=CUS1X*S012-DC2*S022+DC3*S032
      D01V2=(X2-1.0)*DCUPR-DCUSR
C
      DCUPR=CUP1X*S102-DC1*S112
      DCUSR=CUS1X*S102-DC2*S112+DC3*S122
      D10V2=(X2-1.0)*DCUPR-DCUSR
      D1MV3=SM2*(-(X2-3.0)*DCUPR+2.0*DCUSR)
C
      S1=S012-S102/R
      S2=S022-S112/R
      S3=S032-S122/R
      DCUPR=CUP1X*S1-DC1*S2
      DCUSR=CUS1X*S1-DC2*S2+DC3*S3
      DD0V3=SM2*(-(X2-3.0)*DCUPR+2.0*DCUSR)
C
      QXS(0)=QXS(0)+DD0V1
      VXS(0)=VXS(0)+D1MV1/R
      WXS(0)=WXS(0)+D10V2
C
      D1MV1R=DD0V1
      DD0V1R=(-DD0V1+D1MV1/R)/R-D11V1
      D10V2R=D01V2-D10V2/R
C
      QXS(1)=QXS(1)+DD0V1R
      VXS(1)=VXS(1)+(-D1MV1/R+D1MV1R)/R
      WXS(1)=WXS(1)+D10V2R
C
      D1MV1Z=D1MV3/SM2+D10V2
      DD0V1Z=DD0V3/SM2+D01V2-D10V2/R
      D10V2Z=(D10V4-RM2*D11V1)/(RM2+2.D0*SM2) 
C
      QXS(2)=QXS(2)+DD0V1Z
      VXS(2)=VXS(2)+D1MV1Z/R
      WXS(2)=WXS(2)+D10V2Z
C
C   * P-SV (DOWN-GOING REFLECTED FROM UPPER BOUNDARY) *
C
   25 HD1=SDP-DEP1
      ZD1=ODP-DEP1
      DC1=HD1*CDP2X
      DC2=HD1*CDS2X+ZD1*CDP1X
      DC3=ZD1*DC1
      DCDPR=CDP1X*S1M1-DC1*S101
      DCDSR=CDS1X*S1M1-DC2*S101+DC3*S111
      D1MV1=DCDPR+DCDSR
C
      DCDPR=CDP1X*S111-DC1*S121
      DCDSR=CDS1X*S111-DC2*S121+DC3*S131
      D11V1=DCDPR+DCDSR
      D10V4=SM2*((X2-1.0)*DCDPR-2.0*DCDSR)
C
      S1=S001-S1M1/R
      S2=S011-S101/R
      S3=S021-S111/R
      DCDPR=CDP1X*S1-DC1*S2
      DCDSR=CDS1X*S1-DC2*S2+DC3*S3
      DD0V1=DCDPR+DCDSR
C
      DCDPR=CDP1X*S011-DC1*S021
      DCDSR=CDS1X*S011-DC2*S021+DC3*S031
      D01V2=-(X2-1.0)*DCDPR+DCDSR
C
      DCDPR=CDP1X*S101-DC1*S111
      DCDSR=CDS1X*S101-DC2*S111+DC3*S121
      D10V2=-(X2-1.0)*DCDPR+DCDSR
      D1MV3=SM2*((X2-3.0)*DCDPR-2.0*DCDSR)
C
      S1=S011-S101/R
      S2=S021-S111/R
      S3=S031-S121/R
      DCDPR=CDP1X*S1-DC1*S2
      DCDSR=CDS1X*S1-DC2*S2+DC3*S3
      DD0V3=SM2*((X2-3.0)*DCDPR-2.0*DCDSR)
C
      QXS(0)=QXS(0)+DD0V1
      VXS(0)=VXS(0)+D1MV1/R
      WXS(0)=WXS(0)+D10V2
C
      D1MV1R=DD0V1
      DD0V1R=(-DD0V1+D1MV1/R)/R-D11V1
      D10V2R=D01V2-D10V2/R
C
      QXS(1)=QXS(1)+DD0V1R
      VXS(1)=VXS(1)+(-D1MV1/R+D1MV1R)/R
      WXS(1)=WXS(1)+D10V2R
C
      D1MV1Z=D1MV3/SM2+D10V2
      DD0V1Z=DD0V3/SM2+D01V2-D10V2/R
      D10V2Z=(D10V4-RM2*D11V1)/(RM2+2.D0*SM2) 
C
      QXS(2)=QXS(2)+DD0V1Z
      VXS(2)=VXS(2)+D1MV1Z/R
      WXS(2)=WXS(2)+D10V2Z
C
C  ** QZ **
C
      DC1=(ODP-SDP)*CDPZ
      DCDPD=CDPZ*S01
      DCDSD=CDSZ*S01-DC1*S02
      D01V1=DCDPD+DCDSD
      D00V4=SM2*((X2-1.0)*DCDPD-2.0*DCDSD)
C
      DCDPD=CDPZ*S10
      DCDSD=CDSZ*S10-DC1*S11
      D10V1=DCDPD+DCDSD
C
      DCDPD=CDPZ*S00
      DCDSD=CDSZ*S00-DC1*S01
      D00V2=-(X2-1.0)*DCDPD+DCDSD
C
      DCDPD=CDPZ*S11
      DCDSD=CDSZ*S11-DC1*S12
      D11V2=-(X2-1.0)*DCDPD+DCDSD
      D10V3=SM2*((X2-3.0)*DCDPD-2.0*DCDSD)
C
      QZS(0)=D10V1
      WZS(0)=D00V2
C
      D10V1R=D01V1-D10V1/R
      D00V2R=-D11V2
C
      QZS(1)=D10V1R
      WZS(1)=D00V2R
C
      D10V1Z=D10V3/SM2+D11V2
      D00V2Z=(D00V4-RM2*D01V1)/(RM2+2.D0*SM2)
C
      QZS(2)=D10V1Z
      WZS(2)=D00V2Z
C
      IF(LS.EQ.NL)GOTO 27
C
      DC1=D2H*CUP2Z
      DC2=D2H*CUS2Z+D2Z*CUP1Z
      DC3=D2Z*DC1
      DCUPR=CUP1Z*S012-DC1*S022
      DCUSR=CUS1Z*S012-DC2*S022+DC3*S032
      D01V1=DCUPR+DCUSR
      D00V4=SM2*((X2-1.0)*DCUPR-2.0*DCUSR)
C
      DCUPR=CUP1Z*S102-DC1*S112
      DCUSR=CUS1Z*S102-DC2*S112+DC3*S122
      D10V1=DCUPR+DCUSR
C
      DCUPR=CUP1Z*S002-DC1*S012
      DCUSR=CUS1Z*S002-DC2*S012+DC3*S022
      D00V2=(X2-1.0)*DCUPR-DCUSR
C
      DCUPR=CUP1Z*S112-DC1*S122
      DCUSR=CUS1Z*S112-DC2*S122+DC3*S132
      D11V2=(X2-1.0)*DCUPR-DCUSR
      D10V3=SM2*(-(X2-3.0)*DCUPR+2.0*DCUSR)
C
      QZS(0)=QZS(0)+D10V1
      WZS(0)=WZS(0)+D00V2
C
      D10V1R=D01V1-D10V1/R
      D00V2R=-D11V2
C
      QZS(1)=QZS(1)+D10V1R
      WZS(1)=WZS(1)+D00V2R
C
      D10V1Z=D10V3/SM2+D11V2
      D00V2Z=(D00V4-RM2*D01V1)/(RM2+2.D0*SM2)
C
      QZS(2)=QZS(2)+D10V1Z
      WZS(2)=WZS(2)+D00V2Z
C
   27 DC1=HD1*CDP2Z
      DC2=HD1*CDS2Z+ZD1*CDP1Z
      DC3=ZD1*DC1
      DCDPR=CDP1Z*S011-DC1*S021
      DCDSR=CDS1Z*S011-DC2*S021+DC3*S031
      D01V1=DCDPR+DCDSR
      D00V4=SM2*((X2-1.0)*DCDPR-2.0*DCDSR)
C
      DCDPR=CDP1Z*S101-DC1*S111
      DCDSR=CDS1Z*S101-DC2*S111+DC3*S121
      D10V1=DCDPR+DCDSR
C
      DCDPR=CDP1Z*S001-DC1*S011
      DCDSR=CDS1Z*S001-DC2*S011+DC3*S021
      D00V2=-(X2-1.0)*DCDPR+DCDSR
C
      DCDPR=CDP1Z*S111-DC1*S121
      DCDSR=CDS1Z*S111-DC2*S121+DC3*S131
      D11V2=-(X2-1.0)*DCDPR+DCDSR
      D10V3=SM2*((X2-3.0)*DCDPR-2.0*DCDSR)
C
      QZS(0)=QZS(0)+D10V1
      WZS(0)=WZS(0)+D00V2
C
      D10V1R=D01V1-D10V1/R
      D00V2R=-D11V2
C
      QZS(1)=QZS(1)+D10V1R
      WZS(1)=WZS(1)+D00V2R
C
      D10V1Z=D10V3/SM2+D11V2
      D00V2Z=(D00V4-RM2*D01V1)/(RM2+2.D0*SM2)
C
      QZS(2)=QZS(2)+D10V1Z
      WZS(2)=WZS(2)+D00V2Z
C
      GOTO 50
C
C *** SDP < (OR =) ODP < DEP2 ***
C
C  ** QX **
C
C   * SH (TRANSMITTED) *
C
   30 D1MH1=H1T*S1M
      DD0H1=H1T*(S00-S1M/R)
      QXS(0)=D1MH1/R
      VXS(0)=DD0H1
C
      D11H1=H1T*S11
      D1MH1R=DD0H1
      DD0H1R=(-DD0H1+D1MH1/R)/R-D11H1
      QXS(1)=(-D1MH1/R+D1MH1R)/R
      VXS(1)=DD0H1R
C
      D1MH2=H2T*S10
      DD0H2=H2T*(S01-S10/R)
      D1MH1Z=D1MH2/SM3
      DD0H1Z=DD0H2/SM3
      QXS(2)=D1MH1Z/R
      VXS(2)=DD0H1Z
C
C   * P-SV (TRANSMITTED) *
C
      D2H=DEP2-SDP
      ZD2=ODP-DEP2
      DC1=D2H*CDP2X
      DC2=D2H*CDS2X+ZD2*CDP1X
      DC3=ZD2*DC1
      DCDPT=CDP1X*S1M-DC1*S10
      DCDST=CDS1X*S1M-DC2*S10+DC3*S11
      D1MV1=DCDPT+DCDST
C
      DCDPT=CDP1X*S11-DC1*S12
      DCDST=CDS1X*S11-DC2*S12+DC3*S13
      D11V1=DCDPT+DCDST
      D10V4=SM3*((X3-1.0)*DCDPT-2.0*DCDST)
C
      S1=S00-S1M/R
      S2=S01-S10/R
      S3=S02-S11/R
      DCDPT=CDP1X*S1-DC1*S2
      DCDST=CDS1X*S1-DC2*S2+DC3*S3
      DD0V1=DCDPT+DCDST
C
      DCDPT=CDP1X*S01-DC1*S02
      DCDST=CDS1X*S01-DC2*S02+DC3*S03
      D01V2=-(X3-1.0)*DCDPT+DCDST
C
      DCDPT=CDP1X*S10-DC1*S11
      DCDST=CDS1X*S10-DC2*S11+DC3*S12
      D10V2=-(X3-1.0)*DCDPT+DCDST
      D1MV3=SM3*((X3-3.0)*DCDPT-2.0*DCDST)
C
      S1=S01-S10/R
      S2=S02-S11/R
      S3=S03-S12/R
      DCDPT=CDP1X*S1-DC1*S2
      DCDST=CDS1X*S1-DC2*S2+DC3*S3
      DD0V3=SM3*((X3-3.0)*DCDPT-2.0*DCDST)
C
      QXS(0)=QXS(0)+DD0V1
      VXS(0)=VXS(0)+D1MV1/R
      WXS(0)=D10V2
C
      D1MV1R=DD0V1
      DD0V1R=(-DD0V1+D1MV1/R)/R-D11V1
      D10V2R=D01V2-D10V2/R
C
      QXS(1)=QXS(1)+DD0V1R
      VXS(1)=VXS(1)+(-D1MV1/R+D1MV1R)/R
      WXS(1)=D10V2R
C
      D1MV1Z=D1MV3/SM3+D10V2
      DD0V1Z=DD0V3/SM3+D01V2-D10V2/R
      D10V2Z=(D10V4-RM3*D11V1)/(RM3+2.D0*SM3) 
C
      QXS(2)=QXS(2)+DD0V1Z
      VXS(2)=VXS(2)+D1MV1Z/R
      WXS(2)=D10V2Z
C
C  ** QZ **
C
      DC1=D2H*CDP2Z
      DC2=D2H*CDS2Z+ZD2*CDP1Z
      DC3=ZD2*DC1
      DCDPT=CDP1Z*S01-DC1*S02
      DCDST=CDS1Z*S01-DC2*S02+DC3*S03
      D01V1=DCDPT+DCDST
      D00V4=SM3*((X3-1.0)*DCDPT-2.0*DCDST)
C
      DCDPT=CDP1Z*S10-DC1*S11
      DCDST=CDS1Z*S10-DC2*S11+DC3*S12
      D10V1=DCDPT+DCDST
C
      DCDPT=CDP1Z*S00-DC1*S01
      DCDST=CDS1Z*S00-DC2*S01+DC3*S02
      D00V2=-(X3-1.0)*DCDPT+DCDST
C
      DCDPT=CDP1Z*S11-DC1*S12
      DCDST=CDS1Z*S11-DC2*S12+DC3*S13
      D11V2=-(X3-1.0)*DCDPT+DCDST
      D10V3=SM3*((X3-3.0)*DCDPT-2.0*DCDST)
C
      QZS(0)=D10V1
      WZS(0)=D00V2
C
      D10V1R=D01V1-D10V1/R
      D00V2R=-D11V2
C
      QZS(1)=D10V1R
      WZS(1)=D00V2R
C
      D10V1Z=D10V3/SM3+D11V2
      D00V2Z=(D00V4-RM3*D01V1)/(RM3+2.D0*SM3)
C
      QZS(2)=D10V1Z
      WZS(2)=D00V2Z
C
   50 RETURN
      END
C
C
C *** DISP-STRESS VECTOR FOR SH WAVES USING R/T MATRIX ***
C
      SUBROUTINE RTSH(OM,WN,VLX)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VP,VS,CD,CU,VS2,WS2,CMS,SM,EE(2,2)
     *         ,CUSU,CDSU,CDSD,CUSD,DD,DU,VLX(2)
      COMMON /SOLY/SDEP,ODEP,LS,LO
     *    /STRC/DNS(20),VP(20),VS(20),THK(20),NL
      COMMON /DUSH/CD(20),CU(20),CUSU,CDSU,CDSD,CUSD
C
C *** CHECK FOR THE DIMENSION OF CD AND CU ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(*,*)' FROM RTSH: NL.GT.NC ',NL,NC
         STOP
      END IF
C
C *** COMPUTATION OF CD AND CU ***
C
      CALL CDUSH(OM,WN)
C
C *** COMPUTATION OF VLX ***
C
      OM2=OM*OM
      WN2=WN*WN
C
      VS2=VS(LO)*VS(LO)
      WS2=OM2/VS2
      CMS=CDSQRT(WN2-WS2)
      IF(DREAL(CMS).LT.0.D0)CMS=-CMS
C
      SM=DNS(LO)*VS2
      EE(1,1)=(1.D0,0.D0)
      EE(1,2)=(1.D0,0.D0)
      EE(2,1)=-SM*CMS
      EE(2,2)= SM*CMS
C
      IF(LO.EQ.NL)THEN
C        
        IF(LO.NE.LS)THEN
          DD=CDEXP(-CMS*ODEP)
          VLX(1)=EE(1,1)*DD*CD(LO)
          VLX(2)=EE(2,1)*DD*CD(LO)
        ELSE
          IF(ODEP.GE.SDEP)THEN
            DD=CDEXP(-CMS*(ODEP-SDEP))
            VLX(1)=EE(1,1)*DD*CDSD
            VLX(2)=EE(2,1)*DD*CDSD
          ELSE
            DD=CDEXP(-CMS*ODEP)
            DU=CDEXP(-CMS*(SDEP-ODEP))
            VLX(1)=EE(1,1)*DD*CDSU+EE(1,2)*DU*CUSU
            VLX(2)=EE(2,1)*DD*CDSU+EE(2,2)*DU*CUSU
          END IF
        END IF
C
      ELSE
C        
        IF(LO.NE.LS)THEN
          DD=CDEXP(-CMS*ODEP)
          DU=CDEXP(-CMS*(THK(LO)-ODEP))
          VLX(1)=EE(1,1)*DD*CD(LO)+EE(1,2)*DU*CU(LO)
          VLX(2)=EE(2,1)*DD*CD(LO)+EE(2,2)*DU*CU(LO)
        ELSE
          IF(ODEP.GE.SDEP)THEN
            DD=CDEXP(-CMS*(ODEP-SDEP))
            DU=CDEXP(-CMS*(THK(LO)-ODEP))
            VLX(1)=EE(1,1)*DD*CDSD+EE(1,2)*DU*CUSD
            VLX(2)=EE(2,1)*DD*CDSD+EE(2,2)*DU*CUSD
          ELSE
            DD=CDEXP(-CMS*ODEP)
            DU=CDEXP(-CMS*(SDEP-ODEP))
            VLX(1)=EE(1,1)*DD*CDSU+EE(1,2)*DU*CUSU
            VLX(2)=EE(2,1)*DD*CDSU+EE(2,2)*DU*CUSU
          END IF
        END IF
C
      END IF
C
      RETURN
      END
C
C
C *** SH UPGOING/DOWNGOING COEFFICIENTS USING R/T MATRIX ***
C
      SUBROUTINE CDUSH(OM,WN)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VP,VS,CD,CU,VS2,WS2,CMS,SM,EE(2,2),DD,DU
     *      ,CUSU,CDSU,CDSD,CUSD,GTD,GTU,GRD,GRU,AU,BU,AD,BD,CS
      COMMON /SOLY/SDEP,ODEP,LS,LO
     *    /STRC/DNS(20),VP(20),VS(20),THK(20),NL
      COMMON /DUSH/CD(20),CU(20),CUSU,CDSU,CDSD,CUSD 
     *    /GRTSH/GTD(0:20),GTU(0:20),GRD(0:20),GRU(0:20)
C
C *** CHECK FOR THE DIMENSION OF CD, CU, AND GTD-GRU ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(*,*)' FROM CDUSH: NL.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C *** COMPUTATION OF GENERALIZED R/T MATRICES ***
C
      CALL GNRTSH(OM,WN)
C
C *** COMPUTATION OF  CD AND CU ***
C
      OM2=OM*OM
      WN2=WN*WN
C
      VS2=VS(LS)*VS(LS)
      WS2=OM2/VS2
      CMS=CDSQRT(WN2-WS2)
      IF(DREAL(CMS).LT.0.D0)CMS=-CMS
C
      SM=DNS(LS)*VS2
      EE(1,1)=(1.D0,0.D0)
      EE(1,2)=(1.D0,0.D0)
      EE(2,1)=-SM*CMS
      EE(2,2)= SM*CMS
C
      DQ=-WN/(2.D0*3.1415926536D0)
C
      IF((LS.EQ.1).AND.(SDEP.EQ.0.D0))THEN
C
        IF(NL.EQ.1)THEN
          CDSD=DQ/EE(2,1)
          CUSD=(0.D0,0.D0)
        ELSE
          DD=CDEXP(-CMS*THK(1))
          CDSD=DQ/(EE(2,1)+EE(2,2)*DD*GRD(1))
          CUSD=GRD(1)*CDSD
        END IF
C
      ELSE
C
        IF(LS.EQ.NL)THEN 
C
          DD=CDEXP(-CMS*SDEP)
C
          AU=EE(1,1)*DD*GRU(LS-1)+EE(1,2)
          BU=EE(2,1)*DD*GRU(LS-1)+EE(2,2)
          AD=EE(1,1)
          BD=EE(2,1)
          CS=AU/AD
C
          CUSU=DQ/(BD*CS-BU)
          CDSU=GRU(LS-1)*CUSU
          CDSD=CS*CUSU
C
        ELSE
C
          DD=CDEXP(-CMS*SDEP)
          DU=CDEXP(-CMS*(THK(LS)-SDEP))
C
          AU=EE(1,1)*DD*GRU(LS-1)+EE(1,2)
          BU=EE(2,1)*DD*GRU(LS-1)+EE(2,2)
          AD=EE(1,1)+EE(1,2)*DU*GRD(LS)
          BD=EE(2,1)+EE(2,2)*DU*GRD(LS)
          CS=AU/AD
C
          CUSU=DQ/(BD*CS-BU)
          CDSU=GRU(LS-1)*CUSU
          CDSD=CS*CUSU
          CUSD=GRD(LS)*CDSD
C
        ENDIF
C
      END IF
C
      IF(LS.EQ.1)GOTO 10
      CU(LS-1)=GTU(LS-1)*CUSU
      DO 100 I=LS-2,1,-1
        CU(I)=GTU(I)*CU(I+1)
        CD(I+1)=GRU(I)*CU(I+1)
  100 CONTINUE
      CD(1)=GRU(0)*CU(1)
C
   10 IF(LS.EQ.NL)GOTO 20
      CD(LS+1)=GTD(LS)*CDSD
      DO 200 I=LS+1,NL-1
        CD(I+1)=GTD(I)*CD(I)
        CU(I)=GRD(I)*CD(I)
  200 CONTINUE
C
   20 RETURN
      END
C
C
C
C *** SH GENERALIZED R/T MATIX ***
C
      SUBROUTINE GNRTSH(OM,WN)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VP,VS,GTD,GTU,GRD,GRU,TD,TU,RD,RU
      COMMON /SOLY/SDEP,ODEP,LS,LO
     *    /STRC/DNS(20),VP(20),VS(20),THK(20),NL
      COMMON /GRTSH/GTD(0:20),GTU(0:20),GRD(0:20),GRU(0:20)
     *    /MRTSH/TD(0:20),TU(0:20),RD(0:20),RU(0:20)
C
C *** CHECK FOR THE DIMENSION OF GTD, GTU, GRD, AND GRU ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(*,*)' FROM GNRTSH: NL.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C *** MODIFIED R/T MATRIX ***
C
      CALL MDRTSH(OM,WN)
C
C *** GENERALIZED R/T MATRIX ***
C
      GRU(0)=RU(0)
      DO 100 I=1,LS-1
        GTU(I)=TU(I)/(1.D0-RD(I)*GRU(I-1))
        GRU(I)=RU(I)+TD(I)*GRU(I-1)*GTU(I)
  100 CONTINUE
C
      GTD(NL-1)=TD(NL-1)
      GRD(NL-1)=RD(NL-1)
      DO 200 I=NL-2,LS,-1
        GTD(I)=TD(I)/(1.D0-RU(I)*GRD(I+1))
        GRD(I)=RD(I)+TU(I)*GRD(I+1)*GTD(I)
  200 CONTINUE
C
      RETURN
      END
C
C
C
C *** SH MODIFIED R/T MATIX ***
C
      SUBROUTINE MDRTSH(OM,WN)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VP,VS,VS2,WS2,CMS,SM,EA(2,2),EB(2,2)
     *      ,TD,TU,RD,RU,DDA,DUB,EAB
      COMMON /SOLY/SDEP,ODEP,LS,LO
     *    /STRC/DNS(20),VP(20),VS(20),THK(20),NL
      COMMON /MRTSH/TD(0:20),TU(0:20),RD(0:20),RU(0:20)
C
C *** CHECK FOR THE DIMENSION OF TD, TU, RD, AND RU ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(*,*)' FROM MDRTSH: NL.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C *** MODIFIED R/T MATRIX ***
C
      WN2=WN*WN
      OM2=OM*OM
C
      VS2=VS(1)*VS(1)
      WS2=OM2/VS2
      CMS=CDSQRT(WN2-WS2)
      IF(DREAL(CMS).LT.0.D0)CMS=-CMS
C
      SM=DNS(1)*VS2
      EA(1,1)=(1.D0,0.D0)
      EA(1,2)=(1.D0,0.D0)
      EA(2,1)=-SM*CMS
      EA(2,2)= SM*CMS
C
      IF(NL.EQ.1)THEN
        DUB=CDEXP(-CMS*SDEP)
        RU(0)=-EA(2,2)*DUB/EA(2,1)
        RETURN
      END IF
C
      IF(LS.EQ.1)THEN
        DUB=CDEXP(-CMS*SDEP)
        DDA=CDEXP(-CMS*(THK(1)-SDEP))
      ELSE
        DUB=CDEXP(-CMS*THK(1))
        DDA=DUB
      END IF
C
      RU(0)=-EA(2,2)*DUB/EA(2,1)
C
      DO 100 I=1, NL-2
C
        VS2=VS(I+1)*VS(I+1)
        WS2=OM2/VS2
        CMS=CDSQRT(WN2-WS2)
        IF(DREAL(CMS).LT.0.D0)CMS=-CMS
C
        SM=DNS(I+1)*VS2
        EB(1,1)=(1.D0,0.D0)
        EB(1,2)=(1.D0,0.D0)
        EB(2,1)=-SM*CMS
        EB(2,2)= SM*CMS
C
        IF(LS.EQ.(I+1))THEN
          DUB=CDEXP(-CMS*SDEP)
        ELSE
          DUB=CDEXP(-CMS*THK(I+1))
        END IF
C
        EAB=EB(2,1)*EA(1,2)-EB(1,1)*EA(2,2) 
        TD(I)=(EA(1,2)*EA(2,1)-EA(1,1)*EA(2,2))*DDA/EAB
        RU(I)=(EB(1,2)*EA(2,2)-EB(2,2)*EA(1,2))*DUB/EAB
        RD(I)=(EB(1,1)*EA(2,1)-EA(1,1)*EB(2,1))*DDA/EAB
        TU(I)=(EB(1,2)*EB(2,1)-EB(2,2)*EB(1,1))*DUB/EAB
C
        EA(1,1)=EB(1,1)
        EA(1,2)=EB(1,2)
        EA(2,1)=EB(2,1)
        EA(2,2)=EB(2,2)
C
        IF(LS.EQ.(I+1))THEN
          DDA=CDEXP(-CMS*(THK(I+1)-SDEP))
        ELSE          
          DDA=DUB
        END IF
C
  100 CONTINUE
C        
      VS2=VS(NL)*VS(NL)
      WS2=OM2/VS2
      CMS=CDSQRT(WN2-WS2)
      IF(DREAL(CMS).LT.0.D0)CMS=-CMS
C
      SM=DNS(NL)*VS2
      EB(1,1)=(1.D0,0.D0)
      EB(1,2)=(1.D0,0.D0)
      EB(2,1)=-SM*CMS
      EB(2,2)= SM*CMS
C
      EAB=EB(2,1)*EA(1,2)-EB(1,1)*EA(2,2) 
C
      IF(LS.EQ.NL)THEN
C
        DUB=CDEXP(-CMS*SDEP)
        TD(NL-1)=(EA(1,2)*EA(2,1)-EA(1,1)*EA(2,2))*DDA/EAB
        RU(NL-1)=(EB(1,2)*EA(2,2)-EB(2,2)*EA(1,2))*DUB/EAB
        RD(NL-1)=(EB(1,1)*EA(2,1)-EA(1,1)*EB(2,1))*DDA/EAB
        TU(NL-1)=(EB(1,2)*EB(2,1)-EB(2,2)*EB(1,1))*DUB/EAB
C
      ELSE
C
        TD(NL-1)=(EA(1,2)*EA(2,1)-EA(1,1)*EA(2,2))*DDA/EAB
        RD(NL-1)=(EB(1,1)*EA(2,1)-EA(1,1)*EB(2,1))*DDA/EAB
C
      END IF
C
      RETURN
      END
C
C
C *** PSV DISPLACEMENT USING R/T MATRIX ***
C
      SUBROUTINE RTPSV(OM,WN,URX,URZ)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VP,VS,CMP,CMS,EE(4,4),DD(2,2),DU(2,2)
     *    ,CDX,CUX,CUSUX,CDSUX,CDSDX,CUSDX,URX(2,2),E1(2,2)
     *    ,CDZ,CUZ,CUSUZ,CDSUZ,CDSDZ,CUSDZ,URZ(2,2),E2(2,2)
      COMMON /SOLY/SDEP,ODEP,LS,LO
     * /STRC/DNS(20),VP(20),VS(20),THK(20),NL
     * /DUPSVX/CDX(2,20),CUX(2,20),CUSUX(2),CDSUX(2),CDSDX(2),CUSDX(2)
     * /DUPSVZ/CDZ(2,20),CUZ(2,20),CUSUZ(2),CDSUZ(2),CDSDZ(2),CUSDZ(2)
C
C *** CHECK FOR THE DIMENSION OF CD AND CU ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(*,*)' FROM RTPSV: NL.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C *** COMPUTATION OF CD AND CU ***
C
      CALL CDUPSV(OM,WN)
C
C *** COMPUTATION OF URX AND URZ ***
C
      CALL MTXPSV(OM,WN,LO,CMS,CMP,EE)
C
      IF(LO.EQ.NL)THEN
C        
        IF(LO.NE.LS)THEN
C
          CALL DDUPSV(OM,WN,ODEP,CMS,CMP,DD)
C
          E1(1,1)=EE(1,1)*DD(1,1)+EE(1,2)*DD(2,1)
          E1(1,2)=EE(1,1)*DD(1,2)+EE(1,2)*DD(2,2)
          E1(2,1)=EE(2,1)*DD(1,1)+EE(2,2)*DD(2,1)
          E1(2,2)=EE(2,1)*DD(1,2)+EE(2,2)*DD(2,2)
          URX(1,1)=E1(1,1)*CDX(1,LO)+E1(1,2)*CDX(2,LO)
          URX(2,1)=E1(2,1)*CDX(1,LO)+E1(2,2)*CDX(2,LO)
          URZ(1,1)=E1(1,1)*CDZ(1,LO)+E1(1,2)*CDZ(2,LO)
          URZ(2,1)=E1(2,1)*CDZ(1,LO)+E1(2,2)*CDZ(2,LO)
C
          E1(1,1)=EE(3,1)*DD(1,1)+EE(3,2)*DD(2,1)
          E1(1,2)=EE(3,1)*DD(1,2)+EE(3,2)*DD(2,2)
          E1(2,1)=EE(4,1)*DD(1,1)+EE(4,2)*DD(2,1)
          E1(2,2)=EE(4,1)*DD(1,2)+EE(4,2)*DD(2,2)
          URX(1,2)=E1(1,1)*CDX(1,LO)+E1(1,2)*CDX(2,LO)
          URX(2,2)=E1(2,1)*CDX(1,LO)+E1(2,2)*CDX(2,LO)
          URZ(1,2)=E1(1,1)*CDZ(1,LO)+E1(1,2)*CDZ(2,LO)
          URZ(2,2)=E1(2,1)*CDZ(1,LO)+E1(2,2)*CDZ(2,LO)
C
        ELSE
C
          IF(ODEP.GE.SDEP)THEN
C
            CALL DDUPSV(OM,WN,ODEP-SDEP,CMS,CMP,DD)
C
            E1(1,1)=EE(1,1)*DD(1,1)+EE(1,2)*DD(2,1)
            E1(1,2)=EE(1,1)*DD(1,2)+EE(1,2)*DD(2,2)
            E1(2,1)=EE(2,1)*DD(1,1)+EE(2,2)*DD(2,1)
            E1(2,2)=EE(2,1)*DD(1,2)+EE(2,2)*DD(2,2)
            URX(1,1)=E1(1,1)*CDSDX(1)+E1(1,2)*CDSDX(2)
            URX(2,1)=E1(2,1)*CDSDX(1)+E1(2,2)*CDSDX(2)
            URZ(1,1)=E1(1,1)*CDSDZ(1)+E1(1,2)*CDSDZ(2)
            URZ(2,1)=E1(2,1)*CDSDZ(1)+E1(2,2)*CDSDZ(2)
C
            E1(1,1)=EE(3,1)*DD(1,1)+EE(3,2)*DD(2,1)
            E1(1,2)=EE(3,1)*DD(1,2)+EE(3,2)*DD(2,2)
            E1(2,1)=EE(4,1)*DD(1,1)+EE(4,2)*DD(2,1)
            E1(2,2)=EE(4,1)*DD(1,2)+EE(4,2)*DD(2,2)
            URX(1,2)=E1(1,1)*CDSDX(1)+E1(1,2)*CDSDX(2)
            URX(2,2)=E1(2,1)*CDSDX(1)+E1(2,2)*CDSDX(2)
            URZ(1,2)=E1(1,1)*CDSDZ(1)+E1(1,2)*CDSDZ(2)
            URZ(2,2)=E1(2,1)*CDSDZ(1)+E1(2,2)*CDSDZ(2)
C
          ELSE
C
            CALL DDUPSV(OM,WN,ODEP,CMS,CMP,DD)
            CALL DDUPSV(OM,WN,SDEP-ODEP,CMS,CMP,DU)
C
            E1(1,1)=EE(1,1)*DD(1,1)+EE(1,2)*DD(2,1)
            E1(1,2)=EE(1,1)*DD(1,2)+EE(1,2)*DD(2,2)
            E1(2,1)=EE(2,1)*DD(1,1)+EE(2,2)*DD(2,1)
            E1(2,2)=EE(2,1)*DD(1,2)+EE(2,2)*DD(2,2)
            E2(1,1)=EE(1,3)*DU(1,1)+EE(1,4)*DU(2,1)
            E2(1,2)=EE(1,3)*DU(1,2)+EE(1,4)*DU(2,2)
            E2(2,1)=EE(2,3)*DU(1,1)+EE(2,4)*DU(2,1)
            E2(2,2)=EE(2,3)*DU(1,2)+EE(2,4)*DU(2,2)
C
            URX(1,1)=E1(1,1)*CDSUX(1)+E1(1,2)*CDSUX(2)
     *            +E2(1,1)*CUSUX(1)+E2(1,2)*CUSUX(2)
            URX(2,1)=E1(2,1)*CDSUX(1)+E1(2,2)*CDSUX(2)
     *            +E2(2,1)*CUSUX(1)+E2(2,2)*CUSUX(2)
            URZ(1,1)=E1(1,1)*CDSUZ(1)+E1(1,2)*CDSUZ(2)
     *            +E2(1,1)*CUSUZ(1)+E2(1,2)*CUSUZ(2)
            URZ(2,1)=E1(2,1)*CDSUZ(1)+E1(2,2)*CDSUZ(2)
     *            +E2(2,1)*CUSUZ(1)+E2(2,2)*CUSUZ(2)
C
            E1(1,1)=EE(3,1)*DD(1,1)+EE(3,2)*DD(2,1)
            E1(1,2)=EE(3,1)*DD(1,2)+EE(3,2)*DD(2,2)
            E1(2,1)=EE(4,1)*DD(1,1)+EE(4,2)*DD(2,1)
            E1(2,2)=EE(4,1)*DD(1,2)+EE(4,2)*DD(2,2)
            E2(1,1)=EE(3,3)*DU(1,1)+EE(3,4)*DU(2,1)
            E2(1,2)=EE(3,3)*DU(1,2)+EE(3,4)*DU(2,2)
            E2(2,1)=EE(4,3)*DU(1,1)+EE(4,4)*DU(2,1)
            E2(2,2)=EE(4,3)*DU(1,2)+EE(4,4)*DU(2,2)
C
            URX(1,2)=E1(1,1)*CDSUX(1)+E1(1,2)*CDSUX(2)
     *            +E2(1,1)*CUSUX(1)+E2(1,2)*CUSUX(2)
            URX(2,2)=E1(2,1)*CDSUX(1)+E1(2,2)*CDSUX(2)
     *            +E2(2,1)*CUSUX(1)+E2(2,2)*CUSUX(2)
            URZ(1,2)=E1(1,1)*CDSUZ(1)+E1(1,2)*CDSUZ(2)
     *            +E2(1,1)*CUSUZ(1)+E2(1,2)*CUSUZ(2)
            URZ(2,2)=E1(2,1)*CDSUZ(1)+E1(2,2)*CDSUZ(2)
     *            +E2(2,1)*CUSUZ(1)+E2(2,2)*CUSUZ(2)
C
          END IF
C
        END IF
C
      ELSE
C        
        IF(LO.NE.LS)THEN
C
            CALL DDUPSV(OM,WN,ODEP,CMS,CMP,DD)
            CALL DDUPSV(OM,WN,THK(LO)-ODEP,CMS,CMP,DU)
C
            E1(1,1)=EE(1,1)*DD(1,1)+EE(1,2)*DD(2,1)
            E1(1,2)=EE(1,1)*DD(1,2)+EE(1,2)*DD(2,2)
            E1(2,1)=EE(2,1)*DD(1,1)+EE(2,2)*DD(2,1)
            E1(2,2)=EE(2,1)*DD(1,2)+EE(2,2)*DD(2,2)
            E2(1,1)=EE(1,3)*DU(1,1)+EE(1,4)*DU(2,1)
            E2(1,2)=EE(1,3)*DU(1,2)+EE(1,4)*DU(2,2)
            E2(2,1)=EE(2,3)*DU(1,1)+EE(2,4)*DU(2,1)
            E2(2,2)=EE(2,3)*DU(1,2)+EE(2,4)*DU(2,2)
C
            URX(1,1)=E1(1,1)*CDX(1,LO)+E1(1,2)*CDX(2,LO)
     *            +E2(1,1)*CUX(1,LO)+E2(1,2)*CUX(2,LO)
            URX(2,1)=E1(2,1)*CDX(1,LO)+E1(2,2)*CDX(2,LO)
     *            +E2(2,1)*CUX(1,LO)+E2(2,2)*CUX(2,LO)
            URZ(1,1)=E1(1,1)*CDZ(1,LO)+E1(1,2)*CDZ(2,LO)
     *            +E2(1,1)*CUZ(1,LO)+E2(1,2)*CUZ(2,LO)
            URZ(2,1)=E1(2,1)*CDZ(1,LO)+E1(2,2)*CDZ(2,LO)
     *            +E2(2,1)*CUZ(1,LO)+E2(2,2)*CUZ(2,LO)
C
            E1(1,1)=EE(3,1)*DD(1,1)+EE(3,2)*DD(2,1)
            E1(1,2)=EE(3,1)*DD(1,2)+EE(3,2)*DD(2,2)
            E1(2,1)=EE(4,1)*DD(1,1)+EE(4,2)*DD(2,1)
            E1(2,2)=EE(4,1)*DD(1,2)+EE(4,2)*DD(2,2)
            E2(1,1)=EE(3,3)*DU(1,1)+EE(3,4)*DU(2,1)
            E2(1,2)=EE(3,3)*DU(1,2)+EE(3,4)*DU(2,2)
            E2(2,1)=EE(4,3)*DU(1,1)+EE(4,4)*DU(2,1)
            E2(2,2)=EE(4,3)*DU(1,2)+EE(4,4)*DU(2,2)
C
            URX(1,2)=E1(1,1)*CDX(1,LO)+E1(1,2)*CDX(2,LO)
     *            +E2(1,1)*CUX(1,LO)+E2(1,2)*CUX(2,LO)
            URX(2,2)=E1(2,1)*CDX(1,LO)+E1(2,2)*CDX(2,LO)
     *            +E2(2,1)*CUX(1,LO)+E2(2,2)*CUX(2,LO)
            URZ(1,2)=E1(1,1)*CDZ(1,LO)+E1(1,2)*CDZ(2,LO)
     *            +E2(1,1)*CUZ(1,LO)+E2(1,2)*CUZ(2,LO)
            URZ(2,2)=E1(2,1)*CDZ(1,LO)+E1(2,2)*CDZ(2,LO)
     *            +E2(2,1)*CUZ(1,LO)+E2(2,2)*CUZ(2,LO)
C
        ELSE
C
          IF(ODEP.GE.SDEP)THEN
C
            CALL DDUPSV(OM,WN,ODEP-SDEP,CMS,CMP,DD)
            CALL DDUPSV(OM,WN,THK(LO)-ODEP,CMS,CMP,DU)
C
            E1(1,1)=EE(1,1)*DD(1,1)+EE(1,2)*DD(2,1)
            E1(1,2)=EE(1,1)*DD(1,2)+EE(1,2)*DD(2,2)
            E1(2,1)=EE(2,1)*DD(1,1)+EE(2,2)*DD(2,1)
            E1(2,2)=EE(2,1)*DD(1,2)+EE(2,2)*DD(2,2)
            E2(1,1)=EE(1,3)*DU(1,1)+EE(1,4)*DU(2,1)
            E2(1,2)=EE(1,3)*DU(1,2)+EE(1,4)*DU(2,2)
            E2(2,1)=EE(2,3)*DU(1,1)+EE(2,4)*DU(2,1)
            E2(2,2)=EE(2,3)*DU(1,2)+EE(2,4)*DU(2,2)
C
            URX(1,1)=E1(1,1)*CDSDX(1)+E1(1,2)*CDSDX(2)
     *            +E2(1,1)*CUSDX(1)+E2(1,2)*CUSDX(2)
            URX(2,1)=E1(2,1)*CDSDX(1)+E1(2,2)*CDSDX(2)
     *            +E2(2,1)*CUSDX(1)+E2(2,2)*CUSDX(2)
            URZ(1,1)=E1(1,1)*CDSDZ(1)+E1(1,2)*CDSDZ(2)
     *            +E2(1,1)*CUSDZ(1)+E2(1,2)*CUSDZ(2)
            URZ(2,1)=E1(2,1)*CDSDZ(1)+E1(2,2)*CDSDZ(2)
     *            +E2(2,1)*CUSDZ(1)+E2(2,2)*CUSDZ(2)
C
            E1(1,1)=EE(3,1)*DD(1,1)+EE(3,2)*DD(2,1)
            E1(1,2)=EE(3,1)*DD(1,2)+EE(3,2)*DD(2,2)
            E1(2,1)=EE(4,1)*DD(1,1)+EE(4,2)*DD(2,1)
            E1(2,2)=EE(4,1)*DD(1,2)+EE(4,2)*DD(2,2)
            E2(1,1)=EE(3,3)*DU(1,1)+EE(3,4)*DU(2,1)
            E2(1,2)=EE(3,3)*DU(1,2)+EE(3,4)*DU(2,2)
            E2(2,1)=EE(4,3)*DU(1,1)+EE(4,4)*DU(2,1)
            E2(2,2)=EE(4,3)*DU(1,2)+EE(4,4)*DU(2,2)
C
            URX(1,2)=E1(1,1)*CDSDX(1)+E1(1,2)*CDSDX(2)
     *            +E2(1,1)*CUSDX(1)+E2(1,2)*CUSDX(2)
            URX(2,2)=E1(2,1)*CDSDX(1)+E1(2,2)*CDSDX(2)
     *            +E2(2,1)*CUSDX(1)+E2(2,2)*CUSDX(2)
            URZ(1,2)=E1(1,1)*CDSDZ(1)+E1(1,2)*CDSDZ(2)
     *            +E2(1,1)*CUSDZ(1)+E2(1,2)*CUSDZ(2)
            URZ(2,2)=E1(2,1)*CDSDZ(1)+E1(2,2)*CDSDZ(2)
     *            +E2(2,1)*CUSDZ(1)+E2(2,2)*CUSDZ(2)
C
          ELSE
C
            CALL DDUPSV(OM,WN,ODEP,CMS,CMP,DD)
            CALL DDUPSV(OM,WN,SDEP-ODEP,CMS,CMP,DU)
C
            E1(1,1)=EE(1,1)*DD(1,1)+EE(1,2)*DD(2,1)
            E1(1,2)=EE(1,1)*DD(1,2)+EE(1,2)*DD(2,2)
            E1(2,1)=EE(2,1)*DD(1,1)+EE(2,2)*DD(2,1)
            E1(2,2)=EE(2,1)*DD(1,2)+EE(2,2)*DD(2,2)
            E2(1,1)=EE(1,3)*DU(1,1)+EE(1,4)*DU(2,1)
            E2(1,2)=EE(1,3)*DU(1,2)+EE(1,4)*DU(2,2)
            E2(2,1)=EE(2,3)*DU(1,1)+EE(2,4)*DU(2,1)
            E2(2,2)=EE(2,3)*DU(1,2)+EE(2,4)*DU(2,2)
C
            URX(1,1)=E1(1,1)*CDSUX(1)+E1(1,2)*CDSUX(2)
     *            +E2(1,1)*CUSUX(1)+E2(1,2)*CUSUX(2)
            URX(2,1)=E1(2,1)*CDSUX(1)+E1(2,2)*CDSUX(2)
     *            +E2(2,1)*CUSUX(1)+E2(2,2)*CUSUX(2)
            URZ(1,1)=E1(1,1)*CDSUZ(1)+E1(1,2)*CDSUZ(2)
     *            +E2(1,1)*CUSUZ(1)+E2(1,2)*CUSUZ(2)
            URZ(2,1)=E1(2,1)*CDSUZ(1)+E1(2,2)*CDSUZ(2)
     *            +E2(2,1)*CUSUZ(1)+E2(2,2)*CUSUZ(2)
C
            E1(1,1)=EE(3,1)*DD(1,1)+EE(3,2)*DD(2,1)
            E1(1,2)=EE(3,1)*DD(1,2)+EE(3,2)*DD(2,2)
            E1(2,1)=EE(4,1)*DD(1,1)+EE(4,2)*DD(2,1)
            E1(2,2)=EE(4,1)*DD(1,2)+EE(4,2)*DD(2,2)
            E2(1,1)=EE(3,3)*DU(1,1)+EE(3,4)*DU(2,1)
            E2(1,2)=EE(3,3)*DU(1,2)+EE(3,4)*DU(2,2)
            E2(2,1)=EE(4,3)*DU(1,1)+EE(4,4)*DU(2,1)
            E2(2,2)=EE(4,3)*DU(1,2)+EE(4,4)*DU(2,2)
C
            URX(1,2)=E1(1,1)*CDSUX(1)+E1(1,2)*CDSUX(2)
     *            +E2(1,1)*CUSUX(1)+E2(1,2)*CUSUX(2)
            URX(2,2)=E1(2,1)*CDSUX(1)+E1(2,2)*CDSUX(2)
     *            +E2(2,1)*CUSUX(1)+E2(2,2)*CUSUX(2)
            URZ(1,2)=E1(1,1)*CDSUZ(1)+E1(1,2)*CDSUZ(2)
     *            +E2(1,1)*CUSUZ(1)+E2(1,2)*CUSUZ(2)
            URZ(2,2)=E1(2,1)*CDSUZ(1)+E1(2,2)*CDSUZ(2)
     *            +E2(2,1)*CUSUZ(1)+E2(2,2)*CUSUZ(2)
C
          END IF
C
        END IF
C
      END IF
C
      RETURN
      END
C
C
C
C *** P-SV UPGOING/DOWNGOING COEFFICIENTS USING R/T MATRIX ***
C
      SUBROUTINE CDUPSV(OM,WN)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VP,VS,CMP,CMS,EE(4,4),DD(2,2),DU(2,2)
     *    ,CDX,CUX,CUSUX,CDSUX,CDSDX,CUSDX,GTD,GTU,GRD,GRU,DT
     *    ,CDZ,CUZ,CUSUZ,CDSUZ,CDSDZ,CUSDZ
     *    ,AU(2,2),BU(2,2),AD(2,2),BD(2,2),CS(2,2),E1(2,2),E2(2,2)
      COMMON /SOLY/SDEP,ODEP,LS,LO 
     *  /STRC/DNS(20),VP(20),VS(20),THK(20),NL
     *  /DUPSVX/CDX(2,20),CUX(2,20),CUSUX(2),CDSUX(2),CDSDX(2),CUSDX(2) 
     *  /DUPSVZ/CDZ(2,20),CUZ(2,20),CUSUZ(2),CDSUZ(2),CDSDZ(2),CUSDZ(2) 
     *  /GRTPSV/GTD(2,2,0:20),GTU(2,2,0:20),GRD(2,2,0:20),GRU(2,2,0:20)
C
C *** CHECK FOR THE DIMENSION OF CD, CU, AND GTD-GRU ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(*,*)' FROM CDUPSV: NL.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C *** COMPUTATION OF GENERALIZED R/T MATRICES ***
C
      CALL GNRTPSV(OM,WN)
C
C *** COMPUTATION OF CD AND CU ***
C
      CALL MTXPSV(OM,WN,LS,CMS,CMP,EE)
C
      DQX=-WN/(2.D0*3.1415926536D0)
      DQZ= WN/(2.D0*3.1415926536D0)
C
      IF((LS.EQ.1).AND.(SDEP.EQ.0.D0))THEN
C
        IF(NL.EQ.1)THEN
C
          DT=EE(3,1)*EE(4,2)-EE(3,2)*EE(4,1)
          E2(1,1)= EE(4,2)/DT
          E2(1,2)=-EE(3,2)/DT
          E2(2,1)=-EE(4,1)/DT
          E2(2,2)= EE(3,1)/DT
C
          CDSDX(1)=E2(1,1)*DQX
          CDSDX(2)=E2(2,1)*DQX
          CUSDX(1)=(0.D0,0.D0)
          CUSDX(2)=(0.D0,0.D0)
          CDSDZ(1)=E2(1,2)*DQZ
          CDSDZ(2)=E2(2,2)*DQZ
          CUSDZ(1)=(0.D0,0.D0)
          CUSDZ(2)=(0.D0,0.D0)
C
        ELSE
C
          CALL DDUPSV(OM,WN,THK(1),CMS,CMP,DU)
C
          E2(1,1)=EE(3,3)*DU(1,1)+EE(3,4)*DU(2,1)
          E2(1,2)=EE(3,3)*DU(1,2)+EE(3,4)*DU(2,2)
          E2(2,1)=EE(4,3)*DU(1,1)+EE(4,4)*DU(2,1)
          E2(2,2)=EE(4,3)*DU(1,2)+EE(4,4)*DU(2,2)
          E1(1,1)=EE(3,1)+E2(1,1)*GRD(1,1,1)+E2(1,2)*GRD(2,1,1)
          E1(1,2)=EE(3,2)+E2(1,1)*GRD(1,2,1)+E2(1,2)*GRD(2,2,1)
          E1(2,1)=EE(4,1)+E2(2,1)*GRD(1,1,1)+E2(2,2)*GRD(2,1,1)
          E1(2,2)=EE(4,2)+E2(2,1)*GRD(1,2,1)+E2(2,2)*GRD(2,2,1)
          DT=E1(1,1)*E1(2,2)-E1(1,2)*E1(2,1)
          E2(1,1)= E1(2,2)/DT
          E2(1,2)=-E1(1,2)/DT
          E2(2,1)=-E1(2,1)/DT
          E2(2,2)= E1(1,1)/DT
C
          CDSDX(1)=E2(1,1)*DQX
          CDSDX(2)=E2(2,1)*DQX
          CUSDX(1)=GRD(1,1,1)*CDSDX(1)+GRD(1,2,1)*CDSDX(2)
          CUSDX(2)=GRD(2,1,1)*CDSDX(1)+GRD(2,2,1)*CDSDX(2)
          CDSDZ(1)=E2(1,2)*DQZ
          CDSDZ(2)=E2(2,2)*DQZ
          CUSDZ(1)=GRD(1,1,1)*CDSDZ(1)+GRD(1,2,1)*CDSDZ(2)
          CUSDZ(2)=GRD(2,1,1)*CDSDZ(1)+GRD(2,2,1)*CDSDZ(2)
C
        END IF
C
      ELSE
C
        IF(LS.EQ.NL)THEN 
C
          CALL DDUPSV(OM,WN,SDEP,CMS,CMP,DD)
C
          E1(1,1)=DD(1,1)*GRU(1,1,LS-1)+DD(1,2)*GRU(2,1,LS-1)
          E1(1,2)=DD(1,1)*GRU(1,2,LS-1)+DD(1,2)*GRU(2,2,LS-1)
          E1(2,1)=DD(2,1)*GRU(1,1,LS-1)+DD(2,2)*GRU(2,1,LS-1)
          E1(2,2)=DD(2,1)*GRU(1,2,LS-1)+DD(2,2)*GRU(2,2,LS-1)
C
          AU(1,1)=EE(1,1)*E1(1,1)+EE(1,2)*E1(2,1)+EE(1,3)
          AU(1,2)=EE(1,1)*E1(1,2)+EE(1,2)*E1(2,2)+EE(1,4)
          AU(2,1)=EE(2,1)*E1(1,1)+EE(2,2)*E1(2,1)+EE(2,3)
          AU(2,2)=EE(2,1)*E1(1,2)+EE(2,2)*E1(2,2)+EE(2,4)
          BU(1,1)=EE(3,1)*E1(1,1)+EE(3,2)*E1(2,1)+EE(3,3)
          BU(1,2)=EE(3,1)*E1(1,2)+EE(3,2)*E1(2,2)+EE(3,4)
          BU(2,1)=EE(4,1)*E1(1,1)+EE(4,2)*E1(2,1)+EE(4,3)
          BU(2,2)=EE(4,1)*E1(1,2)+EE(4,2)*E1(2,2)+EE(4,4)
          AD(1,1)=EE(1,1)
          AD(1,2)=EE(1,2)
          AD(2,1)=EE(2,1)
          AD(2,2)=EE(2,2)
          BD(1,1)=EE(3,1)
          BD(1,2)=EE(3,2)
          BD(2,1)=EE(4,1)
          BD(2,2)=EE(4,2)
C
          DT=AD(1,1)*AD(2,2)-AD(1,2)*AD(2,1)
          E1(1,1)= AD(2,2)/DT
          E1(1,2)=-AD(1,2)/DT
          E1(2,1)=-AD(2,1)/DT
          E1(2,2)= AD(1,1)/DT
C
          CS(1,1)=E1(1,1)*AU(1,1)+E1(1,2)*AU(2,1)
          CS(1,2)=E1(1,1)*AU(1,2)+E1(1,2)*AU(2,2)
          CS(2,1)=E1(2,1)*AU(1,1)+E1(2,2)*AU(2,1)
          CS(2,2)=E1(2,1)*AU(1,2)+E1(2,2)*AU(2,2)
C
          E1(1,1)=BD(1,1)*CS(1,1)+BD(1,2)*CS(2,1)-BU(1,1)
          E1(1,2)=BD(1,1)*CS(1,2)+BD(1,2)*CS(2,2)-BU(1,2)
          E1(2,1)=BD(2,1)*CS(1,1)+BD(2,2)*CS(2,1)-BU(2,1)
          E1(2,2)=BD(2,1)*CS(1,2)+BD(2,2)*CS(2,2)-BU(2,2)
          DT=E1(1,1)*E1(2,2)-E1(1,2)*E1(2,1)
          E2(1,1)= E1(2,2)/DT
          E2(1,2)=-E1(1,2)/DT
          E2(2,1)=-E1(2,1)/DT
          E2(2,2)= E1(1,1)/DT
C
          CUSUX(1)=E2(1,1)*DQX
          CUSUX(2)=E2(2,1)*DQX
          CDSUX(1)=GRU(1,1,LS-1)*CUSUX(1)+GRU(1,2,LS-1)*CUSUX(2)
          CDSUX(2)=GRU(2,1,LS-1)*CUSUX(1)+GRU(2,2,LS-1)*CUSUX(2)
          CDSDX(1)=CS(1,1)*CUSUX(1)+CS(1,2)*CUSUX(2)
          CDSDX(2)=CS(2,1)*CUSUX(1)+CS(2,2)*CUSUX(2)
C
          CUSUZ(1)=E2(1,2)*DQZ
          CUSUZ(2)=E2(2,2)*DQZ
          CDSUZ(1)=GRU(1,1,LS-1)*CUSUZ(1)+GRU(1,2,LS-1)*CUSUZ(2)
          CDSUZ(2)=GRU(2,1,LS-1)*CUSUZ(1)+GRU(2,2,LS-1)*CUSUZ(2)
          CDSDZ(1)=CS(1,1)*CUSUZ(1)+CS(1,2)*CUSUZ(2)
          CDSDZ(2)=CS(2,1)*CUSUZ(1)+CS(2,2)*CUSUZ(2)
C
        ELSE
C
          CALL DDUPSV(OM,WN,SDEP,CMS,CMP,DD)
          CALL DDUPSV(OM,WN,THK(LS)-SDEP,CMS,CMP,DU)
C
          E1(1,1)=DD(1,1)*GRU(1,1,LS-1)+DD(1,2)*GRU(2,1,LS-1)
          E1(1,2)=DD(1,1)*GRU(1,2,LS-1)+DD(1,2)*GRU(2,2,LS-1)
          E1(2,1)=DD(2,1)*GRU(1,1,LS-1)+DD(2,2)*GRU(2,1,LS-1)
          E1(2,2)=DD(2,1)*GRU(1,2,LS-1)+DD(2,2)*GRU(2,2,LS-1)
C
          AU(1,1)=EE(1,1)*E1(1,1)+EE(1,2)*E1(2,1)+EE(1,3)
          AU(1,2)=EE(1,1)*E1(1,2)+EE(1,2)*E1(2,2)+EE(1,4)
          AU(2,1)=EE(2,1)*E1(1,1)+EE(2,2)*E1(2,1)+EE(2,3)
          AU(2,2)=EE(2,1)*E1(1,2)+EE(2,2)*E1(2,2)+EE(2,4)
          BU(1,1)=EE(3,1)*E1(1,1)+EE(3,2)*E1(2,1)+EE(3,3)
          BU(1,2)=EE(3,1)*E1(1,2)+EE(3,2)*E1(2,2)+EE(3,4)
          BU(2,1)=EE(4,1)*E1(1,1)+EE(4,2)*E1(2,1)+EE(4,3)
          BU(2,2)=EE(4,1)*E1(1,2)+EE(4,2)*E1(2,2)+EE(4,4)
C
          E1(1,1)=DU(1,1)*GRD(1,1,LS)+DU(1,2)*GRD(2,1,LS)
          E1(1,2)=DU(1,1)*GRD(1,2,LS)+DU(1,2)*GRD(2,2,LS)
          E1(2,1)=DU(2,1)*GRD(1,1,LS)+DU(2,2)*GRD(2,1,LS)
          E1(2,2)=DU(2,1)*GRD(1,2,LS)+DU(2,2)*GRD(2,2,LS)
C
          AD(1,1)=EE(1,3)*E1(1,1)+EE(1,4)*E1(2,1)+EE(1,1)
          AD(1,2)=EE(1,3)*E1(1,2)+EE(1,4)*E1(2,2)+EE(1,2)
          AD(2,1)=EE(2,3)*E1(1,1)+EE(2,4)*E1(2,1)+EE(2,1)
          AD(2,2)=EE(2,3)*E1(1,2)+EE(2,4)*E1(2,2)+EE(2,2)
          BD(1,1)=EE(3,3)*E1(1,1)+EE(3,4)*E1(2,1)+EE(3,1)
          BD(1,2)=EE(3,3)*E1(1,2)+EE(3,4)*E1(2,2)+EE(3,2)
          BD(2,1)=EE(4,3)*E1(1,1)+EE(4,4)*E1(2,1)+EE(4,1)
          BD(2,2)=EE(4,3)*E1(1,2)+EE(4,4)*E1(2,2)+EE(4,2)
C
          DT=AD(1,1)*AD(2,2)-AD(1,2)*AD(2,1)
          E1(1,1)= AD(2,2)/DT
          E1(1,2)=-AD(1,2)/DT
          E1(2,1)=-AD(2,1)/DT
          E1(2,2)= AD(1,1)/DT
          CS(1,1)=E1(1,1)*AU(1,1)+E1(1,2)*AU(2,1)
          CS(1,2)=E1(1,1)*AU(1,2)+E1(1,2)*AU(2,2)
          CS(2,1)=E1(2,1)*AU(1,1)+E1(2,2)*AU(2,1)
          CS(2,2)=E1(2,1)*AU(1,2)+E1(2,2)*AU(2,2)
C
          E1(1,1)=BD(1,1)*CS(1,1)+BD(1,2)*CS(2,1)-BU(1,1)
          E1(1,2)=BD(1,1)*CS(1,2)+BD(1,2)*CS(2,2)-BU(1,2)
          E1(2,1)=BD(2,1)*CS(1,1)+BD(2,2)*CS(2,1)-BU(2,1)
          E1(2,2)=BD(2,1)*CS(1,2)+BD(2,2)*CS(2,2)-BU(2,2)
          DT=E1(1,1)*E1(2,2)-E1(1,2)*E1(2,1)
          E2(1,1)= E1(2,2)/DT
          E2(1,2)=-E1(1,2)/DT
          E2(2,1)=-E1(2,1)/DT
          E2(2,2)= E1(1,1)/DT
C
          CUSUX(1)=E2(1,1)*DQX
          CUSUX(2)=E2(2,1)*DQX
          CDSUX(1)=GRU(1,1,LS-1)*CUSUX(1)+GRU(1,2,LS-1)*CUSUX(2)
          CDSUX(2)=GRU(2,1,LS-1)*CUSUX(1)+GRU(2,2,LS-1)*CUSUX(2)
          CDSDX(1)=CS(1,1)*CUSUX(1)+CS(1,2)*CUSUX(2)
          CDSDX(2)=CS(2,1)*CUSUX(1)+CS(2,2)*CUSUX(2)
          CUSDX(1)=GRD(1,1,LS)*CDSDX(1)+GRD(1,2,LS)*CDSDX(2)
          CUSDX(2)=GRD(2,1,LS)*CDSDX(1)+GRD(2,2,LS)*CDSDX(2)
C
          CUSUZ(1)=E2(1,2)*DQZ
          CUSUZ(2)=E2(2,2)*DQZ
          CDSUZ(1)=GRU(1,1,LS-1)*CUSUZ(1)+GRU(1,2,LS-1)*CUSUZ(2)
          CDSUZ(2)=GRU(2,1,LS-1)*CUSUZ(1)+GRU(2,2,LS-1)*CUSUZ(2)
          CDSDZ(1)=CS(1,1)*CUSUZ(1)+CS(1,2)*CUSUZ(2)
          CDSDZ(2)=CS(2,1)*CUSUZ(1)+CS(2,2)*CUSUZ(2)
          CUSDZ(1)=GRD(1,1,LS)*CDSDZ(1)+GRD(1,2,LS)*CDSDZ(2)
          CUSDZ(2)=GRD(2,1,LS)*CDSDZ(1)+GRD(2,2,LS)*CDSDZ(2)
C
        ENDIF
C
      END IF
C
      IF(LS.EQ.1)GOTO 10
      CUX(1,LS-1)=GTU(1,1,LS-1)*CUSUX(1)+GTU(1,2,LS-1)*CUSUX(2)
      CUX(2,LS-1)=GTU(2,1,LS-1)*CUSUX(1)+GTU(2,2,LS-1)*CUSUX(2)
      CUZ(1,LS-1)=GTU(1,1,LS-1)*CUSUZ(1)+GTU(1,2,LS-1)*CUSUZ(2)
      CUZ(2,LS-1)=GTU(2,1,LS-1)*CUSUZ(1)+GTU(2,2,LS-1)*CUSUZ(2)
      DO 100 I=LS-2,1,-1
        CUX(1,I)  =GTU(1,1,I)*CUX(1,I+1)+GTU(1,2,I)*CUX(2,I+1)
        CUX(2,I)  =GTU(2,1,I)*CUX(1,I+1)+GTU(2,2,I)*CUX(2,I+1)
        CDX(1,I+1)=GRU(1,1,I)*CUX(1,I+1)+GRU(1,2,I)*CUX(2,I+1)
        CDX(2,I+1)=GRU(2,1,I)*CUX(1,I+1)+GRU(2,2,I)*CUX(2,I+1)
        CUZ(1,I)  =GTU(1,1,I)*CUZ(1,I+1)+GTU(1,2,I)*CUZ(2,I+1)
        CUZ(2,I)  =GTU(2,1,I)*CUZ(1,I+1)+GTU(2,2,I)*CUZ(2,I+1)
        CDZ(1,I+1)=GRU(1,1,I)*CUZ(1,I+1)+GRU(1,2,I)*CUZ(2,I+1)
        CDZ(2,I+1)=GRU(2,1,I)*CUZ(1,I+1)+GRU(2,2,I)*CUZ(2,I+1)
  100 CONTINUE
      CDX(1,1)=GRU(1,1,0)*CUX(1,1)+GRU(1,2,0)*CUX(2,1)
      CDX(2,1)=GRU(2,1,0)*CUX(1,1)+GRU(2,2,0)*CUX(2,1)
      CDZ(1,1)=GRU(1,1,0)*CUZ(1,1)+GRU(1,2,0)*CUZ(2,1)
      CDZ(2,1)=GRU(2,1,0)*CUZ(1,1)+GRU(2,2,0)*CUZ(2,1)
C
   10 IF(LS.EQ.NL)GOTO 20
      CDX(1,LS+1)=GTD(1,1,LS)*CDSDX(1)+GTD(1,2,LS)*CDSDX(2)
      CDX(2,LS+1)=GTD(2,1,LS)*CDSDX(1)+GTD(2,2,LS)*CDSDX(2)
      CDZ(1,LS+1)=GTD(1,1,LS)*CDSDZ(1)+GTD(1,2,LS)*CDSDZ(2)
      CDZ(2,LS+1)=GTD(2,1,LS)*CDSDZ(1)+GTD(2,2,LS)*CDSDZ(2)
      DO 200 I=LS+1,NL-1
        CDX(1,I+1)=GTD(1,1,I)*CDX(1,I)+GTD(1,2,I)*CDX(2,I)
        CDX(2,I+1)=GTD(2,1,I)*CDX(1,I)+GTD(2,2,I)*CDX(2,I)
        CUX(1,I)  =GRD(1,1,I)*CDX(1,I)+GRD(1,2,I)*CDX(2,I)
        CUX(2,I)  =GRD(2,1,I)*CDX(1,I)+GRD(2,2,I)*CDX(2,I)
        CDZ(1,I+1)=GTD(1,1,I)*CDZ(1,I)+GTD(1,2,I)*CDZ(2,I)
        CDZ(2,I+1)=GTD(2,1,I)*CDZ(1,I)+GTD(2,2,I)*CDZ(2,I)
        CUZ(1,I)  =GRD(1,1,I)*CDZ(1,I)+GRD(1,2,I)*CDZ(2,I)
        CUZ(2,I)  =GRD(2,1,I)*CDZ(1,I)+GRD(2,2,I)*CDZ(2,I)
  200 CONTINUE
C
   20 RETURN
      END
C
C
C
C *** P-SV GENERALIZED R/T MATIX ***
C
      SUBROUTINE GNRTPSV(OM,WN)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VP,VS,GTD,GTU,GRD,GRU,TD,TU,RD,RU,E1(2,2),E2(2,2),DT
      COMMON /SOLY/SDEP,ODEP,LS,LO
     *    /STRC/DNS(20),VP(20),VS(20),THK(20),NL
      COMMON /MRTPSV/TD(2,2,0:20),TU(2,2,0:20),RD(2,2,0:20),RU(2,2,0:20)
     *  /GRTPSV/GTD(2,2,0:20),GTU(2,2,0:20),GRD(2,2,0:20),GRU(2,2,0:20)
C
C *** CHECK FOR THE DIMENSION OF GTD, GTU, GRD, AND GRU ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(*,*)' FROM GNRTPSV: NL.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C *** MODIFIED R/T MATRIX ***
C
      CALL MDRTPSV(OM,WN)
C
C *** GENERALIZED R/T MATRIX ***
C
      GRU(1,1,0)=RU(1,1,0)
      GRU(1,2,0)=RU(1,2,0)
      GRU(2,1,0)=RU(2,1,0)
      GRU(2,2,0)=RU(2,2,0)
C
      DO 100 I=1,LS-1
C
        E1(1,1)=1.D0-RD(1,1,I)*GRU(1,1,I-1)-RD(1,2,I)*GRU(2,1,I-1)
        E1(1,2)=    -RD(1,1,I)*GRU(1,2,I-1)-RD(1,2,I)*GRU(2,2,I-1)
        E1(2,1)=    -RD(2,1,I)*GRU(1,1,I-1)-RD(2,2,I)*GRU(2,1,I-1)
        E1(2,2)=1.D0-RD(2,1,I)*GRU(1,2,I-1)-RD(2,2,I)*GRU(2,2,I-1)
        DT=E1(1,1)*E1(2,2)-E1(1,2)*E1(2,1)
        E2(1,1)= E1(2,2)/DT
        E2(1,2)=-E1(1,2)/DT
        E2(2,1)=-E1(2,1)/DT
        E2(2,2)= E1(1,1)/DT
        GTU(1,1,I)=E2(1,1)*TU(1,1,I)+E2(1,2)*TU(2,1,I)
        GTU(1,2,I)=E2(1,1)*TU(1,2,I)+E2(1,2)*TU(2,2,I)
        GTU(2,1,I)=E2(2,1)*TU(1,1,I)+E2(2,2)*TU(2,1,I)
        GTU(2,2,I)=E2(2,1)*TU(1,2,I)+E2(2,2)*TU(2,2,I)
C
        E1(1,1)=TD(1,1,I)*GRU(1,1,I-1)+TD(1,2,I)*GRU(2,1,I-1)
        E1(1,2)=TD(1,1,I)*GRU(1,2,I-1)+TD(1,2,I)*GRU(2,2,I-1)
        E1(2,1)=TD(2,1,I)*GRU(1,1,I-1)+TD(2,2,I)*GRU(2,1,I-1)
        E1(2,2)=TD(2,1,I)*GRU(1,2,I-1)+TD(2,2,I)*GRU(2,2,I-1)
        E2(1,1)=E1(1,1)*GTU(1,1,I)+E1(1,2)*GTU(2,1,I)
        E2(1,2)=E1(1,1)*GTU(1,2,I)+E1(1,2)*GTU(2,2,I)
        E2(2,1)=E1(2,1)*GTU(1,1,I)+E1(2,2)*GTU(2,1,I)
        E2(2,2)=E1(2,1)*GTU(1,2,I)+E1(2,2)*GTU(2,2,I)
        GRU(1,1,I)=RU(1,1,I)+E2(1,1)
        GRU(1,2,I)=RU(1,2,I)+E2(1,2)
        GRU(2,1,I)=RU(2,1,I)+E2(2,1)
        GRU(2,2,I)=RU(2,2,I)+E2(2,2)
C
  100 CONTINUE
C
      GTD(1,1,NL-1)=TD(1,1,NL-1)
      GTD(1,2,NL-1)=TD(1,2,NL-1)
      GTD(2,1,NL-1)=TD(2,1,NL-1)
      GTD(2,2,NL-1)=TD(2,2,NL-1)
C
      GRD(1,1,NL-1)=RD(1,1,NL-1)
      GRD(1,2,NL-1)=RD(1,2,NL-1)
      GRD(2,1,NL-1)=RD(2,1,NL-1)
      GRD(2,2,NL-1)=RD(2,2,NL-1)
C
      DO 200 I=NL-2,LS,-1
C
        E1(1,1)=1.D0-RU(1,1,I)*GRD(1,1,I+1)-RU(1,2,I)*GRD(2,1,I+1)
        E1(1,2)=    -RU(1,1,I)*GRD(1,2,I+1)-RU(1,2,I)*GRD(2,2,I+1)
        E1(2,1)=    -RU(2,1,I)*GRD(1,1,I+1)-RU(2,2,I)*GRD(2,1,I+1)
        E1(2,2)=1.D0-RU(2,1,I)*GRD(1,2,I+1)-RU(2,2,I)*GRD(2,2,I+1)
        DT=E1(1,1)*E1(2,2)-E1(1,2)*E1(2,1)
        E2(1,1)= E1(2,2)/DT
        E2(1,2)=-E1(1,2)/DT
        E2(2,1)=-E1(2,1)/DT
        E2(2,2)= E1(1,1)/DT
        GTD(1,1,I)=E2(1,1)*TD(1,1,I)+E2(1,2)*TD(2,1,I)
        GTD(1,2,I)=E2(1,1)*TD(1,2,I)+E2(1,2)*TD(2,2,I)
        GTD(2,1,I)=E2(2,1)*TD(1,1,I)+E2(2,2)*TD(2,1,I)
        GTD(2,2,I)=E2(2,1)*TD(1,2,I)+E2(2,2)*TD(2,2,I)
C
        E1(1,1)=TU(1,1,I)*GRD(1,1,I+1)+TU(1,2,I)*GRD(2,1,I+1)
        E1(1,2)=TU(1,1,I)*GRD(1,2,I+1)+TU(1,2,I)*GRD(2,2,I+1)
        E1(2,1)=TU(2,1,I)*GRD(1,1,I+1)+TU(2,2,I)*GRD(2,1,I+1)
        E1(2,2)=TU(2,1,I)*GRD(1,2,I+1)+TU(2,2,I)*GRD(2,2,I+1)
        E2(1,1)=E1(1,1)*GTD(1,1,I)+E1(1,2)*GTD(2,1,I)
        E2(1,2)=E1(1,1)*GTD(1,2,I)+E1(1,2)*GTD(2,2,I)
        E2(2,1)=E1(2,1)*GTD(1,1,I)+E1(2,2)*GTD(2,1,I)
        E2(2,2)=E1(2,1)*GTD(1,2,I)+E1(2,2)*GTD(2,2,I)
        GRD(1,1,I)=RD(1,1,I)+E2(1,1)
        GRD(1,2,I)=RD(1,2,I)+E2(1,2)
        GRD(2,1,I)=RD(2,1,I)+E2(2,1)
        GRD(2,2,I)=RD(2,2,I)+E2(2,2)
C
  200 CONTINUE
C
      RETURN
      END
C
C *** P-SV MODIFIED R/T MATIX ***
C
      SUBROUTINE MDRTPSV(OM,WN)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 CMS,CMP,EA(4,4),EB(4,4),TD,TU,RD,RU
     *       ,DDA(2,2),DUB(2,2),DET,VP,VS
      COMMON /SOLY/SDEP,ODEP,LS,LO
     *    /STRC/DNS(20),VP(20),VS(20),THK(20),NL
     *    /MRTPSV/TD(2,2,0:20),TU(2,2,0:20),RD(2,2,0:20),RU(2,2,0:20)
C
C *** CHECK FOR THE DIMENSION OF COMMON ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(*,*)' FROM MDRTPSV: NL.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C *** MODIFIED R/T MATRIX ***
C
      CALL MTXPSV(OM,WN,1,CMS,CMP,EA)
C
      IF(LS.EQ.1)THEN
        CALL DDUPSV(OM,WN,SDEP,CMS,CMP,DUB)
      ELSE
        CALL DDUPSV(OM,WN,THK(1),CMS,CMP,DUB)
      END IF
C
      DET=EA(3,1)*EA(4,2)-EA(3,2)*EA(4,1)
      RU(1,1,0)=-(DUB(1,1)*(EA(3,3)*EA(4,2)-EA(3,2)*EA(4,3))
     *           +DUB(2,1)*(EA(3,4)*EA(4,2)-EA(3,2)*EA(4,4)))/DET
      RU(1,2,0)=-(DUB(1,2)*(EA(3,3)*EA(4,2)-EA(3,2)*EA(4,3))
     *           +DUB(2,2)*(EA(3,4)*EA(4,2)-EA(3,2)*EA(4,4)))/DET
      RU(2,1,0)=-(DUB(1,1)*(EA(3,1)*EA(4,3)-EA(3,3)*EA(4,1))
     *           +DUB(2,1)*(EA(3,1)*EA(4,4)-EA(3,4)*EA(4,1)))/DET
      RU(2,2,0)=-(DUB(1,2)*(EA(3,1)*EA(4,3)-EA(3,3)*EA(4,1))
     *           +DUB(2,2)*(EA(3,1)*EA(4,4)-EA(3,4)*EA(4,1)))/DET
C
      IF(NL.EQ.1)RETURN
C
      IF(LS.EQ.1)THEN
        CALL DDUPSV(OM,WN,THK(1)-SDEP,CMS,CMP,DDA)
      ELSE
        CALL DDUPSV(OM,WN,THK(1),CMS,CMP,DDA)
      END IF
C
      DO 100 I=1, NL-1
C
        CALL MTXPSV(OM,WN,I+1,CMS,CMP,EB)
C
        IF(LS.EQ.(I+1))THEN
          CALL DDUPSV(OM,WN,SDEP,CMS,CMP,DUB)
        ELSE
          CALL DDUPSV(OM,WN,THK(I+1),CMS,CMP,DUB)
        END IF
C
        CALL MTRPSV(EA,EB,DDA,DUB,I) 
C
        DO 110 J=1,4
          DO 120 K=1,4
            EA(J,K)=EB(J,K)
  120     CONTINUE
  110   CONTINUE
C
        IF(LS.EQ.(I+1))THEN
          IF(LS.NE.NL)THEN
            CALL DDUPSV(OM,WN,THK(I+1)-SDEP,CMS,CMP,DDA)
          END IF
        ELSE
          CALL DDUPSV(OM,WN,THK(I+1),CMS,CMP,DDA)
        END IF
C
  100 CONTINUE
C
      RETURN
      END
C
C
C *** P-SV MODIFIED R/T MATRIX FOR Ith LAYER ***
C
        SUBROUTINE MTRPSV(EA,EB,DDA,DUB,L) 
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 EA(4,4),EB(4,4),DDA(2,2),DUB(2,2),TD,TU,RD,RU
     *          ,E1(4,4),E2(4,4),E3(4,4),VP,VS
     * ,a1(4),a2(4),a3(4),b1(4),b2(4),b3(4)  
     * ,c1(4),c2(4),c3(4),d1(4),d2(4),d3(4)  
     * ,aa(4),EE(4,4)
      COMMON /STRC/DNS(20),VP(20),VS(20),THK(20),NL
     *      /MRTPSV/TD(2,2,0:20),TU(2,2,0:20),RD(2,2,0:20),RU(2,2,0:20)
C
C *** CHECK FOR THE DIMENSION OF COMMON ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(*,*)' FROM MTRPSV: NL.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C *** MATRIX ***
C
      DO 100 I=1,2
        DO 110 J=1,2
          EE(I,J)= EB(I,J)
          EE(I,J+2)=-EA(I,J+2)
          EE(I+2,J)= EB(I+2,J)
          EE(I+2,J+2)=-EA(I+2,J+2)
          E2(I,J)= EA(I,J)
          E2(I,J+2)=-EB(I,J+2)
          E2(I+2,J)= EA(I+2,J)
          E2(I+2,J+2)=-EB(I+2,J+2)
  110   CONTINUE
  100 CONTINUE
C
C     * INVERSE MATRIX OF EE *
C
        DO 222 I=1,4
C
        IF(I.EQ.1)THEN
          J=2
          K=3
          M=4
        ELSE
          IF(I.EQ.2)THEN
            J=3
            K=4
            M=1
          ELSE
            IF(I.EQ.3)THEN
              J=4
              K=1
              M=2
            ELSE
              J=1
              K=2
              M=3
            END IF
          END IF
        END IF
C
        A1(I)=EE(I,4)*EE(J,1)-EE(I,1)*EE(J,4)  
        A2(I)=EE(I,4)*EE(J,2)-EE(I,2)*EE(J,4)  
        A3(I)=EE(I,4)*EE(J,3)-EE(I,3)*EE(J,4)  
        B1(I)=EE(I,3)*EE(J,2)-EE(I,2)*EE(J,3)  
        B2(I)=EE(I,1)*EE(J,3)-EE(I,3)*EE(J,1)  
        B3(I)=EE(I,2)*EE(J,1)-EE(I,1)*EE(J,2)  
        C1(I)=EE(K,3)*EE(M,2)-EE(K,2)*EE(M,3)  
        C2(I)=EE(K,1)*EE(M,3)-EE(K,3)*EE(M,1)  
        C3(I)=EE(K,2)*EE(M,1)-EE(K,1)*EE(M,2)  
        D1(I)=EE(K,4)*EE(M,1)-EE(K,1)*EE(M,4)  
        D2(I)=EE(K,4)*EE(M,2)-EE(K,2)*EE(M,4)  
        D3(I)=EE(K,4)*EE(M,3)-EE(K,3)*EE(M,4)  
C
        AA(I)=A1(I)*C1(I)+A2(I)*C2(I)+A3(I)*C3(I)
     *       +B1(I)*D1(I)+B2(I)*D2(I)+B3(I)*D3(I)
C
        E1(1,I)=(-EE(J,2)*D3(I)+EE(J,3)*D2(I)-EE(J,4)*C1(I))/AA(I)
        E1(2,I)=( EE(J,1)*D3(I)-EE(J,3)*D1(I)-EE(J,4)*C2(I))/AA(I)
        E1(3,I)=(-EE(J,1)*D2(I)+EE(J,2)*D1(I)-EE(J,4)*C3(I))/AA(I)
        E1(4,I)=( EE(J,1)*C1(I)+EE(J,2)*C2(I)+EE(J,3)*C3(I))/AA(I)
C
  222   CONTINUE
C
      DO 200 I=1,4
        DO 210 J=1,4
          E3(I,J)=(0.D0,0.D0)
          DO 220 K=1,4
            E3(I,J)=E3(I,J)+E1(I,K)*E2(K,J) 
  220     CONTINUE
  210   CONTINUE
  200 CONTINUE
C
      TD(1,1,L)=E3(1,1)*DDA(1,1)+E3(1,2)*DDA(2,1)  
      TD(2,1,L)=E3(2,1)*DDA(1,1)+E3(2,2)*DDA(2,1)  
      TD(1,2,L)=E3(1,1)*DDA(1,2)+E3(1,2)*DDA(2,2)  
      TD(2,2,L)=E3(2,1)*DDA(1,2)+E3(2,2)*DDA(2,2)  
      RD(1,1,L)=E3(3,1)*DDA(1,1)+E3(3,2)*DDA(2,1)  
      RD(2,1,L)=E3(4,1)*DDA(1,1)+E3(4,2)*DDA(2,1)  
      RD(1,2,L)=E3(3,1)*DDA(1,2)+E3(3,2)*DDA(2,2)  
      RD(2,2,L)=E3(4,1)*DDA(1,2)+E3(4,2)*DDA(2,2)  
C
      IF(L.EQ.NL)RETURN
C
      RU(1,1,L)=E3(1,3)*DUB(1,1)+E3(1,4)*DUB(2,1)  
      RU(2,1,L)=E3(2,3)*DUB(1,1)+E3(2,4)*DUB(2,1)  
      RU(1,2,L)=E3(1,3)*DUB(1,2)+E3(1,4)*DUB(2,2)  
      RU(2,2,L)=E3(2,3)*DUB(1,2)+E3(2,4)*DUB(2,2)  
      TU(1,1,L)=E3(3,3)*DUB(1,1)+E3(3,4)*DUB(2,1)  
      TU(2,1,L)=E3(4,3)*DUB(1,1)+E3(4,4)*DUB(2,1)  
      TU(1,2,L)=E3(3,3)*DUB(1,2)+E3(3,4)*DUB(2,2)  
      TU(2,2,L)=E3(4,3)*DUB(1,2)+E3(4,4)*DUB(2,2)  
C
      RETURN
      END
C
C
C *** P-SV ELEMENT MATIX ***
C
      SUBROUTINE MTXPSV(OM,WN,L,CMS,CMP,EE)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VP,VS,VS2,WS2,CMS,VP2,WP2,CMP
     *      ,SM,EE(4,4),X0,AA,CMSG,CMPG
      COMMON /STRC/DNS(20),VP(20),VS(20),THK(20),NL
C
C *** CHECK FOR THE DIMENSION OF COMMOM ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(*,*)' FROM MTXPSV: NL.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C *** ELEMENT MATRIX ***
C
      VS2=VS(L)*VS(L)
      VP2=VP(L)*VP(L)
      SM=DNS(L)*VS2
C
      IF(OM.EQ.0.D0)THEN
C  
        X0=(1.D0+VS2/VP2)/(1.D0-VS2/VP2)
        AA=WN*SM
C
        EE(1,1)=(1.D0,0.D0)
        EE(1,2)=(1.D0,0.D0)
        EE(1,3)=(1.D0,0.D0)
        EE(1,4)=(1.D0,0.D0)
        EE(2,1)=-(X0-1.D0)
        EE(2,2)=(1.D0,0.D0)
        EE(2,3)= X0-1.D0
        EE(2,4)=(-1.D0,0.D0)
        EE(3,1)= (X0-3.D0)*AA
        EE(3,2)= -2.D0*AA
        EE(3,3)=-(X0-3.D0)*AA
        EE(3,4)= 2.D0*AA
        EE(4,1)= (X0-1.D0)*AA
        EE(4,2)= -2.D0*AA
        EE(4,3)= (X0-1.D0)*AA
        EE(4,4)=-2.D0*AA
C
      ELSE
C
        WN2=WN*WN
        OM2=OM*OM
C
        WS2=OM2/VS2
        CMS=CDSQRT(WN2-WS2)
        IF(DREAL(CMS).LT.0.D0)CMS=-CMS
        CMSG=CMS/WN
C
        WP2=OM2/VP2
        CMP=CDSQRT(WN2-WP2)
        IF(DREAL(CMP).LT.0.D0)CMP=-CMP
        CMPG=CMP/WN
C
        X0=2.D0*WN-WS2/WN
C
        EE(1,1)=(-1.D0,0.D0)
        EE(1,2)= CMSG
        EE(1,3)=(-1.D0,0.D0)
        EE(1,4)= CMSG
        EE(2,1)=-CMPG
        EE(2,2)=(1.D0,0.D0)
        EE(2,3)= CMPG
        EE(2,4)=(-1.D0,0.D0)
        EE(3,1)= 2.D0*SM*CMP 
        EE(3,2)=-SM*X0 
        EE(3,3)=-2.D0*SM*CMP 
        EE(3,4)= SM*X0 
        EE(4,1)= SM*X0 
        EE(4,2)=-2.D0*SM*CMS 
        EE(4,3)= SM*X0 
        EE(4,4)=-2.D0*SM*CMS 
C
      END IF
C
      RETURN
      END
C
C
C *** P-SV UP/DOWN-GOING MATRIX ***
C
      SUBROUTINE DDUPSV(OM,WN,THK,CMS,CMP,DU)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 DU(2,2),CMS,CMP
C
      IF(OM.EQ.0.D0)THEN
        EW=DEXP(-WN*THK)
        DU(1,1)= EW
        DU(1,2)= (0.D0,0.D0)
        DU(2,1)= -WN*THK*EW
        DU(2,2)= EW
      ELSE
        DU(1,1)= CDEXP(-CMP*THK)
        DU(1,2)= (0.D0,0.D0)
        DU(2,1)= (0.D0,0.D0)
        DU(2,2)= CDEXP(-CMS*THK)
      END IF
C
      RETURN
      END
C
C
C   NAME                - BESJ0 : DOUBLE PRECISION
C
      FUNCTION BESJ0 (X)
C
      IMPLICIT REAL*8 (A-H,O-Z)
      DATA C0,C1,C2,C3,C4,C5,C6,C7,C8,C9,C10
     *   / -2.2270 614            D-9,  7.8594 0046           D-8,
     *     -1.9891 96578 6        D-6,  4.0310 96038 59       D-5,
     *     -6.4499 78225 950      D-4,  7.9012 34182 8039     D-3,
     *     -7.1111 11100 46304    D-2,  4.4444 44444 25769 5  D-1,
     *     -1.7777 77777 77584 46 D 0,  3.9999 99999 99989 66 D 0,
     *     -3.9999 99999 99999 78 D 0                             /
      DATA D0,D1,D2,D3,D4,D5,D6,D7,D8,D9,D10,D11,D12,D13
     *   /  3.3150 18             D-11,  4.6970 63             D-11,
     *     -5.6633 6210           D- 9, -4.8842 9867           D- 9,
     *      6.7810 96622 5        D- 7,  2.1332 90985 8        D- 7,
     *     -5.3795 00918 201      D- 5,  1.1202 21502 384      D- 5,
     *      2.5202 11970 18629    D- 3, -1.7073 87596 90737    D- 3,
     *     -5.6148 69347 44975 4  D- 2,  5.6040 47189 80261 6  D- 2,
     *      3.2757 91375 91465 09 D- 1, -1.7759 67713 14338 34 D- 1 /
      DATA E0,E1,E2,E3,E4,E5,E6,E7,E8,E9,E10,E11,E12,E13,E14
     *   / -1.9865 0              D-12, -1.5178 02          D-11,
     *      3.9649 790            D-10,  2.4051 2241        D- 9,
     *     -5.7071 43476          D- 8, -2.5797 89310 5     D- 7,
     *      5.6607 46146 93       D- 6,  1.7226 12595 802   D- 5,
     *     -3.5284 91002 3231     D- 4, -5.9314 89738 8202  D- 4,
     *      1.1790 12935 71024 6  D- 2,  6.3961 29897 84324 D- 3,
     *     -1.5037 41226 51373 90 D- 1,  4.6828 23482 34590 D- 3,
     *      3.0007 92705 19555 60 D- 1                            /
      DATA A0,A1,A2,A3,A4,A5,A6,A7,A8,A9,A10
     *   /  2.8228 605             D-10, -1.9415 5134       D-9,
     *      6.5160 2613            D- 9, -1.5339 82923      D-8,
     *      3.2676 79187           D- 8, -8.0992 84883      D-8,
     *      2.8872 15007 2         D- 7, -1.7424 98968 02   D-6,
     *      2.1846 78334 973       D- 5, -8.7658 21590 6110 D-4,
     *      0.7978 84560 80286 526 D  0                         /
      DATA B0,B1,B2,B3,B4,B5,B6,B7,B8,B9,B10,B11,B12
     *   / -8.9331 052           D-10,  6.6532 7442       D-9,
     *     -2.3224 79945         D- 8,  5.1878 39807      D-8,
     *     -8.7651 63373         D- 8,  1.3089 75220 0    D-7,
     *     -2.0632 29204 9       D- 7,  4.0856 34631 8    D-7,
     *     -1.1593 30628 42      D- 6,  5.2586 51475 08   D-6,
     *     -4.4239 73781 374     D- 5,  9.1310 64157 0497 D-4,
     *     -9.9735 57010 03581 5 D- 2                         /
      DATA PI4  / 0.78539 81633 97448 31D0 /
      AX=ABS(X)
      IF (AX .GE. 1.0D12) GO TO 40
      IF (AX .GE. 8.0D0) GO TO 30
      IF (AX .GE. 6.0D0) GO TO 20
      IF (AX .GE. 4.0D0) GO TO 10
C
C     !X! < 4
C
      Z=(X*0.25D0)**2
      BESJ0=((((((((((C0*Z +C1)*Z +C2)*Z +C3)*Z +C4)*Z +C5)*Z +C6)*Z
     *       +C7)*Z +C8)*Z +C9)*Z +C10)*Z +1.0D0
      GO TO 9999
C
C     4 <= !X! < 6
C
   10 Z=AX-5.0D0
      BESJ0=((((((((((((D0*Z +D1)*Z +D2)*Z +D3)*Z +D4)*Z +D5)*Z
     *       +D6)*Z +D7)*Z +D8)*Z +D9)*Z +D10)*Z +D11)*Z +D12)*Z +D13
      GO TO 9999
C
C     6 <= !X! < 8
C
   20 Z=AX-7.0D0
      BESJ0=(((((((((((((E0*Z +E1)*Z +E2)*Z +E3)*Z +E4)*Z +E5)*Z
     *       +E6)*Z +E7)*Z +E8)*Z +E9)*Z +E10)*Z +E11)*Z +E12)*Z
     *       +E13)*Z +E14
      GO TO 9999
C
C     8 <= !X! < 1.0D12
C
   30 Z=(8.0D0/AX)**2
      A=(((((((((A0*Z +A1)*Z +A2)*Z +A3)*Z +A4)*Z +A5)*Z +A6)*Z
     *  +A7)*Z +A8)*Z +A9)*Z +A10
      B=((((((((((((B0*Z +B1)*Z +B2)*Z +B3)*Z +B4)*Z +B5)*Z +B6)*Z
     *  +B7)*Z +B8)*Z +B9)*Z +B10)*Z +B11)*Z +B12) /AX
      Q=AX-PI4
      BESJ0=(A*COS(Q) - B*SIN(Q)) /SQRT(AX)
      GO TO 9999
C
C     PARAMETER ERROR
C
   40 BESJ0=0.0D0
C
 9999 RETURN
      END
C
C
C
C   NAME                - BESJ1 : DOUBLE PRECISION
C
      FUNCTION BESJ1 (X)
C
C
      IMPLICIT REAL*8 (A-H,O-Z)
      DATA C0,C1,C2,C3,C4,C5,C6,C7,C8,C9
     *   /  1.2249 4472           D-8, -3.9299 65711          D-7,
     *      8.9514 96872 0        D-6, -1.6124 40490 740      D-4,
     *      2.2574 92590 3842     D-3, -2.3703 70267 63209    D-2,
     *      1.7777 77775 58267 6  D-1, -8.8888 88888 61512 4  D-1,
     *      2.6666 66666 66491 50 D 0, -3.9999 99999 99995 63 D 0 /
      DATA D0,D1,D2,D3,D4,D5,D6,D7,D8,D9,D10,D11,D12,D13,D14
     *   /  2.2193 1              D-12,  4.1163 4              D-12,
     *     -4.3816 515            D-10, -5.7599 652            D-10,
     *      6.2306 13774          D- 8,  4.8857 13666          D- 8,
     *     -6.1029 92681 89       D- 6, -1.7066 40506 78       D- 6,
     *      3.7656 50660 9458     D- 4, -6.7213 28811 703      D- 5,
     *     -1.2601 05985 12045 4  D- 2,  6.8295 50387 40438    D- 3,
     *      1.6844 60804 23507 80 D- 1, -1.1208 09437 96045 28 D- 1,
     *     -3.2757 91375 91465 22 D- 1                              /
      DATA E0,E1,E2,E3,E4,E5,E6,E7,E8,E9,E10,E11,E12,E13,E14
     *   / -1.0746 3             D-12,  2.7780 03             D-11,
     *      2.0080 677           D-10, -4.7578 8198           D- 9,
     *     -2.6460 77937         D- 8,  5.7071 42412 9        D- 7,
     *      2.3218 13149 97      D- 6, -4.5285 96911 743      D- 5,
     *     -1.2058 28825 8765    D- 4,  2.1170 94601 37864    D- 3,
     *      2.9657 44869 54234   D- 3, -4.7160 51742 84081 7  D- 2,
     *     -1.9188 38969 35370 8 D- 2,  3.0074 82453 02747 75 D- 1,
     *     -4.6828 23482 34583   D- 3                              /
      DATA A0,A1,A2,A3,A4,A5,A6,A7,A8,A9,A10
     *   / -3.0207 732            D-10,  2.0806 2912        D-9,
     *     -7.0003 5699           D- 9,  1.6557 68548       D-8,
     *     -3.5582 32648          D- 8,  8.9550 40706       D-8,
     *     -3.2722 78326 5        D- 7,  2.0593 18935 93    D-6,
     *     -2.8088 72167 333      D- 5,  1.4609 70265 11633 D-3,
     *      7.9788 45608 02865 46 D- 1                          /
      DATA B0,B1,B2,B3,B4,B5,B6,B7,B8,B9,B10,B11,B12
     *   /  9.4650 569            D-10, -7.0536 9184        D-9,
     *      2.4646 12368          D- 8, -5.5142 32274       D-8,
     *      9.3432 85420          D- 8, -1.4024 48404 3     D-7,
     *      2.2297 66956 2        D- 7, -4.4752 03358 4     D-7,
     *      1.2957 31867 99       D- 6, -6.0676 75972 38    D-6,
     *      5.4070 79074 732      D- 5, -1.2783 48981 98997 D-3,
     *      2.9920 67103 01074 48 D- 1                          /
      DATA PI4  / 0.78539 81633 97448 31D0 /
      AX=ABS(X)
      IF (AX .GE. 1.0D12) GO TO 50
      IF (AX .GE. 8.0D0) GO TO 30
      IF (AX .GE. 6.0D0) GO TO 20
      IF (AX .GE. 4.0D0) GO TO 10
C
C     !X! < 4
C
      Y=X*0.25D0
      Z=Y*Y
      BESJ1=((((((((((C0*Z +C1)*Z +C2)*Z +C3)*Z +C4)*Z +C5)*Z +C6)*Z
     *       +C7)*Z +C8)*Z +C9)*Z +2.0D0) *Y
      GO TO 9999
C
C     4 <= !X! < 6
C
   10 Z=AX-5.0D0
      BESJ1=(((((((((((((D0*Z +D1)*Z +D2)*Z +D3)*Z +D4)*Z +D5)*Z
     *       +D6)*Z +D7)*Z +D8)*Z +D9)*Z +D10)*Z +D11)*Z +D12)*Z
     *       +D13)*Z +D14
      GO TO 40
C
C     6 <= !X! < 8
C
   20 Z=AX-7.0D0
      BESJ1=(((((((((((((E0*Z +E1)*Z +E2)*Z +E3)*Z +E4)*Z +E5)*Z
     *       +E6)*Z +E7)*Z +E8)*Z +E9)*Z +E10)*Z +E11)*Z +E12)*Z
     *       +E13)*Z +E14
      GO TO 40
C
C     8 <= !X! < 1.0D12
C
   30 Z=(8.0D0/AX)**2
      A=(((((((((A0*Z +A1)*Z +A2)*Z +A3)*Z +A4)*Z +A5)*Z +A6)*Z +A7)*Z
     *  +A8)*Z +A9)*Z +A10
      B=((((((((((((B0*Z +B1)*Z +B2)*Z +B3)*Z +B4)*Z +B5)*Z +B6)*Z
     *  +B7)*Z +B8)*Z +B9)*Z +B10)*Z +B11)*Z +B12) /AX
      Q=AX-PI4
      BESJ1=(A*SIN(Q) + B*COS(Q)) /SQRT(AX)
   40 IF (X .LT. 0.0D0) BESJ1=-BESJ1
      GO TO 9999
C
C     PARAMETER ERROR
C
   50 BESJ1=0.0D0
C
 9999 RETURN
      END
