C
C *** CALCULATION OF Q, V AND Q USING 
C     THE NEWTON'S (TRAPEZOIDAL)INTEGRAL FORMULA ***
C
      SUBROUTINE QVW(A,B,N,DNS,VP,VS,THK,ML,IER,L1,DEP1,L2,DEP2
     *              ,ISO,FQX,FQZ,FVX,FWX,FWZ)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VP(ML),VS(ML),CPX
     *          ,FQX,FQZ,FVX,FWX,FWZ
     *          ,QXY,QZY,VXY,WXY,WZY
     *          ,QXQ,QZQ,VXQ,WXQ,WZQ
     *          ,QXF,QZF,VXF,WXF,WZF
     *          ,WRR(2),URR(2),VLR
      DIMENSION DNS(ML),THK(ML)
      COMMON /FDAT/ OM,NL,RR
C
      IF(N.LT.1)THEN
        IER=1000
        FQX=(0.D0,0.D0)
        FQZ=(0.D0,0.D0)
        FVX=(0.D0,0.D0)
        FWX=(0.D0,0.D0)
        FWZ=(0.D0,0.D0)
        GO TO 9999
      END IF
      IF(A.EQ.B)THEN
        IER=2000
        FQX=(0.D0,0.D0)
        FQZ=(0.D0,0.D0)
        FVX=(0.D0,0.D0)
        FWX=(0.D0,0.D0)
        FWZ=(0.D0,0.D0)
        GO TO 9999
      END IF
C
      CPX=(0.D0,1.D0)
      IF(ISO.EQ.1)THEN
C
        H=(B-A)/N
        CALL UVWRL(A,DNS,VP,VS,THK,ML,L1,DEP1,L2,DEP2,WRR,URR,VLR)
        RK=A*RR
        BJ0=BESJ0(RK)
        BJ1=BESJ1(RK)
        BJD=BJ0-BJ1/RK
C
        QXY= URR(1)*BJD-CPX*VLR*BJ1/RK
        VXY=-URR(1)*BJ1/RK+CPX*VLR*BJD
        WXY=-CPX*WRR(1)*BJ1
        QZY=-URR(2)*BJ1
        WZY=-CPX*WRR(2)*BJ0
C
        QXQ=(0.D0,0.D0)
        VXQ=(0.D0,0.D0)
        WXQ=(0.D0,0.D0)
        QZQ=(0.D0,0.D0)
        WZQ=(0.D0,0.D0)
C
        DO  100 I=1,N
C
          X=A+H*I
          CALL UVWRL(X,DNS,VP,VS,THK,ML,L1,DEP1,L2,DEP2,WRR,URR,VLR)
          RK=X*RR
          BJ0=BESJ0(RK)
          BJ1=BESJ1(RK)
          BJD=BJ0-BJ1/RK
C
          QXF= URR(1)*BJD-CPX*VLR*BJ1/RK
          VXF=-URR(1)*BJ1/RK+CPX*VLR*BJD
          WXF=-CPX*WRR(1)*BJ1
          QZF=-URR(2)*BJ1
          WZF=-CPX*WRR(2)*BJ0
C
          QXQ= QXQ+H*(QXY+QXF)/2.D0
          VXQ= VXQ+H*(VXY+VXF)/2.D0
          WXQ= WXQ+H*(WXY+WXF)/2.D0
          QZQ= QZQ+H*(QZY+QZF)/2.D0
          WZQ= WZQ+H*(WZY+WZF)/2.D0
C
          QXY= QXF
          VXY= VXF
          WXY= WXF
          QZY= QZF
          WZY= WZF
C
  100   CONTINUE
C
      ELSE
C
        H=(B-A)/N
        CALL UVWRL(A,DNS,VP,VS,THK,ML,L1,DEP1,L2,DEP2,WRR,URR,VLR)
        RK=A*RR
        BJ0=BESJ0(RK)
        BJ1=BESJ1(RK)
        BJD=BJ0-BJ1/RK
C
        QXY= URR(1)*BJD-CPX*VLR*BJ1/RK
        VXY=-URR(1)*BJ1/RK+CPX*VLR*BJD
        WXY= URR(2)*BJ1
        QZY= CPX*WRR(1)*BJ1
        WZY=-CPX*WRR(2)*BJ0
C
        QXQ=(0.D0,0.D0)
        VXQ=(0.D0,0.D0)
        WXQ=(0.D0,0.D0)
        QZQ=(0.D0,0.D0)
        WZQ=(0.D0,0.D0)
C
        DO 200 I=1,N
C
          X=A+H*I
          CALL UVWRL(X,DNS,VP,VS,THK,ML,L1,DEP1,L2,DEP2,WRR,URR,VLR)
          RK=X*RR
          BJ0=BESJ0(RK)
          BJ1=BESJ1(RK)
          BJD=BJ0-BJ1/RK
C
          QXF= URR(1)*BJD-CPX*VLR*BJ1/RK
          VXF=-URR(1)*BJ1/RK+CPX*VLR*BJD
          WXF= URR(2)*BJ1
          QZF= CPX*WRR(1)*BJ1
          WZF=-CPX*WRR(2)*BJ0
C
          QXQ= QXQ+H*(QXY+QXF)/2.D0
          VXQ= VXQ+H*(VXY+VXF)/2.D0
          WXQ= WXQ+H*(WXY+WXF)/2.D0
          QZQ= QZQ+H*(QZY+QZF)/2.D0
          WZQ= WZQ+H*(WZY+WZF)/2.D0
C
          QXY= QXF
          VXY= VXF
          WXY= WXF
          QZY= QZF
          WZY= WZF
C
  200   CONTINUE
C
      END IF
C
      FQX= QXQ
      FVX= VXQ
      FWX= WXQ
      FQZ= QZQ
      FWZ= WZQ
C
      IER=0
C
 9999 RETURN
C
      END
