c
c *** This code computes Green's functions from a fault model 
c     with sub-faults, considering 
c     1) the smooth rupture velocity using the Gaussian integration 
c        up to 36(=6*6) points on each sub-fault,
c     2) the smooth ramp function with a triangle slip velocity, and
c     3) time windows up to 6 windows.
c     "phs1.f" is needed for the normal modes before running grflt9.f.
c    (July, 1997)
C
c     References:
C     1) HISADA, Bull. Seismo. Soc. Am., V.84, Oct. 1994, (Part 1)
C     2) HISADA, Bull. Seismo. Soc. Am., V.85, Aug. 1995, (Part 2)
C     3) HISADA, Bull. Seismo. Soc. Am., V.85, Oct. 1995, (comments)
C
C *** FOR BUGS OR COMMENTS, PLEASE E-MAIL to "hisada@cc.kogakuin.ac.jp",
C     OR MAIL TO "Yoshiaki Hisada, Dept. of Architecture,
C     Kogakuin Univeristy, Nishi-Shinjuku 1-24-2, Tokyo 163-91, Japan"
C
C    (Originally made in 1992, and modified in 1995, 1996, 1997, and 2000)
c
C *** This code was modified from GRFAULT.F to save computation time
c     for large number of source points. This can be done by making 
c     the R/T coefficient library, before starting the wavenumber integration.
c
c *** For the numerical wave-number integrations, both Simpson's and Filon's
c     integration fomulas are used; First, the simpson's is used, then swiched
c     to the Filon's after a large wavenumber * distance (k*r > 10). (95/09)
c
c *** Before running this code, you must run "normal.f" to get the location of 
c     poles corresponding to the normal mode solutions; the data stored in
c     "phsl_dat" and "phsr_dat" for the phase velocities of the Love and 
c     and Rayleigh waves, respectively. (95/09)
c
c *** You can choose the ramp function or the smoothed ramp function as a source
c     slip function. The slip velocity of the ramp function is rectangular, and
c     that of the smoothed is triangle. Therefore, note that if you choose the 
c     smoothed, the maximum amplitude of its velocity is up to two times larger	
c     than those of the ramp function. (95/11)
c
C    (Modified from "grfault.f" in Sep., 95, from "grflt4.f" in Nov., 95
c          and from "grflt61.f" in Aug., 96, from "grinv0.f" and "grflt9.f")
c
C  1) This code is modified from grflt941.f to skip zero slip sub-faults
c    (Oct, 1998)    
C  2) This code is modified from grflt951.f to include the exp-type source
c     time function (Aug, 2000)    
C  3) This code is modified from grflt95e.f to include finer integration points
c     around branch points (Aug, 2000)    c
c  4) This code is modified from grflt00.f to include the Phinney method (Aug, 2000)
c  5) This code is modified from grflt00p.f to include the dynamic distribution
c     of the Gaussian integration points on the fault plane (Aug, 2000)
c
c * Parameters : 
c     MOM : Number of Omega
c     MOB : NUmber of Observation Points
c
      PROGRAM grflt06s
C
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(MOM=531,MOB=51)
      COMPLEX*16 GXs(MOM,MOB),GYs(MOM,MOB),GZs(MOM,MOB)
      COMPLEX*16 GX(MOM,MOB),GY(MOM,MOB),GZ(MOM,MOB)
C
      WRITE(*,*)
      WRITE(*,*)
      WRITE(*,*)'grflt06s.f computes strong ground motions considering'
      WRITE(*,*)' the static faulting and the frequency-dependent Q'
      write(*,*)' (Originally made in 1992, and modified in 2002)'
      write(*,*)
      WRITE(*,*)'"grflt06s.in" is the input data'
      WRITE(*,*)' Please Check "fort.10" for troubles.'
      WRITE(*,*)
     *' For questions, e-mail to "hisada@cc.kogakuin.ac.jp".'
      WRITE(*,*)' This code is open for the academic use.'
      WRITE(*,*)
     *' You must run "phs3sq.f" before running this code !!'
C
      CALL INDAT(MOM,MOB)
      CALL PHSBRC(MOM)
      write(*,*)
      write(*,*)'Skip Calculation of Static Terms? (No=0)'
	read(*,*)iskip
	if(iskip.eq.0)CALL HKQVWs(MOM,MOB,GXs,GYs,GZs)
c
      WRITE(*,*)
      write(*,*)'Skip Calculation of Dynamic Terms? (No=0)'
	read(*,*)iskip
	if(iskip.ne.0)then
	  do iom=1,mom
	    do iob=1,mob
	      GX(iom,iob)=GXs(iom,iob)
	      GY(iom,iob)=GYs(iom,iob)
	      GZ(iom,iob)=GZs(iom,iob)
	    end do
	  end do
	else
        CALL HKQVWd(MOM,MOB,GXs,GYs,GZs,GX,GY,GZ)
	end if
c
      CALL OTDAT(MOM,MOB,GX,GY,GZ)
c
      WRITE(*,*)
      WRITE(*,*)' Use "grfftsd.f" for FFT'
      WRITE(*,*)' Use "wavevds.m" for Plotting'
C
      STOP
      END
C
c
C *** DATA INPUT ***
C
      SUBROUTINE INDAT(MOM,MOB)
C
      IMPLICIT REAL*8(A-H,O-Z)
      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
     *    /INTG/VPMAX,VSMAX,RYMIN,CFINAL,N11,N12,M11
     *    /FLT/STR,DIP,RAK(1027,9),FL,FW,VR,FD(1027,9),RT(2,9),TWI
     *        ,xflt,yflt,zflt,delay
     *    /IFLT/NF,NTW,IS,ngauss,ndip,nstr,iramp
     *    /Hypo/xhypo,yhypo,zhypo /outom/ioutom
     *    /SIGN/ISIGN /OTST/IOT,IST
     *    /WVNUMs/WNRs(3),NWN1o,NWN2o,NWN1,NWN2,NWT1,NWT2
c
C *** PARAMETERS IN COMMOM STATEMENT ***
C
      NML =20
      NMOB=51
      NFLT=1027
      IF(MOB.LT.NMOB)NMOB=MOB
C
C *** OMEGA DATA ***
C
      write(*,*)
      OPEN (UNIT=7,FILE='grflt06s.in',STATUS='OLD')
C
      read(7,1100)dt,nt
      read(7,1110)period,OMI
 1100 format(/f8.0,i8)
 1110 format(2f8.0)
      duration=dt*nt
      dom=datan(1.d0)*8.d0/duration
      iom=1
      nom=duration/period
      WRITE(10,*)'FROM indat:'
      WRITE(10,*)
      WRITE(10,1130)nt,dt,period
      WRITE(10,1140)NOM,DOM,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+1).GT.MOM)THEN
        WRITE(10,*)' NOM+1 > MOM !!!!!!!!!!!!!'
        STOP
      END IF
C
C *** MEDIUM DATA ***
C
      READ(7,1200)NL
      DO 100 L=1,NL
        READ(7,*)DNS(L),VP(L),QP(L),FP(L),VS(L),QS(L),FS(L),THK(L)
  100 CONTINUE
 1200 FORMAT(/I10)
C
      WRITE(10,2200)NL
      DO 140 L=1,NL
      WRITE(10,2210)L,DNS(L),VP(L),QP(L),FP(L),VS(L),QS(L),FS(L),THK(L)
  140 CONTINUE
 2200 FORMAT(//' Total Number of Layers =',I5/
     *       /' Layer,Density,Vp(m/s),   Qp  ,   Fp  '
     *        ',Vs(m/s),   Qs  ,   Fs  ,Thickness(m)')
 2210 FORMAT(I5,F9.2,7F8.1,F15.2)
C
      IF(NL.GT.NML)THEN
        WRITE(10,*)' NL>NML !!!!!!!!!!!!!'
        STOP
      END IF
C
C *** MAX Vs and Vp and MIN Rayleigh Phase Velocity ***
C
      VSMAX=VS(NL)
      VPMAX=VP(NL)
      VSMIN=VS(1)
      MINVS=1
      DO 160 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
  160 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,1300)SLENGTH,DWIDTH
c      write(*,*)'SLENGTH,DWIDTH',SLENGTH,DWIDTH
      READ(7,1310)NSTR,NDIP,ngauss
c      write(*,*)'NSTR,NDIP,ngauss',NSTR,NDIP,ngauss
      NF=NSTR*NDIP
      IF(NF.GT.NFLT)THEN
        WRITE(10,*)' NF.GT.NFLT !!!!!!'
        STOP
      END IF
      if(ngauss.gt.6)then
        write(10,*)' ngauss > 6 !!!'
        stop
      end if
      fl=slength/nstr
      fw=dwidth/ndip
      READ(7,1320)DELAY,IRAMP
c      write(*,*)'DELAY,IRAMP',DELAY,IRAMP
      READ(7,1330)XHYPO,YHYPO,ZHYPO
c      write(*,*)'XHYPO,YHYPO,ZHYPO',XHYPO,YHYPO,ZHYPO
      READ(7,1330)xflt,yflt,zflt
c      write(*,*)'xflt,yflt,zflt',xflt,yflt,zflt
      READ(7,1330)STR,DIP
c      write(*,*)'STR,DIP',STR,DIP
      str=str*0.01745329252d0
      dip=dip*0.01745329252d0
      READ(7,1330)VR
      READ(7,1335)NTW,TWI
