!=========================================================================
! 二次元トレーサー移流モデル: 時間増分計算モジュール
!
! 履歴	1997/10/16	小高正嗣
!	1997/10/22	小高正嗣
!	1997/11/13	小高正嗣
!	1997/11/18	小高正嗣
!=========================================================================
MODULE cal_module
USE param_module
PUBLIC 

	CONTAINS

!=========================================================================
SUBROUTINE get_dtrc1_C2_2D( u, v, TRC, DTRC )! 二次中央差分

	REAL,DIMENSION(:,:),INTENT(in)	:: u
	REAL,DIMENSION(:,:),INTENT(in)	:: v
	REAL,DIMENSION(0:,0:),INTENT(in):: TRC
	REAL,DIMENSION(0:,0:),INTENT(out):: DTRC

	INTEGER				:: dim, xdim, ydim
	REAL,DIMENSION(:,:),ALLOCATABLE	:: FLUX_X! フラックス1
	REAL,DIMENSION(:,:),ALLOCATABLE	:: FLUX_Y! フラックス2
	REAL,DIMENSION(:,:),ALLOCATABLE	:: ADVE_X! 移流項1
	REAL,DIMENSION(:,:),ALLOCATABLE	:: ADVE_Y! 移流項2
	REAL,DIMENSION(:,:),ALLOCATABLE	:: TRCM_X! 変数(半整数値1)
	REAL,DIMENSION(:,:),ALLOCATABLE	:: TRCM_Y! 変数(半整数値2)

	xdim = SIZE( TRC, 1 ) - 1
	ydim = SIZE( TRC, 2 ) - 1
	dim  = xdim - 1

	ALLOCATE( TRCM_X(xdim,0:ydim), TRCM_Y(0:xdim,ydim) )
	ALLOCATE( FLUX_X(xdim,ydim), FLUX_Y(xdim,ydim) )
	ALLOCATE( ADVE_X(0:xdim,0:ydim), ADVE_Y(0:xdim,0:ydim) )

	CALL get_mval_2D( TRC, TRCM_X, TRCM_Y )

	FLUX_X = u * TRCM_X(:,1:) 
	FLUX_Y = v * TRCM_Y(1:,:) 

	ADVE_X(1:dim,1:dim) = - ( FLUX_X(2:,:dim) - FLUX_X(:dim,:dim) )/dx
	ADVE_Y(1:dim,1:dim) = - ( FLUX_Y(:dim,2:) - FLUX_Y(:dim,:dim) )/dy

	call bound2( ADVE_X )
	call bound2( ADVE_Y )

	DTRC = dt * ( ADVE_X + ADVE_Y )

	DEALLOCATE( TRCM_X, TRCM_Y, FLUX_X, FLUX_Y, ADVE_X, ADVE_Y )
	
END SUBROUTINE get_dtrc1_C2_2D
!=========================================================================
SUBROUTINE get_dtrc1_C4_2D( u, v, TRC, DTRC )! 四次中央差分

	REAL,DIMENSION(:,:),INTENT(in)	:: u
	REAL,DIMENSION(:,:),INTENT(in)	:: v
	REAL,DIMENSION(0:,0:),INTENT(in):: TRC
	REAL,DIMENSION(0:,0:),INTENT(out):: DTRC

	INTEGER				:: dim, xdim, ydim
	REAL,DIMENSION(:,:),ALLOCATABLE	:: FLUX_X! フラックス1
	REAL,DIMENSION(:,:),ALLOCATABLE	:: FLUX_Y! フラックス2
	REAL,DIMENSION(:,:),ALLOCATABLE	:: ADVE_X! 移流項1
	REAL,DIMENSION(:,:),ALLOCATABLE	:: ADVE_Y! 移流項2
	REAL,DIMENSION(:,:),ALLOCATABLE	:: TRCM_X! 変数(半整数値1)
	REAL,DIMENSION(:,:),ALLOCATABLE	:: TRCM_Y! 変数(半整数値2)

	xdim = SIZE( TRC, 1 ) - 1
	ydim = SIZE( TRC, 2 ) - 1
	dim  = xdim - 1

	ALLOCATE( TRCM_X(xdim,0:ydim), TRCM_Y(0:xdim,ydim) )
	ALLOCATE( FLUX_X(0:xdim+1,ydim), FLUX_Y(xdim,0:ydim+1) )
	ALLOCATE( ADVE_X(0:xdim,0:ydim), ADVE_Y(0:xdim,0:ydim) )

	CALL get_mval_2D( TRC, TRCM_X, TRCM_Y )

	FLUX_X(1:xdim,:) = u * TRCM_X(:,1:) 
	FLUX_Y(:,1:ydim) = v * TRCM_Y(1:,:) 

	FLUX_X(0,:)      = FLUX_X(2,:)
	FLUX_X(xdim+1,:) = FLUX_X(dim,:)
	FLUX_Y(:,0)      = FLUX_Y(:,2)
	FLUX_Y(:,ydim+1) = FLUX_Y(:,dim)

	ADVE_X(1:dim,1:dim) = &
	& ( - 1.125   * ( FLUX_X(2:xdim,1:dim) - FLUX_X(1:dim,1:dim)  ) &
	&   + 0.125/3.* ( FLUX_X(3:,1:dim)     - FLUX_X(:dim-1,1:dim) ) )/dx

	ADVE_Y(1:dim,1:dim) = &
	& ( - 1.125   * ( FLUX_Y(1:dim,2:ydim) - FLUX_Y(1:dim,1:dim)  ) &
	&   + 0.125/3.* ( FLUX_Y(1:dim,3:)     - FLUX_Y(1:dim,:dim-1) ) )/dy

	call bound2( ADVE_X )
	call bound2( ADVE_Y )

	DTRC = dt * ( ADVE_X + ADVE_Y )

	DEALLOCATE( TRCM_X, TRCM_Y, FLUX_X, FLUX_Y, ADVE_X, ADVE_Y )
	
