!= Module Damping
!
! Authors::   SUGIYAMA Ko-ichiro, ODAKA Masatsugu
! Version::   $Id: damping.f90,v 1.5 2011-06-17 19:07:27 sugiyama Exp $
! Tag Name::  $Name: arare5-20110623-2 $
! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview
!
!減衰率とその計算を行うためのパッケージ型モジュール
!  * 音波減衰項の係数
!  * スポンジ層の設定(境界付近で波の反射を抑え吸収するための層)
! 
!== Error Handling
!
!== Bugs
!
!== Note
!
!  * この関数は, 基本場が零な変数(速度, エクスナー関数の擾乱)に
!    適用することを想定している. 
!  * 各格子点に対する関数を定義する必要がある
!
!== Future Plans
!
!

module Damping
  !
  !減衰率とその計算を行うためのパッケージ型モジュール
  !  * 音波減衰項の係数
  !  * スポンジ層の設定(境界付近で波の反射を抑え吸収するための層)
  ! 

  !モジュール読み込み
  use dc_types,   only : DP
  use dc_iounit,  only: FileOpen
  use dc_message, only: MessageNotify
  use timeset, only: DelTimeLong
  use gridset, only: imin,       &! x 方向の配列の下限
    &                imax,       &! x 方向の配列の上限
    &                jmin,       &! y 方向の配列の下限
    &                jmax,       &! y 方向の配列の上限
    &                kmin,       &! z 方向の配列の下限
    &                kmax,       &! z 方向の配列の上限
    &                nx,         &! x 方向の物理領域の上限
    &                ny,         &! y 方向の物理領域の上限
    &                nz           ! z 方向の物理領域の上限
  use axesset, only: x_X,        &!X 座標軸(スカラー格子点)
    &                y_Y,        &!Y 座標軸(スカラー格子点)
    &                z_Z,        &!Z 座標軸(スカラー格子点)
    &                p_X,        &!X 座標軸(フラックス格子点)
    &                q_Y,        &!Y 座標軸(フラックス格子点)
    &                r_Z,        &!Z 座標軸(フラックス格子点)
    &                x_dx, y_dy, z_dz, &! 格子間隔
    &                XMax,          &!X 座標の最大値
    &                YMax,          &!Y 座標の最大値
    &                ZMax            !Z 座標の最大値 
  use namelist_util, only: namelist_filename

  !暗黙の型宣言禁止
  implicit none

  !private 属性を指定
  private 
  
  !関数には public 属性を指定
  public Damping_Init
  public xyz_DampSponge_xyz
  public xyr_DampSponge_xyr
  public pyz_DampSponge_pyz
  public xqz_DampSponge_xqz
  public xyz_DampSponge
  public xyr_DampSponge
  public pyz_DampSponge
  public xqz_DampSponge
  public DampSponge_xyz
  public DampSponge_xyr
  public DampSponge_pyz
  public DampSponge_xqz

  !変数定義
  real(DP), save :: EFTime     = 100.0d0 !スポンジ層の e-folding time
  real(DP), save :: DampDepthH = 0.0d0   !スポンジ層の厚さ(水平方向)
  real(DP), save :: DampDepthV = 0.0d0   !スポンジ層の厚さ(鉛直方向)
  real(DP), allocatable, save :: xyz_DampRateH(:,:,:) !xyz 格子減衰係数(水平方向)
  real(DP), allocatable, save :: xyz_DampRateV(:,:,:) !xyz 格子減衰係数(鉛直方向)
  real(DP), allocatable, save :: pyz_DampRateH(:,:,:) !pyz 格子減衰係数(水平方向)
  real(DP), allocatable, save :: pyz_DampRateV(:,:,:) !pyz 格子減衰係数(鉛直方向)
  real(DP), allocatable, save :: xqz_DampRateH(:,:,:) !xqz 格子減衰係数(水平方向)
  real(DP), allocatable, save :: xqz_DampRateV(:,:,:) !xqz 格子減衰係数(鉛直方向)
  real(DP), allocatable, save :: xyr_DampRateH(:,:,:) !xyr 格子減衰係数(水平方向)
  real(DP), allocatable, save :: xyr_DampRateV(:,:,:) !xyr 格子減衰係数(鉛直方向)

