c
c *** This code computes high-frequency strong groun motions from a fault model
c     with sub-faults. The main features are the use of
c     1) Boore's point sources (1983) on sub-faults,
c     2) Hisada's hybrid phase spectra (zero at low freq, and random at high freq.)
c     2) Irikura's semi-empirical method (1986),
c     3) Green's function of layered half-spaces
c     (August, 2006)
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 *** Before running this code, you must run "phs3x.f" to get the locations of
c     poles corresponding to the normal mode solutions; the data stored in
c     "phsl1_dat" and "phsr1_dat" for the phase velocities of the Love and
c     and Rayleigh waves, respectively.
c
c *** After running this code, you can plot waveforms using a FFT code "grfftps.f".
c
c * Parameters :
c     MOM : Max Number of Omegas (ex., 18212)
c     MOB : Max NUmber of Observation Points
c     MFL : Max Number of Sub-Faults
c     NL  : Max Number of Layers used in common statements (ex., 20)
c     NWN : Max Number of Wavenumbers used in common statements (ex., 29998)
c     NW  : Max Number of Poles and Branch Points used in common statements (ex., 801)
c
      PROGRAM grflto2x
C
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(MOM=18212,MOB=51,MFL=503)
      COMPLEX*16 GX(MOM,MOB),GY(MOM,MOB),GZ(MOM,MOB),Frate(MOM,MFL,2)
     *          ,CPX,EDT,SS,SX,SY
      DIMENSION rate(mom,MFL,2)
      COMMON /FREQ/DOM,IOM,NOM,OMI
     *    /STRC/DNS(20),VP(20),QP(20),FP(20)
     *                 ,VS(20),QS(20),FS(20),THK(20),NL
     *    /FLT/STR,DIP,RAK(503),FL,FW,FD(503),VR,TLc,Pn(503)
     *        ,XF(503),YF(503),TD(503),DP(503),fc(503),fmax(503),SM0L
     *    /IFLT/NF,NSTR,NDIP,NFfunc,ND,NSF,NDF
C
      WRITE(*,*)
      WRITE(*,*)
      WRITE(*,*)'"grflto2x.f" computes broad-band strong motion'
      WRITE(*,*)'   using the Boore source model and the Green function'
      write(*,*)'   of layered half-space'
      write(*,*)' (Originally made in 1992, and modified in 2010)'
      write(*,*)
      WRITE(*,*)'"grflto2x-in.csv" is the input data'
      WRITE(*,*)
     *' You must run "phs3x.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.'
C
      CALL INDAT(MOM,MOB,MFL,Frate,rate)
c
c   * Output of Soure Spectra and Source Time Functiona
c
      CPX=(0.D0,1.D0)
      pai=4.d0*datan(1.d0)
      sNF=Dfloat(ND)
      FDsum=0.d0
      Fmaxsum=0.d0
      do j=1,NF
        FDsum=FDsum+FD(j)
        Fmaxsum=Fmaxsum+Fmax(j)
      end do
      FDavg=FDsum/dfloat(NF)
      Fmaxavg=Fmaxsum/dfloat(NF)
      SLENGTH=FL*dfloat(NSTR)
      DWIDTH=FW*dfloat(NDIP)
c
      VSsum=0.d0
      do 100 id=1,NDIP
        ZS=DP(id)
        DRU=0.d0
        DRL=THK(1)
        DO 120 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)
  120   CONTINUE
        VSsum=VSsum+VS(LS)
  100 continue
      VSavg=VSsum/dfloat(NDIP)
c
      RL=dsqrt(SLENGTH*DWIDTH/pai)/1000.d0
      if(NF.eq.1)then
        fcL=fc(1)
        sdropL=SM0L*(fcL/(4.9d+6*VSavg/1000.d0))**3
      else
        sdropL=7.d0*SM0L/(16.d0*RL*RL*RL*(10.d0**21))
        sfac=1.d0/3.d0
        fcL=(4.9d+6)*(sdropL/SM0L)**sfac*VSavg/1000.d0
      end if
      write(*,*)
      write(*,*)'VSavg(m/s),FDavg(m)=',VSavg,FDavg
      write(*,*)'Fmaxavg(Hz),SM0L(dyn-cm)=',Fmaxavg,SM0L
      write(*,*)'RL(km),Fmax(No.1 Hz)=',RL,Fmax(1)
      write(*,*)'sdropL(bar),fcL(Hz)=',sdropL,fcL
      write(10,*)
      write(10,*)'VSavg(m/s),FDavg(m)=',VSavg,FDavg
      write(10,*)'Fmaxavg(Hz),SM0L(dyn-cm)=',Fmaxavg,SM0L
      write(10,*)'RL(km),Fmax(No.1 Hz)=',RL,Fmax(1)
      write(10,*)'sdropL(bar),fcL(Hz)=',sdropL,fcL
c
      open(21,file='source-spectra-acc.csv')
      write(21,*)'freq.(Hz),moment-theory,moment-X,moment-Y,F-func'
      DO 200 IO=IOM,NOM
        OM=DOM*DFLOAT(IO)
        freq=om/(2.d0*pai)
     	FFc=freq/fcL
        SS=SM0L/(1.d0+FFc*FFc)
        FFmax=freq/Fmaxavg
        SS=SS*OM*OM/dsqrt(1.d0+FFmax**(2.d0*Pn(1)))
        SX=(0.d0,0.d0)
        SY=(0.d0,0.d0)
        OMTd=OM*TLc
        sinc=dsin(OMTd)/OMTd
c
        IS=0
        do 220 id=1,NDIP
          ZS=DP(id)
          DRU=0.d0
          DRL=THK(1)
          DO 230 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)
  230     CONTINUE
          SM=DNS(LS)*VS(LS)*VS(LS)
c
          do 240 j=1,nstr
            IS=IS+1
            if(FD(IS).eq.0.d0)goto 200
            EDT=CDEXP(cpx*OM*TD(IS))
            SM0=SM*FD(IS)*FL*FW/sNF*1.0D+10
c    * Irikura's F-function *
            if(NFfunc.eq.0)then
              EDT=EDT*(1.d0+(sNF-1.d0)*sinc*cdexp(cpx*OMTd))
c    * Exp-type F-function *
            else
              EDT=EDT*(1.d0+(sNF-1.d0)/(1.d0-cpx*OMTd))
            end if
            SX=SX+SM0*EDT*Frate(io+1,is,1)
            SY=SY+SM0*EDT*Frate(io+1,is,2)
  240     continue
  220   continue
        write(21,2000)freq,cdabs(SS),cdabs(SX),cdabs(SY),cdabs(EDT)
  200 continue
 2000 format(f8.4,4(',',e12.4))
      close(21)
c
      CALL PHSBRC(MOM)
      CALL HKQVW(MOM,MOB,MFL,Frate,GX,GY,GZ)
      CALL OTDAT(MOM,MOB,GX,GY,GZ)
      WRITE(*,*)
      WRITE(*,*)'Use FFT program "grfftpx.f" to plot waveforms'
C
      STOP
      END
C
c
C *** DATA INPUT ***
C
      SUBROUTINE INDAT(MOM,MOB,MFL,Frate,rate)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 Frate(MOM,MFL,2),cpx,Ffunc
      DIMENSION rate(mom,MFL,2),sdrop(503)
      character infile*20
      COMMON /FREQ/DOM,IOM,NOM,OMI /OBS/XO(51),YO(51),ZO(51),NOB
     *    /time/dt,period /ntime/nt
     *    /STRC/DNS(20),VP(20),QP(20),FP(20)
     *                 ,VS(20),QS(20),FS(20),THK(20),NL
     *    /INTG/VPMAX,VSMAX,RYMIN,CFINAL,N11,N12,M11
     *    /FLT/STR,DIP,RAK(503),FL,FW,FD(503),VR,TLc,Pn(503)
     *        ,XF(503),YF(503),TD(503),DP(503),fc(503),fmax(503),SM0L
     *    /Hypo/xhypo,yhypo,zhypo /Radiation/fR1,fR2,fran
     *    /IFLT/NF,NSTR,NDIP,NFfunc,ND,NSF,NDF /SIGN/ISIGN
     *    /OTST/IOT,IST /outom/ioutom /Nboore/Niter,idum
c
c *** Check of Dimensions of Common Paramters ***
C
      NML =20
c
      nc=51
      if(nc.lt.MOB)then
        write(10,*)'from Indat:'
        WRITE(10,*)' NC < MOB !!!!!!!!!!!!!'
        write(*,*)'from Indat:'
        WRITE(*,*)' NC < MOB !!!!!!!!!!!!!'
        STOP
      END IF
c
      nc=503
      if(nc.lt.MFL)then
        write(10,*)'from Indat:'
        WRITE(10,*)' NC < MFL !!!!!!!!!!!!!'
        write(*,*)'from Indat:'
        WRITE(*,*)' NC < MFL !!!!!!!!!!!!!'
        STOP
      END IF
C
C *** OMEGA DATA ***
C
      write(*,*)
      write(*,*)'Enter File name for "grflto2x-in.csv"'
      read(*,'(a20)')infile
      OPEN (UNIT=7,FILE=infile,STATUS='OLD')
c      OPEN (UNIT=7,FILE='grflto2x-test1.csv',STATUS='OLD')
C
      pi2=datan(1.d0)*8.d0
      read(7,*)
      read(7,*)
      read(7,*)dt,nt
      write(*,*)
      write(*,*)'delta time (s), and number of times =',dt,nt
      read(7,*)
      read(7,*)period,OMI
      write(*,*)'minimum period (s) =',period
      duration=dt*nt
      dom=pi2/duration
      iom=1
      nom=int(duration/period)
c
      WRITE(10,*)'FROM indat:'
      WRITE(10,*)
      WRITE(10,1130)nt,dt,period
      WRITE(10,1140)NOM,DOM,OMI
 1130 format(/' Num. of Time, Delta Time and Min. Period (s)=',i5,2f8.4)
 1140 format( ' Num. of Omega, dOM, and Imag OM (rad/sec) =',i5,2f10.7)
C
      IF(NOM.GT.MOM)THEN
        write(10,*)'from Indat:'
        WRITE(10,*)' NOM>MOM !!!!!!!!!!!!!'
        write(*,*)'from Indat:'
        WRITE(*,*)' NOM>MOM !!!!!!!!!!!!!'
        stop
      END IF
C
C *** MEDIUM DATA ***
C
      read(7,*)
      read(7,*)
      READ(7,*)NL
      WRITE(10,2200)NL
      WRITE(*,2200)NL
      if(NML.lt.NL)then
        write(10,*)'from Indat:'
        WRITE(10,*)' NC < NL !!!!!!!!!!!!!'
        write(*,*)'from Indat:'
        WRITE(*,*)' NC < NL !!!!!!!!!!!!!'
        STOP
      END IF
      read(7,*)
      DO 100 L=1,NL
        READ(7,*)i,DNS(L),VP(L),QP(L),FP(L),VS(L),QS(L),FS(L),THK(L)
      WRITE(*,2210)L,DNS(L),VP(L),QP(L),FP(L),VS(L),QS(L),FS(L),THK(L)
      WRITE(10,2210)L,DNS(L),VP(L),QP(L),FP(L),VS(L),QS(L),FS(L),THK(L)
  100 CONTINUE
 2200 FORMAT(//' Total Number of Layers =',I5/,' Layer,Density,Vp(real),
     *    Qp  ,  Fp  ,Vs(real),    Qs  ,  Fs  ,Thickness(m)')
 2210 FORMAT(I5,F9.3,2(F9.2,F9.3,F7.3),F12.3)
C
C *** MAX Vs and Vp and MIN Rayleigh Phase Velocity ***
C
      VSMAX=VS(NL)
      VPMAX=VP(NL)
      VSMIN=VS(1)
      MINVS=1
      DO 120 IL=1,NL
        IF(VS(IL).GT.VSMAX)VSMAX=VS(IL)
        IF(VP(IL).GT.VPMAX)VPMAX=VP(IL)
        IF(VS(IL).LT.VSMIN)THEN
          VSMIN=VS(IL)
          MINVS=IL
        END IF
  120 CONTINUE
      RYMIN=FRHF(VP(MINVS),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 *** FAULT DATA ***
C
      read(7,*)
      read(7,*)
      READ(7,*)SLENGTH,DWIDTH,NSTR,NDIP,ND
      write(*,*)
      write(*,*)'Fault Model Parameters:'
      write(*,*)'  Fault L and W (m) =',SLENGTH,DWIDTH
      write(*,*)'  Numbers of sub-fault along L and W =',NSTR,NDIP
      write(*,*)'  ND (=Large EQ Slip/Small EQ Slip)',ND
      NF=NSTR*NDIP
c      sNF=dsqrt(dfloat(NF))
      sNF=dfloat(ND)
      IF(NF.GT.MFL)THEN
        write(10,*)'from Indat:'
        WRITE(10,*)' NF.GT.MFL !!!!!!'
        write(*,*)'from Indat:'
        WRITE(*,*)' NF.GT.MFL !!!!!!'
        STOP
      END IF
c
      FL=slength/nstr
      FW=dwidth/ndip
      read(7,*)
      READ(7,*)DELAY,STR,DIP,VR,dtime,idum
      write(*,*)'  Starting Time (s) and Vr (m/s)=',DELAY,Vr
      write(*,*)'  Strike and Dip Angles (deg) =',STR,DIP
      write(*,*)'  Average Random Delay Time (s) =',dtime
      str=str*pi2/360.d0
      dip=dip*pi2 /360.d0
      read(7,*)
      READ(7,*)XHYPO,YHYPO,ZHYPO
      write(*,*)'  Hypocenter Location (X,Y,Z: m)   =',XHYPO,YHYPO,ZHYPO
      read(7,*)
      READ(7,*)xflt,yflt,zflt
      write(*,*)'  Fault Origin Location (X,Y,Z: m) =',xflt,yflt,zflt
      read(7,*)
c
      dfz=fw*dsin(dip)
      FDsum=0.d0
      SM0=0.d0
      j=1
      do 150 id=1,NDIP
        ZS=zflt-dfz*(dfloat(id-1)+0.5d0)
        DRU=0.d0
        DRL=THK(1)
        DO 160 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)
  160   CONTINUE
        SM=DNS(LS)*VS(LS)*VS(LS)
c        write(*,*)'id,ZS,LS,SM=',id,ZS,LS,SM
        do 150 is=1,NSTR
          READ(7,*)i,FD(j),rak(j),sdrop(j),fc(j),fmax(j),Pn(j)
          rak(j)=rak(j)*pi2/360.d0
          FDsum=FDsum+FD(j)
          sSM0=SM*FD(j)*FL*FW*1.0D+10
c      write(*,*)'i,FD,rak,sdrop,fmax=',i,FD(j),rak(j),sdrop(j),fmax(j)
          SM0=SM0+sSM0
          j=j+1
  150 continue
      SM0L=SM0
      FDavg=FDsum/dfloat(NF)*100.d0
c      SM0avg=SM0L/dfloat(NF)
      write(*,*)
      write(*,*)'  Average Slip (cm)     =',FDavg
      write(*,*)'  Seismic Moment(dyne-cm)=',SM0L
      sMw=(dlog10(SM0L)-16.1d0)/1.5d0
      sMj=(dlog10(SM0L)-17.72d0)/1.17d0
      write(*,*)'  Mw (by Hanks and Kanamori,1979) =',sMw
      write(*,*)'  Mj (by M.Takemura, 1990)        =',sMj
c
      write(*,*)
      sfac=1.d0/3.d0
c      sSM0=SM0avg/sNF
      j=1
      do 170 id=1,NDIP
        if(sdrop(j).eq.0.d0)goto 170
        ZS=zflt-dfz*(dfloat(id-1)+0.5d0)
        DRU=0.d0
        DRL=THK(1)
        DO 180 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)
  180   CONTINUE
        SM=DNS(LS)*VS(LS)*VS(LS)
        do 185 is=1,NSTR
          sSM0=SM*(FD(j)/sNF)*FL*FW*1.0D+10
          fc(j)=(4.9d+6)*(sdrop(j)/sSM0)**sfac*Vs(LS)/1000.d0
          j=j+1
  185   continue
  170 continue
c
c *** Sub-fault parameters ***
c
c   * Dip direction *
      dcs=fw*dcos(dip)
      cst=dcos(str)
      sst=dsin(str)
c
c   * Introduction of Random Rupture Time Delay *
      iflt=0
      ihyp=1
      dismin=999999999999999.d0
      do 200 i=1,ndip
        dc=dcs*(dfloat(i-1)+0.5d0)
        zc=zflt-dc*dtan(dip)
        if(zc.lt.0.d0)then
          write(*,*)'Fault Depth is less than 0!: Dip Number=',i
c          stop
        end if
c     * Strike direction *
        do 220 j=1,nstr
          iflt=iflt+1
          sc=fL*(dfloat(j-1)+0.5d0)
          dp(iflt)=zc
          zz=zc-zhypo
          xf(iflt)=xflt+sc*cst+dc*sst
          yf(iflt)=yflt+sc*sst-dc*cst
	  xx=xf(iflt)-xhypo
	  yy=yf(iflt)-yhypo
          dis=dsqrt(xx*xx+yy*yy+zz*zz)
          if(dis.lt.dismin)then
            ihyp=iflt
            dismin=dis
          end if
          td(iflt)=delay+dis/VR
  220   continue
  200 continue
c
      do iflt=1,NF
        if(iflt.ne.ihyp)then
          td(iflt)=td(iflt)+dtime*dabs(gasdev(idum))
        end if
      end do
C
      WRITE(10,2300)NF
      WRITE(10,2305)STR,DIP
      WRITE(10,2315)slength,dwidth
      WRITE(10,2320)FL,FW
      WRITE(10,2330)FDavg,SM0L
      WRITE(10,2335)sMw,sMj
      write(10,2340)ihyp,dismin
 2300 FORMAT(//'Seismic Source Data:'
     *      ,/'  Total Number of Sources =',I3)
 2305 FORMAT( '  Strike Angle (rd)=',F12.8,',  Dip Angle   (rd)=',F12.8)
 2315 FORMAT( '  Fault Length  (m)=',F12.3,',  Width(m)        =',F12.3)
 2320 FORMAT( '  Sub-Length    (m)=',F12.3,',  Sub-Width(m)    =',F12.3)
 2330 FORMAT( '  Avgerage Slip (m)=',e12.4,
     *        ',  Seismic Moment (dyne-cm)=',e12.4)
 2335 FORMAT( '  Moment Magnitude =',F12.1,',  JMA Magnitude   =',F12.1)
 2340 FORMAT( '  Fault No. of Hyp =',i12,
     *        ',  Distance(Center to Hyp: m) =',F12.4)