c      write(*,*)'Vr,NTW,TWI',Vr,NTW,TWI
      if(ntw.gt.9)then
         write(*,*)'ntw > 9 !!!'
         stop
      end if
      do i=1,NTW
        READ(7,1340)RT(1,i),RT(2,i)
c        write(*,*)'i,RT(1,i),RT(2,i)',i,RT(1,i),RT(2,i)
      end do
      read(7,*)
      do j=1,NF
        READ(7,*)(FD(j,i),i=1,NTW)
c        write(*,*)'j,FD(j,i)',j,(FD(j,i),i=1,NTW)
      end do
      read(7,*)
      do j=1,NF
        READ(7,*)(rak(j,i),i=1,NTW)
c        write(*,*)'j,rak(j,i)',j,(rak(j,i),i=1,NTW)
        do i=1,ntw
          rak(j,i)=rak(j,i)*0.01745329252d0
        end do
      end do
 1300 FORMAT(/2f10.0)
 1310 FORMAT(3I5)
 1320 FORMAT(F10.0,I5)
 1330 format(3f10.0)
 1335 format(i8,f8.0)
 1340 format(2f8.0)
c 1345 format(6f6.0)
C
      WRITE(10,2300)IRAMP,NF
      WRITE(10,2305)STR,DIP
      WRITE(10,2315)slength,dwidth
      WRITE(10,2320)FL,FW
      write(10,2325)xflt,yflt
      WRITE(10,2330)zflt,VR
      WRITE(10,2335)ntw,twi
      WRITE(10,2340)ngauss,ngauss*ngauss
      do i=1,ntw
        WRITE(10,2345)i,RT(1,i),RT(2,i)
      end do
 2300 FORMAT(//' Seismic Source Data:' 
     *        'IRAMP =',I3,' (=0: Ramp Function ;=1: Smoothed Ramp)' 
     *      //'   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)
 2325 FORMAT( '   X Fault Edge (m)=',F12.3,',  Y Fault Edge (m)=',F12.3)
 2330 FORMAT( '   Z Fault Edge (m)=',F12.3,',  Rupture Vel(m/s)=',F12.3)
 2335 FORMAT( '   Num. of Time Windows='i5,',  Interval (sec)  =',F12.4)
 2340 FORMAT( '   Max Num. of Gauss Points =',i5
     *       ,',  Max Num. of Guass Points per Sub-Faults =',i5)
 2345 FORMAT( '   Time Window=',i3,', Half Rise-Time 1 & 2 (s)=',2F10.5)
c      
      II=0
      DO 270 I=1,NDIP
      DO 270 J=1,NSTR
      II=II+1
      WRITE(10,2350)II,I,J
      WRITE(10,2355)(FD(ii,k),k=1,ntw)
      WRITE(10,2356)(rak(ii,k),k=1,ntw)
  270 CONTINUE
 2350 FORMAT(/'   Source Number   =',I4,' (',I2,',',I2,')')
 2355 FORMAT( '   Dislocation for Each Time Windows (m)=',6F9.3) 
 2356 FORMAT( '   Rake Ang. for Each Time Windows (rad)=',6F9.5) 
c
C *** NUMERICAL INTEGRAL DATA ***
C
      READ(7,1400)WNRs(2)
      WNRs(1)=WNRs(2)/100000.0d0
      READ(7,1410)NWN1o
      READ(7,1420)WNRs(3)
      READ(7,1410)NWN2o
 1400 FORMAT(/F10.0)
 1410 FORMAT(I10)
 1420 FORMAT(F10.0)
C
      WRITE(10,2400)
      WRITE(10,2410)WNRs(1),WNRs(2)
      WRITE(10,2420)0.0,WNRs(3)
      WRITE(10,2430)NWN1o,NWN2o
 2400 FORMAT(//' Static Wavenum. Integ. Data for Greenfield Method:')
 2410 FORMAT( /' The first integ. range of w*k on real axis   =',2e12.4)
 2420 FORMAT(  ' The second integ. range of w*k on imag. axis =',2e12.4)
 2430 FORMAT(  ' Initial Integ. Num. of Adaptive Newton-Cotes =',2I12)
C
      READ(7,1450)N11
      READ(7,1460)N12
      READ(7,1470)cfac
      CFINAL=rymin/cfac
 1450 FORMAT(/I10)
 1460 FORMAT(I10)
 1470 FORMAT(F10.0)
C
      WRITE(10,2450)
      WRITE(10,2460)N11
      WRITE(10,2470)N12
      WRITE(10,2480)CFINAL
      WRITE(10,2490)2*(N11+N12)+1
 2450 FORMAT(//' Number of Points for Numerical Integrations:')
 2460 FORMAT( /'   From 0              to 1.2*om/Ray(min) =',I5)
 2470 FORMAT(  '   From 1.2*om/Ry(min) to om/C(final)     =',I5)
 2480 FORMAT(  '   C(final = Minimun Rayleigh Phase Vel.) =',F10.3)
 2490 FORMAT(  '   Initial Number of Points: 2*(N1+N2)+1  =',I5)
c
C *** CHANGE OF SIGNS FOR IMAGINARY PARTS ***
C
      READ(7,1700)ISIGN
      IF(ISIGN.EQ.1)WRITE(10,2700)
      IF(ISIGN.EQ.0)WRITE(10,2710)
 1700 FORMAT(/I10)
 2700 FORMAT(/' SIGNS OF IMAGINARY PARTS FOR RESULTS ARE CHANGED')
 2710 FORMAT(/' SIGNS OF IMAGINARY PARTS FOR RESULTS ARE NOT CHANGED')
C
C *** OP DATA ***
C
      READ(7,1500)NOB
      WRITE(10,2500)NOB
 1500 FORMAT(/I10)
 2500 FORMAT(//' TOTAL NUMBER OF OBSERVATION POINTS =',I6/)
C
      IF(NOB.GT.NMOB)THEN
        WRITE(10,*)' NOB>NMOB !!!!!!!!!!!!!'
        STOP
      END IF
C
      DO 300 I=1,NOB
        READ(7,1520)XO(I),YO(I),ZO(I)
        WRITE(10,2520)I,XO(I),YO(I),ZO(I)
  300 CONTINUE
 1520 FORMAT(3F10.0)
 2520 FORMAT('   OB NUM.=',I5,',     (X,Y,Z) =',3F12.4)
C
      CLOSE(UNIT=7)
C
c *** output of fault and observation data ***
c
      open(9,file='fault_dat',status='unknown')
        write(9,9500)XHYPO,YHYPO,ZHYPO
        write(9,9500)xflt,yflt,zflt
        write(9,9500)STR,DIP,RAK(1,1)
        write(9,9500)slength,dwidth,VR
        write(9,9510)NSTR,NDIP,ngauss
c        do 700 iflt=1,nf
c        do 700 id=1,ngauss
c        do 700 ig=1,ngauss
c        write(9,9500)xf(iflt,id,ig),yf(iflt,id,ig),dp(iflt,id)
c  700   continue
 9500   format(3f12.4)
 9510   format(3i12)
      close(9) 
c
      open(9,file='obs_dat',status='unknown')
        do i=1,nob
        write(9,9500)xo(i),yo(i),zo(i)
        end do
      close(9) 
C
C *** OUTPUT OF INTEGRAND OF U,V,AND W ***
C
   50 IOT=0
      IF(NOB.NE.1)RETURN
      IF(NF.NE.1)RETURN
c      IF((NOM-IOM).NE.0)RETURN
C
c      WRITE(*,7000)
c      READ(*,*)IOT
c 7000 FORMAT(/' Output Integrands?: Enter No=0 or Yes=1: ')
C
      df0=dom/6.283185307d0
      if(iot.eq.1)then
        f1=0.d0
        f2=df0*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+1
          freq=df0*(i-1)
          dfreq=dabs(freq-ofreq)
          if(dfreq.lt.df)then
            ioutom=i
            df=dfreq
          end if
        end do
        write(*,*)'Actual Output Frequency (Hz)=',df0*(ioutom-1)
        write(10,*)'Actual Output Frequency (Hz)=',df0*(ioutom-1)
        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 *** PHASE VELOCITY OF THE RAYLEIGH WAVE IN THE HOMOGENEOUS HALF-SPACE ***
C
      FUNCTION FRHF(VP,VS)
      IMPLICIT REAL*8(A-H,O-Z)
C
      DZ=0.0001D0
C
      SP=VS/VP
      IF(SP.GE.0.0D0) Z= 0.91263D0
      IF(SP.GE.0.2D0) Z= 0.90767D0
      IF(SP.GE.0.4D0) Z= 0.88899D0
      IF(SP.GE.0.6D0) Z= 0.83575D0
      IF(SP.GE.0.7D0) Z= 0.77062D0
      IF(SP.GE.0.8D0) Z= 0.63338D0
C
      F0=FRHF1(SP,Z)
  100 Z=Z-DZ
      F1=FRHF1(SP,Z)
      AF=F0/F1
      F00=F0
      F0 =F1
      IF(AF.GE.0.D0)GOTO 100
      FZ=Z+DZ*F0/(F0-F00)
C
      FRHF=VS*DSQRT(FZ)
C
      RETURN
      END
C
C
      FUNCTION FRHF1(SP,Z)
      IMPLICIT REAL*8(A-H,O-Z)
C
      Z2=Z*Z
      Z3=Z2*Z
C
      F1=16.D0*(Z-1.D0)*SP*SP
      FRHF1=F1-Z3+8.D0*Z2-24.D0*Z+16.D0
C
      RETURN
      END
C
C
C *** DATA OUTPUT ***
C
      SUBROUTINE OTDAT(MOM,MOB,GX,GY,GZ)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 GX(MOM,MOB),GY(MOM,MOB),GZ(MOM,MOB),AA,CC
      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(1027,9),FL,FW,VR,FD(1027,9),RT(2,9),TWI
     *        ,xflt,yflt,zflt,delay 
     *    /IFLT/NF,NTW,IS,ngauss,ndip,nstr,iramp
     *    /SIGN/ISIGN /Hypo/xhypo,yhypo,zhypo 

C
C *** CHECK ***
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(10,*)'FROM OTDAT:NL>NC',NL,NC
         STOP
      END IF
C
      NC=51
      IF(NOB.GT.NC)THEN
         WRITE(10,*)'FROM OTDAT:NOB>NC',NOB,NC
         STOP
      END IF
C
      NC=1027
      IF(NF.GT.NC)THEN
         WRITE(10,*)'FROM OTDAT:NF>NC',NF,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+1
          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.dxyz',STATUS='UNKNOWN')
        WRITE(8,2500)DOM,IOM,NOM,OMI
        WRITE(8,2510)NOB
        DO 500 IOB=1,NOB
          DO 510 K=IOM,NOM+1
            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,F10.7)
 2510 FORMAT(1X,I10,///)
 2520 FORMAT(1X,6E12.5)
C
C *** OUTPUT OF DATA (GQ, GV AND GW)***
C
        DO 600 IOB=1,NOB
          RX=XO(IOB)-xhypo
          RY=YO(IOB)-yhypo
          RR=DSQRT(RX*RX+RY*RY)
          CR=RX/RR
          SR=RY/RR
          DO 610 K=IOM,NOM+1
            AA=GX(K,IOB)
            CC=GY(K,IOB)
            GX(K,IOB)= AA*CR+CC*SR
            GY(K,IOB)=-AA*SR+CC*CR
  610     CONTINUE
  600   CONTINUE
C      
      OPEN(8,FILE='grfault.dqvw',STATUS='UNKNOWN')
        WRITE(8,3500)DOM,IOM,NOM+1,OMI
        WRITE(8,3510)NOB
        DO 800 IOB=1,NOB
          DO 810 K=IOM,NOM+1
            WRITE(8,3520)GX(K,IOB),GY(K,IOB),GZ(K,IOB)
  810     CONTINUE
  800   CONTINUE
      CLOSE(8)
C
 3500 FORMAT(1X,F10.6,2I10,F10.7)
 3510 FORMAT(1X,I10,///)
 3520 FORMAT(1X,6E12.5)
C
      RETURN
      END
C
C *** COMPUTATION OF Static GREEN FUNCTIONS USING THE R/T METHOD ***
C
      SUBROUTINE HKQVWs(MOM,MOB,GXS,GYS,GZS)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 SM,RAM,OM,OM1,im,VSc,VS2,VPc,VP2
     *      ,GXS(MOM,MOB),GYS(MOM,MOB),GZS(MOM,MOB)
     *      ,QXS(0:2),QZS(0:2),VXS(0:2),WXS(0:2),WZS(0:2)
     *      ,QXG(0:2),QZG(0:2),VXG(0:2),WXG(0:2),WZG(0:2)
     *      ,SMO,EDT0,EDT,GX,GY,GZ,QA,QB,VA,VB,WA,WB,GX0,GY0,GZ0
     *      ,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),gw(8,8),gaux(8,8)
      dimension dc(1027),sc(1027),RRm(1027)        
      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(1027,9),FL,FW,VR,FD(1027,9),RT(2,9),TWI
     *     ,xflt,yflt,zflt,delay  
     *    /IFLT/NF,NTW,IS,ngauss,ndip,nstr,iramp 
     *    /Hypo/xhypo,yhypo,zhypo 
     *    /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
     *    /WVNUMs/WNRs(3),NWN1o,NWN2o,NWN1,NWN2,NWT1,NWT2
c
      data gaux(1,1)/0.0d0/
      data gaux(1,2),gaux(2,2)
     * /-0.577350269189626d0, 0.577350269189626d0/
      data gaux(1,3),gaux(2,3),gaux(3,3)
     * /-0.774596669241483d0,0.d0, 0.774596669241483d0/
      data gaux(1,4),gaux(2,4),gaux(3,4),gaux(4,4)
     * /-0.861136311594053d0,-0.339981043584856d0
     * , 0.339981043584856d0, 0.861136311594053d0/
      data gaux(1,5),gaux(2,5),gaux(3,5),gaux(4,5),gaux(5,5)
     * /-0.906179845938664d0,-0.538469310105683d0,0.0d0
     * , 0.538469310105683d0, 0.906179845938664d0/
      data gaux(1,6),gaux(2,6),gaux(3,6),gaux(4,6),gaux(5,6),gaux(6,6)
     * /-0.932469514203152d0,-0.661209386466265d0,-0.238619186083197d0
     * , 0.238619186083197d0, 0.661209386466265d0, 0.932469514203152d0/ 
c
      data gw(1,1) /2.d0/
      data gw(1,2),gw(2,2) /1.d0,1.d0/
      data gw(1,3),gw(2,3),gw(3,3)
     * / 0.555555555555556d0, 0.888888888888889d0, 0.555555555555556d0/
      data gw(1,4),gw(2,4),gw(3,4),gw(4,4)
     * / 0.347854845137454d0, 0.652145154862546d0
     * , 0.652145154862546d0, 0.347854845137454d0/
      data gw(1,5),gw(2,5),gw(3,5),gw(4,5),gw(5,5)
     * / 0.236926885056189d0, 0.478628670499366d0, 0.568888888888889d0
     * , 0.478628670499366d0, 0.236926885056189d0/
      data gw(1,6),gw(2,6),gw(3,6),gw(4,6),gw(5,6),gw(6,6)
     * / 0.171324492379170d0, 0.360761573048139d0, 0.467913934572691d0
     * , 0.467913934572691d0, 0.360761573048139d0, 0.171324492379170d0/
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=1027
      IF(NF.GE.NC)THEN
        WRITE(10,*)' FROM HKQVW:NF>NC !!!!!!!!!!'
        WRITE(10,*)' NF,NC=',NF,NC
        STOP
      END IF
C
c     * Omega Data *
c
          OM=dcmplx(0.d0,OMI)
C
C   * output of integrands for check *
C
      write(*,*)
      write(*,*)'Output Static Integrands for checking? (no=0, yes=1)'
	read(*,*)IOT
c
      IF(IOT.eq.1)THEN
	  write(*,*)'The first OP is used as the SP location!'
	  write(*,2000)xo(1),yo(1),zo(1)
2000    format(' (x,y,z)=',3f12.4)
	  write(*,*)'Enter the location of OP (x,y,z)'
	  read(*,*)xot,yot,zot
        OPEN (20,FILE='sint1_d',STATUS='unknown')
        OPEN (21,FILE='sint2_d',STATUS='unknown')
      END IF
C
c
c *** The fault integration for the static (asymptotic) solutions ***
c
c      write(1,*)'STR,DIP',STR,DIP
c	write(1,*)'FL,FW',FL,FW
c      write(1,*)'xflt,yflt,zflt',xflt,yflt,zflt  

      im=DCMPLX(0.d0,1.d0)
      cdp=dcos(dip)
      sdp=dsin(dip)
      cst=dcos(str)
      sst=dsin(str)
c
      DO 100 IOB=1,NOB
c
        if(iot.eq.1)then
	    dc(1)=1.d0
		sc(1)=1.d0
		RRm(1)=9999999.d0
		goto 25
	  end if
c
c     * Finding the closest distance between OB and sub-fault planes *
c
        DO 200 idip=1,ndip 
c
          d10=FW*(idip-1)
          d20=FW*idip
          zflt1=zflt-d10*sdp
          zflt2=zflt-d20*sdp
c
          DO 220 istr=1,nstr
c
c          write(1,*)
c          write(1,*)'iob,idip,istr  =',iob,idip,istr
c            
            s10=FL*(istr-1)
            s20=FL*istr
            xflt1=xflt+s10*cst+d10*cdp*sst
            xflt2=xflt+s20*cst+d20*cdp*sst
            yflt1=yflt+s10*sst-d10*cdp*cst
            yflt2=yflt+s20*sst-d20*cdp*cst
c
c            write(1,*)'Flt Corner 1: x,y,z',xflt1,yflt1,zflt1
c            write(1,*)'Flt Corner 2: x,y,z',xflt2,yflt2,zflt2
c
            dRR=0.0000001d0
            Rmin=99999999999.d0
            Rold=Rmin
            iloop=1
            d0=d10
            s0=s10
            Dwid=FW
            Slen=FL
c
 10         d1=d0+Dwid/4.d0
            d2=d0+Dwid*3.d0/4.d0
            Z1=zflt-d1*sdp
            Z2=zflt-d2*sdp
c            if(Z2.gt.Z1)then
c              if(Z1.gt.zflt2)Z1=zflt2
c              if(Z1.lt.zflt1)Z1=zflt1
c              if(Z2.gt.zflt2)Z2=zflt2
c              if(Z2.lt.zflt1)Z2=zflt1
c            else
c              if(Z1.gt.zflt1)Z1=zflt1
c              if(Z1.lt.zflt2)Z1=zflt2
c              if(Z2.gt.zflt1)Z2=zflt1
c              if(Z2.lt.zflt2)Z2=zflt2
c            end if
            RZ1=Z1-ZO(IOB)
            RZ2=Z2-ZO(IOB)
c       write(1,*)
c       write(1,*)'1:iloop',iloop
c       write(1,*)'1:d0,s0',d0,s0
c       write(1,*)'dwid,slen',dwid,slen
c       write(1,*)'d1,s1',d1,s1
c       write(1,*)'d2,s2',d2,s2
c       write(1,*)'z1,z2',z1,z2
c
            s1=s0+Slen/4.d0
            s2=s0+Slen*3.d0/4.d0
            X1=xflt+s1*cst+d1*cdp*sst
            X2=xflt+s2*cst+d2*cdp*sst
            Y1=yflt+s1*sst-d1*cdp*cst
            Y2=yflt+s2*sst-d2*cdp*cst
c
            X11=xflt+s1*cst+d1*cdp*sst
            Y11=yflt+s1*sst-d1*cdp*cst
            X12=xflt+s2*cst+d1*cdp*sst
            Y12=yflt+s2*sst-d1*cdp*cst
            X21=xflt+s1*cst+d2*cdp*sst
            Y21=yflt+s1*sst-d2*cdp*cst
            X22=xflt+s2*cst+d2*cdp*sst
            Y22=yflt+s2*sst-d2*cdp*cst
c
c            if(X2.gt.X1)then
c              if(X11.gt.xflt2)X11=xflt2
c              if(X11.lt.xflt1)X11=xflt1
c              if(X12.gt.xflt2)X12=xflt2
c              if(X12.lt.xflt1)X12=xflt1
c              if(X21.gt.xflt2)X21=xflt2
c              if(X21.lt.xflt1)X21=xflt1
c              if(X22.gt.xflt2)X22=xflt2
c              if(X22.lt.xflt1)X22=xflt1
c            else
c              if(X11.gt.xflt1)X11=xflt1
c              if(X11.lt.xflt2)X11=xflt2
c              if(X12.gt.xflt1)X12=xflt1
c              if(X12.lt.xflt2)X12=xflt2
c              if(X21.gt.xflt1)X21=xflt1
c              if(X21.lt.xflt2)X21=xflt2
c              if(X22.gt.xflt1)X22=xflt1
c              if(X22.lt.xflt2)X22=xflt2
c            end if
c            if(Y2.gt.Y1)then
c              if(Y11.gt.Yflt2)Y11=Yflt2
c              if(Y11.lt.Yflt1)Y11=Yflt1
c              if(Y12.gt.Yflt2)Y12=Yflt2
c              if(Y12.lt.Yflt1)Y12=Yflt1
c              if(Y21.gt.Yflt2)Y21=Yflt2
c              if(Y21.lt.Yflt1)Y21=Yflt1
c              if(Y22.gt.Yflt2)Y22=Yflt2
c              if(Y22.lt.Yflt1)Y22=Yflt1
c            else
c              if(Y11.gt.Yflt1)Y11=Yflt1
c              if(Y11.lt.Yflt2)Y11=Yflt2
c              if(Y12.gt.Yflt1)Y12=Yflt1
c              if(Y12.lt.Yflt2)Y12=Yflt2
c              if(Y21.gt.Yflt1)Y21=Yflt1
c              if(Y21.lt.Yflt2)Y21=Yflt2
c              if(Y22.gt.Yflt1)Y22=Yflt1
c              if(Y22.lt.Yflt2)Y22=Yflt2
c            end if
c
c        write(1,*)'x11,y11',x11,y11
c        write(1,*)'x12,y12',x12,y12
c        write(1,*)'x21,y21',x21,y21
c        write(1,*)'x22,y22',x22,y22
c
            RX11=X11-XO(IOB)
            RX12=X12-XO(IOB)
            RX21=X21-XO(IOB)
            RX22=X22-XO(IOB)
            RY11=Y11-YO(IOB)
            RY12=Y12-YO(IOB)
            RY21=Y21-YO(IOB)
            RY22=Y22-YO(IOB)
c      write(1,*)'xo,yo,zo',xo(IOB),yo(IOB),zo(IOB)
c      write(1,*)'rx11,ry11',rx11,ry11
c      write(1,*)'rx12,ry12',rx12,ry12
c      write(1,*)'rx21,ry21',rx21,ry21
c      write(1,*)'rx22,ry22',rx22,ry22
c      write(1,*)'rz1,rz2',rz1,rz2
            R11=DSQRT(RX11*RX11+RY11*RY11+RZ1*RZ1)
            R12=DSQRT(RX12*RX12+RY12*RY12+RZ1*RZ1)
            R21=DSQRT(RX21*RX21+RY21*RY21+RZ2*RZ2)
            R22=DSQRT(RX22*RX22+RY22*RY22+RZ2*RZ2)
c      write(1,*)'R11,R12',R11,R12
c       write(1,*)'R21,R22',R21,R22
c        
            imin=0
            if(R11.lt.Rmin)then
              imin=11
              Rmin=R11*(1.d0-dRR/10000.d0)
            end if
            if(R12.lt.Rmin)then
              imin=12
              Rmin=R12*(1.d0-dRR/10000.d0)
            end if
            if(R21.lt.Rmin)then
              imin=21
              Rmin=R21*(1.d0-dRR/10000.d0)
            end if
            if(R22.lt.Rmin)then
              imin=22
              Rmin=R22*(1.d0-dRR/10000.d0)
            end if
c        
c      write(1,*)'imin',imin
c	if(imin.eq.11)then
c	  write(2,999)iob,idip,istr,imin,x11,y11,z1,r11
c	end if
c	if(imin.eq.12)then
c	  write(2,999)iob,idip,istr,imin,x12,y12,z1,r12
c	end if
c	if(imin.eq.21)then
c	  write(2,999)iob,idip,istr,imin,x21,y21,z2,r21
c	end if
c	if(imin.eq.22)then
c	  write(2,999)iob,idip,istr,imin,x22,y22,z2,r22
c	end if
c999	format(4i5,4f12.4)
c
            if(imin.ne.0)then
c
              if(imin.eq.12)then
                s0=s0+Slen/2.d0
              else
                if(imin.eq.21)then
                  d0=d0+Dwid/2.d0
                else 
                  if(imin.eq.22)then
                    s0=s0+Slen/2.d0
                    d0=d0+Dwid/2.d0
                  end if
                end if
              end if
c
            else
c
              Rmin=999999999999.d0
              jmin=0
              if(R11.lt.Rmin)then
                jmin=11
                Rmin=R11*(1.d0-dRR/10000.d0)
              end if 
              if(R12.lt.Rmin)then
                jmin=12
                Rmin=R12*(1.d0-dRR/10000.d0)
              end if        
              if(R21.lt.Rmin)then
                jmin=21
                Rmin=R21*(1.d0-dRR/10000.d0)
              end if 
              if(R22.lt.Rmin)then
                jmin=22
                Rmin=R22*(1.d0-dRR/10000.d0)
              end if 
c
c       write(1,*)'jmin',jmin
c
              if(jmin.eq.12)then
                s0=s0+Slen/2.d0
              else
                if(jmin.eq.21)then
                  d0=d0+Dwid/2.d0
                else 
                  if(jmin.eq.22)then
                    s0=s0+Slen/2.d0
                    d0=d0+Dwid/2.d0
                  end if
                end if
              end if
c
            end if 
c
c 20         write(1,*)'2:iloop,imin',iloop,imin
c        write(1,*)'Rold,Rmin',Rold,Rmin
c        write(1,*)'2:d0,s0',d0,s0

c      if(dabs(Rmin-Rold).gt.dRR)then
 20     if(dabs(Rmin-Rold).gt.dRR)then
c
              iloop=iloop+1
              Dwid=Dwid/2.d0
              Slen=Slen/2.d0
              if(iloop.gt.100)then
                write(*,*)'iloop > 100 !!!'
                write(*,*)'The observation point might be too '
                write(*,*)'close to the fault plane!: IOB=',IOB
                stop
              end if
              Rold=Rmin    
              goto 10
c
            else
c
              Rmin0=R11
              if(R12.lt.Rmin0)then
                Rmin0=R12
              else
                if(R21.lt.Rmin0)then
                  Rmin0=R21
                else
                  if(R22.lt.Rmin0)then
                    Rmin0=R22
                  end if
                end if
              end if
c
              Rmax0=R11
              if(R12.gt.Rmax0)then
                Rmax0=R12
              else
                if(R21.gt.Rmax0)then
                  Rmax0=R21
                else
                  if(R22.gt.Rmax0)then
                    Rmax0=R22
                  end if
                end if
              end if
c
c      write(1,*)
c      write(1,*)'Rmax0,Rmin0',Rmax0,Rmin0
c
              if((Rmax0-Rmin0).gt.dRR)then
                iloop=iloop+1
                Dwid=Dwid/2.d0
                Slen=Slen/2.d0
                if(iloop.gt.100)then
                  write(*,*)'iloop > 100 !!!'
                  write(*,*)'The observation point might be too '
                  write(*,*)'close to the fault plane!: IOB=',IOB
                  stop
                end if
                Rold=Rmin    
                goto 10
              end if
c
            end if
c
            xc=(x1+x2)/2.d0
            yc=(y1+y2)/2.d0
            zc=(z1+z2)/2.d0
c
            IS=nstr*(idip-1)+istr
            dc(IS)=(d1+d2)/2.d0        
            sc(IS)=(s1+s2)/2.d0
            RRm(IS)=Rmin        

c        write(1,*)'FW,FL',FW,FL
c        write(1,*)
c        write(1,*)'dc,sc',dc(IS),sc(IS)
c        write(1,9990)xc,yc,zc
        write(*,*)'dc,sc',dc(IS),sc(IS)
        write(*,9990)xc,yc,zc
 9990   format(' xc,yc,zc=',3f12.3)
c        write(1,*)'Rmin,iloop',Rmin,iloop
        write(*,*)'Rmin,iloop',Rmin,iloop
c	stop
c
 220    CONTINUE
 200    CONTINUE

        write(*,*)
c
c *** The fault integration for the static solutions ***
C
C
  25    ZS=ZO(IOB)
        if(iob.ge.2)then
          if(zs.eq.zo(iob-1))goto 30
        end if
c
        DRU=0.0
        DRL=THK(1)
        DO 250 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)
 250    CONTINUE
C
        DEP1=0.0
        DO 270 L=1,LS-1
            DEP1=DEP1+THK(L)
 270    CONTINUE
        SDEP=ZS-DEP1
        DEP2=DEP1+THK(LS)
c
c     * Starting Fault Integration *
c
 30     GX0=(0.D0,0.D0)
        GY0=(0.D0,0.D0)
        GZ0=(0.D0,0.D0)
c
        ngdipT=0
        ngstrT=0
        ngT=0
        NWT1=0
        NWT2=0
        write(10,*)
        write(10,*)'from HKQVWs: iob=',iob
        DO 300 idip=1,ndip 
c
            da=FW*(idip-1)
            db=FW*idip
            za=zflt-da*sdp*1.0000001d0
            zb=zflt-db*sdp*0.9999999d0
            if(zb.lt.0.d0)then
              write(*,*)'Fault Depth is less than 0!: Dip Number=',idip
              stop
            end if
            DRU=0.0
            DRL=THK(1)
            do L=1,NL
              IF((za.GE.DRU).AND.(za.LT.DRL))L1=L
              IF((za.GE.DRU).AND.(L.EQ.NL))  L1=L
              DRU=DRL
              DRL=DRL+THK(L+1)
            end do
            DRU=0.0
            DRL=THK(1)
            do L=1,NL
              IF((zb.GE.DRU).AND.(zb.LT.DRL))L2=L
              IF((zb.GE.DRU).AND.(L.EQ.NL))  L2=L
              DRU=DRL
              DRL=DRL+THK(L+1)
            end do
            if(L1.ne.L2)then
              write(*,*)' The L1 Layer is different from the L2 Layer!'
              write(*,*)' Layers 1 and 2 =',L1,L2
c             stop
            end if
c
          DO 320 istr=1,nstr
c
          write(*,*)
          write(*,*)'iob,idip,istr',iob,idip,istr
            IS=nstr*(idip-1)+istr
            rWS=FW/RRm(IS)
            rLS=FL/RRm(IS)
            NWS=rWS
            NLS=rLS
            if(rWS.le.0.8)NWS=2            
            if(rLS.le.0.8)NLS=2         
            if(rWS.le.0.6)NWS=1            
            if(rLS.le.0.6)NLS=1         
            if(NWS.lt.1)NWS=1            
            if(NLS.lt.1)NLS=1         
c         write(*,*)'rWS,rLS',rWS,rLS
         write(*,*)'NWS,NLS',NWS,NLS
            if(NWS.gt.500)then
               write(*,*)'Warning! NWS > 500',NWS
               write(*,*)'The OB might be too close to the fault !'
               write(*,*)'iob,idip,istr',iob,idip,istr
            end if
            if(NLS.gt.500)then
               write(*,*)'Warning! NLS > 500',NLS
               write(*,*)'The OB might be too close to the fault !'
               write(*,*)'iob,idip,istr',iob,idip,istr
            end if
c
            dWS=FW/dfloat(NWS)
            dLS=FL/dfloat(NLS)
            if(dWS.gt.FW)dWS=FW
            if(dLS.gt.FL)dLS=FL
c         write(*,*)'dWS,dLS',dWS,dLS
c
           ngdipS=0
           ngstrS=0
           ngS=0
           do 400 iWS=1,NWS
c
            da=FW*(idip-1)+dWS*(iWS-1)
            db=FW*(idip-1)+dWS*iWS
            dcn=(da+db)/2.d0
            za=zflt-da*sdp*1.0000001d0
            zb=zflt-db*sdp*0.9999999d0
            zcn=(za+zb)/2.d0
c  
            do 420 iLS=1,NLS
c
              sa=FL*(istr-1)+dLS*(iLS-1)
              sb=FL*(istr-1)+dLS*iLS
              scn=(sa+sb)/2.d0
              xcn=xflt+scn*cst+dcn*cdp*sst
              ycn=yflt+scn*sst-dcn*cdp*cst
              xn=xcn-xo(iob)
              yn=ycn-yo(iob)
              zn=zcn-zo(iob)
              Rcn=dsqrt(xn*xn+yn*yn+zn*zn)
c       write(*,*)'iWS,iLS',iWS,iLS
c       write(*,*)'da,db',da,db
c       write(*,*)'sa,sb',sa,sb
c       write(*,*)'Rcn',Rcn
              ratioW=Rcn/dWS 
              ratioL=Rcn/dLS 
c       write(*,*)'ratioW,ratioL',ratioW,ratioL

            ngausW=1
            if(ratioW.lt.100.0d0)then
              ngausW=2
              if(ratioW.lt.10.0d0)then
                ngausW=3
                if(ratioW.lt.4.0d0)then
                  ngausW=4
                  if(ratioW.lt.2.0d0)then
                    ngausW=5
                    if(ratioW.lt.1.0d0)then
                      ngausW=6
                      if(ratioW.lt.0.5d0)then
                  write(*,*)'ngausW > 6 !!!!'
                  write(*,*)'ngausW > 6 -> ngausW = 6'
                        ngausW=6
                      end if
                    end if
                  end if
                end if
              end if
            end if
c
            ngausL=1
            if(ratioL.lt.100.0d0)then
              ngausL=2
              if(ratioL.lt.10.0d0)then
                ngausL=3
                if(ratioL.lt.4.0d0)then
                  ngausL=4
                  if(ratioL.lt.2.0d0)then
                    ngausL=5
                    if(ratioL.lt.1.0d0)then
                      ngausL=6
                      if(ratioL.lt.0.5d0)then
                  write(*,*)'ngausL > 6 !!!!'
                  write(*,*)'ngausL > 6 -> ngausL = 6'
                        ngausL=6
                      end if
                    end if
                  end if
                end if
              end if
            end if
c
        ngdipS=ngdipS+ngausL
        ngstrS=ngstrS+ngausW
        ngS=ngS+ngausL*ngausW
c
c            write(*,1200)idip,za,zb         
c            write(*,1210)ratioW,ngausW
c            write(*,1220)ratioL,ngausL
c            write(10,1200)idip,za,zb         
c            write(10,1210)ratioW,ngausW
c            write(10,1220)ratioL,ngausL
c 1200     format(' Elem Num.for W, and the edge depths =',i3,2f11.2)
c 1210     format(' Dist. Ratio & Gauss Pts per Sub-W =',e12.4,' ->',i3)
c 1220     format(' Dist. Ratio & Gauss Pts per Sub-L =',e12.4,' ->',i3)
c
c           * Gaussian Quadrature (id*ig points) *
c
            do 500 id=1,ngausW
c
              zd=((1.d0-gaux(id,ngausW))*za
     *           +(1.d0+gaux(id,ngausW))*zb)/2.d0
              if(iot.eq.1)zd=zot
c  
              DRU=0.0
              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
              zz=zd-zhypo
c
              LO=L1
              VPc=VP(LO)*DCMPLX(1.D0,-0.5D0/QP(LO))
              VSc=VS(LO)*DCMPLX(1.D0,-0.5D0/QS(LO))
              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
              dd=((1.d0-gaux(id,ngausW))*da
     *           +(1.d0+gaux(id,ngausW))*db)/2.d0
c     
              do 550 ig=1,ngausL
c
                sg=((1.d0-gaux(ig,ngausL))*sa
     *             +(1.d0+gaux(ig,ngausL))*sb)/2.d0
                XF=xflt+sg*cst+dd*cdp*sst
                YF=yflt+sg*sst-dd*cdp*cst
	          if(iot.eq.1)then
	            XF=xot
	            YF=yot
	          end if
c
                xx=XF-xhypo
                yy=YF-yhypo
                xy=dsqrt(xx*xx+yy*yy+zz*zz)
                TD=delay+xy/vr            
                RX=XF-XO(IOB)
                RY=YF-YO(IOB)
                RR=DSQRT(RX*RX+RY*RY)
c
c       write(*,*)'idip,istr,id,ig',idip,istr,id,ig
c       write(*,*)'xhypo,yhypo,zhypo',xhypo,yhypo,zhypo
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 *
c
                NWN1=NWN1o
                NWN2=NWN2o
                CALL QVWRZS1(OM,QXG,QZG,VXG,WXG,WZG)
                CALL QVWRZS2(OM,QXS,QZS,VXS,WXS,WZS)
c
c      write(*,*)'xf,yf,zd',xf,yf,zd
c      write(*,*)'rx,ry,rr',rx,ry,rr
c      write(*,*)'om',om
c      write(*,*)'VX1',VXG(1)
c      write(*,*)'VX2',VXG(2)
c      write(*,*)'VX1',VXS(1),VXG(1)
c      write(*,*)'VX2',VXS(2),VXG(2)
c      write(1,*)'QX0',QXS(0),QXG(0)
c      write(1,*)'QX1',QXS(1),QXG(1)
c      write(1,*)'QX2',QXS(2),QXG(2)
c      write(1,*)'VX1',VXS(1),VXG(1)
c      write(1,*)'VX2',VXS(2),VXG(1)
c	stop
c       write(2,*)'rx,ry,rr',rx,ry,rr
c      write(2,9999)QXS(2),QXG(2),QXS(2)+QXG(2)
c 9999	format(10e12.4)
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
c       write(*,*)'QA',QA
c       write(*,*)'QB',QB
c       write(*,*)'C2',C2
c       write(*,*)'CX',CX
                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               * Time Windows *
c
            IS=nstr*(idip-1)+istr
            DO 600 IO=IOM,NOM+1
              OMR=DOM*DFLOAT(IO-1)
              OM1=dcmplx(OMR,OMI)
              EDT=CDEXP(im*OM1*TD)
                DO 620 itw=1,ntw
                  CALL MO(itw,OM1,en,SMO)
                  SMO=SMO/dfloat(nWS*nLS)
                  EDT0=CDEXP(im*OM1*TWI*(itw-1))
                  GX=(0.D0,0.D0)
                  GY=(0.D0,0.D0)
                  GZ=(0.D0,0.D0)
                  DO I=1,3
                    DO J=1,3
                      GX=GX+EN(I,J)*UX(J,I)
                      GY=GY+EN(I,J)*UY(J,I)
                      GZ=GZ+EN(I,J)*UZ(J,I)
c        write(*,*)'GZ ',GZ
                    end do
                  end do
c       write(*,*)'1 GX ',GX
c        write(*,*)'1 GZ ',GZ
                 GX0=GX*SMO*EDT*EDT0*gw(id,ngausW)*gw(ig,ngausL)/4.d0
                  GY0=GY*SMO*EDT*EDT0*gw(id,ngausW)*gw(ig,ngausL)/4.d0
                  GZ0=GZ*SMO*EDT*EDT0*gw(id,ngausW)*gw(ig,ngausL)/4.d0
c       write(*,*)'2 GX0',GX0
                  GXS(IO,IOB)=GXS(IO,IOB)+GX0
                  GYS(IO,IOB)=GYS(IO,IOB)+GY0
                  GZS(IO,IOB)=GZS(IO,IOB)+GZ0
c	write(*,*)'SMO,EDT,EDT0',SMO,EDT,EDT0
c	write(*,*)'im,OM1,TD',im,OM1,TD
c	if(io.eq.10)stop
C
                IF(IOT.eq.1)THEN
      write(*,*)'GXS',GXS(io,iob)
      write(*,*)'GYS',GYS(io,iob)
      write(*,*)'GZS',GZS(io,iob)
      write(*,*)
      write(*,*)'Use sint1.m for ploting Green functions integrands 1' 
      write(*,*)'Use sint2.m for ploting Green functions integrands 2' 
                  CLOSE(20)
                  CLOSE(21)
	            stop
                END IF
c
  620         CONTINUE   
  600       CONTINUE   
c
c    * id & ig -> ngaussW & ngaussL (Gauss Points) *
  550       CONTINUE    
  500       CONTINUE
c
c    * iWS & iLS -> nWS & nLS (Sub Sub-Faults) *
  420     CONTINUE
  400     CONTINUE
c
        write(*,9992)idip,istr,ngdipS,ngstrS,ngS
 9992   format(' dip & str=',2i4,', Gauss Pts for W & L=',2i6
     *        ,', Total Num=',i7)
        ngdipT=ngdipT+ngdipS
        ngstrT=ngstrT+ngstrS
        ngT=ngT+ngS
c
c    * idip & istr -> ndip & nstr (Sub Faults) *
  320   CONTINUE
  300   CONTINUE
c
        write(*,*)
        write(*,9994)iob,ngdipT,ngstrT,ngT
 9994   format(' iob =',i4,', Gauss Pts for W & L ='
     *,2i7,', Total Pnts Num.=',i8)
	  write(*,*)'Numbers of Integ Pnts for Range 1 & 2 =',NWT1,NWT2
	  write(10,*)'Numbers of Integ Pnts for Range 1 & 2 =',NWT1,NWT2
	  write(*,*)
        write(*,9996)dreal(GXS(1,IOB)),dreal(GYS(1,IOB)),
     *               dreal(GZS(1,IOB))
 9996   format(' Static Displacement for X, Y & Z =',3e12.4)
c
c    * iob -> nob (Observation Points) *
  100 CONTINUE
c
      RETURN
      END
c
C
C
C *** COMPUTATION OF Dynamic GREEN FUNCTIONS USING THE R/T METHOD ***
C
      SUBROUTINE HKQVWd(MOM,MOB,GXS,GYS,GZS,GGX,GGY,GGZ)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VPc,VSc,VS2,VP2,SM,RAM
     *      ,GGX(MOM,MOB),GGY(MOM,MOB),GGZ(MOM,MOB)
     *      ,GXS(MOM,MOB),GYS(MOM,MOB),GZS(MOM,MOB),OM,im
     *      ,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)
     *      ,SMO,EDT0,EDT,GX,GY,GZ,GX0,GY0,GZ0
	* 	  ,QA,QB,VA,VB,WA,WB
     *      ,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),gw(8,8),gaux(8,8)
      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(1027,9),FL,FW,VR,FD(1027,9),RT(2,9),TWI
     *     ,xflt,yflt,zflt,delay  
     *    /IFLT/NF,NTW,IS,ngauss,ndip,nstr,iramp 
     *    /Hypo/xhypo,yhypo,zhypo 
     *    /OTST/IOT,IST  /FCQVW/SM,RAM,RR  /SOLY/SDEP,ODEP,LS,LO
     *    /DPSO /DEP1,DEP2,ZS,ZZO /outom/ioutom
     *    /INTG/VPMAX,VSMAX,RYMIN,CFINAL,N11,N12,M11
     *    /WVNUM/WVN(2001),DWN,WNR(3),NWN
c
      data gaux(1,1)/0.0d0/
      data gaux(1,2),gaux(2,2)
     * /-0.577350269189626d0, 0.577350269189626d0/
      data gaux(1,3),gaux(2,3),gaux(3,3)
     * /-0.774596669241483d0,0.d0, 0.774596669241483d0/
      data gaux(1,4),gaux(2,4),gaux(3,4),gaux(4,4)
     * /-0.861136311594053d0,-0.339981043584856d0
     * , 0.339981043584856d0, 0.861136311594053d0/
      data gaux(1,5),gaux(2,5),gaux(3,5),gaux(4,5),gaux(5,5)
     * /-0.906179845938664d0,-0.538469310105683d0,0.0d0
     * , 0.538469310105683d0, 0.906179845938664d0/
      data gaux(1,6),gaux(2,6),gaux(3,6),gaux(4,6),gaux(5,6),gaux(6,6)
     * /-0.932469514203152d0,-0.661209386466265d0,-0.238619186083197d0
     * , 0.238619186083197d0, 0.661209386466265d0, 0.932469514203152d0/ 
c
      data gw(1,1) /2.d0/
      data gw(1,2),gw(2,2) /1.d0,1.d0/
      data gw(1,3),gw(2,3),gw(3,3)
     * / 0.555555555555556d0, 0.888888888888889d0, 0.555555555555556d0/
      data gw(1,4),gw(2,4),gw(3,4),gw(4,4)
     * / 0.347854845137454d0, 0.652145154862546d0
     * , 0.652145154862546d0, 0.347854845137454d0/
      data gw(1,5),gw(2,5),gw(3,5),gw(4,5),gw(5,5)
     * / 0.236926885056189d0, 0.478628670499366d0, 0.568888888888889d0
     * , 0.478628670499366d0, 0.236926885056189d0/
      data gw(1,6),gw(2,6),gw(3,6),gw(4,6),gw(5,6),gw(6,6)
     * / 0.171324492379170d0, 0.360761573048139d0, 0.467913934572691d0
     * , 0.467913934572691d0, 0.360761573048139d0, 0.171324492379170d0/
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=1027
      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
      write(*,*)
      write(*,*)'Output Dynamic Integrands for checking? (no=0, yes=1)'
	read(*,*)IOT
C
      write(*,*)
      write(*,*)'Please choose the following combinations of integrands'
      write(*,*)' 0 : Regular (= Dynamic Integrand - Static Integrand)'
      write(*,*)' 1 : Check 1 (= Dynamic Integrand, only)'
      write(*,*)' 2 : Check 2 (= Static Integrand, only)'
	read(*,*)IST
c
c
      IF(IOT.eq.1)THEN
	  write(*,*)'The first observation point will be used!'
	  write(*,2000)xo(1),yo(1),zo(1)
2000    format(' (x,y,z)=',3f12.4)
	  write(*,*)'Enter the location of the source (x,y,z)'
	  read(*,*)xot,yot,zot
        OPEN (30,FILE='dint_d',STATUS='unknown')
C
        df0=dom/6.283185307d0
        f1=0.d0
        f2=df0*nom
        write(*,*)
        write(*,*)'Frequency Range (Hz) =',f1,f2
        write(10,*)'Frequency Range (Hz) =',f1,f2
        write(*,*)'Enter Target Frequency (Hz)'
        read(*,*)ofreq
        df=9999.d0
        do i=1,nom+1
          freq=df0*(i-1)
          dfreq=dabs(freq-ofreq)
          if(dfreq.lt.df)then
            ioutom=i
            df=dfreq
          end if
        end do
        write(*,*)'Actual Output Frequency (Hz)=',df0*(ioutom-1)
        write(10,*)'Actual Output Frequency (Hz)=',df0*(ioutom-1)
      END IF
C
c   * Green's function computations *
c
      im=DCMPLX(0.d0,1.d0)
      dcs=FW*dcos(dip)
      cst=dcos(str)
      sst=dsin(str)
c
      DO 100 IO=IOM,NOM+1
C
        if((iot.eq.1).and.(io.ne.ioutom))goto 100
c
	  if(io.eq.1)then
      write(*,*)
          do IOB=1,NOB
      write(*,*)'static, dynamic & total (io,iob)',io,iob
            GGX(IO,IOB)=GXS(IO,IOB)
            GGY(IO,IOB)=GYS(IO,IOB)
            GGZ(IO,IOB)=GZS(IO,IOB)
      write(*,*)'GGX',GGX(io,iob)
	    end do
	    goto 100
        end if
c
        OMR=DOM*DFLOAT(IO-1)
        OM=dcmplx(OMR,OMI)
	  freq=omr/6.2831853072d0
C
c     * Configure wavenumbers for numerical integrations *
c
c       * range 1: from near 0 to near om/Rayleigh(min) *
        WNR(2)=OMR/rymin*1.2
        WNR(1)=WNR(2)/n11/10000.d0
        CALL WVNUMV(MOM,IO,dreal(OM))
c
        NWN=2*(M11+N12)+1
        NC=2001
        IF(NWN.GT.NC)THEN
          WRITE(10,*)'FROM HKQVWd: NWN > 2001 !!!'
        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(*,*)
        write(*,1005)om,nwn
        write(10,*)
        write(10,1005)om,nwn
        write(10,1010)wnr(1),wnr(2),wnr(3)
        write(10,1020)dwn
 1005   format(' Omega, Num. of Wave-number integ. points =',2f10.4,I6)
 1010   format(' Boundaries for the Integration ranges:'/x,3e12.4)
 1020   format(' Integration Increments for Range 2:'/x,e12.4)
C
c     * the observation points *
c
        DO 200 IOB=1,NOB
C
          GGX(IO,IOB)=(0.D0,0.D0)
          GGY(IO,IOB)=(0.D0,0.D0)
          GGZ(IO,IOB)=(0.D0,0.D0)
c
          ZS=ZO(IOB)
          if(iob.ge.2)then
            if(zs.eq.zo(iob-1))goto 10
          end if

          DRU=0.0
          DRL=THK(1)
          DO 220 L=1,NL
            IF((ZS.GE.DRU).AND.(ZS.LT.DRL))LS=L
            IF((ZS.GE.DRU).AND.(L.EQ.NL))LS=L
            DRU=DRL
            DRL=DRL+THK(L+1)
 220     CONTINUE
C
          DEP1=0.0
          DO 230 L=1,LS-1
            DEP1=DEP1+THK(L)
 230      CONTINUE
          SDEP=ZS-DEP1
          DEP2=DEP1+THK(LS)
C
c       * Configure the up/down-going coefficient library *
c
          CALL CUDSH(OM)
          CALL CUDPSV(OM)
          CALL CUDSH0(dcmplx(0.d0,OMI))
          CALL CUDPSV0(dcmplx(0.d0,OMI))
c
c       * fault points and Green's function computation  *
c
 10       DO 300 idip=1,ndip 
 
            ngausL=1
            ngausW=1
            if(iot.eq.1)goto 15

            da=dcs*(idip-1)
            db=dcs*idip
            za=zflt-da*dtan(dip)*1.0000001d0
            zb=zflt-db*dtan(dip)*0.9999999d0
            if(zb.lt.0.d0)then
              write(*,*)'Fault Depth is less than 0!: Dip Number=',idip
              stop
            end if
            DRU=0.0
            DRL=THK(1)
            do L=1,NL
              IF((za.GE.DRU).AND.(za.LT.DRL))L1=L
              IF((za.GE.DRU).AND.(L.EQ.NL))  L1=L
              DRU=DRL
              DRL=DRL+THK(L+1)
            end do
            DRU=0.0
            DRL=THK(1)
            do L=1,NL
              IF((zb.GE.DRU).AND.(zb.LT.DRL))L2=L
              IF((zb.GE.DRU).AND.(L.EQ.NL))  L2=L
              DRU=DRL
              DRL=DRL+THK(L+1)
            end do
            if(L1.ne.L2)then
              write(*,*)' The L1 Layer is different from the L2 Layer!'
              write(*,*)' Layers 1 and 2 =',L1,L2
c             stop
            end if
c            
            Vmin=Vs(L1)
            if(Vs(L2).lt.Vmin)Vmin=Vs(L2)
            if(VR.lt.Vmin)Vmin=VR
c	      if(omr.lt.0.00001d0)then
c              wvlen=Vmin*6.283185307d0/0.000001d0
c	      else
              wvlen=Vmin*6.283185307d0/OMR
c	      end if
c
            gausW=FW*ngauss/wvlen
            ngausW=1
            if(gausW.gt.0.7d0)then
              ngausW=2
              if(gausW.gt.2.0d0)then
                ngausW=3
                if(gausW.gt.3.0d0)then
                  ngausW=4
                  if(gausW.gt.5.0d0)then
                    ngausW=5
                    if(gausW.gt.6.4d0)then
                      ngausW=6
                      if(gausW.gt.6.0d0)then
                  write(*,*)'ngausW > 6 !!!!; Use smaller sub-faults!'
                  write(*,*)'ngausW > 6 -> ngausW = 6'
                        ngausW=6
                      end if
                    end if
                  end if
                end if
              end if
            end if

            gausL=FL*ngauss/wvlen
            ngausL=1
            if(gausL.gt.0.7d0)then
              ngausL=2
              if(gausL.gt.2.0d0)then 
                ngausL=3
                if(gausL.gt.3.0d0)then
                  ngausL=4
                  if(gausL.gt.5.0d0)then
                    ngausL=5
                    if(gausL.gt.6.4d0)then
                      ngausL=6
                      if(gausL.gt.6.0d0)then
                  write(*,*)'ngausL > 6 !!!!; Use smaller sub-faults!'
                  write(*,*)'ngausL > 6 -> ngausL = 6'
                        ngausL=6
                      end if
                    end if
                  end if
                end if
              end if
            end if

   15     if(iob.eq.1)then 
            write(*,1200)idip,za,zb         
            write(*,1210)gausW,ngausW
            write(*,1220)gausL,ngausL
            write(10,1200)idip,za,zb         
            write(10,1210)gausW,ngausW
            write(10,1220)gausL,ngausL
          end if
 1200     format('  Elem Num.for W, and the edge depths =',i3,2f10.2)
 1210     format('  Gauss Points per Sub-Fault along W =',f9.2,' ->',i3)
 1220     format('  Gauss Points per Sub-Fault along L =',f9.2,' ->',i3)

          DO 300 istr=1,nstr
c
            IS=nstr*(idip-1)+istr
            sa=FL*(istr-1)
            sb=FL*istr
c
            do itw=1,ntw
              if(FD(IS,itw).ne.0.d0)goto 20
            end do
            goto 300
 20         continue
c
c           * Gaussian Quadrature (id*ig points) *
c
            do 400 id=1,ngausW
c
              zd=((1.d0-gaux(id,ngausW))*za
     *           +(1.d0+gaux(id,ngausW))*zb)/2.d0
              if(iot.eq.1)zd=zot
  
              DRU=0.0
              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
              zz=zd-zhypo

              LO=L1
              QPe=QP(LO)*(freq)**FP(LO)
              QSe=QS(LO)*(freq)**FS(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
              dd=((1.d0-gaux(id,ngausW))*da
     *           +(1.d0+gaux(id,ngausW))*db)/2.d0
     
	        do 450 ig=1,ngausL

                sg=((1.d0-gaux(ig,ngausL))*sa
     *             +(1.d0+gaux(ig,ngausL))*sb)/2.d0
                XF=xflt+sg*cst+dd*sst
                YF=yflt+sg*sst-dd*cst
	          if(iot.eq.1)then
	            XF=xot
	            YF=yot
	          end if
                xx=XF-xhypo
                yy=YF-yhypo
                xy=dsqrt(xx*xx+yy*yy+zz*zz)
                TD=delay+xy/vr            
                RX=XF-XO(IOB)
                RY=YF-YO(IOB)
                RR=DSQRT(RX*RX+RY*RY)
c
                EDT=CDEXP(im*OM*TD)
C
                DO J=0,2
                  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)
                  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)
    	      end do