contains 
  
!!!------------------------------------------------------------------------!!!
  subroutine Damping_Init
    !
    ! 音波減衰項とスポンジ層の減衰係数の初期化
    ! 

    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(DP)                  :: Time     !
    real(DP)                  :: DepthH   !スポンジ層の厚さ(水平方向)
    real(DP)                  :: DepthV   !スポンジ層の厚さ(鉛直方向)
    real(DP), parameter       :: Pi =3.1415926535897932385d0   !円周率
    integer                  :: unit ! 装置番号
    integer                   :: i, j, k

    !NAMELIST から取得
    NAMELIST /damping_nml/ Time, DepthH, DepthV

    call FileOpen(unit, file=namelist_filename, mode='r')
    read(unit, NML=damping_nml)
    close(unit)

    !初期化
    allocate( &
      & xyz_DampRateH(imin:imax,jmin:jmax,kmin:kmax), &
      & xyz_DampRateV(imin:imax,jmin:jmax,kmin:kmax), &
      & pyz_DampRateH(imin:imax,jmin:jmax,kmin:kmax), &
      & pyz_DampRateV(imin:imax,jmin:jmax,kmin:kmax), &
      & xqz_DampRateH(imin:imax,jmin:jmax,kmin:kmax), &
      & xqz_DampRateV(imin:imax,jmin:jmax,kmin:kmax), &
      & xyr_DampRateH(imin:imax,jmin:jmax,kmin:kmax), &
      & xyr_DampRateV(imin:imax,jmin:jmax,kmin:kmax)    )
    xyz_DampRateH = 0.0d0
    xyz_DampRateV = 0.0d0
    pyz_DampRateH = 0.0d0
    pyz_DampRateV = 0.0d0
    xqz_DampRateH = 0.0d0
    xqz_DampRateV = 0.0d0
    xyr_DampRateH = 0.0d0
    xyr_DampRateV = 0.0d0

    !値の入力
    EFTime     = Time
    DampDepthH = DepthH
    DampDepthV = DepthV
    
    !-----------------------------------------------------------------    
    ! スポンジ層の減衰率
    !-----------------------------------------------------------------
    !水平方向の東側・西側境界
    if ( DampDepthH < x_dx(1) ) then 
      call MessageNotify( "W", &
        & "DampSponge_Init", "DampDepthH is too thin. DelX is %f", d=(/x_dx(1)/))

    else if ( DampDepthH < x_dx(nx) ) then 
      call MessageNotify( "W", &
        & "DampSponge_Init", "DampDepthH is too thin. DelX is %f", d=(/x_dx(nx)/))

    else
      do i = imin, imax
        !スカラー格子点の西側境界
        if ( x_X(i) < DampDepthH) then 
          xyz_DampRateH(i,:,:) = ((1.0d0 - x_X(i) / DampDepthH) ** 3.0d0) / EFTime
        end if
        
        !フラックス格子点の西側境界
        if ( p_X(i) < DampDepthH) then 
          pyz_DampRateH(i,:,:) = ((1.0d0 - p_X(i) / DampDepthH) ** 3.0d0) / EFTime
         end if
        
        !スカラー格子点の東側境界    
        if ( x_X(i) > ( XMax - DampDepthH ) ) then 
          xyz_DampRateH(i,:,:) = &
            & ((1.0d0 - (XMax - x_X(i)) / DampDepthH) ** 3.0d0) / EFTime 
        end if
        
        !フラックス格子点の東側境界    
        if ( p_X(i) > ( XMax - DampDepthH ) ) then 
          pyz_DampRateH(i,:,:) = &
            & ((1.0d0 - (XMax - p_X(i)) / DampDepthH) ** 3.0d0) / EFTime 
        end if
      end do
    end if
    !sf と ss は X 方向に関しては同じ
    xyr_DampRateH  = xyz_DampRateH

    !水平方向の南側・北側境界
    if ( DampDepthH < y_dy(1) ) then 
      call MessageNotify( "W", &
        & "DampSponge_Init", "DampDepthH is too thin. DelY is %f", d=(/x_dx(1)/))

    else if ( DampDepthH < y_dy(ny) ) then 
      call MessageNotify( "W", &
        & "DampSponge_Init", "DampDepthH is too thin. DelY is %f", d=(/y_dy(ny)/))

    else
      do j = jmin, jmax
        !スカラー格子点の西側境界
        if ( y_Y(j) < DampDepthH) then 
          xyz_DampRateH(:,j,:) = ((1.0d0 - y_Y(j) / DampDepthH) ** 3.0d0) / EFTime
        end if
        
        !フラックス格子点の西側境界
        if ( q_Y(j) < DampDepthH) then 
          pyz_DampRateH(:,j,:) = ((1.0d0 - q_Y(j) / DampDepthH) ** 3.0d0) / EFTime
         end if
        
        !スカラー格子点の東側境界    
        if ( y_Y(j) > ( YMax - DampDepthH ) ) then 
          xyz_DampRateH(:,j,:) = &
            & ((1.0d0 - (YMax - y_Y(j)) / DampDepthH) ** 3.0d0) / EFTime 
        end if
        
        !フラックス格子点の東側境界    
        if ( q_Y(j) > ( YMax - DampDepthH ) ) then 
          pyz_DampRateH(:,j,:) = &
            & ((1.0d0 - (YMax - q_Y(j)) / DampDepthH) ** 3.0d0) / EFTime 
        end if
      end do
    end if
    !sf と ss は X 方向に関しては同じ
    xqz_DampRateH  = xyz_DampRateH

    
    !鉛直方向の上部境界    
    if ( DampDepthV < z_dz(nz) ) then 
      call MessageNotify( "W", &
        & "DampSponge_Init", "DampDepthV is too thin. DelZ is %f", d=(/z_dz(nz)/) )
      
    else
      do k = kmin, kmax
        !スカラー格子点
        if ( z_Z(k) >= ( ZMax - DampDepthV ) ) then 
          xyz_DampRateV(:,:,k) =  &
            & (1.0d0 - dcos(Pi * (z_Z(k) - ZMax + DampDepthV) / DampDepthV)) &
            &  / EFTime 
        end if
        
        !フラックス格子点
        if ( r_Z(k) >= ( ZMax - DampDepthV ) ) then 
          xyr_DampRateV(:,:,k) =  &
            & (1.0d0 - dcos(Pi * (r_Z(k) - ZMax + DampDepthV)/ DampDepthV)) &
            &  / EFTime 
        end if
      end do
    end if
    !fs と ss は Z 方向に関しては同じ
    pyz_DampRateV  = xyz_DampRateV
    
    !-----------------------------------------------------------------    
    ! 値の確認
    !-----------------------------------------------------------------

    call MessageNotify( "M", "DampSponge_Init", "EFTime = %f", d=(/EFTime/) )
    call MessageNotify( "M", "DampSponge_Init", "DampDepthH = %f", d=(/DampDepthH/) )
    call MessageNotify( "M", "DampSponge_Init", "DampDepthV = %f", d=(/DampDepthV/) )  

  end subroutine Damping_Init