c
      II=0
      DO 270 I=1,NDIP
        dc=dcs*(dfloat(i-1)+0.5d0)
        zc=zflt-dc*dtan(dip)
        DO 270 J=1,NSTR
          II=II+1
          sc=fL*(dfloat(j-1)+0.5d0)
          zz=zc-zhypo
          xf(ii)=xflt+sc*cst+dc*sst
          yf(ii)=yflt+sc*sst-dc*cst
	  xx=xf(ii)-xhypo
	  yy=yf(ii)-yhypo
          dis=dsqrt(xx*xx+yy*yy+zz*zz)
        WRITE(10,2350)II,I,J
        WRITE(10,2360)FD(ii),rak(ii)
        write(10,2370)XF(II),YF(II),DP(II)
        write(10,2375)sdrop(ii),fc(ii),fmax(ii)
        write(10,2380)TD(II),dis/VR,TD(II)-dis/VR
  270 CONTINUE
 2350 FORMAT(/'   Source Number   =',I4,' (',I2,',',I2,')')
 2360 FORMAT( '   Dislocation (m)=',e12.4,', Rake Ang. (rad)=',F9.5)
 2370 FORMAT( '   Source X, Y, and Z Location (m) =',3f9.1)
 2375 format( '   Stress Drop (bar), fc & fmax (Hz) =',f8.2,f9.3,f12.3)
 2380 FORMAT( '   Total Time Delay, r/Vr & delta Tr (s) =',3f9.3)
c
c *** Input Boore's Source Parameters ***
c
      call Boore(MOM,MFL,Frate,rate)
C
C *** Radiation Pattern for S and P Waves at high Frequencies ***
C
      read(7,*)
      read(7,*)
      READ(7,*)fR1,fR2
      write(*,2410)fR1,fR2
      write(10,2410)fR1,fR2
 2410 FORMAT(/' Transient Corners Frequencies (Hz) = ',2f12.3)
C
C *** Selection of F-function, which adjusts the difference of the slip
c     functions between the whole- and sub-faults ***
C
      read(7,*)
      read(7,*)
      READ(7,*)NFfunc,TLc
c      write(*,*)NFfunc,TLc
      WRITE(10,2430)NFfunc,TLc
 2430 FORMAT(/' Choice of F-function and TLc (s)',i5,f8.2)
C
C *** Selection of High-freq. Source Types ***
C
      read(7,*)
      read(7,*)
      READ(7,*)NSF
c      write(*,*)NSF
      write(10,2450)
      IF(NSF.EQ.0)WRITE(10,2460)
      IF(NSF.EQ.1)WRITE(10,2470)
      IF(NSF.EQ.2)WRITE(10,2480)
      IF(NSF.EQ.3)WRITE(10,2490)
 2450 FORMAT(/' Choice of High-Frequency Source Types:')
 2460 FORMAT(' Full-Wave Source (P, SH, and SV weves)')
 2470 FORMAT(' P-Wave Source')
 2480 FORMAT(' SH-Wave Source')
 2490 FORMAT(' SV-Wave Source')
C
C *** Selection of High-Freq. Distribution of Horizontal Components ***
C
      read(7,*)
      read(7,*)
      READ(7,*)NDF
c      write(*,*)NDF
      write(10,2500)
      IF(NDF.EQ.0)WRITE(10,2510)
      IF(NDF.EQ.1)WRITE(10,2520)
 2500 FORMAT(/' Distribution of High-Frequency Horizontal Amplitudes:')
 2510 FORMAT(' Equal Distribution for X and Y components')
 2520 FORMAT(' Original Distribution for P, SH, and SV Waves')
c
c *** Output of F-function ***
c
      write(*,*)'  F-function (No.1) is stored in F-func.csv'
      open(21,file='F-func(No.1).csv')
      write(21,*)'freq.(Hz),',' F-func'
      CPX=(0.D0,1.D0)
      pai=4.d0*datan(1.d0)
      do io=iom,nom
        OM=DOM*DFLOAT(IO)
        freq=om/(2.d0*pai)
        OMTd=OM*TLc
c   *  Irikura's F-function (1987) *
        if(NFfunc.eq.0)then
          sinc=dsin(OMTd)/OMTd
          Ffunc=1.d0+(sNF-1.d0)*sinc*exp(cpx*OMTd)
c   *  Exp-type Slip Velocity Function *
        else
          Ffunc=1.d0+(sNF-1.d0)/(1.d0-cpx*OMTd)
        end if
        write(21,2550)freq,cdabs(Ffunc)
      end do
      close(21)
 2550 format(f12.4,',',e12.4)
C
C *** NUMERICAL INTEGRAL DATA ***
C
      read(7,*)
      read(7,*)
      READ(7,*)N11
      read(7,*)
      READ(7,*)N12
      read(7,*)
      READ(7,*)cfac
c      write(*,*)N11,N12,cfac
      CFINAL=rymin/cfac
