c
C *** This code is modified from grflt941.f to skip zero slip sub-faults
c    (Oct, 1998)    
c
C *** This code is modified from grflt951.f to include the exp-type source
c     time function (Aug, 2000)    
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, and 1997)
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 * Parameters : 
c     MOM : Number of Omega
c     MOB : NUmber of Observation Points
c
      PROGRAM grflt95p
C
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(MOM=531,MOB=51)
      COMPLEX*16 GX(MOM,MOB),GY(MOM,MOB),GZ(MOM,MOB)
C
      WRITE(*,*)
      WRITE(*,*)
      WRITE(*,*)'"grflt95p.f" constructs Green functions'
      WRITE(*,*)' using the Phinney method'
     *'           (Originally made in 1992, and modified in 2000)'
      write(*,*)
      WRITE(*,*)'"grflt95p.in" is the input data'
      WRITE(*,*)
     *' You must run "phs1.f" before running this code !!'
      WRITE(*,*)' Please Check "fort.10" for troubles.'
      WRITE(*,*)
     *' For questions, e-mail to "hisada@cc.kogakuin.ac.jp".'
      WRITE(*,*)' This code is open for the academic use.'
C
      CALL INDAT(MOM,MOB)
      CALL INPHS(MOM)
      CALL HKQVW(MOM,MOB,GX,GY,GZ)
      CALL OTDAT(MOM,MOB,GX,GY,GZ)
C
      STOP
      END
C
c
C *** DATA INPUT ***
C
      SUBROUTINE INDAT(MOM,MOB)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VP,VS
C      character infile*20
      COMMON /FREQ/DOM,IOM,NOM,OMI /OBS/XO(51),YO(51),ZO(51),NOB
     *    /STRC/DNS(20),VP(20),VS(20),THK(20),NL
     *    /INTG/VPMAX,VSMAX,RYMIN,CFINAL,N11,N12,N13,N14,M13
     *    /FLT/STR,DIP,RAK(203,6),FL,FW,VR,FD(203,6),RT(2,6),TWI
     *     ,XF(203,8,8),YF(203,8,8),TD(203,8,8),DP(203,8)
     *    /IFLT/LF(203,8),NF,NTW,IS,ngauss
     *    /Hypo/xhypo,yhypo,zhypo 
     *    /SIGN/ISIGN /OTST/IOT,IST /ramp/iramp,NSTR,NDIP
      dimension gaux(8,8)
      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
C *** PARAMETERS IN COMMOM STATEMENT ***
C
      NML =20
      NMOB=51
      NFLT=203
      IF(MOB.LT.NMOB)NMOB=MOB
C
C *** OMEGA DATA ***
C
      write(*,*)
C      write(*,*)'Enter File name for "grflt95e.in"'
C      read(*,'(a20)')infile
C      OPEN (UNIT=7,FILE=infile,STATUS='OLD')
      OPEN (UNIT=7,FILE='grflt95p.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.GT.MOM)THEN
        WRITE(10,*)' NOM>MOM !!!!!!!!!!!!!'
        STOP
      END IF
C
C *** MEDIUM DATA ***
C
      READ(7,1200)NL
      DO 100 IL=1,NL
        READ(7,1210)DNS(IL),VPP,QP,VSS,QS,THK(IL)
        VP(IL)=VPP*DCMPLX(1.D0,-0.5D0/QP)
        VS(IL)=VSS*DCMPLX(1.D0,-0.5D0/QS)
  100 CONTINUE
 1200 FORMAT(/I10)
 1210 FORMAT(6E8.0)
C
      WRITE(10,2200)NL
      DO 140 IL=1,NL
        WRITE(10,2210)IL,DNS(IL),VP(IL),VS(IL),THK(IL)
  140 CONTINUE
 2200 FORMAT(//' Total Number of Layers =',I5/
     *       /' Layer, Density , Vp(real), Vp(imag)'
     *        ', Vs(real), Vs(imag),Thickness(m)')
 2210 FORMAT(I5,5F10.3,F12.3)
C
      IF(NL.GT.NML)THEN
        WRITE(10,*)' NL>NML !!!!!!!!!!!!!'
        STOP
      END IF
C
C *** MAX Vs and Vp and MIN Rayleigh Phase Velocity ***
C
      VSMAX=dreal(VS(NL))
      VPMAX=dreal(VP(NL))
      VSMIN=DREAL(VS(1))
      MINVS=1
      DO 160 IL=1,NL
        IF(dreal(VS(IL)).GT.VSMAX)VSMAX=dreal(VS(IL))
        IF(dreal(VP(IL)).GT.VPMAX)VPMAX=dreal(VP(IL))
        IF(dreal(VS(IL)).LT.VSMIN)THEN
          VSMIN=dreal(VS(IL))
          MINVS=IL
        END IF
  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,depth
c      write(*,*)'xflt,yflt,depth',xflt,yflt,depth
      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.6)then
         write(*,*)'ntw > 6 !!!'
         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
      dcs=fw*dcos(dip)
      cst=dcos(str)
      sst=dsin(str)
      iflt=0
							
