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

	REAL,DIMENSION(xdim)	:: xplot	! 格子点座標値
	REAL,DIMENSION(ydim)	:: yplot	! 格子点座標値
	REAL,DIMENSION(xdim,ydim):: u, v	! 流速(格子点値)
!	REAL,DIMENSION(xdim,ydim):: udif, vdif	! 拡散流速(格子点値)

	REAL,DIMENSION(0:xdim,0:ydim):: TRC	! 変数(格子点値)
	REAL,DIMENSION(xdim,ydim):: TRCP	! 変数(半整数値)

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

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

	xplot = dx * (/ (i, i = 0, xdim-1) /)
	yplot = dy * (/ (i, i = 0, ydim-1) /)

!	u = 0.
!	v = 0.1

	DO i = 1,xdim
           DO j = 1,ydim
	      u(i,j) = - omega * ( ( j - 0.5 )*dy - y0 )
              v(i,j) =   omega * ( ( i - 0.5 )*dx - x0 )
           END DO
 	END DO

!- 初期値の設定(円錐分布)

	DO i = 1, xdim
	   DO j = 1, ydim
	      TRC(i,j) = 4.0*( 1. - SQRT((i + 0.5 - xdim0)**2 + &
       &		                 (j + 0.5 - ydim0)**2)/7.) 
	      TRC(i,j) = MAX( TRC(i,j), 0. )
	   END DO 
	END DO

	CALL bound2( TRC )

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

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

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

!- 初期値の描画

	CALL dennou_set
	CALL get_plotval( TRC, TRCP )
	CALL dennou_draw( xplot, yplot, TRCP, header)
	CALL dennou_draw( xplot, yplot, TRCP, header)
	CALL dennou_draw( xplot, yplot, TRCP, header)
	CALL dennou_draw( xplot, yplot, TRCP, header)

!- 時間積分

	DO i = 1, tstep

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

!- 1 次上流差分

	   CALL get_dtrc_u1_2D( u, v, TRC, DTRC1 )
	   CALL get_dtrc_u1_2D( u, v, TRC+DTRC1, DTRC2 )
!	   CALL get_dtrc1_C2_2D( u, v, TRC, DTRC1 )
!	   CALL get_dtrc1_C2_2D( u, v, TRC+DTRC1, DTRC2 )
	
	   TRC = TRC + ( DTRC1 + DTRC2 )/2

	   CALL bound2( TRC )

	   CALL get_plotval( TRC, TRCP )

	   IF ( MOD( i, 155 ) == 0 ) THEN
	     CALL dennou_draw( xplot, yplot, TRCP, header)
	     CALL dennou_draw( xplot, yplot, TRCP, header)
	     CALL dennou_draw( xplot, yplot, TRCP, header)
	     CALL dennou_draw( xplot, yplot, TRCP, header)
	   END IF

	END DO

	CALL dennou_cls

END PROGRAM main



