parameter(nnn=20000)
      dimension ddy(nnn),acc(nnn),vel(nnn),dis(nnn)
      real ddy1(nnn),ddy2(nnn),acc1(nnn),acc2(nnn),vel1(nnn),vel2(nnn)
      real accm(200),velm(200),iv1,iv2,iv
      character fmt*8,dm*1,eqn*6
      data h/0.05/
c     "\twc\eqp.dat"
c     "\era\eq\imvelcns.eq"
c     "\era\eq\imvelcew.eq"
c     "\twc\psi.dat"
c     "\twc\a12.dat"
c     "\twc\spc\be5_imvelc.spc"
      open(1,file='eqp.dat')!■計算する地震動の名前を書いたファイル
      open(7,file='psi.dat')!■提案震度出力ファイル
      open(8,file='a12.dat')!■1-2秒応答出力ファイル
      write(7,'(a25)')'  eqn     Iv1  Iv2  Ip  C'
 1    read(1,'(a6)')eqn
      if(eqn.eq.'000')goto2
      write(6,*)eqn
      open(2,file='\era\eq\'//eqn//'ns.eq')!■地震動データ
      open(3,file='\era\eq\'//eqn//'ew.eq')!■地震動データ
      open(9,file='spc\be5_'//eqn//'.spc') !■２方向ベクトル和スペクトル
      read(2,'(a1)')dm
      read(3,'(a1)')dm
      read(2,'(4f8.0)')dt,stt,edt,fct
      read(3,'(4f8.0)')dt,stt,edt,fct
c     dt: 時間刻み, stt: スタート時刻, edt: 終了時刻, ■
c     fct: 倍率の1/100（0.01だとそのまま、ということ）■
      read(2,'(a8)')fmt! ■地震動データフォーマット(ns)
      read(3,'(a8)')fmt! ■地震動データフォーマット(ew)
      read(2,'(i4)')inum!■地震動データの１行のデータ数(ns)
      read(3,'(i4)')inum!■地震動データの１行のデータ数(ew)
      nn=(edt-stt)/dt
      jj=nn/inum+1
      do j=1,jj
      read(2,fmt)(ddy1(i),i=j*inum-(inum-1),j*inum)
      read(3,fmt)(ddy2(i),i=j*inum-(inum-1),j*inum)
        do i=j*inum-(inum-1),j*inum
        ddy1(i)=ddy1(i)*fct*100.
        ddy2(i)=ddy2(i)*fct*100.
        end do
      end do
      do k=1,200
      t=real(k)/100.
      w=6.283185/t
      call resp(h,w,dt,nn,ddy1,acc1,vel1,dis,nnn,sa1,sv1,sd1)
      call resp(h,w,dt,nn,ddy2,acc2,vel2,dis,nnn,sa2,sv2,sd2)
      accm(k)=0.
      velm(k)=0.
      do m=2,nn
      acc(m)=sqrt(acc1(m)**2+acc2(m)**2)
      vel(m)=sqrt(vel1(m)**2+vel2(m)**2)
      if(acc(m).gt.accm(k))accm(k)=acc(m)
      if(vel(m).gt.velm(k))velm(k)=vel(m)
      end do
      write(9,'(3(1pe10.3))')t,accm(k),velm(k)
      end do!k
c
      v1=0.
c     a1=0.
      do k=100,200
      v1=v1+velm(k)
c     a1=a1+accm(k)
      end do
      v1=v1/101.
c     a1=a1/101.
      write(8,'(a6,1pe10.3)')eqn,v1
      iv1=2.17*log(v1)/log(10.)+1.
      if(iv1.ge.5.5)then
        iv=iv1
        iflg=1
      else
        do k=10,100
        v2=v2+velm(k)
        end do
        v2=v2/91.
        iv2=1.92*log(v2)/log(10.)+2.02
        if(iv2.lt.5.5)then
          v2=0.
          iv=iv2
          iflg=2
        else
          iv=(iv1+iv2)/2.
          iflg=3
        end if
      end if
      write(7,'(a8,3f5.2,i2)')eqn,iv1,iv2,iv,iflg
      goto1
 2    continue
      stop
      end
c
      subroutine resp(h,w,dt,nn,ddy,acc,vel,dis,nd,sa,sv,sd)
      dimension ddy(nd),acc(nd),vel(nd),dis(nd)
      w2=w*w
      hw=h*w
      wd=w*sqrt(1.-h*h)
      wdt=wd*dt
      e=exp(-hw*dt)
      cwdt=cos(wdt)
      swdt=sin(wdt)
      a11=e*(cwdt+hw*swdt/wd)
      a12=e*swdt/wd
      a21=-e*w2*swdt/wd
      a22=e*(cwdt-hw*swdt/wd)
      ss=-hw*swdt-wd*cwdt
      cc=-hw*cwdt+wd*swdt
      s1=(e*ss+wd)/w2
      c1=(e*cc+hw)/w2
      s2=(e*dt*ss+hw*s1+wd*c1)/w2
      c2=(e*dt*cc+hw*c1-wd*s1)/w2
      s3=dt*s1-s2
      c3=dt*c1-c2
      b11=-s2/wdt
      b12=-s3/wdt
      b21=(hw*s2-wd*c2)/wdt
      b22=(hw*s3-wd*c3)/wdt
      acc(1)=2.*h*w*ddy(1)*dt
      vel(1)=-ddy(1)*dt
      dis(1)=0.
      dx=vel(1)
      x=0.
      sa=0.
      sv=0.
      sd=0.
      do 110 m=2,nn
      dxf=dx
      xf=x
c     write(6,*)m,ddy(m)
      ddym=ddy(m)
      ddyf=ddy(m-1)
      x=a12*dxf+a11*xf+b12*ddym+b11*ddyf
      dx=a22*dxf+a21*xf+b22*ddym+b21*ddyf
      ddx=-2.*hw*dx-w2*x
      acc(m)=ddx
      vel(m)=dx
      dis(m)=x
      sa=amax1(sa,abs(ddx))
      sv=amax1(sv,abs(dx))
      sd=amax1(sd,abs(x))
 110  continue
      return
      end