C
C *** THE COMPUTATION OF Q, V, AND W AND THEIR DERIVATIVES 
C     USING Newton-Cotes Quadrature (for Static Range 1)***
C
      SUBROUTINE QVWRZS1(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)
     *          ,QX1(0:2),QZ1(0:2),VX1(0:2),WX1(0:2),WZ1(0:2)
     *          ,QX2(0:2),QZ2(0:2),VX2(0:2),WX2(0:2),WZ2(0:2)
     *          ,QX3(0:2),QZ3(0:2),VX3(0:2),WX3(0:2),WZ3(0:2)
     *          ,QX4(0:2),QZ4(0:2),VX4(0:2),WX4(0:2),WZ4(0:2)
     *          ,QX5(0:2),QZ5(0:2),VX5(0:2),WX5(0:2),WZ5(0:2)
     *          ,QX6(0:2),QZ6(0:2),VX6(0:2),WX6(0:2),WZ6(0:2)
     *          ,QX7(0:2),QZ7(0:2),VX7(0:2),WX7(0:2),WZ7(0:2)
     *          ,QX8(0:2),QZ8(0:2),VX8(0:2),WX8(0:2),WZ8(0:2)
     *          ,QX9(0:2),QZ9(0:2),VX9(0:2),WX9(0:2),WZ9(0:2)
     *          ,QX10(0:2),QZ10(0:2),VX10(0:2),WX10(0:2),WZ10(0:2)
     *          ,QX11(0:2),QZ11(0:2),VX11(0:2),WX11(0:2),WZ11(0:2)
     *          ,QX12(0:2),QZ12(0:2),VX12(0:2),WX12(0:2),WZ12(0:2)
     *          ,QX13(0:2),QZ13(0:2),VX13(0:2),WX13(0:2),WZ13(0:2)
     *          ,QX14(0:2),QZ14(0:2),VX14(0:2),WX14(0:2),WZ14(0:2)
     *          ,QX15(0:2),QZ15(0:2),VX15(0:2),WX15(0:2),WZ15(0:2)
     *          ,QX16(0:2),QZ16(0:2),VX16(0:2),WX16(0:2),WZ16(0:2)
     *          ,QX17(0:2),QZ17(0:2),VX17(0:2),WX17(0:2),WZ17(0:2)
	*          ,OM,SM,RAM,aa
      COMMON /WVNUMs/WNRs(3),NWN1o,NWN2o,NWN1,NWN2,NWT1,NWT2
     *    /OTST/IOT,IST  /FCQVW/SM,RAM,RR  /SOLY/SDEP,ODEP,LS,LO
C
C   *Newton-Cotes Quadrature (for Static Range 1) *
C
      ns=NWN1
      if(ns.eq.0)return
      w2=WNRs(2)
	w1=WNRs(1)
	iter=0
      sImax1=0.d0
      sImax2=0.d0
	nint=0
c      write(*,*)
c
   10 iter=iter+1
c      write(*,*)
c      write(*,*)'iter=',iter
      dw=(w2-w1)/ns
	H=dw/RR
      CALL INTEGS1(OM,W1,QX1,QZ1,VX1,WX1,WZ1)
	nint=nint+1
      IF(IOT.eq.1)THEN
        WRITE(20,1100)W1,QX1(0),VX1(0),WX1(0),QZ1(0),WZ1(0)
	*                  ,QX1(1),VX1(1),WX1(1),QZ1(1),WZ1(1)
	*                  ,QX1(2),VX1(2),WX1(2),QZ1(2),WZ1(2)
      END IF
 1100 FORMAT(31E12.4)
c
      DO  100 I=1,ns
c
c     * 3 Points Newton-Cotes (Simpson) Quadrature *
c     * 1-2-3 *
c
	  index=1
        W3=W1+dw*dfloat(I)
        W2=W3-dw*0.5d0
        CALL INTEGS1(OM,W2,QX2,QZ2,VX2,WX2,WZ2)
        CALL INTEGS1(OM,W3,QX3,QZ3,VX3,WX3,WZ3)
	  nint=nint+2
	  sI11=cdabs((QZ1(2)+4.0*QZ2(2)+QZ3(2)))*H/6.0
	  sI12=cdabs((VX1(2)+4.0*VX2(2)+VX3(2)))*H/6.0
	  sInt1=sI11
	  sInt2=sI12
	  if(sI11.ge.sImax1)sImax1=sI11
	  if(sI12.ge.sImax2)sImax2=sI12
	  if(((sI11/sImax1).lt.0.001d0).and.
     *                       ((sI12/sImax2).lt.0.001d0))then
c	    write(*,*)'i,index=',i,index
c	    write(*,*)'sI11,sI12',sI11,sI12
c	    write(*,*)'sImax1,sImax2',sImax1,sImax2
c	    write(*,*)'sI11/sImax1 < 0.001d !!!: sI11/sImax1=',sI11/sImax1
c	    write(*,*)'sI12/sImax2 < 0.001d !!!: sI12/sImax2=',sI12/sImax2
          DO J=0,2
            GQX(J)= GQX(J)+H*(QX1(J)+4.0*QX2(J)+QX3(J))/6.0
            GVX(J)= GVX(J)+H*(VX1(J)+4.0*VX2(J)+VX3(J))/6.0
            GWX(J)= GWX(J)+H*(WX1(J)+4.0*WX2(J)+WX3(J))/6.0
            GQZ(J)= GQZ(J)+H*(QZ1(J)+4.0*QZ2(J)+QZ3(J))/6.0
            GWZ(J)= GWZ(J)+H*(WZ1(J)+4.0*WZ2(J)+WZ3(J))/6.0
            QX1(J)= QX3(J)
            VX1(J)= VX3(J)
            WX1(J)= WX3(J)
            QZ1(J)= QZ3(J)
            WZ1(J)= WZ3(J)
          end do
          IF(IOT.eq.1)THEN
            WRITE(20,1100)W2,QX2(0),VX2(0),WX2(0),QZ2(0),WZ2(0)
	*                      ,QX2(1),VX2(1),WX2(1),QZ2(1),WZ2(1)
	*                      ,QX2(2),VX2(2),WX2(2),QZ2(2),WZ2(2)
            WRITE(20,1100)W3,QX3(0),VX3(0),WX3(0),QZ3(0),WZ3(0)
	*                      ,QX3(1),VX3(1),WX3(1),QZ3(1),WZ3(1)
	*                      ,QX3(2),VX3(2),WX3(2),QZ3(2),WZ3(2)
          END IF
	    goto 100
	  end if
c
c     * 5 Points Newton-Cotes Quadrature *
c     * 1-4-2-5-3 *
c
	  index=2
        W4=W3-dw*0.75d0
        CALL INTEGS1(OM,W4,QX4,QZ4,VX4,WX4,WZ4)
        W5=W3-dw*0.25d0
        CALL INTEGS1(OM,W5,QX5,QZ5,VX5,WX5,WZ5)
	  nint=nint+2
	  aa=7.0*QZ1(2)+32.0*QZ4(2)+12.0*QZ2(2)+32.0*QZ5(2)+7.0*QZ3(2)
	  sI21=cdabs(aa)*H/90.0
	  aa=7.0*VX1(2)+32.0*VX4(2)+12.0*VX2(2)+32.0*VX5(2)+7.0*VX3(2)
	  sI22=cdabs(aa)*H/90.0
	  sInt1=sI21
	  sInt2=sI22
	  if(sI21.ge.sImax1)sImax1=sI21
	  if(sI22.ge.sImax2)sImax2=sI22
	  er11=dabs((sI11-sI21)/sI21)
	  er12=dabs((sI12-sI22)/sI22)
	  if((er11.lt.0.0001d0).and.(er12.lt.0.0001d0))then
c	    write(*,*)'i,index=',i,index
c	    write(*,*)'sI11,sI12',sI11,sI12
c 	    write(*,*)'sI21,sI22',sI21,sI22
c	    write(*,*)'er11,er12',er11,er12
          do J=0,2
            aa= 7.0*QX1(J)+32.0*QX4(J)+12.0*QX2(J)
            GQX(J)= GQX(J)+H*(aa+32.0*QX5(J)+7.0*QX3(J))/90.0
            aa= 7.0*VX1(J)+32.0*VX4(J)+12.0*VX2(J)
            GVX(J)= GVX(J)+H*(aa+32.0*VX5(J)+7.0*VX3(J))/90.0
            aa= 7.0*WX1(J)+32.0*WX4(J)+12.0*WX2(J)
            GWX(J)= GWX(J)+H*(aa+32.0*WX5(J)+7.0*WX3(J))/90.0
            aa= 7.0*QZ1(J)+32.0*QZ4(J)+12.0*QZ2(J)
            GQZ(J)= GQZ(J)+H*(aa+32.0*QZ5(J)+7.0*QZ3(J))/90.0
            aa= 7.0*WZ1(J)+32.0*WZ4(J)+12.0*WZ2(J)
            GWZ(J)= GWZ(J)+H*(aa+32.0*WZ5(J)+7.0*WZ3(J))/90.0
            QX1(J)= QX3(J)
            VX1(J)= VX3(J)
            WX1(J)= WX3(J)
            QZ1(J)= QZ3(J)
            WZ1(J)= WZ3(J)
          end do
          IF(IOT.eq.1)THEN
            WRITE(20,1100)W4,QX4(0),VX4(0),WX4(0),QZ4(0),WZ4(0)
	*                      ,QX4(1),VX4(1),WX4(1),QZ4(1),WZ4(1)
	*                      ,QX4(2),VX4(2),WX4(2),QZ4(2),WZ4(2)
            WRITE(20,1100)W2,QX2(0),VX2(0),WX2(0),QZ2(0),WZ2(0)
	*                      ,QX2(1),VX2(1),WX2(1),QZ2(1),WZ2(1)
	*                      ,QX2(2),VX2(2),WX2(2),QZ2(2),WZ2(2)
            WRITE(20,1100)W5,QX5(0),VX5(0),WX5(0),QZ5(0),WZ5(0)
	*                      ,QX5(1),VX5(1),WX5(1),QZ5(1),WZ5(1)
	*                      ,QX5(2),VX5(2),WX5(2),QZ5(2),WZ5(2)
            WRITE(20,1100)W3,QX3(0),VX3(0),WX3(0),QZ3(0),WZ3(0)
	*                      ,QX3(1),VX3(1),WX3(1),QZ3(1),WZ3(1)
	*                      ,QX3(2),VX3(2),WX3(2),QZ3(2),WZ3(2)
          END IF
	    goto 100
	  end if
c
c     * 9 Points Newton-Cotes Quadrature *
c     * 1-6-4-7-2-8-5-9-3 *
c
	  index=3
        W6=W3-dw*0.875d0
        CALL INTEGS1(OM,W6,QX6,QZ6,VX6,WX6,WZ6)
        W7=W3-dw*0.625d0
        CALL INTEGS1(OM,W7,QX7,QZ7,VX7,WX7,WZ7)
        W8=W3-dw*0.375d0
        CALL INTEGS1(OM,W8,QX8,QZ8,VX8,WX8,WZ8)
        W9=W3-dw*0.125d0
        CALL INTEGS1(OM,W9,QX9,QZ9,VX9,WX9,WZ9)
	  nint=nint+4
	  aa=989.0*QZ1(2)+5888.0*QZ6(2)-928.0*QZ4(2)+10496.0*QZ7(2)-4540.0
	*   *QZ2(2)+10496.0*QZ8(2)-928.0*QZ5(2)+5888.0*QZ9(2)+989.0*QZ3(2)
	  sI31=cdabs(aa)*H/28350.0
	  aa=989.0*VX1(2)+5888.0*VX6(2)-928.0*VX4(2)+10496.0*VX7(2)-4540.0
	*   *VX2(2)+10496.0*VX8(2)-928.0*VX5(2)+5888.0*VX9(2)+989.0*VX3(2)
	  sI32=cdabs(aa)*H/28350.0
	  sInt1=sI31
	  sInt2=sI32
	  if(sI31.ge.sImax1)sImax1=sI31
	  if(sI32.ge.sImax2)sImax2=sI32
	  er21=dabs((sI21-sI31)/sI31)
	  er22=dabs((sI22-sI32)/sI32)
	  if((er21.lt.0.0001d0).and.(er22.lt.0.0001d0))then
