! 1次元移流方程式
! 時間積分: オイラー法
! 周期境界条件
! 空間積分: 中心差分法
program main
implicit none

! --- 変数の宣言 ---
integer , parameter :: nmax = 30     ! 分割数
real    , parameter :: xmax = 10.0e0 ! 0<x<10
real :: dx = xmax / nmax             ! 格子間隔
real :: dt =0.010e0                  ! 時間間隔
real :: cs = 10.0e0                   ! 速度
real , allocatable :: u_a(:)         !  (n+1)\Delta t
real , allocatable :: u_b(:)         !  n\Delta t
real , allocatable :: x_x(:)         ! 座標
integer :: it, ix, ins                 ! do ループ作業用変数他
integer :: nt                        ! ステップ数
! --- 初期条件 ---
allocate ( u_a(1:nmax+1), u_b(1:nmax+1), x_x(1:nmax+1) )



! --- 初期値 ---
do ix=1, nmax+1
   u_a(ix) = 0.0e0
   u_b(ix) = 0.0e0
   x_x(ix) = 0.0e0
enddo
! --- 座標設定 ---
do ix=1, nmax+1
   x_x(ix) = (ix - 1) * dx
enddo
do ix=1, nmax+1
   if (x_x(ix) <= 1.0e0) then
      u_b(ix) = 1.0e0
   else
      u_b(ix) = 0.0e0
   endif
enddo
u_b(nmax+1)=1.0e0
! 

WRITE(*,*) ' WORKSTATION ID (I)  ? ;'
    call sgpwsn
    read (*,*) ins
    call gropn( ins )
!    CALL GROPN( 2 )                                                                         
    call grfrm 
    call grswnd( 0.0e0, 10.0e0, -1.20e0, 1.20e0 )
    call grsvpt(0.2, 0.8, 0.2, 0.8)
    call grstrf

    call ussttl('x', ' ', 'p', ' ')
    call usdaxs
    call sgplzu(nmax+1, x_x, u_b, 1, 1) 

! --- 計算 ---
write(*,*),'time step input =>'
read(*,*) nt

do it=1, nt  
   u_a(1) = u_b(1) - cs * dt * (u_b(2) - u_b(nmax)) * 0.50e0 / dx
   do ix=2, nmax
      u_a(ix) = u_b(ix) - cs * dt * (u_b(ix+1) - u_b(ix-1)) * 0.50e0 / dx
   enddo
!   u_a(nmax) = u_b(nmax) 
   u_a(nmax+1) = u_b(nmax+1) - cs * dt * (u_b(2) - u_b(nmax)) * 0.50e0 / dx   
   u_b=u_a

! 図を描く
   call grfrm
   call grswnd( 0.0e0, 15.0e0, -1.20e0, 1.20e0 )
   call grsvpt(0.2, 0.8, 0.2, 0.8)
   call grstrf
   
   call ussttl('x', ' ', 'p', ' ')
   call usdaxs
   call sgplzu(nmax+1, x_x, u_b, 1, 1)
!   open(20,file='output.d')
!   write(20,*) x_x, u_b
!   close(20)
enddo
!write(*,*) u_b

call grcls

end program main
