C     
c *** Green's functions using Boor's method
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    (Originally made in 1992, and modified in 1998)
c
c * Parameters : 
c     MOM : Number of Omega
c     MOB : NUmber of Observation Points
c
      PROGRAM grfltHF6
C
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(MOM=8200,MOB=13,MFLT=7)
      COMPLEX*16 GX(MOM,MOB),GY(MOM,MOB),GZ(MOM,MOB),FX(MOM)
      DIMENSION FD(MFLT,MFLT),TD(MFLT,MFLT),XR(MOM)
     *         ,XF(MFLT,MFLT),YF(MFLT,MFLT),ZF(MFLT,MFLT)
C
      WRITE(*,*)
      WRITE(*,*)'"Green functions using the Boore method'
      WRITE(*,*)
     *'           (Originally made in 1992, and modified in 1998)'
      write(*,*)
      WRITE(*,*)'"grfltHF61.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,FD,MFLT)
      CALL BRSAH(MFLT,MOM,MOB,XF,YF,ZF,FD,TD,GX,FX,XR)
      CALL BRSAH(MFLT,MOM,MOB,XF,YF,ZF,FD,TD,GY,FX,XR)
      CALL BRSAV(MFLT,MOM,MOB,XF,YF,ZF,FD,TD,GZ,FX,XR)
      CALL OTDAT(MOM,MOB,GX,GY,GZ)
C
      STOP
      END
C
c
C *** DATA INPUT ***
C
      SUBROUTINE INDAT(MOM,MOB,XF,YF,ZF,TD,FD,MFLT)
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION FD(MFLT,MFLT),TD(MFLT,MFLT)
     *         ,XF(MFLT,MFLT),YF(MFLT,MFLT),ZF(MFLT,MFLT)
      COMMON /FREQ/DOM,IOM,NOM /OBS/XO(13),YO(13),ZO(13),NOB
     *    /STRC/DNS,VP,VS,Qp,Qs,Qpn,Qsn
     *    /FLT/STR,DIP,FL,FW,Vr,Fmax,Sdrop,Rtime
     *    /IFLT/NSTR,NDIP,idum
     *    /Hypo/xhypo,yhypo,zhypo 
     *    /SIGN/ISIGN
c
C *** PARAMETERS IN COMMOM STATEMENT ***
C
      NMOB=13
      IF(MOB.LT.NMOB)NMOB=MOB
