!****************module***********************************************
module para
	implicit none
	integer :: i,j,k
	real(8),parameter :: xmax = 1.d0, &
										 & xmin = 0.d0, &
										 & tmax = 2.d0, &
										 & tmin = 0.d0, &
										 & nyu	= 0.01d0
	real(8) :: dx,dt,pi,t
	integer,parameter :: M = 100, &
										 & N = 10000
	real(8),dimension(0:M) :: x
end module para
!****************main*************************************************
program burgers
	use para
	implicit none
	real(8),dimension(0:M) :: u
	
	call calcpara
	
	call icon(u)
	
	call bcon(u)
	
	do k = 1,N
		
		t = dble(k)*dt
		call calc_u(u)
		
	end do
	
	write(*,*) 't=',t,'dx=',dx,'dt=',dt
	
	stop
	
end program burgers
!**********計算パラメータの計算***************************************
subroutine calcpara
	use para
	implicit none
	
	pi = 4.d0*datan(1.d0)
	
	dx = (xmax-xmin)/dble(M)
	
	dt = (tmax-tmin)/dble(N)
	
	do i = 0,M
		x(i) = dble(i)*dx
	end do

	return
	
end subroutine calcpara
!*******Inicial Condition*********************************************
subroutine icon(u)
	use para
	implicit none
	real(8),dimension(0:M) :: u
	
	do i = 0,M
		u(i) = dsin(2.d0*pi*x(i))
	end do
	
	return
	
end subroutine icon
!******Boundary Condition*********************************************
subroutine bcon(u)
	use para
	implicit none
	real(8),dimension(0:M) :: u
	
	u(0) = 0.d0
	u(M) = 0.d0
	
	return
	
end subroutine bcon
!********Calculation u************************************************
subroutine calc_u(u)
	use para
	implicit none
	real(8),dimension(0:M) :: u
	
	do i = 1,M-1
		u(i) = u(i)-dt/(2.d0*dx)*u(i)*(u(i+1)-u(i-1)) &
				& +nyu*dt/(dx*dx)*(u(i+1)-2.d0*u(i)+u(i-1))
	end do

	return
	
end subroutine calc_u