!= Module NumDiffusion
!
! Authors::   SUGIYAMA Ko-ichiro
! Version::   $Id: numdiffusion.f90,v 1.13 2007-08-24 07:44:16 sugiyama Exp $
! Tag Name::  $Name: arare4-20080627 $
! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview
!
!数値拡散項の計算モジュール
! 
!== Error Handling
!
!== Bugs
!
!== Note
!
!== Future Plans
!
!

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

  !モジュール読み込み
  use dc_message, only: MessageNotify

  use gridset, only:  DimXMin,           &! x 方向の配列の下限
    &                 DimXMax,           &! x 方向の配列の上限
    &                 DimZMin,           &! z 方向の配列の下限
    &                 DimZMax,           &! z 方向の配列の上限
    &                 SpcNum,            &!
    &                 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
  use StorePotTemp, only: StorePotTempDiff
  use StoreMixRt,   only: StoreMixRtDiff

  !暗黙の型宣言禁止
  implicit none

  !属性の指定
  private

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

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

contains

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

    !暗黙の型宣言禁止
    implicit none
    
    ! CReSS マニュアルでは, 2 次精度中心差分の場合, Alpha < 1/8 くらいが適当と述べている. 
    ! deepconv では NuH として 500 以上の値が入っていたそうだ,
    NuH  = Alpha * ( DelX ** 2.0d0 ) / DelTimeLong
    NuV  = Alpha * ( DelZ ** 2.0d0 ) / DelTimeLong
    
    !確認
    call MessageNotify( "M", "NumDiffusion_init", "NuH = %f", d=(/NuH/) )
    call MessageNotify( "M", "NumDiffusion_init", "NuV = %f", d=(/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 ))) 
    
    call StorePotTempDiff( xz_NumDiffScalar )

  end function xz_NumDiffScalar


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

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


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

    !変数定義
    real(8), intent(in) :: xza_Scalar(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                                    !スカラー量
    real(8)             :: xza_NumDiffScalar(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                                    !水平方向の数値拡散
    integer             :: s
    

    do s = 1, SpcNum
      xza_NumDiffScalar(:,:,s) =   &
        &    NuH * (xz_dx_pz(pz_dx_xz( xza_Scalar(:,:,s) ))) &
        &  + NuV * (xz_dz_xr(xr_dz_xz( xza_Scalar(:,:,s) ))) 
    end do

    call StoreMixRtDiff( xza_NumDiffScalar )

  end function xza_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