!========================================================================= ! 一次元トレーサー移流モデル: 時間増分計算モジュール ! ! 履歴 1997/10/16 小高正嗣 ! 1997/10/22 小高正嗣 ! 1997/11/13 小高正嗣 !========================================================================= MODULE cal_module PUBLIC get_dtrc1,get_dtrc2,get_dtrcu,get_mval,bound CONTAINS !========================================================================= SUBROUTINE get_dtrc1( dt, dx, u, TRC, DTRC )! 二次中央差分 REAL,INTENT(in) :: dt REAL,INTENT(in) :: dx REAL,DIMENSION(:),INTENT(in) :: u REAL,DIMENSION(0:),INTENT(in) :: TRC REAL,DIMENSION(0:),INTENT(out) :: DTRC INTEGER :: dim, xdim REAL,DIMENSION(:),ALLOCATABLE :: TRCM ! トレーサー変数(半整数値) REAL,DIMENSION(:),ALLOCATABLE :: FLUX ! フラックス REAL,DIMENSION(:),ALLOCATABLE :: ADVE ! 移流項 xdim = SIZE(TRC) - 1 dim = xdim - 1 ALLOCATE( TRCM(xdim) ) ALLOCATE( FLUX(xdim) ) ALLOCATE( ADVE(0:xdim) ) CALL get_mval( TRC, TRCM ) FLUX = u * TRCM ADVE(1:dim) = - ( FLUX(2:) - FLUX(:dim) )/dx DTRC = dt * ADVE DEALLOCATE( TRCM, FLUX, ADVE ) END SUBROUTINE get_dtrc1 !========================================================================= SUBROUTINE get_dtrc2( dt, dx, u, TRC, DTRC ) ! 四次中央差分 REAL,INTENT(in) :: dt REAL,INTENT(in) :: dx REAL,DIMENSION(:),INTENT(in) :: u REAL,DIMENSION(0:),INTENT(in) :: TRC REAL,DIMENSION(0:),INTENT(out) :: DTRC INTEGER :: dim, xdim REAL,DIMENSION(:),ALLOCATABLE :: TRCM ! トレーサー変数(半整数値) REAL,DIMENSION(:),ALLOCATABLE :: GRAD ! 一階微分 REAL,DIMENSION(:),ALLOCATABLE :: FLUX ! フラックス REAL,DIMENSION(:),ALLOCATABLE :: ADVE ! 移流項 xdim = SIZE(TRC) - 1 dim = xdim - 1 ALLOCATE( TRCM(xdim) ) ALLOCATE( GRAD(0:xdim) ) ALLOCATE( FLUX(xdim) ) ALLOCATE( ADVE(0:xdim) ) CALL get_mval( TRC, TRCM ) GRAD(1:dim) = ( u(2:)*TRCM(2:) - u(:dim)*TRCM(:dim) )/dx CALL bound( GRAD ) FLUX = u*TRCM + ( GRAD(1:) - GRAD(:xdim) )*dx/24. ADVE(1:dim) = - ( FLUX(2:) - FLUX(:dim) )/dx DTRC = dt * ADVE DEALLOCATE( TRCM, FLUX, ADVE, GRAD ) END SUBROUTINE get_dtrc2 !========================================================================= SUBROUTINE get_dtrcu1( dt, dx, u, TRC, DTRC )! 上流差分 REAL,INTENT(in) :: dt REAL,INTENT(in) :: dx REAL,DIMENSION(:),INTENT(in) :: u REAL,DIMENSION(0:),INTENT(in) :: TRC REAL,DIMENSION(0:),INTENT(out) :: DTRC INTEGER :: dim, xdim REAL,DIMENSION(:),ALLOCATABLE :: FLUX ! フラックス REAL,DIMENSION(:),ALLOCATABLE :: ADVE ! 移流項 xdim = SIZE(TRC) - 1 dim = xdim - 1 ALLOCATE( FLUX(xdim) ) ALLOCATE( ADVE(0:xdim) ) FLUX = MAX( u, 0. )*TRC(:dim) + MIN( u, 0. )*TRC(1:) ADVE(1:dim) = - ( FLUX(2:) - FLUX(:dim) )/dx DTRC = dt * ADVE DEALLOCATE( FLUX, ADVE ) END SUBROUTINE get_dtrcu1 !========================================================================= SUBROUTINE get_dtrcu_mpdata( dt, dx, u, TRC, udif, DTRC )! MPDATA 上流差分 REAL,INTENT(in) :: dt REAL,INTENT(in) :: dx REAL,DIMENSION(:),INTENT(in) :: u REAL,DIMENSION(0:),INTENT(in) :: TRC REAL,DIMENSION(:),INTENT(out) :: udif REAL,DIMENSION(0:),INTENT(out) :: DTRC INTEGER :: dim, xdim REAL,DIMENSION(:),ALLOCATABLE :: TRCM ! トレーサー変数(半整数値) REAL,DIMENSION(:),ALLOCATABLE :: FLUX ! フラックス REAL,DIMENSION(:),ALLOCATABLE :: ADVE ! 移流項 REAL,PARAMETER :: epsilon = 1.e-8 xdim = SIZE(TRC) - 1 dim = xdim - 1 ALLOCATE( TRCM(xdim) ) ALLOCATE( FLUX(xdim) ) ALLOCATE( ADVE(0:xdim) ) CALL get_mval( TRC, TRCM ) udif = ( 0.5*dx*ABS(u) - dt*u**2 )*( TRC(1:) - TRC(:dim) ) / & & ( dx*( TRCM + epsilon ) ) ! udif = 0.5*ABS(u)*( TRC(1:) - TRC(:dim) ) / ( TRCM + epsilon ) FLUX = MAX( udif, 0. )*TRC(:dim) + MIN( udif, 0. )*TRC(1:) ADVE(1:dim) = - ( FLUX(2:) - FLUX(:dim) )/dx DTRC = dt * ADVE DEALLOCATE( TRCM, FLUX, ADVE ) END SUBROUTINE get_dtrcu_mpdata !========================================================================= SUBROUTINE get_mval( val, mval ) REAL,DIMENSION(0:),INTENT(in) :: val REAL,DIMENSION(:),INTENT(out) :: mval INTEGER :: xdim xdim = SIZE( mval ) mval = ( val(1:) + val(:xdim) )/2 END SUBROUTINE get_mval !========================================================================= SUBROUTINE bound( val ) REAL,DIMENSION(0:),INTENT(inout):: val INTEGER :: dim,xdim xdim = SIZE( val ) - 1 dim = xdim - 1 val(0) = val(dim) val(xdim) = val(1) END SUBROUTINE bound !========================================================================= END MODULE cal_module