!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2005. All rights reserved.
!---------------------------------------------------------------------
!= Module NumDiffusion
!
!   * Developer: SUGIYAMA Ko-ichiro
!   * Version: $Id: $
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview
!
!数値拡散項の計算モジュール
! 
!== Error Handling
!
!== Bugs
!
!== Note
!
!== Future Plans
!
!

module NumDiffusion
  !
  !数値拡散項の計算モジュール
  !現在は 2 次精度中心差分を利用
  !

  !モジュール読み込み
  use gridset, only:  DimXMin,           &! x 方向の配列の下限
    &                 DimXMax,           &! x 方向の配列の上限
    &                 DimZMin,           &! z 方向の配列の下限
    &                 DimZMax,           &! z 方向の配列の上限
    &                 DelX,              &! X 方向の格子点間隔
    &                 DelZ                ! Z 方向の格子点間隔
  use timeset, only:  DelTimeLong         !長い時間ステップ
  use differentiate_center2, only: xz_dx_pz, xz_dz_xr, &
    &                              xr_dx_pr, xr_dz_xz, &
    &                              pz_dx_xz, pz_dz_pr, &
    &                              pr_dx_xr, pr_dz_pz

  !暗黙の型宣言禁止
  implicit none

  !属性の指定
  private

  !関数を public にする
  public NumDiffusion_Init
  public xz_NumDiffScalar
  public pz_NumDiffVelX
  public xr_NumDiffVelZ

  !変数の定義
  real(8)  :: NuH   = 0.0d0   !数値粘性の係数 (水平方向)
  real(8)  :: NuV   = 0.0d0   !数値粘性の係数 (鉛直方向)
  real(8)  :: Alpha = 1.0d-5
 
  save NuH, NuV, Alpha

contains

!!!-------------------------------------------------------------------!!!
  subroutine NumDiffusion_init()
    !
    ! NumDiffusion モジュールの初期化ルーチン
    !

    !暗黙の型宣言禁止
    implicit none
    
    ! 2 次精度中心差分の場合
    ! CReSS マニュアルでは, Alpha < 1/8 くらいが適当と述べている. 
    NuH  = Alpha * ( DelX ** 2.0d0 ) / DelTimeLong
    NuV  = Alpha * ( DelZ ** 2.0d0 ) / DelTimeLong
    
    !  ! 4 次精度中心差分の場合は以下を利用. 
    !  ! CReSS マニュアルでは, Alpha = 0.001 くらいが適当と述べている. 
    !  NuH  = Alpha * ( DelX ** 4.0d0 ) / DelTimeLong
    !  NuV  = Alpha * ( DelZ ** 4.0d0 ) / DelTimeLong

    !確認
    write(*,*) "NumDiffusion_init, NuH: ", NuH
    write(*,*) "NumDiffusion_init, NuV: ", NuV

  end subroutine NumDiffusion_init

!!!-------------------------------------------------------------------!!!
  function xz_NumDiffScalar(xz_Scalar)
    !
    ! x, z 方向に半格子ずれた点での数値拡散項を評価
    !
    
    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(8), intent(in) :: xz_Scalar(DimXMin:DimXMax, DimZMin:DimZMax)
                                                    !スカラー量
    real(8)             :: xz_NumDiffScalar(DimXMin:DimXMax, DimZMin:DimZMax)
                                                    !水平方向の数値拡散
    
    xz_NumDiffScalar = 0.0d0  !初期化
    xz_NumDiffScalar =   &
      &    NuH * (xz_dx_pz(pz_dx_xz( xz_Scalar ))) &
      &  + NuV * (xz_dz_xr(xr_dz_xz( xz_Scalar ))) 
    
  end function xz_NumDiffScalar


!!!-------------------------------------------------------------------!!!
  function pz_NumDiffVelX(pz_VarX)
    !
    ! z 方向に半格子ずれた点での数値拡散項を評価
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(8), intent(in) :: pz_VarX(DimXMin:DimXMax, DimZMin:DimZMax)
                                                    !物理量
    real(8)             :: pz_NumDiffVelX(DimXMin:DimXMax, DimZMin:DimZMax)
                                                    !数値拡散
    
    pz_NumDiffVelX = 0.0d0  !初期化
    pz_NumDiffVelX = &
      &   NuH * ( pz_dx_xz( xz_dx_pz( pz_VarX ) ) ) &
      & + NuV * ( pz_dz_pr( pr_dz_pz( pz_VarX ) ) )
    
  end function pz_NumDiffVelX
  

!!!-------------------------------------------------------------------!!!  
  function xr_NumDiffVelZ(xr_VarZ)
    !
    ! x 方向に半格子ずれた点での数値拡散項を評価
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(8), intent(in) :: xr_VarZ(DimXMin:DimXMax, DimZMin:DimZMax)
                                                    !物理量
    real(8)             :: xr_NumDiffVelZ(DimXMin:DimXMax, DimZMin:DimZMax)
                                                    !数値拡散

    xr_NumDiffVelZ = 0.0d0 
    xr_NumDiffVelZ = &
      &   NuH * ( xr_dx_pr( pr_dx_xr( xr_VarZ ) ) ) &
      & + NuV * ( xr_dz_xz( xz_dz_xr( xr_VarZ ) ) )
    
  end function xr_NumDiffVelZ

end module NumDiffusion