module para
  implicit none
  integer :: i,j,k
  real(8),parameter :: c = 1.d0, &
                     & x0 = 0.d0, &
                     & xm = 1.d0, &
                     & t0 = 0.d0, &
                     & tn = 0.4d0
  integer,parameter :: m = 161, &
                     & n = 65
  real(8) :: dx,dt,pi
end module para
!*********************************************************************
program wave_equation
  use para
  implicit none
  real(8),dimension(0:m,0:n) :: u,v
  real(8),dimension(0:m)     :: x
  real(8),dimension(0:n)     :: t
  
  call calc_param(x,t)
  call bcon(u)
  call icon(x,u,v)
  call calc_nextu(u,v)
  stop
  
end program wave_equation
!*********************************************************************
subroutine calc_param(x,t)
  use para
  implicit none
  real(8),dimension(0:m) :: x
  real(8),dimension(0:n) :: t
  
  x(0) = x0
  x(m) = xm
  t(0) = t0
  t(n) = tn
  
  dx = (xm-x0)/dble(m-1)
  dt = (tn-t0)/dble(n-1)
  
  do i = 1,m
    x(i) = dble(i-1)*dx
  end do
  
  do i = 1,n
    t(i) = dble(i-1)*dt
  end do
  
  return
end subroutine calc_param
!***********boundary condition****************************************
subroutine bcon(u)
  use para
  implicit none
  real(8),dimension(0:m,0:n) :: u
  
  do i = 0,n
     u(0,i) = 0.d0
     u(m,i) = 0.d0
  end do
  
  return
  
end subroutine bcon
!************initial condition****************************************
subroutine icon(x,u,v)
  use para
  implicit none
  real(8),dimension(0:m) :: x
  real(8),dimension(0:m,0:n) :: u,v
  
  pi = 4.d0*datan(1.d0)
  
  do i = 0, m
     if(3.d0/8.d0 <= x(i) .and. x(i) <= 5.d0/8.d0) then
        
        u(i,0) = 0.5d0*cos(8.d0*pi*(x(i)-0.5d0))+0.5d0
        
     else if(x(i) <= 3.d0/8.d0 .and. 5.d0/8.d0 <= x(i)) then
           
        u(i,0) = 0.d0
        
     end if
  end do
  
  do i = 0,m
     v(i,0) = 0.d0
  end do
  
  return
  
end subroutine icon
!*********************************************************************
subroutine calc_nextu(u,v)
  use para
  implicit none
  real(8),dimension(0:m,0:n) :: u,v
  do j = 0,n
    do i = 1,m-1
       u(i,j+1) = u(i,j)+c*dt*v(i,j)
    end do
    
    do i = 1,m-1
       v(i,j+1) = v(i,j)+c*dt/(dx*dx) &
                & *(u(i+1,j+1)-2.d0*u(i,j+1)+u(i-1,j+1))
    end do
  end do
  
  return
end subroutine calc_nextu