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

	CONTAINS

!=========================================================================
SUBROUTINE get_dtrc_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! 変数(u 格子上)
	REAL,DIMENSION(:,:),ALLOCATABLE	:: TRCM_Y! 変数(v 格子上)

	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 )

!- 境界上での人工減衰

	DTRC(1,:)   = - 5.*dt*TRC(1,:)
	DTRC(dim,:) = - 5.*dt*TRC(dim,:)
	DTRC(:,1)   = - 5.*dt*TRC(:,1)
	DTRC(:,dim) = - 5.*dt*TRC(:,dim)

	DEALLOCATE( TRCM_X, TRCM_Y, FLUX_X, FLUX_Y, ADVE_X, ADVE_Y )
	
END SUBROUTINE get_dtrc_C2_2D
!=========================================================================
SUBROUTINE get_dtrc_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! 変数(u 格子上)
	REAL,DIMENSION(:,:),ALLOCATABLE	:: TRCM_Y! 変数(v 格子上)

	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 )

!- 境界上での人工減衰

	DTRC(1,:)   = - 5.*dt*TRC(1,:)
	DTRC(dim,:) = - 5.*dt*TRC(dim,:)
	DTRC(:,1)   = - 5.*dt*TRC(:,1)
	DTRC(:,dim) = - 5.*dt*TRC(:,dim)

	DEALLOCATE( TRCM_X, TRCM_Y, FLUX_X, FLUX_Y, ADVE_X, ADVE_Y )
	
END SUBROUTINE get_dtrc_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	:: TRCM_X! 変数(u 格子上)
	REAL,DIMENSION(:,:),ALLOCATABLE	:: TRCM_Y! 変数(v 格子上)
	REAL,DIMENSION(:,:),ALLOCATABLE	:: DTRC_X! 変数X方向差分(u 格子上)
	REAL,DIMENSION(:,:),ALLOCATABLE	:: DTRC_Y! 変数Y方向差分(v 格子上)
	REAL,DIMENSION(:,:),ALLOCATABLE	:: DTRC_Xvar! 変数X方向差分(v 格子上)
	REAL,DIMENSION(:,:),ALLOCATABLE	:: DTRC_Yvar! 変数Y方向差分(u 格子上)
	REAL,DIMENSION(:,:),ALLOCATABLE	:: u_var ! u(v格子上)
	REAL,DIMENSION(:,:),ALLOCATABLE	:: v_var ! v(u格子上)
	REAL,PARAMETER			:: epsilon = 1.e-8

	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( DTRC_X(xdim,ydim), DTRC_Y(xdim,ydim) )
	ALLOCATE( DTRC_Xvar(xdim,ydim), DTRC_Yvar(xdim,ydim) )
	ALLOCATE( u_var(xdim,ydim), v_var(xdim,ydim) )

	CALL get_mval_2D( TRC, TRCM_X, TRCM_Y )

	DTRC_X = TRC(1:,1:) - TRC(:dim,1:)
	DTRC_Y = TRC(1:,1:) - TRC(1:,:dim)

	CALL get_mval_utov( DTRC_X, DTRC_Xvar )
	CALL get_mval_vtou( DTRC_Y, DTRC_Yvar )
	CALL get_mval_utov( u, u_var )
	CALL get_mval_vtou( v, v_var )

	udif = 0.5*( ( ABS(u)*dx - dt*u**2 ) * DTRC_X / dx &
	&           - dt*u*v_var*DTRC_Yvar / dy ) / ( TRCM_X(:,1:) + epsilon )

	vdif = 0.5*( ( ABS(v)*dx - dt*v**2 ) * DTRC_Y / dy &
	&           - dt*u_var*v*DTRC_Xvar / dx ) / ( TRCM_Y(1:,:) + epsilon )

	CALL get_dtrc_u1_2D( udif, vdif, TRC, DTRC )	! 上流差分	

	DEALLOCATE( TRCM_X, TRCM_Y, DTRC_X, DTRC_Y, &
	&           DTRC_Xvar, DTRC_Yvar, u_var, v_var )