c   * Dip direction *      
      do 200 i=1,ndip

        da=dcs*(i-1)
        db=dcs*i
        za=depth-da*dtan(dip)
        zb=depth-db*dtan(dip)
        if((za.lt.0.d0).or.(zb.lt.0.d0))then
          write(*,*)'Fault Depth is less than 0!: Dip Number=',i
          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(*,*)'Dip number, L1, and L2 =',i,L1,L2
c          stop
        end if

c     * Strike direction *
        do 220 j=1,nstr   
   
          iflt=iflt+1

          sa=fl*(j-1)
          sb=fl*j
          do 240 id=1,ngauss

            zd=((1.d0-gaux(id,ngauss))*za
     *         +(1.d0+gaux(id,ngauss))*zb)/2.d0

            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

            LF(iflt,id)=L1
            dp(iflt,id)=zd
            zz=zd-zhypo
c
            dd=((1.d0-gaux(id,ngauss))*da
     *         +(1.d0+gaux(id,ngauss))*db)/2.d0

            do ig=1,ngauss
              sg=((1.d0-gaux(ig,ngauss))*sa
     *           +(1.d0+gaux(ig,ngauss))*sb)/2.d0
              xf(iflt,id,ig)=xflt+sg*cst+dd*sst
              yf(iflt,id,ig)=yflt+sg*sst-dd*cst
	      xx=xf(iflt,id,ig)-xhypo
	      yy=yf(iflt,id,ig)-yhypo
              dis=dsqrt(xx*xx+yy*yy+zz*zz)
              td(iflt,id,ig)=delay+dis/vr            
            end do

  240     continue

  220   continue

  200 continue
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)depth,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( '   Num. of Gauss Points =',i5
     *       ,',  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)
      do id=1,ngauss
        write(10,2360)id,DP(II,id),LF(II,id)
        write(10,2365)(XF(II,id,ig),ig=1,ngauss)
        write(10,2370)(YF(II,id,ig),ig=1,ngauss)
        write(10,2375)(TD(II,id,ig),ig=1,ngauss)
      end do
  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) 
 2360 format( '   Gauss Point Number along Dip=',i5
     *        ',  Depth (m) =',f8.2,',   Layer=',i5)
 2365 FORMAT( '   Source X-Axis along Str (m) =',8f8.1)
 2370 FORMAT( '   Source Y-Axis along Str (m) =',8f8.1)
 2375 FORMAT( '   Total Time Delay      (sec) =',8f8.3)
c
C *** NUMERICAL INTEGRAL DATA ***
C
      READ(7,1400)N11
      READ(7,1410)N12
      READ(7,1410)N13
      READ(7,1410)N14
      CFINAL=rymin/10.0
      write(*,*)'CFINAL=rymin/10.0'
 1400 FORMAT(/I10)
 1410 FORMAT(I10)
 1420 FORMAT(F10.0)