c
c         * the wavenumber integration using Simpson's and Filon's fomulas *
                CALL QVWRZS0(OM,QXS,QZS,VXS,WXS,WZS)
                CALL QVWRZF0(OM,QXG,QZG,VXG,WXG,WZG)
c
c      write(*,*)'xf,yf,zd',xf,yf,zd
c      write(*,*)'rx,ry,rr',rx,ry,rr
c      write(*,*)'om',om
c      write(*,*)'QX0',QXS(0),QXG(0)
c      write(*,*)'QX1',QXS(1),QXG(1)
c      write(*,*)'QX2',QXS(2),QXG(2)
c      write(*,*)'VX1',VXS(1),VXG(1)
c      write(*,*)'VX2',VXS(2),VXG(1)
c      stop
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               * Time Windows *
c
                do 470 itw=1,ntw
c
                  CALL MO(itw,OM,en,SMO)
                  EDT0=CDEXP(im*OM*TWI*(itw-1))
                  GX=(0.D0,0.D0)
                  GY=(0.D0,0.D0)
                  GZ=(0.D0,0.D0)
                  DO I=1,3
                  DO J=1,3
                    GX=GX+EN(I,J)*UX(J,I)
                    GY=GY+EN(I,J)*UY(J,I)
                    GZ=GZ+EN(I,J)*UZ(J,I)
