!========================================================================= ! 一次元トレーサー移流モデル ! ! 履歴 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