!  移流方程式 ∂ u/∂ t +c∂ u/∂ x =0 の解. 
!  条件:
!    const=0. c, x: 任意. x方向に周期境界. 
!    u の初期条件: 計算領域内に収まるような, 半波長分のsin波.
!    よって, u(x,0) = sin
!  解析解は u=f(x-ct)
!  時間方向: リープフロッグ, 空間方向: 中央差分による解法を用いた
!  (2016/03/22 須藤康平 作成)

program advection_leapfrog
  use dcl
  implicit none
  integer :: n          !時間ステップ数
  integer :: nmax =1000  !時間ステップ数
  integer :: j          !時間ステップ数
  
  real(4) :: c = 0.1e0  !移流速度 
  real(4) :: pi = 4.0e0 *atan(1.0e0)           !円周率  
  complex(4),parameter :: im = (0.0e0, 1.0e0)  !虚数単位 
  real(4) :: u_x(0:100)     !座標
  real(4) :: delta_t = 0.5e0        !時間間隔
  real(4) :: delta_x                !空間格子間隔
  
  complex(4) :: u_before(0:100) = complex(0.0e0, 0.0e0)   !数値解
  complex(4) :: u_now(0:100) = complex(0.0e0, 0.0e0)   !数値解
  complex(4) :: u_after(0:100) = complex(0.0e0, 0.0e0)   !数値解
  complex(4) :: u_kaiseki(0:100) = complex(0.0e0, 0.0e0)   !数値解
  

!  !格子定数を入力
!  write(*,*) 'imax ='
!  read(*,*) imax
!  !空間間隔を計算
!  delta_x = L/imax

  !空間格子間隔
  delta_x = 3.0e0*pi/100 
  !初期値を設定
  u_x(0) = 0.0e0

  !x 座標を指定
  do j=0,100
    u_x(j) = real(j* delta_x)
  enddo

  !初期条件を指定
  do j=1,100 
      u_before(j) = real(sin(u_x(j) -pi))
    if(real(u_before(j)) .lt. 0) then
      u_before(j) =0
    endif
  enddo

  !解析解を求める
  u_kaiseki(:) = u_before(:)

  !dcl で描画する
  call DclOpenGraphics()
  call DclNewFrame
  call DclDrawScaledGraph( u_x, real(u_kaiseki) )
  call DclDrawScaledGraph( u_x, real(u_before) )


  !n=1 ステップ目を計算 (オイラー法を用いる)
  do j=1,100 
    u_now(j) = u_before(j) -c*delta_t/delta_x*(u_before(j)-u_before(j-1)) 

    u_kaiseki(j) = sin(u_x(j) -c*delta_t)
    if(real(u_kaiseki(j)) .lt. 0) then
      u_kaiseki(j) = 0
    endif
  enddo
  

  call DclNewFrame
  call DclDrawScaledGraph( u_x, real(u_kaiseki) )
  call DclDrawScaledGraph( u_x, real(u_now) )



  !ループ計算
  do n = 2,nmax

    !空間座標毎の数値解を求める
    do j=2,99
      u_after(j) = u_before(j) -c*delta_t/delta_x*(u_now(j+1)-u_now(j-1))

    enddo
    
    u_after(100) = u_before(100) -c*delta_t/delta_x*(u_now(1)-u_now(99))
    u_after(1) = u_before(1) -c*delta_t/delta_x*(u_now(2)-u_now(0))

    !境界条件を計算
    u_after(0) = u_after(100)

    !u_now, u_after を次時間ステップへ進める
    u_before(:) = u_now(:)
    u_now(:) = u_after(:)

    call DclNewFrame
    call DclDrawScaledGraph( u_x, real(u_after) )
  enddo

  call DclCloseGraphics
endprogram

