module defpi
    implicit none
    real(8),parameter::pi=3.14159265358979323846D0
!    real(8),parameter::beta=1.0D0/4.0D0  !Average acceleration method
    real(8),parameter::beta=1.0D0/6.0D0  !Linear acceleration method
    integer,parameter::md=5
end module defpi


program f90_ars_nmk
    !*****************************************************************
    !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::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))
    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 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 NEWMARK(acc2,vel2,dis2,acc1,vel1,dis1,am,ac,ak,dt,ff)
        !*****************************************************************
        !Newmark Beta Method
        !*****************************************************************
        real(8),intent(out)::acc2,vel2,dis2
        real(8),intent(in)::acc1,vel1,dis1
        real(8),intent(in)::am,ac,ak,dt,ff

        acc2=(ff-ac*(vel1+0.5D0*dt*acc1)-ak*(dis1+dt*vel1+(0.5D0-beta)*dt*dt*acc1))&
             /(am+ac*0.5D0*dt+ak*beta*dt*dt)
        vel2=vel1+0.5D0*dt*(acc1+acc2)
        dis2=dis1+dt*vel1+(0.5D0-beta)*dt*dt*acc1+beta*dt*dt*acc2
    end subroutine NEWMARK


    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
        real(8)::acc2,vel2,dis2
        real(8)::acc1,vel1,dis1
        real(8)::am,ac,ak,f1,f2,fm,ddt
        real(8)::racc,rvel,rdis
        real(8)::acmax,vemax,dimax

        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

            acc1=0.0D0
            vel1=0.0D0
            dis1=0.0D0
            acmax=acc1
            vemax=vel1
            dimax=dis1
            do i=1,ndata
                f1=0.0D0
                if(2<=i)f1=-am*acc(i-1)
                f2=-am*acc(i)
                do j=1,md
                    fm=f1+(f2-f1)/dble(md)*dble(j)
                    call NEWMARK(acc2,vel2,dis2,acc1,vel1,dis1,am,ac,ak,ddt,fm)
                    acc1=acc2
                    vel1=vel2
                    dis1=dis2
                end do
                racc=acc2+acc(i)     !Absolute acceleration
                rvel=vel2            !Relative velocity
                rdis=dis2            !Relative displacement
                if(acmax<abs(racc))acmax=abs(racc)
                if(vemax<abs(rvel))vemax=abs(rvel)
                if(dimax<abs(rdis))dimax=abs(rdis)
            end do
            resacc(k)=acmax
            resvel(k)=vemax
            resdis(k)=dimax
        end do
    end subroutine CALRES


    subroutine del_spaces(s)
        character (*), intent (inout) :: s
        character (len=len(s)) tmp
        integer i, j
        j = 1
        do i = 1, len(s)
            if (s(i:i)==' ') cycle
            tmp(j:j) = s(i:i)
            j = j + 1
        end do
        s = tmp(1:j-1)
    end subroutine del_spaces

end program f90_ars_nmk