c	    write(*,*)'i,index=',i,index
c	    write(*,*)'sI11,sI12',sI11,sI12
c 	    write(*,*)'sI21,sI22',sI21,sI22
c 	    write(*,*)'sI31,sI32',sI31,sI32
c	    write(*,*)'er21,er22',er21,er22
          do J=0,2
            aa= 989.0*QX1(J)+5888.0*QX6(J)-928.0*QX4(J)+10496.0*QX7(J)
            aa= aa-4540.0*QX2(J)+10496.0*QX8(J)-928.0*QX5(J)
            GQX(J)= GQX(J)+H*(aa+5888.0*QX9(J)+989.0*QX3(J))/28350.0
            aa= 989.0*VX1(J)+5888.0*VX6(J)-928.0*VX4(J)+10496.0*VX7(J)
            aa= aa-4540.0*VX2(J)+10496.0*VX8(J)-928.0*VX5(J)
            GVX(J)= GVX(J)+H*(aa+5888.0*VX9(J)+989.0*VX3(J))/28350.0
            aa= 989.0*WX1(J)+5888.0*WX6(J)-928.0*WX4(J)+10496.0*WX7(J)
            aa= aa-4540.0*WX2(J)+10496.0*WX8(J)-928.0*WX5(J)
            GWX(J)= GWX(J)+H*(aa+5888.0*WX9(J)+989.0*WX3(J))/28350.0
            aa= 989.0*QZ1(J)+5888.0*QZ6(J)-928.0*QZ4(J)+10496.0*QZ7(J)
            aa= aa-4540.0*QZ2(J)+10496.0*QZ8(J)-928.0*QZ5(J)
            GQZ(J)= GQZ(J)+H*(aa+5888.0*QZ9(J)+989.0*QZ3(J))/28350.0
            aa= 989.0*WZ1(J)+5888.0*WZ6(J)-928.0*WZ4(J)+10496.0*WZ7(J)
            aa= aa-4540.0*WZ2(J)+10496.0*WZ8(J)-928.0*WZ5(J)
            GWZ(J)= GWZ(J)+H*(aa+5888.0*WZ9(J)+989.0*WZ3(J))/28350.0
            QX1(J)= QX3(J)
            VX1(J)= VX3(J)
            WX1(J)= WX3(J)
            QZ1(J)= QZ3(J)
            WZ1(J)= WZ3(J)
          end do
          IF(IOT.eq.1)THEN
            WRITE(20,1100)W6,QX6(0),VX6(0),WX6(0),QZ6(0),WZ6(0)
	*                      ,QX6(1),VX6(1),WX6(1),QZ6(1),WZ6(1)
	*                      ,QX6(2),VX6(2),WX6(2),QZ6(2),WZ6(2)
            WRITE(20,1100)W4,QX4(0),VX4(0),WX4(0),QZ4(0),WZ4(0)
	*                      ,QX4(1),VX4(1),WX4(1),QZ4(1),WZ4(1)
	*                      ,QX4(2),VX4(2),WX4(2),QZ4(2),WZ4(2)
            WRITE(20,1100)W7,QX7(0),VX7(0),WX7(0),QZ7(0),WZ7(0)
	*                      ,QX7(1),VX7(1),WX7(1),QZ7(1),WZ7(1)
	*                      ,QX7(2),VX7(2),WX7(2),QZ7(2),WZ7(2)
            WRITE(20,1100)W2,QX2(0),VX2(0),WX2(0),QZ2(0),WZ2(0)
	*                      ,QX2(1),VX2(1),WX2(1),QZ2(1),WZ2(1)
	*                      ,QX2(2),VX2(2),WX2(2),QZ2(2),WZ2(2)
            WRITE(20,1100)W8,QX8(0),VX8(0),WX8(0),QZ8(0),WZ8(0)
	*                      ,QX8(1),VX8(1),WX8(1),QZ8(1),WZ8(1)
	*                      ,QX8(2),VX8(2),WX8(2),QZ8(2),WZ8(2)
            WRITE(20,1100)W5,QX5(0),VX5(0),WX5(0),QZ5(0),WZ5(0)
	*                      ,QX5(1),VX5(1),WX5(1),QZ5(1),WZ5(1)
	*                      ,QX5(2),VX5(2),WX5(2),QZ5(2),WZ5(2)
            WRITE(20,1100)W9,QX9(0),VX9(0),WX9(0),QZ9(0),WZ9(0)
	*                      ,QX9(1),VX9(1),WX9(1),QZ9(1),WZ9(1)
	*                      ,QX9(2),VX9(2),WX9(2),QZ9(2),WZ9(2)
            WRITE(20,1100)W3,QX3(0),VX3(0),WX3(0),QZ3(0),WZ3(0)
	*                      ,QX3(1),VX3(1),WX3(1),QZ3(1),WZ3(1)
	*                      ,QX3(2),VX3(2),WX3(2),QZ3(2),WZ3(2)
          END IF
	    goto 100
	  end if
c
c     * 2 x 9 Points Newton-Cotes Quadrature *
c     * 1-10-6-11-4-12-7-13-2 *  +  * 2-14-8-15-5-16-9-17-3 *
c
	  index=4
        W10=W3-dw*0.9375d0
        CALL INTEGS1(OM,W10,QX10,QZ10,VX10,WX10,WZ10)
        W11=W3-dw*0.8125d0
        CALL INTEGS1(OM,W11,QX11,QZ11,VX11,WX11,WZ11)
        W12=W3-dw*0.6875d0
        CALL INTEGS1(OM,W12,QX12,QZ12,VX12,WX12,WZ12)
        W13=W3-dw*0.5625d0
        CALL INTEGS1(OM,W13,QX13,QZ13,VX13,WX13,WZ13)
        W14=W3-dw*0.4375d0
        CALL INTEGS1(OM,W14,QX14,QZ14,VX14,WX14,WZ14)
        W15=W3-dw*0.3125d0
        CALL INTEGS1(OM,W15,QX15,QZ15,VX15,WX15,WZ15)
        W16=W3-dw*0.1875d0
        CALL INTEGS1(OM,W16,QX16,QZ16,VX16,WX16,WZ16)
        W17=W3-dw*0.0625d0
        CALL INTEGS1(OM,W17,QX17,QZ17,VX17,WX17,WZ17)
	  nint=nint+8
	  aa=   989.0* QZ1(2)+5888.0*QZ10(2)-  928.0* QZ6(2)
	*     +10496.0*QZ11(2)-4540.0* QZ4(2)+10496.0*QZ12(2)
     *       -928.0* QZ7(2)+5888.0*QZ13(2)+  989.0* QZ2(2)
	  aa=aa+989.0* QZ2(2)+5888.0*QZ14(2)-  928.0* QZ8(2)
     *     +10496.0*QZ15(2)-4540.0* QZ5(2)+10496.0*QZ16(2)
	*       -928.0* QZ9(2)+5888.0*QZ17(2)+  989.0* QZ3(2)
	  sI41=cdabs(aa)*H/56700.0
	  aa=   989.0* VX1(2)+5888.0*VX10(2)-  928.0* VX6(2)
	*     +10496.0*VX11(2)-4540.0* VX4(2)+10496.0*VX12(2)
     *       -928.0* VX7(2)+5888.0*VX13(2)+  989.0* VX2(2)
	  aa=aa+989.0* VX2(2)+5888.0*VX14(2)-  928.0* VX8(2)
     *     +10496.0*VX15(2)-4540.0* VX5(2)+10496.0*VX16(2)
	*       -928.0* VX9(2)+5888.0*VX17(2)+  989.0* VX3(2)
	  sI42=cdabs(aa)*H/56700.0
	  sInt1=sI41
	  sInt2=sI42
	  if(sI41.ge.sImax1)sImax1=sI41
	  if(sI42.ge.sImax2)sImax2=sI42
	  er31=dabs((sI31-sI41)/sI41)
	  er32=dabs((sI32-sI42)/sI42)
c
c	  write(*,*)'i,index=',i,index
c	  write(*,*)'sI11,sI12',sI11,sI12
c 	  write(*,*)'sI21,sI22',sI21,sI22
c 	  write(*,*)'sI31,sI32',sI31,sI32
c 	  write(*,*)'sI41,sI42',sI41,sI42
c	  write(*,*)'er31,er32',er31,er32
	  if((er31.gt.0.0001d0).or.(er32.gt.0.0001d0))then
	    if(((sI41/sImax1).gt.0.0001d0).
     *                      or.((sI42/sImax2).lt.0.0001d0))then
c	      write(*,*)'sI41/sImax1,sI42/sImax2=',sI41/sImax1,sI42/sImax2
c	      write(*,*)'Increase Number of ns for Static Integ 1 !!!!!'
	    end if
	  end if
c
        do J=0,2
	      aa=   989.0* QX1(J)+5888.0*QX10(J)-  928.0* QX6(J)
	*         +10496.0*QX11(J)-4540.0* QX4(J)+10496.0*QX12(J)
     *           -928.0* QX7(J)+5888.0*QX13(J)+  989.0* QX2(J)
	      aa=aa+989.0* QX2(J)+5888.0*QX14(J)-  928.0* QX8(J)
     *         +10496.0*QX15(J)-4540.0* QX5(J)+10496.0*QX16(J)
	*           -928.0* QX9(J)+5888.0*QX17(J)+  989.0* QX3(J)
            GQX(J)= GQX(J)+H*aa/56700.0
	      aa=   989.0* VX1(J)+5888.0*VX10(J)-  928.0* VX6(J)
	*         +10496.0*VX11(J)-4540.0* VX4(J)+10496.0*VX12(J)
     *           -928.0* VX7(J)+5888.0*VX13(J)+  989.0* VX2(J)
	      aa=aa+989.0* VX2(J)+5888.0*VX14(J)-  928.0* VX8(J)
     *         +10496.0*VX15(J)-4540.0* VX5(J)+10496.0*VX16(J)
	*           -928.0* VX9(J)+5888.0*VX17(J)+  989.0* VX3(J)
            GVX(J)= GVX(J)+H*aa/56700.0
	      aa=   989.0* WX1(J)+5888.0*WX10(J)-  928.0* WX6(J)
	*         +10496.0*WX11(J)-4540.0* WX4(J)+10496.0*WX12(J)
     *           -928.0* WX7(J)+5888.0*WX13(J)+  989.0* WX2(J)
	      aa=aa+989.0* WX2(J)+5888.0*WX14(J)-  928.0* WX8(J)
     *         +10496.0*WX15(J)-4540.0* WX5(J)+10496.0*WX16(J)
	*           -928.0* WX9(J)+5888.0*WX17(J)+  989.0* WX3(J)
            GWX(J)= GWX(J)+H*aa/56700.0
	      aa=   989.0* QZ1(J)+5888.0*QZ10(J)-  928.0* QZ6(J)
	*         +10496.0*QZ11(J)-4540.0* QZ4(J)+10496.0*QZ12(J)
     *           -928.0* QZ7(J)+5888.0*QZ13(J)+  989.0* QZ2(J)
	      aa=aa+989.0* QZ2(J)+5888.0*QZ14(J)-  928.0* QZ8(J)
     *         +10496.0*QZ15(J)-4540.0* QZ5(J)+10496.0*QZ16(J)
	*           -928.0* QZ9(J)+5888.0*QZ17(J)+  989.0* QZ3(J)
            GQZ(J)= GQZ(J)+H*aa/56700.0
	      aa=   989.0* WZ1(J)+5888.0*WZ10(J)-  928.0* WZ6(J)
	*         +10496.0*WZ11(J)-4540.0* WZ4(J)+10496.0*WZ12(J)
     *           -928.0* WZ7(J)+5888.0*WZ13(J)+  989.0* WZ2(J)
	      aa=aa+989.0* WZ2(J)+5888.0*WZ14(J)-  928.0* WZ8(J)
     *         +10496.0*WZ15(J)-4540.0* WZ5(J)+10496.0*WZ16(J)
	*           -928.0* WZ9(J)+5888.0*WZ17(J)+  989.0* WZ3(J)
            GWZ(J)= GWZ(J)+H*aa/56700.0
            QX1(J)= QX3(J)
            VX1(J)= VX3(J)
            WX1(J)= WX3(J)
            QZ1(J)= QZ3(J)
            WZ1(J)= WZ3(J)
        end do
          IF(IOT.eq.1)THEN
            WRITE(20,1100)W6,QX6(0),VX6(0),WX6(0),QZ6(0),WZ6(0)
	*                      ,QX6(1),VX6(1),WX6(1),QZ6(1),WZ6(1)
	*                      ,QX6(2),VX6(2),WX6(2),QZ6(2),WZ6(2)
            WRITE(20,1100)W4,QX4(0),VX4(0),WX4(0),QZ4(0),WZ4(0)
	*                      ,QX4(1),VX4(1),WX4(1),QZ4(1),WZ4(1)
	*                      ,QX4(2),VX4(2),WX4(2),QZ4(2),WZ4(2)
            WRITE(20,1100)W7,QX7(0),VX7(0),WX7(0),QZ7(0),WZ7(0)
	*                      ,QX7(1),VX7(1),WX7(1),QZ7(1),WZ7(1)
	*                      ,QX7(2),VX7(2),WX7(2),QZ7(2),WZ7(2)
            WRITE(20,1100)W2,QX2(0),VX2(0),WX2(0),QZ2(0),WZ2(0)
	*                      ,QX2(1),VX2(1),WX2(1),QZ2(1),WZ2(1)
	*                      ,QX2(2),VX2(2),WX2(2),QZ2(2),WZ2(2)
            WRITE(20,1100)W8,QX8(0),VX8(0),WX8(0),QZ8(0),WZ8(0)
	*                      ,QX8(1),VX8(1),WX8(1),QZ8(1),WZ8(1)
	*                      ,QX8(2),VX8(2),WX8(2),QZ8(2),WZ8(2)
            WRITE(20,1100)W5,QX5(0),VX5(0),WX5(0),QZ5(0),WZ5(0)
	*                      ,QX5(1),VX5(1),WX5(1),QZ5(1),WZ5(1)
	*                      ,QX5(2),VX5(2),WX5(2),QZ5(2),WZ5(2)
            WRITE(20,1100)W9,QX9(0),VX9(0),WX9(0),QZ9(0),WZ9(0)
	*                      ,QX9(1),VX9(1),WX9(1),QZ9(1),WZ9(1)
	*                      ,QX9(2),VX9(2),WX9(2),QZ9(2),WZ9(2)
            WRITE(20,1100)W3,QX3(0),VX3(0),WX3(0),QZ3(0),WZ3(0)
	*                      ,QX3(1),VX3(1),WX3(1),QZ3(1),WZ3(1)
	*                      ,QX3(2),VX3(2),WX3(2),QZ3(2),WZ3(2)
          END IF

  100 CONTINUE
c
c	write(*,*)
c	write(*,*)'sInt1       ,sInt2        ',sInt1,sInt2
c	write(*,*)'sImax1      ,sImax2       ',sImax1,sImax2
c	write(*,*)'sInt1/sImax1,sInt2/sImax2=',sInt1/sImax1,sInt2/sImax2
	if(((sInt1/sImax1).lt.0.00001d0)
	*            .and.((sInt2/sImax2).lt.0.00001d0))then
c	  write(*,*)'sInt1/sImax1 & sInt2/sImax2 < 0.00001d !!!'
c	  write(*,*)'NWN2 is set to be 0 !!!'
	  NWN2=0
	  goto 20
	end if
c
	if(((sInt1/sImax1).lt.0.5d0).and.((sInt2/sImax2).lt.0.5d0))then
	  if(iter.gt.20)then
	    write(10,*)'iter > 20 !!'
          write(10,*)'SDEP,ODEP,LS,LO',SDEP,ODEP,LS,LO
          write(10,*)'OM,RR',OM,RR
	    goto 20
	  end if
c	  write(*,*)'sInt1/sImax1 & sInt2/sImax2 < 0.5d !!!' 
c	  write(*,*) 
	  w1=WNRs(2)*iter
	  w2=w1+WNRs(2)
	  if(ns.ge.2)then
	    ns=ns/2
	    sImax1=sImax1*2.d0
	    sImax2=sImax2*2.d0
	  end if
	  goto 10
	end if
c
c   20	write(*,*)
c      write(*,*)'Static Integration 1:'
c	write(*,*)'Number of Integration Range =',iter
c   20 write(*,*)'Number of Integ. Points for Static Range 1=',nint
   20 write(10,*)'Itr. & Intg. Pnt. Num. for Static Range 1=',iter,nint
      NWT1=NWT1+nint
c
      return
	end
C
C
C *** THE INTEGRANDS OF Q, V, AND W AND THEIR DERIVATIVES for Static Range 1 ***
C
      SUBROUTINE INTEGS1(OM,RK,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),OM
     *   ,VLX(2),URX(2,2),URZ(2,2),VLXZ,URXZ(2),URZZ(2),SM,RAM,WN
      DIMENSION BJ0(0:1),BJ1(0:1),BJK(0:1),BJD(0:1)
      COMMON /SOLY/SDEP,ODEP,LS,LO 
     *       /OTST/IOT,IST /FCQVW/SM,RAM,RR
