c
C *** THE COMPUTATION OF Q, V, AND W AND THEIR DERIVATIVES 
C     USING SIMPSON'S INTEGATION FOMULA (before k*r < 8)***
C
      SUBROUTINE QVWRZS0(OM,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),OM
     *          ,QXB(0:2),QZB(0:2),VXB(0:2),WXB(0:2),WZB(0:2),SM
     *          ,QXC(0:2),QZC(0:2),VXC(0:2),WXC(0:2),WZC(0:2),RAM
      COMMON /INTG/VPMAX,VSMAX,RYMIN,CFINAL,N11,N12,M11
     *       /WVNUM/WVN(2001),DWN,WNR(3),NWN
     *       /FCQVW/SM,RAM,RR  /IRANGE/NRNG,IRNG
	*       /sInteg/sImax1,sImax2
C
c *** Dimension Check ***
c
      NC=2001
      IF(NWN.GT.NC)THEN
         WRITE(10,*)' FROM QVWRZS: NWN.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C   * RANGE 1 : From near Zero to near om/Rayleigh(min) *
C
      sImax1=0.d0
      sImax2=0.d0
      NRNG=0
      IRNG=0
C
      WA=WNR(1)
      K=1
      CALL INTEGS0(OM,K,WA,QXA,QZA,VXA,WXA,WZA)
      DO  100 I=1,M11
        K=2*I
        WB=WVN(K)
        CALL INTEGS0(OM,K,WB,QXB,QZB,VXB,WXB,WZB)
        K=K+1
        WC=WVN(K)
        CALL INTEGS0(OM,K,WC,QXC,QZC,VXC,WXC,WZC)
c
        H=WC-WB
        DO 110 J=0,2
         GQX(J)= GQX(J)+H*(QXA(J)+4.0*QXB(J)+QXC(J))/3.0
         GVX(J)= GVX(J)+H*(VXA(J)+4.0*VXB(J)+VXC(J))/3.0
         GWX(J)= GWX(J)+H*(WXA(J)+4.0*WXB(J)+WXC(J))/3.0
         GQZ(J)= GQZ(J)+H*(QZA(J)+4.0*QZB(J)+QZC(J))/3.0
         GWZ(J)= GWZ(J)+H*(WZA(J)+4.0*WZB(J)+WZC(J))/3.0
         QXA(J)= QXC(J)
         VXA(J)= VXC(J)
         WXA(J)= WXC(J)
         QZA(J)= QZC(J)
         WZA(J)= WZC(J)
  110   CONTINUE
c
c     * Check	of Convergence of Integrands *
c
	  sI11=cdabs(QZA(2)+4.0*QZB(2)+QZB(2))/3.0
	  sI12=cdabs(WZA(2)+4.0*WZB(2)+WZB(2))/3.0
	  if(sI11.ge.sImax1)sImax1=sI11
	  if(sI12.ge.sImax2)sImax2=sI12
c	write(2,999)k,si11,simax1,si11/simax1,si12,simax2,si12/simax2
c 999  format(i5,6e12.4)
	  if(((sI11/sImax1).lt.0.0001d0).and.
     *                       ((sI12/sImax2).lt.0.0001d0))then
c	    write(*,*)' Simpson Range 1: k*R & Int.Num.=',WC*RR,K
c	    write(*,*)'sI11,sI12',sI11,sI12
c	    write(*,*)'sImax1,sImax2',sImax1,sImax2
c	    write(*,*)'sI11/sImax1 < 0.0001d !!: sI11/sImax1=',sI11/sImax1
c	    write(*,*)'sI12/sImax2 < 0.0001d !!: sI12/sImax2=',sI12/sImax2
	    goto 10
	  end if
c
        IF(RR*WC.GT.8.D0)THEN
          NRNG=1
          IRNG=I
          IF(I.EQ.M11)THEN
            NRNG=2
            IRNG=0
          END IF
          RETURN
        END IF
c
  100 CONTINUE
c
C   * RANGE 2 : From near om/Rayleigh(min) to near om/C(final) *
C
      H=DWN
      WA=WNR(2)
      DO  200 I=1,N12
        K=2*(M11+I)
        WB=WVN(K)
        CALL INTEGS0(OM,K,WB,QXB,QZB,VXB,WXB,WZB)
        K=K+1
        WC=WVN(K)
        CALL INTEGS0(OM,K,WC,QXC,QZC,VXC,WXC,WZC)
c
        DO 210 J=0,2
         GQX(J)= GQX(J)+H*(QXA(J)+4.0*QXB(J)+QXC(J))/3.0
         GVX(J)= GVX(J)+H*(VXA(J)+4.0*VXB(J)+VXC(J))/3.0
         GWX(J)= GWX(J)+H*(WXA(J)+4.0*WXB(J)+WXC(J))/3.0
         GQZ(J)= GQZ(J)+H*(QZA(J)+4.0*QZB(J)+QZC(J))/3.0
         GWZ(J)= GWZ(J)+H*(WZA(J)+4.0*WZB(J)+WZC(J))/3.0
         QXA(J)= QXC(J)
         VXA(J)= VXC(J)
         WXA(J)= WXC(J)
         QZA(J)= QZC(J)
         WZA(J)= WZC(J)
  210   CONTINUE
c
c     * Check	of Convergence of Integrands *
c
	  sI11=cdabs((QZA(2)+4.0*QZB(2)+QZB(2)))/3.0
	  sI12=cdabs((WZA(2)+4.0*WZB(2)+WZB(2)))/3.0
	  if(sI11.ge.sImax1)sImax1=sI11
	  if(sI12.ge.sImax2)sImax2=sI12
c	write(2,999)k,si11,simax1,si11/simax1,si12,simax2,si12/simax2
	  if(((sI11/sImax1).lt.0.0001d0).and.
     *                       ((sI12/sImax2).lt.0.0001d0))then
c	    write(*,*)' Simpson Range 2: k*R & Int.Num.=',WC*RR,K
c	    write(*,*)'sI11,sI12',sI11,sI12
c	    write(*,*)'sImax1,sImax2',sImax1,sImax2
c	    write(*,*)'sI11/sImax1 < 0.0001d !!: sI11/sImax1=',sI11/sImax1
c	    write(*,*)'sI12/sImax2 < 0.0001d !!: sI12/sImax2=',sI12/sImax2
	    goto 10
	  end if
c
        IF(RR*WC.GT.8.D0)THEN
          NRNG=2
          IRNG=I
          IF(I.EQ.N12)THEN
            NRNG=3
            IRNG=0
          END IF
          RETURN
        END IF

  200 CONTINUE
C
   10 NRNG=3
c       
      RETURN
      END
C
C
C *** THE INTEGRANDS OF Q, V, AND W AND THEIR DERIVATIVES for Simpson's ***
C
      SUBROUTINE INTEGS0(OM,K,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),VLX0(2),URX0(2,2),URZ0(2,2)
     *   ,VLXZ,URXZ(2),URZZ(2),SM,RAM,OM
      DIMENSION BJ0(0:1),BJ1(0:1),BJK(0:1),BJD(0:1)
      COMMON /FREQ/DOM,IOM,NOM,OMI /SOLY/SDEP,ODEP,LS,LO 
     *       /STRC/DNS(20),VP(20),QP(20),FP(20)
     *                    ,VS(20),QS(20),FS(20),THK(20),NL
     *       /OTST/IOT,IST /FCQVW/SM,RAM,RR
C
C   * CHECK OF DIMENSION *
C
      NC=20
      IF(NL.GE.NC)THEN
        WRITE(10,*)' FROM INTEGS:NL>NC !!!!!!!!!!'
        WRITE(10,*)' NL,NC=',NL,NC
        STOP
      END IF
C
c   * the displacement-stress vectors *
c
      if(IST.eq.0)then
        CALL RTSH(OM,K,WN,VLX)
        CALL RTPSV(OM,K,WN,URX,URZ)
        CALL RTSH0(dcmplx(0.d0,OMI),K,WN,VLX0)
        CALL RTPSV0(dcmplx(0.d0,OMI),K,WN,URX0,URZ0)
        VLX(1)=VLX(1)-VLX0(1)
        URX(1,1)=URX(1,1)-URX0(1,1)
        URX(2,1)=URX(2,1)-URX0(2,1)
        URZ(1,1)=URZ(1,1)-URZ0(1,1)
        URZ(2,1)=URZ(2,1)-URZ0(2,1)
        VLX(2)=VLX(2)-VLX0(2)
        URX(1,2)=URX(1,2)-URX0(1,2)
        URX(2,2)=URX(2,2)-URX0(2,2)
        URZ(1,2)=URZ(1,2)-URZ0(1,2)
        URZ(2,2)=URZ(2,2)-URZ0(2,2)
	else
        if(IST.eq.1)then
          CALL RTSH(OM,K,WN,VLX)
          CALL RTPSV(OM,K,WN,URX,URZ)
	  else
          if(IST.eq.2)then
            CALL RTSH0(dcmplx(0.d0,OMI),K,WN,VLX)
            CALL RTPSV0(dcmplx(0.d0,OMI),K,WN,URX,URZ)
c            CALL RTSHc(dcmplx(0.d0,OMI),dcmplx(WN,0.d0),VLX)
c            CALL RTPSVc(dcmplx(0.d0,OMI),dcmplx(WN,0.d0),URX,URZ)
	    end if
	  end if
      end if
C
c	write(*,*)'k,wn',k,wn
c	write(*,*)'vlx(1)',vlx(1)
c	write(*,*)'vlx(2)',vlx(2)
c	write(*,*)'URX(2,1)',URX(2,1)
c	write(*,*)'URX(2,2)',URX(2,2)
c	write(*,*)'vlx(1)',vlx(1)
c	write(*,*)'vlx(2)',vlx(2)
c	write(*,*)'URX(2,1)',URX(2,1)
c	write(*,*)'URX(2,2)',URX(2,2)
c	stop
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.0*SM)
      URZZ(2)=(URZ(2,2)-WN*RAM*URZ(1,1))/(RAM+2.0*SM)
C
c   * the Bessel functions *
c
      RK=WN*RR
      BJ0(0)=BESJ0(RK)
      BJ1(0)=BESJ1(RK)
      IF(RK.EQ.0.0)THEN
        BJK(0)=0.5
      ELSE
        BJK(0)=BJ1(0)/RK
      END IF
      BJD(0)=BJ0(0)-BJK(0)
C
      BJ0(1)=-BJ1(0)*WN
      BJ1(1)= BJD(0)*WN
      IF(RK.EQ.0.0)THEN
        BJK(1)=0.0
      ELSE
        BJK(1)=(BJ1(1)-BJ1(0)/RR)/RK
      END IF
      BJD(1)= BJ0(1)-BJK(1)
C
c   * the integrands *
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(30,1100)RK,QX(0),VX(0),WX(0),QZ(0),WZ(0)
	*                  ,QX(1),VX(1),WX(1),QZ(1),WZ(1)
	*                  ,QX(2),VX(2),WX(2),QZ(2),WZ(2)
      END IF
 1100 FORMAT(31E12.4)
C
      RETURN
      END
