!=========================================================================
! 二次元トレーサー移流モデル
!
! 履歴	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_1	! 変数1(格子点値)
	REAL,DIMENSION(xdim,ydim):: TRCP_1	! 変数1(半整数値)
	REAL,DIMENSION(0:xdim,0:ydim):: TRC_2	! 変数2(格子点値)
	REAL,DIMENSION(xdim,ydim):: TRCP_2	! 変数2(半整数値)
	REAL,DIMENSION(0:xdim,0:ydim):: TRC_3	! 変数3(格子点値)
	REAL,DIMENSION(xdim,ydim):: TRCP_3	! 変数3(半整数値)
	REAL,DIMENSION(0:xdim,0:ydim):: TRC_4	! 変数4(格子点値)
	REAL,DIMENSION(xdim,ydim):: TRCP_4	! 変数4(半整数値)

	INTEGER			:: tstep	! 時間ステップ数
	REAL			:: cx		! FCT 用係数
	REAL			:: round	! 回転数
	CHARACTER(LEN=10)	:: theader	! 時刻ヘッダー
	CHARACTER(LEN=20),DIMENSION(4) :: header! スキームヘッダー
	CHARACTER(LEN=4)	:: time		! 時刻
	REAL,DIMENSION(0:xdim,0:ydim):: DTRC1	! 時間積分用配列1
	REAL,DIMENSION(0:xdim,0:ydim):: DTRC2	! 時間積分用配列2
	REAL,DIMENSION(0:xdim,0:ydim):: DTRC3	! 時間積分用配列3
	REAL,DIMENSION(0:xdim,0:ydim):: DTRC4	! 時間積分用配列4


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

	xplot = dx * (/ (i, i = 0, xdim-1) /)
	yplot = dy * (/ (i, i = 0, ydim-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_1(i,j) = amp*( 1. - SQRT((i + 0.5 - xdim0)**2 + &
       &                                (j + 0.5 - ydim0)**2)/radi ) 
	   TRC_1(i,j) = MAX( TRC_1(i,j), 0. )
	END DO 
	END DO

	CALL bound2( TRC_1 )

	TRC_2 = TRC_1
	TRC_3 = TRC_1
	TRC_4 = TRC_1

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

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

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

!- FCT 用パラメータの読み込み

	WRITE(*,*) 'INPUT CORRECTIVE FACTOR ? (R)'
	READ(*,*) cx

!- 初期値の描画

	round = 0.

	header(1) = 'U P S T R E A M'
	header(2) = 'C E N T E R E D - 2'
!	header(2) = 'I N I T I A L'
	header(3) = 'M P D A T A'
	header(4) = 'F C T'		

	WRITE(time,'(I4.4)') 0
	theader = 'tstep=' // time 

	CALL dennou_set(1)
!	CALL get_plotval( TRC_1, TRCP_1 )
!	CALL get_plotval( TRC_2, TRCP_2 )
!	CALL get_plotval( TRC_3, TRCP_3 )
!	CALL get_plotval( TRC_4, TRCP_4 )
!	CALL dennou_draw( xplot, yplot, TRCP_1, theader, header(1))
!	CALL dennou_draw( xplot, yplot, TRCP_2, theader, header(2))
!	CALL dennou_draw( xplot, yplot, TRCP_3, theader, header(3))
!	CALL dennou_draw( xplot, yplot, TRCP_4, theader, header(4))

!- 時間積分

	DO i = 1, tstep

	   WRITE(time,'(I4.4)') i
	   theader = 'tstep=' // time 

!- 1 次上流差分

!	   CALL get_dtrc_u1_2D( u, v, TRC_1, DTRC1 )
!	   CALL get_dtrc_u1_2D( u, v, TRC_1+DTRC1, DTRC2 )
!	   TRC_1 = TRC_1 + ( DTRC1 + DTRC2 )/2
!	   CALL bound2( TRC_1 )

!- 2 次中心差分

!	   CALL get_dtrc_C2_2D( u, v, TRC_2, DTRC1 )
!	   CALL get_dtrc_C2_2D( u, v, TRC_2+DTRC1, DTRC2 )
!	   TRC_2 = TRC_2 + ( DTRC1 + DTRC2 )/2
!	   CALL bound2( TRC_2 )

!- MPDATA (1 次上流差分)

!	   CALL get_dtrc_u1_2D( u, v, TRC_3, DTRC1 )
!	   TRC_3 = TRC_3 + DTRC1
!	   CALL get_dtrc_u_mpdata_2D( u, v, TRC_3, udif, vdif, DTRC1 )
!	   TRC_3 = TRC_3 + DTRC1
!	   CALL bound2( TRC_3 )

	   
!- FCT ( 1 次上流 + 2 次中心 )

	   CALL get_dtrc_u1_2D( u, v, TRC_4, DTRC1 )
	   CALL get_dtrc_u1_2D( u, v, TRC_4+DTRC1, DTRC2 )
	   CALL get_dtrc_C2_2D( u, v, TRC_4, DTRC3 )
	   CALL get_dtrc_C2_2D( u, v, TRC_4+DTRC3, DTRC4 )
	   TRC_4 = TRC_4 + ( 1. - cx )*( DTRC1 + DTRC2 )/2 + &
	&                           cx*( DTRC3 + DTRC4 )/2
	   CALL bound2( TRC_4 )

!- 描画

!	   IF ( MOD( i, 100 ) == 0 ) THEN
!	     CALL get_plotval( TRC_1, TRCP_1 )
!	     CALL get_plotval( TRC_2, TRCP_2 )
!	     CALL get_plotval( TRC_3, TRCP_3 )
!	     CALL get_plotval( TRC_4, TRCP_4 )
!	     CALL dennou_draw( xplot, yplot, TRCP_1, theader, header(1))
!	     CALL dennou_draw( xplot, yplot, TRCP_2, theader, header(2))
!	     CALL dennou_draw( xplot, yplot, TRCP_3, theader, header(3))
!	     CALL dennou_draw( xplot, yplot, TRCP_4, theader, header(4))
!	   END IF

	END DO

	CALL get_plotval( TRC_4, TRCP_4 )
	CALL dennou_draw( xplot, yplot, TRCP_4, theader, header(4))
	CALL dennou_cls

END PROGRAM main