c      write(*,*)'i,j,UX',i,j,UX(J,I)
                  end do
                  end do
                  GX0=GX*SMO*EDT*EDT0*gw(id,ngausW)*gw(ig,ngausL)/4.d0
                  GY0=GY*SMO*EDT*EDT0*gw(id,ngausW)*gw(ig,ngausL)/4.d0
                  GZ0=GZ*SMO*EDT*EDT0*gw(id,ngausW)*gw(ig,ngausL)/4.d0
c
                  GGX(IO,IOB)=GGX(IO,IOB)+GX0
                  GGY(IO,IOB)=GGY(IO,IOB)+GY0
                  GGZ(IO,IOB)=GGZ(IO,IOB)+GZ0
c
                IF(IOT.eq.1)THEN
                  CLOSE(30)
      write(*,*)'GX0',GX0
      write(*,*)'GY0',GY0
      write(*,*)'GZ0',GZ0
      write(*,*)
      write(*,*)'dint.m plots the Green function integrands'
	stop
                END IF
c
  470           CONTINUE
c
  450         CONTINUE
c    
  400       CONTINUE
c
  300     CONTINUE
c
      write(*,*)'static, dynamic & total (io,iob)',io,iob
      write(*,*)'GXS',GXS(io,iob)
      write(*,*)'GGX',GGX(io,iob)
          GGX(IO,IOB)=GGX(IO,IOB)+GXS(IO,IOB)
          GGY(IO,IOB)=GGY(IO,IOB)+GYS(IO,IOB)
          GGZ(IO,IOB)=GGZ(IO,IOB)+GZS(IO,IOB)
      write(*,*)'GGX',GGX(io,iob)