C
c   * the displacement-stress vectors *
c
	WN=RK/RR
      CALL RTSHc(OM,WN,VLX)
      CALL RTPSVc(OM,WN,URX,URZ)
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
      call bj01(RK,BJ0(0),BJ1(0))
      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)
      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)
      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
      RETURN
      END
C
c
C
C *** THE COMPUTATION OF Q, V, AND W AND THEIR DERIVATIVES 
C     USING Simpson's INTEGATION FOMULA (for Static Range 1)***
C
      SUBROUTINE QVWRZS2(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)
     *          ,QX1(0:2),QZ1(0:2),VX1(0:2),WX1(0:2),WZ1(0:2)
     *          ,QX2(0:2),QZ2(0:2),VX2(0:2),WX2(0:2),WZ2(0:2)
     *          ,QX3(0:2),QZ3(0:2),VX3(0:2),WX3(0:2),WZ3(0:2)
     *          ,QX4(0:2),QZ4(0:2),VX4(0:2),WX4(0:2),WZ4(0:2)
     *          ,QX5(0:2),QZ5(0:2),VX5(0:2),WX5(0:2),WZ5(0:2)
     *          ,QX6(0:2),QZ6(0:2),VX6(0:2),WX6(0:2),WZ6(0:2)
     *          ,QX7(0:2),QZ7(0:2),VX7(0:2),WX7(0:2),WZ7(0:2)
     *          ,QX8(0:2),QZ8(0:2),VX8(0:2),WX8(0:2),WZ8(0:2)
     *          ,QX9(0:2),QZ9(0:2),VX9(0:2),WX9(0:2),WZ9(0:2)
     *          ,QX10(0:2),QZ10(0:2),VX10(0:2),WX10(0:2),WZ10(0:2)
     *          ,QX11(0:2),QZ11(0:2),VX11(0:2),WX11(0:2),WZ11(0:2)
     *          ,QX12(0:2),QZ12(0:2),VX12(0:2),WX12(0:2),WZ12(0:2)
     *          ,QX13(0:2),QZ13(0:2),VX13(0:2),WX13(0:2),WZ13(0:2)
     *          ,QX14(0:2),QZ14(0:2),VX14(0:2),WX14(0:2),WZ14(0:2)
     *          ,QX15(0:2),QZ15(0:2),VX15(0:2),WX15(0:2),WZ15(0:2)
     *          ,QX16(0:2),QZ16(0:2),VX16(0:2),WX16(0:2),WZ16(0:2)
     *          ,QX17(0:2),QZ17(0:2),VX17(0:2),WX17(0:2),WZ17(0:2)
     *          ,OM,SM,RAM,aa
      COMMON /WVNUMs/WNRs(3),NWN1o,NWN2o,NWN1,NWN2,NWT1,NWT2
     *       /OTST/IOT,IST /FCQVW/SM,RAM,RR
C
C   * Simpson's INTEGATION FOMULA (for Static Range 2) *
C
      ns=NWN2
      if(ns.eq.0)return
c
      WR=WNRs(2)
      w2=WNRs(3)
	w1=0.d0
	iter=0
      sImax1=0.d0
      sImax2=0.d0
	nint=0
c      write(*,*)
c
   10 iter=iter+1
c      write(*,*)
c      write(*,*)'iter=',iter
      dw=(w2-w1)/ns
	H=dw/RR
      CALL INTEGS2(OM,WR,W1,QX1,QZ1,VX1,WX1,WZ1)
	nint=nint+1
      IF(IOT.eq.1)THEN
        WRITE(21,1100)W1,QX1(0),VX1(0),WX1(0),QZ1(0),WZ1(0)
	*                  ,QX1(1),VX1(1),WX1(1),QZ1(1),WZ1(1)
	*                  ,QX1(2),VX1(2),WX1(2),QZ1(2),WZ1(2)
      END IF
 1100 FORMAT(31E12.4)
c
      DO  100 I=1,ns
c
c     * 3 Points Newton-Cotes (Simpson) Quadrature *
c     * 1-2-3 *
c
	  index=1
        W3=W1+dw*dfloat(I)
        W2=W3-dw*0.5d0
        CALL INTEGS2(OM,WR,W2,QX2,QZ2,VX2,WX2,WZ2)
        CALL INTEGS2(OM,WR,W3,QX3,QZ3,VX3,WX3,WZ3)
	  nint=nint+2
	  sI11=cdabs((QZ1(2)+4.0*QZ2(2)+QZ3(2)))*H/6.0
	  sI12=cdabs((VX1(2)+4.0*VX2(2)+VX3(2)))*H/6.0
	  sInt1=sI11
	  sInt2=sI12
	  if(sI11.ge.sImax1)sImax1=sI11
	  if(sI12.ge.sImax2)sImax2=sI12
	  if(((sI11/sImax1).lt.0.001d0).and.
     *                       ((sI12/sImax2).lt.0.001d0))then
c	    write(*,*)'i,index=',i,index
c	    write(*,*)'sI11,sI12',sI11,sI12
c	    write(*,*)'sImax1,sImax2',sImax1,sImax2
c	    write(*,*)'sI11/sImax1 < 0.001d !!!: sI11/sImax1=',sI11/sImax1
c	    write(*,*)'sI12/sImax2 < 0.001d !!!: sI12/sImax2=',sI12/sImax2
          DO J=0,2
            GQX(J)= GQX(J)+H*(QX1(J)+4.0*QX2(J)+QX3(J))/6.0
            GVX(J)= GVX(J)+H*(VX1(J)+4.0*VX2(J)+VX3(J))/6.0
            GWX(J)= GWX(J)+H*(WX1(J)+4.0*WX2(J)+WX3(J))/6.0
            GQZ(J)= GQZ(J)+H*(QZ1(J)+4.0*QZ2(J)+QZ3(J))/6.0
            GWZ(J)= GWZ(J)+H*(WZ1(J)+4.0*WZ2(J)+WZ3(J))/6.0
            QX1(J)= QX3(J)
            VX1(J)= VX3(J)
            WX1(J)= WX3(J)
            QZ1(J)= QZ3(J)
            WZ1(J)= WZ3(J)
          end do
          IF(IOT.eq.1)THEN
            WRITE(21,1100)W2,QX2(0),VX2(0),WX2(0),QZ2(0),WZ2(0)
	*                      ,QX2(1),VX2(1),WX2(1),QZ2(1),WZ2(1)
	*                      ,QX2(2),VX2(2),WX2(2),QZ2(2),WZ2(2)
            WRITE(21,1100)W3,QX3(0),VX3(0),WX3(0),QZ3(0),WZ3(0)
	*                      ,QX3(1),VX3(1),WX3(1),QZ3(1),WZ3(1)
	*                      ,QX3(2),VX3(2),WX3(2),QZ3(2),WZ3(2)
          END IF
	    goto 100
	  end if
c
c     * 5 Points Newton-Cotes Quadrature *
c     * 1-4-2-5-3 *
c
	  index=2
        W4=W3-dw*0.75d0
        CALL INTEGS2(OM,WR,W4,QX4,QZ4,VX4,WX4,WZ4)
        W5=W3-dw*0.25d0
        CALL INTEGS2(OM,WR,W5,QX5,QZ5,VX5,WX5,WZ5)
	  nint=nint+2
	  aa=7.0*QZ1(2)+32.0*QZ4(2)+12.0*QZ2(2)+32.0*QZ5(2)+7.0*QZ3(2)
	  sI21=cdabs(aa)*H/90.0
	  aa=7.0*VX1(2)+32.0*VX4(2)+12.0*VX2(2)+32.0*VX5(2)+7.0*VX3(2)
	  sI22=cdabs(aa)*H/90.0
	  sInt1=sI21
	  sInt2=sI22
	  if(sI21.ge.sImax1)sImax1=sI21
	  if(sI22.ge.sImax2)sImax2=sI22
	  er11=dabs((sI11-sI21)/sI21)
	  er12=dabs((sI12-sI22)/sI22)
	  if((er11.lt.0.0001d0).and.(er12.lt.0.0001d0))then
c	    write(*,*)'i,index=',i,index
c	    write(*,*)'sI11,sI12',sI11,sI12
c	    write(*,*)'sI21,sI22',sI21,sI22
c	    write(*,*)'er11,er12',er11,er12
          do J=0,2
            aa= 7.0*QX1(J)+32.0*QX4(J)+12.0*QX2(J)
            GQX(J)= GQX(J)+H*(aa+32.0*QX5(J)+7.0*QX3(J))/90.0
            aa= 7.0*VX1(J)+32.0*VX4(J)+12.0*VX2(J)
            GVX(J)= GVX(J)+H*(aa+32.0*VX5(J)+7.0*VX3(J))/90.0
            aa= 7.0*WX1(J)+32.0*WX4(J)+12.0*WX2(J)
            GWX(J)= GWX(J)+H*(aa+32.0*WX5(J)+7.0*WX3(J))/90.0
            aa= 7.0*QZ1(J)+32.0*QZ4(J)+12.0*QZ2(J)
            GQZ(J)= GQZ(J)+H*(aa+32.0*QZ5(J)+7.0*QZ3(J))/90.0
            aa= 7.0*WZ1(J)+32.0*WZ4(J)+12.0*WZ2(J)
            GWZ(J)= GWZ(J)+H*(aa+32.0*WZ5(J)+7.0*WZ3(J))/90.0
            QX1(J)= QX3(J)
            VX1(J)= VX3(J)
            WX1(J)= WX3(J)
            QZ1(J)= QZ3(J)
            WZ1(J)= WZ3(J)
          end do
          IF(IOT.eq.1)THEN
            WRITE(21,1100)W4,QX4(0),VX4(0),WX4(0),QZ4(0),WZ4(0)
	*                      ,QX4(1),VX4(1),WX4(1),QZ4(1),WZ4(1)
	*                      ,QX4(2),VX4(2),WX4(2),QZ4(2),WZ4(2)
            WRITE(21,1100)W2,QX2(0),VX2(0),WX2(0),QZ2(0),WZ2(0)
	*                      ,QX2(1),VX2(1),WX2(1),QZ2(1),WZ2(1)
	*                      ,QX2(2),VX2(2),WX2(2),QZ2(2),WZ2(2)
            WRITE(21,1100)W5,QX5(0),VX5(0),WX5(0),QZ5(0),WZ5(0)
	*                      ,QX5(1),VX5(1),WX5(1),QZ5(1),WZ5(1)
	*                      ,QX5(2),VX5(2),WX5(2),QZ5(2),WZ5(2)
            WRITE(21,1100)W3,QX3(0),VX3(0),WX3(0),QZ3(0),WZ3(0)
	*                      ,QX3(1),VX3(1),WX3(1),QZ3(1),WZ3(1)
	*                      ,QX3(2),VX3(2),WX3(2),QZ3(2),WZ3(2)
          END IF
	    goto 100
	  end if
c
c     * 9 Points Simpson *
c     * 1-6-4-7-2-8-5-9-3 *
c
	  index=3
        W6=W3-dw*0.875d0
        CALL INTEGS2(OM,WR,W6,QX6,QZ6,VX6,WX6,WZ6)
        W7=W3-dw*0.625d0
        CALL INTEGS2(OM,WR,W7,QX7,QZ7,VX7,WX7,WZ7)
        W8=W3-dw*0.375d0
        CALL INTEGS2(OM,WR,W8,QX8,QZ8,VX8,WX8,WZ8)
        W9=W3-dw*0.125d0
        CALL INTEGS2(OM,WR,W9,QX9,QZ9,VX9,WX9,WZ9)
	  nint=nint+4
	  aa=989.0*QZ1(2)+5888.0*QZ6(2)-928.0*QZ4(2)+10496.0*QZ7(2)-4540.0
	*   *QZ2(2)+10496.0*QZ8(2)-928.0*QZ5(2)+5888.0*QZ9(2)+989.0*QZ3(2)
	  sI31=cdabs(aa)*H/28350.0
	  aa=989.0*VX1(2)+5888.0*VX6(2)-928.0*VX4(2)+10496.0*VX7(2)-4540.0
	*   *VX2(2)+10496.0*VX8(2)-928.0*VX5(2)+5888.0*VX9(2)+989.0*VX3(2)
	  sI32=cdabs(aa)*H/28350.0
	  sInt1=sI31
	  sInt2=sI32
	  if(sI31.ge.sImax1)sImax1=sI31
	  if(sI32.ge.sImax2)sImax2=sI32
	  er21=dabs((sI21-sI31)/sI31)
	  er22=dabs((sI22-sI32)/sI32)
	  if((er21.lt.0.0001d0).and.(er22.lt.0.0001d0))then
