C
c *** Boore's Moment Rate Function (2004/8/20) ***
c     1) Using Green's funtion of homogeneous full-space, strong ground
c        motions are simulated. This is modified from boore3.f
c       (from grfltHF1.f) and boore2a.f.
c     2) Green's function can be the complete function or the far field
c        function (=1/r terms).
c     3) Boore's envelope function is applied to moment rate (velocity
c        in far field), but not to acceleration -> boore-env.csv
c      (Modified: Actually, commented out in this program)
c     4) Random phases are used at the frequencies higher than Fran (Hz)
c        -> phase.csv
c     5) Iterations are carried out to fit Boore's omega
c        squared model -> fourier.csv
c     6) Accelaration, velocity, displacement (moment rate funtion), and
c        moment function are stored in wave.csv -> wave.csv
c     7) Input data for Boore's source model, fc and fmax, and the other
c        parameters, idum, fran, and Niter are read from input data
c        'boore5.in'.
c     8) Strong motions of sub faults are superposed by amplification
c        factors of one at lower frequencies, and sqrt(NL*NW) at higher
c        frequencies.
c     9) Isotropic radiation patterns are used for S and P waves at
c        high frequencies.
c     10) To obtain wave forms, fft pronrams (e.g., grfftp.f) are needed.
c
c *** This program was modified from boore6.f to correct the reduction of
c     amplitudes due to phase differences caused by Vr and Vs among
c     sub-faults using the following techniques.
c      1) The amplification factor for Vr and Vs are used only at coherent
c        (low) frequencies (f < fran). At random (high) frequency range
c        (f > fran), the amplification factors for numbers of sub-faults
c        (sqrt(NL*NW) are used.
c      2) In order to remove the artificial dominant frequencies, random
c          distributions of rupture times are introduced.
c      3) The amplification factor for Vr and Vs are smoothed out to avoid
c        extremely large amplitudes, by using the moving average method.
c       (Aug., 2005)
c
c   * Parameters :
c
c     MOM : Number of Omega
c     MOB : NUmber of Observation Points
c     MFL : Number of Sub-Faults
c
      PROGRAM boor10
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)
     *          ,VTsum(MOM)
      DIMENSION rate(mom,MFL),Vamp(MOM)
C
      WRITE(*,*)
      WRITE(*,*)'"boor10.f" simulates broad-band ground motions using '
      WRITE(*,*)' the Boore moment rate function and homogenous full'
      write(*,*)' space Green function'
      write(*,*)'  (Originally made in 1992, and modified in 2005)'
      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,VTsum,Vamp,GX,GY,GZ)
      CALL OTDAT(MOM,MOB,GX,GY,GZ)
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 VP,VS,Frate(MOM,MFL)
      DIMENSION rate(mom,MFL)
      COMMON /FREQ/DOM,IOM,NOM /OBS/XO(51),YO(51),ZO(51),NOB
     *    /time/dt,period /ntime/nt
     *    /STRC/DNS,VP,VS  /Hypo/xhypo,yhypo,zhypo
     *    /FLT/STR,DIP,RAK(203),FL,FW,FD(203),VR
     *        ,XF(203),YF(203),TD(203),DP(203)
     *    /IFLT/NF,NSTR,NDIP,NCF /SIGN/ISIGN
     *    /Radiation/FsH,FpH,fR1,fR2 /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(*,*)'"boor10.in" is the input data'
      write(*,*)
      OPEN (UNIT=7,FILE='boor10.in',STATUS='OLD')
