!========================================================================= ! 一次元トレーサー移流モデル: 時間増分計算モジュール ! ! 履歴 1997/10/16 小高正嗣 ! 1997/10/22 小高正嗣 ! 1997/11/13 小高正嗣 ! 1998/01/17 小高正嗣 ! 2011/06/27 小高正嗣 ! !========================================================================= MODULE cal_module PUBLIC CONTAINS !========================================================================= SUBROUTINE get_dtrc_C2( 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_dtrc_C2 !========================================================================= SUBROUTINE get_dtrc_C4( 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_dtrc_C4 !========================================================================= SUBROUTINE get_dtrc_u1( 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_dtrc_u1 !========================================================================= SUBROUTINE get_dtrc_u_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_dtrc_u_mpdata !========================================================================= SUBROUTINE get_dtrc_fct_shasta( dt, dx, u, TRCM, AFLUX, CFLUX )! SHASTA REAL,INTENT(in) :: dt REAL,INTENT(in) :: dx REAL,DIMENSION(:),INTENT(in) :: u REAL,DIMENSION(:),INTENT(inout) :: TRCM REAL,DIMENSION(:),INTENT(inout) :: AFLUX! 反拡散フラックス REAL,DIMENSION(:),INTENT(inout) :: CFLUX! 補正フラックス INTEGER :: dim, xdim REAL,DIMENSION(:),ALLOCATABLE :: TRCMT! トレーサー変数(半整数値2) REAL,DIMENSION(:),ALLOCATABLE :: Q1 ! SHASTAパラメータ1 REAL,DIMENSION(:),ALLOCATABLE :: Q2 ! SHASTAパラメータ2 REAL,DIMENSION(:),ALLOCATABLE :: ut ! 速度(仮値) REAL,DIMENSION(:),ALLOCATABLE :: DIF ! トレーサー空間差分 REAL,DIMENSION(:),ALLOCATABLE :: SN ! AFLUX の符合 xdim = SIZE(TRCM) dim = xdim - 1 ALLOCATE( TRCMT(-1:xdim+2) ) ALLOCATE( ut(0:xdim+1) ) ALLOCATE( Q1(xdim) ) ALLOCATE( Q2(xdim) ) ALLOCATE( DIF(0:xdim+2) ) ALLOCATE( SN(xdim+1) ) !-輸送過程(transport stage) TRCMT(1:xdim) = TRCM TRCMT(-1) = TRCM(xdim-2) TRCMT(0) = TRCM(xdim-1) TRCMT(xdim+1) = TRCM(2) TRCMT(xdim+1) = TRCM(3) ut(1:xdim) = u ut(0) = u(xdim-1) ut(xdim+1) = u(2) Q1 = ( 0.5 - ut(1:xdim)*dt/dx )/( 1. + dt*( ut(2:) - ut(1:) )/dx ) Q2 = ( 0.5 + ut(1:xdim)*dt/dx )/( 1. + dt*( ut(1:) - ut(0:) )/dx ) TRCM = 0.5*Q1**2*( TRCMT(2:) - TRCMT(1:) ) & & + 0.5*Q2**2*( TRCMT(0:) - TRCMT(1:) ) & & + ( Q1 + Q2 )*TRCMT(1:) !-補正過程(corrected stage) TRCMT(1:xdim) = TRCM TRCMT(-1) = TRCM(xdim-2) TRCMT(0) = TRCM(xdim-1) TRCMT(xdim+1) = TRCM(2) TRCMT(xdim+1) = TRCM(3) DIF = TRCMT(0:) - TRCMT(:xdim+1) AFLUX = 0.125*DIF(1:) SN = SIGN(1.,AFLUX) CFLUX = SN*MAX(0., MIN( SN*DIF(:xdim), ABS(AFLUX), SN*DIF(2:) ) ) TRCM = TRCM - (CFLUX(2:) - CFLUX(:xdim)) DEALLOCATE( TRCMT, ut, Q1, Q2, DIF, SN ) END SUBROUTINE get_dtrc_fct_shasta !========================================================================= SUBROUTINE get_dtrc_nonfct_shasta( dt, dx, u, TRCM )! non-fct SHASTA REAL,INTENT(in) :: dt REAL,INTENT(in) :: dx REAL,DIMENSION(:),INTENT(in) :: u REAL,DIMENSION(:),INTENT(inout) :: TRCM INTEGER :: dim, xdim REAL,DIMENSION(:),ALLOCATABLE :: TRCMT! トレーサー変数(半整数値2) REAL,DIMENSION(:),ALLOCATABLE :: Q1 ! SHASTAパラメータ1 REAL,DIMENSION(:),ALLOCATABLE :: Q2 ! SHASTAパラメータ2 REAL,DIMENSION(:),ALLOCATABLE :: ut ! 速度(仮値) xdim = SIZE(TRCM) dim = xdim - 1 ALLOCATE( TRCMT(-1:xdim+2) ) ALLOCATE( ut(0:xdim+1) ) ALLOCATE( Q1(xdim) ) ALLOCATE( Q2(xdim) ) !-輸送過程(transport stage) TRCMT(1:xdim) = TRCM TRCMT(-1) = TRCM(xdim-2) TRCMT(0) = TRCM(xdim-1) TRCMT(xdim+1) = TRCM(2) TRCMT(xdim+1) = TRCM(3) ut(1:xdim) = u ut(0) = u(xdim-1) ut(xdim+1) = u(2) Q1 = ( 0.5 - ut(1:xdim)*dt/dx )/( 1. + dt*( ut(2:) - ut(1:) )/dx ) Q2 = ( 0.5 + ut(1:xdim)*dt/dx )/( 1. + dt*( ut(1:) - ut(0:) )/dx ) TRCM = 0.5*Q1**2*( TRCMT(2:) - TRCMT(1:) ) & & + 0.5*Q2**2*( TRCMT(0:) - TRCMT(1:) ) & & + ( Q1 + Q2 )*TRCMT(1:) DEALLOCATE( TRCMT, ut, Q1, Q2 ) END SUBROUTINE get_dtrc_nonfct_shasta !========================================================================= 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