C
      WRITE(10,2600)
      WRITE(10,2610)N11
      WRITE(10,2620)N12
      WRITE(10,2630)CFINAL
      WRITE(10,2640)2*(N11+N12)+1
 2600 FORMAT(//' Number of Points for Numerical Integrations:')
 2610 FORMAT( /'   From 0              to 1.2*om/Ray(min) =',I5)
 2620 FORMAT(  '   From 1.2*om/Ry(min) to om/C(final)     =',I5)
 2630 FORMAT(  '   C(final = Minimun Rayleigh Phase Vel.) =',F10.3)
 2640 FORMAT(  '   Initial Number of Points: 2*(N1+N2)+1  =',I5)
C
C *** DATA FOR THE ASYMPTOTIC ***
C
      read(7,*)
      read(7,*)
      READ(7,*)IST
c      write(*,*)IST
      IF(IST.EQ.1)THEN
        WRITE(10,2650)
      ELSE
        WRITE(10,2660)
      END IF
 2650 FORMAT(//' THE ASYMPTOTIC METHOD IS USED')
 2660 FORMAT(//' THE ASYMPTOTIC METHOD IS NOT USED')
C
C *** CHANGE OF SIGNS FOR IMAGINARY PARTS ***
C
      read(7,*)
      read(7,*)
      READ(7,*)ISIGN
c      write(*,*)ISIGN
      IF(ISIGN.EQ.1)WRITE(10,2700)
      IF(ISIGN.EQ.0)WRITE(10,2710)
 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,*)
      read(7,*)
      READ(7,*)NOB
c      write(*,*)NOB
      WRITE(10,2800)NOB
 2800 FORMAT(//' TOTAL NUMBER OF OBSERVATION POINTS =',I6/)
C
      IF(NOB.GT.MOB)THEN
        WRITE(10,*)' NOB>NMOB !!!!!!!!!!!!!'
        STOP
      END IF
C
      read(7,*)
      DO 300 I=1,NOB
        READ(7,*)j,XO(I),YO(I),ZO(I)
c        write(*,*)j,XO(I),YO(I),ZO(I)
        WRITE(10,2820)I,XO(I),YO(I),ZO(I)
  300 CONTINUE
 2820 FORMAT('   OB NUM.=',I5,',     (X,Y,Z) =',3F12.4)
C
      CLOSE(UNIT=7)
C
c *** output of fault and observation data ***
c
c      open(9,file='fault_dat',status='unknown')
c        write(9,9500)XHYPO,YHYPO,ZHYPO
c        write(9,9500)xflt,yflt,zflt
c        write(9,9500)STR,DIP,RAK(1)
c        write(9,9500)slength,dwidth,VR
c        write(9,9510)NSTR,NDIP
c 9500   format(3f12.4)
c 9510   format(3i12)
c      close(9)
c
c      open(9,file='obs_dat',status='unknown')
c        do i=1,nob
c        write(9,9500)xo(i),yo(i),zo(i)
c        end do
c      close(9)
C
C *** OUTPUT OF INTEGRAND OF U,V,AND W ***
C
      IOT=0
      IF(NOB.NE.1)RETURN
      IF(NF.NE.1)RETURN
c      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
        f1=dom/6.283185307d0
        f2=f1*nom
        write(*,*)
        write(*,*)'Frequency Range (Hz) =',f1,f2
        write(10,*)'Frequency Range (Hz) =',f1,f2
        write(*,*)'Enter Output Frequency (Hz)'
        read(*,*)ofreq
        df=9999.d0
        do i=1,nom
          freq=f1*i
          dfreq=dabs(freq-ofreq)
          if(dfreq.lt.df)then
            ioutom=i
            df=dfreq
          end if
        end do
        write(*,*)'Actual Output Frequency (Hz)=',f1*ioutom
        write(10,*)'Actual Output Frequency (Hz)=',f1*ioutom
        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)vp(i)
            write(9,9000)vs(i)
 500      continue
        close(9)
 9000   format(e15.7)
 9100   format(i5)
c
      end if

      RETURN
      END
C
c
c
c *** Randam Number with a normally distributed deviate with zero mean
c     and unit variance, using ran0 as the source of uniform deviates
c
      function gasdev(idum)
      IMPLICIT REAL*8(A-H,O-Z)
      save iset,gset
      data iset/0/

      if(iset.eq.0)then
   10   v1=2.d0*ran0(idum)-1.d0
        v2=2.d0*ran0(idum)-1.d0
        rsq=v1*v1+v2*v2
        if(rsq.ge.1.d0.or.rsq.eq.0.d0)goto 10
        fac=dsqrt(-2.d0*dlog(rsq)/rsq)
        gset=v1*fac
        gasdev=v2*fac
        iset=1
      else
        gasdev=gset
        iset=0
      end if

      return
      end

c *** Randam Function **
c
      function ran0(idum)
      IMPLICIT REAL*8(A-H,O-Z)
      parameter(IA=16807,IM=2147483647,AM=1.d0/IM,IQ=127773
     *         ,IR=2836,MASK=123459876)
c
      idum=ieor(idum,MASK)
      k=idum/IQ
      idum=IA*(idum-k*IQ)-IR*k
      if(idum.lt.0)idum=idum+IM
      ran0=AM*idum
      idum=ieor(idum,MASK)
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 *** Calculation of Normalized Boore's Moment Rate Functions
c     (Far Field Displacement) at the all sub-faults ***
c
      SUBROUTINE Boore(MOM,MFL,Frate,rate)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 FX(MOM),FX1(MOM),cran,Frate(MOM,MFL,2)
      DIMENSION XR(MOM),bov(mom),bovN(mom),rate(mom,MFL,2),bos(mom)
      COMMON /FREQ/DOM,IOM,NOM,OMI /OBS/XO(51),YO(51),ZO(51),NOB
     *    /time/dt,period /ntime/nt
     *    /FLT/STR,DIP,RAK(503),FL,FW,FD(503),VR,TLc,Pn(503)
     *        ,XF(503),YF(503),TD(503),DP(503),fc(503),fmax(503),SM0L
     *    /Radiation/fR1,fR2,fran
     *    /IFLT/NF,NSTR,NDIP,NFfunc,ND,NSF,NDF /Nboore/Niter,idum
c
      nc=51
      if(nc.lt.NOB)then
        write(10,*)'from Boore:'
        WRITE(10,*)' NC < NOB !!!!!!!!!!!!!'
        write(*,*)'from Boore:'
        WRITE(*,*)' NC < NOB !!!!!!!!!!!!!'
        STOP
      END IF
c
      nc=203
      if(nc.lt.NF)then
        write(10,*)'from Indat:'
        WRITE(10,*)' NC < NF !!!!!!!!!!!!!'
        write(*,*)'from Indat:'
        WRITE(*,*)' NC < NF !!!!!!!!!!!!!'
        STOP
      END IF
c
c   * Input Data *
c
      read(7,*)
      read(7,*)
      read(7,*)fran,Niter,idum
      write(10,*)
      write(10,*)'From subroutine Boore:'
      write(10,*)'  Upper Limit Freq. for Random Phases (Hz) ', Fran
      write(10,*)'  Number of Iterrations =',Niter
      write(10,*)'  Initial Random Number (1-999) =',idum
      if(idum.lt.0)then
      idum0=idum
      idum=-idum
      write(10,*)'  Same random phases are used for all the sub-faults'
      end if
c
c   * Initial data
c
      pi=datan(1.d0)*4.d0
      dur=dt*dfloat(nt)
      nfast=nom*2
c
c   * Spectra of Nomalized Moment Rate Function *
c
      write(*,*)'  One Example of Phase Spectra is stored in phase.csv'
      open(21,file='Source(No.1)-phase.csv')
      write(21,*)'freq.(Hz),','phase,','envelope(+),','envelope(-)'
c
      write(*,*)'  Now calculating moment-rates for all the sub-faults!'
      do 100 iflt=1,NF
        write(*,*)' Fault No. =',iflt
        if((iflt.ge.2).and.(idum0.lt.0))goto 170
c
c   * Boore's Envelope Function *
c     (Tdu: the source duration (see eq.(6) in Boore (1983))
c
        Tdu=1.d0/fc(iflt)
        Tw=2.d0*Tdu
        b=1.2531d0
        c=5.d0*b/Tw
        a=(5.d0*2.718281829d0/Tw)**b
        bov(1)=0.d0
        bovN(1)=0.d0
        Tdu2=Tdu*2.d0
        Tdu3=Tdu*3.d0
        write(*,*)'Boore Envelope Function is stored in boore-env.csv'
        open(20,file='boore-env.csv')
        write(20,*)'time,','boore,','envelope'
        do i=2,nt
          time=dt*(i-1)
          bov(i)=a*(time**b)*dexp(-c*time)
          if(time.gt.Tdu2)then
            if(time.le.Tdu3)then
              bov(i)=bov(i)*dcos(pi*(time-Tdu2)/2.d0/Tdu)
            else
              bov(i)=0.d0
            end if
          end if
          bovN(i)=bov(i)**(1.d0/dfloat(Niter))
          if((bovN(i).lt.1.d-10).or.(Niter.eq.0))bovN(i)=0.d0
          write(20,9999)time,bov(i),bovN(i)
        end do
        close(20)
 9999   format(f12.4,5(',',e15.7))
c
        DO 110 irnd=1,2
c
          DO 120 IO=1,nom+1
            OM=DOM*DFLOAT(IO-1)
            f=om/(2.d0*pi)
c     * Brune's omega-square model *
     	    FFc=f/Fc(iflt)
	    S=1.d0/(1.d0+FFc*FFc)
c     * Fmax by Boore (1983) n=4 *
            FFmax=f/Fmax(iflt)
            P=1.d0/dsqrt(1.d0+FFmax**(2.d0*Pn(iflt)))
c     * Fmax by Kamae et al. (1991) *
c            P=1.d0/(1.d0+FFmax)
            SA=S*P
c     * Double derivatives of Moment Rate Function *
c            FX(io)=SA*cran1(io)*(-om*om)
            bos(io)=SA*om*om
c
c     * Generation of Phases (zero for lower and random for higher freq.)*
            ran=(ran0(idum)-0.5d0)*pi*2.d0
            phsav=pi
            if(f.lt.Fran)then
              amp=f/Fran
              amp=dsin(amp*pi/2.d0)
              amp=amp*amp
              amp=amp*amp
              ran=ran*amp
              phsav=phsav*amp
            end if
            dr=+dcos(ran)
            di=-dsin(ran)
            cran=dcmplx(dr,di)
            if((iflt.eq.1).and.(irnd.eq.1))then
              write(21,*)f,',',ran,',',phsav,',',-phsav
            end if
c
c     * Double Derivative of Moment Rate Function (Far-Field Acceleration) *
            FX(io)=SA*cran*(-om*om)
            FX1(io)=FX(io)
c
  120   CONTINUE
c
c     * Calculating First Waveforms (accelaration) using Boore's envelop funtion *
c
          DO I=2,NFAST/2
            FX(NFAST+2-I)=dCMPLX(dREAL(FX(I)),-DIMAG(FX(I)))
          end do
          CALL dWVFAST(NFAST,FX,MOM,1)
          DO I = 1,NFAST
            XR(I)=dREAL(FX(I))/DUR
            XR(i)=XR(i)*bov(i)
          end do
c
c     * Iterations to obtain Boore's waveform functions *
c
        do 150 iter=1,Niter
c
c       * Correction of Fourier Amplitude *
	  do i=1,nfast
            FX(I)=dCMPLX(XR(I),0.0d0)
	  end do
          CALL dWVFAST(nfast,FX,MOM,-1)
          DO K=1,nom+1
            FX(K)=FX(K)*dt
            FX(K)=FX(K)*cdabs(FX1(K))/cdabs(FX(K))
	  end do
c
c       * Waveform (accelaration) *
          DO I=2,NFAST/2
            FX(NFAST+2-I)=dCMPLX(dREAL(FX(I)),-DIMAG(FX(I)))
          end do
          CALL dWVFAST(NFAST,FX,MOM,1)
          DO 140 I = 1,NFAST
            XR(I)=dREAL(FX(I))/DUR
            XR(i)=XR(i)*bovN(i)
  140     CONTINUE
c
  150   CONTINUE
c
c   * Double derivatives of Momant Rate Function (Far Field acceleration) *
c
          DO I = 1,NFAST
            rate(i,iflt,irnd)=xr(i)
          end do
c
c   * Fourier Amplitude of Double derivatives of Momant Rate Function *
	  do i=1,nfast
            FX(I)=dCMPLX(rate(i,iflt,irnd),0.0d0)
	  end do
          CALL dWVFAST(nfast,FX,MOM,-1)
          DO K=1,nom+1
            Frate(K,iflt,irnd)=dconjg(FX(K))*dt
	  end do
c
  110   continue
c
  170   if(idum0.lt.0)then
          DO 180 irnd=1,2
            DO I = 1,NFAST
              rate(i,iflt,irnd)=rate(i,1,irnd)
            end do
            DO K=1,nom+1
              Frate(K,iflt,irnd)=Frate(K,1,irnd)
	    end do
  180     continue
        end if
c
  100 continue
      close(21)
c
c   * Data Output *
c
      write(*,*)
      write(*,*)'  Two sets of the double derivatives of the normalized'
      write(*,*)'  moment rates are stored in'
      write(*,*)'  Moment-Rate-Acc-1.csv and Moment-Rate-Acc-2.csv'
      open(22,file='Moment-Rate-Acc-1.csv')
        write(22,2000)(iflt,iflt=1,NF)
        do i=1,nfast
	  time=dt*(i-1)
          write(22,2001)time,(rate(i,iflt,1),iflt=1,NF)
        end do
      close(22)
      open(22,file='Moment-Rate-Acc-2.csv')
        write(22,2000)(iflt,iflt=1,NF)
        do i=1,nfast
	  time=dt*(i-1)
          write(22,2001)time,(rate(i,iflt,2),iflt=1,NF)
        end do
      close(22)
 2000 format('time(sec)',300(',Sub-Fault No.',I3))
 2001 format(f8.4,300(',',e17.5))
c
      write(*,*)'  Two-sets of Fourier amp. of moment-rates are stored'
      write(*,*)'  in Moment_Rate-Amp-1.csv and Moment_Rate-Amp-2.csv'
      open(23,file='Moment-Rate-AccAmp-1.csv')
        write(23,2100)(iflt,iflt=1,NF)
        DO K=1,nom+1
          om=dom*(k-1)
          f=om/(2.d0*pi)
c        write(23,2101)f,cdabs(fx1(k)),(cdabs(Frate(k,iflt,1)),iflt=1,NF)
        write(23,2101)f,bos(k),(cdabs(Frate(k,iflt,1)),iflt=1,NF)
        end do
      close(23)
      open(23,file='Moment-Rate-AccAmp-2.csv')
        write(23,2100)(iflt,iflt=1,NF)
        DO K=1,nom+1
          om=dom*(k-1)
          f=om/(2.d0*pi)
        write(23,2101)f,bos(k),(cdabs(Frate(k,iflt,2)),iflt=1,NF)
        end do
      close(23)
 2100 format('freq.(Hz),Theory for the last sub-fault',
     *    300(',Sub-Fault No.',I3))
 2101 format(f8.4,300(',',e15.7))
c
      return
      END
c
C
C *** DATA OUTPUT ***
C
      SUBROUTINE OTDAT(MOM,MOB,GX,GY,GZ)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 GX(MOM,MOB),GY(MOM,MOB),GZ(MOM,MOB)
      COMMON /FREQ/DOM,IOM,NOM,OMI /OBS/XO(51),YO(51),ZO(51),NOB
     *    /SIGN/ISIGN
C
C *** CHECK ***
C
      NC=51
      IF(NOB.GT.NC)THEN
         WRITE(10,*)'FROM OTDAT:NOB>NC',NOB,NC
         STOP
      END IF
C
C *** OUTPUT OF DATA (GX, GY AND GZ)***
C
      IF(ISIGN.EQ.1)THEN
        DO 100 IOB=1,NOB
        DO 110 K=IOM,NOM
          GX(K,IOB)=DCMPLX(DREAL(GX(K,IOB)),-DIMAG(GX(K,IOB)))
          GY(K,IOB)=DCMPLX(DREAL(GY(K,IOB)),-DIMAG(GY(K,IOB)))
          GZ(K,IOB)=DCMPLX(DREAL(GZ(K,IOB)),-DIMAG(GZ(K,IOB)))
  110   CONTINUE
  100   CONTINUE
      END IF
C
      OPEN(UNIT=8,FILE='grfault.axyz',STATUS='UNKNOWN')
        WRITE(8,2500)DOM,IOM,NOM
        WRITE(8,2510)NOB
        DO 500 IOB=1,NOB
          DO 510 K=IOM,NOM
            WRITE(8,2520)GX(K,IOB),GY(K,IOB),GZ(K,IOB)
  510     CONTINUE
  500   CONTINUE
      CLOSE(8)
C
 2500 FORMAT(1X,F10.6,2I10)
 2510 FORMAT(1X,I10,///)
 2520 FORMAT(1X,6E12.5)
C
      RETURN
      END
c
C *** INPUT OF Pole and Branch Points ***
C
      SUBROUTINE PHSBRC(MOM)
C
      IMPLICIT REAL*8(A-H,O-Z)
      dimension CR(18212,0:801),CL(18212,0:801),MDL(18212),MDR(18212)
      COMMON /FREQ/DOM,IOM,NOM,OMI /POLE/CS(18212,801),NW(18212)
     *       /STRC/DNS(20),VP(20),QP(20),FP(20)
     *                    ,VS(20),QS(20),FS(20),THK(20),NL
     *       /TMP1/CC(801)
C
C  ** CHECK ** 
C
      NC=18212
      IF(MOM.GT.NC)THEN
         WRITE(10,*)'FROM PHSBRC: MOM>NC !!'
         STOP
      END IF
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(10,*)'FROM PHSBRC: NL>NC !!'
         STOP
      END IF

c   * Number of Poles in Common Block*
      NC=801
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 PHSBRC: MDL>NC !!'
         STOP
        END IF
        IF(MDR(IO).GT.NC)THEN
         WRITE(10,*)'FROM PHSBRC: MDR>NC !!'
         STOP
        END IF
  100 CONTINUE
      close(22)
      close(23)

      open(25,file='phsl1_dat',status='old')
      open(26,file='phsr1_dat',status='old')
      DO 120 IO=IOM,NOM
        DO 110 MD=0,MDL(io)
          read(25,*)CL(IO,MD)
  110   CONTINUE
        DO 115 MD=0,MDR(io)
          read(26,*)CR(IO,MD)
  115   CONTINUE
  120 CONTINUE
      close(25)
      close(26)
C
C  ** CONSTRUCT PHASE VELOCITY DATA **
C
      write(10,1000)
 1000 FORMAT(//' From PHSBRC:')

      DO 200 IO=IOM,NOM

        IF(NL.EQ.1)THEN
          N=1
          CC(N)=CR(IO,0)
          GOTO 210
        END IF
C
        N=0
        IR=0
        IL=0
c
    5   N=N+1
        if((CR(IO,IR).lt.CL(IO,IL)).and.(CR(IO,IR).ne.0.d0))then
          CC(N)=CR(IO,IR)
          IR=IR+1
        else
          IF((CL(IO,IL).EQ.0.D0).AND.(CR(IO,IR).NE.0.D0))THEN
            CC(N)=CR(IO,IR)
            IR=IR+1
          ELSE
            if(CL(IO,IL).ne.0.d0)then
              CC(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 5
c
 210    NW(IO)=N+NL*2
c
        NC=801
        IF(NW(IO).GE.NC)THEN
          WRITE(10,*)' FROM PHSBRC: NW.GE.NC !!, IO',IO
          STOP
        END IF
C
        DO 220 I=1,N-1
          IF(CC(I).EQ.0.D0)THEN
            WRITE(10,*)' FROM PHSBRC: CS=0 !!, IO',IO
            STOP
          END IF
          IF(CC(I).GT.CC(I+1))THEN
            WRITE(10,*)' FROM PHSBRC: CS(I)>CS(I+1) !!, IO',IO
            STOP
          END IF
  220   CONTINUE
        IF(CC(N).EQ.0.D0)THEN
          WRITE(10,*)' FROM PHSBRC: CS=0 !!, IO',IO
          STOP
        END IF
c
        ib=0
        do 300 IL=1,NL
          VSR=VS(IL)
          VPR=VP(IL)
          is=0
          i=1
          C1=CC(i)
          if(VSR.lt.C1)then
            do ii=NW(IO),i+1,-1
              CC(ii)=CC(ii-1)
             end do
            CC(i)=VSR
            is=1
            ib=ib+1
            goto 20
          end if            
   10     i=i+1
          C2=CC(i)
          if(C2.eq.0.d0)goto 20
          if(VSR.lt.C2)then
            do ii=NW(IO),i+1,-1
              CC(ii)=CC(ii-1)
            end do
            CC(i)=VSR
            is=1
            ib=ib+1
            goto 20
          end if
          goto 10
c          
   20     if(is.eq.0)then
            CC(N+ib+1)=VSR
            ib=ib+1
          end if
c          
          ip=0
          i=1
   30     i=i+1
          C2=CC(i)
          if(C2.eq.0.d0)goto 40
          if(VPR.lt.C2)then
            do ii=NW(IO),i+1,-1
              CC(ii)=CC(ii-1)
            end do
            CC(i)=VPR
            ip=1
            ib=ib+1
            goto 40
          end if
          goto 30
c          
   40     if(ip.eq.0)then
            CC(N+ib+1)=VPR
            ib=ib+1
          end if
          
  300   continue
c          
        do i=1,nw(io)
          CS(IO,i)=CC(i)
          cc(i)=0.d0
        end do
c
c        write(10,1100)io*dom,nw(io)
c        write(10,1110)
c        do 320 i=1,NW(IO),5
c          write(10,1120)CS(IO,i),CS(IO,i+1)
c     *               ,CS(IO,i+2),CS(IO,i+3),CS(IO,i+4)
c  320   continue
c 1100   FORMAT('   Omega, Number of Phase Velocities=',f8.5,I5)
c 1110   FORMAT( '   Phase Velocities (m/sec) =')
c 1120   FORMAT(2x,5D12.4)
C
        DO 340 I=1,NW(IO)-1
          IF(CS(IO,I).EQ.0.D0)THEN
            WRITE(10,*)' FROM PHSBRC: CS=0 !!'
            WRITE(*,*)' FROM PHSBRC: CS=0 !!'
            STOP
          END IF
          IF(CS(IO,I).GT.CS(IO,I+1))THEN
            WRITE(10,*)' FROM PHSBRC: CS(I)>CS(I+1) !!'
            WRITE(*,*)' FROM PHSBRC: CS(I)>CS(I+1) !!'
            STOP
          END IF
  340   CONTINUE
        IF(CS(IO,N).EQ.0.D0)THEN
          WRITE(10,*)' FROM PHSBRC: CS=0 !!'
          WRITE(*,*)' FROM PHSBRC: CS=0 !!'
          STOP
        END IF
c
  200 continue
C
      RETURN
      END
C
C
C *** COMPUTATION OF GREEN FUNCTIONS USING THE R/T METHOD ***
c *** For large wavenumber * distance (k*r > 8), the Filon's integration 
c     formula is used instead of the Simpson's.
C
      SUBROUTINE HKQVW(MOM,MOB,MFL,Frate,GGX,GGY,GGZ)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 Vpc,Vsc,SM,RAM,SM1,XX1,SM2,XX2,SM3,XX3,RM1,RM2,RM3
     *      ,OM,VS2,VP2,Frate(MOM,MFL,2),cpx
     *      ,GGX(MOM,MOB),GGY(MOM,MOB),GGZ(MOM,MOB)
     *      ,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)
     *      ,SM0,SM00,EDT,QA,QB,VA,VB,WA,WB,GX,GY,GZ
     *      ,GXHp,GYHp,GZHp,GXHs,GYHs,GZHs
     *      ,GXHsh,GYHsh,GTHsh,GXHsv,GYHsv,GZHsv
     *      ,GXL,GYL,GZL,GXH,GYH,GZH
     *      ,QX(0:3),QY(0:3),QZ(0:3),VX(0:3),VY(0:3)
     *      ,WX(0:3),WY(0:3),WZ(0:3),UX(3,3),UY(3,3),UZ(3,3)
      DIMENSION EN(3,3),rd(3)
      COMMON /FREQ/DOM,IOM,NOM,OMI /OBS/XO(51),YO(51),ZO(51),NOB
     *    /STRC/DNS(20),VP(20),QP(20),FP(20)
     *                 ,VS(20),QS(20),FS(20),THK(20),NL
     *    /FLT/STR,DIP,RAK(503),FL,FW,FD(503),VR,TLc,Pn(503)
     *        ,XF(503),YF(503),TD(503),DP(503),fc(503),fmax(503),SM0L
      COMMON /IFLT/NF,NSTR,NDIP,NFfunc,ND,NSF,NDF
     *    /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 /outom/ioutom
     *    /INTG/VPMAX,VSMAX,RYMIN,CFINAL,N11,N12,M11
     *    /WVNUM/WVN(29998),DWN,WNR(3),NWN
     *    /Radiation/fR1,fR2,fran
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
      NC=503
      IF(NF.GE.NC)THEN
        WRITE(10,*)' FROM HKQVW:NF>NC !!!!!!!!!!'
        WRITE(10,*)' NF,NC=',NF,NC
        STOP
      END IF
C
C *** GREEN'S FUNCTION ***
C
C   * output of integrands for check *
C
      IF(IOT.eq.1)THEN
      OPEN (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
c *** Various Parameters ***
c
      CPX=(0.D0,1.D0)
      pai=4.d0*datan(1.d0)
      CST= DCOS(STR)
      SST= DSIN(STR)
      CDP= DCOS(DIP)
      SDP= DSIN(DIP)
c      CDP2=CDP*CDP
c      SDP2=SDP*SDP
c      CDP4= DCOS(4.d0*DIP)
c      CDP22=(1.d0+CDP4)/2.d0
c      SDP22=(1.d0-CDP4)/2.d0
c
c
c *** Homogeneous High-Freq. Radiation Coefficients for P and S waves ***
c
      write(10,*)
      FpH=dsqrt(4.d0/15.d0)
      FsH=dsqrt(0.4D0)
      WRITE(10,*)
      write(10,*)'Radiation Coefficients'
      WRITE(10,*)'Homogeneous Coefficients for P and S waves =',FpH,FsH
C
C *** GREEN'S FUNCTION ***
C
      do iob=1,nob
      do io=1,nom
        GGX(IO,IOB)=(0.D0,0.D0)
        GGY(IO,IOB)=(0.D0,0.D0)
        GGZ(IO,IOB)=(0.D0,0.D0)
      end do
      end do
c
c   * Green's function computations *
c
      write(*,*)
      write(*,*)'From subroutine HKQVW:'
      write(*,*)'  Now calculating ground motions for each frequency!'
      write(10,*)
      write(10,*)'Radiation Coefficients from subroutine HKQVW:'
c      sNF=dsqrt(dfloat(NF))
      sNF=dfloat(ND)
      DO 100 IO=IOM,NOM
C
        OMR=DOM*DFLOAT(IO)
        OM=dcmplx(OMR,OMI)
        freq=omr/(2.d0*pai)
        OMTd=OMR*TLc
        sinc=dsin(OMTd)/OMTd
        if((iot.eq.1).and.(io.ne.ioutom))goto 100
C
c     * Configure wavenumbers for numerical integrations *
c
c       * range 1: from near 0 to near om/Rayleigh(min) *
        WNR(2)=OMR/rymin*1.2d0
        WNR(1)=WNR(2)/n11/10000.d0
        CALL WVNUMV(MOM,IO,dreal(OM))
c
        NWN=2*(M11+N12)+1
        NC=29998
        IF(NWN.GT.NC)THEN
          WRITE(10,*)'FROM HKQVW: NWN > NC !!!',NWN,NC
        STOP
        END IF
C
c       * range 2: from near om/Rayleigh(min) to near om/C(final) *
        WNR(3)=OMR/cfinal
        DWN=(WNR(3)-WNR(2))/(2.D0*N12)
        DO 150 K=1,2*N12+1
          WVN(2*M11+K)=WNR(2)+DWN*(K-1)
 150    CONTINUE
C
        write(*,1005)freq,dfloat(io)/dfloat(nom)*100.d0,nwn
        write(10,*)
        write(10,1005)om,nwn
        write(10,1010)wnr(1),wnr(2),wnr(3)
        write(10,1020)dwn
 1005   format('   Freq (%) =',f7.3,' (',f5.1,' %),'
     *        ,' Num. of Wave-number integ. points =',I5)
 1010   format(' Boundaries for the Integration ranges:',3e12.4)
 1020   format(' Integration Increments for Range 2:',e12.4)
C
c     * the observation points *
c
        DO 200 IOB=1,NOB
c
          ZS=ZO(IOB)
          if(iob.ge.2)then
            if(zs.eq.zo(iob-1))goto 10
          end if
c
          DRU=0.d0
          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.d0
          DO 230 L=1,LS-1
            DEP1=DEP1+THK(L)
 230      CONTINUE
          SDEP=ZS-DEP1
          DEP2=DEP1+THK(LS)
C
          QPe=FuncQP(OM,LS)
          QSe=FuncQS(OM,LS)
          VPc=VP(LS)*DCMPLX(1.D0,-0.5D0/QPe)
          VSc=VS(LS)*DCMPLX(1.D0,-0.5D0/QSe)
          VP2= VPc*VPc
          VS2= VSc*VSc
          SM2= DNS(LS)*VS2
          RM2= DNS(LS)*VP2-2.d0*SM2
          XX2 = (1.d0+VS2/VP2)/(1.d0-VS2/VP2)
          IF(LS.GE.2)THEN
            QPe=FuncQP(OM,LS-1)
            QSe=FuncQS(OM,LS-1)
            VPc=VP(LS-1)*DCMPLX(1.D0,-0.5D0/QPe)
            VSc=VS(LS-1)*DCMPLX(1.D0,-0.5D0/QSe)
            VP2= VPc*VPc
            VS2= VSc*VSc
            SM1= DNS(LS-1)*VS2
            RM1= DNS(LS-1)*VP2-2.d0*SM1
            XX1 = (1.d0+VS2/VP2)/(1.d0-VS2/VP2)
          END IF
          IF(LS.NE.NL)THEN
            QPe=FuncQP(OM,LS+1)
            QSe=FuncQS(OM,LS+1)
            VPc=VP(LS+1)*DCMPLX(1.D0,-0.5D0/QPe)
            VSc=VS(LS+1)*DCMPLX(1.D0,-0.5D0/QSe)
            VP2= VPc*VPc
            VS2= VSc*VSc
            SM3= DNS(LS+1)*VS2
            RM3= DNS(LS+1)*VP2-2.d0*SM3
            XX3 = (1.d0+VS2/VP2)/(1.d0-VS2/VP2)
          END IF
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
 10       DO 300 is=1,NF
c
            if(FD(IS).eq.0.d0)goto 300
c
            RX=XF(IS)-XO(IOB)
            RY=YF(IS)-YO(IOB)
            RR=DSQRT(RX*RX+RY*RY)
            RZ=DP(IS)-ZO(IOB)
            RRR=DSQRT(RX*RX+RY*RY+RZ*RZ)
            rd(1)=Rx/RRR
            rd(2)=Ry/RRR
            rd(3)=Rz/RRR
            CRK= DCOS(RAK(is))
            SRK= DSIN(RAK(is))
c            CRK2=CRK*CRK
c            SRK2=SRK*SRK
c
              zd=dp(is)
              DRU=0.d0
              DRL=THK(1)
              do L=1,NL
                IF((zd.GE.DRU).AND.(zd.LT.DRL))L1=L
                IF((zd.GE.DRU).AND.(L.EQ.NL))  L1=L
                DRU=DRL
                DRL=DRL+THK(L+1)
              end do
c
              LO=L1
              QPe=FuncQP(OM,LO)
              QSe=FuncQS(OM,LO)
              VPc=VP(LO)*DCMPLX(1.D0,-0.5D0/QPe)
              VSc=VS(LO)*DCMPLX(1.D0,-0.5D0/QSe)
              VP2= VPc*VPc
              VS2= VSc*VSc
              SM = DNS(LO)*VS2
              RAM= DNS(LO)*VP2-2.D0*SM
C
              ZZO=zd
              ODEP=ZZO
              DO L=2,LO
                ODEP=ODEP-THK(L-1)
              end do
c
c    *      Time Delay *
c
              EDT=CDEXP(cpx*OM*TD(IS))
c
c    *      Seismic Moment *
            SM00=FD(IS)*SM*FL*FW/sNF
c    *      Slip Velocity Function -> Slip Displacement Function *
            SM0=SM00*CPX/OM
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         * see Hisada(1994) eq.(26) *
                CR=RX/RR
                SR=RY/RR
                C2=CR*CR
                S2=SR*SR
                CS=CR*SR
                CX=S2/RR
                SY=C2/RR
                CY=-CS/RR
                SX=CY
C       
c         * Urx, Ury, Urz: see Hisada(1994) eq.(11) *
                QA=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         * 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    * Low Frequency Wave *
c
            XN=-SDP*SST
            YN= SDP*CST
            ZN=-CDP
            EX= CRK*CST+SRK*CDP*SST
            EY= CRK*SST-SRK*CDP*CST
            EZ=-SRK*SDP
C
            EN(1,1)= 2.D0*EX*XN
            EN(1,2)= EX*YN +EY*XN
            EN(1,3)= EX*ZN +EZ*XN
            EN(2,1)= EY*XN +EX*YN
            EN(2,2)= 2.D0*EY*YN
            EN(2,3)= EY*ZN +EZ*YN
            EN(3,1)= EZ*XN +EX*ZN
            EN(3,2)= EZ*YN +EY*ZN
            EN(3,3)= 2.D0*EZ*ZN
c
              GXL=(0.D0,0.D0)
              GYL=(0.D0,0.D0)
              GZL=(0.D0,0.D0)
              DO I=1,3
              DO J=1,3
                GXL=GXL+EN(I,J)*UX(J,I)
                GYL=GYL+EN(I,J)*UY(J,I)
                GZL=GZL+EN(I,J)*UZ(J,I)
              end do
              end do
c
c    * High Frequency Wave *
             GXH=(0.D0,0.D0)
             GYH=(0.D0,0.D0)
             GZH=(0.D0,0.D0)
             GXHp=(0.D0,0.D0)
             GYHp=(0.D0,0.D0)
             GZHp=(0.D0,0.D0)
             GXHs=(0.D0,0.D0)
             GYHs=(0.D0,0.D0)
             GZHs=(0.D0,0.D0)
c
c      * High Frequency P Wave (Explosive Source) *
c
             if((NSF.eq.2).or.(NSF.eq.3))goto 555
c
             do i=1,3
               GXHp=GXHp+UX(i,i)
               GYHp=GYHp+UY(i,i)
               GZHp=GZHp+UZ(i,i)
             end do
             GXHp=FpH*GXHp
             GYHp=FpH*GYHp
             GZHp=FpH*GZHp
c
             if(NSF.eq.1)goto 666
c
c      * Check of the SH and SV Waves Radiation Coefficients *
c
c             FshH=2.d0*CRK2*(1.d0+SDP2)+SRK2*(1.d0+CDP22)
c             FshH=dsqrt(FshH/3.d0)/2.d0
c             FsvH=CRK2*(4.d0+10.d0*CDP2)+SRK2*(14.d0+5.d0*SDP22)
c             FsvH=dsqrt(FsvH/15.d0)/2.d0
c             FsXL=re*XN+(EX-2.d0*re*rd(1))*rn
c             FsYL=re*YN+(EY-2.d0*re*rd(2))*rn
c             FshL= FsXL*rd(2)-FsYL*rd(1)
c             FsvL=-FsXL*rd(1)-FsYL*rd(2)
c             if(FshL.lt.0.d0)FshH=-FshH
c             if(FsvL.lt.0.d0)FsvH=-FsvH
c
c      * High Frequency SH Wave (Fault Plane is perpendicular to the OP-SP line) *
c
  555        GXHsh=(0.D0,0.D0)
             GYHsh=(0.D0,0.D0)
             if(NSF.eq.3)goto 565
c
c      * Unit Vector of the Fault Plane *
             XN=-rd(1)
             YN=-rd(2)
             ZN=-rd(3)
             XY2=dsqrt(XN*XN+YN*YN)
             XNS=XN/XY2
             YNS=YN/XY2
C
c          * Unit Slip Vector of the SH Source *
             EX= YNS
             EY=-XNS
             EZ= 0.d0
c
c      * (Ei x Nj + Ej x Ni) and SH wave Green function *
c
             EN(1,1)= 2.D0*EX*XN
             EN(1,2)= EX*YN +EY*XN
             EN(1,3)= EX*ZN +EZ*XN
             EN(2,1)= EY*XN +EX*YN
             EN(2,2)= 2.D0*EY*YN
             EN(2,3)= EY*ZN +EZ*YN
             EN(3,1)= EZ*XN +EX*ZN
             EN(3,2)= EZ*YN +EY*ZN
             EN(3,3)= 2.D0*EZ*ZN
c
              DO I=1,3
              DO J=1,3
                GXHsh=GXHsh+EN(I,J)*UX(J,I)
                GYHsh=GYHsh+EN(I,J)*UY(J,I)
              end do
              end do
c
c          * Transverse Component and Distribution to X and Y components *
             GTHsh=-YNS*GXHsh+XNS*GYHsh
             GXHsh=FsH*GTHsh/dsqrt(2.d0)
             GYHsh=GXHsh
c
c      * High Frequency SV Wave (Fault Plane is perpendicular to the OP-SP line) *
c
  565        GXHsv=(0.D0,0.D0)
             GYHsv=(0.D0,0.D0)
             GZHsv=(0.D0,0.D0)
             if(NSF.eq.2)goto 575
c
c      * Unit Vector of the Fault Plane *
             XN=-rd(1)
             YN=-rd(2)
             ZN=-rd(3)
             XY2=dsqrt(XN*XN+YN*YN)
             XNS=XN/XY2
             YNS=YN/XY2
c
c          * Unit Slip Vector of the SV Source *
             EX= ZN*XNS
             EY= ZN*YNS
             EZ= -dsqrt(1-ZN*ZN)
c
c      * (Ei x Nj + Ej x Ni) and SH wave Green function *
c
             EN(1,1)= 2.D0*EX*XN
             EN(1,2)= EX*YN +EY*XN
             EN(1,3)= EX*ZN +EZ*XN
             EN(2,1)= EY*XN +EX*YN
             EN(2,2)= 2.D0*EY*YN
             EN(2,3)= EY*ZN +EZ*YN
             EN(3,1)= EZ*XN +EX*ZN
             EN(3,2)= EZ*YN +EY*ZN
             EN(3,3)= 2.D0*EZ*ZN
c
          write(*,*)'0:GXHsv',GXHsv
              DO I=1,3
              DO J=1,3
                GXHsv=GXHsv+EN(I,J)*UX(J,I)
                GYHsv=GYHsv+EN(I,J)*UY(J,I)
                GZHsv=GZHsv+EN(I,J)*UZ(J,I)
              end do
              end do
c
c          * Radial Component and Distribution to X and Y components *
             GXHsv=FsH*GXHsv
             GYHsv=FsH*GYHsv
             GZHsv=FsH*GZHsv
c
c          * Sum of SH and SV and Amplitudes of X and Y are set to be equal *
  575        GXHs=GXHsh+GXHsv
             GYHs=GYHsh+GYHsv
             GZHs=GZHsv
c
c    * Hybrid Wave *
c
  666        GXH=GXHp+GXHs
             GYH=GYHp+GYHs
             GZH=GZHp+GZHs
             if(NDF.eq.0)then
               aGXH=cdabs(GXH)
               aGYH=cdabs(GYH)
               aGHH=sqrt(aGXH*aGXH+aGYH*aGYH)/dsqrt(2.d0)
               GXH=aGHH*GXH/aGXH
               GYH=aGHH*GYH/aGYH
             end if
c
             if(freq.lt.fR1)then
                GX=GXL
                GY=GYL
                GZ=GZL
             else
                if(freq.lt.fR2)then
                  dfreq=(freq-fR1)/(fR2-fR1)
                  GX=GXL+(GXH-GXL)*dfreq
                  GY=GYL+(GYH-GYL)*dfreq
                  GZ=GZL+(GZH-GZL)*dfreq
                else
                  GX=GXH
                  GY=GYH
                  GZ=GZH
                end if
             end if
c
c    * Irikura's F-function *
            if(NFfunc.eq.0)then
              EDT=EDT*(1.d0+(sNF-1.d0)*sinc*exp(cpx*OMTd))
c    * Exp-type F-function *
            else
              EDT=EDT*(1.d0+(sNF-1.d0)/(1.d0-cpx*OMTd))
            end if
c
            GX=GX*SM0*EDT*Frate(io+1,is,1)
            GY=GY*SM0*EDT*Frate(io+1,is,2)
            GZ=GZ*SM0*EDT*(Frate(io+1,is,1)+Frate(io+1,is,2))/2.d0
            GGX(IO,IOB)=GGX(IO,IOB)+GX
            GGY(IO,IOB)=GGY(IO,IOB)+GY
            GGZ(IO,IOB)=GGZ(IO,IOB)+GZ

c
  300     CONTINUE
C
  200   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
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,M11
     *       /WVNUM/WVN(29998),DWN,WNR(3),NWN
     *       /POLE/CS(18212,801),NW(18212) /TMP/WVT(29998)
C
C  ** CHECK ** 
C
      NC=18212
      IF(MOM.GT.NC)THEN
         WRITE(10,*)'FROM WVNUM: MOM > NC !!'
         STOP
      END IF
C
      NC=801
      IF(NW(IO).GT.NC)THEN
         WRITE(10,*)'FROM WVNUM: NW(IO) > NC !!'
         STOP
      END IF
C
      NC=29998
      IF((N11+1+NW(IO)*9).GT.NC)THEN
         WRITE(10,*)'FROM WVNUM: N11+1+NW(IO)*9 > NC !!'
	   write(10,*)'N11,NW(io),io,NC',N11,NW(io),io,NC
         STOP
      END IF
C
C  ** POINTS **
C
      DW1=(WNR(2)-WNR(1))/N11
      DO 100 K=1,N11+1
        WVT(K)=WNR(1)+DW1*(K-1)
 100  CONTINUE
c
      MW=N11
      DW=DW1/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
c
        IF(IPL.EQ.0)THEN
          W1=POLE-DW1/2.D0
          CALL WVNM(W1,MW)
          W2=POLE-DW1/4.D0
          CALL WVNM(W2,MW)
          W3=POLE-DW1/8.D0
          CALL WVNM(W3,MW)
          W4=POLE-DW1/16.D0
          CALL WVNM(W4,MW)
          W5=POLE
          CALL WVNM(W5,MW)
          W6=POLE+DW1/16.D0
          CALL WVNM(W6,MW)
          W7=POLE+DW1/8.D0
          CALL WVNM(W7,MW)
          W8=POLE+DW1/4.D0
          CALL WVNM(W8,MW)
          W9=POLE+DW1/2.D0
          CALL WVNM(W9,MW)
        ELSE
          CALL WVNM(POLE,MW)
        END IF
 200  CONTINUE
c
      DW=DW1/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(2)).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
c
      DO 400 J=1,MW
        WVN(2*J-1)=WVT(J)
        WVN(2*J  )=(WVT(J)+WVT(J+1))/2.D0
 400  CONTINUE
      WVN(2*MW+1)=WVT(MW+1)
c
      M11=MW
c
      return
      end
C
C
      SUBROUTINE WVNM(W,MW)
      implicit real*8(a-h,o-z)
      COMMON /TMP/WVT(29998)
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
c
      RETURN
      END
c
C
C *** THE COMPUTATION OF Q, V, AND W AND THEIR DERIVATIVES 
C     USING SIMPSON'S INTEGATION FOMULA (before k*r < 8)***
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),OM
     *          ,QXB(0:2),QZB(0:2),VXB(0:2),WXB(0:2),WZB(0:2),SM
     *          ,QXC(0:2),QZC(0:2),VXC(0:2),WXC(0:2),WZC(0:2),RAM
      COMMON /INTG/VPMAX,VSMAX,RYMIN,CFINAL,N11,N12,M11
     *       /WVNUM/WVN(29998),DWN,WNR(3),NWN
     *       /FCQVW/SM,RAM,RR  /IRANGE/NRNG,IRNG
C
c *** Dimension Check ***
c
      NC=29998
      IF(NWN.GT.NC)THEN
         WRITE(10,*)' FROM QVWRZS: NWN.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C   * RANGE 1 : From near Zero to near om/Rayleigh(min) *
C
      NRNG=0
      IRNG=0
C
      WA=WNR(1)
      K=1
      CALL INTEGS(OM,K,WA,QXA,QZA,VXA,WXA,WZA)
      DO  100 I=1,M11
        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)
c
        H=WC-WB
        DO 110 J=0,2
         GQX(J)= GQX(J)+H*(QXA(J)+4.d0*QXB(J)+QXC(J))/3.d0
         GVX(J)= GVX(J)+H*(VXA(J)+4.d0*VXB(J)+VXC(J))/3.d0
         GWX(J)= GWX(J)+H*(WXA(J)+4.d0*WXB(J)+WXC(J))/3.d0
         GQZ(J)= GQZ(J)+H*(QZA(J)+4.d0*QZB(J)+QZC(J))/3.d0
         GWZ(J)= GWZ(J)+H*(WZA(J)+4.d0*WZB(J)+WZC(J))/3.d0
         QXA(J)= QXC(J)
         VXA(J)= VXC(J)
         WXA(J)= WXC(J)
         QZA(J)= QZC(J)
         WZA(J)= WZC(J)
  110   CONTINUE
c
        IF(RR*WC.GT.8.D0)THEN
          NRNG=1
          IRNG=I
          IF(I.EQ.M11)THEN
            NRNG=2
            IRNG=0
          END IF
          RETURN
        END IF
c
  100 CONTINUE
c
C   * RANGE 2 : From near om/Rayleigh(min) to near om/C(final) *
C
      H=DWN
      WA=WNR(2)
      DO  200 I=1,N12
        K=2*(M11+I)
        WB=WVN(K)
        CALL 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)
c
        DO 210 J=0,2
         GQX(J)= GQX(J)+H*(QXA(J)+4.d0*QXB(J)+QXC(J))/3.d0
         GVX(J)= GVX(J)+H*(VXA(J)+4.d0*VXB(J)+VXC(J))/3.d0
         GWX(J)= GWX(J)+H*(WXA(J)+4.d0*WXB(J)+WXC(J))/3.d0
         GQZ(J)= GQZ(J)+H*(QZA(J)+4.d0*QZB(J)+QZC(J))/3.d0
         GWZ(J)= GWZ(J)+H*(WZA(J)+4.d0*WZB(J)+WZC(J))/3.d0
         QXA(J)= QXC(J)
         VXA(J)= VXC(J)
         WXA(J)= WXC(J)
         QZA(J)= QZC(J)
         WZA(J)= WZC(J)
  210   CONTINUE
c
        IF(RR*WC.GT.8.D0)THEN
          NRNG=2
          IRNG=I
          IF(I.EQ.N12)THEN
            NRNG=3
            IRNG=0
          END IF
          RETURN
        END IF

  200 CONTINUE
C
      NRNG=3
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,OM
      DIMENSION BJ0(0:1),BJ1(0:1),BJK(0:1),BJD(0:1)
      COMMON /SOLY/SDEP,ODEP,LS,LO 
     *       /STRC/DNS(20),VP(20),QP(20),FP(20)
     *                    ,VS(20),QS(20),FS(20),THK(20),NL
     *       /OTST/IOT,IST /FCQVW/SM,RAM,RR
C
C   * CHECK OF DIMENSION *
C
      NC=20
      IF(NL.GE.NC)THEN
        WRITE(10,*)' FROM INTEGS:NL>NC !!!!!!!!!!'
        WRITE(10,*)' NL,NC=',NL,NC
        STOP
      END IF
C
c   * the displacement-stress vectors *
c
      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.d0*SM)
      URZZ(2)=(URZ(2,2)-WN*RAM*URZ(1,1))/(RAM+2.d0*SM)
C
c   * the Bessel functions *
c
      RK=WN*RR
      call bj01(RK,BJ0(0),BJ1(0))
      IF(RK.EQ.0.0)THEN
        BJK(0)=0.5d0
      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.d0
      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 > 8.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,OM,a1,b1
      COMMON /FCQVW/SM,RAM,RR  /IRANGE/NRNG,IRNG
     *       /INTG/VPMAX,VSMAX,RYMIN,CFINAL,N11,N12,M11
     *       /WVNUM/WVN(29998),DWN,WNR(3),NWN
C
C *** Dimesion Check ***
C
      NC=29998
      IF(NWN.GT.NC)THEN
         WRITE(10,*)' FROM QVWRZF: NWN.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C   * RANGE 1 : From near Zero to near om/Rayleigh(min) *
C
      IF(NRNG.EQ.3)RETURN

      JRNG=0
      IF(NRNG.GT.1)GOTO 10
C
      K=2*IRNG+1
      WA=WVN(K)
c      write(*,*)'QVWRZF(Range 1):WA',WA
      CALL INTEGF(OM,K,WA,QXSA,QZSA,VXSA,WXSA,WZSA
     *                   ,QXCA,QZCA,VXCA,WXCA,WZCA)
c
      DO  300 I=IRNG+1,M11
c
        K=2*I
        WB=WVN(K)
        CALL 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)
