C     
c *** This code computes Green's functions of the homogeneous full-space
c     from a fault model with sub-faults, considering 
c     1) the smooth rupture velocity using the many integration point
c        up to 441(=21*21) 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
c *** The rupture velocity can vary with the Gaussian distributed 
c     deviate in order to generate high-frequency waves:
c     1) Sub faults (NSTR*NDIP) are divided into sub-fault 1 (nsub1*nsub1) 
c        for variable Vr
c     2) Sub-fault 1 are divided into sub-fault 2 (nsub2*nsub2) for integrations
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 1998)
c * Parameters : 
c     MOM : Number of Omega
c     MOB : NUmber of Observation Points
c
      PROGRAM grfltHF
C
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(MOM=551,MOB=51,MFLT=7,MSUB1=201)
      COMPLEX*16 GX(MOM,MOB),GY(MOM,MOB),GZ(MOM,MOB)
      DIMENSION RAK(MFLT,MFLT,6),FD(MFLT,MFLT,6)
     *         ,XF(MFLT,MFLT,MSUB1,MSUB1)
     *	       ,YF(MFLT,MFLT,MSUB1,MSUB1)
     *	       ,ZF(MFLT,MFLT,MSUB1,MSUB1)
     *	       ,TD(MFLT,MFLT,MSUB1,MSUB1)
C
      WRITE(*,*)
      WRITE(*,*)'"grfltHF5.f" constructs Green functions'
      WRITE(*,*)
     *'           (Originally made in 1992, and modified in 1998)'
      write(*,*)
      WRITE(*,*)'"grfltHF5.in" is the input data'
      write(*,*)
      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,XF,YF,ZF,TD,RAK,FD,MFLT,MSUB1)
      CALL HKQVW(MOM,MOB,GX,GY,GZ,XF,YF,ZF,TD,RAK,FD
     *          ,MFLT,MSUB1)
      CALL OTDAT(MOM,MOB,GX,GY,GZ)
C
      STOP
      END
C
c
C *** DATA INPUT ***
C
      SUBROUTINE INDAT(MOM,MOB,XF,YF,ZF,TD,RAK,FD,MFLT,MSUB1)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VP,VS
      DIMENSION RAK(MFLT,MFLT,6),FD(MFLT,MFLT,6)
     *         ,XF(MFLT,MFLT,MSUB1,MSUB1)
     *	       ,YF(MFLT,MFLT,MSUB1,MSUB1)
     *	       ,ZF(MFLT,MFLT,MSUB1,MSUB1)
     *	       ,TD(MFLT,MFLT,MSUB1,MSUB1)
      COMMON /FREQ/DOM,IOM,NOM /OBS/XO(51),YO(51),ZO(51),NOB
     *    /STRC/DNS,VP,VS
     *    /FLT/STR,DIP,FL,FW,RT(2,6),TWI
     *    /IFLT/NSTR,NDIP,nsub1,NTW,iramp
     *    /Hypo/xhypo,yhypo,zhypo 
     *    /SIGN/ISIGN
c
C *** PARAMETERS IN COMMOM STATEMENT ***
C
      NMOB=51
      IF(MOB.LT.NMOB)NMOB=MOB
C
C *** OMEGA DATA ***
C
      write(*,*)
      OPEN (UNIT=7,FILE='grfltHFs.in',STATUS='OLD')
C
      read(7,1100)dt,nt
      read(7,1110)period
 1100 format(/f8.0,i8)
 1110 format(f8.0)
      duration=dt*nt
      dom=datan(1.d0)*8.d0/duration
      iom=1
      nom=duration/period

      WRITE(10,*)'FROM grfltHF3.in:'
      WRITE(10,*)
      WRITE(10,1130)nt,dt,period
      WRITE(10,1140)NOM,DOM
 1130 format(/' Num. of Time, Delta Time and Min. Period (s)=',i5,2f8.4)
 1140 format( ' Num. of Omega, and Delta Omega (rad/sec) =',i5,f10.7)
C
      IF(NOM.GT.MOM)THEN
        WRITE(10,*)' NOM>MOM !!!!!!!!!!!!!'
        STOP
      END IF