c
  200   CONTINUE
C
  100 CONTINUE
C
      RETURN
      END
c
C
c
C *** INPUT OF Pole and Branch Points ***
C
      SUBROUTINE PHSBRC(MOM)
C
      IMPLICIT REAL*8(A-H,O-Z)
      dimension CR(531,0:201),CL(531,0:201),MDL(531),MDR(531)
      COMMON /FREQ/DOM,IOM,NOM,OMI /POLE/CS(531,201),NW(531)
     *       /STRC/DNS(20),VP(20),QP(20),FP(20)
     *                    ,VS(20),QS(20),FS(20),THK(20),NL
     *       /TMP1/CC(201)
C
C  ** CHECK ** 
C
      NC=531
      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=201
c
C  ** INPUT OF PHASE VELOCITIES FOR LOVE AND RAYLEIGH MODES **
C
      open(22,file='mdl_dat',status='old')
      open(23,file='mdr_dat',status='old')
      DO 100 IO=IOM,NOM+1
        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(20,file='phsl_dat',status='old')
      DO 110 I=0,MDL(NOM)
      DO 110 IO=IOM,NOM+1
        read(20,*)CL(IO,I)
  110 CONTINUE
      close(20)

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

      DO 200 IO=IOM,NOM+1

        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=201
        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
        write(10,1100)(io-1)*dom,nw(io)
        write(10,1110)
        do 320 i=1,NW(IO),5
          write(10,1120)CS(IO,i),CS(IO,i+1)
     *               ,CS(IO,i+2),CS(IO,i+3),CS(IO,i+4)
  320   continue
 1100   FORMAT('   Omega, Number of Phase Velocities=',f8.5,I5)
 1110   FORMAT( '   Phase Velocities (m/sec) =')
 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 *** MOMENT TENSOR FOR THE SOURCE ***
