C     
c *** This code computes Green's functions from a point source. 
c     "phs1.f" is needed for the normal modes before running grpnt1.f.
c    (July, 1997)
C
c     References:
C     1) HISADA, Bull. Seismo. Soc. Am., V.84, Oct. 1994, (Part 1)
C     2) HISADA, Bull. Seismo. Soc. Am., V.85, Aug. 1995, (Part 2)
C     3) HISADA, Bull. Seismo. Soc. Am., V.85, Oct. 1995, (comments)
C
C *** FOR BUGS OR COMMENTS, PLEASE E-MAIL to "hisada@cc.kogakuin.ac.jp",
C     OR MAIL TO "Yoshiaki Hisada, Dept. of Architecture,
C     Kogakuin Univeristy, Nishi-Shinjuku 1-24-2, Tokyo 163-91, Japan"
C
C    (Originally made in 1992, and modified in 1995, 1996, and 1997)
C    (Modified from "grpoint.f" and "grflt92.f")
c
c * Parameters : 
c     MOM : Number of Omega
c     MOB : NUmber of Observation Points
c
      PROGRAM grpnt1
C
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(MOM=211,MOB=31)
      COMPLEX*16 GQ(MOM,MOB),GV(MOM,MOB),GW(MOM,MOB)
c     *          ,SQ(MOM,MOB),SV(MOM,MOB),SW(MOM,MOB)
C
      WRITE(*,*)
      WRITE(*,*)'"grpnt1.f" constructs Green functions'
      WRITE(*,*)
     *'           (Originally made in 1992, and modified in 1997)'
      write(*,*)
      WRITE(*,*)'"grpnt1.in" is the input data'
      write(*,*)
      WRITE(*,*)
     *' You must run "phs1.f" before running this code !!'
      WRITE(*,*)' Please Check "fort.10" for troubles.'
      WRITE(*,*)
     *' For questions, e-mail to "hisada@cc.kogakuin.ac.jp".'
      WRITE(*,*)' This code is open for the academic use.'
      WRITE(*,*)
      WRITE(*,*)'Output data are stored in '
      WRITE(*,*)'  "grpoint.dxyz" and "grpoint.dqvw" for displacement' 
c      WRITE(*,*)'  "grpoint.sxyz" and "grpoint.sqvw" for stress.' 
      WRITE(*,*)
      WRITE(*,*)'Please use "grfft.f" to obtain seismograms. '
C
      CALL INDAT(MOM,MOB)
      CALL INPHS(MOM)
c      CALL PNQVW(MOM,MOB,GQ,GV,GW,SQ,SV,SW)
      CALL PNQVW(MOM,MOB,GQ,GV,GW)
      CALL OTDAT(MOM,MOB,GQ,GV,GW,1)
c      CALL OTDAT(MOM,MOB,SQ,SV,SW,2)
C
      STOP
      END
C
c
C *** DATA INPUT ***
C
      SUBROUTINE INDAT(MOM,MOB)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VP,VS
      COMMON /FREQ/DOM,IOM,NOM /OBS/XO(51),YO(51),ZO(51),NOB
     *    /STRC/DNS(20),VP(20),VS(20),THK(20),NL
     *    /INTG/VPMAX,VSMAX,RYMIN,CFINAL,N11,N12,N13,N14,M13
     *    /SIGN/ISIGN /OTST/IOT,IST /SRCK/ZS,TP,TS,QX,QY,QZ,AA,IR 
C
C *** PARAMETERS IN COMMOM STATEMENT ***
C
      NML =20
      NMOB=51
      IF(MOB.LT.NMOB)NMOB=MOB
C
C *** OMEGA DATA ***
C
      OPEN (UNIT=7,FILE='grpnt1.in',STATUS='OLD')
C
      read(7,1100)dt,nt
      read(7,1110)period
 1100 format(/f8.0,i8)
 1110 format(f8.0)
      duration=dt*nt
      dom=datan(1.d0)*8.d0/duration
      iom=1
      nom=duration/period

      WRITE(10,*)'FROM indat:'
      WRITE(10,*)
      WRITE(10,1130)nt,dt,period
      WRITE(10,1140)NOM,DOM
 1130 format(/' Num. of Time, Delta Time and Min. Period (s)=',i5,2f8.4)
 1140 format( ' Num. of Omega, and Delta Omega (rad/sec) =',i5,f10.7)
C
      IF(NOM.GT.MOM)THEN
        WRITE(10,*)' NOM>MOM !!!!!!!!!!!!!'
        STOP
      END IF
C
C *** MEDIUM DATA ***
C
      READ(7,1200)NL
      DO 100 IL=1,NL
        READ(7,1210)DNS(IL),VPP,QP,VSS,QS,THK(IL)
        VP(IL)=VPP*DCMPLX(1.D0,-0.5D0/QP)
        VS(IL)=VSS*DCMPLX(1.D0,-0.5D0/QS)
  100 CONTINUE
 1200 FORMAT(/I10)
 1210 FORMAT(6E8.0)
C
      WRITE(10,2200)NL
      DO 200 IL=1,NL
        WRITE(10,2210)IL,DNS(IL),VP(IL),VS(IL),THK(IL)
  200 CONTINUE
 2200 FORMAT(//' Total Number of Layers =',I5/
     *       /' Layer, Density , Vp(real), Vp(imag)'
     *        ', Vs(real), Vs(imag),Thickness(m)')
 2210 FORMAT(I5,5F10.3,F12.3)
C
      IF(NL.GT.NML)THEN
        WRITE(10,*)' NL>NML !!!!!!!!!!!!!'
        STOP
      END IF
C
C *** MAX Vs and Vp and MIN Rayleigh Phase Velocity ***
C
      VSMAX=dreal(VS(NL))
      VPMAX=dreal(VP(NL))
      VSMIN=DREAL(VS(1))
      MINVS=1
      DO 220 IL=1,NL
        IF(dreal(VS(IL)).GT.VSMAX)VSMAX=dreal(VS(IL))
        IF(dreal(VP(IL)).GT.VPMAX)VPMAX=dreal(VP(IL))
        IF(dreal(VS(IL)).LT.VSMIN)THEN
          VSMIN=dreal(VS(IL))
          MINVS=IL
        END IF
  220 CONTINUE
      RYMIN=FRHF(dreal(VP(MINVS)),dreal(VS(MINVS)))
C
      WRITE(10,2250)VSMAX
      WRITE(10,2260)VPMAX
      WRITE(10,2270)RYMIN
 2250 FORMAT(/'    Max Vs (m/sec)              =',f10.3)
 2260 FORMAT( '    Max Vp (m/sec)              =',f10.3)
 2270 FORMAT( '    Min Rayleigh Phase Velocity =',f10.3)
