!=========================================================================
! 一次元トレーサー移流モデル: 時間増分計算モジュール
!
! 履歴	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
