module defpi implicit none real(8),parameter::pi=3.14159265358979323846D0 end module defpi program f90_EQSP !***************************************************************** !Acceleration response spetcrum (by Dr. Yorihiko Ohsaki) !***************************************************************** use defpi implicit none integer::i,ndata,nnn real(8)::dt,h real(8)::accmax,velmax,dismax real(8),allocatable::acc(:),vel(:),dis(:) real(8),allocatable::tk(:),resacc(:),resvel(:),resdis(:) character::fnameR*50,fnameW*50 character::dummy*50,linebuff*300 h=0.05 !Damping factor nnn=200 !Number of devision of period span allocate(tk(1:nnn)) allocate(resacc(1:nnn)) allocate(resvel(1:nnn)) allocate(resdis(1:nnn)) call getarg(1,fnameR) !Input data file name (time history data of acceleration call getarg(2,fnameW) !Output file name open(11,file=fnameR,status='old') read(11,'(a)') linebuff read(11,*) dummy,dt read(11,*) dummy,ndata allocate(acc(1:ndata)) allocate(vel(1:ndata)) allocate(dis(1:ndata)) do i=1,ndata read(11,*) acc(i) end do close(11) do i=1,nnn tk(i)=10.0D0**(log10(2.0D0*dt)+(log10(10.0D0)-log10(2.0D0*dt))/dble(nnn)*dble(i-1)) end do call IACC(ndata,dt,acc,vel,dis,accmax,velmax,dismax) !Integral of acceleration call CRAC(ndata,dt,acc,vel,dis,accmax,velmax,dismax) !Amendment of base line call ERES(ndata,dt,h,acc,nnn,tk,resacc,resvel,resdis) !calculation of response spectrum open (12, file=fnameW, status='replace') write(12,'(a)') trim(linebuff) write(linebuff,*) 'nnn,',nnn call del_spaces(linebuff) write(12,'(a)') trim(linebuff) write(12,'(a)') 'period(sec),acc(gal),vel(kine),dis(cm)' do i=1,nnn write(linebuff,*) tk(i),',',resacc(i),',',resvel(i),',',resdis(i) call del_spaces(linebuff) write(12,'(a)') trim(linebuff) end do close(12) stop contains subroutine ERES(nn,dt,h,acc,nt,tk,resacc,resvel,resdis) !************************************************************* !Calculation of response spectrum !************************************************************* !nn : Number od time history data of acceleration !dt : Time increment (input data) !h : Damping factor (input data) !acc() : Time history data of acceleration (input data) !nt : Number of period span devision (i=1 to nt) (input data) !tk() : Calculated period (input data) (input data) !resacc(): Acceleration response spectrum (output data) !resvel(): Velocity response spectrum (output data) !resdis(): Displacement response spectrum (output data) integer,intent(in)::nn,nt real(8),intent(in)::dt,h real(8),intent(in)::acc(1:nn) real(8),intent(in)::tk(1:nt) real(8),intent(out)::resacc(1:nt),resvel(1:nt),resdis(1:nt) integer::i,m real(8)::accmax,velmax,dismax real(8)::w,w2,hw,wd,wdt,e,cwdt,swdt real(8)::ss,cc,s1,c1,s2,c2,s3,c3 real(8)::A11,A12,A21,A22,B11,B12,B21,B22 real(8)::dxf,xf,ddym,ddyf,x,dx,ddx do i=1,nt w=2.0D0*pi/tk(i) w2=w*w hw=h*w wd=w*sqrt(1.0D0-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 accmax=2.0D0*hw*abs(acc(1))*dt velmax=abs(acc(1))*dt dismax=0.0 dxf=-acc(1)*dt xf=0.0 do m=2,nn ddym=acc(m) ddyf=acc(m-1) x=A12*dxf+A11*xf+B12*ddym+B11*ddyf dx=A22*dxf+A21*xf+B22*ddym+B21*ddyf ddx=-2.0D0*hw*dx-w2*x if(accmax<=abs(ddx))accmax=abs(ddx) if(velmax<=abs(dx))velmax=abs(dx) if(dismax<=abs(x))dismax=abs(x) dxf=dx xf=x end do resacc(i)=accmax resvel(i)=velmax resdis(i)=dismax end do end subroutine ERES subroutine IACC(nn,dt,acc,vel,dis,accmax,velmax,dismax) !***************************************************************** !Integral of acceleration !***************************************************************** !nn : Number of time history data of acceleration !dt : Time increment (input data) !acc() : Time history data of Acceleration (input data) !vel() : Time history data of Velocity (output data) !dis() : Time history data of Displacement (output data) !accmax: Maximum value of acceleration (output data) !velmax: Maximum value of velocity (output data) !dismax: Maximum value of displacement (output data) integer,intent(in)::nn real(8),intent(in)::dt real(8),intent(in)::acc(1:nn) real(8),intent(out)::vel(1:nn),dis(1:nn) real(8),intent(out)::accmax,velmax,dismax integer::i accmax=0.0 velmax=0.0 dismax=0.0 vel(1)=0.0 dis(1)=0.0 do i=2,nn vel(i)=vel(i-1)+(acc(i-1)+acc(i))*dt/2.0D0 dis(i)=dis(i-1)+vel(i-1)*dt+(acc(i-1)/3.0D0+acc(i)/6.0D0)*dt*dt if(accmax