END SUBROUTINE get_dtrc_u_mpdata_2D
!=========================================================================
SUBROUTINE get_dtrc_fct_2D( u, v, TRC, DTRC )! FCT

	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	:: AFLUX_X! antidifuusiveフラックス1
	REAL,DIMENSION(:,:),ALLOCATABLE	:: AFLUX_Y! antidifuusiveフラックス2
	REAL,DIMENSION(:,:),ALLOCATABLE	:: ADVE_X! 移流項1
	REAL,DIMENSION(:,:),ALLOCATABLE	:: ADVE_Y! 移流項2
	REAL,DIMENSION(:,:),ALLOCATABLE :: TRC_L! 低次スキームによる計算値
	REAL,DIMENSION(:,:),ALLOCATABLE :: TRC_MAX
	REAL,DIMENSION(:,:),ALLOCATABLE :: TRC_MIN
	REAL,DIMENSION(:,:),ALLOCATABLE :: CF_X	! 補正係数1
	REAL,DIMENSION(:,:),ALLOCATABLE :: CF_Y	! 補正係数2

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

	ALLOCATE( ADVE_X(0:xdim,0:ydim), ADVE_Y(0:xdim,0:ydim) )
	ALLOCATE( AFLUX_X(xdim,ydim), AFLUX_Y(xdim,ydim) )
	ALLOCATE( TRC_L(0:xdim,0:xdim) )
	ALLOCATE( TRC_MAX(0:xdim,0:xdim), TRC_MIN(0:xdim,0:xdim) )
	ALLOCATE( CF_X(xdim,xdim), CF_Y(xdim,xdim) )

!- 低次スキームフラックスの計算: 上流差分

	CALL get_dtrc_u1_2D( u, v, TRC, TRC_L )
	TRC_L = TRC + TRC_L

!- antidiffusive フラックスの計算

	call get_aflux( TRC, AFLUX_X, AFLUX_Y )

!- リミッタの設定(1)

	call limit_aflux( TRC_L, AFLUX_X, AFLUX_Y )

!- リミット値の設定

	call get_limitval( TRC, TRC_L, TRC_MAX, TRC_MIN )

!- 補正係数の計算
	
	call get_cfact( TRC_L, TRC_MAX, TRC_MIN, AFLUX_X, AFLUX_Y, CF_X, CF_Y )

!- 補正フラックスの計算

	AFLUX_X = CF_X * AFLUX_X 
	AFLUX_Y = CF_Y * AFLUX_Y 

	ADVE_X(1:dim,1:dim) = - ( AFLUX_X(2:,:dim) - AFLUX_X(:dim,:dim) )/dx
	ADVE_Y(1:dim,1:dim) = - ( AFLUX_Y(:dim,2:) - AFLUX_Y(:dim,:dim) )/dy

	call bound2( ADVE_X )
	call bound2( ADVE_Y )

!- 増分計算

	DTRC = TRC_L - TRC + dt * ( ADVE_X + ADVE_Y )

	DEALLOCATE( ADVE_X, ADVE_Y, AFLUX_X, AFLUX_Y, TRC_L )
	DEALLOCATE( TRC_MAX, TRC_MIN, CF_X, CF_Y )

	CONTAINS