C
C *** MEDIUM DATA ***
C
      READ(7,1200)DNS,VPP,QP,VSS,QS
      VP=VPP*DCMPLX(1.D0,-0.5D0/QP)
      VS=VSS*DCMPLX(1.D0,-0.5D0/QS)
 1200 FORMAT(/6E8.0)
C
      WRITE(10,2200)
      WRITE(10,2210)DNS,VP,VS
 2200 FORMAT(/' Density , Vp(real), Vp(imag)'
     *        ', Vs(real), Vs(imag),Thickness(m)')
 2210 FORMAT(5F10.3,F12.3)
C
C *** FAULT DATA ***
C
      READ(7,1300)SLENGTH,DWIDTH
c      write(*,*)'SLENGTH,DWIDTH',SLENGTH,DWIDTH
      READ(7,1310)NSTR,NDIP
      READ(7,1310)nsub1
      write(*,*)'NSTR,NDIP',NSTR,NDIP
      write(*,*)'nsub1',nsub1
      IF(NSTR.GT.MFLT)THEN
        WRITE(10,*)' NSTR.GT.MFLT !!!!!!'
        WRITE(*,*)' NSTR.GT.MFLT !!!!!!'
        STOP
      END IF
      IF(NDIP.GT.MFLT)THEN
        WRITE(10,*)' NDIP.GT.MFLT !!!!!!'
        WRITE(*,*)' NDIP.GT.MFLT !!!!!!'
        STOP
      END IF
      if(nsub1.gt.msub1)then
        write(10,*)' nsub1 > msub1 !!!'
        write(*,*)' nsub1 > msub1 !!!'
        stop
      end if
      READ(7,1320)DELAY,IRAMP
c      write(*,*)'DELAY,IRAMP',DELAY,IRAMP
      READ(7,1330)SHYPO,DHYPO
c      write(*,*)'SHYPO,DHYPO',SHYPO,DHYPO
      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,1332)VR
      READ(7,1335)NTW,TWI
 1300 FORMAT(/2f10.0)
 1310 FORMAT(2I5)
 1320 FORMAT(F10.0,I10)
 1325 FORMAT(i10)
 1330 format(3f10.0)
 1332 format(2f10.0,i5)
 1335 format(i8,f8.0)
c
c      write(*,*)'NTW,TWI',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 i=1,ndip
      do j=1,nstr
        READ(7,*)(FD(i,j,k),k=1,NTW)
c        write(*,*)'i,j,FD(i,j,k)',i,j,(FD(i,j,k),k=1,NTW)
      end do
      end do
      read(7,*)
      do i=1,ndip
      do j=1,nstr
        READ(7,*)(rak(i,j,k),k=1,NTW)
c        write(*,*)'i,j,rak(i,j,k)',i,j,(rak(i,j,k),k=1,NTW)
        do k=1,ntw
          rak(i,j,k)=rak(i,j,k)*0.01745329252d0
        end do
      end do
      end do
 1340 format(2f8.0)
c 1345 format(6f6.0)
c
c   * Generation of Rupture Velocities at Sub-faults 
c							
c   * Dip direction *      
C
      fl=slength/nstr
      fw=dwidth/ndip

      cdp=dcos(dip)
      sdp=dsin(dip)
      cst=dcos(str)
      sst=dsin(str)
      xhypo=xflt+shypo*cst+dhypo*cdp*sst
      yhypo=yflt+shypo*sst-dhypo*cdp*cst
      zhypo=zflt-dhypo*sdp
c
      open(1,file='vel_dat',status='unknown')
      do 200 i=1,ndip

        da=fw*(i-1)
        dd1=fw/nsub1

c     * Strike direction *
        do 210 j=1,nstr   
   
          sa=fl*(j-1)
          ds1=fl/nsub1

c       * Sub-Fault 1 for Integration
          do 220 id=1,nsub1
          do 230 ig=1,nsub1
c
            da1=da+dd1*(id-1)+dd1/2.d0
            sa1=sa+ds1*(ig-1)+ds1/2.d0