END SUBROUTINE get_dtrc1_C4_2D
!=========================================================================
SUBROUTINE get_dtrc_u1_2D( u, v, TRC, DTRC )! 上流差分

	REAL,DIMENSION(:,:),INTENT(in)	:: u
	REAL,DIMENSION(:,:),INTENT(in)	:: v
	REAL,DIMENSION(0:,0:),INTENT(in):: TRC
	REAL,DIMENSION(0:,0:),INTENT(out):: DTRC

	INTEGER				:: dim, xdim, ydim
	REAL,DIMENSION(:,:),ALLOCATABLE	:: FLUX_X! フラックス1
	REAL,DIMENSION(:,:),ALLOCATABLE	:: FLUX_Y! フラックス2
	REAL,DIMENSION(:,:),ALLOCATABLE	:: ADVE_X! 移流項1
	REAL,DIMENSION(:,:),ALLOCATABLE	:: ADVE_Y! 移流項2

	xdim = SIZE( TRC, 1 ) - 1
	ydim = SIZE( TRC, 2 ) - 1
	dim  = xdim - 1

	ALLOCATE( FLUX_X(xdim,ydim), FLUX_Y(xdim,ydim) )
	ALLOCATE( ADVE_X(0:xdim,0:ydim), ADVE_Y(0:xdim,0:ydim) )
	
	FLUX_X = MAX( u, 0. )*TRC(:dim,1:) + MIN( u, 0. )*TRC(1:,1:)
	FLUX_Y = MAX( v, 0. )*TRC(1:,:dim) + MIN( v, 0. )*TRC(1:,1:)

	ADVE_X(1:dim,1:dim) = - ( FLUX_X(2:,:dim) - FLUX_X(:dim,:dim) )/dx
	ADVE_Y(1:dim,1:dim) = - ( FLUX_Y(:dim,2:) - FLUX_Y(:dim,:dim) )/dy

	call bound2( ADVE_X )
	call bound2( ADVE_Y )

	DTRC = dt * ( ADVE_X + ADVE_Y )

	DEALLOCATE( FLUX_X, FLUX_Y, ADVE_X, ADVE_Y )

