      parameter(MOB=150,MEV=21,MFQ=17)

      implicit real*8(a-h,o-z)

      dimension s1(MOB,MFQ),f(MFQ),fr(MFQ),ref(1),f0(MEV)
      dimension s2(MOB,MFQ),ra(1),rad(MOB),om(MOB,MFQ)
      character INFILE1(MEV)*100,INFILE2(MEV)*100,INFILE3(MEV)*100,
     &          OTFILE1(MEV)*100,OTFILE2(MEV)*100,OTFILE3(MEV)*100,a*10
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c Event 1
c    moment magnitude = 4.9
c
c    m; seismic moment (dyne*cm)
c    fc; corner frequency (Hz)
c
c    JMA    magnitude      
c       f0=1.79
c    Moment magnitude
        f0(1)=1.56d0

c Event 2
c    moment magnitude = 4.5
c
c    m; seismic moment (dyne*cm)
c    f0; corner frequency (Hz)
c
c    Corner Frequency (Hz)
        f0(2)=2.47d0

c Event 3 
c    moment magnitude = 4.7
c
c    m; seismic moment (dyne*cm)
c    f0; corner frequency (Hz)
c 
c    Corner Frequency (Hz)      
        f0(3)=1.96d0

c Event 4 
c    moment magnitude = 5.3
c
c    m; seismic moment (dyne*cm)
c    f0; corner frequency (Hz)
c 
c    Corner Frequency
        f0(4)=1.12d0

c Event 5
c    moment magnitude = 4.9
c
c    m; seismic moment (dyne*cm)
c    f0; corner frequency (Hz)
c 
c    Corner Frequency      
        f0(5)=1.56d0

c Event 6 
c    moment magnitude = 4.6
c
c    m; seismic moment (dyne*cm)
c    f0; corner frequency (Hz)
c 
c    Corner Frequency      
        f0(6)=2.21d0

c Event 7
c    moment magnitude = 4.4
c
c    m; seismic moment (dyne*cm)
c    f0; corner frequency (Hz)
c 
c    Corner Frequency
        f0(7)=2.78d0

c Event 8
c    moment magnitude = 4.6
c
c    m; seismic moment (dyne*cm)
c    f0; corner frequency (Hz)
c
c    Corner Frequency
        f0(8)=2.21d0

c Event 9
c    moment magnitude = 4.6
c
c    m; seismic moment (dyne*cm)
c    f0;corner frequency
c
c    Corner Frequency
        f0(9)=2.21d0

c Event 10
c    moment magnitude = 4.0
c
c    m; seismic moment (dyne*cm)
c    f0;corner frequency
c
c    Corner Frequency
        f0(10)=4.40d0

c Event 11
c    moment magnitude = 5.3
c
c    m; seismic moment (dyne*cm)
c    f0;corner frequency
c
c    Corner Frequency
        f0(11)=0.98d0

c Event 17 (in Theisis, 12)
c    moment magnitude = 5.1
c
c    m; seismic moment (dyne*cm)
c    f0;corner frequency
c
c    Corner Frequency      
        f0(17)=1.24d0

c Event 18 (13)
c    moment magnitude = 4.7
c
c    m; seismic moment (dyne*cm)
c    f0;corner frequency
c
c    Corner Frequency      
        f0(18)=1.97d0

c Event 19
c    moment magnitude = 6.1
c
c    m; seismic moment (dyne*cm)
c    f0;corner frequency
c
c    Corner Frequency
        f0(19)=0.39d0

c Event 20
c    moment magnitude = 4.6
c
c    m; seismic moment (dyne*cm)
c    f0;corner frequency
c
c    Corner Frequency
        f0(20)=2.21d0

c Event 21
c    moment magnitude = 4.8
c
c    m; seismic moment (dyne*cm)
c    f0;corner frequency
c
c    Corner Frequency
        f0(21)=1.75d0

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      NEV=MEV

      Vs=4.5d0

      beta=Vs*1000.d0*100.d0

      ro=3.4d0

      PI=dacos(-1.d0)

      NDATA=2048

      DT=0.005d0

      df=1.d0/(dfloat(NDATA)*DT)

      a='0123456789'

      open(17,file='omega.txt')
      open(18,file='mo16.inp')

      do 100 iEV=1,NEV

cccc ループから外すもの（対象外） cccc
         if (iEV.eq.12) then
            goto 100
         else if (iEV.eq.13) then
            goto 100
         else if (iEV.eq.14) then
            goto 100
         else if (iEV.eq.15) then
            goto 100
         else if (iEV.eq.16) then
            goto 100
         endif