c
        H=WC-WB
        th=RR*H
        t=1.d0/th
        t2=t*t
        t3=t2*t
        si=dsin(th)
        ss=si*si
        ci=dcos(th)
        cc=ci*ci
        s2=dsin(2.d0*th)
c
        a1=t+0.5d0*s2*t2-2.d0*ss*t3
        b1=2.d0*((1.d0+cc)*t2-s2*t3)
        c1=4.d0*(si*t3-ci*t2)
c
        SA=1.D0/DTAN(RR*WA)
        CA=DTAN(RR*WA)
        SC=1.D0/DTAN(RR*WC)
        CC=DTAN(RR*WC)
        DO 350 J=0,2
          QXS1(J)=(QXSA(J)+QXSC(J))*0.5d0
          QZS1(J)=(QZSA(J)+QZSC(J))*0.5d0
          VXS1(J)=(VXSA(J)+VXSC(J))*0.5d0
          WXS1(J)=(WXSA(J)+WXSC(J))*0.5d0
          WZS1(J)=(WZSA(J)+WZSC(J))*0.5d0
          QXC1(J)=(QXCA(J)+QXCC(J))*0.5d0
          QZC1(J)=(QZCA(J)+QZCC(J))*0.5d0
          VXC1(J)=(VXCA(J)+VXCC(J))*0.5d0
          WXC1(J)=(WXCA(J)+WXCC(J))*0.5d0
          WZC1(J)=(WZCA(J)+WZCC(J))*0.5d0
c
          S1=A1*(QXSA(J)*SA-QXSC(J)*SC)+B1*QXS1(J)+C1*QXSB(J)
          S2=A1*(QXCC(J)*CC-QXCA(J)*CA)+B1*QXC1(J)+C1*QXCB(J)
          GQX(J)=GQX(J)+H*(S1+S2)
          S1=A1*(QZSA(J)*SA-QZSC(J)*SC)+B1*QZS1(J)+C1*QZSB(J)
          S2=A1*(QZCC(J)*CC-QZCA(J)*CA)+B1*QZC1(J)+C1*QZCB(J)
          GQZ(J)=GQZ(J)+H*(S1+S2)
          S1=A1*(VXSA(J)*SA-VXSC(J)*SC)+B1*VXS1(J)+C1*VXSB(J)
          S2=A1*(VXCC(J)*CC-VXCA(J)*CA)+B1*VXC1(J)+C1*VXCB(J)
          GVX(J)=GVX(J)+H*(S1+S2)
          S1=A1*(WXSA(J)*SA-WXSC(J)*SC)+B1*WXS1(J)+C1*WXSB(J)
          S2=A1*(WXCC(J)*CC-WXCA(J)*CA)+B1*WXC1(J)+C1*WXCB(J)
          GWX(J)=GWX(J)+H*(S1+S2)
          S1=A1*(WZSA(J)*SA-WZSC(J)*SC)+B1*WZS1(J)+C1*WZSB(J)
          S2=A1*(WZCC(J)*CC-WZCA(J)*CA)+B1*WZC1(J)+C1*WZCB(J)
          GWZ(J)=GWZ(J)+H*(S1+S2)
c
          QXSA(J)=QXSC(J)
          QZSA(J)=QZSC(J)
          VXSA(J)=VXSC(J)
          WXSA(J)=WXSC(J)
          WZSA(J)=WZSC(J)
          QXCA(J)=QXCC(J)
          QZCA(J)=QZCC(J)
          VXCA(J)=VXCC(J)
          WXCA(J)=WXCC(J)
          WZCA(J)=WZCC(J)
