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_blc
    !*****************************************************************
    !Acceleration Base Line Correction
    !*****************************************************************
    use defpi
    implicit none
    integer::ndata,i
    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)::dt                       !Time increment
    character::fnameR*50,fnameW*50    !Input & output file name
    character::dummy*50,linebuff*300

    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)

    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

    open (12, file=fnameW, status='replace')
    write(12,'(a)') trim(linebuff)
    write(12,*) 'dt,',dt
    write(12,*) 'ndata,',ndata
    do i=1,ndata
        write(12,*) acc(i)
    end do
    close(12)

    stop

contains

    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
        real(8)::ddt,a1,a2,am1,am2,vel1,vel2,dis1,dis2
        integer::i,j

        ddt=dt/dble(md)
        accmax=0.0D0
        velmax=0.0D0
        dismax=0.0D0
        vel1=0.0D0
        dis1=0.0D0
        do i=1,nn
            a1=0.0D0
            if(2<=i)a1=acc(i-1)
            a2=acc(i)
            do j=1,md
                am1=a1+(a2-a1)/dble(md)*dble(j-1)
                am2=a1+(a2-a1)/dble(md)*dble(j)
                vel2=vel1+0.5D0*(am1+am2)*ddt
                dis2=dis1+vel1*ddt+(0.5D0-beta)*ddt*ddt*(am1)+beta*ddt*ddt*am2
                vel1=vel2
                dis1=dis2
            end do
            vel(i)=vel2     !ground velocity
            dis(i)=dis2     !ground displacement
            if(accmax<abs(acc(i)))accmax=abs(acc(i))
            if(velmax<abs(vel(i)))velmax=abs(vel(i))
            if(dismax<abs(dis(i)))dismax=abs(dis(i))
        end do
    end subroutine IACC


    subroutine CRAC(nn,dt,acc,vel,dis,accmax,velmax,dismax)
        !*****************************************************************
        !Amendment of base line
        !*****************************************************************
        !nn    : Number of time history data of acceleration
        !dt    : Time increment                    (input data)
        !acc() : Time histori data of acceleration (input & output data)
        !vel() : Time history of velocity          (input data)
        !dis() : Time history data of displacement (input & output data)
        !accmax: Maximum value of acceleration     (input data)
        !velmax: Maximum value of velocity         (input data)
        !dismax: Maximum value of displacement     (input data)

        integer,intent(in)::nn
        real(8),intent(in)::dt,accmax,velmax,dismax
        real(8),intent(in)::vel(1:nn)
        real(8),intent(inout)::dis(1:nn)
        real(8),intent(inout)::acc(1:nn)

        integer::i
        real(8)::tt,t,sum,a0,a1,acmax,coef

        tt=dble(nn-1)*dt
        t=0.0D0
        do i=1,nn
            dis(i)=dis(i)*(3.0D0*tt-2.0D0*t)*t*t
            t=t+dt
        end do
        sum=(dis(1)+dis(nn))/2.0D0
        do i=2,nn-1
            sum=sum+dis(i)
        end do
        sum=sum*dt
        a1=28.0D0/13.0D0/tt/tt*(2.0D0*vel(nn)-15.0D0/tt**5.0D0*sum)
        a0=vel(nn)/tt-a1/2.0D0*tt
        t=0.0D0
        acmax=0.0D0
        do i=1,nn
            acc(i)=acc(i)-a0-a1*t
            if(acmax<abs(acc(i)))acmax=abs(acc(i))
            t=t+dt
        end do
        coef=accmax/acmax
        do i=1,nn
            acc(i)=acc(i)*coef
        end do
    end subroutine CRAC

end program f90_blc

