module defpi implicit none real(8),parameter::pi=3.14159265358979323846D0 end module defpi program f90_ARS_RUNGE !***************************************************************** !Acceleration Response Spectrum (time domein) !***************************************************************** use defpi implicit none integer::ndata,i,k,nnn real(8)::accmax,velmax,dismax !Maximum values of acceleration, velocity and displacement real(8),allocatable::acc(:) !Acceleration (input data) real(8),allocatable::vel(:) !Velocity real(8),allocatable::dis(:) !Displacement real(8),allocatable::resacc(:) !Acceleration real(8),allocatable::resvel(:) !Velocity real(8),allocatable::resdis(:) !Displacement real(8),allocatable::Tk(:) !Natural period real(8)::hh !Damping factor real(8)::dt !Time increment character::fnameR*50,fnameW*50 !Input & output file name character::dummy*50,linebuff*300 nnn=200 hh=0.05 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 k=1,nnn Tk(k)=10.0D0**(log10(2.0*dt)+(log10(10.0)-log10(2.0*dt))/dble(nnn)*dble(k-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 CALRES(ndata,nnn,dt,hh,acc,Tk,resacc,resvel,resdis) !Calculation of response open (12, file=fnameW, status='replace') write(12,'(a)') trim(linebuff) write(linebuff,*) 'nt,',nnn call del_spaces(linebuff) write(12,'(a)') trim(linebuff) write(12,'(a)') 'period(sec),acc(gal),vel(kine),dis(cm)' do k=1,nnn write(linebuff,*) Tk(k),',',resacc(k),',',resvel(k),',',resdis(k) call del_spaces(linebuff) write(12,'(a)') trim(linebuff) end do close(12) stop contains subroutine FUNC(am,ac,ak,ff,x,y,dxdt,dydt) !***************************************************************** !Function for Runge-Kutta method !***************************************************************** real(8),intent(in)::am,ac,ak,ff,x,y real(8),intent(out)::dxdt,dydt dxdt=y dydt=ff/am-ac/am*y-ak/am*x end subroutine FUNC subroutine CALRES(ndata,nnn,dt,hh,acc,Tk,resacc,resvel,resdis) !***************************************************************** !Calculation of response !***************************************************************** integer,intent(in)::ndata,nnn real(8),intent(in)::dt,hh real(8),intent(in)::acc(1:ndata),Tk(1:nnn) real(8),intent(out)::resacc(1:nnn),resvel(1:nnn),resdis(1:nnn) integer::i,j,k,md real(8)::am,ac,ak,f1,f2,fm1,fm2,fmm real(8)::racc,rvel,rdis real(8)::acmax,vemax,dimax real(8)::x1,y1,x2,y2,ddt real(8)::kx1,kx2,kx3,kx4 real(8)::ky1,ky2,ky3,ky4 md=5 ddt=dt/dble(md) do k=1,nnn am=1.0D0 !Mass ak=4.0D0*pi*pi*am/Tk(k)/Tk(k) !Stiffness ac=2.0D0*hh*sqrt(ak*am) !Damping constant racc=0.0D0 rvel=0.0D0 rdis=0.0D0 acmax=racc vemax=rvel dimax=rdis do i=1,ndata f1=0.0D0 if(2<=i)f1=-am*acc(i-1) f2=-am*acc(i) x1=rdis y1=rvel do j=1,md fm1=f1+(f2-f1)*dble(j-1)/dble(md) fm2=f1+(f2-f1)*dble(j)/dble(md) fmm=0.5D0*(fm1+fm2) call FUNC(am,ac,ak,fm1,x1,y1,kx1,ky1) call FUNC(am,ac,ak,fmm,x1+ddt*kx1/2.0D0,y1+ddt*ky1/2.0D0,kx2,ky2) call FUNC(am,ac,ak,fmm,x1+ddt*kx2/2.0D0,y1+ddt*ky2/2.0D0,kx3,ky3) call FUNC(am,ac,ak,fm2,x1+ddt*kx3,y1+ddt*ky3,kx4,ky4) x2=x1+ddt/6.0D0*(kx1+2.0D0*kx2+2.0D0*kx3+kx4) y2=y1+ddt/6.0D0*(ky1+2.0D0*ky2+2.0D0*ky3+ky4) x1=x2 y1=y2 end do racc=f2/am-ac/am*y2-ak/am*x2 racc=racc+acc(i) !Absolute acceleration rvel=y2 !Relative velocity rdis=x2 !Relative displacement if(acmax