C
C
C *** THE COMPUTATION OF Q, V, AND W AND THEIR DERIVATIVES 
C     USING Filon'S INTEGRATION FORMULA (after k*r > 8.0) ***
C
      SUBROUTINE QVWRZF0(OM,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)
     *      ,QXSA(0:2),QZSA(0:2),VXSA(0:2),WXSA(0:2),WZSA(0:2)
     *      ,QXCA(0:2),QZCA(0:2),VXCA(0:2),WXCA(0:2),WZCA(0:2)
     *      ,QXSB(0:2),QZSB(0:2),VXSB(0:2),WXSB(0:2),WZSB(0:2)
     *      ,QXCB(0:2),QZCB(0:2),VXCB(0:2),WXCB(0:2),WZCB(0:2)
     *      ,QXSC(0:2),QZSC(0:2),VXSC(0:2),WXSC(0:2),WZSC(0:2)
     *      ,QXCC(0:2),QZCC(0:2),VXCC(0:2),WXCC(0:2),WZCC(0:2)
     *      ,QXS1(0:2),QZS1(0:2),VXS1(0:2),WXS1(0:2),WZS1(0:2)
     *      ,QXC1(0:2),QZC1(0:2),VXC1(0:2),WXC1(0:2),WZC1(0:2)
     *      ,SM,RAM,S1,S2,OM
      COMMON /FCQVW/SM,RAM,RR  /IRANGE/NRNG,IRNG
     *       /INTG/VPMAX,VSMAX,RYMIN,CFINAL,N11,N12,M11
     *       /WVNUM/WVN(2001),DWN,WNR(3),NWN
	*       /sInteg/sImax1,sImax2
C
C *** Dimesion Check ***
C
      NC=2001
      IF(NWN.GT.NC)THEN
         WRITE(10,*)' FROM QVWRZF: NWN.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C   * RANGE 1 : From near Zero to near om/Rayleigh(min) *
C
c      write(2,*)
c
      IF(NRNG.EQ.3)RETURN

      JRNG=0
      IF(NRNG.GT.1)GOTO 10
C
      K=2*IRNG+1
      WA=WVN(K)
c      write(*,*)'QVWRZF(Range 1):WA*RR',WA*RR
      CALL INTEGF0(OM,K,WA,QXSA,QZSA,VXSA,WXSA,WZSA
     *                    ,QXCA,QZCA,VXCA,WXCA,WZCA)
c
      DO  300 I=IRNG+1,M11
c
        K=2*I
        WB=WVN(K)
        CALL INTEGF0(OM,K,WB,QXSB,QZSB,VXSB,WXSB,WZSB
     *                      ,QXCB,QZCB,VXCB,WXCB,WZCB)
        K=K+1
        WC=WVN(K)
        CALL INTEGF0(OM,K,WC,QXSC,QZSC,VXSC,WXSC,WZSC
     *                      ,QXCC,QZCC,VXCC,WXCC,WZCC)
c
        H=WC-WB
        th=RR*H
        t=1.d0/th
        t2=t*t
        t3=t2*t
        si=dsin(th)
        ss=si*si
        ci=dcos(th)
        cc=ci*ci
        s2=dsin(2.d0*th)
c
        a1=t+0.5d0*s2*t2-2.d0*ss*t3
        b1=2.d0*((1.d0+cc)*t2-s2*t3)
        c1=4.d0*(si*t3-ci*t2)
c
        SA=1.D0/DTAN(RR*WA)
        CA=DTAN(RR*WA)
        SC=1.D0/DTAN(RR*WC)
        CC=DTAN(RR*WC)
        DO 350 J=0,2
          QXS1(J)=(QXSA(J)+QXSC(J))*0.5d0
          QZS1(J)=(QZSA(J)+QZSC(J))*0.5d0
          VXS1(J)=(VXSA(J)+VXSC(J))*0.5d0
          WXS1(J)=(WXSA(J)+WXSC(J))*0.5d0
          WZS1(J)=(WZSA(J)+WZSC(J))*0.5d0
          QXC1(J)=(QXCA(J)+QXCC(J))*0.5d0
          QZC1(J)=(QZCA(J)+QZCC(J))*0.5d0
          VXC1(J)=(VXCA(J)+VXCC(J))*0.5d0
          WXC1(J)=(WXCA(J)+WXCC(J))*0.5d0
          WZC1(J)=(WZCA(J)+WZCC(J))*0.5d0
c
          S1=A1*(QXSA(J)*SA-QXSC(J)*SC)+B1*QXS1(J)+C1*QXSB(J)
          S2=A1*(QXCC(J)*CC-QXCA(J)*CA)+B1*QXC1(J)+C1*QXCB(J)
          GQX(J)=GQX(J)+H*(S1+S2)
          S1=A1*(QZSA(J)*SA-QZSC(J)*SC)+B1*QZS1(J)+C1*QZSB(J)
          S2=A1*(QZCC(J)*CC-QZCA(J)*CA)+B1*QZC1(J)+C1*QZCB(J)
          GQZ(J)=GQZ(J)+H*(S1+S2)
          S1=A1*(VXSA(J)*SA-VXSC(J)*SC)+B1*VXS1(J)+C1*VXSB(J)
          S2=A1*(VXCC(J)*CC-VXCA(J)*CA)+B1*VXC1(J)+C1*VXCB(J)
          GVX(J)=GVX(J)+H*(S1+S2)
          S1=A1*(WXSA(J)*SA-WXSC(J)*SC)+B1*WXS1(J)+C1*WXSB(J)
          S2=A1*(WXCC(J)*CC-WXCA(J)*CA)+B1*WXC1(J)+C1*WXCB(J)
          GWX(J)=GWX(J)+H*(S1+S2)
          S1=A1*(WZSA(J)*SA-WZSC(J)*SC)+B1*WZS1(J)+C1*WZSB(J)
          S2=A1*(WZCC(J)*CC-WZCA(J)*CA)+B1*WZC1(J)+C1*WZCB(J)
          GWZ(J)=GWZ(J)+H*(S1+S2)
c
          QXSA(J)=QXSC(J)
          QZSA(J)=QZSC(J)
          VXSA(J)=VXSC(J)
          WXSA(J)=WXSC(J)
          WZSA(J)=WZSC(J)
          QXCA(J)=QXCC(J)
          QZCA(J)=QZCC(J)
          VXCA(J)=VXCC(J)
          WXCA(J)=WXCC(J)
          WZCA(J)=WZCC(J)
c
          WA=WC
  350   CONTINUE
c
c     * Check	of Convergence of Integrands *
c
        S1=A1*(QZSA(2)*SA-QZSC(2)*SC)+B1*QZS1(2)+C1*QZSB(2)
        S2=A1*(QZCC(2)*CC-QZCA(2)*CA)+B1*QZC1(2)+C1*QZCB(2)
	  sI11=cdabs(S1+S2)
        S1=A1*(WZSA(2)*SA-WZSC(2)*SC)+B1*WZS1(2)+C1*WZSB(2)
        S2=A1*(WZCC(2)*CC-WZCA(2)*CA)+B1*WZC1(2)+C1*WZCB(2)
	  sI12=cdabs(S1+S2)
	  if(sI11.ge.sImax1)sImax1=sI11
	  if(sI12.ge.sImax2)sImax2=sI12
c	write(2,999)k,si11,simax1,si11/simax1,si12,simax2,si12/simax2
c 999  format(i5,6e12.4)
	  if(((sI11/sImax1).lt.0.00001d0).and.
     *                       ((sI12/sImax2).lt.0.00001d0))then
c	    write(*,*)' Filon Range 1: k*R & Int.Num.=',WC*RR,K
c	    write(*,*)'sI11,sI12',sI11,sI12
c	    write(*,*)'sImax1,sImax2',sImax1,sImax2
c	    write(*,*)'sI11/sImax1 < 0.00001d !: sI11/sImax1=',sI11/sImax1
c	    write(*,*)'sI12/sImax2 < 0.00001d !: sI12/sImax2=',sI12/sImax2
	    goto 20
	  end if
c
  300 CONTINUE

c      write(*,*)'f1:gqx0',gqx(0)
C
C   * RANGE 2 : From om/Rayleigh(min) to om/C(final) *
C
      JRNG=1
      IRNG=0
C
   10 H=DWN
      th=RR*H
      t=1.d0/th
      t2=t*t
      t3=t2*t
      si=dsin(th)
      ss=si*si
      ci=dcos(th)
      cc=ci*ci
      s2=dsin(2.d0*th)
      a1=t+0.5d0*s2*t2-2.d0*ss*t3
      b1=2.d0*((1.d0+cc)*t2-s2*t3)
      c1=4.d0*(si*t3-ci*t2)
c
      K=2*(M11+IRNG)+1
      WA=WVN(K)
c      write(*,*)'QVWRZF(Range 2):WA*RR',WA*RR
      IF(JRNG.EQ.1)GOTO 15
      CALL INTEGF0(OM,K,WA,QXSA,QZSA,VXSA,WXSA,WZSA
     *                    ,QXCA,QZCA,VXCA,WXCA,WZCA)
c
   15 DO  400 I=IRNG+1,N12
        K=2*(M11+I)
        WB=WVN(K)
        CALL INTEGF0(OM,K,WB,QXSB,QZSB,VXSB,WXSB,WZSB
     *                      ,QXCB,QZCB,VXCB,WXCB,WZCB)
        K=K+1
        WC=WVN(K)
        CALL INTEGF0(OM,K,WC,QXSC,QZSC,VXSC,WXSC,WZSC
     *                      ,QXCC,QZCC,VXCC,WXCC,WZCC)
c
        SA=1.D0/DTAN(RR*WA)
        CA=DTAN(RR*WA)
        SC=1.D0/DTAN(RR*WC)
        CC=DTAN(RR*WC)
        DO 450 J=0,2
          QXS1(J)=(QXSA(J)+QXSC(J))*0.5d0
          QZS1(J)=(QZSA(J)+QZSC(J))*0.5d0
          VXS1(J)=(VXSA(J)+VXSC(J))*0.5d0
          WXS1(J)=(WXSA(J)+WXSC(J))*0.5d0
          WZS1(J)=(WZSA(J)+WZSC(J))*0.5d0
          QXC1(J)=(QXCA(J)+QXCC(J))*0.5d0
          QZC1(J)=(QZCA(J)+QZCC(J))*0.5d0
          VXC1(J)=(VXCA(J)+VXCC(J))*0.5d0
          WXC1(J)=(WXCA(J)+WXCC(J))*0.5d0
          WZC1(J)=(WZCA(J)+WZCC(J))*0.5d0
c
          S1=A1*(QXSA(J)*SA-QXSC(J)*SC)+B1*QXS1(J)+C1*QXSB(J)
          S2=A1*(QXCC(J)*CC-QXCA(J)*CA)+B1*QXC1(J)+C1*QXCB(J)
          GQX(J)=GQX(J)+H*(S1+S2)
          S1=A1*(QZSA(J)*SA-QZSC(J)*SC)+B1*QZS1(J)+C1*QZSB(J)
          S2=A1*(QZCC(J)*CC-QZCA(J)*CA)+B1*QZC1(J)+C1*QZCB(J)
          GQZ(J)=GQZ(J)+H*(S1+S2)
          S1=A1*(VXSA(J)*SA-VXSC(J)*SC)+B1*VXS1(J)+C1*VXSB(J)
          S2=A1*(VXCC(J)*CC-VXCA(J)*CA)+B1*VXC1(J)+C1*VXCB(J)
          GVX(J)=GVX(J)+H*(S1+S2)
          S1=A1*(WXSA(J)*SA-WXSC(J)*SC)+B1*WXS1(J)+C1*WXSB(J)
          S2=A1*(WXCC(J)*CC-WXCA(J)*CA)+B1*WXC1(J)+C1*WXCB(J)
          GWX(J)=GWX(J)+H*(S1+S2)
          S1=A1*(WZSA(J)*SA-WZSC(J)*SC)+B1*WZS1(J)+C1*WZSB(J)
          S2=A1*(WZCC(J)*CC-WZCA(J)*CA)+B1*WZC1(J)+C1*WZCB(J)
          GWZ(J)=GWZ(J)+H*(S1+S2)