!!!------------------------------------------------------------------------!!!
  function xyz_DampSponge_xyz(xyz_Var)
    !
    ! xyz 格子点に対するスポンジ層
    !

    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(DP), intent(in)   :: xyz_Var(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)               :: xyz_DampSponge_xyz(imin:imax,jmin:jmax,kmin:kmax)
    
    !スポンジ層によるダンピングを計算
    xyz_DampSponge_xyz =  ( xyz_DampRateH + xyz_DampRateV ) * xyz_Var
    
  end function xyz_DampSponge_xyz
!!!------------------------------------------------------------------------!!!
  function xyz_DampSponge(xyz_VarA, xyz_VarB)
    !
    ! xyz 格子点に対するスポンジ層
    !

    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(DP), intent(in)   :: xyz_VarA(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)   :: xyz_VarB(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)               :: xyz_DampSponge(imin:imax,jmin:jmax,kmin:kmax)
    
    !スポンジ層によるダンピングを計算
    xyz_DampSponge =  xyz_VarA - ( xyz_DampRateH + xyz_DampRateV ) * xyz_VarB * (2.0d0 * DelTimeLong)
    
  end function xyz_DampSponge


!!!------------------------------------------------------------------------!!!
  function xyr_DampSponge_xyr(xyr_Var)
    !
    ! xyr 格子点に対するスポンジ層
    !
    
    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(DP), intent(in)   :: xyr_Var(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)               :: xyr_DampSponge_xyr(imin:imax,jmin:jmax,kmin:kmax)

    !スポンジ層によるダンピングを計算  
    xyr_DampSponge_xyr =  ( xyr_DampRateH + xyr_DampRateV )* xyr_Var
    
  end function xyr_DampSponge_xyr