c
          WA=WC
  350   CONTINUE
c
  300 CONTINUE

c      write(*,*)'f1:gqx0',gqx(0)
C
C   * RANGE 2 : From om/Rayleigh(min) to om/C(final) *
C
      JRNG=1
      IRNG=0
C
 10   H=DWN
      K=2*(M11+IRNG)+1
      WA=WVN(K)
c      write(*,*)'QVWRZF(Range 2):WA',WA
      IF(JRNG.EQ.1)GOTO 15
      CALL INTEGF(OM,K,WA,QXSA,QZSA,VXSA,WXSA,WZSA
     *                   ,QXCA,QZCA,VXCA,WXCA,WZCA)

 15   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,N12
        K=2*(M11+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(3))
      CC=DTAN(RR*WNR(3))
      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(*,*)'f2: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,OM
      DIMENSION BJ0S(0:1),BJ1S(0:1),BJKS(0:1),BJDS(0:1)
     *         ,BJ0C(0:1),BJ1C(0:1),BJKC(0:1),BJDC(0:1)
      COMMON /SOLY/SDEP,ODEP,LS,LO 
     *       /STRC/DNS(20),VP(20),QP(20),FP(20)
     *                    ,VS(20),QS(20),FS(20),THK(20),NL
     *       /OTST/IOT,IST /FCQVW/SM,RAM,RR
C
C *** CHECK OF DIMENSION ***
C
      NC=20
      IF(NL.GE.NC)THEN
        WRITE(10,*)' FROM INTEGF:NL>NC !!!!!!!!!!'
        WRITE(10,*)' NL,NC=',NL,NC
        STOP
      END IF
C
c   * the displacement-stress vectors *
c
      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.d0*SM)
      URZZ(2)=(URZ(2,2)-WN*RAM*URZ(1,1))/(RAM+2.d0*SM)
C
c   * SIN and COS parts of the approximate Bessel functions *
c
      PI=4.D0*DATAN(1.D0)
      RK=WN*RR
      A=1.D0/DSQRT(PI*RK)
      SRK=A*DSIN(RK)
      CRK=A*DCOS(RK)
c
      x1=0.125d0/rk
      x2=0.0703125d0/(rk*rk)
      x3=0.073242188d0/(rk*rk*rk)      
      BJ0S(0)= SRK*(1.D0+x1-x2-x3)
      BJ0C(0)= CRK*(1.D0-x1-x2+x3)

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

      BJKS(1)=(BJ1S(1)-BJ1S(0)/RR)/RK
      BJKC(1)=(BJ1C(1)-BJ1C(0)/RR)/RK
      BJDS(1)= BJ0S(1)-BJKS(1)
      BJDC(1)= BJ0C(1)-BJKC(1)
C
c   * the integrands *
c
      QXS(0)= URX(1,1)*BJDS(0)+VLX(1)*BJKS(0)
      QXC(0)= URX(1,1)*BJDC(0)+VLX(1)*BJKC(0)
      VXS(0)= URX(1,1)*BJKS(0)+VLX(1)*BJDS(0)
      VXC(0)= URX(1,1)*BJKC(0)+VLX(1)*BJDC(0)
      WXS(0)= URX(2,1)*BJ1S(0)
      WXC(0)= URX(2,1)*BJ1C(0)
      QZS(0)= URZ(1,1)*BJ1S(0)
      QZC(0)= URZ(1,1)*BJ1C(0)
      WZS(0)= URZ(2,1)*BJ0S(0)
      WZC(0)= URZ(2,1)*BJ0C(0)
C
      QXS(1)= URX(1,1)*BJDS(1)+VLX(1)*BJKS(1)
      QXC(1)= URX(1,1)*BJDC(1)+VLX(1)*BJKC(1)
      VXS(1)= URX(1,1)*BJKS(1)+VLX(1)*BJDS(1)
      VXC(1)= URX(1,1)*BJKC(1)+VLX(1)*BJDC(1)
      WXS(1)= URX(2,1)*BJ1S(1)
      WXC(1)= URX(2,1)*BJ1C(1)
      QZS(1)= URZ(1,1)*BJ1S(1)
      QZC(1)= URZ(1,1)*BJ1C(1)
      WZS(1)= URZ(2,1)*BJ0S(1)
      WZC(1)= URZ(2,1)*BJ0C(1)
C
      QXS(2)= URXZ(1)*BJDS(0)+VLXZ*BJKS(0)
      QXC(2)= URXZ(1)*BJDC(0)+VLXZ*BJKC(0)
      VXS(2)= URXZ(1)*BJKS(0)+VLXZ*BJDS(0)
      VXC(2)= URXZ(1)*BJKC(0)+VLXZ*BJDC(0)
      WXS(2)= URXZ(2)*BJ1S(0)
      WXC(2)= URXZ(2)*BJ1C(0)
      QZS(2)= URZZ(1)*BJ1S(0)
      QZC(2)= URZZ(1)*BJ1C(0)
      WZS(2)= URZZ(2)*BJ0S(0)
      WZC(2)= URZZ(2)*BJ0C(0)
C
      IF(IOT.eq.1)THEN
        WRITE(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 VSc,CD,CU,VS2,WS2,CMS,SM,EE(2,2),DD,DU
     *  ,CUSU,CDSU,CDSD,CUSD,GTD,GTU,GRD,GRU,AU,BU,AD,BD,CS,OM,OM2
      COMMON /SOLY/SDEP,ODEP,LS,LO 
     *       /STRC/DNS(20),VP(20),QP(20),FP(20)
     *                    ,VS(20),QS(20),FS(20),THK(20),NL
     *       /DUSH/CD(29998,20),CU(29998,20)
     *       ,CUSU(29998),CDSU(29998),CDSD(29998),CUSD(29998)
     *       /GRTSH/GTD(0:20),GTU(0:20),GRD(0:20),GRU(0:20)
     *       /WVNUM/WVN(29998),DWN,WNR(3),NWN
C
C *** CHECK FOR THE DIMENSION OF CD, CU, AND GTD-GRU ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(10,*)' FROM CUDSH: NL.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
      NC=29998
      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

      QSe=FuncQS(OM,LS)
      VSc=VS(LS)*DCMPLX(1.D0,-0.5D0/QSe)
      VS2=VSc*VSc
      WS2=OM2/VS2
      SM=DNS(LS)*VS2
      EE(1,1)=(1.d0,0.d0)
      EE(1,2)=(1.d0,0.d0)
C
      DO 100 I=1,NWN
C
        WN=WVN(I)
C
C     * COMPUTATION OF GENERALIZED R/T MATRICES *
C
        CALL GNRTSH(OM,WN)
C
C     * COMPUTATION OF  CD AND CU *
C
        WN2=WN*WN
        CMS=CDSQRT(WN2-WS2)
        IF(dREAL(CMS).LT.0.d0)CMS=-CMS
C
        EE(2,1)=-SM*CMS
        EE(2,2)= SM*CMS
C
        PI=4.D0*DATAN(1.D0)
        DQ=-WN/(2.d0*PI)
C
        IF((LS.EQ.1).AND.(SDEP.EQ.0.d0))THEN
C
          IF(NL.EQ.1)THEN
            CDSD(I)=DQ/EE(2,1)
            CUSD(I)=(0.d0,0.d0)
          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 CMP,CMS,EE(4,4),DD(2,2),DU(2,2)
     *    ,CDX,CUX,CUSUX,CDSUX,CDSDX,CUSDX,GTD,GTU,GRD,GRU,DT
     *    ,CDZ,CUZ,CUSUZ,CDSUZ,CDSDZ,CUSDZ,OM
     *    ,AU(2,2),BU(2,2),AD(2,2),BD(2,2),CS(2,2),E1(2,2),E2(2,2)
      COMMON /SOLY/SDEP,ODEP,LS,LO 
     *  /STRC/DNS(20),VP(20),QP(20),FP(20)
     *               ,VS(20),QS(20),FS(20),THK(20),NL
     *  /DUPSVX/CDX(29998,2,20),CUX(29998,2,20),CUSUX(29998,2)
     *         ,CDSUX(29998,2),CDSDX(29998,2),CUSDX(29998,2)
     *  /DUPSVZ/CDZ(29998,2,20),CUZ(29998,2,20),CUSUZ(29998,2)
     *         ,CDSUZ(29998,2),CDSDZ(29998,2),CUSDZ(29998,2)
     *  /GRTPSV/GTD(2,2,0:20),GTU(2,2,0:20),GRD(2,2,0:20),GRU(2,2,0:20)
     *  /WVNUM/WVN(29998),DWN,WNR(3),NWN
C
C *** CHECK FOR THE DIMENSION OF CD, CU, AND GTD-GRU ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(10,*)' FROM CUDPSV: NL.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
      NC=29998
      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.d0*PI)
        DQZ= WN/(2.d0*PI)
C
        IF((LS.EQ.1).AND.(SDEP.EQ.0.d0))THEN
C
          IF(NL.EQ.1)THEN
C
            DT=EE(3,1)*EE(4,2)-EE(3,2)*EE(4,1)
            E2(1,1)= EE(4,2)/DT
            E2(1,2)=-EE(3,2)/DT
            E2(2,1)=-EE(4,1)/DT
            E2(2,2)= EE(3,1)/DT
C
            CDSDX(K,1)=E2(1,1)*DQX
            CDSDX(K,2)=E2(2,1)*DQX
            CUSDX(K,1)=(0.d0,0.d0)
            CUSDX(K,2)=(0.d0,0.d0)
            CDSDZ(K,1)=E2(1,2)*DQZ
            CDSDZ(K,2)=E2(2,2)*DQZ
            CUSDZ(K,1)=(0.d0,0.d0)
            CUSDZ(K,2)=(0.d0,0.d0)
C
          ELSE
C
            CALL DDUPSV(OM,WN,THK(1),CMS,CMP,DU)
C
            E2(1,1)=EE(3,3)*DU(1,1)+EE(3,4)*DU(2,1)
            E2(1,2)=EE(3,3)*DU(1,2)+EE(3,4)*DU(2,2)
            E2(2,1)=EE(4,3)*DU(1,1)+EE(4,4)*DU(2,1)
            E2(2,2)=EE(4,3)*DU(1,2)+EE(4,4)*DU(2,2)
            E1(1,1)=EE(3,1)+E2(1,1)*GRD(1,1,1)+E2(1,2)*GRD(2,1,1)
            E1(1,2)=EE(3,2)+E2(1,1)*GRD(1,2,1)+E2(1,2)*GRD(2,2,1)
            E1(2,1)=EE(4,1)+E2(2,1)*GRD(1,1,1)+E2(2,2)*GRD(2,1,1)
            E1(2,2)=EE(4,2)+E2(2,1)*GRD(1,2,1)+E2(2,2)*GRD(2,2,1)
            DT=E1(1,1)*E1(2,2)-E1(1,2)*E1(2,1)
            E2(1,1)= E1(2,2)/DT
            E2(1,2)=-E1(1,2)/DT
            E2(2,1)=-E1(2,1)/DT
            E2(2,2)= E1(1,1)/DT
