! 表題  GFD ワークゼミサンプルプログラム1 (2011/01/06)
!
! Authors::   ODAKA Masatsugu
! Version::   $Id: $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2010. All rights reserved.
! License::   
!
!== 説明
!
! 1次元移流方程式
!
!   dU(x,t)/dt + c dU(x,t)/dx = 0  (c = 1.0)
! 
! を, 境界条件は周期境界条件, 初期条件 
!
!   U(x,0) = 1.0 (0 < x < 1.0)
!          = 0.0 (1.0 < x < 5.0)
!
! の下解く. 
! 時間積分法はオイラー法, 空間微分は中心差分で離散化する.
program sampl1_1

  implicit none

!-- 宣言部 --

  integer,parameter :: imax = 30  ! 格子点数
  real(8) :: U_a(imax+1)! 振幅変数 [t=(n+1)Δt]
  real(8) :: U_b(imax+1)! 振幅変数 [t=nΔt]
  real(8) :: x_X(imax+1)! 座標値
  real(8) :: Ca = 1.0d0 ! 移流速度
  real(8) :: dx         ! 空間格子間隔
  real(8) :: dt         ! 時間格子間隔
  real(8) :: Lx = 1.0d1 ! x 方向計算領域
  real(8) :: pi         ! 円周率
  integer :: it         ! DO ループ変数 
  integer :: ix         ! DO ループ変数 
  integer :: nt         ! DO ループステップ数
  integer :: iws = 1    ! dcl 用変数
  logical :: lwait=.true.! dcl 用内部変数

!-- パラメータの設定 --

  pi = 4.0d0*atan(1.0d0)

  dx = Lx/real(imax)
  dt = 0.1d0

  write(*,*) "Curant number is ", Ca*dt/dx

  if ( dx <= Ca * dt ) then
    write(*,*) "[Error] CFL condition is not satisfied."
    write(*,*) "[Error] Program is forcedly terminated."
    stop
  end if

  write(*,*) "Input the number of time step (integer)"
  read(*,*) nt 

!-- 座標値の設定

  x_X = 0.0d0

  do it = 2, imax + 1
    x_X(it) = x_X(it-1) + dx
  end do

!-- 初期値の設定 --
 
  U_a = 0.0d0
  U_b = 0.0d0

!  where (x_X <= 1.0d0)
!    U_b = 1.0d0
!  end where

  U_b = max(0.0d0, cos(2.0d0*pi*x_X/Lx))

  U_b(imax+1) = U_b(1)

!-- DCL 初期化

  call gropn(iws)
  call grfrm
  call grswnd(0.0, real(Lx),  -1.2, 1.2)
  call grsvpt(0.2, 0.8, 0.2, 0.8)
  call grstrn(1)
  call grstrf
  call usdaxs
  call uulin(imax+1, real(x_X), real(U_b))

!-- 時間積分開始 --

  do it = 1, nt

    U_a(1) = U_b(1) - dt*Ca*(U_b(2)-U_b(imax))/dx*0.5d0

    do ix = 2, imax
      U_a(ix) = U_b(ix) - dt*Ca*(U_b(ix+1)-U_b(ix-1))/dx*0.5d0
    end do

!    U_a(imax+1) = U_b(imax+1) - dt*Ca*(U_b(2)-U_b(imax))/dx*0.5d0
    U_a(imax+1) = U_a(1)

    U_b = U_a

  end do 

!-- DCL 終了 
    call grfrm
    call grswnd(0.0, real(Lx),  -1.2, 1.2)
    call grsvpt(0.2, 0.8, 0.2, 0.8)
    call grstrn(1)
    call grstrf
    call usdaxs
    call uulin(imax+1, real(x_X), real(U_b))


  call grcls

end program sampl1_1