c	    write(*,*)'i,index=',i,index
c	    write(*,*)'sI11,sI12',sI11,sI12
c 	    write(*,*)'sI21,sI22',sI21,sI22
c 	    write(*,*)'sI31,sI32',sI31,sI32
c	    write(*,*)'er21,er22',er21,er22
          do J=0,2
            aa= 989.0*QX1(J)+5888.0*QX6(J)-928.0*QX4(J)+10496.0*QX7(J)
            aa= aa-4540.0*QX2(J)+10496.0*QX8(J)-928.0*QX5(J)
            GQX(J)= GQX(J)+H*(aa+5888.0*QX9(J)+989.0*QX3(J))/28350.0
            aa= 989.0*VX1(J)+5888.0*VX6(J)-928.0*VX4(J)+10496.0*VX7(J)
            aa= aa-4540.0*VX2(J)+10496.0*VX8(J)-928.0*VX5(J)
            GVX(J)= GVX(J)+H*(aa+5888.0*VX9(J)+989.0*VX3(J))/28350.0
            aa= 989.0*WX1(J)+5888.0*WX6(J)-928.0*WX4(J)+10496.0*WX7(J)
            aa= aa-4540.0*WX2(J)+10496.0*WX8(J)-928.0*WX5(J)
            GWX(J)= GWX(J)+H*(aa+5888.0*WX9(J)+989.0*WX3(J))/28350.0
            aa= 989.0*QZ1(J)+5888.0*QZ6(J)-928.0*QZ4(J)+10496.0*QZ7(J)
            aa= aa-4540.0*QZ2(J)+10496.0*QZ8(J)-928.0*QZ5(J)
            GQZ(J)= GQZ(J)+H*(aa+5888.0*QZ9(J)+989.0*QZ3(J))/28350.0
            aa= 989.0*WZ1(J)+5888.0*WZ6(J)-928.0*WZ4(J)+10496.0*WZ7(J)
            aa= aa-4540.0*WZ2(J)+10496.0*WZ8(J)-928.0*WZ5(J)
            GWZ(J)= GWZ(J)+H*(aa+5888.0*WZ9(J)+989.0*WZ3(J))/28350.0
            QX1(J)= QX3(J)
            VX1(J)= VX3(J)
            WX1(J)= WX3(J)
            QZ1(J)= QZ3(J)
            WZ1(J)= WZ3(J)
          end do
          IF(IOT.eq.1)THEN
            WRITE(21,1100)W6,QX6(0),VX6(0),WX6(0),QZ6(0),WZ6(0)
	*                      ,QX6(1),VX6(1),WX6(1),QZ6(1),WZ6(1)
	*                      ,QX6(2),VX6(2),WX6(2),QZ6(2),WZ6(2)
            WRITE(21,1100)W4,QX4(0),VX4(0),WX4(0),QZ4(0),WZ4(0)
	*                      ,QX4(1),VX4(1),WX4(1),QZ4(1),WZ4(1)
	*                      ,QX4(2),VX4(2),WX4(2),QZ4(2),WZ4(2)
            WRITE(21,1100)W7,QX7(0),VX7(0),WX7(0),QZ7(0),WZ7(0)
	*                      ,QX7(1),VX7(1),WX7(1),QZ7(1),WZ7(1)
	*                      ,QX7(2),VX7(2),WX7(2),QZ7(2),WZ7(2)
            WRITE(21,1100)W2,QX2(0),VX2(0),WX2(0),QZ2(0),WZ2(0)
	*                      ,QX2(1),VX2(1),WX2(1),QZ2(1),WZ2(1)
	*                      ,QX2(2),VX2(2),WX2(2),QZ2(2),WZ2(2)
            WRITE(21,1100)W8,QX8(0),VX8(0),WX8(0),QZ8(0),WZ8(0)
	*                      ,QX8(1),VX8(1),WX8(1),QZ8(1),WZ8(1)
	*                      ,QX8(2),VX8(2),WX8(2),QZ8(2),WZ8(2)
            WRITE(21,1100)W5,QX5(0),VX5(0),WX5(0),QZ5(0),WZ5(0)
	*                      ,QX5(1),VX5(1),WX5(1),QZ5(1),WZ5(1)
	*                      ,QX5(2),VX5(2),WX5(2),QZ5(2),WZ5(2)
            WRITE(21,1100)W9,QX9(0),VX9(0),WX9(0),QZ9(0),WZ9(0)
	*                      ,QX9(1),VX9(1),WX9(1),QZ9(1),WZ9(1)
	*                      ,QX9(2),VX9(2),WX9(2),QZ9(2),WZ9(2)
            WRITE(21,1100)W3,QX3(0),VX3(0),WX3(0),QZ3(0),WZ3(0)
	*                      ,QX3(1),VX3(1),WX3(1),QZ3(1),WZ3(1)
	*                      ,QX3(2),VX3(2),WX3(2),QZ3(2),WZ3(2)
          END IF
	    goto 100
	  end if
c
c     * 2 x 9 Points Newton-Cotes Quadrature *
c     * 1-10-6-11-4-12-7-13-2 *  +  * 2-14-8-15-5-16-9-17-3 *
c
	  index=4
        W10=W3-dw*0.9375d0
        CALL INTEGS2(OM,WR,W10,QX10,QZ10,VX10,WX10,WZ10)
        W11=W3-dw*0.8125d0
        CALL INTEGS2(OM,WR,W11,QX11,QZ11,VX11,WX11,WZ11)
        W12=W3-dw*0.6875d0
        CALL INTEGS2(OM,WR,W12,QX12,QZ12,VX12,WX12,WZ12)
        W13=W3-dw*0.5625d0
        CALL INTEGS2(OM,WR,W13,QX13,QZ13,VX13,WX13,WZ13)
        W14=W3-dw*0.4375d0
        CALL INTEGS2(OM,WR,W14,QX14,QZ14,VX14,WX14,WZ14)
        W15=W3-dw*0.3125d0
        CALL INTEGS2(OM,WR,W15,QX15,QZ15,VX15,WX15,WZ15)
        W16=W3-dw*0.1875d0
        CALL INTEGS2(OM,WR,W16,QX16,QZ16,VX16,WX16,WZ16)
        W17=W3-dw*0.0625d0
        CALL INTEGS2(OM,WR,W17,QX17,QZ17,VX17,WX17,WZ17)
	  nint=nint+8
	  aa=   989.0* QZ1(2)+5888.0*QZ10(2)-  928.0* QZ6(2)
	*     +10496.0*QZ11(2)-4540.0* QZ4(2)+10496.0*QZ12(2)
     *       -928.0* QZ7(2)+5888.0*QZ13(2)+  989.0* QZ2(2)
	  aa=aa+989.0* QZ2(2)+5888.0*QZ14(2)-  928.0* QZ8(2)
     *     +10496.0*QZ15(2)-4540.0* QZ5(2)+10496.0*QZ16(2)
	*       -928.0* QZ9(2)+5888.0*QZ17(2)+  989.0* QZ3(2)
	  sI41=cdabs(aa)*H/56700.0
	  aa=   989.0* VX1(2)+5888.0*VX10(2)-  928.0* VX6(2)
	*     +10496.0*VX11(2)-4540.0* VX4(2)+10496.0*VX12(2)
     *       -928.0* VX7(2)+5888.0*VX13(2)+  989.0* VX2(2)
	  aa=aa+989.0* VX2(2)+5888.0*VX14(2)-  928.0* VX8(2)
     *     +10496.0*VX15(2)-4540.0* VX5(2)+10496.0*VX16(2)
	*       -928.0* VX9(2)+5888.0*VX17(2)+  989.0* VX3(2)
	  sI42=cdabs(aa)*H/56700.0
	  sInt1=sI41
	  sInt2=sI42
	  if(sI41.ge.sImax1)sImax1=sI41
	  if(sI42.ge.sImax2)sImax2=sI42
	  er31=dabs((sI31-sI41)/sI41)
	  er32=dabs((sI32-sI42)/sI42)
c
c	  write(*,*)'i,index=',i,index
c	  write(*,*)'sI11,sI12',sI11,sI12
c 	  write(*,*)'sI21,sI22',sI21,sI22
c 	  write(*,*)'sI31,sI32',sI31,sI32
c 	  write(*,*)'sI41,sI42',sI41,sI42
c	  write(*,*)'er31,er32',er31,er32
	  if((er31.gt.0.0001d0).or.(er32.gt.0.0001d0))then
	    if(((sI41/sImax1).gt.0.0001d0).
     *                      or.((sI42/sImax2).lt.0.0001d0))then
c	      write(*,*)'sI41/sImax1,sI42/sImax2=',sI41/sImax1,sI42/sImax2
c	      write(*,*)'Increase Number of ns for Static Integ 2 !!!!!'
	    end if
	  end if
c
        do J=0,2
	      aa=   989.0* QX1(J)+5888.0*QX10(J)-  928.0* QX6(J)
	*         +10496.0*QX11(J)-4540.0* QX4(J)+10496.0*QX12(J)
     *           -928.0* QX7(J)+5888.0*QX13(J)+  989.0* QX2(J)
	      aa=aa+989.0* QX2(J)+5888.0*QX14(J)-  928.0* QX8(J)
     *         +10496.0*QX15(J)-4540.0* QX5(J)+10496.0*QX16(J)
	*           -928.0* QX9(J)+5888.0*QX17(J)+  989.0* QX3(J)
            GQX(J)= GQX(J)+H*aa/56700.0
	      aa=   989.0* VX1(J)+5888.0*VX10(J)-  928.0* VX6(J)
	*         +10496.0*VX11(J)-4540.0* VX4(J)+10496.0*VX12(J)
     *           -928.0* VX7(J)+5888.0*VX13(J)+  989.0* VX2(J)
	      aa=aa+989.0* VX2(J)+5888.0*VX14(J)-  928.0* VX8(J)
     *         +10496.0*VX15(J)-4540.0* VX5(J)+10496.0*VX16(J)
	*           -928.0* VX9(J)+5888.0*VX17(J)+  989.0* VX3(J)
            GVX(J)= GVX(J)+H*aa/56700.0
	      aa=   989.0* WX1(J)+5888.0*WX10(J)-  928.0* WX6(J)
	*         +10496.0*WX11(J)-4540.0* WX4(J)+10496.0*WX12(J)
     *           -928.0* WX7(J)+5888.0*WX13(J)+  989.0* WX2(J)
	      aa=aa+989.0* WX2(J)+5888.0*WX14(J)-  928.0* WX8(J)
     *         +10496.0*WX15(J)-4540.0* WX5(J)+10496.0*WX16(J)
	*           -928.0* WX9(J)+5888.0*WX17(J)+  989.0* WX3(J)
            GWX(J)= GWX(J)+H*aa/56700.0
	      aa=   989.0* QZ1(J)+5888.0*QZ10(J)-  928.0* QZ6(J)
	*         +10496.0*QZ11(J)-4540.0* QZ4(J)+10496.0*QZ12(J)
     *           -928.0* QZ7(J)+5888.0*QZ13(J)+  989.0* QZ2(J)
	      aa=aa+989.0* QZ2(J)+5888.0*QZ14(J)-  928.0* QZ8(J)
     *         +10496.0*QZ15(J)-4540.0* QZ5(J)+10496.0*QZ16(J)
	*           -928.0* QZ9(J)+5888.0*QZ17(J)+  989.0* QZ3(J)
            GQZ(J)= GQZ(J)+H*aa/56700.0
	      aa=   989.0* WZ1(J)+5888.0*WZ10(J)-  928.0* WZ6(J)
	*         +10496.0*WZ11(J)-4540.0* WZ4(J)+10496.0*WZ12(J)
     *           -928.0* WZ7(J)+5888.0*WZ13(J)+  989.0* WZ2(J)
	      aa=aa+989.0* WZ2(J)+5888.0*WZ14(J)-  928.0* WZ8(J)
     *         +10496.0*WZ15(J)-4540.0* WZ5(J)+10496.0*WZ16(J)
	*           -928.0* WZ9(J)+5888.0*WZ17(J)+  989.0* WZ3(J)
            GWZ(J)= GWZ(J)+H*aa/56700.0
            QX1(J)= QX3(J)
            VX1(J)= VX3(J)
            WX1(J)= WX3(J)
            QZ1(J)= QZ3(J)
            WZ1(J)= WZ3(J)
        end do
        IF(IOT.eq.1)THEN
            WRITE(21,1100)W10,QX10(0),VX10(0),WX10(0),QZ10(0),WZ10(0)
	*                       ,QX10(1),VX10(1),WX10(1),QZ10(1),WZ10(1)
	*                       ,QX10(2),VX10(2),WX10(2),QZ10(2),WZ10(2)
            WRITE(21,1100)W6,QX6(0),VX6(0),WX6(0),QZ6(0),WZ6(0)
	*                      ,QX6(1),VX6(1),WX6(1),QZ6(1),WZ6(1)
	*                      ,QX6(2),VX6(2),WX6(2),QZ6(2),WZ6(2)
            WRITE(21,1100)W11,QX11(0),VX11(0),WX11(0),QZ11(0),WZ11(0)
	*                       ,QX11(1),VX11(1),WX11(1),QZ11(1),WZ11(1)
	*                       ,QX11(2),VX11(2),WX11(2),QZ11(2),WZ11(2)
            WRITE(21,1100)W4,QX4(0),VX4(0),WX4(0),QZ4(0),WZ4(0)
	*                      ,QX4(1),VX4(1),WX4(1),QZ4(1),WZ4(1)
	*                      ,QX4(2),VX4(2),WX4(2),QZ4(2),WZ4(2)
            WRITE(21,1100)W12,QX12(0),VX12(0),WX12(0),QZ12(0),WZ12(0)
	*                       ,QX12(1),VX12(1),WX12(1),QZ12(1),WZ12(1)
	*                       ,QX12(2),VX12(2),WX12(2),QZ12(2),WZ12(2)
            WRITE(21,1100)W7,QX7(0),VX7(0),WX7(0),QZ7(0),WZ7(0)
	*                      ,QX7(1),VX7(1),WX7(1),QZ7(1),WZ7(1)
	*                      ,QX7(2),VX7(2),WX7(2),QZ7(2),WZ7(2)
            WRITE(21,1100)W13,QX13(0),VX13(0),WX13(0),QZ13(0),WZ13(0)
	*                       ,QX13(1),VX13(1),WX13(1),QZ13(1),WZ13(1)
	*                       ,QX13(2),VX13(2),WX13(2),QZ13(2),WZ13(2)
            WRITE(21,1100)W2,QX2(0),VX2(0),WX2(0),QZ2(0),WZ2(0)
	*                      ,QX2(1),VX2(1),WX2(1),QZ2(1),WZ2(1)
	*                      ,QX2(2),VX2(2),WX2(2),QZ2(2),WZ2(2)
            WRITE(21,1100)W14,QX14(0),VX14(0),WX14(0),QZ14(0),WZ14(0)
	*                       ,QX14(1),VX14(1),WX14(1),QZ14(1),WZ14(1)
	*                       ,QX14(2),VX14(2),WX14(2),QZ14(2),WZ14(2)
            WRITE(21,1100)W8,QX8(0),VX8(0),WX8(0),QZ8(0),WZ8(0)
	*                      ,QX8(1),VX8(1),WX8(1),QZ8(1),WZ8(1)
	*                      ,QX8(2),VX8(2),WX8(2),QZ8(2),WZ8(2)
            WRITE(21,1100)W15,QX15(0),VX15(0),WX15(0),QZ15(0),WZ15(0)
	*                       ,QX15(1),VX15(1),WX15(1),QZ15(1),WZ15(1)
	*                       ,QX15(2),VX15(2),WX15(2),QZ15(2),WZ15(2)
            WRITE(21,1100)W5,QX5(0),VX5(0),WX5(0),QZ5(0),WZ5(0)
	*                      ,QX5(1),VX5(1),WX5(1),QZ5(1),WZ5(1)
	*                      ,QX5(2),VX5(2),WX5(2),QZ5(2),WZ5(2)
            WRITE(21,1100)W16,QX16(0),VX16(0),WX16(0),QZ16(0),WZ16(0)
	*                       ,QX16(1),VX16(1),WX16(1),QZ16(1),WZ16(1)
	*                       ,QX16(2),VX16(2),WX16(2),QZ16(2),WZ16(2)
            WRITE(21,1100)W9,QX9(0),VX9(0),WX9(0),QZ9(0),WZ9(0)
	*                      ,QX9(1),VX9(1),WX9(1),QZ9(1),WZ9(1)
	*                      ,QX9(2),VX9(2),WX9(2),QZ9(2),WZ9(2)
            WRITE(21,1100)W17,QX17(0),VX17(0),WX17(0),QZ17(0),WZ17(0)
	*                       ,QX17(1),VX17(1),WX17(1),QZ17(1),WZ17(1)
	*                       ,QX17(2),VX17(2),WX17(2),QZ17(2),WZ17(2)
            WRITE(21,1100)W3,QX3(0),VX3(0),WX3(0),QZ3(0),WZ3(0)
	*                      ,QX3(1),VX3(1),WX3(1),QZ3(1),WZ3(1)
	*                      ,QX3(2),VX3(2),WX3(2),QZ3(2),WZ3(2)
          END IF