!!!------------------------------------------------------------------------!!!
  function xyr_DampSponge(xyr_VarA, xyr_VarB)
    !
    ! xyr 格子点に対するスポンジ層
    !
    
    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(DP), intent(in)   :: xyr_VarA(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)   :: xyr_VarB(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)               :: xyr_DampSponge(imin:imax,jmin:jmax,kmin:kmax)

    !スポンジ層によるダンピングを計算  
    xyr_DampSponge =  xyr_VarA - ( xyr_DampRateH + xyr_DampRateV )* xyr_VarB * (2.0d0 * DelTimeLong)
    
  end function xyr_DampSponge


!!!------------------------------------------------------------------------!!!
  function pyz_DampSponge_pyz(pyz_Var)
    !
    ! pyz 格子点に対するスポンジ層
    !

    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(DP), intent(in)   :: pyz_Var(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)               :: pyz_DampSponge_pyz(imin:imax,jmin:jmax,kmin:kmax)
    
    !スポンジ層によるダンピングを計算  
    pyz_DampSponge_pyz = ( pyz_DampRateH + pyz_DampRateV ) * pyz_Var
    
  end function pyz_DampSponge_pyz
  

!!!------------------------------------------------------------------------!!!
  function pyz_DampSponge(pyz_VarA, pyz_VarB)
    !
    ! pyz 格子点に対するスポンジ層
    !

    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(DP), intent(in)   :: pyz_VarA(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)   :: pyz_VarB(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)               :: pyz_DampSponge(imin:imax,jmin:jmax,kmin:kmax)
    
    !スポンジ層によるダンピングを計算  
    pyz_DampSponge = pyz_VarA - ( pyz_DampRateH + pyz_DampRateV ) * pyz_VarB * (2.0d0 * DelTimeLong)
    
  end function pyz_DampSponge
!!!------------------------------------------------------------------------!!!
  function xqz_DampSponge_xqz(xqz_Var)
    !
    ! xqz 格子点に対するスポンジ層
    !

    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(DP), intent(in)   :: xqz_Var(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)               :: xqz_DampSponge_xqz(imin:imax,jmin:jmax,kmin:kmax)
    
    !スポンジ層によるダンピングを計算  
    xqz_DampSponge_xqz = ( xqz_DampRateH + xqz_DampRateV ) * xqz_Var
    
  end function xqz_DampSponge_xqz