C
      pi2=datan(1.d0)*8.d0
      read(7,1100)dt,nt
      read(7,1110)period
 1100 format(/f8.0,i8)
 1110 format(f8.0)
      duration=dt*dfloat(nt)
      dom=pi2/duration
      iom=1
      nom=int(duration/period)

      WRITE(10,*)'FROM boore7.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,*)'from Indat:'
        WRITE(10,*)' NOM>MOM !!!!!!!!!!!!!'
        write(*,*)'from Indat:'
        WRITE(*,*)' 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
      write(*,*)'SLENGTH,DWIDTH',SLENGTH,DWIDTH
      READ(7,1310)NSTR,NDIP
      write(*,*)'NSTR,NDIP',NSTR,NDIP
      NF=NSTR*NDIP
      IF(NF.GT.MFL)THEN
        write(10,*)'from Indat:'
        WRITE(10,*)' NF.GT.MFL !!!!!!'
        write(*,*)'from Indat:'
        WRITE(*,*)' NF.GT.MFL !!!!!!'
        STOP
      END IF
      fL=slength/nstr
      fW=dwidth/ndip
      READ(7,1330)DELAY,VR
      write(*,*)'DELAY,VR',DELAY,VR
      READ(7,1330)XHYPO,YHYPO,ZHYPO
      write(*,*)'XHYPO,YHYPO,ZHYPO',XHYPO,YHYPO,ZHYPO
      READ(7,1330)xflt,yflt,zflt
      write(*,*)'xflt,yflt,zflt',xflt,yflt,zflt
      READ(7,1330)STR,DIP
      write(*,*)'STR,DIP',STR,DIP
      str=str*pi2/360.d0
      dip=dip*pi2 /360.d0
      read(7,*)
      do j=1,NF
        READ(7,*)FD(j)
        write(*,*)'j,FD',j,FD(j)
      end do
      read(7,*)
      do j=1,NF
        READ(7,*)rak(j)
        write(*,*)'j,rak(j)',j,rak(j)
        rak(j)=rak(j)*pi2/360.d0
      end do
 1300 FORMAT(/2f10.0)
 1310 FORMAT(3I5)
 1330 format(3f10.0)
c
c *** Input Boore's Source Parameters ***
c
      call Boore(MOM,MFL,Frate,rate)
c
c *** Sub-fault parameters ***
c
c   * Dip direction *
      dcs=fw*dcos(dip)
      cst=dcos(str)
      sst=dsin(str)
c
      iflt=0
      TLW=dsqrt(fL*fL+fW*fW)/(2.d0*VR)
      do 200 i=1,ndip
c
        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)
          ddt=(0.5d0-ran0(idum))*TLW
          if(dis/VR.lt.dabs(ddt))ddt=0.d0
          td(iflt)=delay+dis/VR+ddt
  220   continue
  200 continue
C
      WRITE(10,2300)NF
      WRITE(10,2305)STR,DIP
      WRITE(10,2315)slength,dwidth
      WRITE(10,2320)FL,FW
      write(10,2325)xflt,yflt,zflt
      WRITE(10,2330)VR
 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)
 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)
c
      II=0
      DO 270 I=1,NDIP
      DO 270 J=1,NSTR
        II=II+1
        WRITE(10,2350)II,I,J
        WRITE(10,2360)FD(ii),rak(ii)
        write(10,2370)XF(II),YF(II),DP(II)
        write(10,2380)TD(II)
  270 CONTINUE
 2350 FORMAT(/'   Source Number   =',I4,' (',I2,',',I2,')')
 2360 FORMAT( '   Dislocation (m)=',6F9.3,', Rake Ang. (rad)=',6F9.5)
 2370 FORMAT( '   Source X, Y, and Z Location (m) =',8f8.1)
 2380 FORMAT( '   Total Time Delay (sec) =',8f8.3)
C
C *** Selection of Green's Function of Homogeneous Full Space ***
C
      READ(7,1400)NCF
      IF(NCF.EQ.1)WRITE(10,2400)
      IF(NCF.EQ.0)WRITE(10,2410)
 1400 FORMAT(/I10)
 2400 FORMAT(/' Complete Homogeneous Full Space Green Function is used')
 2410 FORMAT(/' Far Field Term (=1/r Term) is used')
C
C *** Radiation Pattern for S and P Waves at high Frequencies ***
C
      READ(7,1500)FsH,FpH
      READ(7,1510)fR1,fR2
      WRITE(10,2500)FsH,FpH
      write(10,2510)fR1,fR2
 1500 FORMAT(/2f10.0)
 1510 FORMAT(2f10.0)
 2500 FORMAT(/' Radiation Coefficients for S and P Waves = ',2f8.3)
 2510 FORMAT(/' Transient Corners Frequencies (Hz) = ',2f8.3)
C
C *** CHANGE OF SIGNS FOR IMAGINARY PARTS ***
C
      READ(7,1600)ISIGN
      IF(ISIGN.EQ.1)WRITE(10,2600)
      IF(ISIGN.EQ.0)WRITE(10,2610)
 1600 FORMAT(/I10)
 2600 FORMAT(/' SIGNS OF IMAGINARY PARTS FOR RESULTS ARE CHANGED')
 2610 FORMAT(/' SIGNS OF IMAGINARY PARTS FOR RESULTS ARE NOT CHANGED')
