c
c *** This code is applied to the S11 model of the 2009 Benchmark test
c     (This code is the same as omega2009-v2a.f)
c
c   * Parameters :
c
c     MOM : Number of Omega
c     MOB : NUmber of Observation Points
c     MFL : Number of Sub-Faults
c
      PROGRAM omega2x
C
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(MOM=8211,MOB=51,MFL=203)
      COMPLEX*16 GX(MOM,MOB),GY(MOM,MOB),GZ(MOM,MOB),Frate(MOM,MFL,2)
      DIMENSION rate(mom,MFL,2)
C
      WRITE(*,*)
      WRITE(*,*)'"omega2009.f" simulates broad-band ground motions'
      WRITE(*,*)'   using the Boore source model and the Green function'
      write(*,*)'   of the direct S-wave'
      write(*,*)' (Originally made in 1992, and modified in 2009)'
      WRITE(*,*)
     *' For questions, e-mail to "hisada@cc.kogakuin.ac.jp".'
      WRITE(*,*)' This code is open for the academic use only.'
C
      CALL INDAT(MOM,MOB,MFL,Frate,rate)
      CALL HKQVW(MOM,MOB,MFL,Frate,GX,GY,GZ)
      CALL OTDAT(MOM,MOB,GX,GY,GZ)
      WRITE(*,*)
      WRITE(*,*)'Use FFT program "grfftpxA.f" to plot waveforms'
C
      STOP
      END
C
c
C *** DATA INPUT ***
C
      SUBROUTINE INDAT(MOM,MOB,MFL,Frate,rate)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 Frate(MOM,MFL,2),cpx,Ffunc
      DIMENSION rate(mom,MFL,2),sdrop(203)
      CHARACTER INFILE*40
      COMMON /FREQ/DOM,IOM,NOM /OBS/XO(51),YO(51),ZO(51),NOB
     *    /time/dt,period  /ntime/nt  /STRC/DNS,VP,QP,FP,VS,QS,FS
     *    /FLT/STR,DIP,RAK(203),FL,FW,FD(203),VR,TLc,Pn(203)
     *        ,XF(203),YF(203),TD(203),DP(203),fc(203),fmax(203)
     *    /Hypo/xhypo,yhypo,zhypo /Radiation/fR1,fR2,fran
     *    /IFLT/NF,NSTR,NDIP,NCF,NFfunc /SIGN/ISIGN /Nboore/Niter,idum
c
c *** Check of Dimensions of Common Paramters ***
c
      nc=51
      if(nc.lt.MOB)then
        write(10,*)'from Indat:'
        WRITE(10,*)' NC < MOB !!!!!!!!!!!!!'
        write(*,*)'from Indat:'
        WRITE(*,*)' NC < MOB !!!!!!!!!!!!!'
        STOP
      END IF
c
      nc=203
      if(nc.lt.MFL)then
        write(10,*)'from Indat:'
        WRITE(10,*)' NC < MFL !!!!!!!!!!!!!'
        write(*,*)'from Indat:'
        WRITE(*,*)' NC < MFL !!!!!!!!!!!!!'
        STOP
      END IF
C
C *** OMEGA DATA ***
C
      write(*,*)
      WRITE(*,*)'Enter File name for omega2x-in.csv'
      read(*,'(a20)')infile
      OPEN (UNIT=7,FILE=infile,STATUS='OLD')
c      OPEN (UNIT=7,FILE='omega2009-ex1.csv',STATUS='OLD')
C
      pi2=datan(1.d0)*8.d0
      read(7,*)
      read(7,*)
      read(7,*)dt,nt
      write(*,*)
      write(*,*)'delta time (s)and number of times =',dt,nt
      read(7,*)
      read(7,*)period
      write(*,*)'minimum period (s) =',period
      duration=dt*nt
      dom=pi2/duration
      iom=1
      nom=int(duration/period)

      WRITE(10,*)'FROM Indat:'
      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*2).GT.MOM)THEN
        write(10,*)'from Indat:'
        WRITE(10,*)' NOM*2 >M OM !!!!!!!!!!!!!'
        write(*,*)'from Indat:'
        WRITE(*,*)' NOM*2 > MOM !!!!!!!!!!!!!'
        STOP
      END IF
C
C *** MEDIUM DATA ***
C
      read(7,*)
      read(7,*)
      READ(7,*)DNS,VP,QP,FP,VS,QS,FS