c	    goto 100

  100 CONTINUE
c
c	write(*,*)
c	write(*,*)'sInt1       ,sInt2        ',sInt1,sInt2
c	write(*,*)'sImax1      ,sImax2       ',sImax1,sImax2
c	write(*,*)'sInt1/sImax1,sInt2/sImax2=',sInt1/sImax1,sInt2/sImax2
	if(((sInt1/sImax1).lt.0.00001d0)
	*            .and.((sInt2/sImax2).lt.0.00001d0))then
c	  write(*,*)'sInt1/sImax1 & sInt2/sImax2 < 0.00001d !!!'
	  goto 20
	end if
c
	if(((sInt1/sImax1).lt.0.5d0).and.((sInt2/sImax2).lt.0.5d0))then
	  if(iter.gt.20)then
	    write(10,*)'iter > 20 !!'
	    goto 20
	  end if
c	  write(*,*)'sInt1/sImax1 & sInt2/sImax2 < 0.5d !!!' 
c	  write(*,*) 
	  w1=WNRs(3)*iter
	  w2=w1+WNRs(3)
	  if(ns.ge.2)then
	    ns=ns/2
	    sImax1=sImax1*2.d0
	    sImax2=sImax2*2.d0
	  end if
	  goto 10
	end if
c
c   20	write(*,*)
c      write(*,*)'Static Integration 2:'
c	write(*,*)'Number of Integration Range =',iter
c   20 write(*,*)'Number of Integ. Points for Static Range 2=',nint
   20 write(10,*)'Itr. & Intg. Pnt. Num. for Static Range 2=',iter,nint
      NWT2=NWT2+nint
c
      return
	end
C
C
C *** THE INTEGRANDS OF Q, V, AND W AND THEIR DERIVATIVES for Static Range 1 ***
C
      SUBROUTINE INTEGS2(OM,WR,WI,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),OM
     *   ,VLX(2),URX(2,2),URZ(2,2),VLXZ,URXZ(2),URZZ(2),SM,RAM
     *   ,WN,im,RK,BJ0,BJ1,BN0,BN1,HK0(0:1),HK1(0:1),HKK(0:1),HKD(0:1)
      COMMON /SOLY/SDEP,ODEP,LS,LO 
     *       /OTST/IOT,IST /FCQVW/SM,RAM,RR
C
c   * the displacement-stress vectors on the plus ki axis *
c
      im=(0.d0,1.d0)
      RK=dcmplx(WR,+WI)
      WN=RK/RR
      CALL RTSHc(OM,WN,VLX)
      CALL RTPSVc(OM,WN,URX,URZ)
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 complex Bessel & Hankell functions on the plus ki axis *
c
      call CBJN0(RK,BJ0,BN0)
      call CBJN1(RK,BJ1,BN1)
      HK0(0)=+0.5d0*im*(BJ0+im*BN0)
      HK1(0)=+0.5d0*im*(BJ1+im*BN1)
c
      HKK(0)=HK1(0)/RK
      HKD(0)=HK0(0)-HKK(0)
      HK0(1)=-HK1(0)*WN
      HK1(1)= HKD(0)*WN
      HKK(1)=(HK1(1)-HK1(0)/RR)/RK
      HKD(1)= HK0(1)-HKK(1)
C
c   * the integrands *
c
      QX(0)= URX(1,1)*HKD(0)+VLX(1)*HKK(0)
      VX(0)= URX(1,1)*HKK(0)+VLX(1)*HKD(0)
      WX(0)= URX(2,1)*HK1(0)
      QZ(0)= URZ(1,1)*HK1(0)
      WZ(0)= URZ(2,1)*HK0(0)
      QX(1)= URX(1,1)*HKD(1)+VLX(1)*HKK(1)
      VX(1)= URX(1,1)*HKK(1)+VLX(1)*HKD(1)
      WX(1)= URX(2,1)*HK1(1)
      QZ(1)= URZ(1,1)*HK1(1)
      WZ(1)= URZ(2,1)*HK0(1)
      QX(2)= URXZ(1)*HKD(0)+VLXZ*HKK(0)
      VX(2)= URXZ(1)*HKK(0)+VLXZ*HKD(0)	 
      WX(2)= URXZ(2)*HK1(0)
      QZ(2)= URZZ(1)*HK1(0)
      WZ(2)= URZZ(2)*HK0(0)
C
c   * the displacement-stress vectors on the minus ki axis *
c
      RK=dcmplx(WR,-WI)
      WN=RK/RR
      CALL RTSHc(OM,WN,VLX)
      CALL RTPSVc(OM,WN,URX,URZ)
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 complex Bessel & Hankell functions on the minus ki axis *
c
      HK0(0)=dconjg(HK0(0))
      HK1(0)=dconjg(HK1(0))
c
      HKK(0)=HK1(0)/RK
      HKD(0)=HK0(0)-HKK(0)
      HK0(1)=-HK1(0)*WN
      HK1(1)= HKD(0)*WN
      HKK(1)=(HK1(1)-HK1(0)/RR)/RK
      HKD(1)= HK0(1)-HKK(1)
C
c   * the integrands *
c
      QX(0)= QX(0) + URX(1,1)*HKD(0)+VLX(1)*HKK(0)
      VX(0)= VX(0) + URX(1,1)*HKK(0)+VLX(1)*HKD(0)
      WX(0)= WX(0) + URX(2,1)*HK1(0)
      QZ(0)= QZ(0) + URZ(1,1)*HK1(0)
      WZ(0)= WZ(0) + URZ(2,1)*HK0(0)
      QX(1)= QX(1) + URX(1,1)*HKD(1)+VLX(1)*HKK(1)
      VX(1)= VX(1) + URX(1,1)*HKK(1)+VLX(1)*HKD(1)
      WX(1)= WX(1) + URX(2,1)*HK1(1)
      QZ(1)= QZ(1) + URZ(1,1)*HK1(1)
      WZ(1)= WZ(1) + URZ(2,1)*HK0(1)
      QX(2)= QX(2) + URXZ(1)*HKD(0)+VLXZ*HKK(0)
      VX(2)= VX(2) + URXZ(1)*HKK(0)+VLXZ*HKD(0)	 
      WX(2)= WX(2) + URXZ(2)*HK1(0)
      QZ(2)= QZ(2) + URZZ(1)*HK1(0)
      WZ(2)= WZ(2) + URZZ(2)*HK0(0)
C
      RETURN
      END
C
c
c
C *** SH DISPLACEMENT USING R/T MATRIX ***
C
      SUBROUTINE RTSHc(OM,WN,VLX)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VLX(2),WS2,CMS,SM,EE(2,2),OM,OM2,VSc,VS2
     *          ,DD,DU,AU,BU,AD,BD,CS,GTD,GTU,GRD,GRU
     *          ,CD(20),CU(20),CUSU,CDSU,CDSD,CUSD,WN,WN2,DQ
      COMMON /SOLY/SDEP,ODEP,LS,LO
     *    /STRC/DNS(20),VP(20),QP(20),FP(20)
     *                 ,VS(20),QS(20),FS(20),THK(20),NL
     *       /GRTSH/GTD(0:20),GTU(0:20),GRD(0:20),GRU(0:20)
C
C *** CHECK FOR THE DIMENSION OF CD AND CU ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(10,*)' FROM RTSHc: NL.GT.NC ',NL,NC
         STOP
      END IF
C
C *** COMPUTATION OF GENERALIZED R/T MATRICES ***
C
      CALL GNRTSHc(OM,WN)
C
C *** COMPUTATION OF VLX ***
C
      OM2=OM*OM
      WN2=WN*WN
C
      VSc=VS(LS)*DCMPLX(1.D0,-0.5D0/QS(LS))
      VS2=VSc*VSc
      WS2=OM2/VS2
      CMS=CDSQRT(WN2-WS2)
      IF(dREAL(CMS).LT.0.0)CMS=-CMS
C
      SM=DNS(LS)*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
      DQ=-WN/6.2831853072d0
C
c *** Up-/Down-Going Coefficients ***
c
      IF((LS.EQ.1).AND.(SDEP.EQ.0.0))THEN
C
          IF(NL.EQ.1)THEN
            CDSD=DQ/EE(2,1)
            CUSD=(0.0,0.0)
          ELSE
            DD=CDEXP(-CMS*THK(1))
            CDSD=DQ/(EE(2,1)+EE(2,2)*DD*GRD(1))
            CUSD=GRD(1)*CDSD
          END IF
C
      ELSE
C 
          IF(LS.EQ.NL)THEN 
C
            DD=CDEXP(-CMS*SDEP)
            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
            CUSU=DQ/(BD*CS-BU)
            CDSU=GRU(LS-1)*CUSU
            CDSD=CS*CUSU
C
          ELSE
C
            DD=CDEXP(-CMS*SDEP)
            DU=CDEXP(-CMS*(THK(LS)-SDEP))
            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
            CUSU=DQ/(BD*CS-BU)
            CDSU=GRU(LS-1)*CUSU
            CDSD=CS*CUSU
            CUSD=GRD(LS)*CDSD
C
          ENDIF
C
      END IF
C
      IF(LS.EQ.1)GOTO 10
      CU(LS-1)=GTU(LS-1)*CUSU
      DO 100 J=LS-2,1,-1
          CU(J)=GTU(J)*CU(J+1)
          CD(J+1)=GRU(J)*CU(J+1)
  100 CONTINUE
      CD(1)=GRU(0)*CU(1)
C
   10 IF(LS.EQ.NL)GOTO 20
      CD(LS+1)=GTD(LS)*CDSD
      DO 120 J=LS+1,NL-1
          CD(J+1)=GTD(J)*CD(J)
          CU(J)=GRD(J)*CD(J)
  120 CONTINUE
C
c *** Displacement-Stress Vectors for SH wave ***
c
   20 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(LO)
          VLX(2)=EE(2,1)*DD*CD(LO)
        ELSE
          IF(ODEP.GE.SDEP)THEN
            DD=CDEXP(-CMS*(ODEP-SDEP))
            VLX(1)=EE(1,1)*DD*CDSD
            VLX(2)=EE(2,1)*DD*CDSD
          ELSE
            DD=CDEXP(-CMS*ODEP)
            DU=CDEXP(-CMS*(SDEP-ODEP))
            VLX(1)=EE(1,1)*DD*CDSU+EE(1,2)*DU*CUSU
            VLX(2)=EE(2,1)*DD*CDSU+EE(2,2)*DU*CUSU
          END IF
        END IF
C
      ELSE
C        
        IF(LO.NE.LS)THEN
          DD=CDEXP(-CMS*ODEP)
          DU=CDEXP(-CMS*(THK(LO)-ODEP))
          VLX(1)=EE(1,1)*DD*CD(LO)+EE(1,2)*DU*CU(LO)
          VLX(2)=EE(2,1)*DD*CD(LO)+EE(2,2)*DU*CU(LO)
        ELSE
          IF(ODEP.GE.SDEP)THEN
            DD=CDEXP(-CMS*(ODEP-SDEP))
            DU=CDEXP(-CMS*(THK(LO)-ODEP))
            VLX(1)=EE(1,1)*DD*CDSD+EE(1,2)*DU*CUSD
            VLX(2)=EE(2,1)*DD*CDSD+EE(2,2)*DU*CUSD
          ELSE
            DD=CDEXP(-CMS*ODEP)
            DU=CDEXP(-CMS*(SDEP-ODEP))
            VLX(1)=EE(1,1)*DD*CDSU+EE(1,2)*DU*CUSU
            VLX(2)=EE(2,1)*DD*CDSU+EE(2,2)*DU*CUSU
          END IF
        END IF
C
      END IF
C
      RETURN
      END