C
      SUBROUTINE MO(itw,OM,en,SMO)
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION EN(3,3)
      COMPLEX*16 SMO,CPX,SM,RAM,STT,OM,TTK,TK1,TK2
      COMMON/STRC/DNS(20),VP(20),QP(20),FP(20)
     *                   ,VS(20),QS(20),FS(20),THK(20),NL
     *    /FLT/STR,DIP,RAK(1027,9),FL,FW,VR,FD(1027,9),RT(2,9),TWI
     *     ,xflt,yflt,zflt,delay
     *    /IFLT/NF,NTW,IS,ngauss,ndip,nstr,iramp
     *    /FCQVW/SM,RAM,RR 
C
C *** CHECK OF DIMENSION ***
C
      NC=20
      IF(NL.GE.NC)THEN
        WRITE(*,*)' FROM ENMO:NL>NC !!!!!!!!!!'
        WRITE(*,*)' NL,NC=',NL,NC
        STOP
      END IF
C
      NC=1027
      IF(NF.GE.NC)THEN
        WRITE(*,*)' FROM ENMO:NF>NC !!!!!!!!!!'
        WRITE(*,*)' NF,NC=',NF,NC
        STOP
      END IF
C
      CPX=(0.D0,1.D0)
c
C *** EN ***
C
      CST= DCOS(STR)
      SST= DSIN(STR)
