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_acc_nmk
    !*****************************************************************
    !Acceleration Response (time domein)
    !*****************************************************************
    use defpi
    implicit none
    integer::ndata,i,j
    real(8),allocatable::acc(:)       !Acceleration (input data)
    real(8),allocatable::racc(:)      !Response acceleration (output data)
    real(8),allocatable::rvel(:)      !Response velocity (output data)
    real(8),allocatable::rdis(:)      !Response displacement (output data)
    real(8)::hh                       !Damping factor
    real(8)::Tp                       !Period
    real(8)::dt                       !Time increment
    character::fnameR*50              !Input & output file name
    character::dummy*50,linebuff*300
    real(8)::acc1,vel1,dis1,acc2,vel2,dis2
    real(8)::ddt,am,ak,ac,f1,f2,fm

    call getarg(1,dummy);read(dummy,*) hh
    call getarg(2,dummy);read(dummy,*) Tp
    call getarg(3,fnameR)  !Input file name (time history data of acceleration)

    open(11,file=fnameR,status='old')
    read(11,'(a)') linebuff
    read(11,*) dummy,dt
    read(11,*) dummy,ndata
    allocate(acc(1:ndata),racc(1:ndata),rvel(1:ndata),rdis(1:ndata))
    do i=1,ndata
        read(11,*) acc(i)
    end do
    close(11)

    ddt=dt/dble(md)
    am=1.0D0                 !Mass
    ak=4.0D0*pi*pi*am/Tp/Tp  !Stiffness
    ac=2.0D0*hh*sqrt(ak*am)  !Damping constant

    acc1=0.0D0
    vel1=0.0D0
    dis1=0.0D0
    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(i)=acc1+acc(i)
        rvel(i)=vel1
        rdis(i)=dis1
    end do

    write(6,'(a)') trim(linebuff)
    write(linebuff,*) 'dt,',dt; call del_spaces(linebuff);write(6,'(a)') trim(linebuff)
    write(linebuff,*) 'ndata,',ndata; call del_spaces(linebuff);write(6,'(a)') trim(linebuff)
    write(linebuff,*) 'Tp,',Tp; call del_spaces(linebuff);write(6,'(a)') trim(linebuff)
    write(linebuff,*) 'dt,',dt; call del_spaces(linebuff);write(6,'(a)') trim(linebuff)
    write(6,'(a)') 'time(sec),acc(gal),racc(gal),rvel(kine),rdis(cm)'
    do i=1,ndata
        write(linebuff,*) dt*dble(i-1),',',acc(i),',',racc(i),',',rvel(i),',',rdis(i)
        call del_spaces(linebuff)
        write(6,'(a)') trim(linebuff)
    end do

    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 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_acc_nmk