c
          QXSA(J)=QXSC(J)
          QZSA(J)=QZSC(J)
          VXSA(J)=VXSC(J)
          WXSA(J)=WXSC(J)
          WZSA(J)=WZSC(J)
          QXCA(J)=QXCC(J)
          QZCA(J)=QZCC(J)
          VXCA(J)=VXCC(J)
          WXCA(J)=WXCC(J)
          WZCA(J)=WZCC(J)
c
          WA=WC
  450   CONTINUE
c
c     * Check	of Convergence of Integrands *
c
        S1=A1*(QZSA(2)*SA-QZSC(2)*SC)+B1*QZS1(2)+C1*QZSB(2)
        S2=A1*(QZCC(2)*CC-QZCA(2)*CA)+B1*QZC1(2)+C1*QZCB(2)
	  sI11=cdabs(S1+S2)
        S1=A1*(WZSA(2)*SA-WZSC(2)*SC)+B1*WZS1(2)+C1*WZSB(2)
        S2=A1*(WZCC(2)*CC-WZCA(2)*CA)+B1*WZC1(2)+C1*WZCB(2)
	  sI12=cdabs(S1+S2)
	  if(sI11.ge.sImax1)sImax1=sI11
	  if(sI12.ge.sImax2)sImax2=sI12
c	write(2,999)k,si11,simax1,si11/simax1,si12,simax2,si12/simax2
	  if(((sI11/sImax1).lt.0.00001d0).and.
     *                       ((sI12/sImax2).lt.0.00001d0))then
c	    write(*,*)' Filon Range 2: k*R & Int.Num.=',WC*RR,K
c	    write(*,*)'sI11,sI12',sI11,sI12
c	    write(*,*)'sImax1,sImax2',sImax1,sImax2
c	    write(*,*)'sI11/sImax1 < 0.00001d !: sI11/sImax1=',sI11/sImax1
c	    write(*,*)'sI12/sImax2 < 0.00001d !: sI12/sImax2=',sI12/sImax2
	    goto 20
	    goto 20
	  end if
c
  400 CONTINUE

c      write(*,*)'f2:gqx0',gqx(0)
c
   20 RETURN
      END
C
C
C *** THE INTEGRANDS OF Q, V, AND W AND THEIR DERIVATIVES for Filon's ***
C
      SUBROUTINE INTEGF0(OM,K,WN,QXS,QZS,VXS,WXS,WZS
     *                          ,QXC,QZC,VXC,WXC,WZC)
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)
     *   ,QXC(0:2),QZC(0:2),VXC(0:2),WXC(0:2),WZC(0:2)
     *   ,VLX(2),URX(2,2),URZ(2,2),VLX0(2),URX0(2,2),URZ0(2,2)
     *   ,VLXZ,URXZ(2),URZZ(2),SM,RAM,OM
      DIMENSION BJ0S(0:1),BJ1S(0:1),BJKS(0:1),BJDS(0:1)
     *         ,BJ0C(0:1),BJ1C(0:1),BJKC(0:1),BJDC(0:1)
      COMMON /FREQ/DOM,IOM,NOM,OMI /SOLY/SDEP,ODEP,LS,LO 
     *       /STRC/DNS(20),VP(20),QP(20),FP(20)
     *                    ,VS(20),QS(20),FS(20),THK(20),NL
     *       /OTST/IOT,IST /FCQVW/SM,RAM,RR
C
C *** CHECK OF DIMENSION ***
C
      NC=20
      IF(NL.GE.NC)THEN
        WRITE(10,*)' FROM INTEGF:NL>NC !!!!!!!!!!'
        WRITE(10,*)' NL,NC=',NL,NC
        STOP
      END IF
C
c   * the displacement-stress vectors *
c
      if(IST.eq.0)then
        CALL RTSH(OM,K,WN,VLX)
        CALL RTPSV(OM,K,WN,URX,URZ)
        CALL RTSH0(dcmplx(0.d0,OMI),K,WN,VLX0)
        CALL RTPSV0(dcmplx(0.d0,OMI),K,WN,URX0,URZ0)
        VLX(1)=VLX(1)-VLX0(1)
        URX(1,1)=URX(1,1)-URX0(1,1)
        URX(2,1)=URX(2,1)-URX0(2,1)
        URZ(1,1)=URZ(1,1)-URZ0(1,1)
        URZ(2,1)=URZ(2,1)-URZ0(2,1)
        VLX(2)=VLX(2)-VLX0(2)
        URX(1,2)=URX(1,2)-URX0(1,2)
        URX(2,2)=URX(2,2)-URX0(2,2)
        URZ(1,2)=URZ(1,2)-URZ0(1,2)
        URZ(2,2)=URZ(2,2)-URZ0(2,2)
	else
        if(IST.eq.1)then
          CALL RTSH(OM,K,WN,VLX)
          CALL RTPSV(OM,K,WN,URX,URZ)
	  else
          if(IST.eq.2)then
            CALL RTSH0(dcmplx(0.d0,OMI),K,WN,VLX)
            CALL RTPSV0(dcmplx(0.d0,OMI),K,WN,URX,URZ)
	    end if
	  end if
      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.0*SM)
      URZZ(2)=(URZ(2,2)-WN*RAM*URZ(1,1))/(RAM+2.0*SM)
C
c   * SIN and COS parts of the approximate Bessel functions *
c
      PI=4.D0*DATAN(1.D0)
      RK=WN*RR
      A=1.D0/DSQRT(PI*RK)
      SRK=A*DSIN(RK)
      CRK=A*DCOS(RK)
c
      x1=0.125d0/rk
      x2=0.0703125d0/(rk*rk)
      x3=0.073242188d0/(rk*rk*rk)      
      BJ0S(0)= SRK*(1.D0+x1-x2-x3)
      BJ0C(0)= CRK*(1.D0-x1-x2+x3)

      y1=0.375d0/rk
      y2=0.1171875d0/(rk*rk)
      y3=0.102539063d0/(rk*rk*rk)
      BJ1S(0)= SRK*(1.D0+y1+y2-y3)
      BJ1C(0)=-CRK*(1.D0-y1+y2+y3)
C
      BJKS(0)=BJ1S(0)/RK
      BJKC(0)=BJ1C(0)/RK
      BJDS(0)=BJ0S(0)-BJKS(0)
      BJDC(0)=BJ0C(0)-BJKC(0)
C
      BJ0S(1)=-BJ1S(0)*WN
      BJ0C(1)=-BJ1C(0)*WN
      BJ1S(1)= BJDS(0)*WN
      BJ1C(1)= BJDC(0)*WN

      BJKS(1)=(BJ1S(1)-BJ1S(0)/RR)/RK
      BJKC(1)=(BJ1C(1)-BJ1C(0)/RR)/RK
      BJDS(1)= BJ0S(1)-BJKS(1)
      BJDC(1)= BJ0C(1)-BJKC(1)
C
c   * the integrands *
c
      QXS(0)= URX(1,1)*BJDS(0)+VLX(1)*BJKS(0)
      QXC(0)= URX(1,1)*BJDC(0)+VLX(1)*BJKC(0)
      VXS(0)= URX(1,1)*BJKS(0)+VLX(1)*BJDS(0)
      VXC(0)= URX(1,1)*BJKC(0)+VLX(1)*BJDC(0)
      WXS(0)= URX(2,1)*BJ1S(0)
      WXC(0)= URX(2,1)*BJ1C(0)
      QZS(0)= URZ(1,1)*BJ1S(0)
      QZC(0)= URZ(1,1)*BJ1C(0)
      WZS(0)= URZ(2,1)*BJ0S(0)
      WZC(0)= URZ(2,1)*BJ0C(0)
C
      QXS(1)= URX(1,1)*BJDS(1)+VLX(1)*BJKS(1)
      QXC(1)= URX(1,1)*BJDC(1)+VLX(1)*BJKC(1)
      VXS(1)= URX(1,1)*BJKS(1)+VLX(1)*BJDS(1)
      VXC(1)= URX(1,1)*BJKC(1)+VLX(1)*BJDC(1)
      WXS(1)= URX(2,1)*BJ1S(1)
      WXC(1)= URX(2,1)*BJ1C(1)
      QZS(1)= URZ(1,1)*BJ1S(1)
      QZC(1)= URZ(1,1)*BJ1C(1)
      WZS(1)= URZ(2,1)*BJ0S(1)
      WZC(1)= URZ(2,1)*BJ0C(1)
C
      QXS(2)= URXZ(1)*BJDS(0)+VLXZ*BJKS(0)
      QXC(2)= URXZ(1)*BJDC(0)+VLXZ*BJKC(0)
      VXS(2)= URXZ(1)*BJKS(0)+VLXZ*BJDS(0)
      VXC(2)= URXZ(1)*BJKC(0)+VLXZ*BJDC(0)
      WXS(2)= URXZ(2)*BJ1S(0)
      WXC(2)= URXZ(2)*BJ1C(0)
      QZS(2)= URZZ(1)*BJ1S(0)
      QZC(2)= URZZ(1)*BJ1C(0)
      WZS(2)= URZZ(2)*BJ0S(0)
      WZC(2)= URZZ(2)*BJ0C(0)
C
      IF(IOT.eq.1)THEN
        WRITE(30,1100)RK,QXS(0)+QXC(0),VXS(0)+VXC(0),WXS(0)+WXC(0)
	*    ,QZS(0)+QZC(0),WZS(0)+WZC(0),QXS(1)+QXC(1),VXS(1)+VXC(1)
     *    ,WXS(1)+WXC(1),QZS(1)+QZC(1),WZS(1)+WZC(1),QXS(2)+QXC(2)
     *    ,VXS(2)+VXC(2),WXS(2)+WXC(2),QZS(2)+QZC(2),WZS(2)+WZC(2)
      END IF
 1100 FORMAT(31E12.4)
C
      RETURN
      END
C
C
C
C *** Library for SH UPGOING/DOWNGOING COEFFICIENTS ***
C
      SUBROUTINE CUDSH0(OM)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VSc,CD,CU,VS2,WS2,CMS,SM,EE(2,2),DD,DU
     *  ,CUSU,CDSU,CDSD,CUSD,GTD,GTU,GRD,GRU,AU,BU,AD,BD,CS,OM,OM2
      COMMON /SOLY/SDEP,ODEP,LS,LO 
     *       /STRC/DNS(20),VP(20),QP(20),FP(20)
     *                    ,VS(20),QS(20),FS(20),THK(20),NL
     *       /DUSH0/CD(2001,20),CU(2001,20)
     *       ,CUSU(2001),CDSU(2001),CDSD(2001),CUSD(2001) 
     *       /GRTSH/GTD(0:20),GTU(0:20),GRD(0:20),GRU(0:20)
     *       /WVNUM/WVN(2001),DWN,WNR(3),NWN
C
C *** CHECK FOR THE DIMENSION OF CD, CU, AND GTD-GRU ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(10,*)' FROM CUDSH: NL.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
      NC=2001
      IF(NWN.GT.NC)THEN
         WRITE(10,*)' FROM CUDSH: NWN.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C *** DOWN/UP-R/T COEFFICIENTS FOR ALL LAYERS AT ALL WAVE-NUMBERS ***
C
      OM2=OM*OM
      VSc=VS(LS)*DCMPLX(1.D0,-0.5D0/QS(LS))
      VS2=VSc*VSc
      WS2=OM2/VS2
      SM=DNS(LS)*VS2
      EE(1,1)=(1.0,0.0)
      EE(1,2)=(1.0,0.0)