C
C   * SOURCE POINT DATA *
C
      READ(7,1501)ZS
      WRITE(10,2501)ZS
 1501 FORMAT(/F10.0)
 2501 FORMAT(//' SP DEPTH =',F10.3)
C
      READ(7,1511)QX,QY,QZ
      WRITE(10,2512)QX,QY,QZ
 1511 FORMAT(3F10.0)
 2512 FORMAT(/' QX,QY AND QZ =',3F10.0)
C
      READ(7,1521)IR
C
      IF(IR.EQ.2)THEN
          READ(7,*)
          READ(7,*)
      ELSE
          READ(7,1531)TP,TS
          WRITE(10,2521)IR
          WRITE(10,2531)TP,TS
          AA=DNS(1)*VS(1)*VS(1)
          WRITE(10,2541)AA
          READ(7,1541)POW
          WRITE(10,2542)POW
          AA=AA**POW
          WRITE(10,2543)AA
      END IF
C
 1521 FORMAT(I10)
 1531 FORMAT(2F10.0)
 2521 FORMAT(/' RICKER WAVELET (=1) OR ORIGINAL WAVE (=2) :',I5)
 2531 FORMAT( ' TP (DOMINANT PERIOD) AND TS (DELAY TIME) FOR'
     *  ' RICKER WAVELET =',2F10.4)
 1541 FORMAT(F10.0)
 2541 FORMAT(/' original Amplification Factor (AA) = ',2D12.4)
 2542 FORMAT( ' AA**N; N=: ',F10.5)
 2543 FORMAT(/' New Amplification Factor  = ',2D12.4)
C
C *** NUMERICAL INTEGRAL DATA ***
C
      READ(7,1400)N11
      READ(7,1410)N12
      READ(7,1410)N13
      READ(7,1410)N14
      read(7,1420)vfac
      CFINAL=rymin/vfac
 1400 FORMAT(/I10)
 1410 FORMAT(I10)
 1420 FORMAT(F10.0)
C
      WRITE(10,2400)
      WRITE(10,2410)N11
      WRITE(10,2420)N12
      WRITE(10,2430)N13
      WRITE(10,2440)N14
      WRITE(10,2450)vfac,CFINAL
      WRITE(10,2460)2*(N11+N12+N13+N14)+1
 2400 FORMAT(//' Number of Points for Numerical Integrations:')
 2410 FORMAT( /'   From 0               to 0.8*om/Vp(max)   =',I5)
 2420 FORMAT(  '   From 0.8*om/Vp(max)  to 0.95*om/Vs(max)  =',I5)
 2430 FORMAT(  '   From 0.95*om/Vs(max) to 1.05*om/Ray(min) =',I5)
 2440 FORMAT(  '   From 1.05*om/Ry(min) to om/C(final)      =',I5)
 2450 FORMAT(  '   VFAC(Factor for Minimum Phase Velocity)  =',F10.3,
     *        /'   C(final = Min. Rayleigh Phase  Vel./VFAC)=',F10.3)
 2460 FORMAT(  '   Initial Number of Points: 2(N1+N2+N3+N4)+1 =',I5)
C
C *** DATA FOR THE ASYMPTOTIC ***
C
      READ(7,1600)IST
      IF(IST.EQ.1)THEN
        WRITE(10,2600)
      ELSE
        WRITE(10,2610)
      END IF
 1600 FORMAT(/I10)
 2600 FORMAT(//' THE ASYMPTOTIC METHOD IS USED')
 2610 FORMAT(//' THE ASYMPTOTIC METHOD IS NOT USED')
C
C *** CHANGE OF SIGNS FOR IMAGINARY PARTS ***
C
      READ(7,1700)ISIGN
      IF(ISIGN.EQ.1)WRITE(10,2700)
      IF(ISIGN.EQ.0)WRITE(10,2710)
 1700 FORMAT(/I10)
 2700 FORMAT(/' SIGNS OF IMAGINARY PARTS FOR RESULTS ARE CHANGED')
 2710 FORMAT(/' SIGNS OF IMAGINARY PARTS FOR RESULTS ARE NOT CHANGED')
C
C *** OP DATA ***
C
      READ(7,1500)NOB
      WRITE(10,2500)NOB
 1500 FORMAT(/I10)
 2500 FORMAT(//' TOTAL NUMBER OF OBSERVATION POINTS =',I6/)
C
      IF(NOB.GT.NMOB)THEN
        WRITE(10,*)' NOB>NMOB !!!!!!!!!!!!!'
        STOP
      END IF
C
      DO 300 I=1,NOB
        READ(7,1520)XO(I),YO(I),ZO(I)
        WRITE(10,2520)I,XO(I),YO(I),ZO(I)
  300 CONTINUE
 1520 FORMAT(3F10.0)
 2520 FORMAT('   OB NUM.=',I5,',     (X,Y,Z) =',3F12.4)
C
      CLOSE(UNIT=7)
C
C *** OUTPUT OF INTEGRAND OF U,V,AND W ***
C
   50 IOT=0
      IF(NOB.NE.1)RETURN
      IF((NOM-IOM).NE.0)RETURN
C
      WRITE(*,7000)
      READ(*,*)IOT
 7000 FORMAT(/' Output Integrands?: Enter No=0 or Yes=1: ')
C
      if(iot.eq.1)then
        write(*,*)
        write(*,*)'  You can plot the integrands using "integ.m",'
        write(*,*)'  "integr.m", and "integz.m" for "MATLAB"'
c
        open(9,file='in_dat',status='unknown')
          write(9,9000)dom*iom
          write(9,9100)nl
          do 500 i=1,nl
            write(9,9000)dreal(vp(i))
            write(9,9000)dreal(vs(i))
 500      continue
        close(9)
 9000   format(e15.7)
 9100   format(i5)
c
      end if
C
      RETURN
      END
C
C
C *** PHASE VELOCITY OF THE RAYLEIGH WAVE IN THE HOMOGENEOUS HALF-SPACE
C
      FUNCTION FRHF(VP,VS)
      IMPLICIT REAL*8(A-H,O-Z)
C
      DZ=0.0001D0
C
      SP=VS/VP
      IF(SP.GE.0.0D0) Z= 0.91263D0
      IF(SP.GE.0.2D0) Z= 0.90767D0
      IF(SP.GE.0.4D0) Z= 0.88899D0
      IF(SP.GE.0.6D0) Z= 0.83575D0
      IF(SP.GE.0.7D0) Z= 0.77062D0
      IF(SP.GE.0.8D0) Z= 0.63338D0
C
      F0=FRHF1(SP,Z)
  100 Z=Z-DZ
      F1=FRHF1(SP,Z)
      AF=F0/F1
      F00=F0
      F0 =F1
      IF(AF.GE.0.D0)GOTO 100
      FZ=Z+DZ*F0/(F0-F00)
C
      FRHF=VS*DSQRT(FZ)
C
      RETURN
      END
C
C
      FUNCTION FRHF1(SP,Z)
      IMPLICIT REAL*8(A-H,O-Z)
C
      Z2=Z*Z
      Z3=Z2*Z
C
      F1=16.D0*(Z-1.D0)*SP*SP
      FRHF1=F1-Z3+8.D0*Z2-24.D0*Z+16.D0
C
      RETURN
      END
C
C
C *** DATA OUTPUT ***
C
      SUBROUTINE OTDAT(MOM,MOB,GQ,GV,GW,IDS)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VP,VS,GQ(MOM,MOB),GV(MOM,MOB),GW(MOM,MOB),AA,CC
      COMMON /FREQ/DOM,IOM,NOM /SIGN/ISIGN
     *    /STRC/DNS(20),VP(20),VS(20),THK(20),NL
     *    /OBS/XO(51),YO(51),ZO(51),NOB
      CHARACTER FNAME*15
C
C   * OUTPUT OF DATA (GQ, GV AND GW: CYLINDRICAL COORDINATE) *
C
      IF(ISIGN.EQ.1)THEN
        DO 100 IOB=1,NOB
        DO 110 K=IOM,NOM
          GQ(K,IOB)=DCMPLX(DREAL(GQ(K,IOB)),-DIMAG(GQ(K,IOB)))
          GV(K,IOB)=DCMPLX(DREAL(GV(K,IOB)),-DIMAG(GV(K,IOB)))
          GW(K,IOB)=DCMPLX(DREAL(GW(K,IOB)),-DIMAG(GW(K,IOB)))
  110   CONTINUE
  100   CONTINUE
      END IF
C
      FNAME='grpoint.dqvw'
      IF(IDS.EQ.2)FNAME='grpoint.sqvw'
C
      OPEN(8,FILE=FNAME,STATUS='unknown')
        WRITE(8,2500)DOM,IOM,NOM
        WRITE(8,2510)NOB
        DO 500 IOB=1,NOB
          DO 510 K=IOM,NOM
            WRITE(8,2520)GQ(K,IOB),GV(K,IOB),GW(K,IOB)
  510     CONTINUE
  500   CONTINUE
      CLOSE(8)
C
 2500 FORMAT(1X,F10.6,2I10)
 2510 FORMAT(1X,I10,///)
 2520 FORMAT(1X,6E12.5)
C
C   * OUTPUT OF DATA (GX, GY AND GZ: CARTESIAN COORDINATE) *
C
        DO 600 IOB=1,NOB
          RR=DSQRT(XO(IOB)*XO(IOB)+YO(IOB)*YO(IOB))
          CR=XO(IOB)/RR
          SR=YO(IOB)/RR
          DO 610 K=IOM,NOM
            AA=GQ(K,IOB)
            CC=GV(K,IOB)
            GQ(K,IOB)=AA*CR-CC*SR
            GV(K,IOB)=AA*SR+CC*CR
  610     CONTINUE
  600   CONTINUE
C
      FNAME='grpoint.dxyz'
      IF(IDS.EQ.2)FNAME='grpoint.sxyz'
C      
      OPEN(8,FILE=FNAME,STATUS='unknown')
        WRITE(8,3500)DOM,IOM,NOM
        WRITE(8,3510)NOB
        DO 800 IOB=1,NOB
          DO 810 K=IOM,NOM
            WRITE(8,3520)GQ(K,IOB),GV(K,IOB),GW(K,IOB)
  810     CONTINUE
  800   CONTINUE
      CLOSE(8)
C
 3500 FORMAT(1X,F10.6,2I10)
 3510 FORMAT(1X,I10,///)
 3520 FORMAT(1X,6E12.5)
C
      RETURN
      END
c
c
C *** INPUT OF PHASE VELOCITIES AND MAKE POLE DATA ***
C
      SUBROUTINE INPHS(MOM)
C
      IMPLICIT REAL*8(A-H,O-Z)
      dimension CR(211,0:201),CL(211,0:201),MDL(211),MDR(211)
      COMPLEX*16 VP,VS
      COMMON /FREQ/DOM,IOM,NOM /POLE/CS(211,201),NW(211)
     *       /STRC/DNS(20),VP(20),VS(20),THK(20),NL
C
C  ** CHECK ** 
C
      NC=211
      IF(MOM.GT.NC)THEN
         WRITE(10,*)'FROM INPHS: MOM>NC !!'
         STOP
      END IF
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(10,*)'FROM INPHS: NL>NC !!'
         STOP
      END IF

c   * Number of Poles *
      NC=201
c
C  ** INPUT OF PHASE VELOCITIES FOR LOVE AND RAYLEIGH MODES **
C
      open(22,file='mdl_dat',status='old')
      open(23,file='mdr_dat',status='old')
      DO 100 IO=IOM,NOM
        READ(22,*)MDL(IO)
        READ(23,*)MDR(IO)
        IF(MDL(IO).GT.NC)THEN
         WRITE(10,*)'FROM INPHS: MDL>NC !!'
         STOP
        END IF
        IF(MDR(IO).GT.NC)THEN
         WRITE(10,*)'FROM INPHS: MDR>NC !!'
         STOP
        END IF
  100 CONTINUE
      close(22)
      close(23)

      open(20,file='phsl_dat',status='old')
      DO 110 I=0,MDL(NOM)
      DO 110 IO=IOM,NOM
        read(20,*)CL(IO,I)
  110 CONTINUE
      close(20)

      open(21,file='phsr_dat',status='old')
      DO 120 I=0,MDR(NOM)
      DO 120 IO=IOM,NOM
        read(21,*)CR(IO,I)
  120 CONTINUE
      close(21)
C
C  ** CONSTRUCT PHASE VELOCITY DATA **
C
      write(10,1000)
 1000 FORMAT(//' From INPHS:')

      DO 200 IO=IOM,NOM

        IF(NL.EQ.1)THEN
          N=1
          CS(IO,N)=CR(IO,0)
          GOTO 210
        END IF
C
        N=0
        IR=0
        IL=0

   30   N=N+1
        if((CR(IO,IR).lt.CL(IO,IL)).and.(CR(IO,IR).ne.0.d0))then
          CS(IO,N)=CR(IO,IR)
          IR=IR+1
        else
          IF((CL(IO,IL).EQ.0.D0).AND.(CR(IO,IR).NE.0.D0))THEN
            CS(IO,N)=CR(IO,IR)
            IR=IR+1
          ELSE
            if(CL(IO,IL).ne.0.d0)then
              CS(IO,N)=CL(IO,IL)
              IL=IL+1
            else
              N=N-1
              GOTO 210
            end if
          end if
        end if

        if((IR.le.MDR(IO)).or.(IL.le.MDL(IO)))goto 30


 210    NW(IO)=N
c
        NC=201
        IF(N.GE.NC)THEN
          WRITE(10,*)' FROM INPHS: NW.GE.NC !!'
          STOP
        END IF
C
        DO 220 I=1,NW(IO)-1
          IF(CS(IO,I).EQ.0.D0)THEN
            WRITE(10,*)' FROM INPHS: CS=0 !!'
            STOP
          END IF
          IF(CS(IO,I).GT.CS(IO,I+1))THEN
            WRITE(10,*)' FROM INPHS: CS(I)>CS(I+1) !!'
            STOP
          END IF
  220   CONTINUE
        IF(CS(IO,N).EQ.0.D0)THEN
          WRITE(10,*)' FROM INPHS: CS=0 !!'
          STOP
        END IF
c
        write(10,1100)io*dom,nw(io)
        write(10,1110)
        do 230 i=1,NW(IO),5
          write(10,1120)CS(IO,i),CS(IO,i+1)
     *               ,CS(IO,i+2),CS(IO,i+3),CS(IO,i+4)
  230   continue
 1100   FORMAT('   Omega, Number of Phase Velocities=',f8.5,I5)
 1110   FORMAT( '   Phase Velocities (m/sec) =')
 1120   FORMAT(2x,5D12.4)

  200 continue
C
      RETURN
      END
C
C
C *** COMPUTATION OF GREEN FUNCTIONS USING THE R/T METHOD ***
c *** For large wavenumber * distance (k*r > 10), the Filon's integration 
c     formula is used instead of the Simpson's.
C
c      SUBROUTINE PNQVW(MOM,MOB,GQ,GV,GW,SQ,SV,SW)
      SUBROUTINE PNQVW(MOM,MOB,GQ,GV,GW)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VP,VS,SM,RAM,SM1,XX1,SM2,XX2,SM3,XX3,RM1,RM2,RM3
     *      ,GQ(MOM,MOB),GV(MOM,MOB),GW(MOM,MOB),CC
c     *      ,SQ(MOM,MOB),SV(MOM,MOB),SW(MOM,MOB),VS2,VP2
     *      ,QXG(0:2),QZG(0:2),VXG(0:2),WXG(0:2),WZG(0:2)
     *      ,QXS(0:2),QZS(0:2),VXS(0:2),WXS(0:2),WZS(0:2)
     *      ,QA,QB,VA,VB,WA
     *      ,QX(0:3),QY(0:3),QZ(0:3),VX(0:3),VY(0:3)
     *      ,WX(0:3),WY(0:3),WZ(0:3),UX(3,3),UY(3,3),UZ(3,3)
      COMMON /FREQ/DOM,IOM,NOM /OBS/XO(51),YO(51),ZO(51),NOB
     *    /STRC/DNS(20),VP(20),VS(20),THK(20),NL
     *    /OTST/IOT,IST  /FCQVW/SM,RAM,RR  /SOLY/SDEP,ODEP,LS,LO
     *    /ASYMP/SM1,XX1,SM2,XX2,SM3,XX3 /RAM/RM1,RM2,RM3
     *    /DPSO /DEP1,DEP2,ZS,ZZO  /SRCK/ZZS,TP,TS,QQX,QQY,QQZ,GAA,IR 
     *    /INTG/VPMAX,VSMAX,RYMIN,CFINAL,N11,N12,N13,N14,M13
     *    /WVNUM/WVN(2001),DWN(4),WNR(5),NWN
C
C *** CHECK OF DIMENSION ***
C
      NC=20
      IF(NL.GE.NC)THEN
        WRITE(10,*)' FROM HKQVW:NL>NC !!!!!!!!!!'
        WRITE(10,*)' NL,NC=',NL,NC
        STOP
      END IF
C
      NC=51
      IF(NOB.GE.NC)THEN
        WRITE(10,*)' FROM HKQVW:NOB>NC !!!!!!!!!!'
        WRITE(10,*)' NOB,NC=',NOB,NC
        STOP
      END IF
C
      NC=MOB
      IF(NOB.GE.NC)THEN
        WRITE(10,*)' FROM HKQVW:NOB>MOB !!!!!!!!!!'
        WRITE(10,*)' NOB,NC=',NOB,MOB
        STOP
      END IF
C
C *** GREEN'S FUNCTION ***
C
C   * output of integrands for check *
C
      IF(IOT.eq.1)THEN
      OPEN (19,FILE='wn_dat',STATUS='unknown')
      OPEN (20,FILE='qxr_dat',STATUS='unknown')
      OPEN (21,FILE='qxi_dat',STATUS='unknown')
      OPEN (22,FILE='vxr_dat',STATUS='unknown')
      OPEN (23,FILE='vxi_dat',STATUS='unknown')
      OPEN (24,FILE='wxr_dat',STATUS='unknown')
      OPEN (25,FILE='wxi_dat',STATUS='unknown')
      OPEN (26,FILE='qzr_dat',STATUS='unknown')
      OPEN (27,FILE='qzi_dat',STATUS='unknown')
      OPEN (28,FILE='wzr_dat',STATUS='unknown')
      OPEN (29,FILE='wzi_dat',STATUS='unknown')
      OPEN (30,FILE='qxrr_dat',STATUS='unknown')
      OPEN (31,FILE='qxri_dat',STATUS='unknown')
      OPEN (32,FILE='vxrr_dat',STATUS='unknown')
      OPEN (33,FILE='vxri_dat',STATUS='unknown')
      OPEN (34,FILE='wxrr_dat',STATUS='unknown')
      OPEN (35,FILE='wxri_dat',STATUS='unknown')
      OPEN (36,FILE='qzrr_dat',STATUS='unknown')
      OPEN (37,FILE='qzri_dat',STATUS='unknown')
      OPEN (38,FILE='wzrr_dat',STATUS='unknown')
      OPEN (39,FILE='wzri_dat',STATUS='unknown')
      OPEN (40,FILE='qxzr_dat',STATUS='unknown')
      OPEN (41,FILE='qxzi_dat',STATUS='unknown')
      OPEN (42,FILE='vxzr_dat',STATUS='unknown')
      OPEN (43,FILE='vxzi_dat',STATUS='unknown')
      OPEN (44,FILE='wxzr_dat',STATUS='unknown')
      OPEN (45,FILE='wxzi_dat',STATUS='unknown')
      OPEN (46,FILE='qzzr_dat',STATUS='unknown')
      OPEN (47,FILE='qzzi_dat',STATUS='unknown')
      OPEN (48,FILE='wzzr_dat',STATUS='unknown')
      OPEN (49,FILE='wzzi_dat',STATUS='unknown')
      END IF
C
      BB=TP*TP*0.0253302D0
C
c   * Data for the source layer *
c
      ZS=ZZS
      DRU=0.0
      DRL=THK(1)
      DO 220 L=1,NL
         IF((ZS.GE.DRU).AND.(ZS.LT.DRL))LS=L
         IF((ZS.GE.DRU).AND.(L.EQ.NL))LS=L
         DRU=DRL
         DRL=DRL+THK(L+1)
 220  CONTINUE
C
      DEP1=0.0
      DO 230 L=1,LS-1
         DEP1=DEP1+THK(L)
 230  CONTINUE
      SDEP=ZS-DEP1
      DEP2=DEP1+THK(LS)
C
      SM=DNS(LS)*VS(LS)*VS(LS)
      RAM=DNS(LS)*VP(LS)*VP(LS)-2.D0*SM
C
      VS2= VS(LS)*VS(LS)
      SM2= DNS(LS)*VS2
      VP2= VP(LS)*VP(LS)
      RM2= DNS(LS)*VP2-2.0*SM2
      XX2 = (1.0+VS2/VP2)/(1.0-VS2/VP2)
      IF(LS.GE.2)THEN
         VS2= VS(LS-1)*VS(LS-1)
         SM1= DNS(LS-1)*VS2
         VP2= VP(LS-1)*VP(LS-1)
         RM1= DNS(LS-1)*VP2-2.0*SM1
         XX1 = (1.0+VS2/VP2)/(1.0-VS2/VP2)
      END IF
      IF(LS.NE.NL)THEN
         VS2= VS(LS+1)*VS(LS+1)
         SM3= DNS(LS+1)*VS2
         VP2= VP(LS+1)*VP(LS+1)
         RM3= DNS(LS+1)*VP2-2.0*SM3
         XX3 = (1.0+VS2/VP2)/(1.0-VS2/VP2)
      END IF
C
c   * Green's function computations *
c
      DO 100 IO=IOM,NOM
C
        OM=DOM*DFLOAT(IO)
C
        IF(IR.EQ.1)THEN
            OM2=OM*OM
            CC=GAA*OM2*DEXP(-OM2*BB)*CDEXP(DCMPLX(0.D0,+OM*TS))
        ELSE
            CC=(1.D0,0.0)
        END IF
C
c     * Configure wavenumbers for numerical integrations *
c
c       * range 1 (from near 0 to near om/Vp(max) *
        WNR(2)=om/vpmax*0.8
        WNR(1)=WNR(2)/n11/10000.d0
        DWN(1)=(WNR(2)-WNR(1))/(2.D0*N11)
        DO 120 K=1,2*N11
          WVN(K)=WNR(1)+DWN(1)*(K-1)
 120    CONTINUE

c       * range 2 (from near om/Vp(max) to near om/Vs(max) *
        WNR(3)=om/vsmax*0.95
        DWN(2)=(WNR(3)-WNR(2))/(2.D0*N12)
        DO 130 K=1,2*N12
          WVN(2*N11+K)=WNR(2)+DWN(2)*(K-1)
 130    CONTINUE
C
c       * range 3 (from near om/Vs(max) to near om/Rayleigh(min) *
        WNR(4)=om/rymin*1.05 
        CALL WVNUMV(MOM,IO,OM)
c
        NWN=2*(N11+N12+M13+N14)+1
        NC=2001
        IF(NWN.GT.NC)THEN
          WRITE(10,*)'FROM HKQVW: NWN > 2001 !!!'
        STOP
        END IF
C
c       * range 4 (from near om/Rayleigh(min) to near om/C(final) *
        WNR(5)=om/cfinal
        DWN(4)=(WNR(5)-WNR(4))/(2.D0*N14)
        DO 150 K=1,2*N14+1
          WVN(2*(N11+N12+M13)+K)=WNR(4)+DWN(4)*(K-1)
 150    CONTINUE
C
        write(*,1005)om,nwn
c        write(10,*)
c        write(10,1005)om,nwn
c        write(10,1010)wnr(1),wnr(2),wnr(3),wnr(4),wnr(5)
c        write(10,1020)dwn(1),dwn(2),dwn(3),dwn(4)
 1005   format('    Omega, Num. of integration points =',f10.4,I6)
c 1010   format('    Boundaries for the Integration ranges:'/x,5e12.4)
c 1020   format('    Integration Increments for the ranges:'/x,4e12.4)
C
c     * Configure the up/down-going coefficient library *
c
        CALL CUDSH(OM)
        CALL CUDPSV(OM)
c
c     * fault points and Green's function computation  *
c
        DO 300 IOB=1,NOB
C
          RX=XO(IOB)
          RY=YO(IOB)
          RR=DSQRT(RX*RX+RY*RY)
          CR=RX/RR
          SR=RY/RR
C
          DRU=0.0
          DRL=THK(1)
          DO 320 L=1,NL
            IF((ZO(IOB).GE.DRU).AND.(ZO(IOB).LT.DRL))LO=L
            IF((ZO(IOB).GE.DRU).AND.(L.EQ.NL))LO=L
            DRU=DRL
            DRL=DRL+THK(L+1)
 320      CONTINUE
C
          ZZO=ZO(IOB)
          ODEP=ZZO
          DO 330 L=2,LO
            ODEP=ODEP-THK(L-1)
 330      CONTINUE
c
c       * coefficients of the asymptotic solutions *
          IF((IST.EQ.1).AND.(IABS(LS-LO).LE.1))THEN
             CALL CASYM(LS,NL)
          END IF
C
          DO J=0,2
            QXG(J)=(0.D0,0.D0)
            VXG(J)=(0.D0,0.D0)
            WXG(J)=(0.D0,0.D0)
            QZG(J)=(0.D0,0.D0)
            WZG(J)=(0.D0,0.D0)
            QXS(J)=(0.D0,0.D0)
            VXS(J)=(0.D0,0.D0)
            WXS(J)=(0.D0,0.D0)
            QZS(J)=(0.D0,0.D0)
            WZS(J)=(0.D0,0.D0)
    	  end do
c
c    * the wavenumber integration using Simpson's and Filon's fomulas *
          CALL QVWRZS(OM,QXG,QZG,VXG,WXG,WZG)
          CALL QVWRZF(OM,QXG,QZG,VXG,WXG,WZG)
C
c    * the analytical wavenumber integration of the asymptotic solutions *
          IF((IST.EQ.1).AND.(IABS(LS-LO).LE.1))THEN
            CALL AINTRZ(LS,NL,QXS,QZS,VXS,WXS,WZS)
          END IF
C       
c     * Urx, Ury, Urz: see Hisada(1994) eq.(11) *
          QA=QXG(0)+QXS(0)
            QX(0)=QA*CR
            QY(0)=QA*SR
            QZ(0)=-(QZG(0)+QZS(0))
C
c     * Utx, Uty     : see Hisada(1994) eq.(11) *
          VA=VXG(0)+VXS(0)
            VX(0)=-VA*SR
            VY(0)= VA*CR
C	       
c     * Uzx, Uzy, Uzz: see Hisada(1994) eq.(11) *
            WA=WXG(0)+WXS(0)
              WX(0)=-WA*CR
              WY(0)=-WA*SR
              WZ(0)=-(WZG(0)+WZS(0))
c
c    * see Hisada(1994) eq.(26) *
                C2=CR*CR
                S2=SR*SR
                CS=CR*SR
                CX=S2/RR
                SY=C2/RR
                CY=-CS/RR
                SX=CY
c
c     * Urx,x, Ury,x, Urx,y, Ury,y, Urz,x Urz,y: see Hisada(1994) eq.(27) *
                QB=QXG(1)+QXS(1)
                  QX(1)=QB*C2+QA*CX
                  QY(1)=QB*CS+QA*SX
                  QX(2)=QB*CS+QA*CY
                  QY(2)=QB*S2+QA*SY
                  QB=QZG(1)+QZS(1)
                  QZ(1)=-QB*CR
                  QZ(2)=-QB*SR
c
c     * Utx,x, Uty,x, Utx,y, Uty,y: see Hisada(1994) eq.(27) *
                VB=VXG(1)+VXS(1)
                  VX(1)=-(VB*CS+VA*SX)
                  VY(1)=  VB*C2+VA*CX
                  VX(2)=-(VB*S2+VA*SY)
                  VY(2)=  VB*CS+VA*CY
c    
c     * Uzx,x, Uzy,x, Uzx,y, Uzy,y, Uzz,x Uzz,y: see Hisada(1994) eq.(27) 
                WB=WXG(1)+WXS(1)
                  WX(1)=-(WB*C2+WA*CX)
                  WY(1)=-(WB*CS+WA*SX)
                  WX(2)=-(WB*CS+WA*CY)
                  WY(2)=-(WB*S2+WA*SY)
                  WB=WZG(1)+WZS(1)
                  WZ(1)=-WB*CR
                  WZ(2)=-WB*SR
c
c     * Urx,z, Ury,z, Urz,z: see Hisada(1994) eq.(29) *
                QB=QXG(2)+QXS(2)
                  QX(3)= QB*CR
                  QY(3)= QB*SR
                  QZ(3)=-(QZG(2)+QZS(2))
c		 
c     * Utx,z, Uty,z: see Hisada(1994) eq.(29) *
                VB=VXG(2)+VXS(2)
                  VX(3)=-VB*SR
                  VY(3)= VB*CR
c	     
c     * Uzx,z, Uzy,z, Uzz,z: see Hisada(1994) eq.(29) *
                WB=WXG(2)+WXS(2)
                  WX(3)=-WB*CR
                  WY(3)=-WB*SR
                  WZ(3)=-(WZG(2)+WZS(2))
C		 
c     * Uxx,x, Uxx,y, Uxy,x, Uxy,y, Uxz,x, Uxz,y: see Hisada(1994) eq.(24) *
                UX(1,1)=QX(1)*CR+QX(0)*CX-VX(1)*SR-VX(0)*SX
                UX(1,2)=QX(2)*CR+QX(0)*CY-VX(2)*SR-VX(0)*SY
                UY(1,1)=QY(1)*CR+QY(0)*CX-VY(1)*SR-VY(0)*SX
                UY(1,2)=QY(2)*CR+QY(0)*CY-VY(2)*SR-VY(0)*SY
                UZ(1,1)=QZ(1)*CR+QZ(0)*CX
                UZ(1,2)=QZ(2)*CR+QZ(0)*CY
C	 
c     * Uyx,x, Uyx,y, Uyy,x, Uyy,y, Uyz,x, Uyz,y: see Hisada(1994) eq.(24) *
                UX(2,1)=QX(1)*SR+QX(0)*SX+VX(1)*CR+VX(0)*CX
                UX(2,2)=QX(2)*SR+QX(0)*SY+VX(2)*CR+VX(0)*CY
                UY(2,1)=QY(1)*SR+QY(0)*SX+VY(1)*CR+VY(0)*CX
                UY(2,2)=QY(2)*SR+QY(0)*SY+VY(2)*CR+VY(0)*CY
                UZ(2,1)=QZ(1)*SR+QZ(0)*SX
                UZ(2,2)=QZ(2)*SR+QZ(0)*SY
C
c     * Uzx,x, Uzx,y, Uzy,x, Uzy,y, Uzz,x, Uzz,y: see Hisada(1994) eq.(27) *
                UX(3,1)=WX(1)
                UX(3,2)=WX(2)
                UY(3,1)=WY(1)
                UY(3,2)=WY(2)
                UZ(3,1)=WZ(1)
                UZ(3,2)=WZ(2)
C
c     * Uxx,z, Uxy,z, Uxz,z, Uyx,z, Uyy,z, Uyz,z: see Hisada(1994) eq.(25) *
                UX(1,3)=QX(3)*CR-VX(3)*SR
                UY(1,3)=QY(3)*CR-VY(3)*SR
                UZ(1,3)=QZ(3)*CR
                UX(2,3)=QX(3)*SR+VX(3)*CR
                UY(2,3)=QY(3)*SR+VY(3)*CR
                UZ(2,3)=QZ(3)*SR
C
c     * Uzx,z, Uzy,z, Uzz,z: see Hisada(1994) eq.(29) *
                UX(3,3)=WX(3)
                UY(3,3)=WY(3)
                UZ(3,3)=WZ(3)
C
c *** Displacement Green's functions for point sources ***
c
          GQ(IO,IOB)=CC*(QX(0)*QQX+QY(0)*QQY+QZ(0)*QQZ)
          GV(IO,IOB)=CC*(VX(0)*QQX+VY(0)*QQY)
          GW(IO,IOB)=CC*(WX(0)*QQX+WY(0)*QQY+WZ(0)*QQZ)
C
  300   CONTINUE
C
  100 CONTINUE
C
      IF(IOT.eq.1)THEN
      CLOSE(19)
      CLOSE(20)
      CLOSE(21)
      CLOSE(22)
      CLOSE(23)
      CLOSE(24)
      CLOSE(25)
      CLOSE(26)
      CLOSE(27)
      CLOSE(28)
      CLOSE(29)
      CLOSE(30)
      CLOSE(31)
      CLOSE(32)
      CLOSE(33)
      CLOSE(34)
      CLOSE(35)
      CLOSE(36)
      CLOSE(37)
      CLOSE(38)
      CLOSE(39)
      CLOSE(40)
      CLOSE(41)
      CLOSE(42)
      CLOSE(43)
      CLOSE(44)
      CLOSE(45)
      CLOSE(46)
      CLOSE(47)
      CLOSE(48)
      CLOSE(49)
      END IF
C
      RETURN
      END
C
c
C *** CONSTRUCT POINTS FOR WAVENUMBER INTEGRATION ***
C
      SUBROUTINE WVNUMV(MOM,IO,OM)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON /INTG/VPMAX,VSMAX,RYMIN,CFINAL,N11,N12,N13,N14,M13
     *       /WVNUM/WVN(2001),DWN(4),WNR(5),NWN
     *       /POLE/CS(211,201),NW(211) /TMP/WVT(1001)
C
C  ** CHECK ** 
C
      NC=211
      IF(MOM.GT.NC)THEN
         WRITE(10,*)'FROM WVNUM: MOM > NC !!'
         STOP
      END IF
C
      NC=202
      IF(NW(IO).GT.NC)THEN
         WRITE(10,*)'FROM WVNUM: NW(IO) > NC !!'
         STOP
      END IF
C
      NC=1001
      IF((N13+1+NW(IO)*9).GT.NC)THEN
         WRITE(10,*)'FROM WVNUM: N13+1+NW(IO)*9 > NC !!'
         STOP
      END IF
C
C  ** POINTS **
C
      DWN(3)=(WNR(4)-WNR(3))/N13
      DO 100 K=1,N13+1
        WVT(K)=WNR(3)+DWN(3)*(K-1)
 100  CONTINUE

      MW=N13
      DW=DWN(3)/16.D0
      DO 200 J=1,NW(IO)
        POLE=OM/CS(IO,J)
        IPL=0
        DO 210 K=1,NW(IO)
          PL1=OM/CS(IO,K)
          DP1=PL1-POLE
          IF((DABS(DP1).LT.DW).AND.(DP1.GT.0.D0))IPL=1
 210    CONTINUE

        IF(IPL.EQ.0)THEN
          W1=POLE-DWN(3)/2.D0
          CALL WVNM(W1,MW)
          W2=POLE-DWN(3)/4.D0
          CALL WVNM(W2,MW)
          W3=POLE-DWN(3)/8.D0
          CALL WVNM(W3,MW)
          W4=POLE-DWN(3)/16.D0
          CALL WVNM(W4,MW)
          W5=POLE
          CALL WVNM(W5,MW)
          W6=POLE+DWN(3)/16.D0
          CALL WVNM(W6,MW)
          W7=POLE+DWN(3)/8.D0
          CALL WVNM(W7,MW)
          W8=POLE+DWN(3)/4.D0
          CALL WVNM(W8,MW)
          W9=POLE+DWN(3)/2.D0
          CALL WVNM(W9,MW)
        ELSE
          CALL WVNM(POLE,MW)
        END IF
 200  CONTINUE

      DW=DWN(3)/16.D0*0.9999D0
      DR=DW/1000.D0
      DO 300 I=1,MW
        WA=WVT(I)
        WB=WVT(I+1)
        WD=WB-WA
        ISML=0
        IF(WD.LT.DW)THEN
          ISML=1
          DO 310 J=1,NW(IO)
            POLE=OM/CS(IO,J)
            IF((DABS(WA-POLE).LT.DR).OR.(DABS(WB-POLE).LT.DR))ISML=0
 310      CONTINUE
          IF(DABS(WB-WNR(4)).LT.DR)ISML=0
        END IF
        IF(ISML.EQ.1)THEN
          DO 320 K=I+1,MW
            WVT(K)=WVT(K+1)
 320      CONTINUE
          MW=MW-1
        END IF
 300  CONTINUE

      DO 400 J=1,MW
        WVN(2*(N11+N12+J)-1)=WVT(J)
        WVN(2*(N11+N12+J)  )=(WVT(J)+WVT(J+1))/2.D0
 400  CONTINUE
      WVN(2*(N11+N12+MW)+1)=WVT(MW+1)

      M13=MW

      return
      end
C
C
      SUBROUTINE WVNM(W,MW)
      implicit real*8(a-h,o-z)
      COMMON /TMP/WVT(1001) 
C
      DO 100 I=1,MW
        WA=WVT(I)
        WB=WVT(I+1)
        IF((W.GT.WA).AND.(W.LT.WB))THEN
          DO 105 K=MW+1,I+1,-1
            WVT(K+1)=WVT(K)
 105      CONTINUE
          WVT(I+1)=W
          MW=MW+1
          RETURN
        END IF
 100  CONTINUE

      RETURN
      END
C
C
C *** THE COMPUTATION OF Q, V, AND W AND THEIR DERIVATIVES 
C     USING SIMPSON'S INTEGATION FOMULA (before k*r < 10)***
C
      SUBROUTINE QVWRZS(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)
     *          ,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,N13,N14,M13
     *       /WVNUM/WVN(2001),DWN(4),WNR(5),NWN
     *       /FCQVW/SM,RAM,RR  /IRANGE/NRNG,IRNG
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/Vp(max) *
C
      NRNG=0
      IRNG=0
C
      H=DWN(1)
      K=1
      WA=WNR(1)
      CALL INTEGS(OM,K,WA,QXA,QZA,VXA,WXA,WZA)
      DO  100 I=1,N11
        K=2*I
        WB=WVN(K)
        CALL INTEGS(OM,K,WB,QXB,QZB,VXB,WXB,WZB)
        K=K+1
        WC=WVN(K)
        CALL INTEGS(OM,K,WC,QXC,QZC,VXC,WXC,WZC)

        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

        IF(RR*WC.GT.10.D0)THEN
          NRNG=1
          IRNG=I
          IF(I.EQ.N11)THEN
            NRNG=2
            IRNG=0
          END IF
c      write(*,*)'s1:i,n11,rr*wc',i,n11,rr*wc
c      write(*,*)'s1:gqx0',gqx(0)
          RETURN
        END IF

  100 CONTINUE

c      write(*,*)'s1:gqx0',gqx(0)
C
C   * RANGE 2 : From near om/Vp(max) to near om/Vs(max) *
C
      H=DWN(2)
      WA=WNR(2)
      DO  200 I=1,N12
        K=2*(N11+I)
        WB=WVN(K)
        CALL INTEGS(OM,K,WB,QXB,QZB,VXB,WXB,WZB)
        WC=WA+H*2*I
        K=K+1
        WC=WVN(K)
        CALL INTEGS(OM,K,WC,QXC,QZC,VXC,WXC,WZC)

        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

        IF(RR*WC.GT.10.D0)THEN
          NRNG=2
          IRNG=I
          IF(I.EQ.N12)THEN
            NRNG=3
            IRNG=0
          END IF
c      write(*,*)'s2:i,n12,rr*wc',i,n12,rr*wc
c      write(*,*)'s2:gqx0',gqx(0)
          RETURN
        END IF

  200 CONTINUE

c      write(*,*)'s2:gqx0',gqx(0)
C
C   * RANGE 3 : From near om/Vs(max) to near om/Rayleigh(min) *
C
      WA=WNR(3)
      DO  300 I=1,M13
        K=2*(N11+N12+I)
        WB=WVN(K)
        CALL INTEGS(OM,K,WB,QXB,QZB,VXB,WXB,WZB)
        K=K+1
        WC=WVN(K)
        CALL INTEGS(OM,K,WC,QXC,QZC,VXC,WXC,WZC)

        H=WC-WB
        DO 310 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)
  310   CONTINUE

        IF(RR*WC.GT.10.D0)THEN
          NRNG=3
          IRNG=I
          IF(I.EQ.M13)THEN
            NRNG=4
            IRNG=0
c      write(*,*)'s3:i,M13,rr*wc',i,M13,rr*wc
c      write(*,*)'s3:gqx0',gqx(0)
          END IF
          RETURN
        END IF

  300 CONTINUE

c      write(*,*)'s3:gqx0',gqx(0)
C
C   * RANGE 4 : From near om/Rayleigh(min) to om/C(final) *
C
      H=DWN(4)
      WA=WNR(4)
      DO  400 I=1,N14
        K=2*(N11+N12+M13+I)
        WB=WVN(K)
        CALL INTEGS(OM,K,WB,QXB,QZB,VXB,WXB,WZB)
        K=K+1
        WC=WVN(K)
        CALL INTEGS(OM,K,WC,QXC,QZC,VXC,WXC,WZC)

        DO 410 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)
  410   CONTINUE

        IF(RR*WC.GT.10.D0)THEN
          NRNG=4
          IRNG=I
          IF(I.EQ.N14)THEN
            NRNG=0
            IRNG=0
c      write(*,*)'s4:i,n14,rr*wc',i,n14,rr*wc
c      write(*,*)'s4:gqx0',gqx(0)
          END IF
          RETURN
        END IF

  400 CONTINUE

c      write(*,*)'final s4:gqx0',gqx(0)
C
      RETURN
      END
C
C
C *** THE INTEGRANDS OF Q, V, AND W AND THEIR DERIVATIVES for Simpson's ***
C
      SUBROUTINE INTEGS(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),VLXZ,URXZ(2),URZZ(2),SM,RAM
     *   ,H1X,H2X,V1X,V2X,V3X,V4X,V1Z,V2Z,V3Z,V4Z,VP,VS
      DIMENSION BJ0(0:1),BJ1(0:1),BJK(0:1),BJD(0:1)
      COMMON /SOLY/SDEP,ODEP,LS,LO 
     *       /STRC/DNS(20),VP(20),VS(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
      CALL RTSH(OM,K,WN,VLX)
      CALL RTPSV(OM,K,WN,URX,URZ)
C
      IF((IST.EQ.1).AND.(IABS(LS-LO).LE.1))THEN
        CALL ASYM(LS,NL,WN,H1X,H2X,V1X,V2X,V3X,V4X,V1Z,V2Z,V3Z,V4Z)
        VLX(1)=VLX(1)-H1X
        URX(1,1)=URX(1,1)-V1X
        URX(2,1)=URX(2,1)-V2X
        URZ(1,1)=URZ(1,1)-V1Z
        URZ(2,1)=URZ(2,1)-V2Z
        VLX(2)=VLX(2)-H2X
        URX(1,2)=URX(1,2)-V3X
        URX(2,2)=URX(2,2)-V4X
        URZ(1,2)=URZ(1,2)-V3Z
        URZ(2,2)=URZ(2,2)-V4Z
      END IF
C
      VLXZ=VLX(2)/SM
      URXZ(1)=URX(1,2)/SM+WN*URX(2,1)
      URZZ(1)=URZ(1,2)/SM+WN*URZ(2,1)
      URXZ(2)=(URX(2,2)-WN*RAM*URX(1,1))/(RAM+2.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(19,1100)WN
        WRITE(20,1100)dreal(QX(0))
        WRITE(21,1100)dimag(QX(0))
        WRITE(22,1100)dreal(VX(0))
        WRITE(23,1100)dimag(VX(0))
        WRITE(24,1100)dreal(WX(0))
        WRITE(25,1100)dimag(WX(0))
        WRITE(26,1100)dreal(QZ(0))
        WRITE(27,1100)dimag(QZ(0))
        WRITE(28,1100)dreal(WZ(0))
        WRITE(29,1100)dimag(WZ(0))
        WRITE(30,1100)dreal(QX(1))
        WRITE(31,1100)dimag(QX(1))
        WRITE(32,1100)dreal(VX(1))
        WRITE(33,1100)dimag(VX(1))
        WRITE(34,1100)dreal(WX(1))
        WRITE(35,1100)dimag(WX(1))
        WRITE(36,1100)dreal(QZ(1))
        WRITE(37,1100)dimag(QZ(1))
        WRITE(38,1100)dreal(WZ(1))
        WRITE(39,1100)dimag(WZ(1))
        WRITE(40,1100)dreal(QX(2))
        WRITE(41,1100)dimag(QX(2))
        WRITE(42,1100)dreal(VX(2))
        WRITE(43,1100)dimag(VX(2))
        WRITE(44,1100)dreal(WX(2))
        WRITE(45,1100)dimag(WX(2))
        WRITE(46,1100)dreal(QZ(2))
        WRITE(47,1100)dimag(QZ(2))
        WRITE(48,1100)dreal(WZ(2))
        WRITE(49,1100)dimag(WZ(2))
      END IF
 1100 FORMAT(E15.7)
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 > 10.0) ***
C
      SUBROUTINE QVWRZF(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)
     *      ,QXS2(0:2),QZS2(0:2),VXS2(0:2),WXS2(0:2),WZS2(0:2)
     *      ,QXC2(0:2),QZC2(0:2),VXC2(0:2),WXC2(0:2),WZC2(0:2)
     *      ,SM,RAM,S1,S2
      COMMON /FCQVW/SM,RAM,RR  /IRANGE/NRNG,IRNG
     *       /INTG/VPMAX,VSMAX,RYMIN,CFINAL,N11,N12,N13,N14,M13
     *       /WVNUM/WVN(2001),DWN(4),WNR(5),NWN
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/Vp(max) *
C
      IF(NRNG.EQ.0)RETURN

      JRNG=0
      IF(NRNG.GT.1)GOTO 10
C
      H=DWN(1)
      K=2*IRNG+1
      WA=WVN(K)
      CALL INTEGF(OM,K,WA,QXSA,QZSA,VXSA,WXSA,WZSA
     *                   ,QXCA,QZCA,VXCA,WXCA,WZCA)

      DO 100 J=0,2
        QXS1(J)=0.5D0*QXSA(J)
        QZS1(J)=0.5D0*QZSA(J)
        VXS1(J)=0.5D0*VXSA(J)
        WXS1(J)=0.5D0*WXSA(J)
        WZS1(J)=0.5D0*WZSA(J)
        QXC1(J)=0.5D0*QXCA(J)
        QZC1(J)=0.5D0*QZCA(J)
        VXC1(J)=0.5D0*VXCA(J)
        WXC1(J)=0.5D0*WXCA(J)
        WZC1(J)=0.5D0*WZCA(J)
        QXS2(J)=(0.D0,0.D0)
        QZS2(J)=(0.D0,0.D0)
        VXS2(J)=(0.D0,0.D0)
        WXS2(J)=(0.D0,0.D0)
        WZS2(J)=(0.D0,0.D0)
        QXC2(J)=(0.D0,0.D0)
        QZC2(J)=(0.D0,0.D0)
        VXC2(J)=(0.D0,0.D0)
        WXC2(J)=(0.D0,0.D0)
        WZC2(J)=(0.D0,0.D0)
 100  CONTINUE

      DO  120 I=IRNG+1,N11
        K=2*I
        WB=WVN(K)
        CALL INTEGF(OM,K,WB,QXSB,QZSB,VXSB,WXSB,WZSB
     *                     ,QXCB,QZCB,VXCB,WXCB,WZCB)
        K=K+1
        WC=WVN(K)
        CALL INTEGF(OM,K,WC,QXSC,QZSC,VXSC,WXSC,WZSC
     *                     ,QXCC,QZCC,VXCC,WXCC,WZCC)

        DO 120 J=0,2
          QXS1(J)=QXS1(J)+QXSC(J)
          QZS1(J)=QZS1(J)+QZSC(J)
          VXS1(J)=VXS1(J)+VXSC(J)
          WXS1(J)=WXS1(J)+WXSC(J)
          WZS1(J)=WZS1(J)+WZSC(J)
          QXC1(J)=QXC1(J)+QXCC(J)
          QZC1(J)=QZC1(J)+QZCC(J)
          VXC1(J)=VXC1(J)+VXCC(J)
          WXC1(J)=WXC1(J)+WXCC(J)
          WZC1(J)=WZC1(J)+WZCC(J)
          QXS2(J)=QXS2(J)+QXSB(J)
          QZS2(J)=QZS2(J)+QZSB(J)
          VXS2(J)=VXS2(J)+VXSB(J)
          WXS2(J)=WXS2(J)+WXSB(J)
          WZS2(J)=WZS2(J)+WZSB(J)
          QXC2(J)=QXC2(J)+QXCB(J)
          QZC2(J)=QZC2(J)+QZCB(J)
          VXC2(J)=VXC2(J)+VXCB(J)
          WXC2(J)=WXC2(J)+WXCB(J)
          WZC2(J)=WZC2(J)+WZCB(J)
  120 CONTINUE

      DO 150 J=0,2
        QXS1(J)=QXS1(J)-QXSC(J)*0.5d0
        QZS1(J)=QZS1(J)-QZSC(J)*0.5d0
        VXS1(J)=VXS1(J)-VXSC(J)*0.5d0
        WXS1(J)=WXS1(J)-WXSC(J)*0.5d0
        WZS1(J)=WZS1(J)-WZSC(J)*0.5d0
        QXC1(J)=QXC1(J)-QXCC(J)*0.5d0
        QZC1(J)=QZC1(J)-QZCC(J)*0.5d0
        VXC1(J)=VXC1(J)-VXCC(J)*0.5d0
        WXC1(J)=WXC1(J)-WXCC(J)*0.5d0
        WZC1(J)=WZC1(J)-WZCC(J)*0.5d0
  150 CONTINUE

      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)

      SA=1.D0/DTAN(RR*WA)
      CA=DTAN(RR*WA)
      SC=1.D0/DTAN(RR*WNR(2))
      CC=DTAN(RR*WNR(2))
      DO 170 J=0,2
        S1=A1*(QXSA(J)*SA-QXSC(J)*SC)+B1*QXS1(J)+C1*QXS2(J)
        S2=A1*(QXCC(J)*CC-QXCA(J)*CA)+B1*QXC1(J)+C1*QXC2(J)
        GQX(J)=GQX(J)+H*(S1+S2)
        S1=A1*(QZSA(J)*SA-QZSC(J)*SC)+B1*QZS1(J)+C1*QZS2(J)
        S2=A1*(QZCC(J)*CC-QZCA(J)*CA)+B1*QZC1(J)+C1*QZC2(J)
        GQZ(J)=GQZ(J)+H*(S1+S2)
        S1=A1*(VXSA(J)*SA-VXSC(J)*SC)+B1*VXS1(J)+C1*VXS2(J)
        S2=A1*(VXCC(J)*CC-VXCA(J)*CA)+B1*VXC1(J)+C1*VXC2(J)
        GVX(J)=GVX(J)+H*(S1+S2)
        S1=A1*(WXSA(J)*SA-WXSC(J)*SC)+B1*WXS1(J)+C1*WXS2(J)
        S2=A1*(WXCC(J)*CC-WXCA(J)*CA)+B1*WXC1(J)+C1*WXC2(J)
        GWX(J)=GWX(J)+H*(S1+S2)
        S1=A1*(WZSA(J)*SA-WZSC(J)*SC)+B1*WZS1(J)+C1*WZS2(J)
        S2=A1*(WZCC(J)*CC-WZCA(J)*CA)+B1*WZC1(J)+C1*WZC2(J)
        GWZ(J)=GWZ(J)+H*(S1+S2)
        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)
  170 CONTINUE

c      write(*,*)'f1:gqx0',gqx(0)
C
C   * RANGE 2 : From near om/Vp(max)to near om/Vs(max) *
C
      JRNG=1
      IRNG=0
 10   IF(NRNG.GT.2)GOTO 20
C
      H=DWN(2)
      K=2*(N11+IRNG)+1
      WA=WVN(K)
      IF(JRNG.EQ.1)GOTO 15

      CALL INTEGF(OM,K,WA,QXSA,QZSA,VXSA,WXSA,WZSA
     *                   ,QXCA,QZCA,VXCA,WXCA,WZCA)

 15   DO 200 J=0,2
        QXS1(J)=0.5D0*QXSA(J)
        QZS1(J)=0.5D0*QZSA(J)
        VXS1(J)=0.5D0*VXSA(J)
        WXS1(J)=0.5D0*WXSA(J)
        WZS1(J)=0.5D0*WZSA(J)
        QXC1(J)=0.5D0*QXCA(J)
        QZC1(J)=0.5D0*QZCA(J)
        VXC1(J)=0.5D0*VXCA(J)
        WXC1(J)=0.5D0*WXCA(J)
        WZC1(J)=0.5D0*WZCA(J)
        QXS2(J)=(0.D0,0.D0)
        QZS2(J)=(0.D0,0.D0)
        VXS2(J)=(0.D0,0.D0)
        WXS2(J)=(0.D0,0.D0)
        WZS2(J)=(0.D0,0.D0)
        QXC2(J)=(0.D0,0.D0)
        QZC2(J)=(0.D0,0.D0)
        VXC2(J)=(0.D0,0.D0)
        WXC2(J)=(0.D0,0.D0)
        WZC2(J)=(0.D0,0.D0)
 200  CONTINUE

      DO  220 I=IRNG+1,N12
        K=2*(N11+I)
        WB=WVN(K)
        CALL INTEGF(OM,K,WB,QXSB,QZSB,VXSB,WXSB,WZSB
     *                     ,QXCB,QZCB,VXCB,WXCB,WZCB)
        K=K+1
        WC=WVN(K)
        CALL INTEGF(OM,K,WC,QXSC,QZSC,VXSC,WXSC,WZSC
     *                     ,QXCC,QZCC,VXCC,WXCC,WZCC)

        DO 220 J=0,2
          QXS1(J)=QXS1(J)+QXSC(J)
          QZS1(J)=QZS1(J)+QZSC(J)
          VXS1(J)=VXS1(J)+VXSC(J)
          WXS1(J)=WXS1(J)+WXSC(J)
          WZS1(J)=WZS1(J)+WZSC(J)
          QXC1(J)=QXC1(J)+QXCC(J)
          QZC1(J)=QZC1(J)+QZCC(J)
          VXC1(J)=VXC1(J)+VXCC(J)
          WXC1(J)=WXC1(J)+WXCC(J)
          WZC1(J)=WZC1(J)+WZCC(J)
          QXS2(J)=QXS2(J)+QXSB(J)
          QZS2(J)=QZS2(J)+QZSB(J)
          VXS2(J)=VXS2(J)+VXSB(J)
          WXS2(J)=WXS2(J)+WXSB(J)
          WZS2(J)=WZS2(J)+WZSB(J)
          QXC2(J)=QXC2(J)+QXCB(J)
          QZC2(J)=QZC2(J)+QZCB(J)
          VXC2(J)=VXC2(J)+VXCB(J)
          WXC2(J)=WXC2(J)+WXCB(J)
          WZC2(J)=WZC2(J)+WZCB(J)
  220 CONTINUE

      DO 250 J=0,2
        QXS1(J)=QXS1(J)-QXSC(J)*0.5d0
        QZS1(J)=QZS1(J)-QZSC(J)*0.5d0
        VXS1(J)=VXS1(J)-VXSC(J)*0.5d0
        WXS1(J)=WXS1(J)-WXSC(J)*0.5d0
        WZS1(J)=WZS1(J)-WZSC(J)*0.5d0
        QXC1(J)=QXC1(J)-QXCC(J)*0.5d0
        QZC1(J)=QZC1(J)-QZCC(J)*0.5d0
        VXC1(J)=VXC1(J)-VXCC(J)*0.5d0
        WXC1(J)=WXC1(J)-WXCC(J)*0.5d0
        WZC1(J)=WZC1(J)-WZCC(J)*0.5d0
  250 CONTINUE

      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)

      SA=1.D0/DTAN(RR*WA)
      CA=DTAN(RR*WA)
      SC=1.D0/DTAN(RR*WNR(3))
      CC=DTAN(RR*WNR(3))
      DO 270 J=0,2
        S1=A1*(QXSA(J)*SA-QXSC(J)*SC)+B1*QXS1(J)+C1*QXS2(J)
        S2=A1*(QXCC(J)*CC-QXCA(J)*CA)+B1*QXC1(J)+C1*QXC2(J)
        GQX(J)=GQX(J)+H*(S1+S2)
        S1=A1*(QZSA(J)*SA-QZSC(J)*SC)+B1*QZS1(J)+C1*QZS2(J)
        S2=A1*(QZCC(J)*CC-QZCA(J)*CA)+B1*QZC1(J)+C1*QZC2(J)
        GQZ(J)=GQZ(J)+H*(S1+S2)
        S1=A1*(VXSA(J)*SA-VXSC(J)*SC)+B1*VXS1(J)+C1*VXS2(J)
        S2=A1*(VXCC(J)*CC-VXCA(J)*CA)+B1*VXC1(J)+C1*VXC2(J)
        GVX(J)=GVX(J)+H*(S1+S2)
        S1=A1*(WXSA(J)*SA-WXSC(J)*SC)+B1*WXS1(J)+C1*WXS2(J)
        S2=A1*(WXCC(J)*CC-WXCA(J)*CA)+B1*WXC1(J)+C1*WXC2(J)
        GWX(J)=GWX(J)+H*(S1+S2)
        S1=A1*(WZSA(J)*SA-WZSC(J)*SC)+B1*WZS1(J)+C1*WZS2(J)
        S2=A1*(WZCC(J)*CC-WZCA(J)*CA)+B1*WZC1(J)+C1*WZC2(J)
        GWZ(J)=GWZ(J)+H*(S1+S2)
        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)
  270 CONTINUE

c      write(*,*)'f2:gqx0',gqx(0)
C
C   * RANGE 3 : From near om/Vs(max)to near om/Rayleigh(min) *
C
      JRNG=2
      IRNG=0
C
 20   IF(NRNG.GT.3)GOTO 30
C
      K=2*(N11+N12+IRNG)+1
      WA=WVN(K)
      IF(JRNG.EQ.2)GOTO 25

      CALL INTEGF(OM,K,WA,QXSA,QZSA,VXSA,WXSA,WZSA
     *                   ,QXCA,QZCA,VXCA,WXCA,WZCA)

 25   DO  300 I=IRNG+1,M13

        K=2*(N11+N12+I)
        WB=WVN(K)
        CALL INTEGF(OM,K,WB,QXSB,QZSB,VXSB,WXSB,WZSB
     *                     ,QXCB,QZCB,VXCB,WXCB,WZCB)
        K=K+1
        WC=WVN(K)
        CALL INTEGF(OM,K,WC,QXSC,QZSC,VXSC,WXSC,WZSC
     *                     ,QXCC,QZCC,VXCC,WXCC,WZCC)

        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)

        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)

        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

          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)

          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)

          WA=WC
  350   CONTINUE

  300 CONTINUE

c      write(*,*)'f3:gqx0',gqx(0)
C
C   * RANGE 4 : From om/Rayleigh(min) to om/C(final) *
C
      JRNG=3
      IRNG=0
C
 30   H=DWN(4)
      K=2*(N11+N12+M13+IRNG)+1
      WA=WVN(K)
      IF(JRNG.EQ.3)GOTO 35

      CALL INTEGF(OM,K,WA,QXSA,QZSA,VXSA,WXSA,WZSA
     *                   ,QXCA,QZCA,VXCA,WXCA,WZCA)

 35   DO 400 J=0,2
        QXS1(J)=0.5D0*QXSA(J)
        QZS1(J)=0.5D0*QZSA(J)
        VXS1(J)=0.5D0*VXSA(J)
        WXS1(J)=0.5D0*WXSA(J)
        WZS1(J)=0.5D0*WZSA(J)
        QXC1(J)=0.5D0*QXCA(J)
        QZC1(J)=0.5D0*QZCA(J)
        VXC1(J)=0.5D0*VXCA(J)
        WXC1(J)=0.5D0*WXCA(J)
        WZC1(J)=0.5D0*WZCA(J)
        QXS2(J)=(0.D0,0.D0)
        QZS2(J)=(0.D0,0.D0)
        VXS2(J)=(0.D0,0.D0)
        WXS2(J)=(0.D0,0.D0)
        WZS2(J)=(0.D0,0.D0)
        QXC2(J)=(0.D0,0.D0)
        QZC2(J)=(0.D0,0.D0)
        VXC2(J)=(0.D0,0.D0)
        WXC2(J)=(0.D0,0.D0)
        WZC2(J)=(0.D0,0.D0)
 400  CONTINUE

      DO  420 I=IRNG+1,N14
        K=2*(N11+N12+M13+I)
        WB=WVN(K)
        CALL INTEGF(OM,K,WB,QXSB,QZSB,VXSB,WXSB,WZSB
     *                     ,QXCB,QZCB,VXCB,WXCB,WZCB)
        K=K+1
        WC=WVN(K)
        CALL INTEGF(OM,K,WC,QXSC,QZSC,VXSC,WXSC,WZSC
     *                     ,QXCC,QZCC,VXCC,WXCC,WZCC)

        DO 420 J=0,2
          QXS1(J)=QXS1(J)+QXSC(J)
          QZS1(J)=QZS1(J)+QZSC(J)
          VXS1(J)=VXS1(J)+VXSC(J)
          WXS1(J)=WXS1(J)+WXSC(J)
          WZS1(J)=WZS1(J)+WZSC(J)
          QXC1(J)=QXC1(J)+QXCC(J)
          QZC1(J)=QZC1(J)+QZCC(J)
          VXC1(J)=VXC1(J)+VXCC(J)
          WXC1(J)=WXC1(J)+WXCC(J)
          WZC1(J)=WZC1(J)+WZCC(J)
          QXS2(J)=QXS2(J)+QXSB(J)
          QZS2(J)=QZS2(J)+QZSB(J)
          VXS2(J)=VXS2(J)+VXSB(J)
          WXS2(J)=WXS2(J)+WXSB(J)
          WZS2(J)=WZS2(J)+WZSB(J)
          QXC2(J)=QXC2(J)+QXCB(J)
          QZC2(J)=QZC2(J)+QZCB(J)
          VXC2(J)=VXC2(J)+VXCB(J)
          WXC2(J)=WXC2(J)+WXCB(J)
          WZC2(J)=WZC2(J)+WZCB(J)
  420 CONTINUE

      DO 450 J=0,2
        QXS1(J)=QXS1(J)-QXSC(J)*0.5d0
        QZS1(J)=QZS1(J)-QZSC(J)*0.5d0
        VXS1(J)=VXS1(J)-VXSC(J)*0.5d0
        WXS1(J)=WXS1(J)-WXSC(J)*0.5d0
        WZS1(J)=WZS1(J)-WZSC(J)*0.5d0
        QXC1(J)=QXC1(J)-QXCC(J)*0.5d0
        QZC1(J)=QZC1(J)-QZCC(J)*0.5d0
        VXC1(J)=VXC1(J)-VXCC(J)*0.5d0
        WXC1(J)=WXC1(J)-WXCC(J)*0.5d0
        WZC1(J)=WZC1(J)-WZCC(J)*0.5d0
  450 CONTINUE

      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)

      SA=1.D0/DTAN(RR*WA)
      CA=DTAN(RR*WA)
      SC=1.D0/DTAN(RR*WNR(5))
      CC=DTAN(RR*WNR(5))
      DO 470 J=0,2
        S1=A1*(QXSA(J)*SA-QXSC(J)*SC)+B1*QXS1(J)+C1*QXS2(J)
        S2=A1*(QXCC(J)*CC-QXCA(J)*CA)+B1*QXC1(J)+C1*QXC2(J)
        GQX(J)=GQX(J)+H*(S1+S2)
        S1=A1*(QZSA(J)*SA-QZSC(J)*SC)+B1*QZS1(J)+C1*QZS2(J)
        S2=A1*(QZCC(J)*CC-QZCA(J)*CA)+B1*QZC1(J)+C1*QZC2(J)
        GQZ(J)=GQZ(J)+H*(S1+S2)
        S1=A1*(VXSA(J)*SA-VXSC(J)*SC)+B1*VXS1(J)+C1*VXS2(J)
        S2=A1*(VXCC(J)*CC-VXCA(J)*CA)+B1*VXC1(J)+C1*VXC2(J)
        GVX(J)=GVX(J)+H*(S1+S2)
        S1=A1*(WXSA(J)*SA-WXSC(J)*SC)+B1*WXS1(J)+C1*WXS2(J)
        S2=A1*(WXCC(J)*CC-WXCA(J)*CA)+B1*WXC1(J)+C1*WXC2(J)
        GWX(J)=GWX(J)+H*(S1+S2)
        S1=A1*(WZSA(J)*SA-WZSC(J)*SC)+B1*WZS1(J)+C1*WZS2(J)
        S2=A1*(WZCC(J)*CC-WZCA(J)*CA)+B1*WZC1(J)+C1*WZC2(J)
        GWZ(J)=GWZ(J)+H*(S1+S2)
  470 CONTINUE

c      write(*,*)'f4:gqx0',gqx(0)
c
      RETURN
      END
C
C
C *** THE INTEGRANDS OF Q, V, AND W AND THEIR DERIVATIVES for Filon's ***
C
      SUBROUTINE INTEGF(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),VLXZ,URXZ(2),URZZ(2),SM,RAM
     *   ,H1X,H2X,V1X,V2X,V3X,V4X,V1Z,V2Z,V3Z,V4Z,VP,VS
      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 /SOLY/SDEP,ODEP,LS,LO 
     *       /STRC/DNS(20),VP(20),VS(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
      CALL RTSH(OM,K,WN,VLX)
      CALL RTPSV(OM,K,WN,URX,URZ)
C
      IF((IST.EQ.1).AND.(IABS(LS-LO).LE.1))THEN
        CALL ASYM(LS,NL,WN,H1X,H2X,V1X,V2X,V3X,V4X,V1Z,V2Z,V3Z,V4Z)
        VLX(1)=VLX(1)-H1X
        URX(1,1)=URX(1,1)-V1X
        URX(2,1)=URX(2,1)-V2X
        URZ(1,1)=URZ(1,1)-V1Z
        URZ(2,1)=URZ(2,1)-V2Z
        VLX(2)=VLX(2)-H2X
        URX(1,2)=URX(1,2)-V3X
        URX(2,2)=URX(2,2)-V4X
        URZ(1,2)=URZ(1,2)-V3Z
        URZ(2,2)=URZ(2,2)-V4Z
      END IF
C
      VLXZ=VLX(2)/SM
      URXZ(1)=URX(1,2)/SM+WN*URX(2,1)
      URZZ(1)=URZ(1,2)/SM+WN*URZ(2,1)
      URXZ(2)=(URX(2,2)-WN*RAM*URX(1,1))/(RAM+2.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)
      A=1.D0/(8.D0*RK)
C
      BJ0S(0)= SRK*(1.D0+A)
      BJ0C(0)= CRK*(1.D0-A)
      BJ1S(0)= SRK*(1.D0+3.D0*A)
      BJ1C(0)=-CRK*(1.D0-3.D0*A)
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(19,1100)WN
        WRITE(20,1100)dreal(QXS(0)+QXC(0))
        WRITE(21,1100)dimag(QXS(0)+QXC(0))
        WRITE(22,1100)dreal(VXS(0)+VXC(0))
        WRITE(23,1100)dimag(VXS(0)+VXC(0))
        WRITE(24,1100)dreal(WXS(0)+WXC(0))
        WRITE(25,1100)dimag(WXS(0)+WXC(0))
        WRITE(26,1100)dreal(QZS(0)+QZC(0))
        WRITE(27,1100)dimag(QZS(0)+QZC(0))
        WRITE(28,1100)dreal(WZS(0)+WZC(0))
        WRITE(29,1100)dimag(WZS(0)+WZC(0))
        WRITE(30,1100)dreal(QXS(1)+QXC(1))
        WRITE(31,1100)dimag(QXS(1)+QXC(1))
        WRITE(32,1100)dreal(VXS(1)+VXC(1))
        WRITE(33,1100)dimag(VXS(1)+VXC(1))
        WRITE(34,1100)dreal(WXS(1)+WXC(1))
        WRITE(35,1100)dimag(WXS(1)+WXC(1))
        WRITE(36,1100)dreal(QZS(1)+QZC(1))
        WRITE(37,1100)dimag(QZS(1)+QZC(1))
        WRITE(38,1100)dreal(WZS(1)+WZC(1))
        WRITE(39,1100)dimag(WZS(1)+WZC(1))
        WRITE(40,1100)dreal(QXS(2)+QXC(2))
        WRITE(41,1100)dimag(QXS(2)+QXC(2))
        WRITE(42,1100)dreal(VXS(2)+VXC(2))
        WRITE(43,1100)dimag(VXS(2)+VXC(2))
        WRITE(44,1100)dreal(WXS(2)+WXC(2))
        WRITE(45,1100)dimag(WXS(2)+WXC(2))
        WRITE(46,1100)dreal(QZS(2)+QZC(2))
        WRITE(47,1100)dimag(QZS(2)+QZC(2))
        WRITE(48,1100)dreal(WZS(2)+WZC(2))
        WRITE(49,1100)dimag(WZS(2)+WZC(2))
      END IF
 1100 FORMAT(E15.7)
C
      RETURN
      END
C
C
C *** Library for SH UPGOING/DOWNGOING COEFFICIENTS ***
C
      SUBROUTINE CUDSH(OM)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VP,VS,CD,CU,VS2,WS2,CMS,SM,EE(2,2),DD,DU
     *      ,CUSU,CDSU,CDSD,CUSD,GTD,GTU,GRD,GRU,AU,BU,AD,BD,CS
      COMMON /SOLY/SDEP,ODEP,LS,LO 
     *       /STRC/DNS(20),VP(20),VS(20),THK(20),NL
     *       /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(4),WNR(5),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
      VS2=VS(LS)*VS(LS)
      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(REAL(CMS).LT.0.0)CMS=-CMS
C
        EE(2,1)=-SM*CMS
        EE(2,2)= SM*CMS
C
        PI=4.D0*DATAN(1.D0)
        DQ=-WN/(2.0*PI)
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 CUDPSV(OM)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VP,VS,CMP,CMS,EE(4,4),DD(2,2),DU(2,2)
     *    ,CDX,CUX,CUSUX,CDSUX,CDSDX,CUSDX,GTD,GTU,GRD,GRU,DT
     *    ,CDZ,CUZ,CUSUZ,CDSUZ,CDSDZ,CUSDZ
     *    ,AU(2,2),BU(2,2),AD(2,2),BD(2,2),CS(2,2),E1(2,2),E2(2,2)
      COMMON /SOLY/SDEP,ODEP,LS,LO 
     *  /STRC/DNS(20),VP(20),VS(20),THK(20),NL
     *  /DUPSVX/CDX(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(4),WNR(5),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
        PI=4.D0*DATAN(1.D0)
        DQX=-WN/(2.0*PI)
        DQZ= WN/(2.0*PI)
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 RTSH(OM,K,WN,VLX)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VP,VS,CD,CU,VS2,WS2,CMS,SM,EE(2,2)
     *          ,CUSU,CDSU,CDSD,CUSD,DD,DU,VLX(2)
      COMMON /SOLY/SDEP,ODEP,LS,LO
     *       /STRC/DNS(20),VP(20),VS(20),THK(20),NL
     *       /DUSH/CD(2001,20),CU(2001,20)
     *       ,CUSU(2001),CDSU(2001),CDSD(2001),CUSD(2001)
     *       /WVNUM/WVN(2001),DWN(4),WNR(5),NWN
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
      NC=2001
      IF(NWN.GT.NC)THEN
         WRITE(10,*)' FROM RTSH: NWN.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C *** COMPUTATION OF VLX ***
C
      OM2=OM*OM
      WN2=WN*WN
C
      VS2=VS(LO)*VS(LO)
      WS2=OM2/VS2
      CMS=CDSQRT(WN2-WS2)
      IF(REAL(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 RTPSV(OM,K,WN,URX,URZ)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VP,VS,CMP,CMS,EE(4,4),DD(2,2),DU(2,2)
     *    ,CDX,CUX,CUSUX,CDSUX,CDSDX,CUSDX,URX(2,2),E1(2,2)
     *    ,CDZ,CUZ,CUSUZ,CDSUZ,CDSDZ,CUSDZ,URZ(2,2),E2(2,2)
      COMMON /SOLY/SDEP,ODEP,LS,LO
     *       /STRC/DNS(20),VP(20),VS(20),THK(20),NL
     * /DUPSVX/CDX(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) 
     * /WVNUM/WVN(2001),DWN(4),WNR(5),NWN
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
      NC=2001
      IF(NWN.GT.NC)THEN
         WRITE(10,*)' FROM RTPSV: NWN.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 *** COEFFICIENTS FOR ASYMPTOTIC SOLUTIONS ***
C
      SUBROUTINE CASYM(LS,NL)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 SM1,X1,SM2,X2,SM3,X3,BB(4,4),BB1(4,4),BB2(4,4)
     *   ,H1T,H2T,H1D,H2D,RH,H1R1,H2R1,H1R2,H2R2       
     *   ,CUPX,CUSX,CUP1X,CUS1X,CUP2X,CUS2X
     *   ,CDPX,CDSX,CDP1X,CDS1X,CDP2X,CDS2X
     *   ,CUPZ,CUSZ,CUP1Z,CUS1Z,CUP2Z,CUS2Z
     *   ,CDPZ,CDSZ,CDP1Z,CDS1Z,CDP2Z,CDS2Z
      COMMON/ASYMP/SM1,X1,SM2,X2,SM3,X3
     *      /DPSO /DEP1,DEP2,SDP,ODP
     *      /HX  /H1T,H2T,H1D,H2D,H1R1,H2R1,H1R2,H2R2       
     *      /VUX /CUPX,CUSX,CUP1X,CUS1X,CUP2X,CUS2X
     *      /VDX /CDPX,CDSX,CDP1X,CDS1X,CDP2X,CDS2X
     *      /VUZ /CUPZ,CUSZ,CUP1Z,CUS1Z,CUP2Z,CUS2Z
     *      /VDZ /CDPZ,CDSZ,CDP1Z,CDS1Z,CDP2Z,CDS2Z
C
C *** ODP < DEP1 ***
C
      IF(ODP.GE.DEP1) GOTO 10
      IF(LS.EQ.1) GOTO 10
C
C  ** QX **
C
C   * SH (TRANSMITTED) *
C
      H1T=0.15915494309D0/(SM1+SM2)
      H2T=SM1*H1T
C
C   * P-SV (TRANSMITTED) *
C
      CUPX=0.25d0/(3.141592654d0*SM2*(X2+1.0))
      CUSX=CUPX*(X2-1.0)
C
      CALL EINVE(SM1,SM2,X1,X2,BB)
C
      CUP1X=BB(3,3)*CUPX+BB(3,4)*CUSX
      CUS1X=BB(4,3)*CUPX+BB(4,4)*CUSX
      CUP2X=BB(3,4)*CUPX
      CUS2X=BB(4,4)*CUPX
C
C  ** QZ **
C
C   * P-SV (TRANSMITTED) *
C
      CUSZ=0.25d0/(3.141592654d0*SM2*(X2+1.0))
      CUPZ=-CUSZ
C
      CUP1Z=BB(3,3)*CUPZ+BB(3,4)*CUSZ
      CUS1Z=BB(4,3)*CUPZ+BB(4,4)*CUSZ
      CUP2Z=BB(3,4)*CUPZ
      CUS2Z=BB(4,4)*CUPZ
C
      GOTO 50
C
C *** DEP1 < (OR =) ODP < SDP ***
C
   10 IF(ODP.GE.SDP) GOTO 20
C
C  ** QX **
C
C   * SH *
C
      H2D=0.07957747155D0
      H1D=H2D/SM2
      RH=(SM2-SM1)/(SM2+SM1)
      IF(LS.EQ.1)RH=(1.D0,0.D0)
      H1R1= RH*H1D
      H2R1=-RH*H2D
      RH=(SM3-SM2)/(SM3+SM2)
      IF(LS.EQ.NL)RH=(0.D0,0.D0)
      H1R2=-RH*H1D
      H2R2=-RH*H2D
C
C   * P-SV (UP-GOING DIRECT) *
C
      CUPX=0.25d0/(3.141592654d0*SM2*(X2+1.0))
      CUSX=CUPX*(X2-1.0)
C
C   * P-SV (DOWN-GOING REFLECTED FROM UPPER BOUNDARY) *
C
      IF(LS.EQ.1)THEN
        BB1(1,3)=-(X2-2.0)
        BB1(1,4)=2.D0
        BB1(2,3)=-0.5D0*(X2-1.D0)*(X2-3.D0)
        BB1(2,4)=X2-2.D0
      ELSE 
        CALL EINVE(SM1,SM2,X1,X2,BB1)
      END IF
C
      CDP1X=BB1(1,3)*CUPX+BB1(1,4)*CUSX
      CDS1X=BB1(2,3)*CUPX+BB1(2,4)*CUSX
      CDP2X=BB1(1,4)*CUPX
      CDS2X=BB1(2,4)*CUPX
C
C   * P-SV (UP-GOING REFLECTED FROM LOWER BOUNDARY) *
C
      IF(LS.EQ.NL)GOTO 15
C
      CDPX=CUPX
      CDSX=CUSX
C
      CALL EINVE(SM2,SM3,X2,X3,BB2)
C
      CUP1X=BB2(3,1)*CDPX+BB2(3,2)*CDSX
      CUS1X=BB2(4,1)*CDPX+BB2(4,2)*CDSX
      CUP2X=BB2(3,2)*CDPX
      CUS2X=BB2(4,2)*CDPX
C
C  ** QZ **
C
   15 CUSZ=0.25d0/(3.141592654d0*SM2*(X2+1.0))
      CUPZ=-CUSZ
C
      CDP1Z=BB1(1,3)*CUPZ+BB1(1,4)*CUSZ
      CDS1Z=BB1(2,3)*CUPZ+BB1(2,4)*CUSZ
      CDP2Z=BB1(1,4)*CUPZ
      CDS2Z=BB1(2,4)*CUPZ
C
      IF(LS.EQ.NL)GOTO 50
C
      CDPZ=CUSZ
      CDSZ=-CDPZ
C
      CUP1Z=BB2(3,1)*CDPZ+BB2(3,2)*CDSZ
      CUS1Z=BB2(4,1)*CDPZ+BB2(4,2)*CDSZ
      CUP2Z=BB2(3,2)*CDPZ
      CUS2Z=BB2(4,2)*CDPZ
C
      GOTO 50
C
C *** SDP < (OR =) ODP < DEP2 ***
C
   20 IF((LS.NE.NL).AND.(ODP.GE.DEP2)) GOTO 30
C
C  ** QX **
C
C   * SH *
C
      H2D=0.07957747155D0
      H1D=H2D/SM2
      RH=(SM2-SM1)/(SM2+SM1)
      IF(LS.EQ.1)RH=(1.D0,0.D0)
      H1R1= RH*H1D
      H2R1=-RH*H2D
      RH=(SM3-SM2)/(SM3+SM2)
      IF(LS.EQ.NL)RH=(0.D0,0.D0)
      H1R2=-RH*H1D
      H2R2=-RH*H2D
C
C   * P-SV (DOWN-GOING DIRECT) *
C
      CDPX=0.25d0/(3.141592654d0*SM2*(X2+1.0))
      CDSX=CDPX*(X2-1.0)
C
C   * P-SV (UP-GOING REFLECTED FROM LOWER BOUNDARY) *
C
      IF(LS.EQ.NL)GOTO 25
C
      CALL EINVE(SM2,SM3,X2,X3,BB2)
C
      CUP1X=BB2(3,1)*CDPX+BB2(3,2)*CDSX
      CUS1X=BB2(4,1)*CDPX+BB2(4,2)*CDSX
      CUP2X=BB2(3,2)*CDPX
      CUS2X=BB2(4,2)*CDPX
C
C   * P-SV (DOWN-GOING REFLECTED FROM UPPER BOUNDARY) *
C
   25 CUPX=CDPX
      CUSX=CDSX
C
      IF(LS.EQ.1)THEN
        BB1(1,3)=-(X2-2.0)
        BB1(1,4)=2.D0
        BB1(2,3)=-0.5D0*(X2-1.D0)*(X2-3.D0)
        BB1(2,4)=X2-2.D0
      ELSE 
        CALL EINVE(SM1,SM2,X1,X2,BB1)
      END IF
C
      CDP1X=BB1(1,3)*CUPX+BB1(1,4)*CUSX
      CDS1X=BB1(2,3)*CUPX+BB1(2,4)*CUSX
      CDP2X=BB1(1,4)*CUPX
      CDS2X=BB1(2,4)*CUPX
C
C  ** QZ **
C
      CDPZ=0.25d0/(3.141592654d0*SM2*(X2+1.0))
      CDSZ=-CDPZ
C
      IF(LS.EQ.NL)GOTO 27
C
      CUP1Z=BB2(3,1)*CDPZ+BB2(3,2)*CDSZ
      CUS1Z=BB2(4,1)*CDPZ+BB2(4,2)*CDSZ
      CUP2Z=BB2(3,2)*CDPZ
      CUS2Z=BB2(4,2)*CDPZ
C
   27 CUPZ=-CDPZ
      CUSZ=-CDSZ
C
      CDP1Z=BB1(1,3)*CUPZ+BB1(1,4)*CUSZ
      CDS1Z=BB1(2,3)*CUPZ+BB1(2,4)*CUSZ
      CDP2Z=BB1(1,4)*CUPZ
      CDS2Z=BB1(2,4)*CUPZ
C
      GOTO 50
C
C *** SDP < DEP2 < (OR =) ODP ***
C
C  ** QX **
C
C   * SH (TRANSMITTED) *
C
   30 H1T=0.15915494309D0/(SM3+SM2)
      H2T=-SM3*H1T
C
C   * P-SV (TRANSMITTED) *
C
      CDPX=0.25d0/(3.141592654d0*SM2*(X2+1.0))
      CDSX=CDPX*(X2-1.0)
C
      CALL EINVE(SM2,SM3,X2,X3,BB)
C
      CDP1X=BB(1,1)*CDPX+BB(1,2)*CDSX
      CDS1X=BB(2,1)*CDPX+BB(2,2)*CDSX
      CDP2X=BB(1,2)*CDPX
      CDS2X=BB(2,2)*CDPX
C
C  ** QZ **
C
      CDPZ=0.25d0/(3.141592654d0*SM2*(X2+1.0))
      CDSZ=-CDPZ
C
      CALL EINVE(SM2,SM3,X2,X3,BB)
C
      CDP1Z=BB(1,1)*CDPZ+BB(1,2)*CDSZ
      CDS1Z=BB(2,1)*CDPZ+BB(2,2)*CDSZ
      CDP2Z=BB(1,2)*CDPZ
      CDS2Z=BB(2,2)*CDPZ
C
   50 RETURN
C
      END
C
C
C *** PRODUCT OF INV E(L+1) * E(L) (SEE APPENDIX B IN HISADA, 1995)***
C
      SUBROUTINE EINVE(SM1,SM2,X1,X2,BB)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 BB(4,4),SM1,SM2,X1,X2,A,A2,B
C
      A=SM1/SM2
      A2=A*A
      B=0.5/(X1+A2*X2+A*(X1*X2+1.0))
C
      BB(1,1)=2.0*(X1+1.0)*A*(A+X1)
      BB(2,1)=(X1+1.0)*A*(X1*X2-3.0*X1+2.0-A*(X1*X2-3.0*X2+2.0))
      BB(3,1)=2.0*(X1-2.0)*(1.0-A)*(1.0+A*X2)
      BB(4,1)=A*(X2-1.0)*(X1*X1-3.0*X1+4.0)
     *         -(X1-1.0)*(4.0+A2*(X1-3.0)*X2)
C
      BB(1,2)=(0.0,0.0)
      BB(2,2)=2.0*(X1+1.0)*A*(1.0+A*X2)
      BB(3,2)=4.0*(A-1.0)*(1.0+A*X2)
      BB(4,2)=-BB(3,1)
C
      BB(1,3)=2.0*(X2-2.0)*(A-1.0)*(A+X1)
      BB(2,3)=A*(X1-1.0)*(X2*X2-3.0*X2+4.0)
     *         -(X2-1.0)*(4.0*A2+(X2-3.0)*X1)
      BB(3,3)=2.0*(X2+1.0)*(1.0+A*X2)
      BB(4,3)=(X2+1.0)*(A*(X1*X2-3.0*X2+2.0)-X1*X2+3.0*X1-2.0)
C
      BB(1,4)=4.0*(1.0-A)*(A+X1)
      BB(2,4)=-BB(1,3)
      BB(3,4)=(0.0,0.0)
      BB(4,4)=2.0*(X2+1.0)*(A+X1)
C
      DO 100 J=1,4
      DO 100 I=1,4
      BB(I,J)=BB(I,J)*B
  100 CONTINUE
C
      RETURN
      END
C
C
C *** ASTMPTOTIC SOLUTION ***
C
      SUBROUTINE ASYM(LS,NL,WN,H1X,H2X
     *               ,V1X,V2X,V3X,V4X,V1Z,V2Z,V3Z,V4Z)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 SM1,SM2,SM3,X1,X2,X3,A,B
     *   ,H1T,H2T,H1D,H2D,H1R1,H2R1,H1R2,H2R2       
     *   ,CUPX,CUSX,CUP1X,CUS1X,CUP2X,CUS2X
     *   ,CDPX,CDSX,CDP1X,CDS1X,CDP2X,CDS2X
     *   ,CUPZ,CUSZ,CUP1Z,CUS1Z,CUP2Z,CUS2Z
     *   ,CDPZ,CDSZ,CDP1Z,CDS1Z,CDP2Z,CDS2Z
     *   ,CSD,CPD,CST,CPT,CSR1,CPR1,CSR2,CPR2
     *   ,V1X,V2X,V3X,V4X,V1Z,V2Z,V3Z,V4Z,H1X,H2X
      COMMON /ASYMP/SM1,X1,SM2,X2,SM3,X3
     *       /DPSO /DEP1,DEP2,SDP,ODP
     *      /HX  /H1T,H2T,H1D,H2D,H1R1,H2R1,H1R2,H2R2       
     *      /VUX /CUPX,CUSX,CUP1X,CUS1X,CUP2X,CUS2X
     *      /VDX /CDPX,CDSX,CDP1X,CDS1X,CDP2X,CDS2X
     *      /VUZ /CUPZ,CUSZ,CUP1Z,CUS1Z,CUP2Z,CUS2Z
     *      /VDZ /CDPZ,CDSZ,CDP1Z,CDS1Z,CDP2Z,CDS2Z
C
C *** DISP-STRESS VECTOR DUE TO Qx FOR ODP < DEP1 ***
C
      IF(ODP.GE.DEP1)GOTO 10
      IF(LS.EQ.1)GOTO 10
C
C   * UP-GOING WAVE DUE TO TRANSMITTED WAVE *
C
      WHZ=WN*(SDP-ODP)
      WHZ1=WN*(SDP-DEP1)
      WZ1Z=WN*(DEP1-ODP)
      EX=DEXP(-WHZ)
C
      CPT=EX*(CUP1X-WHZ1*CUP2X)
      CST=EX*(CUS1X-WZ1Z*CUP1X-WHZ1*CUS2X+WZ1Z*WHZ1*CUP2X)
C
      H1X=EX*H1T
      H2X=WN*EX*H2T
C
      A=X1-1.0
      B=X1-3.0
      V1X=CPT+CST
      V2X=A*CPT-CST
      V3X=WN*SM1*(-B*CPT+2.0*CST)
      V4X=WN*SM1*( A*CPT-2.0*CST)
C
C   * DISP-STRESS VECTOR DUE TO Qz *
C
      CPT=EX*(CUP1Z-WHZ1*CUP2Z)
      CST=EX*(CUS1Z-WZ1Z*CUP1Z-WHZ1*CUS2Z+WZ1Z*WHZ1*CUP2Z)
C
      V1Z=CPT+CST
      V2Z=A*CPT-CST
      V3Z=WN*SM1*(-B*CPT+2.0*CST)
      V4Z=WN*SM1*( A*CPT-2.0*CST)
C
      GOTO 50
C
C *** DISP-STRESS VECTOR DUE TO Qx FOR DEP1 < (OR =) ODP < SDP ***
C
   10 IF(ODP.GE.SDP)GOTO 20
C
C   * UP-GOING WAVE DUE TO DIRECT WAVE *
C
      WHZ=WN*(SDP-ODP)
      EXD=DEXP(-WHZ)
C
      CPD=EXD*CUPX
      CSD=EXD*(CUSX-WHZ*CUPX)
C
C   * DOWN-GOING WAVE DUE TO REFLECT WAVE FROM UPPER BOUNDARY *
C
      WHZ1=WN*(SDP-DEP1)
      WZZ1=WN*(ODP-DEP1)
      EXR1=DEXP(-WHZ1-WZZ1)
C
      CPR1=EXR1*(CDP1X-WHZ1*CDP2X)
      CSR1=EXR1*(CDS1X-WHZ1*CDS2X-WZZ1*CDP1X+WZZ1*WHZ1*CDP2X)
C
C   * UP-GOING WAVE DUE TO REFLECT WAVE FROM LOWER BOUNDARY *
C
      IF(LS.EQ.NL)THEN
        CPR2=(0.D0,0.D0)
        CSR2=(0.D0,0.D0)
      ELSE
        WZ2H=WN*(DEP2-SDP)
        WZ2Z=WN*(DEP2-ODP)
        EXR2=DEXP(-WZ2H-WZ2Z)
        CPR2=EXR2*(CUP1X-WZ2H*CUP2X)
        CSR2=EXR2*(CUS1X-WZ2H*CUS2X-WZ2Z*CUP1X+WZ2Z*WZ2H*CUP2X)
      END IF
C
C   * DISP-STRESS VECTORS *
C
      H1X=H1D*EXD+H1R1*EXR1+H1R2*EXR2
      H2X=WN*(H2D*EXD+H2R1*EXR1+H2R2*EXR2)
C
      A=X2-1.0
      B=X2-3.0
      V1X=CPR1+CSR1+CPD+CSD+CPR2+CSR2
      V2X=-A*CPR1+CSR1+A*(CPD+CPR2)-CSD-CSR2
      V3X=WN*SM2*(B*CPR1-2.0*CSR1-B*(CPD+CPR2)+2.0*(CSD+CSR2))
      V4X=WN*SM2*(A*CPR1-2.0*CSR1+A*(CPD+CPR2)-2.0*(CSD+CSR2))
C
C   * DISP-STRESS VECTOR DUE TO Qz *
C
      CPD=EXD*CUPZ
      CSD=EXD*(CUSZ-WHZ*CUPZ)
C
      CPR1=EXR1*(CDP1Z-WHZ1*CDP2Z)
      CSR1=EXR1*(CDS1Z-WHZ1*CDS2Z-WZZ1*CDP1Z+WZZ1*WHZ1*CDP2Z)
C
      IF(LS.EQ.NL)THEN
        CPR2=(0.D0,0.D0)
        CSR2=(0.D0,0.D0)
      ELSE
        CPR2=EXR2*(CUP1Z-WZ2H*CUP2Z)
        CSR2=EXR2*(CUS1Z-WZ2H*CUS2Z-WZ2Z*CUP1Z+WZ2Z*WZ2H*CUP2Z)
      END IF
C
      V1Z=CPR1+CSR1+CPD+CSD+CPR2+CSR2
      V2Z=-A*CPR1+CSR1+A*(CPD+CPR2)-CSD-CSR2
      V3Z=WN*SM2*(B*CPR1-2.0*CSR1-B*(CPD+CPR2)+2.0*(CSD+CSR2))
      V4Z=WN*SM2*(A*CPR1-2.0*CSR1+A*(CPD+CPR2)-2.0*(CSD+CSR2))
C
      GOTO 50
C
C *** DISP-STRESS VECTOR DUE TO Qx FOR SDP < (OR =) ODP < DEP2 ***
C
   20 IF((LS.NE.NL).AND.(ODP.GE.DEP2))GOTO 30
C
C   * DOWN-GOING WAVE DUE TO DIRECT WAVE *
C
      WZH=WN*(ODP-SDP)
      EXD=DEXP(-WZH)
C
      CPD=EXD*CDPX
      CSD=EXD*(CDSX-WZH*CDPX)
C
C   * UP-GOING WAVE DUE TO REFLECT WAVE FROM LOWER BOUNDARY *
C
      IF(LS.EQ.NL)THEN
        CPR2=(0.D0,0.D0)
        CSR2=(0.D0,0.D0)
      ELSE
        WZ2H=WN*(DEP2-SDP)
        WZ2Z=WN*(DEP2-ODP)
        EXR2=DEXP(-WZ2H-WZ2Z)
        CPR2=EXR2*(CUP1X-WZ2H*CUP2X)
        CSR2=EXR2*(CUS1X-WZ2H*CUS2X-WZ2Z*CUP1X+WZ2Z*WZ2H*CUP2X)
      END IF
C
C   * DOWN-GOING WAVE DUE TO REFLECT WAVE FROM UPPER BOUNDARY *
C
      WHZ1=WN*(SDP-DEP1)
      WZZ1=WN*(ODP-DEP1)
      EXR1=DEXP(-WHZ1-WZZ1)
C
      CPR1=EXR1*(CDP1X-WHZ1*CDP2X)
      CSR1=EXR1*(CDS1X-WHZ1*CDS2X-WZZ1*CDP1X+WZZ1*WHZ1*CDP2X)
C
C   * DISP-STRESS VECTORS *
C
      H1X=H1D*EXD+H1R1*EXR1+H1R2*EXR2
      H2X=WN*(-H2D*EXD+H2R1*EXR1+H2R2*EXR2)
C
      A=X2-1.0
      B=X2-3.0
      V1X=CPD+CSD+CPR1+CSR1+CPR2+CSR2
      V2X=-A*(CPD+CPR1)+CSD+CSR1+A*CPR2-CSR2
      V3X=WN*SM2*(B*(CPD+CPR1)-2.0*(CSD+CSR1)-B*CPR2+2.0*CSR2)
      V4X=WN*SM2*(A*(CPD+CPR1)-2.0*(CSD+CSR1)+A*CPR2-2.0*CSR2)
C
C   * DISP-STRESS VECTOR DUE TO Qz *
C
      CPD=EXD*CDPZ
      CSD=EXD*(CDSZ-WZH*CDPZ)
C
      IF(LS.EQ.NL)THEN
        CPR2=(0.D0,0.D0)
        CSR2=(0.D0,0.D0)
      ELSE
        CPD=EXD*CDPZ
        CSD=EXD*(CDSZ-WZH*CDPZ)
        CPR2=EXR2*(CUP1Z-WZ2H*CUP2Z)
        CSR2=EXR2*(CUS1Z-WZ2H*CUS2Z-WZ2Z*CUP1Z+WZ2Z*WZ2H*CUP2Z)
      END IF
C
      CPR1=EXR1*(CDP1Z-WHZ1*CDP2Z)
      CSR1=EXR1*(CDS1Z-WHZ1*CDS2Z-WZZ1*CDP1Z+WZZ1*WHZ1*CDP2Z)
C
      V1Z=CPD+CSD+CPR1+CSR1+CPR2+CSR2
      V2Z=-A*(CPD+CPR1)+CSD+CSR1+A*CPR2-CSR2
      V3Z=WN*SM2*(B*(CPD+CPR1)-2.0*(CSD+CSR1)-B*CPR2+2.0*CSR2)
      V4Z=WN*SM2*(A*(CPD+CPR1)-2.0*(CSD+CSR1)+A*CPR2-2.0*CSR2)
C
      GOTO 50
C
C *** DISP-STRESS VECTOR DUE TO Qx FOR DEP2 < (OR =) ODP ***
C
   30 WZH=WN*(ODP-SDP)
      WZ2H=WN*(DEP2-SDP)
      WZZ2=WN*(ODP-DEP2)
      EX=DEXP(-WZH)
C
      CPT=EX*(CDP1X-WZ2H*CDP2X)
      CST=EX*(CDS1X-WZZ2*CDP1X-WZ2H*CDS2X+WZZ2*WZ2H*CDP2X)
C
      H1X=EX*H1T
      H2X=WN*EX*H2T
C
      A=X3-1.0
      B=X3-3.0
      V1X=CPT+CST
      V2X=-A*CPT+CST
      V3X=WN*SM3*(B*CPT-2.0*CST)
      V4X=WN*SM3*(A*CPT-2.0*CST)
C
C   * DISP-STRESS VECTOR DUE TO Qz *
C
      CPT=EX*(CDP1Z-WZ2H*CDP2Z)
      CST=EX*(CDS1Z-WZZ2*CDP1Z-WZ2H*CDS2Z+WZZ2*WZ2H*CDP2Z)
C
      V1Z=CPT+CST
      V2Z=-A*CPT+CST
      V3Z=WN*SM3*(B*CPT-2.0*CST)
      V4Z=WN*SM3*(A*CPT-2.0*CST)
C
   50 RETURN
      END
C
C
C *** ANALITICAL INTEGRATIONS OF ASYMPTOTIC SOLUTIONS ***
C
      SUBROUTINE AINTRZ(LS,NL,QXS,QZS,VXS,WXS,WZS)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 QXS(0:2),QZS(0:2),VXS(0:2),WXS(0:2),WZS(0:2)
     *   ,SM1,X1,SM2,X2,SM3,X3,RM1,RM2,RM3,SM,RM,DC1,DC2,DC3
     *   ,H1T,H2T,H1D,H2D,H1R1,H2R1,H1R2,H2R2       
     *   ,CUPX,CUSX,CUP1X,CUS1X,CUP2X,CUS2X
     *   ,CDPX,CDSX,CDP1X,CDS1X,CDP2X,CDS2X
     *   ,CUPZ,CUSZ,CUP1Z,CUS1Z,CUP2Z,CUS2Z
     *   ,CDPZ,CDSZ,CDP1Z,CDS1Z,CDP2Z,CDS2Z
     *   ,DCUPT,DCUST,DCUPD,DCUSD,DCUPR,DCUSR
     *   ,DCDPT,DCDST,DCDPD,DCDSD,DCDPR,DCDSR
     *   ,D1MH1,DD0H1,D11H1,D1MH2,DD0H2,D1MH1R,DD0H1R,D1MH1Z,DD0H1Z
     *   ,D1MV1,DD0V1,D11V1,D01V2,D10V2,D1MV3,DD0V3,D10V4
     *   ,D1MV1R,DD0V1R,D10V2R,D1MV1Z,DD0V1Z,D10V2Z
     *   ,D01V1,D10V1,D00V2,D11V2,D10V3,D00V4
     *   ,D10V1R,D00V2R,D10V1Z,D00V2Z
      COMMON/ASYMP/SM1,X1,SM2,X2,SM3,X3 /RAM/RM1,RM2,RM3
     *      /DPSO /DEP1,DEP2,SDP,ODP /FCQVW/SM,RM,R
     *      /HX  /H1T,H2T,H1D,H2D,H1R1,H2R1,H1R2,H2R2       
     *      /VUX /CUPX,CUSX,CUP1X,CUS1X,CUP2X,CUS2X
     *      /VDX /CDPX,CDSX,CDP1X,CDS1X,CDP2X,CDS2X
     *      /VUZ /CUPZ,CUSZ,CUP1Z,CUS1Z,CUP2Z,CUS2Z
     *      /VDZ /CDPZ,CDSZ,CDP1Z,CDS1Z,CDP2Z,CDS2Z
C
C *** ANALITICAL INTEGRATIONS ***
C
C   * INTEGRATIONS FOR DIRECT WAVE *
C
      ZH=DABS(SDP-ODP)
      ZH2=ZH*ZH
      RR=DSQRT(R*R+ZH2)
      R2=RR*RR
      R3=R2*RR
      R5=R3*R2
C
      S00=1.D0/RR
      S01=ZH/R3
      S02=(3.D0*ZH2/R2-1.D0)/R3
      S03=3.D0*ZH*(5.D0*ZH2/R2-3.D0)/R5
C
      S1M=(RR-ZH)/R
      S10=S1M/RR
      S11=R/R3
      S12=3.D0*R*ZH/R5
      S13=3.D0*R*(5.D0*ZH2/R2-1.D0)/R5
C
C   * INTEGRATIONS FOR REFLECTED WAVE FROM UPPER BOUNDARY *
C
      ZH=SDP+ODP-2.D0*DEP1
      ZH2=ZH*ZH
      RR=DSQRT(R*R+ZH2)
      R2=RR*RR
      R3=R2*RR
      R5=R3*R2
C
      S001=1.D0/RR
      S011=ZH/R3
      S021=(3.D0*ZH2/R2-1.D0)/R3
      S031=3.D0*ZH*(5.D0*ZH2/R2-3.D0)/R5
C
      S1M1=(RR-ZH)/R
      S101=S1M1/RR
      S111=R/R3
      S121=3.D0*R*ZH/R5
      S131=3.D0*R*(5.D0*ZH2/R2-1.D0)/R5
C
C   * INTEGRATIONS FOR REFLECTED WAVE FROM LOWER BOUNDARY *
C
      ZH=2.D0*DEP2-SDP-ODP
      ZH2=ZH*ZH
      RR=DSQRT(R*R+ZH2)
      R2=RR*RR
      R3=R2*RR
      R5=R3*R2
C
      S002=1.D0/RR
      S012=ZH/R3
      S022=(3.D0*ZH2/R2-1.D0)/R3
      S032=3.D0*ZH*(5.D0*ZH2/R2-3.D0)/R5
C
      S1M2=(RR-ZH)/R
      S102=S1M2/RR
      S112=R/R3
      S122=3.D0*R*ZH/R5
      S132=3.D0*R*(5.D0*ZH2/R2-1.D0)/R5
C
C *** ODP < DEP1 ***
C
      IF(ODP.GE.DEP1) GOTO 10
      IF(LS.EQ.1) GOTO 10
C
C  ** QX **
C
C   * SH (TRANSMITTED) *
C
      D1MH1=H1T*S1M
      DD0H1=H1T*(S00-S1M/R)
      QXS(0)=D1MH1/R
      VXS(0)=DD0H1
C
      D11H1=H1T*S11
      D1MH1R=DD0H1
      DD0H1R=(-DD0H1+D1MH1/R)/R-D11H1
      QXS(1)=(-D1MH1/R+D1MH1R)/R
      VXS(1)=DD0H1R
C
      D1MH2=H2T*S10
      DD0H2=H2T*(S01-S10/R)
      D1MH1Z=D1MH2/SM1
      DD0H1Z=DD0H2/SM1
      QXS(2)=D1MH1Z/R
      VXS(2)=DD0H1Z
C
C   * P-SV (TRANSMITTED) *
C
      HD1=SDP-DEP1
      D1Z=DEP1-ODP
      DC1=HD1*CUP2X
      DC2=HD1*CUS2X+D1Z*CUP1X
      DC3=D1Z*DC1
      DCUPT=CUP1X*S1M-DC1*S10
      DCUST=CUS1X*S1M-DC2*S10+DC3*S11
      D1MV1=DCUPT+DCUST
C
      DCUPT=CUP1X*S11-DC1*S12
      DCUST=CUS1X*S11-DC2*S12+DC3*S13
      D11V1=DCUPT+DCUST
      D10V4=SM1*((X1-1.0)*DCUPT-2.0*DCUST)
C
      S1=S00-S1M/R
      S2=S01-S10/R
      S3=S02-S11/R
      DCUPT=CUP1X*S1-DC1*S2
      DCUST=CUS1X*S1-DC2*S2+DC3*S3
      DD0V1=DCUPT+DCUST
C
      DCUPT=CUP1X*S01-DC1*S02
      DCUST=CUS1X*S01-DC2*S02+DC3*S03
      D01V2=(X1-1.0)*DCUPT-DCUST
C
      DCUPT=CUP1X*S10-DC1*S11
      DCUST=CUS1X*S10-DC2*S11+DC3*S12
      D10V2=(X1-1.0)*DCUPT-DCUST
      D1MV3=SM1*(-(X1-3.0)*DCUPT+2.0*DCUST)
C
      S1=S01-S10/R
      S2=S02-S11/R
      S3=S03-S12/R
      DCUPT=CUP1X*S1-DC1*S2
      DCUST=CUS1X*S1-DC2*S2+DC3*S3
      DD0V3=SM1*(-(X1-3.0)*DCUPT+2.0*DCUST)
C
      QXS(0)=QXS(0)+DD0V1
      VXS(0)=VXS(0)+D1MV1/R
      WXS(0)=D10V2
C
      D1MV1R=DD0V1
      DD0V1R=(-DD0V1+D1MV1/R)/R-D11V1
      D10V2R=D01V2-D10V2/R
C
      QXS(1)=QXS(1)+DD0V1R
      VXS(1)=VXS(1)+(-D1MV1/R+D1MV1R)/R
      WXS(1)=D10V2R
C
      D1MV1Z=D1MV3/SM1+D10V2
      DD0V1Z=DD0V3/SM1+D01V2-D10V2/R
      D10V2Z=(D10V4-RM1*D11V1)/(RM1+2.D0*SM1) 
C
      QXS(2)=QXS(2)+DD0V1Z
      VXS(2)=VXS(2)+D1MV1Z/R
      WXS(2)=D10V2Z
C
C  ** QZ **
C
C   * P-SV (TRANSMITTED) *
C
      DC1=HD1*CUP2Z
      DC2=HD1*CUS2Z+D1Z*CUP1Z
      DC3=D1Z*DC1
      DCUPT=CUP1Z*S01-DC1*S02
      DCUST=CUS1Z*S01-DC2*S02+DC3*S03
      D01V1=DCUPT+DCUST
      D00V4=SM1*((X1-1.0)*DCUPT-2.0*DCUST)
C
      DCUPT=CUP1Z*S10-DC1*S11
      DCUST=CUS1Z*S10-DC2*S11+DC3*S12
      D10V1=DCUPT+DCUST
C
      DCUPT=CUP1Z*S00-DC1*S01
      DCUST=CUS1Z*S00-DC2*S01+DC3*S02
      D00V2=(X1-1.0)*DCUPT-DCUST
C
      DCUPT=CUP1Z*S11-DC1*S12
      DCUST=CUS1Z*S11-DC2*S12+DC3*S13
      D11V2=(X1-1.0)*DCUPT-DCUST
      D10V3=SM1*(-(X1-3.0)*DCUPT+2.0*DCUST)
C
      QZS(0)=D10V1
      WZS(0)=D00V2
C
      D10V1R=D01V1-D10V1/R
      D00V2R=-D11V2
C
      QZS(1)=D10V1R
      WZS(1)=D00V2R
C
      D10V1Z=D10V3/SM1+D11V2
      D00V2Z=(D00V4-RM1*D01V1)/(RM1+2.D0*SM1)
C
      QZS(2)=D10V1Z
      WZS(2)=D00V2Z
C
      GOTO 50
C
C *** DEP1 < (OR =) ODP < SDP ***
C
   10 IF(ODP.GE.SDP) GOTO 20
C
C  ** QX **
C
C   * SH *
C
      D1MH1=H1D*S1M+H1R1*S1M1+H1R2*S1M2
      DD0H1=H1D*(S00-S1M/R)+H1R1*(S001-S1M1/R)+H1R2*(S002-S1M2/R)
      QXS(0)=D1MH1/R
      VXS(0)=DD0H1
C
      D11H1=H1D*S11+H1R1*S111+H1R2*S112
      D1MH1R=DD0H1
      DD0H1R=(-DD0H1+D1MH1/R)/R-D11H1
      QXS(1)=(-D1MH1/R+D1MH1R)/R
      VXS(1)=DD0H1R
C
      D1MH2=H2D*S10+H2R1*S101+H2R2*S102
      DD0H2=H2D*(S01-S10/R)+H2R1*(S011-S101/R)+H2R2*(S012-S102/R)
      D1MH1Z=D1MH2/SM2
      DD0H1Z=DD0H2/SM2
      QXS(2)=D1MH1Z/R
      VXS(2)=DD0H1Z
C
C   * P-SV (UP-GOING DIRECT) *
C
      DC1=(SDP-ODP)*CUPX
      DCUPD=CUPX*S1M
      DCUSD=CUSX*S1M-DC1*S10
      D1MV1=DCUPD+DCUSD
C
      DCUPD=CUPX*S11
      DCUSD=CUSX*S11-DC1*S12
      D11V1=DCUPD+DCUSD
      D10V4=SM2*((X2-1.0)*DCUPD-2.0*DCUSD)
C
      S1=S00-S1M/R
      S2=S01-S10/R
      DCUPD=CUPX*S1
      DCUSD=CUSX*S1-DC1*S2
      DD0V1=DCUPD+DCUSD
C
      DCUPD=CUPX*S01
      DCUSD=CUSX*S01-DC1*S02
      D01V2=(X2-1.0)*DCUPD-DCUSD
C
      DCUPD=CUPX*S10
      DCUSD=CUSX*S10-DC1*S11
      D10V2=(X2-1.0)*DCUPD-DCUSD
      D1MV3=SM2*(-(X2-3.0)*DCUPD+2.0*DCUSD)
C
      S1=S01-S10/R
      S2=S02-S11/R
      DCUPD=CUPX*S1
      DCUSD=CUSX*S1-DC1*S2
      DD0V3=SM2*(-(X2-3.0)*DCUPD+2.0*DCUSD)
C
      QXS(0)=QXS(0)+DD0V1
      VXS(0)=VXS(0)+D1MV1/R
      WXS(0)=D10V2
C
      D1MV1R=DD0V1
      DD0V1R=(-DD0V1+D1MV1/R)/R-D11V1
      D10V2R=D01V2-D10V2/R
C
      QXS(1)=QXS(1)+DD0V1R
      VXS(1)=VXS(1)+(-D1MV1/R+D1MV1R)/R
      WXS(1)=D10V2R
C
      D1MV1Z=D1MV3/SM2+D10V2
      DD0V1Z=DD0V3/SM2+D01V2-D10V2/R
      D10V2Z=(D10V4-RM2*D11V1)/(RM2+2.D0*SM2) 
C
      QXS(2)=QXS(2)+DD0V1Z
      VXS(2)=VXS(2)+D1MV1Z/R
      WXS(2)=D10V2Z
C
C   * P-SV (DOWN-GOING REFLECTED FROM UPPER BOUNDARY) *
C
      HD1=SDP-DEP1
      ZD1=ODP-DEP1
      DC1=HD1*CDP2X
      DC2=HD1*CDS2X+ZD1*CDP1X
      DC3=ZD1*DC1
      DCDPR=CDP1X*S1M1-DC1*S101
      DCDSR=CDS1X*S1M1-DC2*S101+DC3*S111
      D1MV1=DCDPR+DCDSR
C
      DCDPR=CDP1X*S111-DC1*S121
      DCDSR=CDS1X*S111-DC2*S121+DC3*S131
      D11V1=DCDPR+DCDSR
      D10V4=SM2*((X2-1.0)*DCDPR-2.0*DCDSR)
C
      S1=S001-S1M1/R
      S2=S011-S101/R
      S3=S021-S111/R
      DCDPR=CDP1X*S1-DC1*S2
      DCDSR=CDS1X*S1-DC2*S2+DC3*S3
      DD0V1=DCDPR+DCDSR
C
      DCDPR=CDP1X*S011-DC1*S021
      DCDSR=CDS1X*S011-DC2*S021+DC3*S031
      D01V2=-(X2-1.0)*DCDPR+DCDSR
C
      DCDPR=CDP1X*S101-DC1*S111
      DCDSR=CDS1X*S101-DC2*S111+DC3*S121
      D10V2=-(X2-1.0)*DCDPR+DCDSR
      D1MV3=SM2*((X2-3.0)*DCDPR-2.0*DCDSR)
C
      S1=S011-S101/R
      S2=S021-S111/R
      S3=S031-S121/R
      DCDPR=CDP1X*S1-DC1*S2
      DCDSR=CDS1X*S1-DC2*S2+DC3*S3
      DD0V3=SM2*((X2-3.0)*DCDPR-2.0*DCDSR)
C
      QXS(0)=QXS(0)+DD0V1
      VXS(0)=VXS(0)+D1MV1/R
      WXS(0)=WXS(0)+D10V2
C
      D1MV1R=DD0V1
      DD0V1R=(-DD0V1+D1MV1/R)/R-D11V1
      D10V2R=D01V2-D10V2/R
C
      QXS(1)=QXS(1)+DD0V1R
      VXS(1)=VXS(1)+(-D1MV1/R+D1MV1R)/R
      WXS(1)=WXS(1)+D10V2R
C
      D1MV1Z=D1MV3/SM2+D10V2
      DD0V1Z=DD0V3/SM2+D01V2-D10V2/R
      D10V2Z=(D10V4-RM2*D11V1)/(RM2+2.D0*SM2) 
C
      QXS(2)=QXS(2)+DD0V1Z
      VXS(2)=VXS(2)+D1MV1Z/R
      WXS(2)=WXS(2)+D10V2Z
C
C   * P-SV (UP-GOING REFLECTED FROM LOWER BOUNDARY) *
C
      IF(LS.EQ.NL)GOTO 15
C
      D2H=DEP2-SDP
      D2Z=DEP2-ODP
      DC1=D2H*CUP2X
      DC2=D2H*CUS2X+D2Z*CUP1X
      DC3=D2Z*DC1
      DCUPR=CUP1X*S1M2-DC1*S102
      DCUSR=CUS1X*S1M2-DC2*S102+DC3*S112
      D1MV1=DCUPR+DCUSR
C
      DCUPR=CUP1X*S112-DC1*S122
      DCUSR=CUS1X*S112-DC2*S122+DC3*S132
      D11V1=DCUPR+DCUSR
      D10V4=SM2*((X2-1.0)*DCUPR-2.0*DCUSR)
C
      S1=S002-S1M2/R
      S2=S012-S102/R
      S3=S022-S112/R
      DCUPR=CUP1X*S1-DC1*S2
      DCUSR=CUS1X*S1-DC2*S2+DC3*S3
      DD0V1=DCUPR+DCUSR
C
      DCUPR=CUP1X*S012-DC1*S022
      DCUSR=CUS1X*S012-DC2*S022+DC3*S032
      D01V2=(X2-1.0)*DCUPR-DCUSR
C
      DCUPR=CUP1X*S102-DC1*S112
      DCUSR=CUS1X*S102-DC2*S112+DC3*S122
      D10V2=(X2-1.0)*DCUPR-DCUSR
      D1MV3=SM2*(-(X2-3.0)*DCUPR+2.0*DCUSR)
C
      S1=S012-S102/R
      S2=S022-S112/R
      S3=S032-S122/R
      DCUPR=CUP1X*S1-DC1*S2
      DCUSR=CUS1X*S1-DC2*S2+DC3*S3
      DD0V3=SM2*(-(X2-3.0)*DCUPR+2.0*DCUSR)
C
      QXS(0)=QXS(0)+DD0V1
      VXS(0)=VXS(0)+D1MV1/R
      WXS(0)=WXS(0)+D10V2
C
      D1MV1R=DD0V1
      DD0V1R=(-DD0V1+D1MV1/R)/R-D11V1
      D10V2R=D01V2-D10V2/R
C
      QXS(1)=QXS(1)+DD0V1R
      VXS(1)=VXS(1)+(-D1MV1/R+D1MV1R)/R
      WXS(1)=WXS(1)+D10V2R
C
      D1MV1Z=D1MV3/SM2+D10V2
      DD0V1Z=DD0V3/SM2+D01V2-D10V2/R
      D10V2Z=(D10V4-RM2*D11V1)/(RM2+2.D0*SM2) 
C
      QXS(2)=QXS(2)+DD0V1Z
      VXS(2)=VXS(2)+D1MV1Z/R
      WXS(2)=WXS(2)+D10V2Z
C
C  ** QZ **
C
   15 DC1=(SDP-ODP)*CUPZ
      DCUPD=CUPZ*S01
      DCUSD=CUSZ*S01-DC1*S02
      D01V1=DCUPD+DCUSD
      D00V4=SM2*((X2-1.0)*DCUPD-2.0*DCUSD)
C
      DCUPD=CUPZ*S10
      DCUSD=CUSZ*S10-DC1*S11
      D10V1=DCUPD+DCUSD
C
      DCUPD=CUPZ*S00
      DCUSD=CUSZ*S00-DC1*S01
      D00V2=(X2-1.0)*DCUPD-DCUSD
C
      DCUPD=CUPZ*S11
      DCUSD=CUSZ*S11-DC1*S12
      D11V2=(X2-1.0)*DCUPD-DCUSD
      D10V3=SM2*(-(X2-3.0)*DCUPD+2.0*DCUSD)
C
      QZS(0)=D10V1
      WZS(0)=D00V2
C
      D10V1R=D01V1-D10V1/R
      D00V2R=-D11V2
C
      QZS(1)=D10V1R
      WZS(1)=D00V2R
C
      D10V1Z=D10V3/SM2+D11V2
      D00V2Z=(D00V4-RM2*D01V1)/(RM2+2.D0*SM2)
C
      QZS(2)=D10V1Z
      WZS(2)=D00V2Z
C
      DC1=HD1*CDP2Z
      DC2=HD1*CDS2Z+ZD1*CDP1Z
      DC3=ZD1*DC1
      DCDPR=CDP1Z*S011-DC1*S021
      DCDSR=CDS1Z*S011-DC2*S021+DC3*S031
      D01V1=DCDPR+DCDSR
      D00V4=SM2*((X2-1.0)*DCDPR-2.0*DCDSR)
C
      DCDPR=CDP1Z*S101-DC1*S111
      DCDSR=CDS1Z*S101-DC2*S111+DC3*S121
      D10V1=DCDPR+DCDSR
C
      DCDPR=CDP1Z*S001-DC1*S011
      DCDSR=CDS1Z*S001-DC2*S011+DC3*S021
      D00V2=-(X2-1.0)*DCDPR+DCDSR
C
      DCDPR=CDP1Z*S111-DC1*S121
      DCDSR=CDS1Z*S111-DC2*S121+DC3*S131
      D11V2=-(X2-1.0)*DCDPR+DCDSR
      D10V3=SM2*((X2-3.0)*DCDPR-2.0*DCDSR)
C
      QZS(0)=QZS(0)+D10V1
      WZS(0)=WZS(0)+D00V2
C
      D10V1R=D01V1-D10V1/R
      D00V2R=-D11V2
C
      QZS(1)=QZS(1)+D10V1R
      WZS(1)=WZS(1)+D00V2R
C
      D10V1Z=D10V3/SM2+D11V2
      D00V2Z=(D00V4-RM2*D01V1)/(RM2+2.D0*SM2)
C
      QZS(2)=QZS(2)+D10V1Z
      WZS(2)=WZS(2)+D00V2Z
C
      IF(LS.EQ.NL)GOTO 50
C
      DC1=D2H*CUP2Z
      DC2=D2H*CUS2Z+D2Z*CUP1Z
      DC3=D2Z*DC1
      DCUPR=CUP1Z*S012-DC1*S022
      DCUSR=CUS1Z*S012-DC2*S022+DC3*S032
      D01V1=DCUPR+DCUSR
      D00V4=SM2*((X2-1.0)*DCUPR-2.0*DCUSR)
C
      DCUPR=CUP1Z*S102-DC1*S112
      DCUSR=CUS1Z*S102-DC2*S112+DC3*S122
      D10V1=DCUPR+DCUSR
C
      DCUPR=CUP1Z*S002-DC1*S012
      DCUSR=CUS1Z*S002-DC2*S012+DC3*S022
      D00V2=(X2-1.0)*DCUPR-DCUSR
C
      DCUPR=CUP1Z*S112-DC1*S122
      DCUSR=CUS1Z*S112-DC2*S122+DC3*S132
      D11V2=(X2-1.0)*DCUPR-DCUSR
      D10V3=SM2*(-(X2-3.0)*DCUPR+2.0*DCUSR)
C
      QZS(0)=QZS(0)+D10V1
      WZS(0)=WZS(0)+D00V2
C
      D10V1R=D01V1-D10V1/R
      D00V2R=-D11V2
C
      QZS(1)=QZS(1)+D10V1R
      WZS(1)=WZS(1)+D00V2R
C
      D10V1Z=D10V3/SM2+D11V2
      D00V2Z=(D00V4-RM2*D01V1)/(RM2+2.D0*SM2)
C
      QZS(2)=QZS(2)+D10V1Z
      WZS(2)=WZS(2)+D00V2Z
C
      GOTO 50
C
C *** SDP < (OR =) ODP < DEP2 ***
C
   20 IF((LS.NE.NL).AND.(ODP.GE.DEP2)) GOTO 30
C
C  ** QX **
C
C   * SH *
C
      D1MH1=H1D*S1M+H1R1*S1M1+H1R2*S1M2
      DD0H1=H1D*(S00-S1M/R)+H1R1*(S001-S1M1/R)+H1R2*(S002-S1M2/R)
      QXS(0)=D1MH1/R
      VXS(0)=DD0H1
C
      D11H1=H1D*S11+H1R1*S111+H1R2*S112
      D1MH1R=DD0H1
      DD0H1R=(-DD0H1+D1MH1/R)/R-D11H1
      QXS(1)=(-D1MH1/R+D1MH1R)/R
      VXS(1)=DD0H1R
C
      D1MH2=-H2D*S10+H2R1*S101+H2R2*S102
      DD0H2=-H2D*(S01-S10/R)+H2R1*(S011-S101/R)+H2R2*(S012-S102/R)
      D1MH1Z=D1MH2/SM2
      DD0H1Z=DD0H2/SM2
      QXS(2)=D1MH1Z/R
      VXS(2)=DD0H1Z
C
C   * P-SV (DOWN-GOING DIRECT) *
C
      DC1=(ODP-SDP)*CDPX
      DCDPD=CDPX*S1M
      DCDSD=CDSX*S1M-DC1*S10
      D1MV1=DCDPD+DCDSD
C
      DCDPD=CDPX*S11
      DCDSD=CDSX*S11-DC1*S12
      D11V1=DCDPD+DCDSD
      D10V4=SM2*((X2-1.0)*DCDPD-2.0*DCDSD)
C
      S1=S00-S1M/R
      S2=S01-S10/R
      DCDPD=CDPX*S1
      DCDSD=CDSX*S1-DC1*S2
      DD0V1=DCDPD+DCDSD
C
      DCDPD=CDPX*S01
      DCDSD=CDSX*S01-DC1*S02
      D01V2=-(X2-1.0)*DCDPD+DCDSD
C
      DCDPD=CDPX*S10
      DCDSD=CDSX*S10-DC1*S11
      D10V2=-(X2-1.0)*DCDPD+DCDSD
      D1MV3=SM2*((X2-3.0)*DCDPD-2.0*DCDSD)
C
      S1=S01-S10/R
      S2=S02-S11/R
      DCDPD=CDPX*S1
      DCDSD=CDSX*S1-DC1*S2
      DD0V3=SM2*((X2-3.0)*DCDPD-2.0*DCDSD)
C
      QXS(0)=QXS(0)+DD0V1
      VXS(0)=VXS(0)+D1MV1/R
      WXS(0)=D10V2
C
      D1MV1R=DD0V1
      DD0V1R=(-DD0V1+D1MV1/R)/R-D11V1
      D10V2R=D01V2-D10V2/R
C
      QXS(1)=QXS(1)+DD0V1R
      VXS(1)=VXS(1)+(-D1MV1/R+D1MV1R)/R
      WXS(1)=D10V2R
C
      D1MV1Z=D1MV3/SM2+D10V2
      DD0V1Z=DD0V3/SM2+D01V2-D10V2/R
      D10V2Z=(D10V4-RM2*D11V1)/(RM2+2.D0*SM2) 
C
      QXS(2)=QXS(2)+DD0V1Z
      VXS(2)=VXS(2)+D1MV1Z/R
      WXS(2)=D10V2Z
C
C   * P-SV (UP-GOING REFLECTED FROM LOWER BOUNDARY) *
C
      IF(LS.EQ.NL)GOTO 25
C
      D2H=DEP2-SDP
      D2Z=DEP2-ODP
      DC1=D2H*CUP2X
      DC2=D2H*CUS2X+D2Z*CUP1X
      DC3=D2Z*DC1
      DCUPR=CUP1X*S1M2-DC1*S102
      DCUSR=CUS1X*S1M2-DC2*S102+DC3*S112
      D1MV1=DCUPR+DCUSR
C
      DCUPR=CUP1X*S112-DC1*S122
      DCUSR=CUS1X*S112-DC2*S122+DC3*S132
      D11V1=DCUPR+DCUSR
      D10V4=SM2*((X2-1.0)*DCUPR-2.0*DCUSR)
C
      S1=S002-S1M2/R
      S2=S012-S102/R
      S3=S022-S112/R
      DCUPR=CUP1X*S1-DC1*S2
      DCUSR=CUS1X*S1-DC2*S2+DC3*S3
      DD0V1=DCUPR+DCUSR
C
      DCUPR=CUP1X*S012-DC1*S022
      DCUSR=CUS1X*S012-DC2*S022+DC3*S032
      D01V2=(X2-1.0)*DCUPR-DCUSR
C
      DCUPR=CUP1X*S102-DC1*S112
      DCUSR=CUS1X*S102-DC2*S112+DC3*S122
      D10V2=(X2-1.0)*DCUPR-DCUSR
      D1MV3=SM2*(-(X2-3.0)*DCUPR+2.0*DCUSR)
C
      S1=S012-S102/R
      S2=S022-S112/R
      S3=S032-S122/R
      DCUPR=CUP1X*S1-DC1*S2
      DCUSR=CUS1X*S1-DC2*S2+DC3*S3
      DD0V3=SM2*(-(X2-3.0)*DCUPR+2.0*DCUSR)
C
      QXS(0)=QXS(0)+DD0V1
      VXS(0)=VXS(0)+D1MV1/R
      WXS(0)=WXS(0)+D10V2
C
      D1MV1R=DD0V1
      DD0V1R=(-DD0V1+D1MV1/R)/R-D11V1
      D10V2R=D01V2-D10V2/R
C
      QXS(1)=QXS(1)+DD0V1R
      VXS(1)=VXS(1)+(-D1MV1/R+D1MV1R)/R
      WXS(1)=WXS(1)+D10V2R
C
      D1MV1Z=D1MV3/SM2+D10V2
      DD0V1Z=DD0V3/SM2+D01V2-D10V2/R
      D10V2Z=(D10V4-RM2*D11V1)/(RM2+2.D0*SM2) 
C
      QXS(2)=QXS(2)+DD0V1Z
      VXS(2)=VXS(2)+D1MV1Z/R
      WXS(2)=WXS(2)+D10V2Z
C
C   * P-SV (DOWN-GOING REFLECTED FROM UPPER BOUNDARY) *
C
   25 HD1=SDP-DEP1
      ZD1=ODP-DEP1
      DC1=HD1*CDP2X
      DC2=HD1*CDS2X+ZD1*CDP1X
      DC3=ZD1*DC1
      DCDPR=CDP1X*S1M1-DC1*S101
      DCDSR=CDS1X*S1M1-DC2*S101+DC3*S111
      D1MV1=DCDPR+DCDSR
C
      DCDPR=CDP1X*S111-DC1*S121
      DCDSR=CDS1X*S111-DC2*S121+DC3*S131
      D11V1=DCDPR+DCDSR
      D10V4=SM2*((X2-1.0)*DCDPR-2.0*DCDSR)
C
      S1=S001-S1M1/R
      S2=S011-S101/R
      S3=S021-S111/R
      DCDPR=CDP1X*S1-DC1*S2
      DCDSR=CDS1X*S1-DC2*S2+DC3*S3
      DD0V1=DCDPR+DCDSR
C
      DCDPR=CDP1X*S011-DC1*S021
      DCDSR=CDS1X*S011-DC2*S021+DC3*S031
      D01V2=-(X2-1.0)*DCDPR+DCDSR
C
      DCDPR=CDP1X*S101-DC1*S111
      DCDSR=CDS1X*S101-DC2*S111+DC3*S121
      D10V2=-(X2-1.0)*DCDPR+DCDSR
      D1MV3=SM2*((X2-3.0)*DCDPR-2.0*DCDSR)
C
      S1=S011-S101/R
      S2=S021-S111/R
      S3=S031-S121/R
      DCDPR=CDP1X*S1-DC1*S2
      DCDSR=CDS1X*S1-DC2*S2+DC3*S3
      DD0V3=SM2*((X2-3.0)*DCDPR-2.0*DCDSR)
C
      QXS(0)=QXS(0)+DD0V1
      VXS(0)=VXS(0)+D1MV1/R
      WXS(0)=WXS(0)+D10V2
C
      D1MV1R=DD0V1
      DD0V1R=(-DD0V1+D1MV1/R)/R-D11V1
      D10V2R=D01V2-D10V2/R
C
      QXS(1)=QXS(1)+DD0V1R
      VXS(1)=VXS(1)+(-D1MV1/R+D1MV1R)/R
      WXS(1)=WXS(1)+D10V2R
C
      D1MV1Z=D1MV3/SM2+D10V2
      DD0V1Z=DD0V3/SM2+D01V2-D10V2/R
      D10V2Z=(D10V4-RM2*D11V1)/(RM2+2.D0*SM2) 
C
      QXS(2)=QXS(2)+DD0V1Z
      VXS(2)=VXS(2)+D1MV1Z/R
      WXS(2)=WXS(2)+D10V2Z
C
C  ** QZ **
C
      DC1=(ODP-SDP)*CDPZ
      DCDPD=CDPZ*S01
      DCDSD=CDSZ*S01-DC1*S02
      D01V1=DCDPD+DCDSD
      D00V4=SM2*((X2-1.0)*DCDPD-2.0*DCDSD)
C
      DCDPD=CDPZ*S10
      DCDSD=CDSZ*S10-DC1*S11
      D10V1=DCDPD+DCDSD
C
      DCDPD=CDPZ*S00
      DCDSD=CDSZ*S00-DC1*S01
      D00V2=-(X2-1.0)*DCDPD+DCDSD
C
      DCDPD=CDPZ*S11
      DCDSD=CDSZ*S11-DC1*S12
      D11V2=-(X2-1.0)*DCDPD+DCDSD
      D10V3=SM2*((X2-3.0)*DCDPD-2.0*DCDSD)
C
      QZS(0)=D10V1
      WZS(0)=D00V2
C
      D10V1R=D01V1-D10V1/R
      D00V2R=-D11V2
C
      QZS(1)=D10V1R
      WZS(1)=D00V2R
C
      D10V1Z=D10V3/SM2+D11V2
      D00V2Z=(D00V4-RM2*D01V1)/(RM2+2.D0*SM2)
C
      QZS(2)=D10V1Z
      WZS(2)=D00V2Z
C
      IF(LS.EQ.NL)GOTO 27
C
      DC1=D2H*CUP2Z
      DC2=D2H*CUS2Z+D2Z*CUP1Z
      DC3=D2Z*DC1
      DCUPR=CUP1Z*S012-DC1*S022
      DCUSR=CUS1Z*S012-DC2*S022+DC3*S032
      D01V1=DCUPR+DCUSR
      D00V4=SM2*((X2-1.0)*DCUPR-2.0*DCUSR)
C
      DCUPR=CUP1Z*S102-DC1*S112
      DCUSR=CUS1Z*S102-DC2*S112+DC3*S122
      D10V1=DCUPR+DCUSR
C
      DCUPR=CUP1Z*S002-DC1*S012
      DCUSR=CUS1Z*S002-DC2*S012+DC3*S022
      D00V2=(X2-1.0)*DCUPR-DCUSR
C
      DCUPR=CUP1Z*S112-DC1*S122
      DCUSR=CUS1Z*S112-DC2*S122+DC3*S132
      D11V2=(X2-1.0)*DCUPR-DCUSR
      D10V3=SM2*(-(X2-3.0)*DCUPR+2.0*DCUSR)
C
      QZS(0)=QZS(0)+D10V1
      WZS(0)=WZS(0)+D00V2
C
      D10V1R=D01V1-D10V1/R
      D00V2R=-D11V2
C
      QZS(1)=QZS(1)+D10V1R
      WZS(1)=WZS(1)+D00V2R
C
      D10V1Z=D10V3/SM2+D11V2
      D00V2Z=(D00V4-RM2*D01V1)/(RM2+2.D0*SM2)
C
      QZS(2)=QZS(2)+D10V1Z
      WZS(2)=WZS(2)+D00V2Z
C
   27 DC1=HD1*CDP2Z
      DC2=HD1*CDS2Z+ZD1*CDP1Z
      DC3=ZD1*DC1
      DCDPR=CDP1Z*S011-DC1*S021
      DCDSR=CDS1Z*S011-DC2*S021+DC3*S031
      D01V1=DCDPR+DCDSR
      D00V4=SM2*((X2-1.0)*DCDPR-2.0*DCDSR)
C
      DCDPR=CDP1Z*S101-DC1*S111
      DCDSR=CDS1Z*S101-DC2*S111+DC3*S121
      D10V1=DCDPR+DCDSR
C
      DCDPR=CDP1Z*S001-DC1*S011
      DCDSR=CDS1Z*S001-DC2*S011+DC3*S021
      D00V2=-(X2-1.0)*DCDPR+DCDSR
C
      DCDPR=CDP1Z*S111-DC1*S121
      DCDSR=CDS1Z*S111-DC2*S121+DC3*S131
      D11V2=-(X2-1.0)*DCDPR+DCDSR
      D10V3=SM2*((X2-3.0)*DCDPR-2.0*DCDSR)
C
      QZS(0)=QZS(0)+D10V1
      WZS(0)=WZS(0)+D00V2
C
      D10V1R=D01V1-D10V1/R
      D00V2R=-D11V2
C
      QZS(1)=QZS(1)+D10V1R
      WZS(1)=WZS(1)+D00V2R
C
      D10V1Z=D10V3/SM2+D11V2
      D00V2Z=(D00V4-RM2*D01V1)/(RM2+2.D0*SM2)
C
      QZS(2)=QZS(2)+D10V1Z
      WZS(2)=WZS(2)+D00V2Z
C
      GOTO 50
C
C *** SDP < (OR =) ODP < DEP2 ***
C
C  ** QX **
C
C   * SH (TRANSMITTED) *
C
   30 D1MH1=H1T*S1M
      DD0H1=H1T*(S00-S1M/R)
      QXS(0)=D1MH1/R
      VXS(0)=DD0H1
C
      D11H1=H1T*S11
      D1MH1R=DD0H1
      DD0H1R=(-DD0H1+D1MH1/R)/R-D11H1
      QXS(1)=(-D1MH1/R+D1MH1R)/R
      VXS(1)=DD0H1R
C
      D1MH2=H2T*S10
      DD0H2=H2T*(S01-S10/R)
      D1MH1Z=D1MH2/SM3
      DD0H1Z=DD0H2/SM3
      QXS(2)=D1MH1Z/R
      VXS(2)=DD0H1Z
C
C   * P-SV (TRANSMITTED) *
C
      D2H=DEP2-SDP
      ZD2=ODP-DEP2
      DC1=D2H*CDP2X
      DC2=D2H*CDS2X+ZD2*CDP1X
      DC3=ZD2*DC1
      DCDPT=CDP1X*S1M-DC1*S10
      DCDST=CDS1X*S1M-DC2*S10+DC3*S11
      D1MV1=DCDPT+DCDST
C
      DCDPT=CDP1X*S11-DC1*S12
      DCDST=CDS1X*S11-DC2*S12+DC3*S13
      D11V1=DCDPT+DCDST
      D10V4=SM3*((X3-1.0)*DCDPT-2.0*DCDST)
C
      S1=S00-S1M/R
      S2=S01-S10/R
      S3=S02-S11/R
      DCDPT=CDP1X*S1-DC1*S2
      DCDST=CDS1X*S1-DC2*S2+DC3*S3
      DD0V1=DCDPT+DCDST
C
      DCDPT=CDP1X*S01-DC1*S02
      DCDST=CDS1X*S01-DC2*S02+DC3*S03
      D01V2=-(X3-1.0)*DCDPT+DCDST
C
      DCDPT=CDP1X*S10-DC1*S11
      DCDST=CDS1X*S10-DC2*S11+DC3*S12
      D10V2=-(X3-1.0)*DCDPT+DCDST
      D1MV3=SM3*((X3-3.0)*DCDPT-2.0*DCDST)
C
      S1=S01-S10/R
      S2=S02-S11/R
      S3=S03-S12/R
      DCDPT=CDP1X*S1-DC1*S2
      DCDST=CDS1X*S1-DC2*S2+DC3*S3
      DD0V3=SM3*((X3-3.0)*DCDPT-2.0*DCDST)
C
      QXS(0)=QXS(0)+DD0V1
      VXS(0)=VXS(0)+D1MV1/R
      WXS(0)=D10V2
C
      D1MV1R=DD0V1
      DD0V1R=(-DD0V1+D1MV1/R)/R-D11V1
      D10V2R=D01V2-D10V2/R
C
      QXS(1)=QXS(1)+DD0V1R
      VXS(1)=VXS(1)+(-D1MV1/R+D1MV1R)/R
      WXS(1)=D10V2R
C
      D1MV1Z=D1MV3/SM3+D10V2
      DD0V1Z=DD0V3/SM3+D01V2-D10V2/R
      D10V2Z=(D10V4-RM3*D11V1)/(RM3+2.D0*SM3) 
C
      QXS(2)=QXS(2)+DD0V1Z
      VXS(2)=VXS(2)+D1MV1Z/R
      WXS(2)=D10V2Z
C
C  ** QZ **
C
      DC1=D2H*CDP2Z
      DC2=D2H*CDS2Z+ZD2*CDP1Z
      DC3=ZD2*DC1
      DCDPT=CDP1Z*S01-DC1*S02
      DCDST=CDS1Z*S01-DC2*S02+DC3*S03
      D01V1=DCDPT+DCDST
      D00V4=SM3*((X3-1.0)*DCDPT-2.0*DCDST)
C
      DCDPT=CDP1Z*S10-DC1*S11
      DCDST=CDS1Z*S10-DC2*S11+DC3*S12
      D10V1=DCDPT+DCDST
C
      DCDPT=CDP1Z*S00-DC1*S01
      DCDST=CDS1Z*S00-DC2*S01+DC3*S02
      D00V2=-(X3-1.0)*DCDPT+DCDST
C
      DCDPT=CDP1Z*S11-DC1*S12
      DCDST=CDS1Z*S11-DC2*S12+DC3*S13
      D11V2=-(X3-1.0)*DCDPT+DCDST
      D10V3=SM3*((X3-3.0)*DCDPT-2.0*DCDST)
C
      QZS(0)=D10V1
      WZS(0)=D00V2
C
      D10V1R=D01V1-D10V1/R
      D00V2R=-D11V2
C
      QZS(1)=D10V1R
      WZS(1)=D00V2R
C
      D10V1Z=D10V3/SM3+D11V2
      D00V2Z=(D00V4-RM3*D01V1)/(RM3+2.D0*SM3)
C
      QZS(2)=D10V1Z
      WZS(2)=D00V2Z
C
   50 RETURN
      END
C
C
C *** SH GENERALIZED R/T MATIX ***
C
      SUBROUTINE GNRTSH(OM,WN)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VP,VS,GTD,GTU,GRD,GRU,TD,TU,RD,RU
      COMMON /SOLY/SDEP,ODEP,LS,LO
     *    /STRC/DNS(20),VP(20),VS(20),THK(20),NL
      COMMON /GRTSH/GTD(0:20),GTU(0:20),GRD(0:20),GRU(0:20)
     *    /MRTSH/TD(0:20),TU(0:20),RD(0:20),RU(0:20)
C
C *** CHECK FOR THE DIMENSION OF GTD, GTU, GRD, AND GRU ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(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 VP,VS,VS2,WS2,CMS,SM,EA(2,2),EB(2,2)
     *      ,TD,TU,RD,RU,DDA,DUB,EAB
      COMMON /SOLY/SDEP,ODEP,LS,LO
     *    /STRC/DNS(20),VP(20),VS(20),THK(20),NL
      COMMON /MRTSH/TD(0:20),TU(0:20),RD(0:20),RU(0:20)
C
C *** CHECK FOR THE DIMENSION OF TD, TU, RD, AND RU ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(10,*)' FROM MDRTSH: NL.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C *** MODIFIED R/T MATRIX ***
C
      WN2=WN*WN
      OM2=OM*OM
C
      VS2=VS(1)*VS(1)
      WS2=OM2/VS2
      CMS=CDSQRT(WN2-WS2)
      IF(DREAL(CMS).LT.0.D0)CMS=-CMS
C
      SM=DNS(1)*VS2
      EA(1,1)=(1.D0,0.D0)
      EA(1,2)=(1.D0,0.D0)
      EA(2,1)=-SM*CMS
      EA(2,2)= SM*CMS
C
      IF(NL.EQ.1)THEN
        DUB=CDEXP(-CMS*SDEP)
        RU(0)=-EA(2,2)*DUB/EA(2,1)
        RETURN
      END IF
C
      IF(LS.EQ.1)THEN
        DUB=CDEXP(-CMS*SDEP)
        DDA=CDEXP(-CMS*(THK(1)-SDEP))
      ELSE
        DUB=CDEXP(-CMS*THK(1))
        DDA=DUB
      END IF
C
      RU(0)=-EA(2,2)*DUB/EA(2,1)
C
      DO 100 I=1, NL-2
C
        VS2=VS(I+1)*VS(I+1)
        WS2=OM2/VS2
        CMS=CDSQRT(WN2-WS2)
        IF(DREAL(CMS).LT.0.D0)CMS=-CMS
C
        SM=DNS(I+1)*VS2
        EB(1,1)=(1.D0,0.D0)
        EB(1,2)=(1.D0,0.D0)
        EB(2,1)=-SM*CMS
        EB(2,2)= SM*CMS
C
        IF(LS.EQ.(I+1))THEN
          DUB=CDEXP(-CMS*SDEP)
        ELSE
          DUB=CDEXP(-CMS*THK(I+1))
        END IF
C
        EAB=EB(2,1)*EA(1,2)-EB(1,1)*EA(2,2) 
        TD(I)=(EA(1,2)*EA(2,1)-EA(1,1)*EA(2,2))*DDA/EAB
        RU(I)=(EB(1,2)*EA(2,2)-EB(2,2)*EA(1,2))*DUB/EAB
        RD(I)=(EB(1,1)*EA(2,1)-EA(1,1)*EB(2,1))*DDA/EAB
        TU(I)=(EB(1,2)*EB(2,1)-EB(2,2)*EB(1,1))*DUB/EAB
C
        EA(1,1)=EB(1,1)
        EA(1,2)=EB(1,2)
        EA(2,1)=EB(2,1)
        EA(2,2)=EB(2,2)
C
        IF(LS.EQ.(I+1))THEN
          DDA=CDEXP(-CMS*(THK(I+1)-SDEP))
        ELSE          
          DDA=DUB
        END IF
C
  100 CONTINUE
C        
      VS2=VS(NL)*VS(NL)
      WS2=OM2/VS2
      CMS=CDSQRT(WN2-WS2)
      IF(DREAL(CMS).LT.0.D0)CMS=-CMS
C
      SM=DNS(NL)*VS2
      EB(1,1)=(1.D0,0.D0)
      EB(1,2)=(1.D0,0.D0)
      EB(2,1)=-SM*CMS
      EB(2,2)= SM*CMS
C
      EAB=EB(2,1)*EA(1,2)-EB(1,1)*EA(2,2) 
C
      IF(LS.EQ.NL)THEN
C
        DUB=CDEXP(-CMS*SDEP)
        TD(NL-1)=(EA(1,2)*EA(2,1)-EA(1,1)*EA(2,2))*DDA/EAB
        RU(NL-1)=(EB(1,2)*EA(2,2)-EB(2,2)*EA(1,2))*DUB/EAB
        RD(NL-1)=(EB(1,1)*EA(2,1)-EA(1,1)*EB(2,1))*DDA/EAB
        TU(NL-1)=(EB(1,2)*EB(2,1)-EB(2,2)*EB(1,1))*DUB/EAB
C
      ELSE
C
        TD(NL-1)=(EA(1,2)*EA(2,1)-EA(1,1)*EA(2,2))*DDA/EAB
        RD(NL-1)=(EB(1,1)*EA(2,1)-EA(1,1)*EB(2,1))*DDA/EAB
C
      END IF
C
      RETURN
      END
C
C
C *** P-SV GENERALIZED R/T MATIX ***
C
      SUBROUTINE GNRTPSV(OM,WN)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VP,VS,GTD,GTU,GRD,GRU,TD,TU,RD,RU,E1(2,2),E2(2,2),DT
      COMMON /SOLY/SDEP,ODEP,LS,LO
     *    /STRC/DNS(20),VP(20),VS(20),THK(20),NL
      COMMON /MRTPSV/TD(2,2,0:20),TU(2,2,0:20),RD(2,2,0:20),RU(2,2,0:20)
     *  /GRTPSV/GTD(2,2,0:20),GTU(2,2,0:20),GRD(2,2,0:20),GRU(2,2,0:20)
C
C *** CHECK FOR THE DIMENSION OF GTD, GTU, GRD, AND GRU ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(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,VP,VS
      COMMON /SOLY/SDEP,ODEP,LS,LO
     *    /STRC/DNS(20),VP(20),VS(20),THK(20),NL
     *    /MRTPSV/TD(2,2,0:20),TU(2,2,0:20),RD(2,2,0:20),RU(2,2,0:20)
C
C *** CHECK FOR THE DIMENSION OF COMMON ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(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 *** P-SV MODIFIED R/T MATRIX FOR Ith LAYER ***
C
        SUBROUTINE MTRPSV(EA,EB,DDA,DUB,L) 
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 EA(4,4),EB(4,4),DDA(2,2),DUB(2,2),TD,TU,RD,RU
     *          ,E1(4,4),E2(4,4),E3(4,4),VP,VS
     * ,a1(4),a2(4),a3(4),b1(4),b2(4),b3(4)  
     * ,c1(4),c2(4),c3(4),d1(4),d2(4),d3(4)  
     * ,aa(4),EE(4,4)
      COMMON /STRC/DNS(20),VP(20),VS(20),THK(20),NL
     *      /MRTPSV/TD(2,2,0:20),TU(2,2,0:20),RD(2,2,0:20),RU(2,2,0:20)
C
C *** CHECK FOR THE DIMENSION OF COMMON ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(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 *** P-SV ELEMENT MATIX ***
C
      SUBROUTINE MTXPSV(OM,WN,L,CMS,CMP,EE)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VP,VS,VS2,WS2,CMS,VP2,WP2,CMP
     *      ,SM,EE(4,4),X0,AA,CMSG,CMPG
      COMMON /STRC/DNS(20),VP(20),VS(20),THK(20),NL
C
C *** CHECK FOR THE DIMENSION OF COMMOM ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(10,*)' FROM MTXPSV: NL.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C *** ELEMENT MATRIX ***
C
      VS2=VS(L)*VS(L)
      VP2=VP(L)*VP(L)
      SM=DNS(L)*VS2
C
      IF(OM.EQ.0.D0)THEN
C  
        X0=(1.D0+VS2/VP2)/(1.D0-VS2/VP2)
        AA=WN*SM
C
        EE(1,1)=(1.D0,0.D0)
        EE(1,2)=(1.D0,0.D0)
        EE(1,3)=(1.D0,0.D0)
        EE(1,4)=(1.D0,0.D0)
        EE(2,1)=-(X0-1.D0)
        EE(2,2)=(1.D0,0.D0)
        EE(2,3)= X0-1.D0
        EE(2,4)=(-1.D0,0.D0)
        EE(3,1)= (X0-3.D0)*AA
        EE(3,2)= -2.D0*AA
        EE(3,3)=-(X0-3.D0)*AA
        EE(3,4)= 2.D0*AA
        EE(4,1)= (X0-1.D0)*AA
        EE(4,2)= -2.D0*AA
        EE(4,3)= (X0-1.D0)*AA
        EE(4,4)=-2.D0*AA
C
      ELSE
C
        WN2=WN*WN
        OM2=OM*OM
C
        WS2=OM2/VS2
        CMS=CDSQRT(WN2-WS2)
        IF(DREAL(CMS).LT.0.D0)CMS=-CMS
        CMSG=CMS/WN
C
        WP2=OM2/VP2
        CMP=CDSQRT(WN2-WP2)
        IF(DREAL(CMP).LT.0.D0)CMP=-CMP
        CMPG=CMP/WN
C
        X0=2.D0*WN-WS2/WN
C
        EE(1,1)=(-1.D0,0.D0)
        EE(1,2)= CMSG
        EE(1,3)=(-1.D0,0.D0)
        EE(1,4)= CMSG
        EE(2,1)=-CMPG
        EE(2,2)=(1.D0,0.D0)
        EE(2,3)= CMPG
        EE(2,4)=(-1.D0,0.D0)
        EE(3,1)= 2.D0*SM*CMP 
        EE(3,2)=-SM*X0 
        EE(3,3)=-2.D0*SM*CMP 
        EE(3,4)= SM*X0 
        EE(4,1)= SM*X0 
        EE(4,2)=-2.D0*SM*CMS 
        EE(4,3)= SM*X0 
        EE(4,4)=-2.D0*SM*CMS 
C
      END IF
C
      RETURN
      END
C
C
C *** P-SV UP/DOWN-GOING MATRIX ***
C
      SUBROUTINE DDUPSV(OM,WN,THK,CMS,CMP,DU)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 DU(2,2),CMS,CMP
C
      IF(OM.EQ.0.D0)THEN
        EW=DEXP(-WN*THK)
        DU(1,1)= EW
        DU(1,2)= (0.D0,0.D0)
        DU(2,1)= -WN*THK*EW
        DU(2,2)= EW
      ELSE
        DU(1,1)= CDEXP(-CMP*THK)
        DU(1,2)= (0.D0,0.D0)
        DU(2,1)= (0.D0,0.D0)
        DU(2,2)= CDEXP(-CMS*THK)
      END IF
C
      RETURN
      END
C
C
C   NAME                - BESJ0 : DOUBLE PRECISION
C
      FUNCTION BESJ0 (X)
C
      IMPLICIT REAL*8 (A-H,O-Z)
      DATA C0,C1,C2,C3,C4,C5,C6,C7,C8,C9,C10
     *   / -2.2270 614            D-9,  7.8594 0046           D-8,
     *     -1.9891 96578 6        D-6,  4.0310 96038 59       D-5,
     *     -6.4499 78225 950      D-4,  7.9012 34182 8039     D-3,
     *     -7.1111 11100 46304    D-2,  4.4444 44444 25769 5  D-1,
     *     -1.7777 77777 77584 46 D 0,  3.9999 99999 99989 66 D 0,
     *     -3.9999 99999 99999 78 D 0                             /
      DATA D0,D1,D2,D3,D4,D5,D6,D7,D8,D9,D10,D11,D12,D13
     *   /  3.3150 18             D-11,  4.6970 63             D-11,
     *     -5.6633 6210           D- 9, -4.8842 9867           D- 9,
     *      6.7810 96622 5        D- 7,  2.1332 90985 8        D- 7,
     *     -5.3795 00918 201      D- 5,  1.1202 21502 384      D- 5,
     *      2.5202 11970 18629    D- 3, -1.7073 87596 90737    D- 3,
     *     -5.6148 69347 44975 4  D- 2,  5.6040 47189 80261 6  D- 2,
     *      3.2757 91375 91465 09 D- 1, -1.7759 67713 14338 34 D- 1 /
      DATA E0,E1,E2,E3,E4,E5,E6,E7,E8,E9,E10,E11,E12,E13,E14
     *   / -1.9865 0              D-12, -1.5178 02          D-11,
     *      3.9649 790            D-10,  2.4051 2241        D- 9,
     *     -5.7071 43476          D- 8, -2.5797 89310 5     D- 7,
     *      5.6607 46146 93       D- 6,  1.7226 12595 802   D- 5,
     *     -3.5284 91002 3231     D- 4, -5.9314 89738 8202  D- 4,
     *      1.1790 12935 71024 6  D- 2,  6.3961 29897 84324 D- 3,
     *     -1.5037 41226 51373 90 D- 1,  4.6828 23482 34590 D- 3,
     *      3.0007 92705 19555 60 D- 1                            /
      DATA A0,A1,A2,A3,A4,A5,A6,A7,A8,A9,A10
     *   /  2.8228 605             D-10, -1.9415 5134       D-9,
     *      6.5160 2613            D- 9, -1.5339 82923      D-8,
     *      3.2676 79187           D- 8, -8.0992 84883      D-8,
     *      2.8872 15007 2         D- 7, -1.7424 98968 02   D-6,
     *      2.1846 78334 973       D- 5, -8.7658 21590 6110 D-4,
     *      0.7978 84560 80286 526 D  0                         /
      DATA B0,B1,B2,B3,B4,B5,B6,B7,B8,B9,B10,B11,B12
     *   / -8.9331 052           D-10,  6.6532 7442       D-9,
     *     -2.3224 79945         D- 8,  5.1878 39807      D-8,
     *     -8.7651 63373         D- 8,  1.3089 75220 0    D-7,
     *     -2.0632 29204 9       D- 7,  4.0856 34631 8    D-7,
     *     -1.1593 30628 42      D- 6,  5.2586 51475 08   D-6,
     *     -4.4239 73781 374     D- 5,  9.1310 64157 0497 D-4,
     *     -9.9735 57010 03581 5 D- 2                         /
      DATA PI4  / 0.78539 81633 97448 31D0 /
      AX=ABS(X)
      IF (AX .GE. 1.0D12) GO TO 40
      IF (AX .GE. 8.0D0) GO TO 30
      IF (AX .GE. 6.0D0) GO TO 20
      IF (AX .GE. 4.0D0) GO TO 10
C
C     !X! < 4
C
      Z=(X*0.25D0)**2
      BESJ0=((((((((((C0*Z +C1)*Z +C2)*Z +C3)*Z +C4)*Z +C5)*Z +C6)*Z
     *       +C7)*Z +C8)*Z +C9)*Z +C10)*Z +1.0D0
      GO TO 9999
C
C     4 <= !X! < 6
C
   10 Z=AX-5.0D0
      BESJ0=((((((((((((D0*Z +D1)*Z +D2)*Z +D3)*Z +D4)*Z +D5)*Z
     *       +D6)*Z +D7)*Z +D8)*Z +D9)*Z +D10)*Z +D11)*Z +D12)*Z +D13
      GO TO 9999
C
C     6 <= !X! < 8
C
   20 Z=AX-7.0D0
      BESJ0=(((((((((((((E0*Z +E1)*Z +E2)*Z +E3)*Z +E4)*Z +E5)*Z
     *       +E6)*Z +E7)*Z +E8)*Z +E9)*Z +E10)*Z +E11)*Z +E12)*Z
     *       +E13)*Z +E14
      GO TO 9999
C
C     8 <= !X! < 1.0D12
C
   30 Z=(8.0D0/AX)**2
      A=(((((((((A0*Z +A1)*Z +A2)*Z +A3)*Z +A4)*Z +A5)*Z +A6)*Z
     *  +A7)*Z +A8)*Z +A9)*Z +A10
      B=((((((((((((B0*Z +B1)*Z +B2)*Z +B3)*Z +B4)*Z +B5)*Z +B6)*Z
     *  +B7)*Z +B8)*Z +B9)*Z +B10)*Z +B11)*Z +B12) /AX
      Q=AX-PI4
      BESJ0=(A*COS(Q) - B*SIN(Q)) /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