!-------------------------------------------------------------------------
	SUBROUTINE get_aflux( TRC, AFLUX_X, AFLUX_Y )

	REAL,DIMENSION(0:,0:),INTENT(in):: TRC	
	REAL,DIMENSION(:,:),INTENT(out) :: AFLUX_X, AFLUX_Y

	REAL,DIMENSION(:,:),ALLOCATABLE	:: FLUX_XH! フラックス1(高次スキーム)
	REAL,DIMENSION(:,:),ALLOCATABLE	:: FLUX_YH! フラックス2(高次スキーム)
	REAL,DIMENSION(:,:),ALLOCATABLE	:: FLUX_XL! フラックス1(低次スキーム)
	REAL,DIMENSION(:,:),ALLOCATABLE	:: FLUX_YL! フラックス2(低次スキーム)
	REAL,DIMENSION(:,:),ALLOCATABLE	:: TRCM_X! 変数(u 格子上)
	REAL,DIMENSION(:,:),ALLOCATABLE	:: TRCM_Y! 変数(v 格子上)

	ALLOCATE( FLUX_XL(xdim,ydim), FLUX_YL(xdim,ydim) )
	ALLOCATE( FLUX_XH(xdim,ydim), FLUX_YH(xdim,ydim) )
	ALLOCATE( TRCM_X(xdim,0:ydim), TRCM_Y(0:xdim,ydim) )

!- 低次スキームフラックスの計算: 上流差分

	FLUX_XL = MAX( u, 0. )*TRC(:dim,1:) + MIN( u, 0. )*TRC(1:,1:)
	FLUX_YL = MAX( v, 0. )*TRC(1:,:dim) + MIN( v, 0. )*TRC(1:,1:)

!- 高次スキームフラックスの計算: 二次中央差分

	CALL get_mval_2D( TRC, TRCM_X, TRCM_Y )
	FLUX_XH = u * TRCM_X(:,1:) 
	FLUX_YH = v * TRCM_Y(1:,:) 

!- antidiffusive フラックスの計算

	AFLUX_X = FLUX_XH - FLUX_XL
	AFLUX_Y = FLUX_YH - FLUX_YL

	DEALLOCATE( FLUX_XL, FLUX_YL, FLUX_XH, FLUX_YH )

	END SUBROUTINE get_aflux
!-------------------------------------------------------------------------
	SUBROUTINE limit_aflux( TRC_L, AFLUX_X, AFLUX_Y )

	REAL,DIMENSION(0:,0:),INTENT(in):: TRC_L	
	REAL,DIMENSION(:,:),INTENT(inout) :: AFLUX_X, AFLUX_Y

	INTEGER :: i,j
	REAL,DIMENSION(:,:),ALLOCATABLE	:: DIFF_X! 差分1
	REAL,DIMENSION(:,:),ALLOCATABLE	:: DIFF_Y! 差分2


	ALLOCATE( DIFF_X(0:xdim+1,ydim), DIFF_Y(xdim,0:ydim+1) )

	DIFF_X = 0.
	DIFF_Y = 0.
	DIFF_X(1:xdim,:) = TRC_L(1:,1:) - TRC_L(:dim,1:)
	DIFF_Y(:,1:ydim) = TRC_L(1:,1:) - TRC_L(1:,:dim)

	DO i = 1,xdim
	DO j = 1,ydim
	   IF ( AFLUX_X(i,j)*DIFF_X(i,j) < 0. .AND. &
	&       AFLUX_X(i,j)*DIFF_X(i-1,j) < 0. .OR. &
	&       AFLUX_X(i,j)*DIFF_X(i+1,j) < 0. ) AFLUX_X(i,j) = 0.

	   IF ( AFLUX_Y(i,j)*DIFF_Y(i,j) < 0. .AND. &
	&       AFLUX_Y(i,j)*DIFF_Y(i,j-1) < 0. .OR. &
	&       AFLUX_Y(i,j)*DIFF_Y(i,j+1) < 0. ) AFLUX_Y(i,j) = 0.
	END DO
	END DO

	DEALLOCATE( DIFF_X, DIFF_Y )

	END SUBROUTINE limit_aflux