C
      DO 100 I=1,NWN
C
        WN=WVN(I)
C
C     * COMPUTATION OF GENERALIZED R/T MATRICES *
C
        CALL GNRTSH(OM,WN)
C
C     * COMPUTATION OF  CD AND CU *
C
        WN2=WN*WN
        CMS=CDSQRT(WN2-WS2)
        IF(dREAL(CMS).LT.0.0)CMS=-CMS
C
        EE(2,1)=-SM*CMS
        EE(2,2)= SM*CMS
C
        DQ=-WN/6.2831853072d0
C
        IF((LS.EQ.1).AND.(SDEP.EQ.0.0))THEN
C
          IF(NL.EQ.1)THEN
            CDSD(I)=DQ/EE(2,1)
            CUSD(I)=(0.0,0.0)
          ELSE
            DD=CDEXP(-CMS*THK(1))
            CDSD(I)=DQ/(EE(2,1)+EE(2,2)*DD*GRD(1))
            CUSD(I)=GRD(1)*CDSD(I)
          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(I)=DQ/(BD*CS-BU)
            CDSU(I)=GRU(LS-1)*CUSU(I)
            CDSD(I)=CS*CUSU(I)
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(I)=DQ/(BD*CS-BU)
            CDSU(I)=GRU(LS-1)*CUSU(I)
            CDSD(I)=CS*CUSU(I)
            CUSD(I)=GRD(LS)*CDSD(I)
C
          ENDIF
C
        END IF
C
        IF(LS.EQ.1)GOTO 10
        CU(I,LS-1)=GTU(LS-1)*CUSU(I)
        DO 110 J=LS-2,1,-1
          CU(I,J)=GTU(J)*CU(I,J+1)
          CD(I,J+1)=GRU(J)*CU(I,J+1)
  110   CONTINUE
        CD(I,1)=GRU(0)*CU(I,1)
C
   10   IF(LS.EQ.NL)GOTO 100
        CD(I,LS+1)=GTD(LS)*CDSD(I)
        DO 120 J=LS+1,NL-1
          CD(I,J+1)=GTD(J)*CD(I,J)
          CU(I,J)=GRD(J)*CD(I,J)
  120   CONTINUE
C
  100 CONTINUE
C
      RETURN
      END
C
C *** Library for SH UPGOING/DOWNGOING COEFFICIENTS ***
C
      SUBROUTINE CUDSH(OM)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VSc,CD,CU,VS2,WS2,CMS,SM,EE(2,2),DD,DU
     *  ,CUSU,CDSU,CDSD,CUSD,GTD,GTU,GRD,GRU,AU,BU,AD,BD,CS,OM,OM2
      COMMON /SOLY/SDEP,ODEP,LS,LO 
     *       /STRC/DNS(20),VP(20),QP(20),FP(20)
     *                    ,VS(20),QS(20),FS(20),THK(20),NL
     *       /DUSH/CD(2001,20),CU(2001,20)
     *       ,CUSU(2001),CDSU(2001),CDSD(2001),CUSD(2001) 
     *       /GRTSH/GTD(0:20),GTU(0:20),GRD(0:20),GRU(0:20)
     *       /WVNUM/WVN(2001),DWN,WNR(3),NWN
C
C *** CHECK FOR THE DIMENSION OF CD, CU, AND GTD-GRU ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(10,*)' FROM CUDSH: NL.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
      NC=2001
      IF(NWN.GT.NC)THEN
         WRITE(10,*)' FROM CUDSH: NWN.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C *** DOWN/UP-R/T COEFFICIENTS FOR ALL LAYERS AT ALL WAVE-NUMBERS ***
C
      OM2=OM*OM
      freq=om/6.2831853072d0
      QSe=QS(LS)*(freq)**FS(LS)
      VSc=VS(LS)*DCMPLX(1.D0,-0.5D0/QSe)
      VS2=VSc*VSc
      WS2=OM2/VS2
      SM=DNS(LS)*VS2
      EE(1,1)=(1.0,0.0)
      EE(1,2)=(1.0,0.0)
C
      DO 100 I=1,NWN
C
        WN=WVN(I)
C
C     * COMPUTATION OF GENERALIZED R/T MATRICES *
C
        CALL GNRTSH(OM,WN)
C
C     * COMPUTATION OF  CD AND CU *
C
        WN2=WN*WN
        CMS=CDSQRT(WN2-WS2)
        IF(dREAL(CMS).LT.0.0)CMS=-CMS
C
        EE(2,1)=-SM*CMS
        EE(2,2)= SM*CMS
C
        DQ=-WN/6.2831853072d0
C
        IF((LS.EQ.1).AND.(SDEP.EQ.0.0))THEN
C
          IF(NL.EQ.1)THEN
            CDSD(I)=DQ/EE(2,1)
            CUSD(I)=(0.0,0.0)
          ELSE
            DD=CDEXP(-CMS*THK(1))
            CDSD(I)=DQ/(EE(2,1)+EE(2,2)*DD*GRD(1))
            CUSD(I)=GRD(1)*CDSD(I)
          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(I)=DQ/(BD*CS-BU)
            CDSU(I)=GRU(LS-1)*CUSU(I)
            CDSD(I)=CS*CUSU(I)
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(I)=DQ/(BD*CS-BU)
            CDSU(I)=GRU(LS-1)*CUSU(I)
            CDSD(I)=CS*CUSU(I)
            CUSD(I)=GRD(LS)*CDSD(I)
C
          ENDIF
C
        END IF
C
        IF(LS.EQ.1)GOTO 10
        CU(I,LS-1)=GTU(LS-1)*CUSU(I)
        DO 110 J=LS-2,1,-1
          CU(I,J)=GTU(J)*CU(I,J+1)
          CD(I,J+1)=GRU(J)*CU(I,J+1)
  110   CONTINUE
        CD(I,1)=GRU(0)*CU(I,1)
C
   10   IF(LS.EQ.NL)GOTO 100
        CD(I,LS+1)=GTD(LS)*CDSD(I)
        DO 120 J=LS+1,NL-1
          CD(I,J+1)=GTD(J)*CD(I,J)
          CU(I,J)=GRD(J)*CD(I,J)
  120   CONTINUE
C
  100 CONTINUE
C
      RETURN
      END
C
C
C *** Library for P-SV UPGOING/DOWNGOING COEFFICIENTS ***
C
      SUBROUTINE CUDPSV0(OM)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 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,OM
     *    ,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),QP(20),FP(20)
     *                 ,VS(20),QS(20),FS(20),THK(20),NL
     *  /DUPSVX0/CDX(2001,2,20),CUX(2001,2,20),CUSUX(2001,2)
     *         ,CDSUX(2001,2),CDSDX(2001,2),CUSDX(2001,2) 
     *  /DUPSVZ0/CDZ(2001,2,20),CUZ(2001,2,20),CUSUZ(2001,2)
     *         ,CDSUZ(2001,2),CDSDZ(2001,2),CUSDZ(2001,2) 
     *  /GRTPSV/GTD(2,2,0:20),GTU(2,2,0:20),GRD(2,2,0:20),GRU(2,2,0:20)
     *  /WVNUM/WVN(2001),DWN,WNR(3),NWN
C
C *** CHECK FOR THE DIMENSION OF CD, CU, AND GTD-GRU ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(10,*)' FROM CUDPSV: NL.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
      NC=2001
      IF(NWN.GT.NC)THEN
         WRITE(10,*)' FROM CUDPSV: NWN.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C *** COMPUTATION OF GENERALIZED R/T MATRICES ***
C
      DO 100 K=1,NWN
C
        WN=WVN(K)
        CALL GNRTPSV(OM,WN)
C
C *** COMPUTATION OF CD AND CU ***
C
        CALL MTXPSV(OM,WN,LS,CMS,CMP,EE)
C
        DQX=-WN/6.2831853072d0
        DQZ= WN/6.2831853072d0