END SUBROUTINE get_dtrc_u1_2D
!=========================================================================
SUBROUTINE get_dtrc_u_mpdata_2D( u, v, TRC, udif, vdif, DTRC )! MPDATA 上流差分

	REAL,DIMENSION(:,:),INTENT(in)	:: u
	REAL,DIMENSION(:,:),INTENT(in)	:: v
	REAL,DIMENSION(0:,0:),INTENT(in):: TRC
	REAL,DIMENSION(:,:),INTENT(out)	:: udif
	REAL,DIMENSION(:,:),INTENT(out)	:: vdif
	REAL,DIMENSION(0:,0:),INTENT(out):: DTRC

	INTEGER				:: dim, xdim, ydim
	REAL,DIMENSION(:,:),ALLOCATABLE	:: FLUX_X! フラックス1
	REAL,DIMENSION(:,:),ALLOCATABLE	:: FLUX_Y! フラックス2
	REAL,DIMENSION(:,:),ALLOCATABLE	:: ADVE_X! 移流項1
	REAL,DIMENSION(:,:),ALLOCATABLE	:: ADVE_Y! 移流項2
	REAL,DIMENSION(:,:),ALLOCATABLE	:: TRCM_X! 変数(半整数値1)
	REAL,DIMENSION(:,:),ALLOCATABLE	:: TRCM_Y! 変数(半整数値2)
	REAL,PARAMETER			:: epsilon = 1.e-4

	xdim = SIZE( TRC, 1 ) - 1
	ydim = SIZE( TRC, 2 ) - 1
	dim  = xdim - 1

	ALLOCATE( FLUX_X(xdim,ydim), FLUX_Y(xdim,ydim) )
	ALLOCATE( ADVE_X(0:xdim,0:ydim), ADVE_Y(0:xdim,0:ydim) )
	ALLOCATE( TRCM_X(xdim,0:ydim), TRCM_Y(0:xdim,ydim) )

	CALL get_mval_2D( TRC, TRCM_X, TRCM_Y )

	udif = 0.
	vdif = 0.

	DO i = 1, xdim
	DO j = 1, ydim
	   IF ( TRCM_X(i,j) /= 0. ) THEN
	      udif(i,j) = 0.5*ABS(u(i,j)) * ( TRC(i,j) - TRC(i-1,j) ) &
	&                       / ( TRCM_X(i,j) ) 
	   END IF
	   IF ( TRCM_Y(i,j) /= 0. ) THEN
	      vdif(i,j) = 0.5*ABS(v(i,j)) * ( TRC(i,j) - TRC(i,j-1) )  &
	&                       / ( TRCM_Y(i,j) ) 
	   END IF
	END DO
	END DO

	FLUX_X = MAX( udif, 0. )*TRC(:dim,1:) + MIN( udif, 0. )*TRC(1:,1:)
	FLUX_Y = MAX( vdif, 0. )*TRC(1:,:dim) + MIN( vdif, 0. )*TRC(1:,1:)

	ADVE_X(1:dim,1:dim) = - ( FLUX_X(2:,:dim) - FLUX_X(:dim,:dim) )/dx
	ADVE_Y(1:dim,1:dim) = - ( FLUX_Y(:dim,2:) - FLUX_Y(:dim,:dim) )/dy

	call bound2( ADVE_X )
	call bound2( ADVE_Y )

	DTRC = dt * ( ADVE_X + ADVE_Y )

	DEALLOCATE( FLUX_X, FLUX_Y, ADVE_X, ADVE_Y )

END SUBROUTINE get_dtrc_u_mpdata_2D
!=========================================================================
SUBROUTINE get_mval_2D( val, mval_x, mval_y )

	REAL,DIMENSION(0:,0:),INTENT(in):: val
	REAL,DIMENSION(:,0:),INTENT(out):: mval_x
	REAL,DIMENSION(0:,:),INTENT(out):: mval_y
	INTEGER				:: xdim,ydim

	xdim = SIZE( mval_x, 1 )
	ydim = SIZE( mval_y, 2 )

	mval_x = ( val(1:,:) + val(:xdim-1,:) )/2
	mval_y = ( val(:,1:) + val(:,:ydim-1) )/2

END SUBROUTINE get_mval_2D
!=========================================================================
SUBROUTINE get_plotval( val, pval )

	REAL,DIMENSION(0:,0:),INTENT(in):: val
	REAL,DIMENSION(:,:),INTENT(out)	:: pval
	INTEGER				:: xdim,ydim

	xdim = SIZE( pval, 1 )
	ydim = SIZE( pval, 2 )

	pval = (val(1:,1:) + val(1:,:ydim-1) + &
       &	val(:xdim-1,1:) + val(:xdim-1,:ydim-1) )/4.


END SUBROUTINE get_plotval
!=========================================================================
SUBROUTINE bound2( val )

	REAL,DIMENSION(0:,0:),INTENT(inout):: val
	INTEGER				:: dim,xdim,ydim
	
	xdim = SIZE( val, 1 ) - 1
	ydim = SIZE( val, 2 ) - 1
	dim  = xdim - 1

!-cyclic

!	val(0,:)    = val(dim,:)
!	val(xdim,:) = val(1,:)
!	val(:,0)    = val(:,dim)
!	val(:,ydim) = val(:,1)

!	val(0,0)    = val(dim,dim)
!	val(0,ydim) = val(dim,1)
!	val(xdim,0) = val(1,dim)
!	val(xdim,ydim) = val(1,1)

!-fixed boundary (symmetric)

	val(0,:)    = val(1,:)
	val(xdim,:) = val(dim,:)
	val(:,0)    = val(:,1)
	val(:,ydim) = val(:,dim)

	val(0,0)    = val(1,1)
	val(0,ydim) = val(1,dim)
	val(xdim,0) = val(dim,1)
	val(xdim,ydim) = val(dim,dim)
!	
END SUBROUTINE bound2
!=========================================================================
END MODULE cal_module