C
C *** grfft.f is a FORTRAN code to obtain seismograms in the time domain,
C     which are computed from grpoint.f and grfault.f
C
      PROGRAM grfftsp
C
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(MW=23,ML=44)
      dimension slp35(50,ML,MW),slp125(50,ML,MW)
      dimension slip35(ML,MW),slip125(ML,MW),slip(ML,MW),rake(ML,MW)
      dimension a1(5,ML,MW),a2(5,ML,MW),dis(ML,MW)
c
      NW=22
      NL=43
      OPEN(11,FILE='rateA.dat',STATUS='OLD')
      do 100 iW=1,NW
      do 100 iL=1,NL
        READ(11,*)a1(1,iL,iW),a1(2,iL,iW)
     *           ,a1(3,iL,iW),a1(4,iL,iW),a1(5,iL,iW)
        READ(11,*)a2(1,iL,iW),a2(2,iL,iW)
     *           ,a2(3,iL,iW),a2(4,iL,iW),a2(5,iL,iW)
        do i=0,3
          jj=i*10
          READ(11,*)(slp35(j+jj,iL,iW),j=1,10)
        end do
        READ(11,*)(slp35(j+40,iL,iW),j=1,9)
c
        READ(11,*)
        READ(11,*)
        do i=0,3
          jj=i*10
          READ(11,*)(slp125(j+jj,iL,iW),j=1,10)
        end do
        READ(11,*)(slp125(j+40,iL,iW),j=1,9)
  100 continue
      CLOSE(11)
c
      pi=datan(1.d0)*4.d0
      do 200 iW=1,NW
      do 200 iL=1,NL
        slip35(iL,iW)=0.d0
        slip125(iL,iW)=0.d0
        do k=1,49
          slip35(iL,iW)=slip35(iL,iW)+slp35(k,iL,iW)
          slip125(iL,iW)=slip125(iL,iW)+slp125(k,iL,iW)
        end do
        slip35(iL,iW)=slip35(iL,iW)/2.d0
        slip125(iL,iW)=slip125(iL,iW)/2.d0
        s35=slip35(iL,iW)
        s125=slip125(iL,iW)
        slip(iL,iW)=dsqrt(s35*s35+s125*s125)
        if(s35.ne.0.d0)then
          rake(iL,iW)=datan(s125/s35)*180.d0/pi+35.d0
        else
          rake(iL,iW)=90.d0
        end if
  200 continue
c
      OPEN(12,FILE='rateA-out4.csv',STATUS='unknown')
c
      write(12,*)'Slip (m)'
      write(12,1000)(iL,iL=1,NL)
 1000 format('W/L',100(',L',i2))
      do 300 iW=1,NW
        write(12,1100)iW,(slip(iL,iW),iL=1,NL)
  300 continue
 1100 format('W',i2,100(',',f10.5))
c
      write(12,*)
      write(12,*)'Rake (deg)'
      write(12,1000)(iL,iL=1,NL)
      do 310 iW=1,NW
        write(12,1100)iW,(rake(iL,iW),iL=1,NL)
  310 continue
c
      write(12,*)
      write(12,*)'Slip of Rake = 35 degree (m)'
      write(12,1000)(iL,iL=1,NL)
      do 320 iW=1,NW
        write(12,1100)iW,(slip35(iL,iW),iL=1,NL)
  320 continue
c
      write(12,*)
      write(12,*)'Slip of Rake = 125 degree (m)'
      write(12,1000)(iL,iL=1,NL)
      do 340 iW=1,NW
        write(12,1100)iW,(slip125(iL,iW),iL=1,NL)
  340 continue

      write(12,*)
      write(12,*)'Rupture Time (s)'
      write(12,1000)(iL,iL=1,NL)
      do 360 iW=1,NW
        write(12,1100)iW,(a2(5,iL,iW),iL=1,NL)
  360 continue

      write(12,*)
      write(12,*)'Distance from Epicenter (km)'
      write(12,1000)(iL,iL=1,NL)
      do 380 iW=1,NW
        do iL=1,NL
       dis(iL,iW)=dsqrt(a1(1,iL,iW)*a1(1,iL,iW)+a1(2,iL,iW)*a1(2,iL,iW))
        end do
        write(12,1100)iW,(dis(iL,iW),iL=1,NL)
  380 continue

      write(12,*)
      write(12,*)'Longitude (deg)'
      write(12,1000)(iL,iL=1,NL)
      do 400 iW=1,NW
        write(12,1200)iW,(a1(3,iL,iW),iL=1,NL)
  400 continue

      write(12,*)
      write(12,*)'Latitude (deg)'
      write(12,1000)(iL,iL=1,NL)
      do 420 iW=1,NW
        write(12,1200)iW,(a1(4,iL,iW),iL=1,NL)
  420 continue
 1200 format('W',i2,100(',',f12.7))

      CLOSE(12)
c
      stop
      END