C
C *** OMEGA DATA ***
C
      write(*,*)
      OPEN (UNIT=7,FILE='grfltHF61.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 grfltHF6.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,VP,QP,VS,QS
c      WRITE(*,2210)DNS,VP,Qp,VS,Qs
      READ(7,1210)Qpn,Qsn
c      WRITE(*,2220)Qpn,Qsn
 1200 FORMAT(/6f8.0)
 1210 FORMAT(2f8.0)
C
      WRITE(10,2200)
      WRITE(10,2210)DNS,VP,Qp,VS,Qs
      WRITE(10,2220)Qpn,Qsn
 2200 FORMAT(/' Density , Vp(m/s), Qp0     '
     *        ', Vs(m/s) , Qs0    ,Thickness(m)')
 2210 FORMAT(5F10.3,F12.3)
 2220 format(' n of Q*f^n for Qp and Qs =',2f8.2)
C
C *** FAULT DATA ***
C
      READ(7,1300)SLENGTH,DWIDTH
c      write(*,*)'SLENGTH,DWIDTH',SLENGTH,DWIDTH
      READ(7,1310)NSTR,NDIP
c      write(*,*)'NSTR,NDIP',NSTR,NDIP
      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
      READ(7,1320)DELAY,Sdrop
c      write(*,*)'DELAY,Sdrop',DELAY,Sdrop
      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,Rtime
c      write(*,*)'STR,DIP',STR,DIP
      str=str*0.01745329252d0
      dip=dip*0.01745329252d0
      READ(7,1340)VR,Fmax,idum
c      write(*,*)'Vr,Fmax,idum',Vr,Fmax,idum
 1300 FORMAT(/2f10.0)
 1310 FORMAT(2I5)
 1320 FORMAT(2F10.0)
 1330 format(3f10.0)
 1340 format(2f10.0,i10)
c
      read(7,*)
      do i=1,ndip
      do j=1,nstr
        READ(7,*)FD(i,j)
c        write(*,*)'i,j,FD(i,j,k)',i,j,FD(i,j)
      end do
      end do
c
c   * Locations of Sub_Fault Centers *  
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							
c   * Dip direction *     
c
      do 200 i=1,ndip

        da=fw*(i-1)+fw/2
	da1=da-dhypo

c     * Strike direction *
        do 210 j=1,nstr   
   
          sa=fl*(j-1)+fl/2
	  sa1=sa-shypo
          dis=dsqrt(da1*da1+sa1*sa1)
          ta=dis/Vr
          xa=xflt+sa*cst+da*cdp*sst
          ya=yflt+sa*sst-da*cdp*cst
          za=zflt-da*sdp
          if(za.lt.0.d0)then
             write(*,*)'za < 0 !!'
             stop
          end if
          xf(i,j)=xa
          yf(i,j)=ya
          ZF(i,j)=za
          td(i,j)=delay+ta

  210   continue
  200 continue
C
      WRITE(10,2300)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)VR,Fmax,Sdrop,Rtime
 2300 FORMAT(//' Seismic Source Data:' 
     *      //'   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( '   Vr         (m/s)=',F12.3,',  Fmax        (Hz)=',F12.3)
 2335 FORMAT( '   Num. of Time Windows='i5,',  Interval (sec)  =',F12.4
     *       /'   Stress Drop(bar)=',F12.3,',  Rise Time  (sec)=',F12.3)
c      
      DO 300 I=1,NDIP
      DO 300 J=1,NSTR
      WRITE(10,2350)I,J
      write(10,2355)FD(i,j),TD(i,j)
      write(10,2360)XF(i,j),YF(i,j),ZF(i,j)
  300 continue
 2350 FORMAT(/'   Sub_Fault 0 Number   = (',I2,',',I2,')')
 2355 format( '   Dislocation (m) and Delay Time (s) =',2F10.2) 
 2360 format( '   Locations of Center: X, Y, Z   (m)=',3f10.2)
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
      open(9,file='fault5_dat',status='unknown')
        write(9,9500)XHYPO,YHYPO,ZHYPO
        write(9,9500)xflt,yflt,zflt
        write(9,9500)STR,DIP,0
        write(9,9500)slength,dwidth,VR
        write(9,9510)NSTR,NDIP,0
        write(9,9510)1,1,0
        do 700 i=1,ndip
        do 700 j=1,nstr
        write(9,9500)xf(i,j),yf(i,j),ZF(i,j)
  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) 
c
      open(9,file='time5_dat',status='unknown')
        write(9,9600)SLENGTH,DWIDTH,0
        write(9,9610)NSTR,NDIP,0
        write(9,9610)1,1,0
        do 800 i=1,ndip
	do 800 j=1,nstr
          da=fw*(i-1)+fw/2
          sa=fl*(j-1)+fl/2
          write(9,9620)sa,da,TD(i,j)
  800   continue
      close(9) 
 9600 format(2f10.1,i5)
 9610 format(3i5)
 9620 format(2f10.1,f10.4)

      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)
      COMMON /FREQ/DOM,IOM,NOM /OBS/XO(13),YO(13),ZO(13),NOB
     *    /STRC/DNS,VP,VS,Qp,Qs,Qpn,Qsn
     *    /FLT/STR,DIP,FL,FW,Vr,Fmax,Sdrop,Rtime
     *    /IFLT/NSTR,NDIP,idum
     *    /Hypo/xhypo,yhypo,zhypo 
     *    /SIGN/ISIGN