C
C *** OP DATA ***
C
      READ(7,1700)NOB
      WRITE(10,2700)NOB
 1700 FORMAT(/I10)
 2700 FORMAT(//' TOTAL NUMBER OF OBSERVATION POINTS =',I6/)
C
      IF(NOB.GT.MOB)THEN
        write(*,*)' from Indata:'
        WRITE(10,*)' NOB>MOB !!!!!!!!!!!!!'
        STOP
      END IF
C
      DO 300 I=1,NOB
        READ(7,1720)XO(I),YO(I),ZO(I)
        WRITE(10,2720)I,XO(I),YO(I),ZO(I)
  300 CONTINUE
 1720 FORMAT(3F10.0)
 2720 FORMAT('   OB NUM.=',I5,',     (X,Y,Z) =',3F12.4)
c
      CLOSE(UNIT=7)
c
      RETURN
      END
C
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 VP,VS,FX(MOM),FX1(MOM),cran,Frate(MOM,MFL)
      DIMENSION XR(MOM),bov(mom),bovSQ(mom),rate(mom,MFL)
      COMMON /FREQ/DOM,IOM,NOM /OBS/XO(51),YO(51),ZO(51),NOB
     *    /time/dt,period /ntime/nt /STRC/DNS,VP,VS
     *    /FLT/STR,DIP,RAK(203),FL,FW,FD(203),VR
     *        ,XF(203),YF(203),TD(203),DP(203)
     *    /IFLT/NF,NSTR,NDIP,NCF
     *    /FBoore/fc, fmax, fran /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,*)fc,fmax,fran
      write(*,*)
      write(*,*)'Corner Frequency fc (Hz) =',fc
      write(*,*)'fmax(Hz) =',fmax
      write(*,*)'Upper Limit Freq. for Random Phases (Hz) ', Fran
      read(7,*)Niter,idum
      write(*,*)'Number of Iterrations =',Niter
      write(*,*)'Initial Random Number (1-999) =',idum
c      write(*,*)'Skip Rising and Tapering of Waveform? (yes=1)'
c      read(*,*)iskip
c
c   * Initial data
c
      pi=datan(1.d0)*4.d0
      dur=dt*dfloat(nt)
      nfast=nom*2
      nc=int(2.d0*pi*fc/dom)
c
c   * Boore's Envelope Function *
c
      Tdu=1/fc
      Tw=2.d0*Tdu
      write(*,*)'Tdu,Tw (sec)= ',Tdu,Tw
c      b=1.2531d0
c      c=5.d0*b/Tw
c      a=(5.d0*2.718281829d0/Tw)**b
c      bov(1)=0.d0
c      bovSQ(1)=0.d0
c      write(*,*)'Boore Envelope Function is stored in boore-env.csv'
c      open(20,file='boore-env.csv')
c      write(20,*)'time,','envelope,','sqrt(envelope)'
c      do i=2,nt
c        time=dt*(i-1)
c        bov(i)=a*(time**b)*dexp(-c*time)
c        bovSQ(i)=dsqrt(bov(i))
c        write(20,9999)time,bov(i),bovSQ(i)
c      end do
c      close(20)
 9999 format(f12.4,5(',',e12.4))
c
c   * Spectra of Nomalized Moment Rate Function *
c
      write(*,*)'One Example of Phase Spectra is stored in phase.csv'
      open(21,file='phase.csv')
      write(21,*)'frequency,','phase,','average phase'
      do 100 iflt=1,NF
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
	    S=1.d0/(1.d0+FFc*FFc)
c     * Fmax by Boore (1983) *
            FFmax=f/Fmax
            P=1.d0/dsqrt(1.d0+FFmax*FFmax)
c     * Fmax by Kamae et al. (1991) *
c            P=1.d0/(1.d0+FFmax)
            SA=S*P
c
c     * Random Phases *
            ran=ran0(idum)*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)write(21,9999)f,ran,phsav
c            write(21,9999)f,ran,phsav
c
c     * Moment Rate Function (Displacement in Far Field) *
            FX(io)=SA*cdexp(dcmplx(0.d0,-om/fc))*cran
            FX1(io)=FX(io)
  120   CONTINUE
c
c     * Waveform considering Boore's accelaration envelope function *
c       : Applied to Moment Rate (Displacement) Function
          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
c            XR(i)=XR(i)*bovSQ(i)
          end do
c
c     * Iterations to obtain Boore's waveform functions *
        do 150 iter=1,Niter
	  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
          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
	    time=dt*(i-1)