c
            xx=da1-DHYPO
            yy=sa1-SHYPO
            dis=dsqrt(xx*xx+yy*yy)
            ta1=dis/Vr
            xa1=xflt+sa1*cst+da1*cdp*sst
            ya1=yflt+sa1*sst-da1*cdp*cst
            za1=zflt-da1*sdp
            if(za1.lt.0.d0)then
               write(*,*)'za1 < 0 !!'
               stop
            end if
            xf(i,j,id,ig)=xa1
            yf(i,j,id,ig)=ya1
            ZF(i,j,id,ig)=za1
            td(i,j,id,ig)=delay+ta1
c
  230     continue
  220     continue
c
  210   continue
  200 continue
C
      WRITE(10,2300)IRAMP,nstr,ndip
      WRITE(10,2305)STR,DIP
      WRITE(10,2315)slength,dwidth
      WRITE(10,2320)FL,FW
      write(10,2322)shypo,dhypo
      write(10,2324)xhypo,yhypo,zhypo
      write(10,2325)xflt,yflt,zflt
      WRITE(10,2330)VRavg, Vrstd
      WRITE(10,2335)ntw,twi
      WRITE(10,2340)nsub1
      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)' 
     *      //'   Numbers of Sub_Fault 0 (NSTR,NDIP)=',2I5)
 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)
 2322 FORMAT( '   Str_Hypo.    (m)=',F12.3,',  Dip_Hypo     (m)=',F12.3)
 2324 FORMAT( '   X_Hypocenter (m)=',F12.3,',  Y_Hypocenter (m)=',F12.3
     *       /'   Z_Hypocenter (m)=',F12.3)
 2325 FORMAT( '   X Fault Edge (m)=',F12.3,',  Y Fault Edge (m)=',F12.3
     *       /'   Z Fault Edge (m)=',F12.3)
 2330 FORMAT( '   Average Vr (m/s)=',F12.3,',  Std. Dev.Vr(m/s)=',F12.3)
 2335 FORMAT( '   Num. of Time Windows='i5,',  Interval (sec)  =',F12.4)
 2340 FORMAT( '   Num. of Sub-Fault 1 for Integration =',i5)
 2345 FORMAT( '   Time Window=',i3,', Half Rise-Time 1 & 2 (s)=',2F10.5)
c      
      DO 300 I=1,NDIP
      DO 300 J=1,NSTR
      WRITE(10,2350)I,J
      WRITE(10,2355)(FD(i,j,k),k=1,ntw)
      WRITE(10,2356)(rak(i,j,k),k=1,ntw)
 2350 FORMAT(/'   Sub_Fault 0 Number   = (',I2,',',I2,')')
 2355 FORMAT( '   Dislocation for Each Time Windows (m)=',6F9.3) 
 2356 FORMAT( '   Rake Ang. for Each Time Windows (rad)=',6F9.5) 
      do id=1,nsub1
        write(10,2360)id
c          write(10,2362)(XF(i,j,id,ig),ig=1,nsub1)
c          write(10,2365)(YF(i,j,id,ig),ig=1,nsub1)
c          write(10,2370)(ZF(i,j,id,ig),ig=1,nsub1)
c          write(10,2375)(TD(i,j,id,ig),ig=1,nsub1)
      end do
  300 CONTINUE
 2360 format( '   Sub-Fault Numbers along Dip         =',2i5)
 2362 FORMAT( '     Sub_Fault 2: X-Axis along Str (m) =',5f9.1)
 2365 FORMAT( '     Sub_Fault 2: Y-Axis along Str (m) =',5f9.1)
 2370 FORMAT( '     Sub_Fault 2: Z-Axis along Str (m) =',5f9.1)
 2375 FORMAT( '     Total Time Delay      (sec)       =',5f9.3)
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 400 I=1,NOB
        READ(7,1520)XO(I),YO(I),ZO(I)
        WRITE(10,2520)I,XO(I),YO(I),ZO(I)
  400 CONTINUE
 1520 FORMAT(3F10.0)
 2520 FORMAT('   OB NUM.=',I5,',     (X,Y,Z) =',3F12.2)
C
      CLOSE(UNIT=7)