C
      CDP= DCOS(DIP)
      SDP= DSIN(DIP)
C
      CRK= DCOS(RAK(is,itw))
      SRK= DSIN(RAK(is,itw))
C
      XN=-SDP*SST
      YN= SDP*CST
      ZN=-CDP
C
      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
C *** SMO ***
C
      SMO= FD(IS,itw)*SM*FL*FW
C
      IF(dreal(OM).EQ.0.D0)THEN      
        STT= (1.D0,0.D0)      
      ELSE          
        TTK= OM*(RT(1,itw)+RT(2,itw))/2.D0
        if(iramp.eq.0)then
          STT= CDSIN(TTK)*CDEXP(CPX*TTK)/TTK
        else
          if(iramp.eq.1)then
            TK1= OM*RT(1,itw)/2.D0
            TK2= OM*RT(2,itw)/2.D0
            STT= CDSIN(TK1)/TK1-CDSIN(TK2)/TK2*CDEXP(CPX*TTK)
            STT= CPX*STT*CDEXP(CPX*TK1)/TTK
          else
          T0=RT(1,itw)+RT(2,itw)
          STT=1-CPX*OM*T0
          STT=1/(STT*STT)
          end if
        end if                  
      END IF
C
c      SMO=CPX*SMO*STT/OM
c      write(*,*)'stt',stt
      SMO=SMO*STT
