!= Module NumDiffusion_3D
!
! Authors::   SUGIYAMA Ko-ichiro, ODAKA Masatsugu 
! Version::   $Id: numdiffusion_3d.f90,v 1.2 2007/10/01 09:01:13 odakker Exp $
! Tag Name::  $Name: arare4-20071012 $
! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview
!
!数値拡散項の計算モジュール
! 
!== Error Handling
!
!== Bugs
!
!== Note
!
!== Future Plans
!
!

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

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

  use gridset_3d, only:  DimXMin,           &! x 方向の配列の下限
    &                 DimXMax,           &! x 方向の配列の上限
    &                 DimYMin,           &! y 方向の配列の下限
    &                 DimYMax,           &! y 方向の配列の上限
    &                 DimZMin,           &! z 方向の配列の下限
    &                 DimZMax,           &! z 方向の配列の上限
    &                 SpcNum,            &!
    &                 x_dx,              &! X 方向の格子点間隔
    &                 y_dy,              &! X 方向の格子点間隔
    &                 z_dz                ! Z 方向の格子点間隔
  use timeset, only:  DelTimeLong         !長い時間ステップ
  use xyz_deriv_module, only: xyz_dx_pyz, xyz_dy_xqz, xyz_dz_xyr, &
    &                         pyz_dx_xyz, xqz_dy_xyz, xyr_dz_xyz, &
    &                         xyr_dx_pyr, pyz_dy_pqz, pyz_dz_pyr, & 
    &                         pyr_dx_xyr, pqz_dy_pyz, pyr_dz_pyz, &
    &                         xqz_dx_pqz, xyr_dy_xqr, xqz_dz_xqr, &
    &                         pqz_dx_xqz, xqr_dy_xyr, xqr_dz_xqz
  use StorePotTemp_3d, only: StorePotTempDiff


  !暗黙の型宣言禁止
  implicit none

  !属性の指定
  private

  !関数を public にする
  public NumDiffusion_Init
  public xyz_NumDiffScalar
!  public xyza_NumDiffScalar
  public pyz_NumDiffVelX
  public xqz_NumDiffVelY
  public xyr_NumDiffVelZ

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

contains

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

    !暗黙の型宣言禁止
    implicit none

    real(DP) :: DelXMin, DelYMin, DelZMin

    DelXMin = minval(x_dx)
    DelYMin = minval(y_dy)
    DelZMin = minval(z_dz)
    
    ! CReSS マニュアルでは, 2 次精度中心差分の場合, Alpha < 1/8 くらいが適当と述べている. 
    ! deepconv では NuH として 500 以上の値が入っていたそうだ,
    NuH  = Alpha * ( SQRT(DelXMin*DelYMin) ** 2.0d0 ) / DelTimeLong
    NuV  = Alpha * ( DelZMin ** 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 xyz_NumDiffScalar(xyz_Scalar)
    !
    ! x, y, z 方向に半格子ずれた点での数値拡散項を評価
    !
    
    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(DP), intent(in) :: xyz_Scalar &
      &                     (DimXMin:DimXMax,DimYMin:DImYMax,DimZMin:DimZMax)
                                                    !スカラー量
    real(DP)             :: xyz_NumDiffScalar &
      &                     (DimXMin:DimXMax,DimYMin:DImYMax,DimZMin:DimZMax)
                                                    !水平方向の数値拡散
    
!    xz_NumDiffScalar = 0.0d0  !初期化
    xyz_NumDiffScalar =   &
      &    NuH * (xyz_dx_pyz(pyz_dx_xyz( xyz_Scalar ))) &
      &  + NuH * (xyz_dy_xqz(xqz_dy_xyz( xyz_Scalar ))) &
      &  + NuV * (xyz_dz_xyr(xyr_dz_xyz( xyz_Scalar ))) 
    
    call StorePotTempDiff( xyz_NumDiffScalar )

  end function xyz_NumDiffScalar

!!!-------------------------------------------------------------------!!!
!  function xyza_NumDiffScalar(xyza_Scalar)
!    !
!    ! x, z 方向に半格子ずれた点での数値拡散項を評価
!    !
!    
!    !暗黙の型宣言禁止
!    implicit none!
!
!    !変数定義
!    real(DP), intent(in) :: xyza_Scalar &
!      &                     (DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax,SpcNum)
                                                    !スカラー量
!    real(DP)             :: xyza_NumDiffScalar &
!      &                     (DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax,SpcNum)
                                                    !水平方向の数値拡散
!    integer             :: s
    

!    do s = 1, SpcNum
!      xyza_NumDiffScalar(:,:,:,s) =   &
!        &    NuH * (xyz_dx_pyz(pyz_dx_xyz( xyza_Scalar(:,:,:,s) ))) &
!        &  + NuH * (xyz_dy_xqz(xqz_dy_xyz( xyza_Scalar(:,:,:,s) ))) &
!        &  + NuV * (xyz_dz_xyr(xyr_dz_xyz( xyza_Scalar(:,:,:,s) ))) 
!    end do

!    call Store2Diff( xyza_NumDiffScalar )

!  end function xyza_NumDiffScalar


!!!-------------------------------------------------------------------!!!
  function pyz_NumDiffVelX(pyz_VarX)
    !
    ! z 方向に半格子ずれた点での数値拡散項を評価
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP), intent(in) :: pyz_VarX &
      &                     (DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                                    !物理量
    real(DP)             :: pyz_NumDiffVelX &
      &                     (DimXMin:DimXMax,DimYMin:DImYMax,DimZMin:DimZMax)
                                                    !数値拡散
    
    pyz_NumDiffVelX = 0.0d0  !初期化
    pyz_NumDiffVelX = &
      &   NuH * ( pyz_dx_xyz( xyz_dx_pyz( pyz_VarX ) ) ) &
      & + NuH * ( pyz_dy_pqz( pqz_dy_pyz( pyz_VarX ) ) ) &
      & + NuV * ( pyz_dz_pyr( pyr_dz_pyz( pyz_VarX ) ) )
    
  end function pyz_NumDiffVelX
  

!!!-------------------------------------------------------------------!!!
  function xqz_NumDiffVelY(xqz_VarY)
    !
    ! z 方向に半格子ずれた点での数値拡散項を評価
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP), intent(in) :: xqz_VarY &
      &                     (DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                                    !物理量
    real(DP)             :: xqz_NumDiffVelY &
      &                     (DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                                    !数値拡散
    
!   xqz_NumDiffVelY = 0.0d0  !初期化
    xqz_NumDiffVelY = &
      &   NuH * ( xqz_dx_pqz( pqz_dx_xqz( xqz_VarY ) ) ) &
      & + NuH * ( xqz_dy_xyz( xyz_dy_xqz( xqz_VarY ) ) ) &
      & + NuV * ( xqz_dz_xqr( xqr_dz_xqz( xqz_VarY ) ) )
    
  end function xqz_NumDiffVelY


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

!   xyr_NumDiffVelZ = 0.0d0 
    xyr_NumDiffVelZ = &
      &   NuH * ( xyr_dx_pyr( pyr_dx_xyr( xyr_VarZ ) ) ) &
      & + NuH * ( xyr_dy_xqr( xqr_dy_xyr( xyr_VarZ ) ) ) & 
      & + NuV * ( xyr_dz_xyz( xyz_dz_xyr( xyr_VarZ ) ) )
    
  end function xyr_NumDiffVelZ

end module NumDiffusion_3d