cccccccccccccccccccccccccccccccccccccc

         if (iEV.lt.10) then
            INFILE1(iEV)='f1.txt'
            INFILE2(iEV)='ref.inp'
            INFILE3(iEV)='rad'//a(iEV+1:iEV+1)//'.inp'
            OTFILE1(iEV)='mom'//a(iEV+1:iEV+1)//'.out'
            OTFILE2(iEV)='momout'//a(iEV+1:iEV+1)
            OTFILE3(iEV)='omega'//a(iEV+1:iEV+1)
         else if (iEV.ge.10.and.iEV.lt.20) then
            INFILE1(iEV)='f1.txt'
            INFILE2(iEV)='ref.inp'
            INFILE3(iEV)='rad1'//a(iEV-9:iEV-9)//'.inp'
            OTFILE1(iEV)='mom1'//a(iEV-9:iEV-9)//'.out'
            OTFILE2(iEV)='momout1'//a(iEV-9:iEV-9)
            OTFILE3(iEV)='omega1'//a(iEV-9:iEV-9)
         else if (iEV.ge.20.and.iEV.lt.30) then
            INFILE1(iEV)='f1.txt'
            INFILE2(iEV)='ref.inp'
            INFILE3(iEV)='rad2'//a(iEV-19:iEV-19)//'.inp'
            OTFILE1(iEV)='mom2'//a(iEV-19:iEV-19)//'.out'
            OTFILE2(iEV)='momout2'//a(iEV-19:iEV-19)
            OTFILE3(iEV)='omega2'//a(iEV-19:iEV-19)
         endif
         call CountNumberOfChracters(INFILE1,iEV,1,iend1,MEV)
         call CountNumberOfChracters(INFILE2,iEV,1,iend6,MEV)
         call CountNumberOfChracters(INFILE3,iEV,1,iend2,MEV)
         call CountNumberOfChracters(OTFILE1,iEV,1,iend3,MEV)
         call CountNumberOfChracters(OTFILE2,iEV,1,iend4,MEV)
         call CountNumberOfChracters(OTFILE3,iEV,1,iend5,MEV)

         open(10,file=INFILE1(iEV)(1:iend1),status='old')
         open(11,file=INFILE2(iEV)(1:iend6),status='old')
         open(12,file=INFILE3(iEV)(1:iend2),status='old')
         open(14,file=OTFILE1(iEV)(1:iend3))
         open(15,file=OTFILE2(iEV)(1:iend4))
         open(16,file=OTFILE3(iEV)(1:iend5))



      read(10,*)(fr(j),j=1,MFQ)
      close(10)

      read(11,*)ref(1)
      read(11,*)ra(1)
      close(11)

      sou=4.d0*PI*ro*(beta**3)*ref(1)*(1.d0*10.d0**5)
c      sou=4.d0*PI*ro*(Vs**3)*ref(1)*(1.d0*10.d0**5)
c      sou=4.d0*PI*ro*(Vs**3)*ref(1)

      i=1
3     read(12,*,end=2)rad(i)
      i=i+1
      goto 3
2     close(12)
      N=i-1


      do 1 j=1,MFQ
              f(j)=df*fr(j)
 1        continue
       do 11 k=1,MFQ
       do 12 l=1,N
       S1(l,k)=-dlog10((ra(1)*((2.d0*PI*f(k))**2))/(sou*(1.d0+
     *    (f(k)/f0(iEV))**2)))


 12    continue
 11    continue

       do 21 k=1,MFQ
       do 22 l=1,N
       s2(l,k)=-dlog10((rad(l)*((2.d0*PI*f(k))**2))
     &         /(sou*(1.d0+(f(k)/f0(iEV))**2)))
c      s2(l,k)=-(dlog10((rad(l)*2.d0*dsqrt(2.d0)*((2.d0*PI*f(k))**2))
c     &         /(sou*(1.d0+(f(k)/f0(iEV))**2))))
 22   continue
 21   continue

       do 31 k=1,MFQ
       do 32 l=1,N
       om(l,k)=dlog10((rad(l)*((2.d0*PI*f(k))**2))/(sou*(1.d0+
     *    (f(k)/f0(iEV))**2)))+23.59d0
c       om(l,k)=dlog10((rad(l)*((2.d0*PI*f(k))**2))/(sou*(1.d0+
c     *    (f(k)/f0(iEV))**2)))
 32   continue
 31   continue


       write(6,*)iEV,s2(1,1),ref(1),sou,om(1,1)

       do 130 l=1,N
          write(14,1000)(S1(l,k),k=1,MFQ)
130    continue
       close(14)

       do 110 l=1,N
          write(15,1000)(S2(l,k),k=1,MFQ)
          write(18,1000)(S2(l,k),k=1,MFQ)
110    continue
       close(15)

       do 120 m=1,N
          write(16,1000)(om(m,n),n=1,MFQ)
          write(17,1000)(om(m,n),n=1,MFQ)
120    continue
       close(16)

 1000  format(17f10.5)


100   continue

      close(17)
      close(18)

      stop
      end
      
      

C***********************************************************************
      subroutine CountNumberOfChracters(a,n,istart,iend,MEV)
C***********************************************************************

      implicit real*8(a-h,o-z)

      character a(MEV)*100

      do 100 i=istart,100
         if (a(n)(i:i).eq.',') then
            go to 1
         else if (a(n)(i:i).eq.'') then
            go to 1
         end if
100   continue
1     iend=i-1

      return
      end







