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
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 grfltHF3
C
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(MOM=351,MOB=51)
      COMPLEX*16 GX(MOM,MOB),GY(MOM,MOB),GZ(MOM,MOB)
C
      WRITE(*,*)
      WRITE(*,*)'"grfltHF3.f" constructs Green functions'
      WRITE(*,*)
     *'           (Originally made in 1992, and modified in 1998)'
      write(*,*)
      WRITE(*,*)'"grfltHF3.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)
      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
      COMMON /FREQ/DOM,IOM,NOM /OBS/XO(51),YO(51),ZO(51),NOB
     *    /STRC/DNS,VP,VS
     *    /FLT/STR,DIP,RAK(203,6),FL,FW,FD(203,6),RT(2,6),TWI
     *     ,XF(203,21,21),YF(203,21,21),TD(203,21,21),DP(203,21)
     *    /IFLT/NF,NTW,IS,nsub
     *    /Hypo/xhypo,yhypo,zhypo 
     *    /SIGN/ISIGN /ramp/iramp,NSTR,NDIP
c
C *** PARAMETERS IN COMMOM STATEMENT ***
C
      NMOB=51
      NFLT=203
      IF(MOB.LT.NMOB)NMOB=MOB
C
C *** OMEGA DATA ***
C
      write(*,*)
      OPEN (UNIT=7,FILE='grfltHF3.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,nsub
c      write(*,*)'NSTR,NDIP,nsub',NSTR,NDIP,nsub
      NF=NSTR*NDIP
      IF(NF.GT.NFLT)THEN
        WRITE(10,*)' NF.GT.NFLT !!!!!!'
        STOP
      END IF
      if(nsub.gt.21)then
        write(10,*)' nsub > 21 !!!'
        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,1332)VRavg,TMstd,idum
      READ(7,1335)NTW,TWI
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 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,I10)
 1330 format(3f10.0)
 1332 format(2f10.0,i5)
 1335 format(i8,f8.0)
 1340 format(2f8.0)
c 1345 format(6f6.0)
c
c   * Generation of Rupture Velocities at Sub-faults using the Gaussian 
c     distributed deviate 
c							
c   * Dip direction *      
C
      dcs=fw*dcos(dip)
      cst=dcos(str)
      sst=dsin(str)
c
      iflt=0
      do 200 i=1,ndip

        da=dcs*(i-1)
        db=dcs*i
        dd=dcs/nsub
        df=da+dd/2.d0
        za=depth-da*dtan(dip)
        zb=depth-db*dtan(dip)
        dz=(za-zb)/nsub
        zf=za-dz/2.d0

        if((za.lt.0.d0).or.(zb.lt.0.d0))then
          write(*,*)'Fault Depth is less than 0!: Dip Number=',i
          stop
        end if

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

          sa=fl*(j-1)
          ds=fl/nsub
          sf=sa+ds/2.d0
          do 240 id=1,nsub

            zd=zf-dz*(id-1)
            dp(iflt,id)=zd
            zz=zd-zhypo
            ddd=df+dd*(id-1)

            do ig=1,nsub
              sg=sf+ds*(ig-1)
              xf(iflt,id,ig)=xflt+sg*cst+ddd*sst
              yf(iflt,id,ig)=yflt+sg*sst-ddd*cst
	      xx=xf(iflt,id,ig)-xhypo
	      yy=yf(iflt,id,ig)-yhypo
              dis=dsqrt(xx*xx+yy*yy+zz*zz)
              dev=gasdev(idum)
              td(iflt,id,ig)=delay+dis/VRavg+TMstd*dev
              if(td(iflt,id,ig).lt.0.d0)td(iflt,id,ig)=0.d0
            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,depth
      WRITE(10,2330)VRavg, TMstd
      WRITE(10,2335)ntw,twi
      WRITE(10,2340)nsub,nsub*nsub
      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
     *        ',  Z Fault Edge (m)=',F12.3)
 2330 FORMAT( '   Average Vr (m/s)=',F12.3,',  Std. Dev.Time(s)=',F12.3)
 2335 FORMAT( '   Num. of Time Windows='i5,',  Interval (sec)  =',F12.4)
 2340 FORMAT( '   Num. of Integration Points =',i5
     *       ,',  Num. of Integration 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,nsub
        write(10,2360)id,DP(II,id)
        write(10,2365)(XF(II,id,ig),ig=1,nsub)
        write(10,2370)(YF(II,id,ig),ig=1,nsub)
        write(10,2375)(TD(II,id,ig),ig=1,nsub)
      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( '   Integration Point Number along Dip=',i5
     *        ',  Depth (m) =',f8.2)
 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 *** 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,VRavg
        write(9,9510)NSTR,NDIP,nsub
        do 700 iflt=1,nf
        do 700 id=1,nsub
        do 700 ig=1,nsub
        write(9,9500)xf(iflt,id,ig),yf(iflt,id,ig),dp(iflt,id)
  700   continue
 9500   format(3f12.3)
 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) 

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

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

      return
      end

