!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Module ArareSet
!
!   * Developer: SUGIYAMA Ko-ichiro (sugiyama@gfd-dennou.org)
!   * Version: $Id: arareset.f90,v 1.2 2005/04/22 15:01:52 sugiyama Exp $ 
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview 
!
!arare を実行するために必要なパラメタを NAMELIST から取得する. 
!NAMELIST 名は外部から取得する. 
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!== Future Plans
!
                                                                 !=end

module arareset

  implicit none

  !--- save 属性
  save
                                                                 !=begin
  !== Public Interface 
  !
  real(8)  :: alpha  !音波減衰項の係数
  real(8)  :: beta   !クランクニコルソン法なら 0.5, 完全陰解放なら 1
  real(8)  :: Cm     !乱流エネルギー診断式の係数 
  real(8)  :: tfil   !アッセリンの時間フィルタの係数
  real(8)  :: NuH    !数値粘性の係数 (水平方向)
  real(8)  :: NuV    !数値粘性の係数 (鉛直方向)
  real(8)  :: Min    !混合比の打ち切り
  real(8)  :: EFTime  !スポンジ層の e-folding time
  real(8)  :: SDepthH !スポンジ層の厚さ(水平方向)
  real(8)  :: SDepthV    !スポンジ層の厚さ(鉛直方向)
  real(8), allocatable :: ss_SDumpH(:,:) !スポンジ層の厚さ(水平方向)
  real(8), allocatable :: ss_SDumpV(:,:) !スポンジ層の厚さ(水平方向)
  real(8), allocatable :: fs_SDumpH(:,:) !スポンジ層の厚さ(水平方向)
  real(8), allocatable :: fs_SDumpV(:,:) !スポンジ層の厚さ(水平方向)
  real(8), allocatable :: sf_SDumpH(:,:) !スポンジ層の厚さ(水平方向)
  real(8), allocatable :: sf_SDumpV(:,:) !スポンジ層の厚さ(水平方向)
                                                                 !=end

contains
                                                                 !=begin
  !== Procedure Interface 
  !
  !=== Initialize module and acquire NAMELIST
  !
  !NameList ファイルから arare 実行に必要なパラメタを取得する
  !
  subroutine arareset_init(cfgfile)
    !=== Dependency
    use dc_trace, only: BeginSub, EndSub 
    use physset,  only: Pi
    use gridset,  only: DelX, DelZ, XMax, ZMax, s_X, f_X, s_Z, f_Z, &
         &              DimXMin,DimXMax, DimZMin, DimZMax
    use timeset,  only: DelTimeShort, DelTimeLong
    
    !=== Input 
    character(*), intent(in) :: cfgfile
                                                                 !=end
    !--- 内部変数
    real(8)  ::  AlphaH
    real(8)  ::  AlphaV
    integer  ::  i, k
                                                                 !=begin
    !=== NAMELIST
    NAMELIST /arareset/ alpha, beta, Cm, tfil, AlphaH, AlphaV, &
      &         EFTime, SDepthH, SDepthV, Min 
                                                                 !=end

    call BeginSub("arareset_init", &
      &                 fmt="%c",       &
      &                 c1="Initialize parameters.")

    !--- 変数の取り出し
    open (10, FILE=cfgfile)
    read(10, NML=arareset)
    close(10)
    
    !--- 音波減衰項の係数を決める
    alpha = alpha * ( DelX ** 2.0d0 ) / DelTimeShort

    !--- 数値粘性の係数を決める
    NuH  = AlphaH * ( DelX ** 2.0d0 ) / DelTimeLong
    NuV  = AlphaV * ( DelZ ** 2.0d0 ) / DelTimeLong


    !--- もしもスポンジ層の厚さが格子点間隔より大きくない場合には警告
    !----- X 方向の厚さのチェック
    !----- Z 方向の厚さのチェック
    
    
    !--- 格子点毎のダンピング係数を決める
    !----- 初期化
    allocate( &
         & ss_SDumpH(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_SDumpV(DimXMin:DimXMax, DimZMin:DimZMax),   &
         & fs_SDumpH(DimXMin:DimXMax, DimZMin:DimZMax), &
         & fs_SDumpV(DimXMin:DimXMax, DimZMin:DimZMax),   &
         & sf_SDumpH(DimXMin:DimXMax, DimZMin:DimZMax), &
         & sf_SDumpV(DimXMin:DimXMax, DimZMin:DimZMax)    )
    ss_SDumpH = 0.0d0 
    ss_SDumpV = 0.0d0 
    fs_SDumpH = 0.0d0 
    fs_SDumpV = 0.0d0 
    sf_SDumpH = 0.0d0 
    sf_SDumpV = 0.0d0 

    !----- 水平方向の西側境界
    if ( SDepthH < DelX ) then 
      !警告を出すだけで何もしない. 
      write(*,*) "SDepthH is too thin. DelX is ", DelX

    else
      do i = DimXMin, DimXMax
        !----- スカラー格子点の西側境界
        if ( s_X(i) < SDepthH) then 
          ss_SDumpH(i,:) = ((1.0d0 - s_X(i) / SDepthH) ** 3.0d0) / EFTime
        end if
        
        !----- フラックス格子点の西側境界
        if ( f_X(i) < SDepthH) then 
          fs_SDumpH(i,:) = ((1.0d0 - f_X(i) / SDepthH) ** 3.0d0) / EFTime 
        end if
        
        !----- スカラー格子点の東側境界    
        if ( s_X(i) > ( XMax - SDepthH ) ) then 
          ss_SDumpH(i,:) = ((1.0d0 - (XMax - s_X(i)) / SDepthH) ** 3.0d0) &
            &              / EFTime 
        end if
        
        !----- フラックス格子点の東側境界    
        if ( f_X(i) > ( XMax - SDepthH ) ) then 
          fs_SDumpH(i,:) = ((1.0d0 - (XMax - f_X(i)) / SDepthH) ** 3.0d0) &
            &              / EFTime 
        end if
      end do
    end if
    !--- sf と ss は X 方向に関しては同じ
    sf_SDumpH  = ss_SDumpH
    
    !----- 鉛直方向の上部境界    
    if ( SDepthV < DelZ ) then 
      !警告を出すだけで何もしない. 
      write(*,*) "SDepthV is too thin. DelZ is ", DelZ
      
    else
      do k = DimZMin, DimZMax
        !----- スカラー格子点
        if ( s_Z(k) >= ( ZMax - SDepthV ) ) then 
          ss_SDumpV(:,k) =  &
            & (1.0d0 - dcos(Pi * (s_Z(k) - ZMax + SDepthV) / SDepthV)) &
            &  / EFTime 
        end if
        
        !----- フラックス格子点
        if ( f_Z(k) >= ( ZMax - SDepthV ) ) then 
          sf_SDumpV(:,k) =  &
            & (1.0d0 - dcos(Pi * (f_Z(k) - ZMax + SDepthV)/ SDepthV)) &
            &  / EFTime 
        end if
      end do
    end if
    !--- fs と ss は Z 方向に関しては同じ
    fs_SDumpV  = ss_SDumpV
    
    !==== 確認
!    write(*,*) "alpha", alpha
!    write(*,*) "beta", beta
!    write(*,*) "Cm", Cm
!    write(*,*) "tfil", tfil
!    write(*,*) "NuH", NuH
!    write(*,*) "NuV", NuV

    call EndSub("arareset_init")
    
  end subroutine arareset_init
  
end module arareset
