!=========================================================================
!			一次元トレーサー移流モデル
!
! 履歴	1997/10/16	小高正嗣	2次精度
!
!=========================================================================
PROGRAM main

	INTEGER,PARAMETER	:: dim  = 128	! 格子点数
	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(dim)	:: TRC		! 変数(格子点値)
	REAL,DIMENSION(xdim)	:: TRCP		! 変数(半整数値)

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


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

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

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

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

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


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

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

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


!- 初期値の描画

	CALL dennou_set
	CALL get_mval( TRC, TRCP )
	CALL dennou_draw( xplot,  TRCP , header)

!- 時間積分

	DO i = 1, tstep

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

	   CALL get_dtrc( TRC, DTRC1 )
	   CALL get_dtrc( TRC+DTRC1, DTRC2 )

	   TRC = TRC + ( DTRC1 + DTRC2 )/2

	   CALL get_mval( TRC, TRCP )
	   CALL dennou_draw( xplot,  TRCP , header)

	END DO


!- 後始末

!	CALL GRCLS

	CONTAINS

!=========================================================================
SUBROUTINE get_dtrc( TRC, DTRC )

	REAL,DIMENSION(dim),INTENT(in)	:: TRC
	REAL,DIMENSION(dim),INTENT(out)	:: DTRC
	REAL,DIMENSION(xdim)	:: TRCM		! トレーサー変数(半整数値)
	REAL,DIMENSION(xdim)	:: FLUX		! フラックス
	REAL,DIMENSION(dim)	:: ADVE		! 移流項
	INTEGER			:: i

	CALL get_mval( TRC, TRCM )

	FLUX = u * TRCM 

	DO i = 1, dim    
	   ADVE(i) = - ( FLUX(i+1) - FLUX(i) )/dx
	END DO

	DTRC = dt * ADVE
	
END SUBROUTINE get_dtrc
!=========================================================================
SUBROUTINE get_mval( val, mval )

	REAL,DIMENSION(:),INTENT(in)	:: val
	REAL,DIMENSION(:),INTENT(out)	:: mval
	INTEGER				:: i

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

	DO i = 2, dim
	   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 PROGRAM main