c *** Randam Function **
c
      function ran0(idum)
      IMPLICIT REAL*8(A-H,O-Z)
      parameter(IA=16807,IM=2147483647,AM=1.d0/IM,IQ=127773
     *         ,IR=2836,MASK=123459876)
c
      idum=ieor(idum,MASK)
      k=idum/IQ
      idum=IA*(idum-k*IQ)-IR*k
      if(idum.lt.0)idum=idum+IM
      ran0=AM*idum
      idum=ieor(idum,MASK)
c
      return
      end
C
C
C *** 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,RAK(203,6),FL,FW,FD(203,6),RT(2,6),TWI
     *     ,XF(203,21,21),YF(203,21,21),TD(203,21,21),DP(203,21) 
     *    /IFLT/NF,NTW,IS,nsub 
     *    /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
      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
        WRITE(8,2510)NOB
        DO 500 IOB=1,NOB
          DO 510 K=IOM,NOM
            WRITE(8,2520)GX(K,IOB),GY(K,IOB),GZ(K,IOB)
  510     CONTINUE
  500   CONTINUE
      CLOSE(8)
C
 2500 FORMAT(1X,F10.6,2I10)
 2510 FORMAT(1X,I10,///)
 2520 FORMAT(1X,6E12.5)
C
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
        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)
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)
      COMMON /FREQ/DOM,IOM,NOM /OBS/XO(51),YO(51),ZO(51),NOB
     *    /STRC/DNS,VP,VS
     *    /FLT/STR,DIP,RAK(203,6),FL,FW,FD(203,6),RT(2,6),TWI
     *     ,XF(203,21,21),YF(203,21,21),TD(203,21,21),DP(203,21)  
     *    /IFLT/NF,NTW,IS,nsub 
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
      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   * Green's function computations *
c
      DO 100 IO=IOM,NOM
C
        OM=DOM*DFLOAT(IO)
        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       * fault points and Green's function computation  *
c
          DO 300 IS=1,NF
c
c           * Gaussian Quadrature (id*ig points) *
c
            do 400 id=1,nsub
c
              SM = DNS*VS*VS
              RZ=DP(IS,id)-ZO(IOB)
C
	      do 450 ig=1,nsub

                RX=XF(IS,id,ig)-XO(IOB)
                RY=YF(IS,id,ig)-YO(IOB)
                RR=DSQRT(RX*RX+RY*RY+RZ*RZ)
c
                EDT=CDEXP(DCMPLX(0.D0,OM*TD(IS,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
c
                  CALL MO(itw,OM,en,SM,SMO)
                  EDT0=CDEXP(DCMPLX(0.D0,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/(nsub*nsub) 
                  GY0=GY*SMO*EDT*EDT0/(nsub*nsub) 
                  GZ0=GZ*SMO*EDT*EDT0/(nsub*nsub) 
                  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
      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(itw,OM,en,SM,SMO)
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION EN(3,3)
      COMPLEX*16 VS,VP,SMO,CPX,SM,STT
      COMMON/STRC/DNS,VP,VS
     *    /FLT/STR,DIP,RAK(203,6),FL,FW,FD(203,6),RT(2,6),TWI
     *     ,XF(203,21,21),YF(203,21,21),TD(203,21,21),DP(203,21)
     *    /IFLT/NF,NTW,IS,nsub /ramp/iramp,NSTR,NDIP
C
C *** CHECK OF DIMENSION ***
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
      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=CPX*SMO*STT/OM
C
      RETURN
      END
