parameter(nnn=100000)
      dimension ddy(nnn),acc(nnn),vel(nnn),dis(nnn),vels(nnn)
      real ddy1(nnn),ddy2(nnn),acc1(nnn),acc2(nnn),vel1(nnn),vel2(nnn)
      real ddy3(nnn),acc3(nnn),vel3(nnn)
      real velm(200),velms(200),iv1,iv2,iv3,iv
      character fmt*8,dm*1,eqn*16,FN1*20,FN2*20,FN3*20
      data h/0.05/
      data a1/1.936/b1/2.011/a2/2.030/b2/1.251/a3/2.171/b3/1.002/
c     data a1/1.930/b1/2.023/a2/2.030/b2/1.251/a3/2.171/b3/1.002/
      data bs/5.0/us/5.5/bl/5.5/ul/6.0/
c     "\twc\eqk.dat"
c     "\twc\ksi.dat"
      open(1,file='eqk.dat')!■計算する地震動の名前を書いたファイル
      open(7,file='ksi.dat')!■提案震度出力ファイル!
      write(7,*)'eqn             0.1-1 0.5-1   1-2   Ip'
 1    read(1,'(a16,x,i1)')eqn,j2
      write(6,*)eqn
      if(eqn.eq.'000')goto2
      if(j2.ne.2)then
      FN1=eqn//'.ns'
      FN2=eqn//'.ew'
      FN3=eqn//'.ud'
      else
      FN1=eqn//'.ns2'
      FN2=eqn//'.ew2'
      FN3=eqn//'.ud2'
      end if
      call KNET_IO(FN1,dt,ddy1,nnn,nn)
      call KNET_IO(FN2,dt,ddy2,nnn,nn)
      call KNET_IO(FN3,dt,ddy3,nnn,nn)!3
c
      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)
      call resp(h,w,dt,nn,ddy3,acc3,vel3,dis,nnn,sa3,sv3,sd3)!3
      velm(k)=0.
      velms(k)=0.!3
      do m=2,nn
      vel(m)=sqrt(vel1(m)**2+vel2(m)**2)
      vels(m)=sqrt(vel1(m)**2+vel2(m)**2+vel3(m)**2)!3
      if(vel(m).gt.velm(k))velm(k)=vel(m)
      if(vels(m).gt.velms(k))velms(k)=vels(m)!3
      end do
      end do!k
c     ■短周期（0.1-1秒）
      v1=0.
      do k=10,100
      v1=v1+velms(k)!3
      end do
      v1=v1/91.
      iv1=a1*log(v1)/log(10.)+b1
c     ■中周期（0.5-1秒）
      v2=0.
      do k=50,100
      v2=v2+velm(k)
      end do
      v2=v2/51.
      iv2=a2*log(v2)/log(10.)+b2
c     ■長周期（1-2秒）
      v3=0.
      do k=100,200
      v3=v3+velm(k)
      end do
      v3=v3/101.
      iv3=a3*log(v3)/log(10.)+b3
c     ■ここから
      if(iv3.ge.ul)then
       iv=iv3
       ic=1
       goto 3
      end if
      if(iv3.ge.bl.and.iv3.lt.ul)then
       iv=(iv3*(iv3-bl)+iv2*(ul-iv3))*2!■比例配分の領域の幅が0.5なので*2
       ic=11
       goto 3
      end if
c
      if(iv1.le.bs)then
       iv=iv1
       ic=2
       goto 3
      end if
      if(iv1.gt.bs.and.iv1.le.us)then
       iv=(iv1*(us-iv1)+iv2*(iv1-bs))*2!■比例配分の領域の幅が0.5なので*2
       ic=21
       goto 3
      end if
      iv=iv2
      ic=3
c     ■ここまで
 3    write(7,'(a16,4f6.2,i3)')eqn,iv1,iv2,iv3,iv,ic
      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
c
	SUBROUTINE KNET_IO(FN,DT,D,NDIM,N)
	DIMENSION ND(NDIM),D(NDIM)
	CHARACTER FN*20,SF*54,dm*1

	OPEN(100,FILE='keq\'//FN)
	do i=1,10
	read(100,'(a1)')dm
	end do
	  READ(100,'(18X,I3   )') ISAMP
	  READ(100,'(18X,F7.0 )') DUT
	read(100,'(a1)')dm
	  READ(100,'(18X,A54  )') SF
	do i=1,3
	read(100,'(a1)')dm
	end do

	  DT=1./REAL(ISAMP)
	  N=NINT(DUT*REAL(ISAMP))
	  OPEN(99,FILE='KNET.TMP')
	    WRITE(99,*)SF(1:INDEX(SF,'(')-1)
	    WRITE(99,*)SF(INDEX(SF,'/')+1:54)
	  CLOSE(99)
	  OPEN(99,FILE='KNET.TMP')
	    READ(99,*)SF1
	    READ(99,*)SF2
	  CLOSE(99)
	  READ(100,'(8I9)')(ND(I),I=1,N)
	  CLOSE(100)
c	write(6,*)sf1/sf2
	  SUM=0.
	  DO I=1,N
	    D(I)=ND(I)*SF1/SF2
	    SUM=SUM+D(I)
	  END DO
	  SUM=SUM/REAL(N)
	  DO I=1,N
	    D(I)=D(I)-SUM
	  END DO
	RETURN
	END