c
C
C
C *** PSV DISPLACEMENT USING R/T MATRIX ***
C
      SUBROUTINE RTPSVc(OM,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
     *    ,AU(2,2),BU(2,2),AD(2,2),BD(2,2),CS(2,2),E1(2,2),E2(2,2)
	*    ,GTD,GTU,GRD,GRU,DT,URX(2,2),URZ(2,2),WN,DQX,DQZ
     *    ,CDX(2,20),CUX(2,20),CUSUX(2),CDSUX(2),CDSDX(2),CUSDX(2) 
     *    ,CDZ(2,20),CUZ(2,20),CUSUZ(2),CDSUZ(2),CDSDZ(2),CUSDZ(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
     * /GRTPSV/GTD(2,2,0:20),GTU(2,2,0:20),GRD(2,2,0:20),GRU(2,2,0:20)
C
C *** CHECK FOR THE DIMENSION OF CD AND CU ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(10,*)' FROM RTPSV: NL.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C *** COMPUTATION OF GENERALIZED R/T MATRICES ***
C
      CALL GNRTPSVc(OM,WN)
C
C *** COMPUTATION OF CD AND CU ***
C
      CALL MTXPSVc(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
            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
            CDSDX(1)=E2(1,1)*DQX
            CDSDX(2)=E2(2,1)*DQX
            CUSDX(1)=(0.0,0.0)
            CUSDX(2)=(0.0,0.0)
            CDSDZ(1)=E2(1,2)*DQZ
            CDSDZ(2)=E2(2,2)*DQZ
            CUSDZ(1)=(0.0,0.0)
            CUSDZ(2)=(0.0,0.0)
        ELSE
            CALL DDUPSVc(OM,WN,THK(1),CMS,CMP,DU)
            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
            CDSDX(1)=E2(1,1)*DQX
            CDSDX(2)=E2(2,1)*DQX
            CUSDX(1)=GRD(1,1,1)*CDSDX(1)+GRD(1,2,1)*CDSDX(2)
            CUSDX(2)=GRD(2,1,1)*CDSDX(1)+GRD(2,2,1)*CDSDX(2)
            CDSDZ(1)=E2(1,2)*DQZ
            CDSDZ(2)=E2(2,2)*DQZ
            CUSDZ(1)=GRD(1,1,1)*CDSDZ(1)+GRD(1,2,1)*CDSDZ(2)
            CUSDZ(2)=GRD(2,1,1)*CDSDZ(1)+GRD(2,2,1)*CDSDZ(2)
        END IF
C 
      ELSE
C
        IF(LS.EQ.NL)THEN 
            CALL DDUPSVc(OM,WN,SDEP,CMS,CMP,DD)
            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)  
            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)
            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)
            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
            CUSUX(1)=E2(1,1)*DQX
            CUSUX(2)=E2(2,1)*DQX
            CDSUX(1)=GRU(1,1,LS-1)*CUSUX(1)+GRU(1,2,LS-1)*CUSUX(2)
            CDSUX(2)=GRU(2,1,LS-1)*CUSUX(1)+GRU(2,2,LS-1)*CUSUX(2)
            CDSDX(1)=CS(1,1)*CUSUX(1)+CS(1,2)*CUSUX(2)
            CDSDX(2)=CS(2,1)*CUSUX(1)+CS(2,2)*CUSUX(2)
            CUSUZ(1)=E2(1,2)*DQZ
            CUSUZ(2)=E2(2,2)*DQZ
            CDSUZ(1)=GRU(1,1,LS-1)*CUSUZ(1)+GRU(1,2,LS-1)*CUSUZ(2)
            CDSUZ(2)=GRU(2,1,LS-1)*CUSUZ(1)+GRU(2,2,LS-1)*CUSUZ(2)
            CDSDZ(1)=CS(1,1)*CUSUZ(1)+CS(1,2)*CUSUZ(2)
            CDSDZ(2)=CS(2,1)*CUSUZ(1)+CS(2,2)*CUSUZ(2)
        ELSE
            CALL DDUPSVc(OM,WN,SDEP,CMS,CMP,DD)
            CALL DDUPSVc(OM,WN,THK(LS)-SDEP,CMS,CMP,DU)
            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)
            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)
            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)
            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)
            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)
            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
            CUSUX(1)=E2(1,1)*DQX
            CUSUX(2)=E2(2,1)*DQX
            CDSUX(1)=GRU(1,1,LS-1)*CUSUX(1)+GRU(1,2,LS-1)*CUSUX(2)
            CDSUX(2)=GRU(2,1,LS-1)*CUSUX(1)+GRU(2,2,LS-1)*CUSUX(2)
            CDSDX(1)=CS(1,1)*CUSUX(1)+CS(1,2)*CUSUX(2)
            CDSDX(2)=CS(2,1)*CUSUX(1)+CS(2,2)*CUSUX(2)
            CUSDX(1)=GRD(1,1,LS)*CDSDX(1)+GRD(1,2,LS)*CDSDX(2)
            CUSDX(2)=GRD(2,1,LS)*CDSDX(1)+GRD(2,2,LS)*CDSDX(2)
            CUSUZ(1)=E2(1,2)*DQZ
            CUSUZ(2)=E2(2,2)*DQZ
            CDSUZ(1)=GRU(1,1,LS-1)*CUSUZ(1)+GRU(1,2,LS-1)*CUSUZ(2)
            CDSUZ(2)=GRU(2,1,LS-1)*CUSUZ(1)+GRU(2,2,LS-1)*CUSUZ(2)
            CDSDZ(1)=CS(1,1)*CUSUZ(1)+CS(1,2)*CUSUZ(2)
            CDSDZ(2)=CS(2,1)*CUSUZ(1)+CS(2,2)*CUSUZ(2)
            CUSDZ(1)=GRD(1,1,LS)*CDSDZ(1)+GRD(1,2,LS)*CDSDZ(2)
            CUSDZ(2)=GRD(2,1,LS)*CDSDZ(1)+GRD(2,2,LS)*CDSDZ(2)
        ENDIF
C
      END IF
C
      IF(LS.EQ.1)GOTO 10
        CUX(1,LS-1)=GTU(1,1,LS-1)*CUSUX(1)+GTU(1,2,LS-1)*CUSUX(2)
        CUX(2,LS-1)=GTU(2,1,LS-1)*CUSUX(1)+GTU(2,2,LS-1)*CUSUX(2)
        CUZ(1,LS-1)=GTU(1,1,LS-1)*CUSUZ(1)+GTU(1,2,LS-1)*CUSUZ(2)
        CUZ(2,LS-1)=GTU(2,1,LS-1)*CUSUZ(1)+GTU(2,2,LS-1)*CUSUZ(2)
        DO 100 I=LS-2,1,-1
          CUX(1,I)  =GTU(1,1,I)*CUX(1,I+1)+GTU(1,2,I)*CUX(2,I+1)
          CUX(2,I)  =GTU(2,1,I)*CUX(1,I+1)+GTU(2,2,I)*CUX(2,I+1)
          CDX(1,I+1)=GRU(1,1,I)*CUX(1,I+1)+GRU(1,2,I)*CUX(2,I+1)
          CDX(2,I+1)=GRU(2,1,I)*CUX(1,I+1)+GRU(2,2,I)*CUX(2,I+1)
          CUZ(1,I)  =GTU(1,1,I)*CUZ(1,I+1)+GTU(1,2,I)*CUZ(2,I+1)
          CUZ(2,I)  =GTU(2,1,I)*CUZ(1,I+1)+GTU(2,2,I)*CUZ(2,I+1)
          CDZ(1,I+1)=GRU(1,1,I)*CUZ(1,I+1)+GRU(1,2,I)*CUZ(2,I+1)
          CDZ(2,I+1)=GRU(2,1,I)*CUZ(1,I+1)+GRU(2,2,I)*CUZ(2,I+1)
  100   CONTINUE
        CDX(1,1)=GRU(1,1,0)*CUX(1,1)+GRU(1,2,0)*CUX(2,1)
        CDX(2,1)=GRU(2,1,0)*CUX(1,1)+GRU(2,2,0)*CUX(2,1)
        CDZ(1,1)=GRU(1,1,0)*CUZ(1,1)+GRU(1,2,0)*CUZ(2,1)
        CDZ(2,1)=GRU(2,1,0)*CUZ(1,1)+GRU(2,2,0)*CUZ(2,1)
C
   10 IF(LS.EQ.NL)GOTO 20
        CDX(1,LS+1)=GTD(1,1,LS)*CDSDX(1)+GTD(1,2,LS)*CDSDX(2)
        CDX(2,LS+1)=GTD(2,1,LS)*CDSDX(1)+GTD(2,2,LS)*CDSDX(2)
        CDZ(1,LS+1)=GTD(1,1,LS)*CDSDZ(1)+GTD(1,2,LS)*CDSDZ(2)
        CDZ(2,LS+1)=GTD(2,1,LS)*CDSDZ(1)+GTD(2,2,LS)*CDSDZ(2)
        DO 120 I=LS+1,NL-1
          CDX(1,I+1)=GTD(1,1,I)*CDX(1,I)+GTD(1,2,I)*CDX(2,I)
          CDX(2,I+1)=GTD(2,1,I)*CDX(1,I)+GTD(2,2,I)*CDX(2,I)
          CUX(1,I)  =GRD(1,1,I)*CDX(1,I)+GRD(1,2,I)*CDX(2,I)
          CUX(2,I)  =GRD(2,1,I)*CDX(1,I)+GRD(2,2,I)*CDX(2,I)
          CDZ(1,I+1)=GTD(1,1,I)*CDZ(1,I)+GTD(1,2,I)*CDZ(2,I)
          CDZ(2,I+1)=GTD(2,1,I)*CDZ(1,I)+GTD(2,2,I)*CDZ(2,I)
          CUZ(1,I)  =GRD(1,1,I)*CDZ(1,I)+GRD(1,2,I)*CDZ(2,I)
          CUZ(2,I)  =GRD(2,1,I)*CDZ(1,I)+GRD(2,2,I)*CDZ(2,I)
  120   CONTINUE
C
C *** COMPUTATION OF URX AND URZ ***
C
  20  CALL MTXPSVc(OM,WN,LO,CMS,CMP,EE)
c
      IF(LO.EQ.NL)THEN
C        
        IF(LO.NE.LS)THEN
C
          CALL DDUPSVc(OM,WN,ODEP,CMS,CMP,DD)
C
          E1(1,1)=EE(1,1)*DD(1,1)+EE(1,2)*DD(2,1)
          E1(1,2)=EE(1,1)*DD(1,2)+EE(1,2)*DD(2,2)
          E1(2,1)=EE(2,1)*DD(1,1)+EE(2,2)*DD(2,1)
          E1(2,2)=EE(2,1)*DD(1,2)+EE(2,2)*DD(2,2)
          URX(1,1)=E1(1,1)*CDX(1,LO)+E1(1,2)*CDX(2,LO)
          URX(2,1)=E1(2,1)*CDX(1,LO)+E1(2,2)*CDX(2,LO)
          URZ(1,1)=E1(1,1)*CDZ(1,LO)+E1(1,2)*CDZ(2,LO)
          URZ(2,1)=E1(2,1)*CDZ(1,LO)+E1(2,2)*CDZ(2,LO)
C
          E1(1,1)=EE(3,1)*DD(1,1)+EE(3,2)*DD(2,1)
          E1(1,2)=EE(3,1)*DD(1,2)+EE(3,2)*DD(2,2)
          E1(2,1)=EE(4,1)*DD(1,1)+EE(4,2)*DD(2,1)
          E1(2,2)=EE(4,1)*DD(1,2)+EE(4,2)*DD(2,2)
          URX(1,2)=E1(1,1)*CDX(1,LO)+E1(1,2)*CDX(2,LO)
          URX(2,2)=E1(2,1)*CDX(1,LO)+E1(2,2)*CDX(2,LO)
          URZ(1,2)=E1(1,1)*CDZ(1,LO)+E1(1,2)*CDZ(2,LO)
          URZ(2,2)=E1(2,1)*CDZ(1,LO)+E1(2,2)*CDZ(2,LO)
C
        ELSE
C
          IF(ODEP.GE.SDEP)THEN
C
            CALL DDUPSVc(OM,WN,ODEP-SDEP,CMS,CMP,DD)
C
            E1(1,1)=EE(1,1)*DD(1,1)+EE(1,2)*DD(2,1)
            E1(1,2)=EE(1,1)*DD(1,2)+EE(1,2)*DD(2,2)
            E1(2,1)=EE(2,1)*DD(1,1)+EE(2,2)*DD(2,1)
            E1(2,2)=EE(2,1)*DD(1,2)+EE(2,2)*DD(2,2)
            URX(1,1)=E1(1,1)*CDSDX(1)+E1(1,2)*CDSDX(2)
            URX(2,1)=E1(2,1)*CDSDX(1)+E1(2,2)*CDSDX(2)
            URZ(1,1)=E1(1,1)*CDSDZ(1)+E1(1,2)*CDSDZ(2)
            URZ(2,1)=E1(2,1)*CDSDZ(1)+E1(2,2)*CDSDZ(2)
C
            E1(1,1)=EE(3,1)*DD(1,1)+EE(3,2)*DD(2,1)
            E1(1,2)=EE(3,1)*DD(1,2)+EE(3,2)*DD(2,2)
            E1(2,1)=EE(4,1)*DD(1,1)+EE(4,2)*DD(2,1)
            E1(2,2)=EE(4,1)*DD(1,2)+EE(4,2)*DD(2,2)
            URX(1,2)=E1(1,1)*CDSDX(1)+E1(1,2)*CDSDX(2)
            URX(2,2)=E1(2,1)*CDSDX(1)+E1(2,2)*CDSDX(2)
            URZ(1,2)=E1(1,1)*CDSDZ(1)+E1(1,2)*CDSDZ(2)
            URZ(2,2)=E1(2,1)*CDSDZ(1)+E1(2,2)*CDSDZ(2)
C
          ELSE
C
            CALL DDUPSVc(OM,WN,ODEP,CMS,CMP,DD)
            CALL DDUPSVc(OM,WN,SDEP-ODEP,CMS,CMP,DU)
C
            E1(1,1)=EE(1,1)*DD(1,1)+EE(1,2)*DD(2,1)
            E1(1,2)=EE(1,1)*DD(1,2)+EE(1,2)*DD(2,2)
            E1(2,1)=EE(2,1)*DD(1,1)+EE(2,2)*DD(2,1)
            E1(2,2)=EE(2,1)*DD(1,2)+EE(2,2)*DD(2,2)
            E2(1,1)=EE(1,3)*DU(1,1)+EE(1,4)*DU(2,1)
            E2(1,2)=EE(1,3)*DU(1,2)+EE(1,4)*DU(2,2)
            E2(2,1)=EE(2,3)*DU(1,1)+EE(2,4)*DU(2,1)
            E2(2,2)=EE(2,3)*DU(1,2)+EE(2,4)*DU(2,2)
C
            URX(1,1)=E1(1,1)*CDSUX(1)+E1(1,2)*CDSUX(2)
     *              +E2(1,1)*CUSUX(1)+E2(1,2)*CUSUX(2)
            URX(2,1)=E1(2,1)*CDSUX(1)+E1(2,2)*CDSUX(2)
     *              +E2(2,1)*CUSUX(1)+E2(2,2)*CUSUX(2)
            URZ(1,1)=E1(1,1)*CDSUZ(1)+E1(1,2)*CDSUZ(2)
     *              +E2(1,1)*CUSUZ(1)+E2(1,2)*CUSUZ(2)
            URZ(2,1)=E1(2,1)*CDSUZ(1)+E1(2,2)*CDSUZ(2)
     *              +E2(2,1)*CUSUZ(1)+E2(2,2)*CUSUZ(2)
C
            E1(1,1)=EE(3,1)*DD(1,1)+EE(3,2)*DD(2,1)
            E1(1,2)=EE(3,1)*DD(1,2)+EE(3,2)*DD(2,2)
            E1(2,1)=EE(4,1)*DD(1,1)+EE(4,2)*DD(2,1)
            E1(2,2)=EE(4,1)*DD(1,2)+EE(4,2)*DD(2,2)
            E2(1,1)=EE(3,3)*DU(1,1)+EE(3,4)*DU(2,1)
            E2(1,2)=EE(3,3)*DU(1,2)+EE(3,4)*DU(2,2)
            E2(2,1)=EE(4,3)*DU(1,1)+EE(4,4)*DU(2,1)
            E2(2,2)=EE(4,3)*DU(1,2)+EE(4,4)*DU(2,2)