!-------------------------------------------------------------------------
	SUBROUTINE get_limitval( TRC, TRC_L, TRC_MAX, TRC_MIN )

	REAL,DIMENSION(0:,0:),INTENT(in):: TRC, TRC_L	
	REAL,DIMENSION(0:,0:),INTENT(out) :: TRC_MAX, TRC_MIN

	REAL,DIMENSION(:,:),ALLOCATABLE :: TRC_A
	REAL,DIMENSION(:,:),ALLOCATABLE :: TRC_B

	ALLOCATE( TRC_A(0:xdim,0:xdim), TRC_B(0:xdim,0:xdim) )

	TRC_A = MAX( TRC, TRC_L )
	TRC_B = MIN( TRC, TRC_L )

	TRC_MAX(1:dim,1:dim) = MAX( TRC_A(1:dim,  1:dim),   &
	&                           TRC_A(0:dim-1,1:dim),   & 
	&                           TRC_A(2:xdim, 1:dim),   & 
	&                           TRC_A(1:dim,  0:dim-1), & 
	&                           TRC_A(1:dim,  2:ydim ) )

	TRC_MIN(1:dim,1:dim) = MIN( TRC_B(1:dim,  1:dim),   &
	&                           TRC_B(0:dim-1,1:dim),   & 
	&                           TRC_B(2:xdim, 1:dim),   & 
	&                           TRC_B(1:dim,  0:dim-1), & 
	&                           TRC_B(1:dim,  2:ydim ) )

	call bound2( TRC_MAX )
	call bound2( TRC_MIN )

	DEALLOCATE( TRC_A, TRC_B )

	END SUBROUTINE get_limitval
!-------------------------------------------------------------------------
	SUBROUTINE get_cfact( TRC_L, TRC_MAX, TRC_MIN, AFLUX_X, AFLUX_Y, CF_X, CF_Y )

	REAL,DIMENSION(0:,0:),INTENT(in):: TRC_L
	REAL,DIMENSION(0:,0:),INTENT(in) :: TRC_MAX, TRC_MIN
	REAL,DIMENSION(:,:),INTENT(in) :: AFLUX_X, AFLUX_Y
	REAL,DIMENSION(:,:),INTENT(out):: CF_X, CF_Y

	REAL,DIMENSION(:,:),ALLOCATABLE	:: SN
	REAL,DIMENSION(:,:),ALLOCATABLE	:: PP, QP, RP
	REAL,DIMENSION(:,:),ALLOCATABLE	:: PN, QN, RN

	ALLOCATE( SN(xdim,ydim) )
	ALLOCATE( PP(0:xdim,0:ydim), QP(0:xdim,0:ydim), RP(0:xdim,0:ydim) )
	ALLOCATE( PN(0:xdim,0:ydim), QN(0:xdim,0:ydim), RN(0:xdim,0:ydim) )

	PP(1:dim,1:dim) = MAX(0., AFLUX_X(:dim,:))*dy &
	&               - MIN(0., AFLUX_X(2:,:))*dy   &
	&               + MAX(0., AFLUX_Y(:,:dim))*dx &
	&               - MIN(0., AFLUX_Y(:,2:))*dx 

	PP(0,:)    = - MIN( 0., AFLUX_X(1,:) )*dy 
	PP(xdim,:) = MAX( 0., AFLUX_X(xdim,:) )*dy
	PP(:,0)    = - MIN( 0., AFLUX_Y(:,1) )*dx
	PP(:,ydim) = MAX( 0., AFLUX_X(:,ydim) )*dx

	QP = ( TRC_MAX - TRC_L )*dx*dy

	RP = 0. 
	RP = MIN( 1., QP/(PP + 1.e-8))

	PN(1:dim,1:dim) = MAX(0., AFLUX_X(2:,:))*dy &
	&               - MIN(0., AFLUX_X(:dim,:))*dy &
	&               + MAX(0., AFLUX_Y(:,2:))*dx &
	&               - MIN(0., AFLUX_Y(:,:dim))*dx

	PN(0,:)    = MAX( 0., AFLUX_X(1,:) )*dy 
	PN(xdim,:) = - MIN( 0., AFLUX_X(xdim,:) )*dy 
	PN(:,0)    = MAX( 0., AFLUX_Y(:,1) )*dx 
	PN(:,ydim) = - MIN( 0., AFLUX_X(:,ydim) )*dx 

	QN = ( TRC_L - TRC_MIN )*dx*dy

	RN = 0. 
	RN = MIN( 1., QP/(PN + 1.e-8))

	SN = INT( SIGN( 1., AFLUX_X ) )
	CF_X = MIN( MAX( 0.5*( 1 + SN )*RP(1:,1:), 0.5*( 1 - SN )*RN(1:,1:) ) &
	&     , MAX( 0.5*( 1 + SN )*RN(:dim,1:), 0.5*( 1 - SN )*RP(:dim,1:) ) )

	SN = INT( SIGN( 1., AFLUX_Y ) )
	CF_Y = MIN( MAX( 0.5*( 1 + SN )*RP(1:,1:), 0.5*( 1 - SN )*RN(1:,1:) ) &
	&     , MAX( 0.5*( 1 + SN )*RN(1:,:dim), 0.5*( 1 - SN )*RP(1:,:dim) ) )

	DEALLOCATE( SN, PP, QP, RP, PN, QN, RN )

	END SUBROUTINE get_cfact