C
      WRITE(10,2400)
      WRITE(10,2410)N11
      WRITE(10,2420)N12
      WRITE(10,2430)N13
      WRITE(10,2440)N14
      WRITE(10,2450)CFINAL
      WRITE(10,2460)2*(N11+N12+N13+N14)+1
 2400 FORMAT(//' Number of Points for Numerical Integrations:')
 2410 FORMAT( /'   From 0               to 0.8*om/Vp(max)   =',I5)
 2420 FORMAT(  '   From 0.8*om/Vp(max)  to 0.95*om/Vs(max)  =',I5)
 2430 FORMAT(  '   From 0.95*om/Vs(max) to 1.05*om/Ray(min) =',I5)
 2440 FORMAT(  '   From 1.05*om/Ry(min) to om/C(final)      =',I5)
 2450 FORMAT(  '   C(final = Minimun Rayleigh Phase Vel./10)=',F10.3)
 2460 FORMAT(  '   Initial Number of Points: 2(N1+N2+N3+N4)+1 =',I5)
C
C *** DATA FOR THE ASYMPTOTIC ***
C
      READ(7,1600)IST
      IF(IST.EQ.1)THEN
        WRITE(10,2600)
      ELSE
        WRITE(10,2610)
      END IF
 1600 FORMAT(/I10)
 2600 FORMAT(//' THE ASYMPTOTIC METHOD IS USED')
 2610 FORMAT(//' THE ASYMPTOTIC METHOD IS NOT USED')
C
C *** CHANGE OF SIGNS FOR IMAGINARY PARTS ***
C
      READ(7,1700)ISIGN
      IF(ISIGN.EQ.1)WRITE(10,2700)
      IF(ISIGN.EQ.0)WRITE(10,2710)
 1700 FORMAT(/I10)
 2700 FORMAT(/' SIGNS OF IMAGINARY PARTS FOR RESULTS ARE CHANGED')
 2710 FORMAT(/' SIGNS OF IMAGINARY PARTS FOR RESULTS ARE NOT CHANGED')
C
C *** OP DATA ***
C
      READ(7,1500)NOB
      WRITE(10,2500)NOB
 1500 FORMAT(/I10)
 2500 FORMAT(//' TOTAL NUMBER OF OBSERVATION POINTS =',I6/)
C
      IF(NOB.GT.NMOB)THEN
        WRITE(10,*)' NOB>NMOB !!!!!!!!!!!!!'
        STOP
      END IF
C
      DO 300 I=1,NOB
        READ(7,1520)XO(I),YO(I),ZO(I)
        WRITE(10,2520)I,XO(I),YO(I),ZO(I)
  300 CONTINUE
 1520 FORMAT(3F10.0)
 2520 FORMAT('   OB NUM.=',I5,',     (X,Y,Z) =',3F12.4)
C
      CLOSE(UNIT=7)
C
c *** output of fault and observation data ***
c
      open(9,file='fault_dat',status='unknown')
        write(9,9500)XHYPO,YHYPO,ZHYPO
        write(9,9500)xflt,yflt,depth
        write(9,9500)STR,DIP,RAK(1,1)
        write(9,9500)slength,dwidth,VR
        write(9,9510)NSTR,NDIP,ngauss
        do 700 iflt=1,nf
        do 700 id=1,ngauss
        do 700 ig=1,ngauss
        write(9,9500)xf(iflt,id,ig),yf(iflt,id,ig),dp(iflt,id)
  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
      IF((NOM-IOM).NE.0)RETURN
C
      WRITE(*,7000)
      READ(*,*)IOT
 7000 FORMAT(/' Output Integrands?: Enter No=0 or Yes=1: ')
C
      if(iot.eq.1)then
        write(*,*)
        write(*,*)'  You can plot the integrands using "integ.m",'
        write(*,*)'  "integr.m", and "integz.m" for "MATLAB"'
c
        open(9,file='in_dat',status='unknown')
          write(9,9000)dom*iom
          write(9,9100)nl
          do 500 i=1,nl
            write(9,9000)dreal(vp(i))
            write(9,9000)dreal(vs(i))
 500      continue
        close(9)
 9000   format(e15.7)
 9100   format(i5)
c
      end if

      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 VP,VS,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),VS(20),THK(20),NL
     *    /FLT/STR,DIP,RAK(203,6),FL,FW,VR,FD(203,6),RT(2,6),TWI
     *     ,XF(203,8,8),YF(203,8,8),TD(203,8,8),DP(203,8) 
     *    /IFLT/LF(203,8),NF,NTW,IS,ngauss 
     *    /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=203
      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
          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
            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
            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,OMI
        WRITE(8,3510)NOB
        DO 800 IOB=1,NOB
          DO 810 K=IOM,NOM
            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 *** INPUT OF PHASE VELOCITIES AND MAKE POLE DATA ***
C
      SUBROUTINE INPHS(MOM)
C
      IMPLICIT REAL*8(A-H,O-Z)
      dimension CR(531,0:201),CL(531,0:201),MDL(531),MDR(531)
      COMPLEX*16 VP,VS
      COMMON /FREQ/DOM,IOM,NOM,OMI /POLE/CS(531,201),NW(531)
     *       /STRC/DNS(20),VP(20),VS(20),THK(20),NL
C
C  ** CHECK ** 
C
      NC=531
      IF(MOM.GT.NC)THEN
         WRITE(10,*)'FROM INPHS: MOM>NC !!'
         STOP
      END IF
C
      NC=20
      IF(NL.GT.NC)THEN
         WRITE(10,*)'FROM INPHS: NL>NC !!'
         STOP
      END IF

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

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

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

      DO 200 IO=IOM,NOM

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

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

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


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

  200 continue
C
      RETURN
      END
C
C
C *** COMPUTATION OF GREEN FUNCTIONS USING THE R/T METHOD ***
c *** For large wavenumber * distance (k*r > 10), the Filon's integration 
c     formula is used instead of the Simpson's.
C
      SUBROUTINE HKQVW(MOM,MOB,GGX,GGY,GGZ)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VP,VS,SM,RAM,SM1,XX1,SM2,XX2,SM3,XX3,RM1,RM2,RM3
     *      ,GGX(MOM,MOB),GGY(MOM,MOB),GGZ(MOM,MOB),VS2,VP2
     *      ,QXG(0:2),QZG(0:2),VXG(0:2),WXG(0:2),WZG(0:2)
     *      ,QXS(0:2),QZS(0:2),VXS(0:2),WXS(0:2),WZS(0:2)
     *      ,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),OM,im
      DIMENSION EN(3,3),gw(8,8)
      COMMON /FREQ/DOM,IOM,NOM,OMI /OBS/XO(51),YO(51),ZO(51),NOB
     *    /STRC/DNS(20),VP(20),VS(20),THK(20),NL
     *    /FLT/STR,DIP,RAK(203,6),FL,FW,VR,FD(203,6),RT(2,6),TWI
     *     ,XF(203,8,8),YF(203,8,8),TD(203,8,8),DP(203,8)  
     *    /IFLT/LF(203,8),NF,NTW,IS,ngauss 
     *    /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
     *    /INTG/VPMAX,VSMAX,RYMIN,CFINAL,N11,N12,N13,N14,M13
     *    /WVNUM/WVN(2001),DWN(4),WNR(5),NWN
      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=203
      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   * Green's function computations *
c
      im=DCMPLX(0.d0,1.d0)
      DO 100 IO=IOM,NOM
C
        OMR=DOM*DFLOAT(IO)
        OM=dcmplx(OMR,OMI)
C
c     * Configure wavenumbers for numerical integrations *
c
c       * range 1 (from near 0 to near om/Vp(max) *
        WNR(2)=OMR/vpmax*0.8
        WNR(1)=WNR(2)/n11/10000.d0
        DWN(1)=(WNR(2)-WNR(1))/(2.D0*N11)
        DO 120 K=1,2*N11
          WVN(K)=WNR(1)+DWN(1)*(K-1)
 120    CONTINUE

c       * range 2 (from near om/Vp(max) to near om/Vs(max) *
        WNR(3)=OMR/vsmax*0.95
        DWN(2)=(WNR(3)-WNR(2))/(2.D0*N12)
        DO 130 K=1,2*N12
          WVN(2*N11+K)=WNR(2)+DWN(2)*(K-1)
 130    CONTINUE
C
c       * range 3 (from near om/Vs(max) to near om/Rayleigh(min) *
        WNR(4)=OMR/rymin*1.05 
        CALL WVNUMV(MOM,IO,OM)
c
        NWN=2*(N11+N12+M13+N14)+1
        NC=2001
        IF(NWN.GT.NC)THEN
          WRITE(10,*)'FROM HKQVW: NWN > 2001 !!!'
        STOP
        END IF
C
c       * range 4 (from near om/Rayleigh(min) to near om/C(final) *
        WNR(5)=OMR/cfinal
        DWN(4)=(WNR(5)-WNR(4))/(2.D0*N14)
        DO 150 K=1,2*N14+1
          WVN(2*(N11+N12+M13)+K)=WNR(4)+DWN(4)*(K-1)
 150    CONTINUE
C
        write(*,1005)om,nwn
        write(10,*)
        write(10,1005)om,nwn
        write(10,1010)wnr(1),wnr(2),wnr(3),wnr(4),wnr(5)
        write(10,1020)dwn(1),dwn(2),dwn(3),dwn(4)
 1005   format('    Omega, Num. of integration points =',2f10.4,I6)
 1010   format('    Boundaries for the Integration ranges:'/x,5e12.4)
 1020   format('    Integration Increments for the ranges:'/x,4e12.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
          VS2= VS(LS)*VS(LS)
          SM2= DNS(LS)*VS2
          VP2= VP(LS)*VP(LS)
          RM2= DNS(LS)*VP2-2.0*SM2
          XX2 = (1.0+VS2/VP2)/(1.0-VS2/VP2)
          IF(LS.GE.2)THEN
            VS2= VS(LS-1)*VS(LS-1)
            SM1= DNS(LS-1)*VS2
            VP2= VP(LS-1)*VP(LS-1)
            RM1= DNS(LS-1)*VP2-2.0*SM1
            XX1 = (1.0+VS2/VP2)/(1.0-VS2/VP2)
          END IF
          IF(LS.NE.NL)THEN
            VS2= VS(LS+1)*VS(LS+1)
            SM3= DNS(LS+1)*VS2
            VP2= VP(LS+1)*VP(LS+1)
            RM3= DNS(LS+1)*VP2-2.0*SM3
            XX3 = (1.0+VS2/VP2)/(1.0-VS2/VP2)
          END IF
C
c       * 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
            do itw=1,ntw
              if(FD(IS,itw).ne.0.d0)goto 20
            end do
            goto 300
c 20     write(*,*)'Not 0 Slip:IS,itw,FD(IS,itw)',IS,itw,FD(IS,itw)
 20         continue
c
c           * Gaussian Quadrature (id*ig points) *
c
            do 400 id=1,ngauss
c
              LO=LF(IS,id)
              SM = DNS(LO)*VS(LO)*VS(LO)
              RAM= DNS(LO)*VP(LO)*VP(LO)-2.D0*SM
C
              ZZO=DP(IS,id)
              ODEP=ZZO
              DO L=2,LO
                ODEP=ODEP-THK(L-1)
              end do
C
	      do 450 ig=1,ngauss

                RX=XF(IS,id,ig)-XO(IOB)
                RY=YF(IS,id,ig)-YO(IOB)
                RR=DSQRT(RX*RX+RY*RY)
c
                EDT=CDEXP(im*OM*TD(IS,id,ig))
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               * Time Windows *
c
                do itw=1,ntw
c
                  CALL MO(itw,OM,en,SMO)
                  EDT0=CDEXP(im*OM*TWI*(itw-1))
C
                  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)
                  end do
		  end do

                  GX0=GX*SMO*EDT*EDT0*gw(id,ngauss)*gw(ig,ngauss)/4.d0
                  GY0=GY*SMO*EDT*EDT0*gw(id,ngauss)*gw(ig,ngauss)/4.d0
                  GZ0=GZ*SMO*EDT*EDT0*gw(id,ngauss)*gw(ig,ngauss)/4.d0
                  GGX(IO,IOB)=GGX(IO,IOB)+GX0
                  GGY(IO,IOB)=GGY(IO,IOB)+GY0
                  GGZ(IO,IOB)=GGZ(IO,IOB)+GZ0
c
		end do
c
  450         CONTINUE
c    
  400       CONTINUE
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 *** 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 VS,VP,SMO,CPX,SM,RAM,STT,OM,TTK,TK1,TK2
      COMMON/STRC/DNS(20),VP(20),VS(20),THK(20),NL
     *    /FLT/STR,DIP,RAK(203,6),FL,FW,VR,FD(203,6),RT(2,6),TWI
     *     ,XF(203,8,8),YF(203,8,8),TD(203,8,8),DP(203,8)
     *    /IFLT/LF(203,8),NF,NTW,IS,ngauss
     *    /FCQVW/SM,RAM,RR /ramp/iramp,NSTR,NDIP
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=203
      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
      SMO=CPX*SMO*STT/OM
C
      RETURN
      END
c
c
C *** CONSTRUCT POINTS FOR WAVENUMBER INTEGRATION ***
C
      SUBROUTINE WVNUMV(MOM,IO,OM)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON /INTG/VPMAX,VSMAX,RYMIN,CFINAL,N11,N12,N13,N14,M13
     *       /WVNUM/WVN(2001),DWN(4),WNR(5),NWN
     *       /POLE/CS(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=202
      IF(NW(IO).GT.NC)THEN
         WRITE(10,*)'FROM WVNUM: NW(IO) > NC !!'
         STOP
      END IF
C
      NC=1001
      IF((N13+1+NW(IO)*9).GT.NC)THEN
         WRITE(10,*)'FROM WVNUM: N13+1+NW(IO)*9 > NC !!'
         STOP
      END IF
C
C  ** POINTS **
C
      DWN(3)=(WNR(4)-WNR(3))/N13
      DO 100 K=1,N13+1
        WVT(K)=WNR(3)+DWN(3)*(K-1)
 100  CONTINUE

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

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

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

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

      M13=MW

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

      RETURN
      END