c      write(32,*)iflt,iter,time,xr(i)
c      if(iskip.eq.1)goto 180
            tmax=dfloat(iter)/(2.d0*fc*dfloat(Niter))
            if(time.lt.tmax)then
              XR(i)=XR(i)*dsin(time*pi/(2.d0*tmax))
            else
              if(time.gt.Tdu)then
                 if(time.lt.dfloat(iter+1)*Tdu)then
                   rnum=(time-Tdu)/(dfloat(iter)*Tdu)
                   XR(i)=XR(i)*dcos(rnum*pi/2.d0)
c      write(30,*)iter,time,rnum,dcos(rnum*pi/2.d0)
                 else
                   XR(i)=0.d0
                 end if
              end if
            end if
c      write(33,*)iflt,iter,time,xr(i)
  180     CONTINUE
  150   CONTINUE
c
c   * Normalized Momant Rate Function (Far Field Displacement) *
          disp=xr(1)
          DO I = 2,NFAST
            disp=disp+(xr(i-1)+xr(i))*dt/2.d0
          end do
          Srate=1.d0/disp
          rate(1,iflt)=xr(1)*Srate
          disp=0.d0
          DO I = 2,NFAST
	    time=dt*(i-1)
            rate(i,iflt)=xr(i)*Srate
            disp=disp+(xr(i-1)+xr(i))*dt/2.d0
          end do
c
c   *  Fourier Amplitude *
	  do i=1,nfast
            FX(I)=dCMPLX(rate(i,iflt),0.0d0)
	  end do
          CALL dWVFAST(nfast,FX,MOM,-1)
c   * Time is delayed by 0.5/fc (sec) *
          DO K=1,nom+1
            om=dom*(k-1)
            f=om/(2.d0*pi)
            Frate(K,iflt)=dconjg(FX(K))*dt
            Frate(K,iflt)=Frate(K,iflt)*cdexp(dcmplx(0.d0,-om/fc/2.d0))
	  end do
c
  100 continue
      close(21)
c
c   * Data Output *
c
      write(*,*)'Normalized Moment Rates are stored in Bwave.csv'
      open(22,file='Bwave.csv')
        do i=1,nfast
	  time=dt*(i-1)
          write(22,2000)time,(rate(i,iflt),iflt=1,NF)
        end do
      close(22)
 2000 format(f8.3,100(',',e12.4))
c
      write(*,*)'Fourier Amplitudes of Slips are stored in Bfour.csv'
      open(23,file='Bfour.csv')
        DO K=1,nom+1
          om=dom*(k-1)
          f=om/(2.d0*pi)
          write(23,2100)f,cdabs(fx1(k)),(cdabs(Frate(k,iflt)),iflt=1,NF)
        end do
      close(23)
 2100 format(f8.3,100(',',e12.4))
c
      return
      END
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 *** 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.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
      RETURN
      END
C
C
C *** COMPUTATION OF GREEN FUNCTIONS of the Homogeneous Full-Space ***
C
      SUBROUTINE HKQVW(MOM,MOB,MFL,Frate,VTsum,Vamp,GGX,GGY,GGZ)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VP,VS,SM,Frate(MOM,MFL),cpx
     *      ,GGX(MOM,MOB),GGY(MOM,MOB),GGZ(MOM,MOB)
     *      ,SMO,SM00,EDT,EDTS,GX,GY,GZ,GXL,GYL,GZL
     *      ,GXHs,GYHs,GZHs,GXHp,GYHp,GZHp
     *      ,UX(3,3),UY(3,3),UZ(3,3),VTsum(MOM)
     *      ,Ukij,kai,kair,fair,esr,epr
     *      ,ss,sr,sr2,sr3,pp,pr,pr2,pr3,ps2
      DIMENSION EN(3,3),rd(3),Vamp(MOM)
      COMMON /FREQ/DOM,IOM,NOM /OBS/XO(51),YO(51),ZO(51),NOB
     *    /STRC/DNS,VP,VS /IFLT/NF,NSTR,NDIP,NCF
     *    /FLT/STR,DIP,RAK(203),FL,FW,FD(203),VR
     *        ,XF(203),YF(203),TD(203),DP(203)
     *    /FBoore/fc, fmax, fran /Nboore/Niter,idum
     *    /Radiation/FsH,FpH,fR1,fR2
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)
      SM=DNS*VS*VS
      CST= DCOS(STR)
      SST= DSIN(STR)
      CDP= DCOS(DIP)
      SDP= DSIN(DIP)
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
      write(*,*)'Amplitude Correction Data are stored in Amp-Data.csv'
      open(21,file='Amp-Data.csv')
      write(21,*)'iob,Freq, Low-Cut Fil,High-Cut Fil, Amplitude Factor'
      DO 100 IOB=1,NOB