!-------------------------------------------------------------------------
END SUBROUTINE get_dtrc_fct_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,dim

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

!- 2 点平均
!	mval_x = ( val(1:,:) + val(:xdim-1,:) )/2
!	mval_y = ( val(:,1:) + val(:,:ydim-1) )/2

!- 6 点平均
	mval_x(:,1:dim) = ( val(1:,0:dim-1)   + val(:dim,0:dim-1)   + &
	&                   2.0*val(1:,1:dim) + 2.0*val(:dim,1:dim) + &
	&                   val(1:,2:ydim)    + val(:dim,2:ydim)    )/8. 

	mval_y(1:dim,:) = ( val(0:dim-1,1:)   + val(0:dim-1,:dim)   + &
	&                   2.0*val(1:dim,1:) + 2.0*val(1:dim,:dim) + &
	&                   val(2:xdim,1:)    + val(2:xdim,:dim)    )/8.

END SUBROUTINE get_mval_2D
!=========================================================================
SUBROUTINE get_mval_utov( val_u, val_v )

	REAL,DIMENSION(:,:),INTENT(in)	:: val_u	! u 格子上変数値
	REAL,DIMENSION(:,:),INTENT(out)	:: val_v	! v 格子上変数値
	INTEGER				:: xdim,ydim,dim

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

	val_v(:dim,2:) = ( val_u(:dim,2:)   + val_u(2:xdim,2:) + &
	&                  val_u(:dim,1:dim)+ val_u(2:xdim,1:dim) )/4.

	val_v(:dim,1)  = ( val_u(:dim,1)    + val_u(2:xdim,1) )/2.
	val_v(xdim,2:) = ( val_u(xdim,:dim) + val_u(xdim,2:)  )/2.
	val_v(xdim,1)  =   val_u(xdim,1)

END SUBROUTINE get_mval_utov
!=========================================================================
SUBROUTINE get_mval_vtou( val_v, val_u )

	REAL,DIMENSION(:,:),INTENT(in)	:: val_v	! v 格子上変数値
	REAL,DIMENSION(:,:),INTENT(out)	:: val_u	! u 格子上変数値
	INTEGER				:: xdim,ydim,dim

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

	val_u(2:,:dim) = ( val_v(2:,:dim)   + val_v(2:,2:ydim) + &
	&                  val_v(1:dim,:dim)+ val_v(1:dim,2:ydim) )/4.

	val_u(1,:dim)  = ( val_v(1,:dim)    + val_v(1,2:ydim) )/2.
	val_u(2:,ydim) = ( val_v(:dim,ydim) + val_v(1:,ydim)  )/2.
	val_u(1,ydim)  =   val_v(1,ydim)

END SUBROUTINE get_mval_vtou
!=========================================================================
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 )

!- 4 点平均
	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

!- 周期境界条件

!	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)

!- 壁で対称(1階微分が0)

	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






