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

	INTEGER,PARAMETER	:: dim  = 64	! 格子点数
	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(xdim)	:: udif		! 拡散流速(格子点値)
	REAL,DIMENSION(0:xdim)	:: TRC1		! 変数(格子点値,2次)
	REAL,DIMENSION(xdim)	:: TRC1P	! 変数(半整数値,2次)
	REAL,DIMENSION(0:xdim)	:: TRC2		! 変数(格子点値,4次)
	REAL,DIMENSION(xdim)	:: TRC2P	! 変数(半整数値,4次)
	REAL,DIMENSION(0:xdim)	:: TRC3		! 変数(格子点値,上流)
	REAL,DIMENSION(xdim)	:: TRC3P	! 変数(半整数値,上流)
	REAL,DIMENSION(0:xdim)	:: TRC4		! 変数(格子点値,上流)
	REAL,DIMENSION(xdim)	:: TRC4P	! 変数(半整数値,上流)

	REAL			:: a		! 初期値用パラメータ
	INTEGER			:: tstep	! 時間ステップ数
	CHARACTER(LEN=9)	:: header	! 時刻ヘッダー
	CHARACTER(LEN=3)	:: time		! 時刻
	REAL,DIMENSION(0:xdim)	:: DTRC1	! 時間積分用配列1
	REAL,DIMENSION(0:xdim)	:: DTRC2	! 時間積分用配列2


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

	dx  = ( xmax - xmin )/dim 
	u   = 1.
	dt  = 0.5 * dx/MAXVAL(u)
	a   = 0.1

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

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

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

	CALL bound( TRC1 )

	TRC2 = TRC1
	TRC3 = TRC1
	TRC4 = TRC1

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

	WRITE(*,*) 'INPUT TIME STEP NUMBER ? (I)'
	READ(*,*) tstep

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

!- 初期値の描画

	CALL dennou_set
	CALL get_mval( TRC1, TRC1P )
	CALL get_mval( TRC2, TRC2P )
	CALL get_mval( TRC3, TRC3P )
	CALL get_mval( TRC4, TRC4P )
	CALL dennou_draw( xplot, TRC1P, TRC2P, TRC3P, TRC4P, header)

!- 時間積分

	DO i = 1, tstep

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

!- 2 次中心差分

	   CALL get_dtrc1( dt, dx, u, TRC1, DTRC1 )
	   CALL get_dtrc1( dt, dx, u, TRC1+DTRC1, DTRC2 )

	   TRC1 = TRC1 + ( DTRC1 + DTRC2 )/2

	   CALL bound( TRC1 )

!- 1 次上流差分

	   CALL get_dtrcu1( dt, dx, u, TRC2, DTRC1 )
	   CALL get_dtrcu1( dt, dx, u, TRC2+DTRC1, DTRC2 )

!	   TRC2 = TRC2 + ( DTRC1 + DTRC2 )/2
	   TRC2 = TRC2 + DTRC1

	   CALL bound( TRC2 )

!- MPDATA (1 回反復)

	   CALL get_dtrcu1( dt, dx, u, TRC3, DTRC1 )
	   CALL get_dtrcu1( dt, dx, u, TRC3+DTRC1, DTRC2 )

!	   TRC3 = TRC3 + ( DTRC1 + DTRC2 )/2
	   TRC3 = TRC3 + DTRC1

	   CALL bound( TRC3 )

	   CALL get_dtrcu_mpdata( dt, dx, u, TRC3, udif, DTRC1 )
	   CALL get_dtrcu_mpdata( dt, dx, u, TRC3+DTRC1, udif, DTRC2 )

!	   TRC3 = TRC3 + ( DTRC1 + DTRC2 )/2
	   TRC3 = TRC3 + DTRC1

	   CALL bound( TRC3 )

!- MPDATA (2 回反復)

	   CALL get_dtrcu1( dt, dx, u, TRC4, DTRC1 )
	   CALL get_dtrcu1( dt, dx, u, TRC4+DTRC1, DTRC2 )

!	   TRC4 = TRC4 + ( DTRC1 + DTRC2 )/2
	   TRC4 = TRC4 + DTRC1

	   CALL bound( TRC4 )

	   CALL get_dtrcu_mpdata( dt, dx, u, TRC4, udif, DTRC1 )
	   CALL get_dtrcu_mpdata( dt, dx, u, TRC4+DTRC1, udif, DTRC2 )

!	   TRC4 = TRC4 + ( DTRC1 + DTRC2 )/2
	   TRC4 = TRC4 + DTRC1

	   CALL get_dtrcu_mpdata( dt, dx, udif, TRC4, udif, DTRC1 )
	   CALL get_dtrcu_mpdata( dt, dx, udif, TRC4+DTRC1, udif, DTRC2 )

!	   TRC4 = TRC4 + ( DTRC1 + DTRC2 )/2
	   TRC4 = TRC4 + DTRC1

	   CALL bound( TRC4 )

	   CALL get_mval( TRC1, TRC1P )
	   CALL get_mval( TRC2, TRC2P )
	   CALL get_mval( TRC3, TRC3P )
	   CALL get_mval( TRC4, TRC4P )

	   CALL dennou_draw( xplot, TRC1P, TRC2P, TRC3P, TRC4P, header)

	END DO

	CALL dennou_cls

END PROGRAM main