C
            CDSDX(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 VSc,CD,CU,VS2,WS2,CMS,SM,EE(2,2)
     *          ,CUSU,CDSU,CDSD,CUSD,DD,DU,VLX(2),OM,OM2
      COMMON /SOLY/SDEP,ODEP,LS,LO
     *       /STRC/DNS(20),VP(20),QP(20),FP(20)
     *                    ,VS(20),QS(20),FS(20),THK(20),NL
     *       /DUSH/CD(29998,20),CU(29998,20)
     *       ,CUSU(29998),CDSU(29998),CDSD(29998),CUSD(29998)
C
C *** CHECK FOR THE DIMENSION OF CD AND CU ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(10,*)' FROM RTSH: NL.GT.NC ',NL,NC
         STOP
      END IF
C
C *** COMPUTATION OF VLX ***
C
      OM2=OM*OM
      WN2=WN*WN
C
      QSe=FuncQS(OM,LO)
      VSc=VS(LO)*DCMPLX(1.D0,-0.5D0/QSe)
      VS2=VSc*VSc
      WS2=OM2/VS2
      CMS=CDSQRT(WN2-WS2)
      IF(dREAL(CMS).LT.0.d0)CMS=-CMS
C
      SM=DNS(LO)*VS2
      EE(1,1)=(1.d0,0.d0)
      EE(1,2)=(1.d0,0.d0)
      EE(2,1)=-SM*CMS
      EE(2,2)= SM*CMS
C
      IF(LO.EQ.NL)THEN
C        
        IF(LO.NE.LS)THEN
          DD=CDEXP(-CMS*ODEP)
          VLX(1)=EE(1,1)*DD*CD(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 CMP,CMS,EE(4,4),DD(2,2),DU(2,2),OM
     *    ,CDX,CUX,CUSUX,CDSUX,CDSDX,CUSDX,URX(2,2),E1(2,2)
     *    ,CDZ,CUZ,CUSUZ,CDSUZ,CDSDZ,CUSDZ,URZ(2,2),E2(2,2)
      COMMON /SOLY/SDEP,ODEP,LS,LO
     * /STRC/DNS(20),VP(20),QP(20),FP(20)
     *              ,VS(20),QS(20),FS(20),THK(20),NL
     * /DUPSVX/CDX(29998,2,20),CUX(29998,2,20),CUSUX(29998,2)
     *        ,CDSUX(29998,2),CDSDX(29998,2),CUSDX(29998,2)
     * /DUPSVZ/CDZ(29998,2,20),CUZ(29998,2,20),CUSUZ(29998,2)
     *        ,CDSUZ(29998,2),CDSDZ(29998,2),CUSDZ(29998,2)
C
C *** CHECK FOR THE DIMENSION OF CD AND CU ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(10,*)' FROM RTPSV: NL.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C *** COMPUTATION OF URX AND URZ ***
C
      CALL MTXPSV(OM,WN,LO,CMS,CMP,EE)
C
      IF(LO.EQ.NL)THEN
C        
        IF(LO.NE.LS)THEN
C
          CALL DDUPSV(OM,WN,ODEP,CMS,CMP,DD)
C
          E1(1,1)=EE(1,1)*DD(1,1)+EE(1,2)*DD(2,1)
          E1(1,2)=EE(1,1)*DD(1,2)+EE(1,2)*DD(2,2)
          E1(2,1)=EE(2,1)*DD(1,1)+EE(2,2)*DD(2,1)
          E1(2,2)=EE(2,1)*DD(1,2)+EE(2,2)*DD(2,2)
          URX(1,1)=E1(1,1)*CDX(K,1,LO)+E1(1,2)*CDX(K,2,LO)
          URX(2,1)=E1(2,1)*CDX(K,1,LO)+E1(2,2)*CDX(K,2,LO)
          URZ(1,1)=E1(1,1)*CDZ(K,1,LO)+E1(1,2)*CDZ(K,2,LO)
          URZ(2,1)=E1(2,1)*CDZ(K,1,LO)+E1(2,2)*CDZ(K,2,LO)
C
          E1(1,1)=EE(3,1)*DD(1,1)+EE(3,2)*DD(2,1)
          E1(1,2)=EE(3,1)*DD(1,2)+EE(3,2)*DD(2,2)
          E1(2,1)=EE(4,1)*DD(1,1)+EE(4,2)*DD(2,1)
          E1(2,2)=EE(4,1)*DD(1,2)+EE(4,2)*DD(2,2)
          URX(1,2)=E1(1,1)*CDX(K,1,LO)+E1(1,2)*CDX(K,2,LO)
          URX(2,2)=E1(2,1)*CDX(K,1,LO)+E1(2,2)*CDX(K,2,LO)
          URZ(1,2)=E1(1,1)*CDZ(K,1,LO)+E1(1,2)*CDZ(K,2,LO)
          URZ(2,2)=E1(2,1)*CDZ(K,1,LO)+E1(2,2)*CDZ(K,2,LO)
C
        ELSE
C
          IF(ODEP.GE.SDEP)THEN
C
            CALL DDUPSV(OM,WN,ODEP-SDEP,CMS,CMP,DD)
C
            E1(1,1)=EE(1,1)*DD(1,1)+EE(1,2)*DD(2,1)
            E1(1,2)=EE(1,1)*DD(1,2)+EE(1,2)*DD(2,2)
            E1(2,1)=EE(2,1)*DD(1,1)+EE(2,2)*DD(2,1)
            E1(2,2)=EE(2,1)*DD(1,2)+EE(2,2)*DD(2,2)
            URX(1,1)=E1(1,1)*CDSDX(K,1)+E1(1,2)*CDSDX(K,2)
            URX(2,1)=E1(2,1)*CDSDX(K,1)+E1(2,2)*CDSDX(K,2)
            URZ(1,1)=E1(1,1)*CDSDZ(K,1)+E1(1,2)*CDSDZ(K,2)
            URZ(2,1)=E1(2,1)*CDSDZ(K,1)+E1(2,2)*CDSDZ(K,2)
C
            E1(1,1)=EE(3,1)*DD(1,1)+EE(3,2)*DD(2,1)
            E1(1,2)=EE(3,1)*DD(1,2)+EE(3,2)*DD(2,2)
            E1(2,1)=EE(4,1)*DD(1,1)+EE(4,2)*DD(2,1)
            E1(2,2)=EE(4,1)*DD(1,2)+EE(4,2)*DD(2,2)
            URX(1,2)=E1(1,1)*CDSDX(K,1)+E1(1,2)*CDSDX(K,2)
            URX(2,2)=E1(2,1)*CDSDX(K,1)+E1(2,2)*CDSDX(K,2)
            URZ(1,2)=E1(1,1)*CDSDZ(K,1)+E1(1,2)*CDSDZ(K,2)
            URZ(2,2)=E1(2,1)*CDSDZ(K,1)+E1(2,2)*CDSDZ(K,2)
C
          ELSE
C
            CALL DDUPSV(OM,WN,ODEP,CMS,CMP,DD)
            CALL DDUPSV(OM,WN,SDEP-ODEP,CMS,CMP,DU)
C
            E1(1,1)=EE(1,1)*DD(1,1)+EE(1,2)*DD(2,1)
            E1(1,2)=EE(1,1)*DD(1,2)+EE(1,2)*DD(2,2)
            E1(2,1)=EE(2,1)*DD(1,1)+EE(2,2)*DD(2,1)
            E1(2,2)=EE(2,1)*DD(1,2)+EE(2,2)*DD(2,2)
            E2(1,1)=EE(1,3)*DU(1,1)+EE(1,4)*DU(2,1)
            E2(1,2)=EE(1,3)*DU(1,2)+EE(1,4)*DU(2,2)
            E2(2,1)=EE(2,3)*DU(1,1)+EE(2,4)*DU(2,1)
            E2(2,2)=EE(2,3)*DU(1,2)+EE(2,4)*DU(2,2)
C
            URX(1,1)=E1(1,1)*CDSUX(K,1)+E1(1,2)*CDSUX(K,2)
     *              +E2(1,1)*CUSUX(K,1)+E2(1,2)*CUSUX(K,2)
            URX(2,1)=E1(2,1)*CDSUX(K,1)+E1(2,2)*CDSUX(K,2)
     *              +E2(2,1)*CUSUX(K,1)+E2(2,2)*CUSUX(K,2)
            URZ(1,1)=E1(1,1)*CDSUZ(K,1)+E1(1,2)*CDSUZ(K,2)
     *              +E2(1,1)*CUSUZ(K,1)+E2(1,2)*CUSUZ(K,2)
            URZ(2,1)=E1(2,1)*CDSUZ(K,1)+E1(2,2)*CDSUZ(K,2)
     *              +E2(2,1)*CUSUZ(K,1)+E2(2,2)*CUSUZ(K,2)
C
            E1(1,1)=EE(3,1)*DD(1,1)+EE(3,2)*DD(2,1)
            E1(1,2)=EE(3,1)*DD(1,2)+EE(3,2)*DD(2,2)
            E1(2,1)=EE(4,1)*DD(1,1)+EE(4,2)*DD(2,1)
            E1(2,2)=EE(4,1)*DD(1,2)+EE(4,2)*DD(2,2)
            E2(1,1)=EE(3,3)*DU(1,1)+EE(3,4)*DU(2,1)
            E2(1,2)=EE(3,3)*DU(1,2)+EE(3,4)*DU(2,2)
            E2(2,1)=EE(4,3)*DU(1,1)+EE(4,4)*DU(2,1)
            E2(2,2)=EE(4,3)*DU(1,2)+EE(4,4)*DU(2,2)
C
            URX(1,2)=E1(1,1)*CDSUX(K,1)+E1(1,2)*CDSUX(K,2)
     *              +E2(1,1)*CUSUX(K,1)+E2(1,2)*CUSUX(K,2)
            URX(2,2)=E1(2,1)*CDSUX(K,1)+E1(2,2)*CDSUX(K,2)
     *              +E2(2,1)*CUSUX(K,1)+E2(2,2)*CUSUX(K,2)
            URZ(1,2)=E1(1,1)*CDSUZ(K,1)+E1(1,2)*CDSUZ(K,2)
     *              +E2(1,1)*CUSUZ(K,1)+E2(1,2)*CUSUZ(K,2)
            URZ(2,2)=E1(2,1)*CDSUZ(K,1)+E1(2,2)*CDSUZ(K,2)
     *              +E2(2,1)*CUSUZ(K,1)+E2(2,2)*CUSUZ(K,2)
C
          END IF
C
        END IF
C
      ELSE
C        
        IF(LO.NE.LS)THEN
C
            CALL DDUPSV(OM,WN,ODEP,CMS,CMP,DD)
            CALL DDUPSV(OM,WN,THK(LO)-ODEP,CMS,CMP,DU)
C
            E1(1,1)=EE(1,1)*DD(1,1)+EE(1,2)*DD(2,1)
            E1(1,2)=EE(1,1)*DD(1,2)+EE(1,2)*DD(2,2)
            E1(2,1)=EE(2,1)*DD(1,1)+EE(2,2)*DD(2,1)
            E1(2,2)=EE(2,1)*DD(1,2)+EE(2,2)*DD(2,2)
            E2(1,1)=EE(1,3)*DU(1,1)+EE(1,4)*DU(2,1)
            E2(1,2)=EE(1,3)*DU(1,2)+EE(1,4)*DU(2,2)
            E2(2,1)=EE(2,3)*DU(1,1)+EE(2,4)*DU(2,1)
            E2(2,2)=EE(2,3)*DU(1,2)+EE(2,4)*DU(2,2)
C
            URX(1,1)=E1(1,1)*CDX(K,1,LO)+E1(1,2)*CDX(K,2,LO)
     *              +E2(1,1)*CUX(K,1,LO)+E2(1,2)*CUX(K,2,LO)
            URX(2,1)=E1(2,1)*CDX(K,1,LO)+E1(2,2)*CDX(K,2,LO)
     *              +E2(2,1)*CUX(K,1,LO)+E2(2,2)*CUX(K,2,LO)
            URZ(1,1)=E1(1,1)*CDZ(K,1,LO)+E1(1,2)*CDZ(K,2,LO)
     *              +E2(1,1)*CUZ(K,1,LO)+E2(1,2)*CUZ(K,2,LO)
            URZ(2,1)=E1(2,1)*CDZ(K,1,LO)+E1(2,2)*CDZ(K,2,LO)
     *              +E2(2,1)*CUZ(K,1,LO)+E2(2,2)*CUZ(K,2,LO)
C
            E1(1,1)=EE(3,1)*DD(1,1)+EE(3,2)*DD(2,1)
            E1(1,2)=EE(3,1)*DD(1,2)+EE(3,2)*DD(2,2)
            E1(2,1)=EE(4,1)*DD(1,1)+EE(4,2)*DD(2,1)
            E1(2,2)=EE(4,1)*DD(1,2)+EE(4,2)*DD(2,2)
            E2(1,1)=EE(3,3)*DU(1,1)+EE(3,4)*DU(2,1)
            E2(1,2)=EE(3,3)*DU(1,2)+EE(3,4)*DU(2,2)
            E2(2,1)=EE(4,3)*DU(1,1)+EE(4,4)*DU(2,1)
            E2(2,2)=EE(4,3)*DU(1,2)+EE(4,4)*DU(2,2)
C
            URX(1,2)=E1(1,1)*CDX(K,1,LO)+E1(1,2)*CDX(K,2,LO)
     *              +E2(1,1)*CUX(K,1,LO)+E2(1,2)*CUX(K,2,LO)
            URX(2,2)=E1(2,1)*CDX(K,1,LO)+E1(2,2)*CDX(K,2,LO)
     *              +E2(2,1)*CUX(K,1,LO)+E2(2,2)*CUX(K,2,LO)
            URZ(1,2)=E1(1,1)*CDZ(K,1,LO)+E1(1,2)*CDZ(K,2,LO)
     *              +E2(1,1)*CUZ(K,1,LO)+E2(1,2)*CUZ(K,2,LO)
            URZ(2,2)=E1(2,1)*CDZ(K,1,LO)+E1(2,2)*CDZ(K,2,LO)
     *              +E2(2,1)*CUZ(K,1,LO)+E2(2,2)*CUZ(K,2,LO)
C
        ELSE
C
          IF(ODEP.GE.SDEP)THEN
C
            CALL DDUPSV(OM,WN,ODEP-SDEP,CMS,CMP,DD)
            CALL DDUPSV(OM,WN,THK(LO)-ODEP,CMS,CMP,DU)
C
            E1(1,1)=EE(1,1)*DD(1,1)+EE(1,2)*DD(2,1)
            E1(1,2)=EE(1,1)*DD(1,2)+EE(1,2)*DD(2,2)
            E1(2,1)=EE(2,1)*DD(1,1)+EE(2,2)*DD(2,1)
            E1(2,2)=EE(2,1)*DD(1,2)+EE(2,2)*DD(2,2)
            E2(1,1)=EE(1,3)*DU(1,1)+EE(1,4)*DU(2,1)
            E2(1,2)=EE(1,3)*DU(1,2)+EE(1,4)*DU(2,2)
            E2(2,1)=EE(2,3)*DU(1,1)+EE(2,4)*DU(2,1)
            E2(2,2)=EE(2,3)*DU(1,2)+EE(2,4)*DU(2,2)
C
            URX(1,1)=E1(1,1)*CDSDX(K,1)+E1(1,2)*CDSDX(K,2)
     *              +E2(1,1)*CUSDX(K,1)+E2(1,2)*CUSDX(K,2)
            URX(2,1)=E1(2,1)*CDSDX(K,1)+E1(2,2)*CDSDX(K,2)
     *              +E2(2,1)*CUSDX(K,1)+E2(2,2)*CUSDX(K,2)
            URZ(1,1)=E1(1,1)*CDSDZ(K,1)+E1(1,2)*CDSDZ(K,2)
     *              +E2(1,1)*CUSDZ(K,1)+E2(1,2)*CUSDZ(K,2)
            URZ(2,1)=E1(2,1)*CDSDZ(K,1)+E1(2,2)*CDSDZ(K,2)
     *              +E2(2,1)*CUSDZ(K,1)+E2(2,2)*CUSDZ(K,2)
C
            E1(1,1)=EE(3,1)*DD(1,1)+EE(3,2)*DD(2,1)
            E1(1,2)=EE(3,1)*DD(1,2)+EE(3,2)*DD(2,2)
            E1(2,1)=EE(4,1)*DD(1,1)+EE(4,2)*DD(2,1)
            E1(2,2)=EE(4,1)*DD(1,2)+EE(4,2)*DD(2,2)
            E2(1,1)=EE(3,3)*DU(1,1)+EE(3,4)*DU(2,1)
            E2(1,2)=EE(3,3)*DU(1,2)+EE(3,4)*DU(2,2)
            E2(2,1)=EE(4,3)*DU(1,1)+EE(4,4)*DU(2,1)
            E2(2,2)=EE(4,3)*DU(1,2)+EE(4,4)*DU(2,2)
C
            URX(1,2)=E1(1,1)*CDSDX(K,1)+E1(1,2)*CDSDX(K,2)
     *              +E2(1,1)*CUSDX(K,1)+E2(1,2)*CUSDX(K,2)
            URX(2,2)=E1(2,1)*CDSDX(K,1)+E1(2,2)*CDSDX(K,2)
     *              +E2(2,1)*CUSDX(K,1)+E2(2,2)*CUSDX(K,2)
            URZ(1,2)=E1(1,1)*CDSDZ(K,1)+E1(1,2)*CDSDZ(K,2)
     *              +E2(1,1)*CUSDZ(K,1)+E2(1,2)*CUSDZ(K,2)
            URZ(2,2)=E1(2,1)*CDSDZ(K,1)+E1(2,2)*CDSDZ(K,2)
     *              +E2(2,1)*CUSDZ(K,1)+E2(2,2)*CUSDZ(K,2)
C
          ELSE
C
            CALL DDUPSV(OM,WN,ODEP,CMS,CMP,DD)
            CALL DDUPSV(OM,WN,SDEP-ODEP,CMS,CMP,DU)
C
            E1(1,1)=EE(1,1)*DD(1,1)+EE(1,2)*DD(2,1)
            E1(1,2)=EE(1,1)*DD(1,2)+EE(1,2)*DD(2,2)
            E1(2,1)=EE(2,1)*DD(1,1)+EE(2,2)*DD(2,1)
            E1(2,2)=EE(2,1)*DD(1,2)+EE(2,2)*DD(2,2)
            E2(1,1)=EE(1,3)*DU(1,1)+EE(1,4)*DU(2,1)
            E2(1,2)=EE(1,3)*DU(1,2)+EE(1,4)*DU(2,2)
            E2(2,1)=EE(2,3)*DU(1,1)+EE(2,4)*DU(2,1)
            E2(2,2)=EE(2,3)*DU(1,2)+EE(2,4)*DU(2,2)
C
            URX(1,1)=E1(1,1)*CDSUX(K,1)+E1(1,2)*CDSUX(K,2)
     *              +E2(1,1)*CUSUX(K,1)+E2(1,2)*CUSUX(K,2)
            URX(2,1)=E1(2,1)*CDSUX(K,1)+E1(2,2)*CDSUX(K,2)
     *              +E2(2,1)*CUSUX(K,1)+E2(2,2)*CUSUX(K,2)
            URZ(1,1)=E1(1,1)*CDSUZ(K,1)+E1(1,2)*CDSUZ(K,2)
     *              +E2(1,1)*CUSUZ(K,1)+E2(1,2)*CUSUZ(K,2)
            URZ(2,1)=E1(2,1)*CDSUZ(K,1)+E1(2,2)*CDSUZ(K,2)
     *              +E2(2,1)*CUSUZ(K,1)+E2(2,2)*CUSUZ(K,2)
C
            E1(1,1)=EE(3,1)*DD(1,1)+EE(3,2)*DD(2,1)
            E1(1,2)=EE(3,1)*DD(1,2)+EE(3,2)*DD(2,2)
            E1(2,1)=EE(4,1)*DD(1,1)+EE(4,2)*DD(2,1)
            E1(2,2)=EE(4,1)*DD(1,2)+EE(4,2)*DD(2,2)
            E2(1,1)=EE(3,3)*DU(1,1)+EE(3,4)*DU(2,1)
            E2(1,2)=EE(3,3)*DU(1,2)+EE(3,4)*DU(2,2)
            E2(2,1)=EE(4,3)*DU(1,1)+EE(4,4)*DU(2,1)
            E2(2,2)=EE(4,3)*DU(1,2)+EE(4,4)*DU(2,2)
C
            URX(1,2)=E1(1,1)*CDSUX(K,1)+E1(1,2)*CDSUX(K,2)
     *              +E2(1,1)*CUSUX(K,1)+E2(1,2)*CUSUX(K,2)
            URX(2,2)=E1(2,1)*CDSUX(K,1)+E1(2,2)*CDSUX(K,2)
     *              +E2(2,1)*CUSUX(K,1)+E2(2,2)*CUSUX(K,2)
            URZ(1,2)=E1(1,1)*CDSUZ(K,1)+E1(1,2)*CDSUZ(K,2)
     *              +E2(1,1)*CUSUZ(K,1)+E2(1,2)*CUSUZ(K,2)
            URZ(2,2)=E1(2,1)*CDSUZ(K,1)+E1(2,2)*CDSUZ(K,2)
     *              +E2(2,1)*CUSUZ(K,1)+E2(2,2)*CUSUZ(K,2)
C
          END IF
C
        END IF
C
      END IF
C
      RETURN
      END
C
C
C *** 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.d0))
      CUSX=CUPX*(X2-1.d0)
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.d0))
      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.d0))
      CUSX=CUPX*(X2-1.d0)
C
C   * P-SV (DOWN-GOING REFLECTED FROM UPPER BOUNDARY) *
C
      IF(LS.EQ.1)THEN
        BB1(1,3)=-(X2-2.d0)
        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.d0))
      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.d0))
      CDSX=CDPX*(X2-1.d0)
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.d0)
        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.d0))
      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.d0))
      CDSX=CDPX*(X2-1.d0)
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.d0))
      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.5d0/(X1+A2*X2+A*(X1*X2+1.d0))
C
      BB(1,1)=2.d0*(X1+1.d0)*A*(A+X1)
      BB(2,1)=(X1+1.d0)*A*(X1*X2-3.d0*X1+2.d0-A*(X1*X2-3.d0*X2+2.d0))
      BB(3,1)=2.d0*(X1-2.d0)*(1.d0-A)*(1.d0+A*X2)
      BB(4,1)=A*(X2-1.d0)*(X1*X1-3.d0*X1+4.d0)
     *         -(X1-1.d0)*(4.d0+A2*(X1-3.d0)*X2)
C
      BB(1,2)=(0.d0,0.d0)
      BB(2,2)=2.d0*(X1+1.d0)*A*(1.d0+A*X2)
      BB(3,2)=4.d0*(A-1.d0)*(1.d0+A*X2)
      BB(4,2)=-BB(3,1)
C
      BB(1,3)=2.d0*(X2-2.d0)*(A-1.d0)*(A+X1)
      BB(2,3)=A*(X1-1.d0)*(X2*X2-3.d0*X2+4.d0)
     *         -(X2-1.d0)*(4.d0*A2+(X2-3.d0)*X1)
      BB(3,3)=2.d0*(X2+1.d0)*(1.d0+A*X2)
      BB(4,3)=(X2+1.d0)*(A*(X1*X2-3.d0*X2+2.d0)-X1*X2+3.d0*X1-2.d0)
C
      BB(1,4)=4.d0*(1.d0-A)*(A+X1)
      BB(2,4)=-BB(1,3)
      BB(3,4)=(0.d0,0.d0)
      BB(4,4)=2.d0*(X2+1.d0)*(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.d0
      B=X1-3.d0
      V1X=CPT+CST
      V2X=A*CPT-CST
      V3X=WN*SM1*(-B*CPT+2.d0*CST)
      V4X=WN*SM1*( A*CPT-2.d0*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.d0*CST)
      V4Z=WN*SM1*( A*CPT-2.d0*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.d0
      B=X2-3.d0
      V1X=CPR1+CSR1+CPD+CSD+CPR2+CSR2
      V2X=-A*CPR1+CSR1+A*(CPD+CPR2)-CSD-CSR2
      V3X=WN*SM2*(B*CPR1-2.d0*CSR1-B*(CPD+CPR2)+2.d0*(CSD+CSR2))
      V4X=WN*SM2*(A*CPR1-2.d0*CSR1+A*(CPD+CPR2)-2.d0*(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.d0*CSR1-B*(CPD+CPR2)+2.d0*(CSD+CSR2))
      V4Z=WN*SM2*(A*CPR1-2.d0*CSR1+A*(CPD+CPR2)-2.d0*(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.d0
      B=X2-3.d0
      V1X=CPD+CSD+CPR1+CSR1+CPR2+CSR2
      V2X=-A*(CPD+CPR1)+CSD+CSR1+A*CPR2-CSR2
      V3X=WN*SM2*(B*(CPD+CPR1)-2.d0*(CSD+CSR1)-B*CPR2+2.d0*CSR2)
      V4X=WN*SM2*(A*(CPD+CPR1)-2.d0*(CSD+CSR1)+A*CPR2-2.d0*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.d0*(CSD+CSR1)-B*CPR2+2.d0*CSR2)
      V4Z=WN*SM2*(A*(CPD+CPR1)-2.d0*(CSD+CSR1)+A*CPR2-2.d0*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.d0
      B=X3-3.d0
      V1X=CPT+CST
      V2X=-A*CPT+CST
      V3X=WN*SM3*(B*CPT-2.d0*CST)
      V4X=WN*SM3*(A*CPT-2.d0*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.d0*CST)
      V4Z=WN*SM3*(A*CPT-2.d0*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.d0)*DCUPT-2.d0*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.d0)*DCUPT-DCUST
C
      DCUPT=CUP1X*S10-DC1*S11
      DCUST=CUS1X*S10-DC2*S11+DC3*S12
      D10V2=(X1-1.d0)*DCUPT-DCUST
      D1MV3=SM1*(-(X1-3.d0)*DCUPT+2.d0*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.d0)*DCUPT+2.d0*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.d0)*DCUPT-2.d0*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.d0)*DCUPT-DCUST
C
      DCUPT=CUP1Z*S11-DC1*S12
      DCUST=CUS1Z*S11-DC2*S12+DC3*S13
      D11V2=(X1-1.d0)*DCUPT-DCUST
      D10V3=SM1*(-(X1-3.d0)*DCUPT+2.d0*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.d0)*DCUPD-2.d0*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.d0)*DCUPD-DCUSD
C
      DCUPD=CUPX*S10
      DCUSD=CUSX*S10-DC1*S11
      D10V2=(X2-1.d0)*DCUPD-DCUSD
      D1MV3=SM2*(-(X2-3.d0)*DCUPD+2.d0*DCUSD)
C
      S1=S01-S10/R
      S2=S02-S11/R
      DCUPD=CUPX*S1
      DCUSD=CUSX*S1-DC1*S2
      DD0V3=SM2*(-(X2-3.d0)*DCUPD+2.d0*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.d0)*DCDPR-2.d0*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.d0)*DCDPR+DCDSR
C
      DCDPR=CDP1X*S101-DC1*S111
      DCDSR=CDS1X*S101-DC2*S111+DC3*S121
      D10V2=-(X2-1.d0)*DCDPR+DCDSR
      D1MV3=SM2*((X2-3.d0)*DCDPR-2.d0*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.d0)*DCDPR-2.d0*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.d0)*DCUPR-2.d0*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.d0)*DCUPR-DCUSR