C
        DO 200 IS=1,NF
c
          write(*,*)'iob & is =',iob,is
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
C
          CRK= DCOS(RAK(is))
          SRK= DSIN(RAK(is))
c
c    * Moment Rate Function *
          SM00=FD(IS)*SM*FL*FW
c
          DO 300 IO=IOM,NOM
C
            OM=DOM*DFLOAT(IO)
            freq=om/(2.d0*pai)
c            if(mod(io,50).eq.1)write(*,1000)io,om,freq
c 1000       format(' Omega and frequency =',i5,2f10.4)
c
c    * Slip Velocity Function -> Slip Displacement Function *
            SMO=SM00*Frate(io+1,is)*CPX/OM
c
            EDT=CDEXP(cpx*OM*TD(IS))
            EDTS=CDEXP(cpx*OM*(TD(IS)+RR/dreal(Vs)))
            VTsum(io)=VTsum(io)+EDTS
c
c    * fai, kai, fai,r, and kai,r *
            ss=+cpx*om/vs
            pp=+cpx*om/vp
	    sr=ss*rr
	    pr=pp*rr
	    ps2=pp*pp/(ss*ss)
	    esr=cdexp(sr)/rr
	    epr=cdexp(pr)/rr
c
            if(NCF.eq.1)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 1/r terms (far field terms) *
c
              kai=(0.d0,0.d0)
              kair=esr*ss
              kair=kair-ps2*epr*pp
              fair=esr*ss
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      * Check of the P Wave Radiation Pattern *
c
             re=rd(1)*EX+rd(2)*EY+rd(3)*EZ
             rn=rd(1)*XN+rd(2)*YN+rd(3)*ZN
             FpRL=2.d0*re*rn
             if(FpRL.ge.0.d0)then
               FpRH= FpH
             else
               FpRH=-FpH
             end if
c
c      * High Frequency P Wave *
c
             GXHp=(0.D0,0.D0)
             GYHp=(0.D0,0.D0)
             GZHp=(0.D0,0.D0)
             do i=1,3
               GXHp=GXHp+UX(i,i)
               GYHp=GYHp+UY(i,i)
               GZHp=GZHp+UZ(i,i)
             end do
c
             ABGX=cdabs(GXHp)
             ABGY=cdabs(GYHp)
             ABGZ=cdabs(GZHp)
             ABG3=dsqrt(ABGX*ABGX+ABGY*ABGY+ABGZ*ABGZ)/dsqrt(3.d0)
             if(ABGX.lt.ABG3/1000.d0)then
               if(ABGY.lt.ABG3/1000.d0)then
                 GXHp=FpRH*ABG3*GZHp/ABGZ
                 GYHp=GXHp
                 GZHp=GXHp
               else
                 if(ABGZ.lt.ABG3/1000.d0)then
                   GYHp=FpRH*ABG3*GYHp/ABGY
                   GXHp=GYHp
                   GZHp=GYHp
                 else
                   GZHp=FpRH*GZHp
                   GYHp=FpRH*GYHp/dsqrt(2.d0)
                   GXHp=GYHp
                 end if
               end if
             else
               if(ABGY.lt.ABG3/1000.d0)then
                 if(ABGZ.lt.ABG3/1000.d0)then
                   GXHp=FpRH*ABG3*GXHp/ABGX
                   GYHp=GXHp
                   GZHp=GXHp
                 else
                   GZHp=FpRH*GZHp
                   GXHp=FpRH*GXHp/dsqrt(2.d0)
                   GYHp=GXHp
                 end if
               else
                 if(ABGZ.lt.ABG3/1000.d0)then
                   GXHp=FpRH*GXHp
                   GYHp=FpRH*GYHp
                   GZHp=FpRH*ABG3*(GXHp/ABGX+GYHp/ABGY)
                 else
                   GXHp=FpRH*GXHp
                   GYHp=FpRH*GYHp
                   GZHp=FpRH*GZHp
                 end if
               end if
             end if
