!= Module Damping
!
! Authors::   SUGIYAMA Ko-ichiro, ODAKA Masatsugu
! Version::   $Id: damping.f90,v 1.7 2011-10-05 03:35:23 sugiyama Exp $
! Tag Name::  $Name: arare5-20111010 $
! 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

  !暗黙の型宣言禁止
  implicit none

  !private 属性を指定
  private

  !関数には public 属性を指定
  public Damping_Init
  public SpongeLayer_forcing

  !変数定義
  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_Gamma(:,:,:) !xyz 格子減衰係数(水平方向)
  real(DP), allocatable, save :: pyz_Gamma(:,:,:) !pyz 格子減衰係数(鉛直方向)
  real(DP), allocatable, save :: xqz_Gamma(:,:,:) !xqz 格子減衰係数(鉛直方向)
  real(DP), allocatable, save :: xyr_Gamma(:,:,:) !xyr 格子減衰係数(鉛直方向)

contains

!!!------------------------------------------------------------------------!!!
  subroutine Damping_Init
    !
    ! 音波減衰項とスポンジ層の減衰係数の初期化
    !
    use dc_iounit,  only: FileOpen
    use dc_message, only: MessageNotify
    use gtool_historyauto, only: HistoryAutoAddVariable
    use mpi_wrapper, only: myrank
    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

    !変数定義
    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_Gamma(imin:imax,jmin:jmax,kmin:kmax), &
      & pyz_Gamma(imin:imax,jmin:jmax,kmin:kmax), &
      & xqz_Gamma(imin:imax,jmin:jmax,kmin:kmax), &
      & xyr_Gamma(imin:imax,jmin:jmax,kmin:kmax)    )
    xyz_Gamma = 0.0d0
    pyz_Gamma = 0.0d0
    xqz_Gamma = 0.0d0
    xyr_Gamma = 0.0d0

    !値の入力
    EFTime     = Time
    DampDepthH = DepthH
    DampDepthV = DepthV

    !-----------------------------------------------------------------
    ! スポンジ層の減衰率
    !
    !水平方向の東側・西側境界
    if ( DampDepthH < x_dx(1) ) then
      if (myrank == 0) &
        & call MessageNotify( "W", "Damping_init", "DampDepthH is too thin. DelX is %f", d=(/x_dx(1)/))

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

    else
      do i = imin, imax
        !スカラー格子点の西側境界
        if ( x_X(i) < DampDepthH) then
          xyz_Gamma(i,:,:) = ((1.0d0 - x_X(i) / DampDepthH) ** 3.0d0) / EFTime
        end if

        !フラックス格子点の西側境界
        if ( p_X(i) < DampDepthH) then
          pyz_Gamma(i,:,:) = ((1.0d0 - p_X(i) / DampDepthH) ** 3.0d0) / EFTime
        end if

        !スカラー格子点の東側境界
        if ( x_X(i) > ( XMax - DampDepthH ) ) then
          xyz_Gamma(i,:,:) = &
            & ((1.0d0 - (XMax - x_X(i)) / DampDepthH) ** 3.0d0) / EFTime
        end if

        !フラックス格子点の東側境界
        if ( p_X(i) > ( XMax - DampDepthH ) ) then
          pyz_Gamma(i,:,:) = &
            & ((1.0d0 - (XMax - p_X(i)) / DampDepthH) ** 3.0d0) / EFTime
        end if
      end do
    end if

    ! x 方向には同じ
    !
    xyr_Gamma  = xyz_Gamma
    xqz_Gamma  = xyz_Gamma

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

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

    else
      do j = jmin, jmax
        !スカラー格子点の西側境界
        if ( y_Y(j) < DampDepthH) then
          xyz_Gamma(:,j,:) = ((1.0d0 - y_Y(j) / DampDepthH) ** 3.0d0) / EFTime
        end if

        !フラックス格子点の西側境界
        if ( q_Y(j) < DampDepthH) then
          xqz_Gamma(:,j,:) = ((1.0d0 - q_Y(j) / DampDepthH) ** 3.0d0) / EFTime
         end if

        !スカラー格子点の東側境界
        if ( y_Y(j) > ( YMax - DampDepthH ) ) then
          xyz_Gamma(:,j,:) = &
            & ((1.0d0 - (YMax - y_Y(j)) / DampDepthH) ** 3.0d0) / EFTime
        end if

        !フラックス格子点の東側境界
        if ( q_Y(j) > ( YMax - DampDepthH ) ) then
          xqz_Gamma(:,j,:) = &
            & ((1.0d0 - (YMax - q_Y(j)) / DampDepthH) ** 3.0d0) / EFTime
        end if
      end do
    end if

    ! y 方向には同じ
    !
    pyz_Gamma  = xyz_Gamma
    xyr_Gamma  = xyz_Gamma

    !鉛直方向の上部境界
    if ( DampDepthV < z_dz(nz) ) then
      if (myrank == 0) &
        & call MessageNotify( "W", "Damping_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_Gamma(:,:,k) =  &
            & (1.0d0 - dcos(Pi * (z_Z(k) - ZMax + DampDepthV) / DampDepthV)) &
            &  / EFTime
        end if

        !フラックス格子点
        if ( r_Z(k) >= ( ZMax - DampDepthV ) ) then
          xyr_Gamma(:,:,k) =  &
            & (1.0d0 - dcos(Pi * (r_Z(k) - ZMax + DampDepthV)/ DampDepthV)) &
            &  / EFTime
        end if
      end do
    end if

    ! z 方向には同じ
    !
    pyz_Gamma  = xyz_Gamma
    xqz_Gamma  = xyz_Gamma

    !-----------------------------------------------------------------
    ! 値の確認
    !
    if (myrank == 0) then
      call MessageNotify( "M", "Damping_init", "EFTime = %f", d=(/EFTime/) )
      call MessageNotify( "M", "Damping_init", "DampDepthH = %f", d=(/DampDepthH/) )
      call MessageNotify( "M", "Damping_init", "DampDepthV = %f", d=(/DampDepthV/) )
    end if

    !-----------------------------------------------------------------
    ! 出力
    !
    call HistoryAutoAddVariable(                          &
      & varname='PTempDamp',                              &
      & dims=(/'x','y','z','t'/),                         &
      & longname='Damping term of potential temperature', &
      & units='K.s-1',                                    &
      & xtype='float')

    call HistoryAutoAddVariable(                    &
      & varname='ExnerDamp',                        &
      & dims=(/'x','y','z','t'/),                   &
      & longname='Damping term of exner function',  &
      & units='s-1',                                &
      & xtype='float')

    call HistoryAutoAddVariable(          &
      & varname='VelXDamp',               &
      & dims=(/'x','y','z','t'/),         &
      & longname='Damping term of VelX',  &
      & units='m.s-1',                    &
      & xtype='float')

    call HistoryAutoAddVariable(          &
      & varname='VelYDamp',               &
      & dims=(/'x','y','z','t'/),         &
      & longname='Damping term of VelY',  &
      & units='m.s-1',                    &
      & xtype='float')

    call HistoryAutoAddVariable(          &
      & varname='VelZDamp',               &
      & dims=(/'x','y','z','t'/),         &
      & longname='Damping term of VelZ',  &
      & units='m.s-1',                    &
      & xtype='float')

  end subroutine Damping_Init

  subroutine SpongeLayer_forcing(                                         &
    & pyz_VelXBl,  xqz_VelYBl,  xyr_VelZBl,  xyz_PTempBl,  xyz_ExnerBl,   & !(in)
    & pyz_VelXAl,  xqz_VelYAl,  xyr_VelZAl,  xyz_PTempAl,  xyz_ExnerAl )    !(inout)

    use gtool_historyauto, only: HistoryAutoPut
    use timeset, only: TimeN, DelTimeLong
    use gridset, only: imin,       &! x 方向の配列の下限
      &                imax,       &! x 方向の配列の上限
      &                jmin,       &! y 方向の配列の下限
      &                jmax,       &! y 方向の配列の上限
      &                kmin,       &! z 方向の配列の下限
      &                kmax,       &! z 方向の配列の上限
      &                nx,         &! x 方向の物理領域の上限
      &                ny,         &! y 方向の物理領域の上限
      &                nz           ! z 方向の物理領域の上限

    !暗黙の型宣言禁止
    implicit none

    real(8), intent(in)    :: pyz_VelXBl(imin:imax, jmin:jmax, kmin:kmax)
    real(8), intent(in)    :: xqz_VelYBl(imin:imax, jmin:jmax, kmin:kmax)
    real(8), intent(in)    :: xyr_VelZBl(imin:imax, jmin:jmax, kmin:kmax)
    real(8), intent(in)    :: xyz_PTempBl(imin:imax, jmin:jmax, kmin:kmax)
    real(8), intent(in)    :: xyz_ExnerBl(imin:imax, jmin:jmax, kmin:kmax)
    real(8), intent(inout) :: pyz_VelXAl(imin:imax, jmin:jmax, kmin:kmax)
    real(8), intent(inout) :: xqz_VelYAl(imin:imax, jmin:jmax, kmin:kmax)
    real(8), intent(inout) :: xyr_VelZAl(imin:imax, jmin:jmax, kmin:kmax)
    real(8), intent(inout) :: xyz_PTempAl(imin:imax, jmin:jmax, kmin:kmax)
    real(8), intent(inout) :: xyz_ExnerAl(imin:imax, jmin:jmax, kmin:kmax)
    real(8)                :: pyz_VelXAl0(imin:imax, jmin:jmax, kmin:kmax)
    real(8)                :: xqz_VelYAl0(imin:imax, jmin:jmax, kmin:kmax)
    real(8)                :: xyr_VelZAl0(imin:imax, jmin:jmax, kmin:kmax)
    real(8)                :: xyz_PTempAl0(imin:imax, jmin:jmax, kmin:kmax)
    real(8)                :: xyz_ExnerAl0(imin:imax, jmin:jmax, kmin:kmax)
    real(8)                :: pyz_DampVelX(imin:imax, jmin:jmax, kmin:kmax)
    real(8)                :: xqz_DampVelY(imin:imax, jmin:jmax, kmin:kmax)
    real(8)                :: xyr_DampVelZ(imin:imax, jmin:jmax, kmin:kmax)
    real(8)                :: xyz_DampPTemp(imin:imax, jmin:jmax, kmin:kmax)
    real(8)                :: xyz_DampExner(imin:imax, jmin:jmax, kmin:kmax)

    pyz_VelXAl0  =   pyz_VelXAl
    pyz_DampVelX = - xyz_Gamma * pyz_VelXBl
    pyz_VelXAl   =   pyz_VelXAl0 + (2.0d0 * DelTimeLong) * pyz_DampVelX

    xqz_VelYAl0  =   xqz_VelYAl
    xqz_DampVelY = - xqz_Gamma * xqz_VelYBl
    xqz_VelYAl   =   xqz_VelYAl0 + (2.0d0 * DelTimeLong) * xqz_DampVelY

    xyr_VelZAl0  =   xyr_VelZAl
    xyr_DampVelZ = - xyr_Gamma * xyr_VelZBl
    xyr_VelZAl   =   xyr_VelZAl0 + (2.0d0 * DelTimeLong) * xyr_DampVelZ

    xyz_PTempAl0  =   xyz_PTempAl
    xyz_DampPTemp = - xyz_Gamma * xyz_PTempBl
    xyz_PTempAl   =   xyz_PTempAl0 + (2.0d0 * DelTimeLong) * xyz_DampPTemp

    xyz_ExnerAl0  =   xyz_ExnerAl
    xyz_DampExner = - xyz_Gamma * xyz_ExnerBl
    xyz_ExnerAl   =   xyz_ExnerAl0 + (2.0d0 * DelTimeLong) * xyz_DampExner

    call HistoryAutoPut(TimeN, 'VelXDamp',  pyz_DampVelX(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelYDamp',  xqz_DampVelY(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelZDamp',  xyr_DampVelZ(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'PTempDamp', xyz_DampPTemp(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'ExnerDamp', xyz_DampExner(1:nx,1:ny,1:nz))

  end subroutine SpongeLayer_forcing

end module Damping