C
C *** CHECK ***
C
      NC=13
      IF(NOB.GT.NC)THEN
         WRITE(10,*)'FROM OTDAT:NOB>NC',NOB,NC
         STOP
      END IF
C
C *** OUTPUT OF DATA (GX, GY AND GZ)***
C
      IF(ISIGN.EQ.1)THEN
        DO 100 IOB=1,NOB
        DO 110 K=IOM,NOM
          GX(K,IOB)=DCMPLX(DREAL(GX(K,IOB)),-DIMAG(GX(K,IOB)))
          GY(K,IOB)=DCMPLX(DREAL(GY(K,IOB)),-DIMAG(GY(K,IOB)))
          GZ(K,IOB)=DCMPLX(DREAL(GZ(K,IOB)),-DIMAG(GZ(K,IOB)))
  110   CONTINUE
  100   CONTINUE
      END IF
C      
      vfac=0.6d0
      write(*,*)'Amplitude of V over is assumed as',vfac
      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)*vfac
  510     CONTINUE
  500   CONTINUE
      CLOSE(8)
C
 2500 FORMAT(1X,F10.6,2I10)
 2510 FORMAT(1X,I10,///)
 2520 FORMAT(1X,6E12.5)
C
      RETURN
      END
C
C
C *** COMPUTATION OF GREEN FUNCTIONS of the Homogeneous Full-Space ***
C
      SUBROUTINE BRSAH(MFLT,MOM,MOB,XF,YF,ZF,FD,TD,GX,FX,XR)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 GX(MOM,MOB),FX(MOM),EDT
      DIMENSION FD(MFLT,MFLT),TD(MFLT,MFLT),XR(MOM)
     *         ,XF(MFLT,MFLT),YF(MFLT,MFLT),ZF(MFLT,MFLT)
      COMMON /FREQ/DOM,IOM,NOM /OBS/XO(13),YO(13),ZO(13),NOB
     *    /STRC/DNS,VP,VS,Qp,Qs,Qpn,Qsn
     *    /FLT/STR,DIP,FL,FW,Vr,Fmax,Sdrop,Rtime
     *    /IFLT/NSTR,NDIP,idum
     *    /SIGN/ISIGN
C
C *** CHECK OF DIMENSION ***
C
      NC=13
      IF(NOB.GE.NC)THEN
        WRITE(10,*)' FROM BRQVW:NOB>NC !!!!!!!!!!'
        WRITE(10,*)' NOB,NC=',NOB,NC
        STOP
      END IF
C
      NC=MOB
      IF(NOB.GE.NC)THEN
        WRITE(10,*)' FROM BRQVW:NOB>MOB !!!!!!!!!!'
        WRITE(10,*)' NOB,NC=',NOB,MOB
        STOP
      END IF
C
C   * Find closest (but larger) 2**J value to the input data dimension *
C
      nom1=nom+1
      DO J = 1,999
        N = 2**J
        IF(NOM1.LE.N)THEN
             NUM = N
             GO TO 10
        END IF
      END DO
C
   10 IF(NUM*2.GT.MOM)THEN
        WRITE(*,*)
        WRITE(*,*)' NUM*2.GT.MOM !!!!!'
        WRITE(*,*)' Increase MOM in the parameter statement!'
        STOP
      ENDIF
c
      pi=datan(1.d0)*4.d0
      DUR = pi*2.d0/DOM
      NFAST = NUM*2
      DT = DUR/dFLOAT(NFAST)
      write(*,*)
      write(*,*)'dur,nfast,dt',dur,nfast,dt
      write(*,*)'New NOM for FFT, and Final F',num,dom*num/(pi*2)
      write(*,*)
c
      Tdu=sqrt(FL*FL+FW*FW)/Vr
      Tw=2.d0*Tdu
      Afac=dur/Tw
      write(*,*)'Tw,Afac',Tw,Afac
C
c   * Source Spectra *
c
c
      aND=dsqrt(dfloat(ndip*nstr))
      Rpatn=0.63d0
      PRTITN=1.d0/dsqrt(2.d0)
      DO 200 IOB=1,NOB

	do i=1,mom
          GX(i,iob)=(0.d0,0.d0)
	end do
c
        do 200 is=1,ndip
        do 200 js=1,nstr

          RX=XF(is,js)-XO(IOB)
          RY=YF(is,js)-YO(IOB)
          RZ=ZF(is,js)-ZO(IOB)
          RR=DSQRT(RX*RX+RY*RY+RZ*RZ)
          slip=FD(is,js)/aND
c
          SMo=Dns*Vs*Vs*slip*FL*FW*10000000000.0d0
          if(SMo.eq.0.d0)then
            Fc=9999.d0
          else
            Fc=4900*Vs*(Sdrop/SMo)**(1.d0/3.d0)
	  end if
c
          FX(1)=(0.d0,0.d0)
	  CM=Rpatn*PRTITN*slip*FL*FW/(4.d0*pi*Vs)
          DO 100 IO=IOM,Num      
            OM=DOM*DFLOAT(IO)
            f=om/(2.d0*pi)
	    FFc=f/Fc
	    S=om*om/(1.d0+FFc*FFc)
            P=1.d0/(1.d0+f/Fmax)
              Q=Qs*f**Qsn
            E=dexp(-pi*f*RR/(Vs*Q))
            SA=CM*S*P*E/RR

            box=om*Rtime/2.d0
            EDT=CDEXP(DCMPLX(0.D0,box))

            ran=ran0(idum)*pi*2.d0
            dr=SA*dcos(ran)
            di=SA*dsin(ran)
          FX(io+1)=dcmplx(dr,di)*(1.d0+(aND-1.d0)*dsin(box)*EDT/box)
c          write(3,9999)om,idum,ran,fx(io+1)*100
  100     CONTINUE
 9999 format(f12.4,i12,4f12.5)
c
c   * Waveform using the inverse FFT ***
C
          IF(ISIGN.EQ.1)THEN
            DO K=1,NOM1
              FX(K)=DCMPLX(DREAL(FX(K)),-DIMAG(FX(K)))
	    end do
          END IF
C
          DO io = NOM1+1,NUM+1
            FX(io) = DCMPLX(0.0,0.0)
          end do
C
          DO I=2,NFAST/2
            FX(NFAST+2-I)=dCMPLX(dREAL(FX(I)),-DIMAG(FX(I)))
          end do
C
          CALL dWVFAST(NFAST,FX,MOM,1)
c
          DO I = 1,NFAST
            XR(I)=dREAL(FX(I))/DUR
          end do
c
c   * RE-construction of waveform considering Boore's waveform functions *
c
	  b=1.2531d0
          c=6.266/Tw
          a=(2.718281829d0/(0.2*Tw))**b
          do i=1,nfast
	    time=dt*(i-1)
c            xro=XR(i)
            XR(i)=XR(i)*Afac*a*(time**b)*dexp(-c*time)
c       write(2,9998)time,xro*100,xr(i)*100
	  end do
 9998 format(f12.4,4f12.5)
c
c  * Fourier Transform and *
c
	  do i=1,nfast
            FX(I)=dCMPLX(XR(I),0.0d0)
	  end do
c
          CALL dWVFAST(nfast,FX,MOM,-1)
C
          IF(ISIGN.EQ.1)THEN
            DO K=1,nfast
              FX(K)=FX(K)*dt
              FX(K)=DCMPLX(DREAL(FX(K)),-DIMAG(FX(K)))
c        write(2,9999)dom*(k-1),fx(k)
	    end do
          END IF
C
          DO I=1,num
	    om=dom*i
            EDT=CDEXP(DCMPLX(0.D0,OM*(TD(is,js)+RR/Vs)))
            GX(i,iob)=GX(i,iob)-FX(i+1)*EDT/(om*om)
          end do
C
  200 CONTINUE
c
      RETURN
      END
c
C
C *** COMPUTATION OF GREEN FUNCTIONS of the Homogeneous Full-Space ***
C
      SUBROUTINE BRSAV(MFLT,MOM,MOB,XF,YF,ZF,FD,TD,GZ,FX,XR)
c
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 GZ(MOM,MOB),FX(MOM),EDT
      DIMENSION FD(MFLT,MFLT),TD(MFLT,MFLT),XR(MOM)
     *         ,XF(MFLT,MFLT),YF(MFLT,MFLT),ZF(MFLT,MFLT)
      COMMON /FREQ/DOM,IOM,NOM /OBS/XO(13),YO(13),ZO(13),NOB
     *    /STRC/DNS,VP,VS,Qp,Qs,Qpn,Qsn
     *    /FLT/STR,DIP,FL,FW,Vr,Fmax,Sdrop,Rtime
     *    /IFLT/NSTR,NDIP,idum
     *    /SIGN/ISIGN
C
C *** CHECK OF DIMENSION ***
C
      NC=13
      IF(NOB.GE.NC)THEN
        WRITE(10,*)' FROM BRQVW:NOB>NC !!!!!!!!!!'
        WRITE(10,*)' NOB,NC=',NOB,NC
        STOP
      END IF
C
      NC=MOB
      IF(NOB.GE.NC)THEN
        WRITE(10,*)' FROM BRQVW:NOB>MOB !!!!!!!!!!'
        WRITE(10,*)' NOB,NC=',NOB,MOB
        STOP
      END IF
C
C   * Find closest (but larger) 2**J value to the input data dimension *
C
      nom1=nom+1
      DO J = 1,999
        N = 2**J
        IF(NOM1.LE.N)THEN
             NUM = N
             GO TO 10
        END IF
      END DO
C
   10 IF(NUM*2.GT.MOM)THEN
        WRITE(*,*)
        WRITE(*,*)' NUM*2.GT.MOM !!!!!'
        WRITE(*,*)' Increase MOM in the parameter statement!'
        STOP
      ENDIF
c
      pi=datan(1.d0)*4.d0
      DUR = pi*2.d0/DOM
      NFAST = NUM*2
      DT = DUR/dFLOAT(NFAST)
      write(*,*)
      write(*,*)'dur,nfast,dt',dur,nfast,dt
      write(*,*)'New NOM for FFT, and Final F',num,dom*num/(pi*2)
      write(*,*)
c
      Tdu=sqrt(FL*FL+FW*FW)/Vr
      Tw=2.d0*Tdu
      Afac=dur/Tw
      write(*,*)'Tw,Afac',Tw,Afac
C
c   * Source Spectra *
c
c
      aND=dsqrt(dfloat(ndip*nstr))
      Rpatn=0.63d0
      PRTITN=1.d0/dsqrt(2.d0)
      DO 200 IOB=1,NOB

	do i=1,mom
          GZ(i,iob)=(0.d0,0.d0)
	end do
c
        do 200 is=1,ndip
        do 200 js=1,nstr

          RX=XF(is,js)-XO(IOB)
          RY=YF(is,js)-YO(IOB)
          RZ=ZF(is,js)-ZO(IOB)
          RR=DSQRT(RX*RX+RY*RY+RZ*RZ)
          slip=FD(is,js)/aND
c
          SMo=Dns*Vs*Vs*slip*FL*FW*10000000000.0d0
          if(SMo.eq.0.d0)then
            Fc=9999.d0
          else
            Fc=4900*Vs*(Sdrop/SMo)**(1.d0/3.d0)
	  end if
c
          FX(1)=(0.d0,0.d0)
	  CM=Rpatn*PRTITN*slip*FL*FW/(4.d0*pi*Vs)
          DO 100 IO=IOM,Num      
            OM=DOM*DFLOAT(IO)
            f=om/(2.d0*pi)
	    FFc=f/Fc
	    S=om*om/(1.d0+FFc*FFc)
            P=1.d0/(1.d0+f/Fmax)
              Q=Qs*f**Qsn
            E=dexp(-pi*f*RR/(Vs*Q))
            SA=CM*S*P*E/RR

            box=om*Rtime/2.d0
            EDT=CDEXP(DCMPLX(0.D0,box))

            ran=ran0(idum)*pi*2.d0
            dr=SA*dcos(ran)
            di=SA*dsin(ran)
          FX(io+1)=dcmplx(dr,di)*(1.d0+(aND-1.d0)*dsin(box)*EDT/box)
c          write(3,9999)om,idum,ran,fx(io+1)*100
  100     CONTINUE
 9999 format(f12.4,i12,4f12.5)
c
c   * Waveform using the inverse FFT ***
C
          IF(ISIGN.EQ.1)THEN
            DO K=1,NOM1
              FX(K)=DCMPLX(DREAL(FX(K)),-DIMAG(FX(K)))
	    end do
          END IF
C
          DO io = NOM1+1,NUM+1
            FX(io) = DCMPLX(0.0,0.0)
          end do
C
          DO I=2,NFAST/2
            FX(NFAST+2-I)=dCMPLX(dREAL(FX(I)),-DIMAG(FX(I)))
          end do
C
          CALL dWVFAST(NFAST,FX,MOM,1)
c
          DO I = 1,NFAST
            XR(I)=dREAL(FX(I))/DUR
          end do
c
c   * RE-construction of waveform considering Boore's waveform functions *
c
	  b=1.2531d0
          c=6.266/Tw
          a=(2.718281829d0/(0.2*Tw))**b
          do i=1,nfast
	    time=dt*(i-1)
            XR(i)=XR(i)*Afac*a*(time**b)*dexp(-c*time)
c       write(1,9999)time,xr(i),xr(i),xr(i)
	  end do
c
c  * Fourier Transform and *
c
	  do i=1,nfast
            FX(I)=dCMPLX(XR(I),0.0d0)
	  end do
c
          CALL dWVFAST(nfast,FX,MOM,-1)
C
          IF(ISIGN.EQ.1)THEN
            DO K=1,nfast
              FX(K)=FX(K)*dt
              FX(K)=DCMPLX(DREAL(FX(K)),-DIMAG(FX(K)))
c        write(2,9999)dom*(k-1),fx(k)
	    end do
          END IF
C
          DO I=1,num
	    om=dom*i
            EDT=CDEXP(DCMPLX(0.D0,OM*(TD(is,js)+RR/Vp)))
            GZ(i,iob)=GZ(i,iob)-FX(i+1)*EDT/(om*om)
          end do
C
  200 CONTINUE
c
      RETURN
      END
C
c
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 *** FOURIER TRANSFORM OR INVERSE FOURIER TRANSFORM
C
      SUBROUTINE dWVFAST(N,A,MOM,IND)
c
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 A(MOM),TEMP,THETA
C
C * FFT and inverse FFT *
C
      PAI=3.1415926540d0
      J=1
      DO 140 I=1,N
        IF(I.GE.J) GO TO 110
        TEMP=A(J)
        A(J)=A(I)
        A(I)=TEMP
  110   M=N/2
  120   IF(J.LE.M) GO TO 130
        J=J-M
        M=M/2
        IF(M.GE.2) GO TO 120
  130   J=J+M
  140 CONTINUE
      KMAX=1
  150 IF(KMAX.GE.N) RETURN
      ISTEP=KMAX*2
      DO 170 K=1,KMAX
        THETA=dCMPLX(0.d0,PAI*dFLOAT(IND*(K-1))/dFLOAT(KMAX))
        DO 160 I=K,N,ISTEP
          J=I+KMAX
          TEMP=A(J)*cdEXP(THETA)
          A(J)=A(I)-TEMP
          A(I)=A(I)+TEMP
  160   CONTINUE
  170 CONTINUE
      KMAX=ISTEP
      GO TO 150
C
      END