c
c      * Check of the S Wave Radiation Pattern *
c
             FsXL=re*XN+(EX-2.d0*re*rd(1))*rn
             FsYL=re*YN+(EY-2.d0*re*rd(2))*rn
             FsTL=-FsXL*rd(2)+FsYL*rd(1)
             if(FsTL.ge.0.d0)then
               FsTH= FsH
             else
               FsTH=-FsH
             end if
C
c      * High Frequency S Wave *
c
             XN=-rd(1)
             YN=-rd(2)
             ZN=-rd(3)
             XY2=dsqrt(XN*XN+YN*YN)
c             EX= YN/XY2
c             EY=-XN/XY2
             EX=-YN/XY2
             EY= XN/XY2
             EZ=0.d0
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
             GXHs=(0.D0,0.D0)
             GYHs=(0.D0,0.D0)
             GZHs=(0.D0,0.D0)
              DO I=1,3
              DO J=1,3
                GXHs=GXHs+EN(I,J)*UX(J,I)
                GYHs=GYHs+EN(I,J)*UY(J,I)
                GZHs=GZHs+EN(I,J)*UZ(J,I)
              end do
              end do
c
              ABGX=cdabs(GXHs)
              ABGY=cdabs(GYHs)
              ABG2=dsqrt(ABGX*ABGX+ABGY*ABGY)/dsqrt(2.d0)
              GXHs=FsTH*ABG2*GXHs/ABGX
              GYHs=FsTH*ABG2*GYHs/ABGY
              GZHs=FsTH*GZHs
c
c
c    * Hybrid Wave *
c
              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
            GX=GX*SMO*EDT
            GY=GY*SMO*EDT
            GZ=GZ*SMO*EDT
            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
c *** Amplitude Correction ***
c     1) freq < Fr2: Corrected by Phase differences caused by Vr and Vs
c     2) freq > Fr2: Corrected by sqrt(NW*NL)
c
c   * Smoothing of Amplitude Correction Factors by Phase differences *
        Vamp(iom)=dfloat(NF)/cdabs(VTsum(iom))
        Vamp(iom+1)=dfloat(NF)/cdabs(VTsum(iom+1))
        do io=iom+2,nom-2
          VT1=cdabs(VTsum(io-2))
          VT2=cdabs(VTsum(io-1))
          VT3=cdabs(VTsum(io))
          VT4=cdabs(VTsum(io+1))
          VT5=cdabs(VTsum(io+2))
          VTa=(VT1+VT2+VT3+VT4+VT5)/5.d0
          Vamp(io)=dfloat(NF)/VTa
        end do
        VT1=cdabs(VTsum(nom-2))
        VT2=cdabs(VTsum(nom-1))
        VT3=cdabs(VTsum(nom))
        VTa=(VT1+VT2+VT3)/3.d0
        Vamp(nom-1)=dfloat(NF)/VTa
        VTa=(VT2+VT3)/2.d0
        Vamp(nom)=dfloat(NF)/VTa
c
      Fq2=Fran*0.5d0
c      Fq2=Fran
      write(*,*)'Upper Limit Frequency for Superpositions (Hz) ',Fq2
c      write(*,*)'Enter Upper Limit Frequency for Superpositions (Hz)'
c      read(*,*)Fq2
c
        Sratio=dsqrt(dfloat(NF))
        do io=iom,nom
          OM=DOM*DFLOAT(IO)
          freq=om/(2.d0*pai)
          if(freq.lt.Fq2)then
            ff=freq/Fq2
c            amp=dsin(ff*pai/Fq2)
            amp0=ff
            amp0=amp0*amp0
            amp0=amp0*amp0
            amp1=1.d0-amp0
            Samp=amp0*(Sratio-1.d0)+1.d0
            amp=Samp*amp0+Vamp(io)*amp1
          else
            amp=Sratio
          end if
       write(21,1000)iob,freq,amp0,amp1,amp
          GGX(IO,IOB)=GGX(IO,IOB)*amp
          GGY(IO,IOB)=GGY(IO,IOB)*amp
          GGZ(IO,IOB)=GGZ(IO,IOB)*amp
        end do
C
  100 CONTINUE
c
      close(21)
 1000 format(i3,6(',',e12.4))
C
c      do iob=1,nob
c      do iom=1,nom
c      write(88,*)iob,dom*(iom-1)/(2*pai),cdabs(GGX(iom,iob))
c      end do
c      end do
c
      RETURN
      END
c
c
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
