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

