!=========================================================================
! 一次元トレーサー移流モデル
!
! 履歴	1997/10/16	小高正嗣
!	1997/10/22	小高正嗣
!	1997/11/13	小高正嗣
!	1998/01/17	小高正嗣
!	2011/06/27	小高正嗣
! 
!=========================================================================
PROGRAM main
USE cal_module
USE dcl_module
USE dcl_module2

	INTEGER,PARAMETER	:: dim  = 100	! 格子点数
	INTEGER,PARAMETER	:: xdim = dim+1	! 変数次元数
	REAL,PARAMETER		:: xmin = 0.	! 領域下限
	REAL,PARAMETER		:: xmax = 1.0	! 領域上限

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

	REAL,DIMENSION(xdim)	:: xplot	! 格子点座標値
	REAL,DIMENSION(xdim)	:: u		! 流速(格子点値)
	REAL,DIMENSION(0:xdim)	:: TRC1		! 変数(格子点値,上流差分)
	REAL,DIMENSION(xdim)	:: TRC1P	! 変数(半整数値,上流差分)
	REAL,DIMENSION(0:xdim)	:: TRC2		! 変数(格子点値,nonFCT-SHASTA)
	REAL,DIMENSION(xdim)	:: TRC2P	! 変数(半整数値,nonFCT-SHASTA)
	REAL,DIMENSION(0:xdim)	:: TRC3		! 変数(格子点値,FCT-SHASTA)
	REAL,DIMENSION(xdim)	:: TRC3P	! 変数(半整数値,FCT-SHASTA)
	REAL,DIMENSION(0:xdim)	:: INIT		! 変数(格子点値,初期値)
	REAL,DIMENSION(xdim)	:: INITP	! 変数(半整数値,初期値)

	INTEGER			:: prep		!
	REAL			:: a		! 初期値用パラメータ
	INTEGER			:: tstep	! 時間ステップ数
	CHARACTER(LEN=9)	:: theader	! 時刻ヘッダー
	CHARACTER(LEN=3)	:: time		! 時刻
	CHARACTER(LEN=10)	:: header1	! スキームヘッダー1
	CHARACTER(LEN=14)	:: header2	! スキームヘッダー2
	CHARACTER(LEN=10)	:: header3	! スキームヘッダー3
	REAL,DIMENSION(0:xdim)	:: DTRC1	! 時間積分用配列1
	REAL,DIMENSION(xdim+1)	:: AFLUX	! 反拡散フラックス
	REAL,DIMENSION(xdim+1)	:: CFLUX	! 補正フラックス

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

	dx  = ( xmax - xmin )/dim 
	u   = 1.
	dt  = 0.2 * dx/MAXVAL(u)

	DO i = 1, xdim
	   xplot(i) = dx * ( i - 1 )
	END DO

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

!       a = 0.1
!	DO i = 1, dim
!	   TRC1(i) = 1.0*EXP( - ( (i - dim/2.)*dx/a)**2 )
!	END DO

!- 初期値の設定(矩形波)

	TRC1 = 0.5 
	DO i = 1, dim/5
	   TRC1(i) = 2.
	END DO

	CALL bound( TRC1 )

	TRC2 = TRC1
	TRC3 = TRC1

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

	WRITE(*,*) 'INPUT TIME STEP NUMBER ? (I)'
	READ(*,*) tstep
	i = 0
	WRITE(time,'(I3.3)') i
	theader = 'tstep=' // time 

	header1 = 'upstream'
	header2 = 'non-FCT SHASTA'
	header3 = 'FCT SHASTA'

!- 初期値の描画

	CALL dennou_set_2
	CALL get_mval( TRC1, TRC1P )
	CALL get_mval( TRC2, TRC2P )
	CALL get_mval( TRC3, TRC3P )
	CALL dennou_draw_2( xplot, TRC1P, TRC1P, TRC2P, TRC3P, &
	&                  theader, header1, header2, header3 )

!- 時間積分

	DO i = 1, tstep

	   WRITE(time,'(I3.3)') i
	   theader = 'tstep=' // time

!- 1 次上流差分

	   CALL get_dtrc_u1( dt, dx, u, TRC1, DTRC1 )
	   TRC1 = TRC1 +  DTRC1
	   CALL bound( TRC1 )

!- non-FCT SHASTA
          
	   CALL get_dtrc_nonfct_shasta( dt, dx, u, TRC2P )

!- SHASTA

	   CALL get_dtrc_fct_shasta( dt, dx, u, TRC3P, AFLUX, CFLUX )

!- 描画

	   IF ( i == 100 .or. i == 800 ) THEN

	   prep = INT(i*0.2)
	   write(*,*) prep
	   IF ( prep .gt. dim ) THEN
	      prep = prep - dim
	   END IF

	   INIT = 0.5 
	   DO j = prep, prep+dim/5
	      INIT(j) = 2.
	   END DO
	   CALL bound( INIT )

	   CALL get_mval( INIT, INITP )
	   CALL get_mval( TRC1, TRC1P )

	   CALL dennou_draw_2( xplot, INITP, TRC1P, TRC2P, TRC3P, &
	&                     theader, header1, header2, header3 )
	   END IF


	END DO

	CALL dennou_cls

END PROGRAM main