C
            URX(1,2)=E1(1,1)*CDSUX(1)+E1(1,2)*CDSUX(2)
     *              +E2(1,1)*CUSUX(1)+E2(1,2)*CUSUX(2)
            URX(2,2)=E1(2,1)*CDSUX(1)+E1(2,2)*CDSUX(2)
     *              +E2(2,1)*CUSUX(1)+E2(2,2)*CUSUX(2)
            URZ(1,2)=E1(1,1)*CDSUZ(1)+E1(1,2)*CDSUZ(2)
     *              +E2(1,1)*CUSUZ(1)+E2(1,2)*CUSUZ(2)
            URZ(2,2)=E1(2,1)*CDSUZ(1)+E1(2,2)*CDSUZ(2)
     *              +E2(2,1)*CUSUZ(1)+E2(2,2)*CUSUZ(2)
C
          END IF
C
        END IF
C
      ELSE
C        
        IF(LO.NE.LS)THEN
C
            CALL DDUPSVc(OM,WN,ODEP,CMS,CMP,DD)
            CALL DDUPSVc(OM,WN,THK(LO)-ODEP,CMS,CMP,DU)
C
            E1(1,1)=EE(1,1)*DD(1,1)+EE(1,2)*DD(2,1)
            E1(1,2)=EE(1,1)*DD(1,2)+EE(1,2)*DD(2,2)
            E1(2,1)=EE(2,1)*DD(1,1)+EE(2,2)*DD(2,1)
            E1(2,2)=EE(2,1)*DD(1,2)+EE(2,2)*DD(2,2)
            E2(1,1)=EE(1,3)*DU(1,1)+EE(1,4)*DU(2,1)
            E2(1,2)=EE(1,3)*DU(1,2)+EE(1,4)*DU(2,2)
            E2(2,1)=EE(2,3)*DU(1,1)+EE(2,4)*DU(2,1)
            E2(2,2)=EE(2,3)*DU(1,2)+EE(2,4)*DU(2,2)
C
            URX(1,1)=E1(1,1)*CDX(1,LO)+E1(1,2)*CDX(2,LO)
     *              +E2(1,1)*CUX(1,LO)+E2(1,2)*CUX(2,LO)
            URX(2,1)=E1(2,1)*CDX(1,LO)+E1(2,2)*CDX(2,LO)
     *              +E2(2,1)*CUX(1,LO)+E2(2,2)*CUX(2,LO)
            URZ(1,1)=E1(1,1)*CDZ(1,LO)+E1(1,2)*CDZ(2,LO)
     *              +E2(1,1)*CUZ(1,LO)+E2(1,2)*CUZ(2,LO)
            URZ(2,1)=E1(2,1)*CDZ(1,LO)+E1(2,2)*CDZ(2,LO)
     *              +E2(2,1)*CUZ(1,LO)+E2(2,2)*CUZ(2,LO)
C
            E1(1,1)=EE(3,1)*DD(1,1)+EE(3,2)*DD(2,1)
            E1(1,2)=EE(3,1)*DD(1,2)+EE(3,2)*DD(2,2)
            E1(2,1)=EE(4,1)*DD(1,1)+EE(4,2)*DD(2,1)
            E1(2,2)=EE(4,1)*DD(1,2)+EE(4,2)*DD(2,2)
            E2(1,1)=EE(3,3)*DU(1,1)+EE(3,4)*DU(2,1)
            E2(1,2)=EE(3,3)*DU(1,2)+EE(3,4)*DU(2,2)
            E2(2,1)=EE(4,3)*DU(1,1)+EE(4,4)*DU(2,1)
            E2(2,2)=EE(4,3)*DU(1,2)+EE(4,4)*DU(2,2)
C
            URX(1,2)=E1(1,1)*CDX(1,LO)+E1(1,2)*CDX(2,LO)
     *              +E2(1,1)*CUX(1,LO)+E2(1,2)*CUX(2,LO)
            URX(2,2)=E1(2,1)*CDX(1,LO)+E1(2,2)*CDX(2,LO)
     *              +E2(2,1)*CUX(1,LO)+E2(2,2)*CUX(2,LO)
            URZ(1,2)=E1(1,1)*CDZ(1,LO)+E1(1,2)*CDZ(2,LO)
     *              +E2(1,1)*CUZ(1,LO)+E2(1,2)*CUZ(2,LO)
            URZ(2,2)=E1(2,1)*CDZ(1,LO)+E1(2,2)*CDZ(2,LO)
     *              +E2(2,1)*CUZ(1,LO)+E2(2,2)*CUZ(2,LO)
C
        ELSE
C
          IF(ODEP.GE.SDEP)THEN
C
            CALL DDUPSVc(OM,WN,ODEP-SDEP,CMS,CMP,DD)
            CALL DDUPSVc(OM,WN,THK(LO)-ODEP,CMS,CMP,DU)
C
            E1(1,1)=EE(1,1)*DD(1,1)+EE(1,2)*DD(2,1)
            E1(1,2)=EE(1,1)*DD(1,2)+EE(1,2)*DD(2,2)
            E1(2,1)=EE(2,1)*DD(1,1)+EE(2,2)*DD(2,1)
            E1(2,2)=EE(2,1)*DD(1,2)+EE(2,2)*DD(2,2)
            E2(1,1)=EE(1,3)*DU(1,1)+EE(1,4)*DU(2,1)
            E2(1,2)=EE(1,3)*DU(1,2)+EE(1,4)*DU(2,2)
            E2(2,1)=EE(2,3)*DU(1,1)+EE(2,4)*DU(2,1)
            E2(2,2)=EE(2,3)*DU(1,2)+EE(2,4)*DU(2,2)
C
            URX(1,1)=E1(1,1)*CDSDX(1)+E1(1,2)*CDSDX(2)
     *              +E2(1,1)*CUSDX(1)+E2(1,2)*CUSDX(2)
            URX(2,1)=E1(2,1)*CDSDX(1)+E1(2,2)*CDSDX(2)
     *              +E2(2,1)*CUSDX(1)+E2(2,2)*CUSDX(2)
            URZ(1,1)=E1(1,1)*CDSDZ(1)+E1(1,2)*CDSDZ(2)
     *              +E2(1,1)*CUSDZ(1)+E2(1,2)*CUSDZ(2)
            URZ(2,1)=E1(2,1)*CDSDZ(1)+E1(2,2)*CDSDZ(2)
     *              +E2(2,1)*CUSDZ(1)+E2(2,2)*CUSDZ(2)
C
            E1(1,1)=EE(3,1)*DD(1,1)+EE(3,2)*DD(2,1)
            E1(1,2)=EE(3,1)*DD(1,2)+EE(3,2)*DD(2,2)
            E1(2,1)=EE(4,1)*DD(1,1)+EE(4,2)*DD(2,1)
            E1(2,2)=EE(4,1)*DD(1,2)+EE(4,2)*DD(2,2)
            E2(1,1)=EE(3,3)*DU(1,1)+EE(3,4)*DU(2,1)
            E2(1,2)=EE(3,3)*DU(1,2)+EE(3,4)*DU(2,2)
            E2(2,1)=EE(4,3)*DU(1,1)+EE(4,4)*DU(2,1)
            E2(2,2)=EE(4,3)*DU(1,2)+EE(4,4)*DU(2,2)
C
            URX(1,2)=E1(1,1)*CDSDX(1)+E1(1,2)*CDSDX(2)
     *              +E2(1,1)*CUSDX(1)+E2(1,2)*CUSDX(2)
            URX(2,2)=E1(2,1)*CDSDX(1)+E1(2,2)*CDSDX(2)
     *              +E2(2,1)*CUSDX(1)+E2(2,2)*CUSDX(2)
            URZ(1,2)=E1(1,1)*CDSDZ(1)+E1(1,2)*CDSDZ(2)
     *              +E2(1,1)*CUSDZ(1)+E2(1,2)*CUSDZ(2)
            URZ(2,2)=E1(2,1)*CDSDZ(1)+E1(2,2)*CDSDZ(2)
     *              +E2(2,1)*CUSDZ(1)+E2(2,2)*CUSDZ(2)
C
          ELSE
C
            CALL DDUPSVc(OM,WN,ODEP,CMS,CMP,DD)
            CALL DDUPSVc(OM,WN,SDEP-ODEP,CMS,CMP,DU)
C
            E1(1,1)=EE(1,1)*DD(1,1)+EE(1,2)*DD(2,1)
            E1(1,2)=EE(1,1)*DD(1,2)+EE(1,2)*DD(2,2)
            E1(2,1)=EE(2,1)*DD(1,1)+EE(2,2)*DD(2,1)
            E1(2,2)=EE(2,1)*DD(1,2)+EE(2,2)*DD(2,2)
            E2(1,1)=EE(1,3)*DU(1,1)+EE(1,4)*DU(2,1)
            E2(1,2)=EE(1,3)*DU(1,2)+EE(1,4)*DU(2,2)
            E2(2,1)=EE(2,3)*DU(1,1)+EE(2,4)*DU(2,1)
            E2(2,2)=EE(2,3)*DU(1,2)+EE(2,4)*DU(2,2)
C
            URX(1,1)=E1(1,1)*CDSUX(1)+E1(1,2)*CDSUX(2)
     *              +E2(1,1)*CUSUX(1)+E2(1,2)*CUSUX(2)
            URX(2,1)=E1(2,1)*CDSUX(1)+E1(2,2)*CDSUX(2)
     *              +E2(2,1)*CUSUX(1)+E2(2,2)*CUSUX(2)
            URZ(1,1)=E1(1,1)*CDSUZ(1)+E1(1,2)*CDSUZ(2)
     *              +E2(1,1)*CUSUZ(1)+E2(1,2)*CUSUZ(2)
            URZ(2,1)=E1(2,1)*CDSUZ(1)+E1(2,2)*CDSUZ(2)
     *              +E2(2,1)*CUSUZ(1)+E2(2,2)*CUSUZ(2)
C
            E1(1,1)=EE(3,1)*DD(1,1)+EE(3,2)*DD(2,1)
            E1(1,2)=EE(3,1)*DD(1,2)+EE(3,2)*DD(2,2)
            E1(2,1)=EE(4,1)*DD(1,1)+EE(4,2)*DD(2,1)
            E1(2,2)=EE(4,1)*DD(1,2)+EE(4,2)*DD(2,2)
            E2(1,1)=EE(3,3)*DU(1,1)+EE(3,4)*DU(2,1)
            E2(1,2)=EE(3,3)*DU(1,2)+EE(3,4)*DU(2,2)
            E2(2,1)=EE(4,3)*DU(1,1)+EE(4,4)*DU(2,1)
            E2(2,2)=EE(4,3)*DU(1,2)+EE(4,4)*DU(2,2)
C
            URX(1,2)=E1(1,1)*CDSUX(1)+E1(1,2)*CDSUX(2)
     *              +E2(1,1)*CUSUX(1)+E2(1,2)*CUSUX(2)
            URX(2,2)=E1(2,1)*CDSUX(1)+E1(2,2)*CDSUX(2)
     *              +E2(2,1)*CUSUX(1)+E2(2,2)*CUSUX(2)
            URZ(1,2)=E1(1,1)*CDSUZ(1)+E1(1,2)*CDSUZ(2)
     *              +E2(1,1)*CUSUZ(1)+E2(1,2)*CUSUZ(2)
            URZ(2,2)=E1(2,1)*CDSUZ(1)+E1(2,2)*CDSUZ(2)
     *              +E2(2,1)*CUSUZ(1)+E2(2,2)*CUSUZ(2)
C
          END IF
C
        END IF
C
      END IF
C
      RETURN
      END
C
c
C *** SH GENERALIZED R/T MATIX ***
C
      SUBROUTINE GNRTSHc(OM,WN)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 GTD,GTU,GRD,GRU,TD,TU,RD,RU,OM,WN
      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 MDRTSHc(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 *** SH MODIFIED R/T MATIX ***
C
      SUBROUTINE MDRTSHc(OM,WN)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 WS2,CMS,SM,EA(2,2),EB(2,2),VSc,VS2
     *      ,TD,TU,RD,RU,DDA,DUB,EAB,OM,OM2,WN,WN2
      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
      VSc=VS(1)*DCMPLX(1.D0,-0.5D0/QS(1))
      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)
c        RU(0)=-EA(2,2)*DUB/EA(2,1)
        RU(0)=DUB
        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
c      RU(0)=-EA(2,2)*DUB/EA(2,1)
      RU(0)=DUB
C
      DO 100 I=1, NL-2