!!!------------------------------------------------------------------------!!!
  function xqz_DampSponge(xqz_VarA, xqz_VarB)
    !
    ! xqz 格子点に対するスポンジ層
    !

    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(DP), intent(in)   :: xqz_VarA(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)   :: xqz_VarB(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)               :: xqz_DampSponge(imin:imax,jmin:jmax,kmin:kmax)
    
    !スポンジ層によるダンピングを計算  
    xqz_DampSponge = xqz_VarA - ( xqz_DampRateH + xqz_DampRateV ) * xqz_VarB * (2.0d0 * DelTimeLong)
    
  end function xqz_DampSponge

!!!------------------------------------------------------------------------!!!
  subroutine DampSponge_xyz(xyz_VarA, xyz_VarB)
    !
    ! ss 格子点に対するスポンジ層
    !

    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(8), intent(inout):: xyz_VarA(imin:imax, jmin:jmax, kmin:kmax)
    real(8), intent(in)   :: xyz_VarB(imin:imax, jmin:jmax, kmin:kmax)
    real(8)               :: xyz_Var(imin:imax, jmin:jmax, kmin:kmax)
    
    !スポンジ層によるダンピングを計算
    xyz_Var  = xyz_VarA - ( xyz_DampRateH + xyz_DampRateV ) * xyz_VarB * (2.0d0 * DelTimeLong)
    xyz_VarA = xyz_Var
    
  end subroutine DampSponge_xyz


!!!------------------------------------------------------------------------!!!
  subroutine DampSponge_xyr(xyr_VarA, xyr_VarB)
    !
    ! sf 格子点に対するスポンジ層
    !
    
    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(8), intent(inout):: xyr_VarA(imin:imax, jmin:jmax, kmin:kmax)
    real(8), intent(in)   :: xyr_VarB(imin:imax, jmin:jmax, kmin:kmax)
    real(8)               :: xyr_Var(imin:imax, jmin:jmax, kmin:kmax)

    !スポンジ層によるダンピングを計算  
    xyr_Var  = xyr_VarA - ( xyr_DampRateH + xyr_DampRateV )* xyr_VarB * (2.0d0 * DelTimeLong)
    xyr_VarA = xyr_Var
    
  end subroutine DampSponge_xyr
  

!!!------------------------------------------------------------------------!!!
  subroutine DampSponge_pyz(pyz_VarA, pyz_VarB)
    !
    ! fs 格子点に対するスポンジ層
    !

    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(8), intent(inout):: pyz_VarA(imin:imax, jmin:jmax, kmin:kmax)
    real(8), intent(in)   :: pyz_VarB(imin:imax, jmin:jmax, kmin:kmax)
    real(8)               :: pyz_Var(imin:imax, jmin:jmax, kmin:kmax)
    
    !スポンジ層によるダンピングを計算  
    pyz_Var  = pyz_VarA - ( pyz_DampRateH + pyz_DampRateV ) * pyz_VarB * (2.0d0 * DelTimeLong)
    pyz_VarA = pyz_Var
    
  end subroutine DampSponge_pyz
  
!!!------------------------------------------------------------------------!!!
  subroutine DampSponge_xqz(xqz_VarA, xqz_VarB)
    !
    ! fs 格子点に対するスポンジ層
    !

    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(8), intent(inout):: xqz_VarA(imin:imax, jmin:jmax, kmin:kmax)
    real(8), intent(in)   :: xqz_VarB(imin:imax, jmin:jmax, kmin:kmax)
    real(8)               :: xqz_Var(imin:imax, jmin:jmax, kmin:kmax)
    
    !スポンジ層によるダンピングを計算  
    xqz_Var  = xqz_VarA - ( xqz_DampRateH + xqz_DampRateV ) * xqz_VarB * (2.0d0 * DelTimeLong)
    xqz_VarA = xqz_Var
    
  end subroutine DampSponge_xqz
  
  
end module Damping