C
      DCUPR=CUP1X*S102-DC1*S112
      DCUSR=CUS1X*S102-DC2*S112+DC3*S122
      D10V2=(X2-1.d0)*DCUPR-DCUSR
      D1MV3=SM2*(-(X2-3.d0)*DCUPR+2.d0*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.d0)*DCUPR+2.d0*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.d0)*DCUPD-2.d0*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.d0)*DCUPD-DCUSD
C
      DCUPD=CUPZ*S11
      DCUSD=CUSZ*S11-DC1*S12
      D11V2=(X2-1.d0)*DCUPD-DCUSD
      D10V3=SM2*(-(X2-3.d0)*DCUPD+2.d0*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.d0)*DCDPR-2.d0*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.d0)*DCDPR+DCDSR
C
      DCDPR=CDP1Z*S111-DC1*S121
      DCDSR=CDS1Z*S111-DC2*S121+DC3*S131
      D11V2=-(X2-1.d0)*DCDPR+DCDSR
      D10V3=SM2*((X2-3.d0)*DCDPR-2.d0*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.d0)*DCUPR-2.d0*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.d0)*DCUPR-DCUSR
C
      DCUPR=CUP1Z*S112-DC1*S122
      DCUSR=CUS1Z*S112-DC2*S122+DC3*S132
      D11V2=(X2-1.d0)*DCUPR-DCUSR
      D10V3=SM2*(-(X2-3.d0)*DCUPR+2.d0*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.d0)*DCDPD-2.d0*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.d0)*DCDPD+DCDSD
C
      DCDPD=CDPX*S10
      DCDSD=CDSX*S10-DC1*S11
      D10V2=-(X2-1.d0)*DCDPD+DCDSD
      D1MV3=SM2*((X2-3.d0)*DCDPD-2.d0*DCDSD)
C
      S1=S01-S10/R
      S2=S02-S11/R
      DCDPD=CDPX*S1
      DCDSD=CDSX*S1-DC1*S2
      DD0V3=SM2*((X2-3.d0)*DCDPD-2.d0*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.d0)*DCUPR-2.d0*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.d0)*DCUPR-DCUSR
C
      DCUPR=CUP1X*S102-DC1*S112
      DCUSR=CUS1X*S102-DC2*S112+DC3*S122
      D10V2=(X2-1.d0)*DCUPR-DCUSR
      D1MV3=SM2*(-(X2-3.d0)*DCUPR+2.d0*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.d0)*DCUPR+2.d0*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.d0)*DCDPR-2.d0*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.d0)*DCDPR+DCDSR
C
      DCDPR=CDP1X*S101-DC1*S111
      DCDSR=CDS1X*S101-DC2*S111+DC3*S121
      D10V2=-(X2-1.d0)*DCDPR+DCDSR
      D1MV3=SM2*((X2-3.d0)*DCDPR-2.d0*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.d0)*DCDPR-2.d0*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.d0)*DCDPD-2.d0*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.d0)*DCDPD+DCDSD
C
      DCDPD=CDPZ*S11
      DCDSD=CDSZ*S11-DC1*S12
      D11V2=-(X2-1.d0)*DCDPD+DCDSD
      D10V3=SM2*((X2-3.d0)*DCDPD-2.d0*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.d0)*DCUPR-2.d0*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.d0)*DCUPR-DCUSR
C
      DCUPR=CUP1Z*S112-DC1*S122
      DCUSR=CUS1Z*S112-DC2*S122+DC3*S132
      D11V2=(X2-1.d0)*DCUPR-DCUSR
      D10V3=SM2*(-(X2-3.d0)*DCUPR+2.d0*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.d0)*DCDPR-2.d0*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.d0)*DCDPR+DCDSR
C
      DCDPR=CDP1Z*S111-DC1*S121
      DCDSR=CDS1Z*S111-DC2*S121+DC3*S131
      D11V2=-(X2-1.d0)*DCDPR+DCDSR
      D10V3=SM2*((X2-3.d0)*DCDPR-2.d0*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.d0)*DCDPT-2.d0*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.d0)*DCDPT+DCDST
C
      DCDPT=CDP1X*S10-DC1*S11
      DCDST=CDS1X*S10-DC2*S11+DC3*S12
      D10V2=-(X3-1.d0)*DCDPT+DCDST
      D1MV3=SM3*((X3-3.d0)*DCDPT-2.d0*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.d0)*DCDPT-2.d0*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.d0)*DCDPT-2.d0*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.d0)*DCDPT+DCDST