C
c *** output of fault and observation data ***
c
c      open(9,file='fault5_dat',status='unknown')
c        write(9,9500)XHYPO,YHYPO,ZHYPO
c        write(9,9500)xflt,yflt,zflt
c        write(9,9500)STR,DIP,0
c        write(9,9500)slength,dwidth,VRavg
c        write(9,9510)NSTR,NDIP,0
c        write(9,9510)nsub1,nsub1,0
c        do 700 i=1,ndip
c        do 700 j=1,nstr
c        do 700 id=1,nsub1
c        do 700 ig=1,nsub1
c        write(9,9500)xf(i,j,id,ig),yf(i,j,id,ig),ZF(i,j,id,ig)
c  700   continue
c 9500   format(3f12.3)
c 9510   format(3i12)
c      close(9) 
c
c      open(9,file='obs_dat',status='unknown')
c        do i=1,nob
c        write(9,9500)xo(i),yo(i),zo(i)
c        end do
c      close(9) 
c
c      open(9,file='time5_dat',status='unknown')
c        write(9,9600)SLENGTH,DWIDTH,0
c        write(9,9610)NSTR,NDIP,0
c        write(9,9610)nsub1,nsub1,0
c        do 800 i=1,ndip
c        do 800 id=1,nsub1
c 	do 800 j=1,nstr
c        do 800 ig=1,nsub1
c         da=fw*(i-1)
c          dd1=fw/nsub1
c          sa=fl*(j-1)
c          ds1=fl/nsub1
c            da1=da+dd1*(id-1)+dd1/2.0
c            sa1=sa+ds1*(ig-1)+ds1/2.0
c            write(9,9620)sa1,da1,TD(i,j,id,ig)
c  800   continue
c      close(9) 
c 9600 format(2f10.1,i5)
c 9610 format(3i5)
c 9620 format(2f10.1,f10.4)

      RETURN
      END
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 /OBS/XO(51),YO(51),ZO(51),NOB
     *    /STRC/DNS,VP,VS
     *    /FLT/STR,DIP,FL,FW,RT(2,6),TWI
     *    /IFLT/NSTR,NDIP,nsub1,NTW,iramp
     *    /SIGN/ISIGN /Hypo/xhypo,yhypo,zhypo 

C
C *** CHECK ***
C
      NC=51
      IF(NOB.GT.NC)THEN
         WRITE(10,*)'FROM OTDAT:NOB>NC',NOB,NC
         STOP
      END IF