C
        IF((LS.EQ.1).AND.(SDEP.EQ.0.0))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(K,1)=E2(1,1)*DQX
            CDSDX(K,2)=E2(2,1)*DQX
            CUSDX(K,1)=(0.0,0.0)
            CUSDX(K,2)=(0.0,0.0)
            CDSDZ(K,1)=E2(1,2)*DQZ
            CDSDZ(K,2)=E2(2,2)*DQZ
            CUSDZ(K,1)=(0.0,0.0)
            CUSDZ(K,2)=(0.0,0.0)
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(K,1)=E2(1,1)*DQX
            CDSDX(K,2)=E2(2,1)*DQX
            CUSDX(K,1)=GRD(1,1,1)*CDSDX(K,1)
     *                  +GRD(1,2,1)*CDSDX(K,2)
            CUSDX(K,2)=GRD(2,1,1)*CDSDX(K,1)
     *                  +GRD(2,2,1)*CDSDX(K,2)
            CDSDZ(K,1)=E2(1,2)*DQZ
            CDSDZ(K,2)=E2(2,2)*DQZ
            CUSDZ(K,1)=GRD(1,1,1)*CDSDZ(K,1)
     *                  +GRD(1,2,1)*CDSDZ(K,2)
            CUSDZ(K,2)=GRD(2,1,1)*CDSDZ(K,1)
     *                  +GRD(2,2,1)*CDSDZ(K,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(K,1)=E2(1,1)*DQX
            CUSUX(K,2)=E2(2,1)*DQX
            CDSUX(K,1)=GRU(1,1,LS-1)*CUSUX(K,1)
     *                  +GRU(1,2,LS-1)*CUSUX(K,2)
            CDSUX(K,2)=GRU(2,1,LS-1)*CUSUX(K,1)
     *                  +GRU(2,2,LS-1)*CUSUX(K,2)
            CDSDX(K,1)=CS(1,1)*CUSUX(K,1)+CS(1,2)*CUSUX(K,2)
            CDSDX(K,2)=CS(2,1)*CUSUX(K,1)+CS(2,2)*CUSUX(K,2)
C
            CUSUZ(K,1)=E2(1,2)*DQZ
            CUSUZ(K,2)=E2(2,2)*DQZ
            CDSUZ(K,1)=GRU(1,1,LS-1)*CUSUZ(K,1)
     *                  +GRU(1,2,LS-1)*CUSUZ(K,2)
            CDSUZ(K,2)=GRU(2,1,LS-1)*CUSUZ(K,1)
     *                  +GRU(2,2,LS-1)*CUSUZ(K,2)
            CDSDZ(K,1)=CS(1,1)*CUSUZ(K,1)+CS(1,2)*CUSUZ(K,2)
            CDSDZ(K,2)=CS(2,1)*CUSUZ(K,1)+CS(2,2)*CUSUZ(K,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(K,1)=E2(1,1)*DQX
            CUSUX(K,2)=E2(2,1)*DQX
            CDSUX(K,1)=GRU(1,1,LS-1)*CUSUX(K,1)
     *                  +GRU(1,2,LS-1)*CUSUX(K,2)
            CDSUX(K,2)=GRU(2,1,LS-1)*CUSUX(K,1)
     *                  +GRU(2,2,LS-1)*CUSUX(K,2)
            CDSDX(K,1)=CS(1,1)*CUSUX(K,1)+CS(1,2)*CUSUX(K,2)
            CDSDX(K,2)=CS(2,1)*CUSUX(K,1)+CS(2,2)*CUSUX(K,2)
            CUSDX(K,1)=GRD(1,1,LS)*CDSDX(K,1)
     *                  +GRD(1,2,LS)*CDSDX(K,2)
            CUSDX(K,2)=GRD(2,1,LS)*CDSDX(K,1)
     *                  +GRD(2,2,LS)*CDSDX(K,2)
C
            CUSUZ(K,1)=E2(1,2)*DQZ
            CUSUZ(K,2)=E2(2,2)*DQZ
            CDSUZ(K,1)=GRU(1,1,LS-1)*CUSUZ(K,1)
     *                  +GRU(1,2,LS-1)*CUSUZ(K,2)
            CDSUZ(K,2)=GRU(2,1,LS-1)*CUSUZ(K,1)
     *                  +GRU(2,2,LS-1)*CUSUZ(K,2)
            CDSDZ(K,1)=CS(1,1)*CUSUZ(K,1)+CS(1,2)*CUSUZ(K,2)
            CDSDZ(K,2)=CS(2,1)*CUSUZ(K,1)+CS(2,2)*CUSUZ(K,2)
            CUSDZ(K,1)=GRD(1,1,LS)*CDSDZ(K,1)
     *                  +GRD(1,2,LS)*CDSDZ(K,2)
            CUSDZ(K,2)=GRD(2,1,LS)*CDSDZ(K,1)
     *                  +GRD(2,2,LS)*CDSDZ(K,2)
C
          ENDIF
C
        END IF
C
        IF(LS.EQ.1)GOTO 10
          CUX(K,1,LS-1)=GTU(1,1,LS-1)*CUSUX(K,1)
     *                   +GTU(1,2,LS-1)*CUSUX(K,2)
          CUX(K,2,LS-1)=GTU(2,1,LS-1)*CUSUX(K,1)
     *                   +GTU(2,2,LS-1)*CUSUX(K,2)
          CUZ(K,1,LS-1)=GTU(1,1,LS-1)*CUSUZ(K,1)
     *                   +GTU(1,2,LS-1)*CUSUZ(K,2)
          CUZ(K,2,LS-1)=GTU(2,1,LS-1)*CUSUZ(K,1)
     *                   +GTU(2,2,LS-1)*CUSUZ(K,2)
        DO 110 I=LS-2,1,-1
          CUX(K,1,I)  =GTU(1,1,I)*CUX(K,1,I+1)
     *                  +GTU(1,2,I)*CUX(K,2,I+1)
          CUX(K,2,I)  =GTU(2,1,I)*CUX(K,1,I+1)
     *                  +GTU(2,2,I)*CUX(K,2,I+1)
          CDX(K,1,I+1)=GRU(1,1,I)*CUX(K,1,I+1)
     *                  +GRU(1,2,I)*CUX(K,2,I+1)
          CDX(K,2,I+1)=GRU(2,1,I)*CUX(K,1,I+1)
     *                  +GRU(2,2,I)*CUX(K,2,I+1)
          CUZ(K,1,I)  =GTU(1,1,I)*CUZ(K,1,I+1)
     *                  +GTU(1,2,I)*CUZ(K,2,I+1)
          CUZ(K,2,I)  =GTU(2,1,I)*CUZ(K,1,I+1)
     *                  +GTU(2,2,I)*CUZ(K,2,I+1)
          CDZ(K,1,I+1)=GRU(1,1,I)*CUZ(K,1,I+1)
     *                  +GRU(1,2,I)*CUZ(K,2,I+1)
          CDZ(K,2,I+1)=GRU(2,1,I)*CUZ(K,1,I+1)
     *                  +GRU(2,2,I)*CUZ(K,2,I+1)
  110   CONTINUE
        CDX(K,1,1)=GRU(1,1,0)*CUX(K,1,1)+GRU(1,2,0)*CUX(K,2,1)
        CDX(K,2,1)=GRU(2,1,0)*CUX(K,1,1)+GRU(2,2,0)*CUX(K,2,1)
        CDZ(K,1,1)=GRU(1,1,0)*CUZ(K,1,1)+GRU(1,2,0)*CUZ(K,2,1)
        CDZ(K,2,1)=GRU(2,1,0)*CUZ(K,1,1)+GRU(2,2,0)*CUZ(K,2,1)
C
   10   IF(LS.EQ.NL)GOTO 100
        CDX(K,1,LS+1)=GTD(1,1,LS)*CDSDX(K,1)
     *                 +GTD(1,2,LS)*CDSDX(K,2)
        CDX(K,2,LS+1)=GTD(2,1,LS)*CDSDX(K,1)
     *                 +GTD(2,2,LS)*CDSDX(K,2)
        CDZ(K,1,LS+1)=GTD(1,1,LS)*CDSDZ(K,1)
     *                 +GTD(1,2,LS)*CDSDZ(K,2)
        CDZ(K,2,LS+1)=GTD(2,1,LS)*CDSDZ(K,1)
     *                 +GTD(2,2,LS)*CDSDZ(K,2)
        DO 120 I=LS+1,NL-1
          CDX(K,1,I+1)=GTD(1,1,I)*CDX(K,1,I)
     *                  +GTD(1,2,I)*CDX(K,2,I)
          CDX(K,2,I+1)=GTD(2,1,I)*CDX(K,1,I)
     *                  +GTD(2,2,I)*CDX(K,2,I)
          CUX(K,1,I)  =GRD(1,1,I)*CDX(K,1,I)
     *                  +GRD(1,2,I)*CDX(K,2,I)
          CUX(K,2,I)  =GRD(2,1,I)*CDX(K,1,I)
     *                  +GRD(2,2,I)*CDX(K,2,I)
          CDZ(K,1,I+1)=GTD(1,1,I)*CDZ(K,1,I)
     *                  +GTD(1,2,I)*CDZ(K,2,I)
          CDZ(K,2,I+1)=GTD(2,1,I)*CDZ(K,1,I)
     *                  +GTD(2,2,I)*CDZ(K,2,I)
          CUZ(K,1,I)  =GRD(1,1,I)*CDZ(K,1,I)
     *                  +GRD(1,2,I)*CDZ(K,2,I)
          CUZ(K,2,I)  =GRD(2,1,I)*CDZ(K,1,I)
     *                  +GRD(2,2,I)*CDZ(K,2,I)
  120   CONTINUE
C
  100 CONTINUE
C
      RETURN
      END
C
C
C *** Library for P-SV UPGOING/DOWNGOING COEFFICIENTS ***
C
      SUBROUTINE CUDPSV(OM)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 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,OM
     *    ,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),QP(20),FP(20)
     *                 ,VS(20),QS(20),FS(20),THK(20),NL
     *  /DUPSVX/CDX(2001,2,20),CUX(2001,2,20),CUSUX(2001,2)
     *         ,CDSUX(2001,2),CDSDX(2001,2),CUSDX(2001,2) 
     *  /DUPSVZ/CDZ(2001,2,20),CUZ(2001,2,20),CUSUZ(2001,2)
     *         ,CDSUZ(2001,2),CDSDZ(2001,2),CUSDZ(2001,2) 
     *  /GRTPSV/GTD(2,2,0:20),GTU(2,2,0:20),GRD(2,2,0:20),GRU(2,2,0:20)
     *  /WVNUM/WVN(2001),DWN,WNR(3),NWN
C
C *** CHECK FOR THE DIMENSION OF CD, CU, AND GTD-GRU ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(10,*)' FROM CUDPSV: NL.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
      NC=2001
      IF(NWN.GT.NC)THEN
         WRITE(10,*)' FROM CUDPSV: NWN.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C *** COMPUTATION OF GENERALIZED R/T MATRICES ***
C
      DO 100 K=1,NWN
C
        WN=WVN(K)
        CALL GNRTPSV(OM,WN)
C
C *** COMPUTATION OF CD AND CU ***
C
        CALL MTXPSV(OM,WN,LS,CMS,CMP,EE)
C
        DQX=-WN/6.2831853072d0
        DQZ= WN/6.2831853072d0
C
        IF((LS.EQ.1).AND.(SDEP.EQ.0.0))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(K,1)=E2(1,1)*DQX
            CDSDX(K,2)=E2(2,1)*DQX
            CUSDX(K,1)=(0.0,0.0)
            CUSDX(K,2)=(0.0,0.0)
            CDSDZ(K,1)=E2(1,2)*DQZ
            CDSDZ(K,2)=E2(2,2)*DQZ
            CUSDZ(K,1)=(0.0,0.0)
            CUSDZ(K,2)=(0.0,0.0)
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(K,1)=E2(1,1)*DQX
            CDSDX(K,2)=E2(2,1)*DQX
            CUSDX(K,1)=GRD(1,1,1)*CDSDX(K,1)
     *                  +GRD(1,2,1)*CDSDX(K,2)
            CUSDX(K,2)=GRD(2,1,1)*CDSDX(K,1)
     *                  +GRD(2,2,1)*CDSDX(K,2)
            CDSDZ(K,1)=E2(1,2)*DQZ
            CDSDZ(K,2)=E2(2,2)*DQZ
            CUSDZ(K,1)=GRD(1,1,1)*CDSDZ(K,1)
     *                  +GRD(1,2,1)*CDSDZ(K,2)
            CUSDZ(K,2)=GRD(2,1,1)*CDSDZ(K,1)
     *                  +GRD(2,2,1)*CDSDZ(K,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(K,1)=E2(1,1)*DQX
            CUSUX(K,2)=E2(2,1)*DQX
            CDSUX(K,1)=GRU(1,1,LS-1)*CUSUX(K,1)
     *                +GRU(1,2,LS-1)*CUSUX(K,2)
            CDSUX(K,2)=GRU(2,1,LS-1)*CUSUX(K,1)
     *                +GRU(2,2,LS-1)*CUSUX(K,2)
            CDSDX(K,1)=CS(1,1)*CUSUX(K,1)+CS(1,2)*CUSUX(K,2)
            CDSDX(K,2)=CS(2,1)*CUSUX(K,1)+CS(2,2)*CUSUX(K,2)
C
            CUSUZ(K,1)=E2(1,2)*DQZ
            CUSUZ(K,2)=E2(2,2)*DQZ
            CDSUZ(K,1)=GRU(1,1,LS-1)*CUSUZ(K,1)
     *                +GRU(1,2,LS-1)*CUSUZ(K,2)
            CDSUZ(K,2)=GRU(2,1,LS-1)*CUSUZ(K,1)
     *                +GRU(2,2,LS-1)*CUSUZ(K,2)
            CDSDZ(K,1)=CS(1,1)*CUSUZ(K,1)+CS(1,2)*CUSUZ(K,2)
            CDSDZ(K,2)=CS(2,1)*CUSUZ(K,1)+CS(2,2)*CUSUZ(K,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(K,1)=E2(1,1)*DQX
            CUSUX(K,2)=E2(2,1)*DQX
            CDSUX(K,1)=GRU(1,1,LS-1)*CUSUX(K,1)
     *                  +GRU(1,2,LS-1)*CUSUX(K,2)
            CDSUX(K,2)=GRU(2,1,LS-1)*CUSUX(K,1)
     *                  +GRU(2,2,LS-1)*CUSUX(K,2)
            CDSDX(K,1)=CS(1,1)*CUSUX(K,1)+CS(1,2)*CUSUX(K,2)
            CDSDX(K,2)=CS(2,1)*CUSUX(K,1)+CS(2,2)*CUSUX(K,2)
            CUSDX(K,1)=GRD(1,1,LS)*CDSDX(K,1)
     *                  +GRD(1,2,LS)*CDSDX(K,2)
            CUSDX(K,2)=GRD(2,1,LS)*CDSDX(K,1)
     *                  +GRD(2,2,LS)*CDSDX(K,2)
C
            CUSUZ(K,1)=E2(1,2)*DQZ
            CUSUZ(K,2)=E2(2,2)*DQZ
            CDSUZ(K,1)=GRU(1,1,LS-1)*CUSUZ(K,1)
     *                  +GRU(1,2,LS-1)*CUSUZ(K,2)
            CDSUZ(K,2)=GRU(2,1,LS-1)*CUSUZ(K,1)
     *                  +GRU(2,2,LS-1)*CUSUZ(K,2)
            CDSDZ(K,1)=CS(1,1)*CUSUZ(K,1)+CS(1,2)*CUSUZ(K,2)
            CDSDZ(K,2)=CS(2,1)*CUSUZ(K,1)+CS(2,2)*CUSUZ(K,2)
            CUSDZ(K,1)=GRD(1,1,LS)*CDSDZ(K,1)
     *                  +GRD(1,2,LS)*CDSDZ(K,2)
            CUSDZ(K,2)=GRD(2,1,LS)*CDSDZ(K,1)
     *                  +GRD(2,2,LS)*CDSDZ(K,2)
C
          ENDIF
C
        END IF
C
        IF(LS.EQ.1)GOTO 10
          CUX(K,1,LS-1)=GTU(1,1,LS-1)*CUSUX(K,1)
     *                   +GTU(1,2,LS-1)*CUSUX(K,2)
          CUX(K,2,LS-1)=GTU(2,1,LS-1)*CUSUX(K,1)
     *                   +GTU(2,2,LS-1)*CUSUX(K,2)
          CUZ(K,1,LS-1)=GTU(1,1,LS-1)*CUSUZ(K,1)
     *                   +GTU(1,2,LS-1)*CUSUZ(K,2)
          CUZ(K,2,LS-1)=GTU(2,1,LS-1)*CUSUZ(K,1)
     *                   +GTU(2,2,LS-1)*CUSUZ(K,2)
        DO 110 I=LS-2,1,-1
          CUX(K,1,I)  =GTU(1,1,I)*CUX(K,1,I+1)
     *                  +GTU(1,2,I)*CUX(K,2,I+1)
          CUX(K,2,I)  =GTU(2,1,I)*CUX(K,1,I+1)
     *                  +GTU(2,2,I)*CUX(K,2,I+1)
          CDX(K,1,I+1)=GRU(1,1,I)*CUX(K,1,I+1)
     *                  +GRU(1,2,I)*CUX(K,2,I+1)
          CDX(K,2,I+1)=GRU(2,1,I)*CUX(K,1,I+1)
     *                  +GRU(2,2,I)*CUX(K,2,I+1)
          CUZ(K,1,I)  =GTU(1,1,I)*CUZ(K,1,I+1)
     *                  +GTU(1,2,I)*CUZ(K,2,I+1)
          CUZ(K,2,I)  =GTU(2,1,I)*CUZ(K,1,I+1)
     *                  +GTU(2,2,I)*CUZ(K,2,I+1)
          CDZ(K,1,I+1)=GRU(1,1,I)*CUZ(K,1,I+1)
     *                  +GRU(1,2,I)*CUZ(K,2,I+1)
          CDZ(K,2,I+1)=GRU(2,1,I)*CUZ(K,1,I+1)
     *                  +GRU(2,2,I)*CUZ(K,2,I+1)
  110   CONTINUE
        CDX(K,1,1)=GRU(1,1,0)*CUX(K,1,1)+GRU(1,2,0)*CUX(K,2,1)
        CDX(K,2,1)=GRU(2,1,0)*CUX(K,1,1)+GRU(2,2,0)*CUX(K,2,1)
        CDZ(K,1,1)=GRU(1,1,0)*CUZ(K,1,1)+GRU(1,2,0)*CUZ(K,2,1)
        CDZ(K,2,1)=GRU(2,1,0)*CUZ(K,1,1)+GRU(2,2,0)*CUZ(K,2,1)
C
   10   IF(LS.EQ.NL)GOTO 100
        CDX(K,1,LS+1)=GTD(1,1,LS)*CDSDX(K,1)
     *                 +GTD(1,2,LS)*CDSDX(K,2)
        CDX(K,2,LS+1)=GTD(2,1,LS)*CDSDX(K,1)
     *                 +GTD(2,2,LS)*CDSDX(K,2)
        CDZ(K,1,LS+1)=GTD(1,1,LS)*CDSDZ(K,1)
     *                 +GTD(1,2,LS)*CDSDZ(K,2)
        CDZ(K,2,LS+1)=GTD(2,1,LS)*CDSDZ(K,1)
     *                 +GTD(2,2,LS)*CDSDZ(K,2)
        DO 120 I=LS+1,NL-1
          CDX(K,1,I+1)=GTD(1,1,I)*CDX(K,1,I)
     *                  +GTD(1,2,I)*CDX(K,2,I)
          CDX(K,2,I+1)=GTD(2,1,I)*CDX(K,1,I)
     *                  +GTD(2,2,I)*CDX(K,2,I)
          CUX(K,1,I)  =GRD(1,1,I)*CDX(K,1,I)
     *                  +GRD(1,2,I)*CDX(K,2,I)
          CUX(K,2,I)  =GRD(2,1,I)*CDX(K,1,I)
     *                  +GRD(2,2,I)*CDX(K,2,I)
          CDZ(K,1,I+1)=GTD(1,1,I)*CDZ(K,1,I)
     *                  +GTD(1,2,I)*CDZ(K,2,I)
          CDZ(K,2,I+1)=GTD(2,1,I)*CDZ(K,1,I)
     *                  +GTD(2,2,I)*CDZ(K,2,I)
          CUZ(K,1,I)  =GRD(1,1,I)*CDZ(K,1,I)
     *                  +GRD(1,2,I)*CDZ(K,2,I)
          CUZ(K,2,I)  =GRD(2,1,I)*CDZ(K,1,I)
     *                  +GRD(2,2,I)*CDZ(K,2,I)
  120   CONTINUE
C
  100 CONTINUE
C
      RETURN
      END
C
c
C *** SH DISPLACEMENT USING R/T MATRIX ***
C
      SUBROUTINE RTSH0(OM,K,WN,VLX)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VSc,CD,CU,VS2,WS2,CMS,SM,EE(2,2)
     *          ,CUSU,CDSU,CDSD,CUSD,DD,DU,VLX(2),OM,OM2
      COMMON /SOLY/SDEP,ODEP,LS,LO
     *       /STRC/DNS(20),VP(20),QP(20),FP(20)
     *                    ,VS(20),QS(20),FS(20),THK(20),NL
     *       /DUSH0/CD(2001,20),CU(2001,20)
     *       ,CUSU(2001),CDSU(2001),CDSD(2001),CUSD(2001)
C
C *** CHECK FOR THE DIMENSION OF CD AND CU ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(10,*)' FROM RTSH: NL.GT.NC ',NL,NC
         STOP
      END IF
C
C *** COMPUTATION OF VLX ***
C
      OM2=OM*OM
      WN2=WN*WN
C
      VSc=VS(LO)*DCMPLX(1.D0,-0.5D0/QS(LO))
      VS2=VSc*VSc
      WS2=OM2/VS2
      CMS=CDSQRT(WN2-WS2)
      IF(dREAL(CMS).LT.0.0)CMS=-CMS
C
      SM=DNS(LO)*VS2
      EE(1,1)=(1.0,0.0)
      EE(1,2)=(1.0,0.0)
      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(K,LO)
          VLX(2)=EE(2,1)*DD*CD(K,LO)
        ELSE
          IF(ODEP.GE.SDEP)THEN
            DD=CDEXP(-CMS*(ODEP-SDEP))
            VLX(1)=EE(1,1)*DD*CDSD(K)
            VLX(2)=EE(2,1)*DD*CDSD(K)
          ELSE
            DD=CDEXP(-CMS*ODEP)
            DU=CDEXP(-CMS*(SDEP-ODEP))
            VLX(1)=EE(1,1)*DD*CDSU(K)+EE(1,2)*DU*CUSU(K)
            VLX(2)=EE(2,1)*DD*CDSU(K)+EE(2,2)*DU*CUSU(K)
          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(K,LO)+EE(1,2)*DU*CU(K,LO)
          VLX(2)=EE(2,1)*DD*CD(K,LO)+EE(2,2)*DU*CU(K,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(K)+EE(1,2)*DU*CUSD(K)
            VLX(2)=EE(2,1)*DD*CDSD(K)+EE(2,2)*DU*CUSD(K)
          ELSE
            DD=CDEXP(-CMS*ODEP)
            DU=CDEXP(-CMS*(SDEP-ODEP))
            VLX(1)=EE(1,1)*DD*CDSU(K)+EE(1,2)*DU*CUSU(K)
            VLX(2)=EE(2,1)*DD*CDSU(K)+EE(2,2)*DU*CUSU(K)
          END IF
        END IF
C
      END IF
C
      RETURN
      END
C
c
c
C *** SH DISPLACEMENT USING R/T MATRIX ***
C
      SUBROUTINE RTSH(OM,K,WN,VLX)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VSc,CD,CU,VS2,WS2,CMS,SM,EE(2,2)
     *          ,CUSU,CDSU,CDSD,CUSD,DD,DU,VLX(2),OM,OM2
      COMMON /SOLY/SDEP,ODEP,LS,LO
     *       /STRC/DNS(20),VP(20),QP(20),FP(20)
     *                    ,VS(20),QS(20),FS(20),THK(20),NL
     *       /DUSH/CD(2001,20),CU(2001,20)
     *       ,CUSU(2001),CDSU(2001),CDSD(2001),CUSD(2001)
C
C *** CHECK FOR THE DIMENSION OF CD AND CU ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(10,*)' FROM RTSH: NL.GT.NC ',NL,NC
         STOP
      END IF
C
C *** COMPUTATION OF VLX ***
C
      OM2=OM*OM
      WN2=WN*WN
C
      freq=om/6.2831853072d0
      QSe=QS(LO)*(freq)**FS(LO)
      VSc=VS(LO)*DCMPLX(1.D0,-0.5D0/QSe)
      VS2=VSc*VSc
      WS2=OM2/VS2
      CMS=CDSQRT(WN2-WS2)
      IF(dREAL(CMS).LT.0.0)CMS=-CMS
C
      SM=DNS(LO)*VS2
      EE(1,1)=(1.0,0.0)
      EE(1,2)=(1.0,0.0)
      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(K,LO)
          VLX(2)=EE(2,1)*DD*CD(K,LO)
        ELSE
          IF(ODEP.GE.SDEP)THEN
            DD=CDEXP(-CMS*(ODEP-SDEP))
            VLX(1)=EE(1,1)*DD*CDSD(K)
            VLX(2)=EE(2,1)*DD*CDSD(K)
          ELSE
            DD=CDEXP(-CMS*ODEP)
            DU=CDEXP(-CMS*(SDEP-ODEP))
            VLX(1)=EE(1,1)*DD*CDSU(K)+EE(1,2)*DU*CUSU(K)
            VLX(2)=EE(2,1)*DD*CDSU(K)+EE(2,2)*DU*CUSU(K)
          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(K,LO)+EE(1,2)*DU*CU(K,LO)
          VLX(2)=EE(2,1)*DD*CD(K,LO)+EE(2,2)*DU*CU(K,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(K)+EE(1,2)*DU*CUSD(K)
            VLX(2)=EE(2,1)*DD*CDSD(K)+EE(2,2)*DU*CUSD(K)
          ELSE
            DD=CDEXP(-CMS*ODEP)
            DU=CDEXP(-CMS*(SDEP-ODEP))
            VLX(1)=EE(1,1)*DD*CDSU(K)+EE(1,2)*DU*CUSU(K)
            VLX(2)=EE(2,1)*DD*CDSU(K)+EE(2,2)*DU*CUSU(K)
          END IF
        END IF
C
      END IF
C
      RETURN
      END
C
C
C *** PSV DISPLACEMENT USING R/T MATRIX ***
C
      SUBROUTINE RTPSV0(OM,K,WN,URX,URZ)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 CMP,CMS,EE(4,4),DD(2,2),DU(2,2),OM
     *    ,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),QP(20),FP(20)
     *                 ,VS(20),QS(20),FS(20),THK(20),NL
     * /DUPSVX0/CDX(2001,2,20),CUX(2001,2,20),CUSUX(2001,2)
     *        ,CDSUX(2001,2),CDSDX(2001,2),CUSDX(2001,2) 
     * /DUPSVZ0/CDZ(2001,2,20),CUZ(2001,2,20),CUSUZ(2001,2)
     *        ,CDSUZ(2001,2),CDSDZ(2001,2),CUSDZ(2001,2) 
C
C *** CHECK FOR THE DIMENSION OF CD AND CU ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(10,*)' FROM RTPSV: NL.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
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(K,1,LO)+E1(1,2)*CDX(K,2,LO)
          URX(2,1)=E1(2,1)*CDX(K,1,LO)+E1(2,2)*CDX(K,2,LO)
          URZ(1,1)=E1(1,1)*CDZ(K,1,LO)+E1(1,2)*CDZ(K,2,LO)
          URZ(2,1)=E1(2,1)*CDZ(K,1,LO)+E1(2,2)*CDZ(K,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(K,1,LO)+E1(1,2)*CDX(K,2,LO)
          URX(2,2)=E1(2,1)*CDX(K,1,LO)+E1(2,2)*CDX(K,2,LO)
          URZ(1,2)=E1(1,1)*CDZ(K,1,LO)+E1(1,2)*CDZ(K,2,LO)
          URZ(2,2)=E1(2,1)*CDZ(K,1,LO)+E1(2,2)*CDZ(K,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(K,1)+E1(1,2)*CDSDX(K,2)
            URX(2,1)=E1(2,1)*CDSDX(K,1)+E1(2,2)*CDSDX(K,2)
            URZ(1,1)=E1(1,1)*CDSDZ(K,1)+E1(1,2)*CDSDZ(K,2)
            URZ(2,1)=E1(2,1)*CDSDZ(K,1)+E1(2,2)*CDSDZ(K,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(K,1)+E1(1,2)*CDSDX(K,2)
            URX(2,2)=E1(2,1)*CDSDX(K,1)+E1(2,2)*CDSDX(K,2)
            URZ(1,2)=E1(1,1)*CDSDZ(K,1)+E1(1,2)*CDSDZ(K,2)
            URZ(2,2)=E1(2,1)*CDSDZ(K,1)+E1(2,2)*CDSDZ(K,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(K,1)+E1(1,2)*CDSUX(K,2)
     *              +E2(1,1)*CUSUX(K,1)+E2(1,2)*CUSUX(K,2)
            URX(2,1)=E1(2,1)*CDSUX(K,1)+E1(2,2)*CDSUX(K,2)
     *              +E2(2,1)*CUSUX(K,1)+E2(2,2)*CUSUX(K,2)
            URZ(1,1)=E1(1,1)*CDSUZ(K,1)+E1(1,2)*CDSUZ(K,2)
     *              +E2(1,1)*CUSUZ(K,1)+E2(1,2)*CUSUZ(K,2)
            URZ(2,1)=E1(2,1)*CDSUZ(K,1)+E1(2,2)*CDSUZ(K,2)
     *              +E2(2,1)*CUSUZ(K,1)+E2(2,2)*CUSUZ(K,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(K,1)+E1(1,2)*CDSUX(K,2)
     *              +E2(1,1)*CUSUX(K,1)+E2(1,2)*CUSUX(K,2)
            URX(2,2)=E1(2,1)*CDSUX(K,1)+E1(2,2)*CDSUX(K,2)
     *              +E2(2,1)*CUSUX(K,1)+E2(2,2)*CUSUX(K,2)
            URZ(1,2)=E1(1,1)*CDSUZ(K,1)+E1(1,2)*CDSUZ(K,2)
     *              +E2(1,1)*CUSUZ(K,1)+E2(1,2)*CUSUZ(K,2)
            URZ(2,2)=E1(2,1)*CDSUZ(K,1)+E1(2,2)*CDSUZ(K,2)
     *              +E2(2,1)*CUSUZ(K,1)+E2(2,2)*CUSUZ(K,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(K,1,LO)+E1(1,2)*CDX(K,2,LO)
     *              +E2(1,1)*CUX(K,1,LO)+E2(1,2)*CUX(K,2,LO)
            URX(2,1)=E1(2,1)*CDX(K,1,LO)+E1(2,2)*CDX(K,2,LO)
     *              +E2(2,1)*CUX(K,1,LO)+E2(2,2)*CUX(K,2,LO)
            URZ(1,1)=E1(1,1)*CDZ(K,1,LO)+E1(1,2)*CDZ(K,2,LO)
     *              +E2(1,1)*CUZ(K,1,LO)+E2(1,2)*CUZ(K,2,LO)
            URZ(2,1)=E1(2,1)*CDZ(K,1,LO)+E1(2,2)*CDZ(K,2,LO)
     *              +E2(2,1)*CUZ(K,1,LO)+E2(2,2)*CUZ(K,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(K,1,LO)+E1(1,2)*CDX(K,2,LO)
     *              +E2(1,1)*CUX(K,1,LO)+E2(1,2)*CUX(K,2,LO)
            URX(2,2)=E1(2,1)*CDX(K,1,LO)+E1(2,2)*CDX(K,2,LO)
     *              +E2(2,1)*CUX(K,1,LO)+E2(2,2)*CUX(K,2,LO)
            URZ(1,2)=E1(1,1)*CDZ(K,1,LO)+E1(1,2)*CDZ(K,2,LO)
     *              +E2(1,1)*CUZ(K,1,LO)+E2(1,2)*CUZ(K,2,LO)
            URZ(2,2)=E1(2,1)*CDZ(K,1,LO)+E1(2,2)*CDZ(K,2,LO)
     *              +E2(2,1)*CUZ(K,1,LO)+E2(2,2)*CUZ(K,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(K,1)+E1(1,2)*CDSDX(K,2)
     *              +E2(1,1)*CUSDX(K,1)+E2(1,2)*CUSDX(K,2)
            URX(2,1)=E1(2,1)*CDSDX(K,1)+E1(2,2)*CDSDX(K,2)
     *              +E2(2,1)*CUSDX(K,1)+E2(2,2)*CUSDX(K,2)
            URZ(1,1)=E1(1,1)*CDSDZ(K,1)+E1(1,2)*CDSDZ(K,2)
     *              +E2(1,1)*CUSDZ(K,1)+E2(1,2)*CUSDZ(K,2)
            URZ(2,1)=E1(2,1)*CDSDZ(K,1)+E1(2,2)*CDSDZ(K,2)
     *              +E2(2,1)*CUSDZ(K,1)+E2(2,2)*CUSDZ(K,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(K,1)+E1(1,2)*CDSDX(K,2)
     *              +E2(1,1)*CUSDX(K,1)+E2(1,2)*CUSDX(K,2)
            URX(2,2)=E1(2,1)*CDSDX(K,1)+E1(2,2)*CDSDX(K,2)
     *              +E2(2,1)*CUSDX(K,1)+E2(2,2)*CUSDX(K,2)
            URZ(1,2)=E1(1,1)*CDSDZ(K,1)+E1(1,2)*CDSDZ(K,2)
     *              +E2(1,1)*CUSDZ(K,1)+E2(1,2)*CUSDZ(K,2)
            URZ(2,2)=E1(2,1)*CDSDZ(K,1)+E1(2,2)*CDSDZ(K,2)
     *              +E2(2,1)*CUSDZ(K,1)+E2(2,2)*CUSDZ(K,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(K,1)+E1(1,2)*CDSUX(K,2)
     *              +E2(1,1)*CUSUX(K,1)+E2(1,2)*CUSUX(K,2)
            URX(2,1)=E1(2,1)*CDSUX(K,1)+E1(2,2)*CDSUX(K,2)
     *              +E2(2,1)*CUSUX(K,1)+E2(2,2)*CUSUX(K,2)
            URZ(1,1)=E1(1,1)*CDSUZ(K,1)+E1(1,2)*CDSUZ(K,2)
     *              +E2(1,1)*CUSUZ(K,1)+E2(1,2)*CUSUZ(K,2)
            URZ(2,1)=E1(2,1)*CDSUZ(K,1)+E1(2,2)*CDSUZ(K,2)
     *              +E2(2,1)*CUSUZ(K,1)+E2(2,2)*CUSUZ(K,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(K,1)+E1(1,2)*CDSUX(K,2)
     *              +E2(1,1)*CUSUX(K,1)+E2(1,2)*CUSUX(K,2)
            URX(2,2)=E1(2,1)*CDSUX(K,1)+E1(2,2)*CDSUX(K,2)
     *              +E2(2,1)*CUSUX(K,1)+E2(2,2)*CUSUX(K,2)
            URZ(1,2)=E1(1,1)*CDSUZ(K,1)+E1(1,2)*CDSUZ(K,2)
     *              +E2(1,1)*CUSUZ(K,1)+E2(1,2)*CUSUZ(K,2)
            URZ(2,2)=E1(2,1)*CDSUZ(K,1)+E1(2,2)*CDSUZ(K,2)
     *              +E2(2,1)*CUSUZ(K,1)+E2(2,2)*CUSUZ(K,2)
C
          END IF
C
        END IF
C
      END IF
C
      RETURN
      END
C
C
C *** PSV DISPLACEMENT USING R/T MATRIX ***
C
      SUBROUTINE RTPSV(OM,K,WN,URX,URZ)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 CMP,CMS,EE(4,4),DD(2,2),DU(2,2),OM
     *    ,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),QP(20),FP(20)
     *                 ,VS(20),QS(20),FS(20),THK(20),NL
     * /DUPSVX/CDX(2001,2,20),CUX(2001,2,20),CUSUX(2001,2)
     *        ,CDSUX(2001,2),CDSDX(2001,2),CUSDX(2001,2) 
     * /DUPSVZ/CDZ(2001,2,20),CUZ(2001,2,20),CUSUZ(2001,2)
     *        ,CDSUZ(2001,2),CDSDZ(2001,2),CUSDZ(2001,2) 
C
C *** CHECK FOR THE DIMENSION OF CD AND CU ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(10,*)' FROM RTPSV: NL.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
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(K,1,LO)+E1(1,2)*CDX(K,2,LO)
          URX(2,1)=E1(2,1)*CDX(K,1,LO)+E1(2,2)*CDX(K,2,LO)
          URZ(1,1)=E1(1,1)*CDZ(K,1,LO)+E1(1,2)*CDZ(K,2,LO)
          URZ(2,1)=E1(2,1)*CDZ(K,1,LO)+E1(2,2)*CDZ(K,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(K,1,LO)+E1(1,2)*CDX(K,2,LO)
          URX(2,2)=E1(2,1)*CDX(K,1,LO)+E1(2,2)*CDX(K,2,LO)
          URZ(1,2)=E1(1,1)*CDZ(K,1,LO)+E1(1,2)*CDZ(K,2,LO)
          URZ(2,2)=E1(2,1)*CDZ(K,1,LO)+E1(2,2)*CDZ(K,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(K,1)+E1(1,2)*CDSDX(K,2)
            URX(2,1)=E1(2,1)*CDSDX(K,1)+E1(2,2)*CDSDX(K,2)
            URZ(1,1)=E1(1,1)*CDSDZ(K,1)+E1(1,2)*CDSDZ(K,2)
            URZ(2,1)=E1(2,1)*CDSDZ(K,1)+E1(2,2)*CDSDZ(K,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(K,1)+E1(1,2)*CDSDX(K,2)
            URX(2,2)=E1(2,1)*CDSDX(K,1)+E1(2,2)*CDSDX(K,2)
            URZ(1,2)=E1(1,1)*CDSDZ(K,1)+E1(1,2)*CDSDZ(K,2)
            URZ(2,2)=E1(2,1)*CDSDZ(K,1)+E1(2,2)*CDSDZ(K,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(K,1)+E1(1,2)*CDSUX(K,2)
     *              +E2(1,1)*CUSUX(K,1)+E2(1,2)*CUSUX(K,2)
            URX(2,1)=E1(2,1)*CDSUX(K,1)+E1(2,2)*CDSUX(K,2)
     *              +E2(2,1)*CUSUX(K,1)+E2(2,2)*CUSUX(K,2)
            URZ(1,1)=E1(1,1)*CDSUZ(K,1)+E1(1,2)*CDSUZ(K,2)
     *              +E2(1,1)*CUSUZ(K,1)+E2(1,2)*CUSUZ(K,2)
            URZ(2,1)=E1(2,1)*CDSUZ(K,1)+E1(2,2)*CDSUZ(K,2)
     *              +E2(2,1)*CUSUZ(K,1)+E2(2,2)*CUSUZ(K,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(K,1)+E1(1,2)*CDSUX(K,2)
     *              +E2(1,1)*CUSUX(K,1)+E2(1,2)*CUSUX(K,2)
            URX(2,2)=E1(2,1)*CDSUX(K,1)+E1(2,2)*CDSUX(K,2)
     *              +E2(2,1)*CUSUX(K,1)+E2(2,2)*CUSUX(K,2)
            URZ(1,2)=E1(1,1)*CDSUZ(K,1)+E1(1,2)*CDSUZ(K,2)
     *              +E2(1,1)*CUSUZ(K,1)+E2(1,2)*CUSUZ(K,2)
            URZ(2,2)=E1(2,1)*CDSUZ(K,1)+E1(2,2)*CDSUZ(K,2)
     *              +E2(2,1)*CUSUZ(K,1)+E2(2,2)*CUSUZ(K,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(K,1,LO)+E1(1,2)*CDX(K,2,LO)
     *              +E2(1,1)*CUX(K,1,LO)+E2(1,2)*CUX(K,2,LO)
            URX(2,1)=E1(2,1)*CDX(K,1,LO)+E1(2,2)*CDX(K,2,LO)
     *              +E2(2,1)*CUX(K,1,LO)+E2(2,2)*CUX(K,2,LO)
            URZ(1,1)=E1(1,1)*CDZ(K,1,LO)+E1(1,2)*CDZ(K,2,LO)
     *              +E2(1,1)*CUZ(K,1,LO)+E2(1,2)*CUZ(K,2,LO)
            URZ(2,1)=E1(2,1)*CDZ(K,1,LO)+E1(2,2)*CDZ(K,2,LO)
     *              +E2(2,1)*CUZ(K,1,LO)+E2(2,2)*CUZ(K,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(K,1,LO)+E1(1,2)*CDX(K,2,LO)
     *              +E2(1,1)*CUX(K,1,LO)+E2(1,2)*CUX(K,2,LO)
            URX(2,2)=E1(2,1)*CDX(K,1,LO)+E1(2,2)*CDX(K,2,LO)
     *              +E2(2,1)*CUX(K,1,LO)+E2(2,2)*CUX(K,2,LO)
            URZ(1,2)=E1(1,1)*CDZ(K,1,LO)+E1(1,2)*CDZ(K,2,LO)
     *              +E2(1,1)*CUZ(K,1,LO)+E2(1,2)*CUZ(K,2,LO)
            URZ(2,2)=E1(2,1)*CDZ(K,1,LO)+E1(2,2)*CDZ(K,2,LO)
     *              +E2(2,1)*CUZ(K,1,LO)+E2(2,2)*CUZ(K,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(K,1)+E1(1,2)*CDSDX(K,2)
     *              +E2(1,1)*CUSDX(K,1)+E2(1,2)*CUSDX(K,2)
            URX(2,1)=E1(2,1)*CDSDX(K,1)+E1(2,2)*CDSDX(K,2)
     *              +E2(2,1)*CUSDX(K,1)+E2(2,2)*CUSDX(K,2)
            URZ(1,1)=E1(1,1)*CDSDZ(K,1)+E1(1,2)*CDSDZ(K,2)
     *              +E2(1,1)*CUSDZ(K,1)+E2(1,2)*CUSDZ(K,2)
            URZ(2,1)=E1(2,1)*CDSDZ(K,1)+E1(2,2)*CDSDZ(K,2)
     *              +E2(2,1)*CUSDZ(K,1)+E2(2,2)*CUSDZ(K,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(K,1)+E1(1,2)*CDSDX(K,2)
     *              +E2(1,1)*CUSDX(K,1)+E2(1,2)*CUSDX(K,2)
            URX(2,2)=E1(2,1)*CDSDX(K,1)+E1(2,2)*CDSDX(K,2)
     *              +E2(2,1)*CUSDX(K,1)+E2(2,2)*CUSDX(K,2)
            URZ(1,2)=E1(1,1)*CDSDZ(K,1)+E1(1,2)*CDSDZ(K,2)
     *              +E2(1,1)*CUSDZ(K,1)+E2(1,2)*CUSDZ(K,2)
            URZ(2,2)=E1(2,1)*CDSDZ(K,1)+E1(2,2)*CDSDZ(K,2)
     *              +E2(2,1)*CUSDZ(K,1)+E2(2,2)*CUSDZ(K,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(K,1)+E1(1,2)*CDSUX(K,2)
     *              +E2(1,1)*CUSUX(K,1)+E2(1,2)*CUSUX(K,2)
            URX(2,1)=E1(2,1)*CDSUX(K,1)+E1(2,2)*CDSUX(K,2)
     *              +E2(2,1)*CUSUX(K,1)+E2(2,2)*CUSUX(K,2)
            URZ(1,1)=E1(1,1)*CDSUZ(K,1)+E1(1,2)*CDSUZ(K,2)
     *              +E2(1,1)*CUSUZ(K,1)+E2(1,2)*CUSUZ(K,2)
            URZ(2,1)=E1(2,1)*CDSUZ(K,1)+E1(2,2)*CDSUZ(K,2)
     *              +E2(2,1)*CUSUZ(K,1)+E2(2,2)*CUSUZ(K,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(K,1)+E1(1,2)*CDSUX(K,2)
     *              +E2(1,1)*CUSUX(K,1)+E2(1,2)*CUSUX(K,2)
            URX(2,2)=E1(2,1)*CDSUX(K,1)+E1(2,2)*CDSUX(K,2)
     *              +E2(2,1)*CUSUX(K,1)+E2(2,2)*CUSUX(K,2)
            URZ(1,2)=E1(1,1)*CDSUZ(K,1)+E1(1,2)*CDSUZ(K,2)
     *              +E2(1,1)*CUSUZ(K,1)+E2(1,2)*CUSUZ(K,2)
            URZ(2,2)=E1(2,1)*CDSUZ(K,1)+E1(2,2)*CDSUZ(K,2)
     *              +E2(2,1)*CUSUZ(K,1)+E2(2,2)*CUSUZ(K,2)
C
          END IF
C
        END IF
C
      END IF
C
      RETURN
      END
c
C
C *** SH GENERALIZED R/T MATIX ***
C
      SUBROUTINE GNRTSH(OM,WN)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 GTD,GTU,GRD,GRU,TD,TU,RD,RU,OM
      COMMON /SOLY/SDEP,ODEP,LS,LO
     *    /STRC/DNS(20),VP(20),QP(20),FP(20)
     *                 ,VS(20),QS(20),FS(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(10,*)' 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 VSc,VS2,WS2,CMS,SM,EA(2,2),EB(2,2)
     *      ,TD,TU,RD,RU,DDA,DUB,EAB,OM,OM2
      COMMON /SOLY/SDEP,ODEP,LS,LO
     *    /STRC/DNS(20),VP(20),QP(20),FP(20)
     *                 ,VS(20),QS(20),FS(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(10,*)' FROM MDRTSH: NL.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C *** MODIFIED R/T MATRIX ***
C
      WN2=WN*WN
      OM2=OM*OM
C
      if(dreal(om).eq.0.d0)then
        VSc=VS(1)*DCMPLX(1.D0,-0.5D0/QS(1))
	else
        freq=om/6.2831853072d0
        QSe=QS(1)*(freq)**FS(1)
        VSc=VS(1)*DCMPLX(1.D0,-0.5D0/QSe)
      end if
      VS2=VSc*VSc
      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
        if(dreal(om).eq.0.d0)then
          VSc=VS(I+1)*DCMPLX(1.D0,-0.5D0/QS(I+1))
	  else
          QSe=QS(I+1)*(freq)**FS(I+1)
          VSc=VS(I+1)*DCMPLX(1.D0,-0.5D0/QSe)
        end if
        VS2=VSc*VSc
        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        
      if(dreal(om).eq.0.d0)then
        VSc=VS(NL)*DCMPLX(1.D0,-0.5D0/QS(NL))
	else
        QSe=QS(NL)*(freq)**FS(NL)
        VSc=VS(NL)*DCMPLX(1.D0,-0.5D0/QSe)
      end if
      VS2=VSc*VSc
      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 *** P-SV GENERALIZED R/T MATIX ***
C
      SUBROUTINE GNRTPSV(OM,WN)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 GTD,GTU,GRD,GRU,TD,TU,RD,RU,E1(2,2),E2(2,2),DT,OM
      COMMON /SOLY/SDEP,ODEP,LS,LO
     *    /STRC/DNS(20),VP(20),QP(20),FP(20)
     *                 ,VS(20),QS(20),FS(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(10,*)' 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,OM
      COMMON /SOLY/SDEP,ODEP,LS,LO
     *    /STRC/DNS(20),VP(20),QP(20),FP(20)
     *                 ,VS(20),QS(20),FS(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(10,*)' 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
C *** P-SV ELEMENT MATIX ***
C
      SUBROUTINE MTXPSV(OM,WN,L,CMS,CMP,EE)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VPc,VSc,VS2,WS2,CMS,VP2,WP2,CMP
     *      ,SM,EE(4,4),X0,AA,CMSG,CMPG,OM,OM2
      COMMON/STRC/DNS(20),VP(20),QP(20),FP(20)
     *                   ,VS(20),QS(20),FS(20),THK(20),NL
C
C *** CHECK FOR THE DIMENSION OF COMMOM ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(10,*)' FROM MTXPSV: NL.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C *** ELEMENT MATRIX ***
C
      if(dreal(om).eq.0.d0)then
        VPc=VP(L)*DCMPLX(1.D0,-0.5D0/QP(L))
        VSc=VS(L)*DCMPLX(1.D0,-0.5D0/QS(L))
	else
        freq=om/6.2831853072d0
        QPe=QP(L)*(freq)**FP(L)
        QSe=QS(L)*(freq)**FS(L)
        VPc=VP(L)*DCMPLX(1.D0,-0.5D0/QPe)
        VSc=VS(L)*DCMPLX(1.D0,-0.5D0/QSe)
      end if
      VP2=VPc*VPc
      VS2=VSc*VSc
      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,OM
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)) /DSQRT(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
      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)) /DSQRT(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
