!=========================================================================
!			一次元トレーサー移流モデル
!
! 履歴	1997/10/16	小高正嗣	2次精度
!
!=========================================================================
module cal_module
public :: get_dtrc,get_mval

  contains
!=========================================================================
subroutine get_dtrc( u, dt, dx, TRC, DTRC )

  real,dimension(:),intent(in)  :: u
  real,intent(in)  :: dt, dx
  real,dimension(:),intent(in)  :: TRC 
  real,dimension(:),intent(out) :: DTRC
  real,dimension(:),allocatable :: TRCM ! トレーサー変数(半整数値)
  real,dimension(:),allocatable :: FLUX ! フラックス
  real,dimension(:),allocatable  :: ADVE ! 移流項
  integer              :: i
  integer              :: dimm

  dimm = size(TRC)

  allocate( TRCM(dimm+1) )
  allocate( FLUX(dimm+1) )
  allocate( ADVE(dimm) )

  call get_mval( TRC, TRCM )

  FLUX = u * TRCM 

  do i = 1, dimm    
     ADVE(i) = - ( FLUX(i+1) - FLUX(i) )/dx
  end do

  DTRC = dt * ADVE

  deallocate( TRCM, FLUX, ADVE )

end subroutine get_dtrc
!=========================================================================
subroutine get_mval( val, mval )

  real,dimension(:),intent(in)  :: val
  real,dimension(:),intent(out) :: mval
  integer                       :: i
  integer                       :: dimm, xdim

  dimm = size(val)
  xdim = size(mval)

  mval(1)    = ( val(1) + val(dimm) )/2

  do i = 2, dimm
     mval(i) = ( val(i) + val(i-1) )/2
  end do

  mval(xdim) = mval(1)

end subroutine get_mval
!=========================================================================
!subroutine dennou_set
!
!call SGLSET( 'LCORNER', .FALSE. )
!call SWISET( 'IPOSX', 50 )
!call SWISET( 'IPOSY', 50 )
!call SWLSET( 'LALT', .TRUE. )
!call SWLSET( 'LWAIT', .FALSE. )
!call GROPN( 1 )
!
!end subroutine dennou_set
!=========================================================================
!subroutine dennou_draw( xplot, yplot, header )
!
!real,dimension(:),intent(in):: xplot
!real,dimension(:),intent(in):: yplot
!CHARACTER(LEN=9),intent(in):: header
!
!call GRFRM
!call USSTTL( 'X', '', '', '' ) 
!call USGRPH( xdim, xplot, yplot )
!
!end subroutine dennou_draw
!=========================================================================
end module cal_module

program main
use cal_module
implicit none

  integer,parameter :: dimm  = 128       !格子点数
  integer,parameter :: xdim = dimm+1     ! 変数次元数
  real,parameter    :: xmin = 0.0        ! 領域下限
  real,parameter    :: xmax = 1.0       ! 領域上限

  real :: dx ! 格子間隔
  real :: dt ! 時間間隔

  real,dimension(xdim) :: xplot ! 格子点座標値
  real,dimension(xdim) :: u     ! 流速(格子点値)
  real,dimension(dimm) :: TRC   ! 変数(格子点値)
  real,dimension(xdim) :: TRCP  ! 変数(半整数値)

  real                 :: a      ! 初期値用パラメータ
  integer              :: tstep
!  character(len=9)     :: header ! 時刻ヘッダー
!  character(len=3)     :: time   ! 時刻
  real,dimension(dimm) :: DTRC1  ! 時間積分用配列1
  real,dimension(dimm) :: DTRC2  ! 時間積分用配列2
  integer              :: i

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

  dx  = ( xmax - xmin )/dimm 
  u   = 1.0
  dt  = 0.8 * dx/maxval(u)
  a   = 0.1

  do i = 1, xdim
     xplot(i) = dx * ( i - 1 )
  end do

!- 初期値の設定(ガウス分布)

  do i = 1, dimm
     TRC(i) = 1.0*exp( - ( (i - dimm/2)*dx/a)**2 )
  end do


!- 時間ステップ数の読み込み

  write(unit=6,fmt=*) "INPUT TIME STEP NUMBER ? (I)"
  read(unit=5,fmt=*) tstep

!  write(unit=time,fmt=*) tstep
!  header = "tstep=" // time 


!- 初期値の描画

!call dennou_set
!call get_mval( TRC, TRCP )
!call dennou_draw( xplot,  TRCP , header)

!- 時間積分

  do i = 1, tstep

!     write(unit=time,fmt=*) i
!     header = "tstep=" // time 
     write(unit=6,fmt=*) i

     call get_dtrc( u, dt, dx, TRC, DTRC1 )
     call get_dtrc( u, dt, dx, TRC+DTRC1, DTRC2 )

     TRC = TRC + ( DTRC1 + DTRC2 )/2

     call get_mval( TRC, TRCP )
!   call dennou_draw( xplot,  TRCP , header)

  end do


!- 後始末

!call GRCLS

end program main