C
C *** OUTPUT OF DATA (GX, GY AND GZ)***
C
      IF(ISIGN.EQ.1)THEN
        DO 100 IOB=1,NOB
        DO 110 K=IOM,NOM+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
        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)
 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
        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)
 3510 FORMAT(1X,I10,///)
 3520 FORMAT(1X,6E12.5)
C
      RETURN
      END
C
C
C *** COMPUTATION OF GREEN FUNCTIONS of the Homogeneous Full-Space ***
C
      SUBROUTINE HKQVW(MOM,MOB,GGX,GGY,GGZ,XF,YF,ZF,TD,RAK,FD
     *          ,MFLT,MSUB1)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VP,VS,SM
     *      ,GGX(MOM,MOB),GGY(MOM,MOB),GGZ(MOM,MOB)
     *      ,SMO,EDT0,EDT,GX,GY,GZ,GX0,GY0,GZ0
     *      ,UX(3,3),UY(3,3),UZ(3,3)
     *      ,Ukij,kai,kair,fair,esr,epr
     *      ,ss,sr,sr2,sr3,pp,pr,pr2,pr3,ps2
      DIMENSION EN(3,3),rd(3),RAK(MFLT,MFLT,6),FD(MFLT,MFLT,6)
     *         ,XF(MFLT,MFLT,MSUB1,MSUB1)
     *	       ,YF(MFLT,MFLT,MSUB1,MSUB1)
     *	       ,ZF(MFLT,MFLT,MSUB1,MSUB1)
     *	       ,TD(MFLT,MFLT,MSUB1,MSUB1)
      COMMON /FREQ/DOM,IOM,NOM /OBS/XO(51),YO(51),ZO(51),NOB
     *    /STRC/DNS,VP,VS
     *    /FLT/STR,DIP,FL,FW,RT(2,6),TWI
     *    /IFLT/NSTR,NDIP,nsub1,NTW,iramp
C
C *** CHECK OF DIMENSION ***
C
      NC=51
      IF(NOB.GE.NC)THEN
        WRITE(10,*)' FROM HKQVW:NOB>NC !!!!!!!!!!'
        WRITE(10,*)' NOB,NC=',NOB,NC
        STOP
      END IF
C
      NC=MOB
      IF(NOB.GE.NC)THEN
        WRITE(10,*)' FROM HKQVW:NOB>MOB !!!!!!!!!!'
        WRITE(10,*)' NOB,NC=',NOB,MOB
        STOP
      END IF
C
C *** GREEN'S FUNCTION ***
C
c   * Green's function computations *
c
      SM = DNS*VS*VS
      NN=nsub1*nsub1
c
      DO 100 IO=IOM,NOM+1
C
        OM=DOM*DFLOAT(IO-1)
        if(io.eq.1)om=0.001d0
        write(*,1000)om
 1000   format('    Omega =',f10.4,I6)

        ss=+DCMPLX(0.D0,1.d0)*om/vs
        pp=+DCMPLX(0.D0,1.d0)*om/vp
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
c       * Sub_Fault 0 and Green's function computation  *
c
          do 300 is=1,ndip
          do 300 js=1,nstr
c
c         * Sub_Fault 1 (Randam Rupture Velocities) *
c
            do 400 id=1,nsub1
	    do 400 ig=1,nsub1

                RX=XF(is,js,id,ig)-XO(IOB)
                RY=YF(is,js,id,ig)-YO(IOB)
                RZ=ZF(is,js,id,ig)-ZO(IOB)
                RR=DSQRT(RX*RX+RY*RY+RZ*RZ)
c
                EDT=CDEXP(DCMPLX(0.D0,OM*TD(is,js,id,ig)))
c
c	      * fai, kai, fai,r, and kai,r *
c
		sr=ss*rr
		sr2=sr*sr
		sr3=sr*sr2
		pr=pp*rr
		pr2=pr*pr
		pr3=pr*pr2
		ps2=pp*pp/(ss*ss)
c
		esr=cdexp(sr)/rr
		epr=cdexp(pr)/rr

                kai=(3.d0/sr2-3.d0/sr+1.d0)*esr
                kai=kai-ps2*(3.d0/pr2-3.d0/pr+1.d0)*epr
                kair=(-9.d0/sr3+9.d0/sr2-4.d0/sr+1.d0)*esr*ss
                kair=kair-ps2*(-9.d0/pr3
     *               +9.d0/pr2-4.d0/pr+1.d0)*epr*pp
                fair=(-3.d0/sr3+3.d0/sr2-2.d0/sr+1.d0)*esr*ss
                fair=fair-ps2*(-3.d0/pr3+3.d0/pr2-1.d0/pr)*epr*pp
c
c             * Uki,j *
c
                rd(1)=Rx/RR
                rd(2)=Ry/RR
                rd(3)=Rz/RR
                pai=4.d0*datan(1.d0)
c
c             * Uki,j=Uxi,j=UX(i,j) *
                k=1
                do i=1,3
                  Dki=0.d0
	          if(k.eq.i)Dki=1.d0
                  do j=1,3
                    Dkj=0.d0
	            if(k.eq.j)Dkj=1.d0
                    Dij=0.d0
	            if(i.eq.j)Dij=1.d0
                    call sUkij(Ukij,RR,rd(k),rd(i),rd(j)
     *                        ,Dki,Dkj,Dij,kai,fair,kair)
                    UX(i,j)=Ukij/(4.d0*pai*SM)
	          end do
	        end do
c
c             * Uki,j=Uyi,j=UY(i,j) *
                k=2
                do i=1,3
                  Dki=0.d0
	          if(k.eq.i)Dki=1.d0
                  do j=1,3
                    Dkj=0.d0
	            if(k.eq.j)Dkj=1.d0
                    Dij=0.d0
	            if(i.eq.j)Dij=1.d0
                    call sUkij(Ukij,RR,rd(k),rd(i),rd(j)
     *                        ,Dki,Dkj,Dij,kai,fair,kair)
                    UY(i,j)=Ukij/(4.d0*pai*SM)
	          end do
	        end do
c
c             * Uki,j=Uzi,j=UZ(i,j) *
                k=3
                do i=1,3
                  Dki=0.d0
	          if(k.eq.i)Dki=1.d0
                  do j=1,3
                    Dkj=0.d0
	            if(k.eq.j)Dkj=1.d0
                    Dij=0.d0
	            if(i.eq.j)Dij=1.d0
                    call sUkij(Ukij,RR,rd(k),rd(i),rd(j)
     *                        ,Dki,Dkj,Dij,kai,fair,kair)
                    UZ(i,j)=Ukij/(4.d0*pai*SM)
	          end do
	        end do
c
c               * Time Windows *
c
                do itw=1,ntw
                  if(io.eq.1)om=0.d0
                  CALL MO(is,js,itw,OM,en,SM,SMO,RAK,FD,MFLT)
                  if(io.eq.1)om=0.001d0
                  EDT0=CDEXP(DCMPLX(0.D0,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)
                  end do
		  end do
                  GX0=GX*SMO*EDT*EDT0/NN 
                  GY0=GY*SMO*EDT*EDT0/NN 
                  GZ0=GZ*SMO*EDT*EDT0/NN 
                  GGX(IO,IOB)=GGX(IO,IOB)+GX0
                  GGY(IO,IOB)=GGY(IO,IOB)+GY0
                  GGZ(IO,IOB)=GGZ(IO,IOB)+GZ0
		end do
c    
  400       CONTINUE
c
  300     CONTINUE

      write(*,*)'GX,(io,iob)',io,iob
      write(*,*)GGX(io,iob)
C
  200   CONTINUE
C
  100 CONTINUE
C
      RETURN
      END
C
c
c *** Green Funtion Uki,j ***
c
      subroutine sUkij(Ukij,RR,Rk,Ri,Rj,Dki,Dkj,Dij,kai,fair,kair)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 Ukij,kai,kair,fair

      Ukij=fair*Rj*Dki-kair*Ri*Rj*Rk
      Ukij=Ukij-(Ri*Dkj+Rk*Dij-2*Ri*Rj*Rk)*kai/RR
C
      RETURN
      END
c
c
c
C *** MOMENT TENSOR FOR THE SOURCE ***
C
      SUBROUTINE MO(i,j,itw,OM,en,SM,SMO,RAK,FD,MFLT)
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION EN(3,3),RAK(MFLT,MFLT,6),FD(MFLT,MFLT,6)
      COMPLEX*16 VS,VP,SMO,CPX,SM,STT
      COMMON/STRC/DNS,VP,VS
     *    /FLT/STR,DIP,FL,FW,RT(2,6),TWI
     *    /IFLT/NSTR,NDIP,nsub1,NTW,iramp
c
C *** EN ***
C
      CPX=(0.D0,1.D0)
c
      CST= DCOS(STR)
      SST= DSIN(STR)
C
      CDP= DCOS(DIP)
      SDP= DSIN(DIP)
C
      CRK= DCOS(RAK(i,j,itw))
      SRK= DSIN(RAK(i,j,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(i,j,itw)*SM*FL*FW
C
      TTK= OM*(RT(1,itw)+RT(2,itw))/2.D0
      if(iramp.eq.0)then
        IF(om.EQ.0.D0)THEN
          STT= (1.D0,0.D0)
        ELSE
          STT= DSIN(TTK)*CDEXP(CPX*TTK)/TTK
        END IF
      else
        IF(OM.EQ.0.D0)THEN
          STT= (1.D0,0.D0)
        ELSE
          TK1= OM*RT(1,itw)/2.D0
          TK2= OM*RT(2,itw)/2.D0
          STT= DSIN(TK1)/TK1-DSIN(TK2)/TK2*CDEXP(CPX*TTK)
          STT= CPX*STT*CDEXP(CPX*TK1)/TTK
        END IF
      end if
C
      SMO=SMO*STT
C
      RETURN
      END