C
      RETURN
      END
c
c
C *** CONSTRUCT POINTS FOR WAVENUMBER INTEGRATION ***
C
      SUBROUTINE WVNUMV(MOM,IO,OM)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON /INTG/VPMAX,VSMAX,RYMIN,CFINAL,N11,N12,M11
     *       /WVNUM/WVN(2001),DWN,WNR(3),NWN
     *       /POLE/CS(531,201),NW(531) /TMP/WVT(1001)
C
C  ** CHECK ** 
C
      NC=531
      IF(MOM.GT.NC)THEN
         WRITE(10,*)'FROM WVNUM: MOM > NC !!'
         STOP
      END IF
C
      NC=201
      IF(NW(IO).GT.NC)THEN
         WRITE(10,*)'FROM WVNUM: NW(IO) > NC !!'
         STOP
      END IF
C
      NC=1001
      IF((N11+1+NW(IO)*9).GT.NC)THEN
         WRITE(10,*)'FROM WVNUM: N11+1+NW(IO)*9 > 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(1001) 
C
      DO 100 I=1,MW
        WA=WVT(I)
        WB=WVT(I+1)
        IF((W.GT.WA).AND.(W.LT.WB))THEN
          DO 105 K=MW+1,I+1,-1
            WVT(K+1)=WVT(K)
 105      CONTINUE
          WVT(I+1)=W
          MW=MW+1
          RETURN
        END IF
 100  CONTINUE
c
      RETURN
      END