C
      DCDPT=CDP1Z*S11-DC1*S12
      DCDST=CDS1Z*S11-DC2*S12+DC3*S13
      D11V2=-(X3-1.d0)*DCDPT+DCDST
      D10V3=SM3*((X3-3.d0)*DCDPT-2.d0*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 GTD,GTU,GRD,GRU,TD,TU,RD,RU,OM
      COMMON /SOLY/SDEP,ODEP,LS,LO
     *    /STRC/DNS(20),VP(20),QP(20),FP(20)
     *                 ,VS(20),QS(20),FS(20),THK(20),NL
      COMMON /GRTSH/GTD(0:20),GTU(0:20),GRD(0:20),GRU(0:20)
     *    /MRTSH/TD(0:20),TU(0:20),RD(0:20),RU(0:20)
C
C *** CHECK FOR THE DIMENSION OF GTD, GTU, GRD, AND GRU ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(10,*)' FROM GNRTSH: NL.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C *** MODIFIED R/T MATRIX ***
C
      CALL MDRTSH(OM,WN)
C
C *** GENERALIZED R/T MATRIX ***
C
      GRU(0)=RU(0)
      DO 100 I=1,LS-1
        GTU(I)=TU(I)/(1.D0-RD(I)*GRU(I-1))
        GRU(I)=RU(I)+TD(I)*GRU(I-1)*GTU(I)
  100 CONTINUE
C
      GTD(NL-1)=TD(NL-1)
      GRD(NL-1)=RD(NL-1)
      DO 200 I=NL-2,LS,-1
        GTD(I)=TD(I)/(1.D0-RU(I)*GRD(I+1))
        GRD(I)=RD(I)+TU(I)*GRD(I+1)*GTD(I)
  200 CONTINUE
C
      RETURN
      END
C
C
C
C *** SH MODIFIED R/T MATIX ***
C
      SUBROUTINE MDRTSH(OM,WN)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VSc,VS2,WS2,CMS,SM,EA(2,2),EB(2,2)
     *      ,TD,TU,RD,RU,DDA,DUB,EAB,OM,OM2
      COMMON /SOLY/SDEP,ODEP,LS,LO
     *    /STRC/DNS(20),VP(20),QP(20),FP(20)
     *                 ,VS(20),QS(20),FS(20),THK(20),NL
      COMMON /MRTSH/TD(0:20),TU(0:20),RD(0:20),RU(0:20)
C
C *** CHECK FOR THE DIMENSION OF TD, TU, RD, AND RU ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(10,*)' FROM MDRTSH: NL.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C *** MODIFIED R/T MATRIX ***
C
      WN2=WN*WN
      OM2=OM*OM
C
      QSe=FuncQS(OM,1)
      VSc=VS(1)*DCMPLX(1.D0,-0.5D0/QSe)
      VS2=VSc*VSc
      WS2=OM2/VS2
      CMS=CDSQRT(WN2-WS2)
      IF(DREAL(CMS).LT.0.D0)CMS=-CMS
C
      SM=DNS(1)*VS2
      EA(1,1)=(1.D0,0.D0)
      EA(1,2)=(1.D0,0.D0)
      EA(2,1)=-SM*CMS
      EA(2,2)= SM*CMS
C
      IF(NL.EQ.1)THEN
        DUB=CDEXP(-CMS*SDEP)
c        RU(0)=-EA(2,2)*DUB/EA(2,1)
        RU(0)=DUB
        RETURN
      END IF
C
      IF(LS.EQ.1)THEN
        DUB=CDEXP(-CMS*SDEP)
        DDA=CDEXP(-CMS*(THK(1)-SDEP))
      ELSE
        DUB=CDEXP(-CMS*THK(1))
        DDA=DUB
      END IF
C
c      RU(0)=-EA(2,2)*DUB/EA(2,1)
      RU(0)=DUB
C
      DO 100 I=1, NL-2
C
        QSe=FuncQS(OM,I+1)
        VSc=VS(I+1)*DCMPLX(1.D0,-0.5D0/QSe)
        VS2=VSc*VSc
        WS2=OM2/VS2
        CMS=CDSQRT(WN2-WS2)
        IF(DREAL(CMS).LT.0.D0)CMS=-CMS
C
        SM=DNS(I+1)*VS2
        EB(1,1)=(1.D0,0.D0)
        EB(1,2)=(1.D0,0.D0)
        EB(2,1)=-SM*CMS
        EB(2,2)= SM*CMS
C
        IF(LS.EQ.(I+1))THEN
          DUB=CDEXP(-CMS*SDEP)
        ELSE
          DUB=CDEXP(-CMS*THK(I+1))
        END IF
C
        EAB=EB(2,1)*EA(1,2)-EB(1,1)*EA(2,2) 
        TD(I)=(EA(1,2)*EA(2,1)-EA(1,1)*EA(2,2))*DDA/EAB
        RU(I)=(EB(1,2)*EA(2,2)-EB(2,2)*EA(1,2))*DUB/EAB
        RD(I)=(EB(1,1)*EA(2,1)-EA(1,1)*EB(2,1))*DDA/EAB
        TU(I)=(EB(1,2)*EB(2,1)-EB(2,2)*EB(1,1))*DUB/EAB
C
        EA(1,1)=EB(1,1)
        EA(1,2)=EB(1,2)
        EA(2,1)=EB(2,1)
        EA(2,2)=EB(2,2)
C
        IF(LS.EQ.(I+1))THEN
          DDA=CDEXP(-CMS*(THK(I+1)-SDEP))
        ELSE          
          DDA=DUB
        END IF
C
  100 CONTINUE
C        
      QSe=FuncQS(OM,NL)
      VSc=VS(NL)*DCMPLX(1.D0,-0.5D0/QSe)
      VS2=VSc*VSc
      WS2=OM2/VS2
      CMS=CDSQRT(WN2-WS2)
      IF(DREAL(CMS).LT.0.D0)CMS=-CMS
C
      SM=DNS(NL)*VS2
      EB(1,1)=(1.D0,0.D0)
      EB(1,2)=(1.D0,0.D0)
      EB(2,1)=-SM*CMS
      EB(2,2)= SM*CMS
C
      EAB=EB(2,1)*EA(1,2)-EB(1,1)*EA(2,2) 
C
      IF(LS.EQ.NL)THEN
C
        DUB=CDEXP(-CMS*SDEP)
        TD(NL-1)=(EA(1,2)*EA(2,1)-EA(1,1)*EA(2,2))*DDA/EAB
        RU(NL-1)=(EB(1,2)*EA(2,2)-EB(2,2)*EA(1,2))*DUB/EAB
        RD(NL-1)=(EB(1,1)*EA(2,1)-EA(1,1)*EB(2,1))*DDA/EAB
        TU(NL-1)=(EB(1,2)*EB(2,1)-EB(2,2)*EB(1,1))*DUB/EAB
C
      ELSE
C
        TD(NL-1)=(EA(1,2)*EA(2,1)-EA(1,1)*EA(2,2))*DDA/EAB
        RD(NL-1)=(EB(1,1)*EA(2,1)-EA(1,1)*EB(2,1))*DDA/EAB
C
      END IF
C
      RETURN
      END
C
C
C *** P-SV GENERALIZED R/T MATIX ***
C
      SUBROUTINE GNRTPSV(OM,WN)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 GTD,GTU,GRD,GRU,TD,TU,RD,RU,E1(2,2),E2(2,2),DT,OM
      COMMON /SOLY/SDEP,ODEP,LS,LO
     *    /STRC/DNS(20),VP(20),QP(20),FP(20)
     *                 ,VS(20),QS(20),FS(20),THK(20),NL
      COMMON /MRTPSV/TD(2,2,0:20),TU(2,2,0:20),RD(2,2,0:20),RU(2,2,0:20)
     *  /GRTPSV/GTD(2,2,0:20),GTU(2,2,0:20),GRD(2,2,0:20),GRU(2,2,0:20)
C
C *** CHECK FOR THE DIMENSION OF GTD, GTU, GRD, AND GRU ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(10,*)' FROM GNRTPSV: NL.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C *** MODIFIED R/T MATRIX ***
C
      CALL MDRTPSV(OM,WN)
C
C *** GENERALIZED R/T MATRIX ***
C
      GRU(1,1,0)=RU(1,1,0)
      GRU(1,2,0)=RU(1,2,0)
      GRU(2,1,0)=RU(2,1,0)
      GRU(2,2,0)=RU(2,2,0)
C
      DO 100 I=1,LS-1
C
        E1(1,1)=1.D0-RD(1,1,I)*GRU(1,1,I-1)-RD(1,2,I)*GRU(2,1,I-1)
        E1(1,2)=    -RD(1,1,I)*GRU(1,2,I-1)-RD(1,2,I)*GRU(2,2,I-1)
        E1(2,1)=    -RD(2,1,I)*GRU(1,1,I-1)-RD(2,2,I)*GRU(2,1,I-1)
        E1(2,2)=1.D0-RD(2,1,I)*GRU(1,2,I-1)-RD(2,2,I)*GRU(2,2,I-1)
        DT=E1(1,1)*E1(2,2)-E1(1,2)*E1(2,1)
        E2(1,1)= E1(2,2)/DT
        E2(1,2)=-E1(1,2)/DT
        E2(2,1)=-E1(2,1)/DT
        E2(2,2)= E1(1,1)/DT
        GTU(1,1,I)=E2(1,1)*TU(1,1,I)+E2(1,2)*TU(2,1,I)
        GTU(1,2,I)=E2(1,1)*TU(1,2,I)+E2(1,2)*TU(2,2,I)
        GTU(2,1,I)=E2(2,1)*TU(1,1,I)+E2(2,2)*TU(2,1,I)
        GTU(2,2,I)=E2(2,1)*TU(1,2,I)+E2(2,2)*TU(2,2,I)
C
        E1(1,1)=TD(1,1,I)*GRU(1,1,I-1)+TD(1,2,I)*GRU(2,1,I-1)
        E1(1,2)=TD(1,1,I)*GRU(1,2,I-1)+TD(1,2,I)*GRU(2,2,I-1)
        E1(2,1)=TD(2,1,I)*GRU(1,1,I-1)+TD(2,2,I)*GRU(2,1,I-1)
        E1(2,2)=TD(2,1,I)*GRU(1,2,I-1)+TD(2,2,I)*GRU(2,2,I-1)
        E2(1,1)=E1(1,1)*GTU(1,1,I)+E1(1,2)*GTU(2,1,I)
        E2(1,2)=E1(1,1)*GTU(1,2,I)+E1(1,2)*GTU(2,2,I)
        E2(2,1)=E1(2,1)*GTU(1,1,I)+E1(2,2)*GTU(2,1,I)
        E2(2,2)=E1(2,1)*GTU(1,2,I)+E1(2,2)*GTU(2,2,I)
        GRU(1,1,I)=RU(1,1,I)+E2(1,1)
        GRU(1,2,I)=RU(1,2,I)+E2(1,2)
        GRU(2,1,I)=RU(2,1,I)+E2(2,1)
        GRU(2,2,I)=RU(2,2,I)+E2(2,2)
C
  100 CONTINUE
C
      GTD(1,1,NL-1)=TD(1,1,NL-1)
      GTD(1,2,NL-1)=TD(1,2,NL-1)
      GTD(2,1,NL-1)=TD(2,1,NL-1)
      GTD(2,2,NL-1)=TD(2,2,NL-1)
C
      GRD(1,1,NL-1)=RD(1,1,NL-1)
      GRD(1,2,NL-1)=RD(1,2,NL-1)
      GRD(2,1,NL-1)=RD(2,1,NL-1)
      GRD(2,2,NL-1)=RD(2,2,NL-1)
C
      DO 200 I=NL-2,LS,-1
C
        E1(1,1)=1.D0-RU(1,1,I)*GRD(1,1,I+1)-RU(1,2,I)*GRD(2,1,I+1)
        E1(1,2)=    -RU(1,1,I)*GRD(1,2,I+1)-RU(1,2,I)*GRD(2,2,I+1)
        E1(2,1)=    -RU(2,1,I)*GRD(1,1,I+1)-RU(2,2,I)*GRD(2,1,I+1)
        E1(2,2)=1.D0-RU(2,1,I)*GRD(1,2,I+1)-RU(2,2,I)*GRD(2,2,I+1)
        DT=E1(1,1)*E1(2,2)-E1(1,2)*E1(2,1)
        E2(1,1)= E1(2,2)/DT
        E2(1,2)=-E1(1,2)/DT
        E2(2,1)=-E1(2,1)/DT
        E2(2,2)= E1(1,1)/DT
        GTD(1,1,I)=E2(1,1)*TD(1,1,I)+E2(1,2)*TD(2,1,I)
        GTD(1,2,I)=E2(1,1)*TD(1,2,I)+E2(1,2)*TD(2,2,I)
        GTD(2,1,I)=E2(2,1)*TD(1,1,I)+E2(2,2)*TD(2,1,I)
        GTD(2,2,I)=E2(2,1)*TD(1,2,I)+E2(2,2)*TD(2,2,I)
C
        E1(1,1)=TU(1,1,I)*GRD(1,1,I+1)+TU(1,2,I)*GRD(2,1,I+1)
        E1(1,2)=TU(1,1,I)*GRD(1,2,I+1)+TU(1,2,I)*GRD(2,2,I+1)
        E1(2,1)=TU(2,1,I)*GRD(1,1,I+1)+TU(2,2,I)*GRD(2,1,I+1)
        E1(2,2)=TU(2,1,I)*GRD(1,2,I+1)+TU(2,2,I)*GRD(2,2,I+1)
        E2(1,1)=E1(1,1)*GTD(1,1,I)+E1(1,2)*GTD(2,1,I)
        E2(1,2)=E1(1,1)*GTD(1,2,I)+E1(1,2)*GTD(2,2,I)
        E2(2,1)=E1(2,1)*GTD(1,1,I)+E1(2,2)*GTD(2,1,I)
        E2(2,2)=E1(2,1)*GTD(1,2,I)+E1(2,2)*GTD(2,2,I)
        GRD(1,1,I)=RD(1,1,I)+E2(1,1)
        GRD(1,2,I)=RD(1,2,I)+E2(1,2)
        GRD(2,1,I)=RD(2,1,I)+E2(2,1)
        GRD(2,2,I)=RD(2,2,I)+E2(2,2)
C
  200 CONTINUE
C
      RETURN
      END
C
C *** P-SV MODIFIED R/T MATIX ***
C
      SUBROUTINE MDRTPSV(OM,WN)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 CMS,CMP,EA(4,4),EB(4,4),TD,TU,RD,RU
     *       ,DDA(2,2),DUB(2,2),DET,OM
      COMMON /SOLY/SDEP,ODEP,LS,LO
     *    /STRC/DNS(20),VP(20),QP(20),FP(20)
     *                 ,VS(20),QS(20),FS(20),THK(20),NL
     *    /MRTPSV/TD(2,2,0:20),TU(2,2,0:20),RD(2,2,0:20),RU(2,2,0:20)
C
C *** CHECK FOR THE DIMENSION OF COMMON ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(10,*)' FROM MDRTPSV: NL.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C *** MODIFIED R/T MATRIX ***
C
      CALL MTXPSV(OM,WN,1,CMS,CMP,EA)
C
      IF(LS.EQ.1)THEN
        CALL DDUPSV(OM,WN,SDEP,CMS,CMP,DUB)
      ELSE
        CALL DDUPSV(OM,WN,THK(1),CMS,CMP,DUB)
      END IF
C
      DET=EA(3,1)*EA(4,2)-EA(3,2)*EA(4,1)
      RU(1,1,0)=-(DUB(1,1)*(EA(3,3)*EA(4,2)-EA(3,2)*EA(4,3))
     *           +DUB(2,1)*(EA(3,4)*EA(4,2)-EA(3,2)*EA(4,4)))/DET
      RU(1,2,0)=-(DUB(1,2)*(EA(3,3)*EA(4,2)-EA(3,2)*EA(4,3))
     *           +DUB(2,2)*(EA(3,4)*EA(4,2)-EA(3,2)*EA(4,4)))/DET
      RU(2,1,0)=-(DUB(1,1)*(EA(3,1)*EA(4,3)-EA(3,3)*EA(4,1))
     *           +DUB(2,1)*(EA(3,1)*EA(4,4)-EA(3,4)*EA(4,1)))/DET
      RU(2,2,0)=-(DUB(1,2)*(EA(3,1)*EA(4,3)-EA(3,3)*EA(4,1))
     *           +DUB(2,2)*(EA(3,1)*EA(4,4)-EA(3,4)*EA(4,1)))/DET
C
      IF(NL.EQ.1)RETURN
C
      IF(LS.EQ.1)THEN
        CALL DDUPSV(OM,WN,THK(1)-SDEP,CMS,CMP,DDA)
      ELSE
        CALL DDUPSV(OM,WN,THK(1),CMS,CMP,DDA)
      END IF
C
      DO 100 I=1, NL-1
C
        CALL MTXPSV(OM,WN,I+1,CMS,CMP,EB)
C
        IF(LS.EQ.(I+1))THEN
          CALL DDUPSV(OM,WN,SDEP,CMS,CMP,DUB)
        ELSE
          CALL DDUPSV(OM,WN,THK(I+1),CMS,CMP,DUB)
        END IF
C
        CALL MTRPSV(EA,EB,DDA,DUB,I) 
C
        DO 110 J=1,4
          DO 120 K=1,4
            EA(J,K)=EB(J,K)
  120     CONTINUE
  110   CONTINUE
C
        IF(LS.EQ.(I+1))THEN
          IF(LS.NE.NL)THEN
            CALL DDUPSV(OM,WN,THK(I+1)-SDEP,CMS,CMP,DDA)
          END IF
        ELSE
          CALL DDUPSV(OM,WN,THK(I+1),CMS,CMP,DDA)
        END IF
C
  100 CONTINUE
C
      RETURN
      END
C
C
C *** P-SV MODIFIED R/T MATRIX FOR Ith LAYER ***
C
        SUBROUTINE MTRPSV(EA,EB,DDA,DUB,L) 
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 EA(4,4),EB(4,4),DDA(2,2),DUB(2,2),TD,TU,RD,RU
     *          ,E1(4,4),E2(4,4),E3(4,4)
     * ,a1(4),a2(4),a3(4),b1(4),b2(4),b3(4)  
     * ,c1(4),c2(4),c3(4),d1(4),d2(4),d3(4)  
     * ,aa(4),EE(4,4)
      COMMON /MRTPSV/TD(2,2,0:20),TU(2,2,0:20),RD(2,2,0:20),RU(2,2,0:20)
     *    /STRC/DNS(20),VP(20),QP(20),FP(20)
     *                 ,VS(20),QS(20),FS(20),THK(20),NL
C
C *** CHECK FOR THE DIMENSION OF COMMON ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(10,*)' FROM MTRPSV: NL.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C *** MATRIX ***
C
      DO 100 I=1,2
        DO 110 J=1,2
          EE(I,J)= EB(I,J)
          EE(I,J+2)=-EA(I,J+2)
          EE(I+2,J)= EB(I+2,J)
          EE(I+2,J+2)=-EA(I+2,J+2)
          E2(I,J)= EA(I,J)
          E2(I,J+2)=-EB(I,J+2)
          E2(I+2,J)= EA(I+2,J)
          E2(I+2,J+2)=-EB(I+2,J+2)
  110   CONTINUE
  100 CONTINUE
C
C     * INVERSE MATRIX OF EE *
C
        DO 222 I=1,4
C
        IF(I.EQ.1)THEN
          J=2
          K=3
          M=4
        ELSE
          IF(I.EQ.2)THEN
            J=3
            K=4
            M=1
          ELSE
            IF(I.EQ.3)THEN
              J=4
              K=1
              M=2
            ELSE
              J=1
              K=2
              M=3
            END IF
          END IF
        END IF
C
        A1(I)=EE(I,4)*EE(J,1)-EE(I,1)*EE(J,4)  
        A2(I)=EE(I,4)*EE(J,2)-EE(I,2)*EE(J,4)  
        A3(I)=EE(I,4)*EE(J,3)-EE(I,3)*EE(J,4)  
        B1(I)=EE(I,3)*EE(J,2)-EE(I,2)*EE(J,3)  
        B2(I)=EE(I,1)*EE(J,3)-EE(I,3)*EE(J,1)  
        B3(I)=EE(I,2)*EE(J,1)-EE(I,1)*EE(J,2)  
        C1(I)=EE(K,3)*EE(M,2)-EE(K,2)*EE(M,3)  
        C2(I)=EE(K,1)*EE(M,3)-EE(K,3)*EE(M,1)  
        C3(I)=EE(K,2)*EE(M,1)-EE(K,1)*EE(M,2)  
        D1(I)=EE(K,4)*EE(M,1)-EE(K,1)*EE(M,4)  
        D2(I)=EE(K,4)*EE(M,2)-EE(K,2)*EE(M,4)  
        D3(I)=EE(K,4)*EE(M,3)-EE(K,3)*EE(M,4)  
C
        AA(I)=A1(I)*C1(I)+A2(I)*C2(I)+A3(I)*C3(I)
     *       +B1(I)*D1(I)+B2(I)*D2(I)+B3(I)*D3(I)
C
        E1(1,I)=(-EE(J,2)*D3(I)+EE(J,3)*D2(I)-EE(J,4)*C1(I))/AA(I)
        E1(2,I)=( EE(J,1)*D3(I)-EE(J,3)*D1(I)-EE(J,4)*C2(I))/AA(I)
        E1(3,I)=(-EE(J,1)*D2(I)+EE(J,2)*D1(I)-EE(J,4)*C3(I))/AA(I)
        E1(4,I)=( EE(J,1)*C1(I)+EE(J,2)*C2(I)+EE(J,3)*C3(I))/AA(I)
C
  222   CONTINUE
C
      DO 200 I=1,4
        DO 210 J=1,4
          E3(I,J)=(0.D0,0.D0)
          DO 220 K=1,4
            E3(I,J)=E3(I,J)+E1(I,K)*E2(K,J) 
  220     CONTINUE
  210   CONTINUE
  200 CONTINUE
C
      TD(1,1,L)=E3(1,1)*DDA(1,1)+E3(1,2)*DDA(2,1)  
      TD(2,1,L)=E3(2,1)*DDA(1,1)+E3(2,2)*DDA(2,1)  
      TD(1,2,L)=E3(1,1)*DDA(1,2)+E3(1,2)*DDA(2,2)  
      TD(2,2,L)=E3(2,1)*DDA(1,2)+E3(2,2)*DDA(2,2)  
      RD(1,1,L)=E3(3,1)*DDA(1,1)+E3(3,2)*DDA(2,1)  
      RD(2,1,L)=E3(4,1)*DDA(1,1)+E3(4,2)*DDA(2,1)  
      RD(1,2,L)=E3(3,1)*DDA(1,2)+E3(3,2)*DDA(2,2)  
      RD(2,2,L)=E3(4,1)*DDA(1,2)+E3(4,2)*DDA(2,2)  
C
      IF(L.EQ.NL)RETURN
C
      RU(1,1,L)=E3(1,3)*DUB(1,1)+E3(1,4)*DUB(2,1)  
      RU(2,1,L)=E3(2,3)*DUB(1,1)+E3(2,4)*DUB(2,1)  
      RU(1,2,L)=E3(1,3)*DUB(1,2)+E3(1,4)*DUB(2,2)  
      RU(2,2,L)=E3(2,3)*DUB(1,2)+E3(2,4)*DUB(2,2)  
      TU(1,1,L)=E3(3,3)*DUB(1,1)+E3(3,4)*DUB(2,1)  
      TU(2,1,L)=E3(4,3)*DUB(1,1)+E3(4,4)*DUB(2,1)  
      TU(1,2,L)=E3(3,3)*DUB(1,2)+E3(3,4)*DUB(2,2)  
      TU(2,2,L)=E3(4,3)*DUB(1,2)+E3(4,4)*DUB(2,2)  
C
      RETURN
      END
C
C
C *** P-SV ELEMENT MATIX ***
C
      SUBROUTINE MTXPSV(OM,WN,L,CMS,CMP,EE)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VPc,VSc,VS2,WS2,CMS,VP2,WP2,CMP
     *      ,SM,EE(4,4),X0,AA,CMSG,CMPG,OM,OM2
      COMMON/STRC/DNS(20),VP(20),QP(20),FP(20)
     *                   ,VS(20),QS(20),FS(20),THK(20),NL
C
C *** CHECK FOR THE DIMENSION OF COMMOM ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(10,*)' FROM MTXPSV: NL.GT.NC !!!!!!!!!!!!!!!'
         STOP
      END IF
C
C *** ELEMENT MATRIX ***
C
      QPe=FuncQP(OM,L)
      QSe=FuncQS(OM,L)
      VPc=VP(L)*DCMPLX(1.D0,-0.5D0/QPe)
      VSc=VS(L)*DCMPLX(1.D0,-0.5D0/QSe)
      VP2=VPc*VPc
      VS2=VSc*VSc
      SM=DNS(L)*VS2
C
      IF(OM.EQ.0.D0)THEN
C  
        X0=(1.D0+VS2/VP2)/(1.D0-VS2/VP2)
        AA=WN*SM
C
        EE(1,1)=(1.D0,0.D0)
        EE(1,2)=(1.D0,0.D0)
        EE(1,3)=(1.D0,0.D0)
        EE(1,4)=(1.D0,0.D0)
        EE(2,1)=-(X0-1.D0)
        EE(2,2)=(1.D0,0.D0)
        EE(2,3)= X0-1.D0
        EE(2,4)=(-1.D0,0.D0)
        EE(3,1)= (X0-3.D0)*AA
        EE(3,2)= -2.D0*AA
        EE(3,3)=-(X0-3.D0)*AA
        EE(3,4)= 2.D0*AA
        EE(4,1)= (X0-1.D0)*AA
        EE(4,2)= -2.D0*AA
        EE(4,3)= (X0-1.D0)*AA
        EE(4,4)=-2.D0*AA
C
      ELSE
C
        WN2=WN*WN
        OM2=OM*OM
C
        WS2=OM2/VS2
        CMS=CDSQRT(WN2-WS2)
        IF(DREAL(CMS).LT.0.D0)CMS=-CMS
        CMSG=CMS/WN
C
        WP2=OM2/VP2
        CMP=CDSQRT(WN2-WP2)
        IF(DREAL(CMP).LT.0.D0)CMP=-CMP
        CMPG=CMP/WN
C
        X0=2.D0*WN-WS2/WN
C
        EE(1,1)=(-1.D0,0.D0)
        EE(1,2)= CMSG
        EE(1,3)=(-1.D0,0.D0)
        EE(1,4)= CMSG
        EE(2,1)=-CMPG
        EE(2,2)=(1.D0,0.D0)
        EE(2,3)= CMPG
        EE(2,4)=(-1.D0,0.D0)
        EE(3,1)= 2.D0*SM*CMP 
        EE(3,2)=-SM*X0 
        EE(3,3)=-2.D0*SM*CMP 
        EE(3,4)= SM*X0 
        EE(4,1)= SM*X0 
        EE(4,2)=-2.D0*SM*CMS 
        EE(4,3)= SM*X0 
        EE(4,4)=-2.D0*SM*CMS 
C
      END IF
C
      RETURN
      END
C
C
C *** P-SV UP/DOWN-GOING MATRIX ***
C
      SUBROUTINE DDUPSV(OM,WN,THK,CMS,CMP,DU)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 DU(2,2),CMS,CMP,OM
C
      IF(OM.EQ.0.D0)THEN
        EW=DEXP(-WN*THK)
        DU(1,1)= EW
        DU(1,2)= (0.D0,0.D0)
        DU(2,1)= -WN*THK*EW
        DU(2,2)= EW
      ELSE
        DU(1,1)= CDEXP(-CMP*THK)
        DU(1,2)= (0.D0,0.D0)
        DU(2,1)= (0.D0,0.D0)
        DU(2,2)= CDEXP(-CMS*THK)
      END IF
C
      RETURN
      END
c
c
c
c *** Bessel of the 0th and 1st orders ***
c
      subroutine BJ01 (X,BJ0,BJ1)
c
      IMPLICIT REAL*8 (A-H,O-Z)
      DATA a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15
     *   /  2.5         D-1,  1.111111111 D-1,
     *      6.25        D-2,  4.0         D-2,
     *      2.777777778 D-2,  2.040816327 D-2,
     *      1.5625      D-2,  1.234567901 D-2,
     *      1.0         D-2,  8.26446281  D-3,
     *      6.944444444 D-3,  5.917159763 D-3,
     *      5.102040816 D-3,  4.444444445 D-3,
     *      3.90625     D-3 /
      DATA b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,b15
     *   /  5.0         D-1,  1.666666667 D-1,
     *      8.333333333 D-2,  5.0         D-2,
     *      3.333333333 D-2,  2.380952381 D-2,
     *      1.785714286 D-2,  1.388888889 D-2,
     *      1.111111111 D-2,  9.090909091 D-3,
     *      7.575757576 D-3,  6.41025641  D-3,
     *      5.494505495 D-3,  4.761904762 D-3,
     *      4.166666666 D-3 /
      DATA pi / 3.141592654 D0 /
c
      if(dabs(x).gt.9.d0)goto 10
c
      a=1.d0
      z1=x*0.5d0
      z=z1*z1
      BJ0=a-z*(a-z*a1*(a-z*a2*(a-z*a3*(a-z*a4*(a-z*a5*(a-z*a6*
     *   (a-z*a7*(a-z*a8*(a-z*a9*(a-z*a10*(a-z*a11*(a-z*a12*
     *   (a-z*a13*(a-z*a14*(a-z*a15)))))))))))))))
      BJ1=z1*(a-z*b1*(a-z*b2*(a-z*b3*(a-z*b4*(a-z*b5*(a-z*b6*
     *   (a-z*b7*(a-z*b8*(a-z*b9*(a-z*b10*(a-z*b11*(a-z*b12*
     *   (a-z*b13*(a-z*b14*(a-z*b15)))))))))))))))
      goto 20
c
   10 aa=dsqrt(pi*x)
      co=dcos(x)
      si=dsin(x)
      x1=0.125d0/x
      x2=0.0703125d0/(x*x)
      x3=0.073242188d0/(x*x*x)
      BJ0=(co*(1.d0-x1-x2+x3)+si*(1.d0+x1-x2-x3))/aa
      x1=0.375d0/x
      x2=0.1171875d0/(x*x)
      x3=0.102539063d0/(x*x*x)
      BJ1= (si*(1.d0+x1+x2-x3)-co*(1.d0-x1+x2+x3))/aa
c
   20 return
      end
c
c
c
C
C *** FOURIER TRANSFORM OR INVERSE FOURIER TRANSFORM
C
      SUBROUTINE dWVFAST(N,A,MOM,IND)
c
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 A(MOM),TEMP,THETA
C
C * FFT and inverse FFT *
C
      PAI=3.1415926540d0
      J=1
      DO 140 I=1,N
        IF(I.GE.J) GO TO 110
        TEMP=A(J)
        A(J)=A(I)
        A(I)=TEMP
  110   M=N/2
  120   IF(J.LE.M) GO TO 130
        J=J-M
        M=M/2
        IF(M.GE.2) GO TO 120
  130   J=J+M
  140 CONTINUE
      KMAX=1
  150 IF(KMAX.GE.N) RETURN
      ISTEP=KMAX*2
      DO 170 K=1,KMAX
        THETA=dCMPLX(0.d0,PAI*dFLOAT(IND*(K-1))/dFLOAT(KMAX))
        DO 160 I=K,N,ISTEP
          J=I+KMAX
          TEMP=A(J)*cdEXP(THETA)
          A(J)=A(I)-TEMP
          A(I)=A(I)+TEMP
  160   CONTINUE
  170 CONTINUE
      KMAX=ISTEP
      GO TO 150
C
      END
c
c
c
c *** Function for Quality Factor of Vs ***
c
      Function FuncQS(OM,L)
c
      IMPLICIT REAL*8 (A-H,O-Z)
      complex*16 OM
      COMMON /STRC/DNS(20),VP(20),QP(20),FP(20)
     *            ,VS(20),QS(20),FS(20),THK(20),NL
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(10,*)'FROM OTDAT:NL>NC',NL,NC
         STOP
      END IF
c
      OMR=dreal(OM)
c      if(OMR.eq.0.d0)then
c        FuncQS=QS(L)
c      else
        freq=OMR/6.2831853072d0
        FuncQS=QS(L)*freq**FS(L)
        if(FuncQS.lt.QS(L))FuncQS=QS(L)
c      end if
c
      return
      end
c
c
c
c *** Function for Quality Factor of Vp ***
c
      Function FuncQP(OM,L)
c
      IMPLICIT REAL*8 (A-H,O-Z)
      complex*16 OM
      COMMON /STRC/DNS(20),VP(20),QP(20),FP(20)
     *            ,VS(20),QS(20),FS(20),THK(20),NL
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(10,*)'FROM OTDAT:NL>NC',NL,NC
         STOP
      END IF
c
      OMR=dreal(OM)
c      if(OMR.eq.0.d0)then
c        FuncQP=QP(L)
c      else
        freq=OMR/6.2831853072d0
        FuncQP=QP(L)*freq**FP(L)
        if(FuncQP.lt.QP(L))FuncQP=QP(L)
c      end if
c
      return
      end