C
      WRITE(*,2200)
      WRITE(10,2200)
      WRITE(*,2210)DNS,VP,QP,FP,VS,QS,FS
      WRITE(10,2210)DNS,VP,QP,FP,VS,QS,FS
 2200 FORMAT(//' Density,Vp(real),   Qp  ,  Fp  ,Vs(real),   Qs  ,  Fs')
 2210 FORMAT(F9.3,2(F9.2,F9.3,F7.3))
C
C *** FAULT DATA ***
C
      read(7,*)
      read(7,*)
      READ(7,*)SLENGTH,DWIDTH,NSTR,NDIP
      write(*,*)
      write(*,*)'Fault Model Parameters:'
      write(*,*)'  Fault L and W (m) =',SLENGTH,DWIDTH
      write(*,*)'  Numbers of sub-fault along L and W =',NSTR,NDIP
      NF=NSTR*NDIP
      sNF=dsqrt(dfloat(NF))
      IF(NF.GT.MFL)THEN
        write(10,*)'from Indat:'
        WRITE(10,*)' NF.GT.MFL !!!!!!'
        write(*,*)'from Indat:'
        WRITE(*,*)' NF.GT.MFL !!!!!!'
        STOP
      END IF
c
      fl=slength/nstr
      fw=dwidth/ndip
      read(7,*)
c      READ(7,*)DELAY,STR,DIP,VR,TTd,irs
c      Tfc=1.d0/TTd
      READ(7,*)DELAY,STR,DIP,VR,dtime,idum
      write(*,*)'  Starting Time (s) and Vr (m/s)=',DELAY,Vr
      write(*,*)'  Strike and Dip Angles (deg) =',STR,DIP
c      write(*,*)'  Total Faulting Duration (s) and fc (Hz) =',TTd,Tfc
      write(*,*)'  Average Random Delay Time (s) =',dtime
      str=str*pi2/360.d0
      dip=dip*pi2 /360.d0
      read(7,*)
      READ(7,*)XHYPO,YHYPO,ZHYPO
      write(*,*)'  Hypocenter Location (X,Y,Z: m)   =',XHYPO,YHYPO,ZHYPO
      read(7,*)
      READ(7,*)xflt,yflt,zflt
      write(*,*)'  Fault Origin Location (X,Y,Z: m) =',xflt,yflt,zflt
      read(7,*)
c
      FDsum=0.d0
      do 150 j=1,NF
        READ(7,*)i,FD(j),rak(j),sdrop(j),fc(j),fmax(j),Pn(j)
        rak(j)=rak(j)*pi2/360.d0
        FDsum=FDsum+FD(j)
  150 continue
      FDavg=FDsum/dfloat(NF)
      SM=DNS*VS*VS
      SM0=SM*FDavg*SLENGTH*DWIDTH*1.0D+10
      SM0avg=SM0/dfloat(NF)
      write(*,*)
      write(*,*)'  Average Slip (cm)     =',FDavg*100.d0
      write(*,*)'  Seismic Moment(dyne-cm)=',SM0
      sMw=(dlog10(SM0)-16.1d0)/1.5d0
      sMj=(dlog10(SM0)-17.72d0)/1.17d0
      write(*,*)'  Mw (by Hanks and Kanamori,1979) =',sMw
      write(*,*)'  Mj (by M.Takemura, 1990)        =',sMj
c
      write(*,*)
      sfac=1.d0/3.d0
      sSM0=SM0avg/sNF
      do 170 j=1,NF
        if(sdrop(j).eq.0.d0)goto 170
        fc(j)=(4.9d+6)*(sdrop(j)/sSM0)**sfac*Vs/1000.d0
  170 continue
c
c *** Sub-fault parameters ***
c
c   * Dip direction *
      dcs=fw*dcos(dip)
      cst=dcos(str)
      sst=dsin(str)
c
      iflt=0
      ihyp=1
      dismin=999999999999999.d0
      do 200 i=1,ndip
        dc=dcs*(dfloat(i-1)+0.5d0)
        zc=zflt-dc*dtan(dip)
        if(zc.lt.0.d0)then
          write(*,*)'Fault Depth is less than 0!: Dip Number=',i
c          stop
        end if
c     * Strike direction *
        do 220 j=1,nstr   
          iflt=iflt+1
          sc=fL*(dfloat(j-1)+0.5d0)
          dp(iflt)=zc
          zz=zc-zhypo
          xf(iflt)=xflt+sc*cst+dc*sst
          yf(iflt)=yflt+sc*sst-dc*cst
	  xx=xf(iflt)-xhypo
	  yy=yf(iflt)-yhypo
          dis=dsqrt(xx*xx+yy*yy+zz*zz)
          if(dis.lt.dismin)then
            ihyp=iflt
            dismin=dis
          end if
          td(iflt)=delay+dis/VR
  220   continue
  200 continue
c
      do iflt=1,NF
        if(iflt.ne.ihyp)then
          td(iflt)=td(iflt)+dtime*dabs(gasdev(idum))
        end if
      end do
C
      WRITE(10,2300)NF
      WRITE(10,2305)STR,DIP
      WRITE(10,2315)slength,dwidth
      WRITE(10,2320)FL,FW
      WRITE(10,2330)FDavg,SM0
      WRITE(10,2335)sMw,sMj
 2300 FORMAT(//'Seismic Source Data:'
     *      ,/'  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)
 2330 FORMAT( '  Avgerage Slip (m)=',F12.3,
     *        ',  Seismic Moment (dyne-cm)=',e12.4)
 2335 FORMAT( '  Moment Magnitude =',F12.1,',  JMA Magnitude   =',F12.1)
c
      II=0
      DO 270 I=1,NDIP
        dc=dcs*(dfloat(i-1)+0.5d0)
        zc=zflt-dc*dtan(dip)
        DO 270 J=1,NSTR
          II=II+1
          sc=fL*(dfloat(j-1)+0.5d0)
          zz=zc-zhypo
          xf(ii)=xflt+sc*cst+dc*sst
          yf(ii)=yflt+sc*sst-dc*cst
	  xx=xf(ii)-xhypo
	  yy=yf(ii)-yhypo
          dis=dsqrt(xx*xx+yy*yy+zz*zz)
        WRITE(10,2350)II,I,J
        WRITE(10,2360)FD(ii),rak(ii)
        write(10,2370)XF(II),YF(II),DP(II)
        write(10,2375)sdrop(ii),fc(ii),fmax(ii)
        write(10,2380)TD(II),dis/VR,TD(II)-dis/VR
  270 CONTINUE
 2350 FORMAT(/'   Source Number   =',I4,' (',I2,',',I2,')')
 2360 FORMAT( '   Dislocation (m)=',F12.3,', Rake Ang. (rad)=',F9.5)
 2370 FORMAT( '   Source X, Y, and Z Location (m) =',3f9.1)
 2375 format( '   Stress Drop (bar), fc & fmax (Hz) =',f8.2,2f9.3)
 2380 FORMAT( '   Total Time Delay, r/Vr & delta Tr (s) =',3f9.3)
c
c *** Input and compute Boore's Source Parameters ***
c
      call Boore(MOM,MFL,Frate,rate)
C
C *** Radiation Pattern for S and P Waves at high Frequencies ***
C
      read(7,*)
      read(7,*)
      READ(7,*)fR1,fR2
c      write(*,*)fR1,fR2
      write(10,2410)fR1,fR2
 2410 FORMAT(/' Transient Corners Frequencies (Hz) = ',2E12.4)
C
C *** Selection of F-function, which adjusts the difference of the slip
c     functions between the whole- and sub-faults ***
C
      read(7,*)
      read(7,*)
      READ(7,*)NFfunc,TLc
c      write(*,*)NFfunc,TLc
      WRITE(10,2500)NFfunc,TLc
 2500 FORMAT(/' Choice of F-function and TLc (s)',i5,f8.2)
c
c *** Output of F-function ***
c
      write(*,*)'  F-function (No.1) is stored in F-func.csv'
      open(21,file='F-func.csv')
      write(21,*)'freq.(Hz),',' F-func'
      CPX=(0.D0,1.D0)
      pai=4.d0*datan(1.d0)
      do io=iom,nom
        OM=DOM*DFLOAT(IO)
        freq=om/(2.d0*pai)
        OMTd=OM*TLc
c   *  Irikura's F-function (1987) *
        if(NFfunc.eq.0)then
          sinc=dsin(OMTd)/OMTd
          Ffunc=1.d0+(sNF-1.d0)*sinc*exp(cpx*OMTd)
c   *  Exp-type Slip Velocity Function *
        else
          Ffunc=1.d0+(sNF-1.d0)/(1.d0-cpx*OMTd)
        end if
        write(21,2250)freq,cdabs(Ffunc)
      end do
      close(21)
 2250 format(f12.4,',',e12.4)
C
C *** Selection of Green's Function of Homogeneous Full Space ***
C
      read(7,*)
      read(7,*)
      READ(7,*)NCF
c      write(*,*)NCF
      write(10,2600)
      IF(NCF.EQ.0)WRITE(10,2610)
      IF(NCF.EQ.1)WRITE(10,2620)
      IF(NCF.EQ.2)WRITE(10,2630)
      IF(NCF.EQ.3)WRITE(10,2640)
 2600 FORMAT(/' Choice of Green Function:')
 2610 FORMAT(' Complete Homogeneous Full-Space Green Function')
 2620 FORMAT(' Far-Field P & S waves Term (=1/r Term)')
 2630 FORMAT(' Far-Field S waves Term (=1/r Term)')
 2640 FORMAT(' Far-Field P waves Term (=1/r Term)')
C
C *** CHANGE OF SIGNS FOR IMAGINARY PARTS ***
C
      read(7,*)
      read(7,*)
      READ(7,*)ISIGN
c      write(*,*)ISIGN
      IF(ISIGN.EQ.1)WRITE(10,2700)
      IF(ISIGN.EQ.0)WRITE(10,2710)
 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,*)
      read(7,*)
      READ(7,*)NOB
c      write(*,*)NOB
      WRITE(10,2800)NOB
 2800 FORMAT(//' TOTAL NUMBER OF OBSERVATION POINTS =',I6/)
C
      IF(NOB.GT.MOB)THEN
        WRITE(10,*)' NOB>NMOB !!!!!!!!!!!!!'
        STOP
      END IF
C
      read(7,*)
      DO 300 I=1,NOB
        READ(7,*)j,XO(I),YO(I),ZO(I)
        write(*,*)j,XO(I),YO(I),ZO(I)
        WRITE(10,2820)I,XO(I),YO(I),ZO(I)
  300 CONTINUE
 2820 FORMAT('   OB NUM.=',I5,',     (X,Y,Z) =',3F12.4)
c
      CLOSE(UNIT=7)
c
      RETURN
      END
C
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 *** Calculation of Normalized Boore's Moment Rate Functions
c     (Far Field Displacement) at the all sub-faults ***
c
      SUBROUTINE Boore(MOM,MFL,Frate,rate)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 FX(MOM),FX1(MOM),cran,Frate(MOM,MFL,2)
c      COMPLEX*16 FX(MOM),Frate(MOM,MFL,2),cran1(MOM)
      DIMENSION XR(MOM),bov(mom),bovN(mom),rate(mom,MFL,2),bos(mom)
c      DIMENSION XR(MOM),rate(mom,MFL,2)
      COMMON /FREQ/DOM,IOM,NOM /OBS/XO(51),YO(51),ZO(51),NOB
     *    /time/dt,period /ntime/nt
     *    /FLT/STR,DIP,RAK(203),FL,FW,FD(203),VR,TLc,Pn(203)
     *        ,XF(203),YF(203),TD(203),DP(203),fc(203),fmax(203)
     *    /Radiation/fR1,fR2,fran
     *    /IFLT/NF,NSTR,NDIP,NCF,NFfunc /Nboore/Niter,idum
c
      nc=51
      if(nc.lt.NOB)then
        write(10,*)'from Boore:'
        WRITE(10,*)' NC < NOB !!!!!!!!!!!!!'
        write(*,*)'from Boore:'
        WRITE(*,*)' NC < NOB !!!!!!!!!!!!!'
        STOP
      END IF
c
      nc=203
      if(nc.lt.NF)then
        write(10,*)'from Indat:'
        WRITE(10,*)' NC < NF !!!!!!!!!!!!!'
        write(*,*)'from Indat:'
        WRITE(*,*)' NC < NF !!!!!!!!!!!!!'
        STOP
      END IF
c
c   * Input Data *
c
      read(7,*)
      read(7,*)
      read(7,*)fran,Niter,idum
      write(10,*)
      write(10,*)'From subroutine Boore:'
      write(10,*)'  Upper Limit Freq. for Random Phases (Hz) ', Fran
      write(10,*)'  Number of Iterrations =',Niter
      write(10,*)'  Initial Random Number (1-999) =',idum
c
c   * Initial data
c
      pi=datan(1.d0)*4.d0
      dur=dt*dfloat(nt)
      nfast=nom*2
c
c   * Spectra of Nomalized Moment Rate Function *
c
      write(*,*)'  One Example of Phase Spectra is stored in phase.csv'
      open(21,file='Source(No.1)-phase.csv')
      write(21,*)'freq.(Hz),','phase,','envelope(+),','envelope(-)'
c
      write(*,*)'  Now calculating moment-rates for all the sub-faults!'
      do 100 iflt=1,NF
        write(*,*)' Fault No. =',iflt
c
c   * Boore's Envelope Function *
c     (Tdu: the source duration (see eq.(6) in Boore (1983))
c
        Tdu=1.d0/fc(iflt)
        Tw=2.d0*Tdu
        b=1.2531d0
        c=5.d0*b/Tw
        a=(5.d0*2.718281829d0/Tw)**b
        bov(1)=0.d0
        bovN(1)=0.d0
        Tdu2=Tdu*2.d0
        Tdu3=Tdu*3.d0
        pi=datan(1.d0)*4.d0
        write(*,*)'Boore Envelope Function is stored in boore-env.csv'
        open(20,file='boore-env.csv')
        write(20,*)'time,','boore,','envelope'
        do i=2,nt
          time=dt*(i-1)
          bov(i)=a*(time**b)*dexp(-c*time)
          if(time.gt.Tdu2)then
            if(time.le.Tdu3)then
              bov(i)=bov(i)*dcos(pi*(time-Tdu2)/2.d0/Tdu)
            else
              bov(i)=0.d0
            end if
          end if
          bovN(i)=bov(i)**(1.d0/dfloat(Niter))
          if((bovN(i).lt.1.d-10).or.(Niter.eq.0))bovN(i)=0.d0
          write(20,9999)time,bov(i),bovN(i)
        end do
        close(20)
 9999   format(f12.4,5(',',e15.7))
c
        DO 110 irnd=1,2
c
          DO 120 IO=1,nom+1
            OM=DOM*DFLOAT(IO-1)
            f=om/(2.d0*pi)
c     * Brune's omega-square model *
     	    FFc=f/Fc(iflt)
	    S=1.d0/(1.d0+FFc*FFc)
c     * Fmax by Boore (1983) n=4 *
            FFmax=f/Fmax(iflt)
            P=1.d0/dsqrt(1.d0+FFmax**(2.d0*Pn(iflt)))
c     * Fmax by Kamae et al. (1991) *
c            P=1.d0/(1.d0+FFmax)
            SA=S*P
c     * Double derivatives of Moment Rate Function *
c            FX(io)=SA*cran1(io)*(-om*om)
            bos(io)=SA*om*om
c
c     * Generation of Phases (zero for lower and random for higher freq.)*
            ran=(ran0(idum)-0.5d0)*pi*2.d0
            phsav=pi
            if(f.lt.Fran)then
              amp=f/Fran
              amp=dsin(amp*pi/2.d0)
              amp=amp*amp
              amp=amp*amp
              ran=ran*amp
              phsav=phsav*amp
            end if
            dr=+dcos(ran)
            di=-dsin(ran)
            cran=dcmplx(dr,di)
            if((iflt.eq.1).and.(irnd.eq.1))then
              write(21,*)f,',',ran,',',phsav,',',-phsav
            end if
c
c     * Double Derivative of Moment Rate Function (Far-Field Acceleration) *
            FX(io)=SA*cran*(-om*om)
            FX1(io)=FX(io)
c
  120   CONTINUE
c
c     * Calculating First Waveforms (accelaration) using Boore's envelop funtion *
c
          DO I=2,NFAST/2
            FX(NFAST+2-I)=dCMPLX(dREAL(FX(I)),-DIMAG(FX(I)))
          end do
          CALL dWVFAST(NFAST,FX,MOM,1)
          DO I = 1,NFAST
            XR(I)=dREAL(FX(I))/DUR
            XR(i)=XR(i)*bov(i)
          end do
c
c     * Iterations to obtain Boore's waveform functions *
c
        do 150 iter=1,Niter
c
c       * Correction of Fourier Amplitude *
	  do i=1,nfast
            FX(I)=dCMPLX(XR(I),0.0d0)
	  end do
          CALL dWVFAST(nfast,FX,MOM,-1)
          DO K=1,nom+1
            FX(K)=FX(K)*dt
            FX(K)=FX(K)*cdabs(FX1(K))/cdabs(FX(K))
	  end do
c
c       * Waveform (accelaration) *
          DO I=2,NFAST/2
            FX(NFAST+2-I)=dCMPLX(dREAL(FX(I)),-DIMAG(FX(I)))
          end do
          CALL dWVFAST(NFAST,FX,MOM,1)
          DO 180 I = 1,NFAST
            XR(I)=dREAL(FX(I))/DUR
            XR(i)=XR(i)*bovN(i)
  180     CONTINUE
c
  150   CONTINUE
c
c   * Double derivatives of Momant Rate Function (Far Field acceleration) *
c
          DO I = 1,NFAST
            rate(i,iflt,irnd)=xr(i)
          end do
c
c   * Fourier Amplitude of Double derivatives of Momant Rate Function *
c
	  do i=1,nfast
            FX(I)=dCMPLX(rate(i,iflt,irnd),0.0d0)
	  end do
          CALL dWVFAST(nfast,FX,MOM,-1)
          DO K=1,nom+1
            om=dom*(k-1)
            f=om/(2.d0*pi)
            Frate(K,iflt,irnd)=dconjg(FX(K))*dt
	  end do
c
  110   continue
c
  100 continue
      close(21)
c
c   * Data Output *
c
      write(*,*)
      write(*,*)'  Two sets of the double derivatives of the normalized'
      write(*,*)'  moment rates are stored in'
      write(*,*)'  Moment-Rate-Acc-1.csv and Moment-Rate-Acc-2.csv'
      open(22,file='Moment-Rate-Acc-1.csv')
        write(22,2000)(iflt,iflt=1,NF)
        do i=1,nfast
	  time=dt*(i-1)
          write(22,2001)time,(rate(i,iflt,1),iflt=1,NF)
        end do
      close(22)
      open(22,file='Moment-Rate-Acc-2.csv')
        write(22,2000)(iflt,iflt=1,NF)
        do i=1,nfast
	  time=dt*(i-1)
          write(22,2001)time,(rate(i,iflt,2),iflt=1,NF)
        end do
      close(22)
 2000 format('time(sec)',300(',Sub-Fault No.',I3))
 2001 format(f8.4,300(',',e17.5))
c
      write(*,*)'  Two-sets of Fourier amp. of moment-rates are stored'
      write(*,*)'  in Moment_Rate-Amp-1.csv and Moment_Rate-Amp-2.csv'
      open(23,file='Moment-Rate-AccAmp-1.csv')
        write(23,2100)(iflt,iflt=1,NF)
        DO K=1,nom+1
          om=dom*(k-1)
          f=om/(2.d0*pi)
c        write(23,2101)f,cdabs(fx1(k)),(cdabs(Frate(k,iflt,1)),iflt=1,NF)
        write(23,2101)f,bos(k),(cdabs(Frate(k,iflt,1)),iflt=1,NF)
        end do
      close(23)
      open(23,file='Moment-Rate-AccAmp-2.csv')
        write(23,2100)(iflt,iflt=1,NF)
        DO K=1,nom+1
          om=dom*(k-1)
          f=om/(2.d0*pi)
        write(23,2101)f,bos(k),(cdabs(Frate(k,iflt,2)),iflt=1,NF)
        end do
      close(23)
 2100 format('freq.(Hz),Theory for the last sub-fault',
     *    300(',Sub-Fault No.',I3))
 2101 format(f8.4,300(',',e15.7))
c
      return
      END
c
C
C *** DATA OUTPUT ***
C
      SUBROUTINE OTDAT(MOM,MOB,GX,GY,GZ)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 GX(MOM,MOB),GY(MOM,MOB),GZ(MOM,MOB)
      COMMON /FREQ/DOM,IOM,NOM /OBS/XO(51),YO(51),ZO(51),NOB
     *    /SIGN/ISIGN
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
          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.axyz',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
      RETURN
      END
C
C
C *** COMPUTATION OF GREEN FUNCTIONS of the Homogeneous Full-Space ***
C
      SUBROUTINE HKQVW(MOM,MOB,MFL,Frate,GGX,GGY,GGZ)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VPc,VSc,SM,Frate(MOM,MFL,2),cpx
     *      ,GGX(MOM,MOB),GGY(MOM,MOB),GGZ(MOM,MOB)
     *      ,SMO,SM00,EDT,GX,GY,GZ,GXL,GYL,GZL
     *      ,GXHp,GYHp,GZHp,GRHp,GXHs,GYHs,GZHs
     *      ,GXHsh,GYHsh,GTHsh
     *      ,UX(3,3),UY(3,3),UZ(3,3)
     *      ,Ukij,kai,kair,fair,esr,epr
     *      ,ss,sr,sr2,sr3,pp,pr,pr2,pr3,ps2
c     *      ,GXHsh,GYHsh,GTHsh,GXHsv,GYHsv,GZHsv,GRHsv
      DIMENSION EN(3,3),rd(3)
      COMMON /FREQ/DOM,IOM,NOM /OBS/XO(51),YO(51),ZO(51),NOB
     *    /STRC/DNS,VP,QP,FP,VS,QS,FS /IFLT/NF,NSTR,NDIP,NCF,NFfunc
     *    /FLT/STR,DIP,RAK(203),FL,FW,FD(203),VR,TLc,Pn(203)
     *        ,XF(203),YF(203),TD(203),DP(203),fc(203),fmax(203)
     *    /Radiation/fR1,fR2,fran
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=203
      IF(NF.GE.NC)THEN
        WRITE(10,*)' FROM HKQVW:NF>NC !!!!!!!!!!'
        WRITE(10,*)' NF,NC=',NF,NC
        STOP
      END IF
C
c *** Various Parameters ***
c
      CPX=(0.D0,1.D0)
      pai=4.d0*datan(1.d0)
      CST= DCOS(STR)
      SST= DSIN(STR)
      CDP= DCOS(DIP)
      SDP= DSIN(DIP)
c      CDP2=CDP*CDP
c      SDP2=SDP*SDP
c      CDP4= DCOS(4.d0*DIP)
c      CDP22=(1.d0+CDP4)/2.d0
c      SDP22=(1.d0-CDP4)/2.d0
c
c *** Homogeneous High-Freq. Radiation Coefficients for P and S waves ***
c
      write(10,*)
      FpH=dsqrt(4.d0/15.d0)
      FsH=dsqrt(0.4D0)
      WRITE(10,*)
      write(10,*)'Radiation Coefficients'
      WRITE(10,*)'Homogeneous Coefficients for P and S waves =',FpH,FsH
C
C *** GREEN'S FUNCTION ***
C
      do iob=1,nob
      do io=1,nom
        GGX(IO,IOB)=(0.D0,0.D0)
        GGY(IO,IOB)=(0.D0,0.D0)
        GGZ(IO,IOB)=(0.D0,0.D0)
      end do
      end do
c
c   * Green's function computations *
c
      sNF=dsqrt(dfloat(NF))
      DO 100 IO=IOM,NOM
C
        OM=DOM*DFLOAT(IO)
        freq=om/(2.d0*pai)
        OMTd=OM*TLc
        sinc=dsin(OMTd)/OMTd
        write(*,1005)freq,dfloat(io)/dfloat(nom)*100.d0
 1005   format('   Freq (%) =',f7.3,' (',f5.1,' %),')
c
        DO 200 IOB=1,NOB
C
          DO 300 IS=1,NF
c
            if(FD(IS).eq.0.d0)goto 300
c
            RX=XF(IS)-XO(IOB)
            RY=YF(IS)-YO(IOB)
            RZ=DP(IS)-ZO(IOB)
            RR=DSQRT(RX*RX+RY*RY+RZ*RZ)
            rd(1)=Rx/RR
            rd(2)=Ry/RR
            rd(3)=Rz/RR
            CRK= DCOS(RAK(is))
            SRK= DSIN(RAK(is))
c            CRK2=CRK*CRK
c            SRK2=SRK*SRK
c
            EDT=CDEXP(cpx*OM*TD(IS))
c
c    * Complex VP, VS, and, Moment Rate Function *
            call QPQS(freq,QPe,QSe)
            if(QP.eq.0.d0)then
              VPc=VP*DCMPLX(1.D0,0.d0)
            else
              VPc=VP*DCMPLX(1.D0,-0.5D0/QPe)
            end if
            if(QS.eq.0.d0)then
              VSc=VS*DCMPLX(1.D0,0.d0)
            else
              VSc=VS*DCMPLX(1.D0,-0.5D0/QSe)
            end if
            SM=DNS*VSc*VSc
c    * Seismic Moment *
            SM00=FD(IS)*SM*FL*FW/sNF
c    * Slip Velocity Function -> Slip Displacement Function *
            SMO=SM00*CPX/OM
c
c    * fai, kai, fai,r, and kai,r *
            ss=+cpx*om/VSc
            pp=+cpx*om/VPc
	    sr=ss*rr
	    pr=pp*rr
	    ps2=pp*pp/(ss*ss)
	    esr=cdexp(sr)/rr
	    epr=cdexp(pr)/rr
c
            if(NCF.eq.0)then
c
c    * Complete Green's Functions of Homogeneous Full-Space *
c
	      sr2=sr*sr
	      sr3=sr*sr2
	      pr2=pr*pr
	      pr3=pr*pr2
              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
            else
c
c    * Green's Function of the far field S & P terms *
c
              kai=(0.d0,0.d0)
              if(NCF.eq.1)then
                kair=esr*ss
                kair=kair-ps2*epr*pp
                fair=esr*ss
              else
                if(NCF.eq.2)then
                  kair=esr*ss
                  fair=esr*ss
                else
                  if(NCF.eq.3)then
                    kair=(0.d0,0.d0)
                    kair=-ps2*epr*pp
                    fair=(0.d0,0.d0)
                  end if
                end if
              end if
c
            end if
c
c    * Green's Functions: Uki,j=Uxi,j=UX(i,j) *
              do i=1,3
                D1i=0.d0
	        if(1.eq.i)D1i=1.d0
                D2i=0.d0
	        if(2.eq.i)D2i=1.d0
                D3i=0.d0
	        if(3.eq.i)D3i=1.d0
                do j=1,3
                  D1j=0.d0
	          if(1.eq.j)D1j=1.d0
                  D2j=0.d0
	          if(2.eq.j)D2j=1.d0
                  D3j=0.d0
	          if(3.eq.j)D3j=1.d0
                  Dij=0.d0
	          if(i.eq.j)Dij=1.d0
                    rdij1=rd(i)*rd(j)*rd(1)
                    Ukij=fair*rd(j)*D1i-kair*rdij1
                    Ukij=Ukij-(rd(i)*D1j+rd(1)*Dij-2.d0*rdij1)*kai/RR
                    UX(i,j)=Ukij/(4.d0*pai*SM)
c
                    rdij2=rd(i)*rd(j)*rd(2)
                    Ukij=fair*rd(j)*D2i-kair*rdij2
                    Ukij=Ukij-(rd(i)*D2j+rd(2)*Dij-2.d0*rdij2)*kai/RR
                    UY(i,j)=Ukij/(4.d0*pai*SM)
c
                    rdij3=rd(i)*rd(j)*rd(3)
                    Ukij=fair*rd(j)*D3i-kair*rdij3
                    Ukij=Ukij-(rd(i)*D3j+rd(3)*Dij-2.d0*rdij3)*kai/RR
                    UZ(i,j)=Ukij/(4.d0*pai*SM)
	        end do
	      end do
c
c    * Low Frequency Wave *
c
            XN=-SDP*SST
            YN= SDP*CST
            ZN=-CDP
            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

              GXL=(0.D0,0.D0)
              GYL=(0.D0,0.D0)
              GZL=(0.D0,0.D0)
              DO I=1,3
              DO J=1,3
                GXL=GXL+EN(I,J)*UX(J,I)
                GYL=GYL+EN(I,J)*UY(J,I)
                GZL=GZL+EN(I,J)*UZ(J,I)
              end do
              end do
c
c    * High Frequency Wave *
c
c      * High Frequency P Wave (Explosive Source) *
c
             GXHp=(0.D0,0.D0)
             GYHp=(0.D0,0.D0)
             GZHp=(0.D0,0.D0)
             if(NCF.eq.2)goto 555
c
             do i=1,3
               GXHp=GXHp+UX(i,i)
               GYHp=GYHp+UY(i,i)
               GZHp=GZHp+UZ(i,i)
             end do
c
c          * Unit Vector of the SP-OP line *
             XN=-rd(1)
             YN=-rd(2)
             XY2=dsqrt(XN*XN+YN*YN)
             XNS=XN/XY2
             YNS=YN/XY2
c
c          * Radial Component and Distribution to X and Y components *
             GRHp= XNS*GXHp+YNS*GYHp
             GXHp=FpH*GRHp/dsqrt(2.d0)
             GYHp=GXHp
             GZHp=FpH*GZHp
c
c      * High Frequency SH Wave (Fault Plane is perpendicular to the OP-SP line) *
c
  555        GXHsh=(0.D0,0.D0)
             GYHsh=(0.D0,0.D0)
c             GXHsv=(0.D0,0.D0)
c             GYHsv=(0.D0,0.D0)
c             GZHsv=(0.D0,0.D0)
             if(NCF.eq.3)goto 666
c
c      *
c          * Unit Vector of the Fault Plane *
             XN=-rd(1)
             YN=-rd(2)
             ZN=-rd(3)
             XY2=dsqrt(XN*XN+YN*YN)
             XNS=XN/XY2
             YNS=YN/XY2
C
c      * High Frequency SH Wave *
c          * Unit Slip Vector of the SH Source *
             EX= YNS
             EY=-XNS
             EZ= 0.d0
c
c      * (Ei x Nj + Ej x Ni) and SH wave Green function *
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
              DO I=1,3
              DO J=1,3
                GXHsh=GXHsh+EN(I,J)*UX(J,I)
                GYHsh=GYHsh+EN(I,J)*UY(J,I)
              end do
              end do
c
c      * Transverse Component and Distribution to X and Y components *
             GTHsh=-YNS*GXHsh+XNS*GYHsh
             GXHs=FsH*GTHsh/dsqrt(2.d0)
             GYHs=GXHs
             GZHs=(0.d0,0.d0)
c
c    * Hybrid Wave *
c
  666        if(freq.lt.fR1)then
                GX=GXL
                GY=GYL
                GZ=GZL
             else
                if(freq.lt.fR2)then
                  dfreq=(freq-fR1)/(fR2-fR1)
                  GX=GXL+(GXHp+GXHs-GXL)*dfreq
                  GY=GYL+(GYHp+GYHs-GYL)*dfreq
                  GZ=GZL+(GZHp+GZHs-GZL)*dfreq
                else
                  GX=GXHp+GXHs
                  GY=GYHp+GYHs
                  GZ=GZHp+GZHs
                end if
             end if
c
c    * Irikura's F-function *
            if(NFfunc.eq.0)then
              EDT=EDT*(1.d0+(sNF-1.d0)*sinc*exp(cpx*OMTd))
c    * Exp-type F-function *
            else
              EDT=EDT*(1.d0+(sNF-1.d0)/(1.d0-cpx*OMTd))
            end if
c
            GX=GX*SMO*EDT*Frate(io+1,is,1)
            GY=GY*SMO*EDT*Frate(io+1,is,2)
            GZ=GZ*SMO*EDT*(Frate(io+1,is,1)+Frate(io+1,is,2))/2.d0
            GGX(IO,IOB)=GGX(IO,IOB)+GX
            GGY(IO,IOB)=GGY(IO,IOB)+GY
            GGZ(IO,IOB)=GGZ(IO,IOB)+GZ
c
  300     CONTINUE
C
  200   CONTINUE
C
  100 CONTINUE
c
      RETURN
      END
c
c
c
C
C *** FOURIER TRANSFORM OR INVERSE FOURIER TRANSFORM
C
      SUBROUTINE dWVFAST(N,A,MOM,IND)
C
C *** N   : NUMBER OF FOURIER COEFF. OR WAVE ( = 2**I )
C *** A(I): FOURIER COEFFICIENTS ( N/2+1 ) OR WAVE ( REAL PART )
C    ( FOURIER SPECTRUM = FOURIER COEFFICIENTS * DURATION/2.0 )
C *** MOM : N IN MAIN DIMMENSION
C         : IF MWV > 2**13(=8192), CAN NOT USE !!!
C *** IND : = -1 : TRANS.,   = 1 : INVERSE TRANS.
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 A(MOM),TEMP,THETA
C
C *** CALUCURATION OF 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
c
c
c *** Quality Factors of Vp and Vs ***
c
      subroutine QPQS(freq,QPe,QSe)
c
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON /STRC/DNS,VP,QP,FP,VS,QS,FS

c      if(freq.eq.0.d0)then
c        QPe=QP
c        QSe=QS
c      else
        QPe=QP*freq**FP
        if(QPe.lt.QP)QPe=QP
        QSe=QS*freq**FS
        if(QSe.lt.QS)QSe=QS
c      end if
c
      return
      end