C
        VSc=VS(I+1)*DCMPLX(1.D0,-0.5D0/QS(I+1))
        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        
      VSc=VS(NL)*DCMPLX(1.D0,-0.5D0/QS(NL))
      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 GNRTPSVc(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)
      COMPLEX*16 DT,OM,WN
      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 MDRTPSVc(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
C
C *** P-SV MODIFIED R/T MATIX ***
C
      SUBROUTINE MDRTPSVc(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,WN
      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 MDRTPSVc: NL.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C *** MODIFIED R/T MATRIX ***
C
      CALL MTXPSVc(OM,WN,1,CMS,CMP,EA)
C
      IF(LS.EQ.1)THEN
        CALL DDUPSVc(OM,WN,SDEP,CMS,CMP,DUB)
      ELSE
        CALL DDUPSVc(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 DDUPSVc(OM,WN,THK(1)-SDEP,CMS,CMP,DDA)
      ELSE
        CALL DDUPSVc(OM,WN,THK(1),CMS,CMP,DDA)
      END IF
C
      DO 100 I=1, NL-1
C
        CALL MTXPSVc(OM,WN,I+1,CMS,CMP,EB)
C
        IF(LS.EQ.(I+1))THEN
          CALL DDUPSVc(OM,WN,SDEP,CMS,CMP,DUB)
        ELSE
          CALL DDUPSVc(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 DDUPSVc(OM,WN,THK(I+1)-SDEP,CMS,CMP,DDA)
          END IF
        ELSE
          CALL DDUPSVc(OM,WN,THK(I+1),CMS,CMP,DDA)
        END IF
C
  100 CONTINUE
C
      RETURN
      END
C
C
C *** P-SV ELEMENT MATIX ***
C
      SUBROUTINE MTXPSVc(OM,WN,L,CMS,CMP,EE)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 WS2,CMS,WP2,CMP,EE(4,4),X0,AA,CMSG,CMPG
	*          ,OM,OM2,WN,WN2,VSc,VS2,VPc,VP2,SM
      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 MTXPSVc: NL.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C *** ELEMENT MATRIX ***
C
      VPc=VP(L)*DCMPLX(1.D0,-0.5D0/QP(L))
      VSc=VS(L)*DCMPLX(1.D0,-0.5D0/QS(L))
      VP2=VPc*VPc
      VS2=VSc*VSc
      SM=DNS(L)*VS2
C
c      IF(OM.EQ.0.D0)THEN
      IF(dreal(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 DDUPSVc(OM,WN,THK,CMS,CMP,DU)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 DU(2,2),CMS,CMP,OM,WN,EW
C
c      IF(OM.EQ.0.D0)THEN
      IF(dreal(OM).EQ.0.D0)THEN
        EW=cDEXP(-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 *** P-SV MODIFIED R/T MATRIX FOR Ith LAYER ***
C
        SUBROUTINE MTRPSV(EA,EB,DDA,DUB,L) 
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 EA(4,4),EB(4,4),DDA(2,2),DUB(2,2),TD,TU,RD,RU
     *          ,E1(4,4),E2(4,4),E3(4,4)
     *          ,a1(4),a2(4),a3(4),b1(4),b2(4),b3(4)  
     *          ,c1(4),c2(4),c3(4),d1(4),d2(4),d3(4)  
     *          ,aa(4),EE(4,4)
      COMMON /MRTPSV/TD(2,2,0:20),TU(2,2,0:20),RD(2,2,0:20),RU(2,2,0:20)
     *    /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 COMMON ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(10,*)' FROM MTRPSV: NL.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C *** MATRIX ***
C
      DO 100 I=1,2
        DO 110 J=1,2
          EE(I,J)= EB(I,J)
          EE(I,J+2)=-EA(I,J+2)
          EE(I+2,J)= EB(I+2,J)
          EE(I+2,J+2)=-EA(I+2,J+2)
          E2(I,J)= EA(I,J)
          E2(I,J+2)=-EB(I,J+2)
          E2(I+2,J)= EA(I+2,J)
          E2(I+2,J+2)=-EB(I+2,J+2)
  110   CONTINUE
  100 CONTINUE
C
C     * INVERSE MATRIX OF EE *
C
        DO 222 I=1,4
C
        IF(I.EQ.1)THEN
          J=2
          K=3
          M=4
        ELSE
          IF(I.EQ.2)THEN
            J=3
            K=4
            M=1
          ELSE
            IF(I.EQ.3)THEN
              J=4
              K=1
              M=2
            ELSE
              J=1
              K=2
              M=3
            END IF
          END IF
        END IF
C
        A1(I)=EE(I,4)*EE(J,1)-EE(I,1)*EE(J,4)  
        A2(I)=EE(I,4)*EE(J,2)-EE(I,2)*EE(J,4)  
        A3(I)=EE(I,4)*EE(J,3)-EE(I,3)*EE(J,4)  
        B1(I)=EE(I,3)*EE(J,2)-EE(I,2)*EE(J,3)  
        B2(I)=EE(I,1)*EE(J,3)-EE(I,3)*EE(J,1)  
        B3(I)=EE(I,2)*EE(J,1)-EE(I,1)*EE(J,2)  
        C1(I)=EE(K,3)*EE(M,2)-EE(K,2)*EE(M,3)  
        C2(I)=EE(K,1)*EE(M,3)-EE(K,3)*EE(M,1)  
        C3(I)=EE(K,2)*EE(M,1)-EE(K,1)*EE(M,2)  
        D1(I)=EE(K,4)*EE(M,1)-EE(K,1)*EE(M,4)  
        D2(I)=EE(K,4)*EE(M,2)-EE(K,2)*EE(M,4)  
        D3(I)=EE(K,4)*EE(M,3)-EE(K,3)*EE(M,4)  
C
        AA(I)=A1(I)*C1(I)+A2(I)*C2(I)+A3(I)*C3(I)
     *       +B1(I)*D1(I)+B2(I)*D2(I)+B3(I)*D3(I)
C
        E1(1,I)=(-EE(J,2)*D3(I)+EE(J,3)*D2(I)-EE(J,4)*C1(I))/AA(I)
        E1(2,I)=( EE(J,1)*D3(I)-EE(J,3)*D1(I)-EE(J,4)*C2(I))/AA(I)
        E1(3,I)=(-EE(J,1)*D2(I)+EE(J,2)*D1(I)-EE(J,4)*C3(I))/AA(I)
        E1(4,I)=( EE(J,1)*C1(I)+EE(J,2)*C2(I)+EE(J,3)*C3(I))/AA(I)
C
  222   CONTINUE
C
      DO 200 I=1,4
        DO 210 J=1,4
          E3(I,J)=(0.D0,0.D0)
          DO 220 K=1,4
            E3(I,J)=E3(I,J)+E1(I,K)*E2(K,J) 
  220     CONTINUE
  210   CONTINUE
  200 CONTINUE
C
      TD(1,1,L)=E3(1,1)*DDA(1,1)+E3(1,2)*DDA(2,1)  
      TD(2,1,L)=E3(2,1)*DDA(1,1)+E3(2,2)*DDA(2,1)  
      TD(1,2,L)=E3(1,1)*DDA(1,2)+E3(1,2)*DDA(2,2)  
      TD(2,2,L)=E3(2,1)*DDA(1,2)+E3(2,2)*DDA(2,2)  
      RD(1,1,L)=E3(3,1)*DDA(1,1)+E3(3,2)*DDA(2,1)  
      RD(2,1,L)=E3(4,1)*DDA(1,1)+E3(4,2)*DDA(2,1)  
      RD(1,2,L)=E3(3,1)*DDA(1,2)+E3(3,2)*DDA(2,2)  
      RD(2,2,L)=E3(4,1)*DDA(1,2)+E3(4,2)*DDA(2,2)  
C
      IF(L.EQ.NL)RETURN
C
      RU(1,1,L)=E3(1,3)*DUB(1,1)+E3(1,4)*DUB(2,1)  
      RU(2,1,L)=E3(2,3)*DUB(1,1)+E3(2,4)*DUB(2,1)  
      RU(1,2,L)=E3(1,3)*DUB(1,2)+E3(1,4)*DUB(2,2)  
      RU(2,2,L)=E3(2,3)*DUB(1,2)+E3(2,4)*DUB(2,2)  
      TU(1,1,L)=E3(3,3)*DUB(1,1)+E3(3,4)*DUB(2,1)  
      TU(2,1,L)=E3(4,3)*DUB(1,1)+E3(4,4)*DUB(2,1)  
      TU(1,2,L)=E3(3,3)*DUB(1,2)+E3(3,4)*DUB(2,2)  
      TU(2,2,L)=E3(4,3)*DUB(1,2)+E3(4,4)*DUB(2,2)  
C
      RETURN
      END
c
c
c *** Complex Bessel and Neumann Functions of the 1st order ***
c
	subroutine CBJN1 (X,BJ1,BN1)
      
      IMPLICIT REAL*8 (A-H,O-Z)           
      complex*16 X,BJ1,BN1,z,z1,aa,co,si,X1,X2,X3           
      DATA a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15
     *   /  5.0         D-1,  1.666666667 D-1,
     *      8.333333333 D-2,  5.0         D-2,
     *      3.333333333 D-2,  2.380952381 D-2,
     *      1.785714286 D-2,  1.388888889 D-2,
     *      1.111111111 D-2,  9.090909091 D-3,
     *      7.575757576 D-3,  6.41025641  D-3,
     *      5.494505495 D-3,  4.761904762 D-3,
     *      4.166666666 D-3 /
      DATA b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,b15
     * /0.50000000000000 D0, 1.25000000000000 D0,    
     *  1.66666666666666 D0, 1.95833333333333 D0,     
     *  2.18333333333333 D0, 2.36666666666667 D0,     
     *  2.52142857142857 D0, 2.65535714285714 D0,     
     *  2.77341269841270 D0, 2.87896825396825 D0,     
     *  2.97442279942280 D0, 3.06154401154401 D0,     
     *  3.14167221667222 D0, 3.21584804084805 D0,     
     *  3.28489565989566 D0 /    
      DATA pi / 3.141592654 D0 /
      DATA eu / 0.5772157 D0 /
c
      if(cdabs(x).gt.9.d0)goto 10
c
      a=1.d0
      z1=x*0.5d0
      z=z1*z1                       
      BJ1=z1*(a-z*a1*(a-z*a2*(a-z*a3*(a-z*a4*(a-z*a5*(a-z*a6*
     *   (a-z*a7*(a-z*a8*(a-z*a9*(a-z*a10*(a-z*a11*(a-z*a12* 
     *   (a-z*a13*(a-z*a14*(a-z*a15)))))))))))))))
      BN1=z1*(b1-z*a1*(b2-z*a2*(b3-z*a3*(b4-z*a4*(b5-z*a5*(b6-z*a6*
     *   (b7-z*a7*(b8-z*a8*(b9-z*a9*(b10-z*a10*(b11-z*a11*(b12-z*
     *   a12*(b13-z*a13*(b14-z*a14*(b15-z*a15)))))))))))))))
      BN1=2.d0*((cdlog(z1)+eu)*BJ1-1.d0/x-BN1)/pi
      goto 20
c
   10 aa=cdsqrt(pi*x)
      co=cdcos(x)
      si=cdsin(x)
      x1=0.375d0/x
      x2=0.1171875d0/(x*x)
      x3=0.102539063d0/(x*x*x)
      BJ1= (si*(1.d0+x1+x2-x3)-co*(1.d0-x1+x2+x3))/aa 
      BN1=-(co*(1.d0+x1+x2-x3)+si*(1.d0-x1+x2+x3))/aa 
c 
   20 return
      end      
c
c
c
c *** Complex Bessel and Neumann Functions of the 0th order ***
c
      subroutine CBJN0 (X,BJ0,BN0)
c
      IMPLICIT REAL*8 (A-H,O-Z)
      complex*16 X,BJ0,BN0,z,z1,aa,co,si,X1,X2,X3           
      DATA a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15
     *   /  2.5         D-1,  1.111111111 D-1,
     *      6.25        D-2,  4.0         D-2,
     *      2.777777778 D-2,  2.040816327 D-2,
     *      1.5625      D-2,  1.234567901 D-2,
     *      1.0         D-2,  8.26446281  D-3,
     *      6.944444444 D-3,  5.917159763 D-3,
     *      5.102040816 D-3,  4.444444445 D-3,
     *      3.90625     D-3 /
      DATA b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,b15
     * /1.00000000000000 D0, 1.50000000000000 D0,    
     *  1.83333333333333 D0, 2.08333333333333 D0,     
     *  2.28333333333333 D0, 2.45000000000000 D0,     
     *  2.59285714285714 D0, 2.71785714285714 D0,     
     *  2.82896825396825 D0, 2.92896825396825 D0,     
     *  3.01987734487734 D0, 3.10321067821068 D0,     
     *  3.18013375513376 D0, 3.25156232656233 D0,     
     *  3.31822899322899 D0 /    
      DATA pi / 3.141592654 D0 /
      DATA eu / 0.5772157 D0 /
c      
      if(cdabs(x).gt.9.d0)goto 10
c      
      a=1.d0
      z1=x*0.5d0           
      z=z1*z1           
      BJ0=a-z*(a-z*a1*(a-z*a2*(a-z*a3*(a-z*a4*(a-z*a5*(a-z*a6*
     *   (a-z*a7*(a-z*a8*(a-z*a9*(a-z*a10*(a-z*a11*(a-z*a12* 
     *   (a-z*a13*(a-z*a14*(a-z*a15)))))))))))))))
      BN0=z*(b1-z*a1*(b2-z*a2*(b3-z*a3*(b4-z*a4*(b5-z*a5*(b6-z*a6*
     *   (b7-z*a7*(b8-z*a8*(b9-z*a9*(b10-z*a10*(b11-z*a11*(b12-z*
     *   a12*(b13-z*a13*(b14-z*a14*(b15-z*a15)))))))))))))))
      BN0=2.d0*((cdlog(z1)+eu)*BJ0+BN0)/pi
      goto 20
c
   10 aa=cdsqrt(pi*x)
      co=cdcos(x)
      si=cdsin(x)
      x1=0.125d0/x
      x2=0.0703125d0/(x*x)
      x3=0.073242188d0/(x*x*x)
      BJ0=(co*(1.d0-x1-x2+x3)+si*(1.d0+x1-x2-x3))/aa 
      BN0=(si*(1.d0-x1-x2+x3)-co*(1.d0+x1-x2-x3))/aa 
c 
   20 return   
      end
c
c
c *** Bessel of the 0th and 1st orders ***
c
	subroutine BJ01 (X,BJ0,BJ1)
c
      IMPLICIT REAL*8 (A-H,O-Z)
      DATA a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15
     *   /  2.5         D-1,  1.111111111 D-1,
     *      6.25        D-2,  4.0         D-2,
     *      2.777777778 D-2,  2.040816327 D-2,
     *      1.5625      D-2,  1.234567901 D-2,
     *      1.0         D-2,  8.26446281  D-3,
     *      6.944444444 D-3,  5.917159763 D-3,
     *      5.102040816 D-3,  4.444444445 D-3,
     *      3.90625     D-3 /
      DATA b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,b15
     *   /  5.0         D-1,  1.666666667 D-1,
     *      8.333333333 D-2,  5.0         D-2,
     *      3.333333333 D-2,  2.380952381 D-2,
     *      1.785714286 D-2,  1.388888889 D-2,
     *      1.111111111 D-2,  9.090909091 D-3,
     *      7.575757576 D-3,  6.41025641  D-3,
     *      5.494505495 D-3,  4.761904762 D-3,
     *      4.166666666 D-3 /
      DATA pi / 3.141592654 D0 /
c
      if(dabs(x).gt.9.d0)goto 10
c
      a=1.d0
      z1=x*0.5d0
      z=z1*z1                       
      BJ0=a-z*(a-z*a1*(a-z*a2*(a-z*a3*(a-z*a4*(a-z*a5*(a-z*a6*
     *   (a-z*a7*(a-z*a8*(a-z*a9*(a-z*a10*(a-z*a11*(a-z*a12* 
     *   (a-z*a13*(a-z*a14*(a-z*a15)))))))))))))))
      BJ1=z1*(a-z*b1*(a-z*b2*(a-z*b3*(a-z*b4*(a-z*b5*(a-z*b6*
     *   (a-z*b7*(a-z*b8*(a-z*b9*(a-z*b10*(a-z*b11*(a-z*b12* 
     *   (a-z*b13*(a-z*b14*(a-z*b15)))))))))))))))
      goto 20
c
   10 aa=dsqrt(pi*x)
      co=dcos(x)
      si=dsin(x)
      x1=0.125d0/x
      x2=0.0703125d0/(x*x)
      x3=0.073242188d0/(x*x*x)
      BJ0=(co*(1.d0-x1-x2+x3)+si*(1.d0+x1-x2-x3))/aa 
      x1=0.375d0/x
      x2=0.1171875d0/(x*x)
      x3=0.102539063d0/(x*x*x)
      BJ1= (si*(1.d0+x1+x2-x3)-co*(1.d0-x1+x2+x3))/aa 
c 
   20 return
      end      
