!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Module DisturbSet_calc
!
!   * Developer: SUGIYAMA Ko-ichiro (sugiyama@gfd-dennou.org)
!   * Version: $Id: disturbset.f90,v 1.2.2.1 2005/07/07 10:36:12 kitamo Exp $ 
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview 
!
!擾乱場の設定を読み込み, 擾乱場の値を計算する. 
!
!典型的な擾乱場として以下が設定可能. 
!   * 無次元圧力: 擾乱なし or X 方向ガウシアン分布 or Z 方向ガウシアン分布 or XZ 方向ガウシアン分布 
!   * 温位: 擾乱なし or X 方向ガウシアン分布 or Z 方向ガウシアン分布 or XZ 方向ガウシアン分布 
!   * 速度: 擾乱なし or X 方向速度一定 or Z 方向速度一定 
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!== Future Plans
!
!  * 設定可能な擾乱のタイプを増やす
!  * エラー処理に gf4f90io を利用
!
                                                                 !=end

module disturbset
  !== 暗黙の型宣言禁止
  implicit none

  !== save 属性
  save
                                                                 !=begin
  !== Public Interface
  real(8), allocatable  :: ss_PotTempDisturb(:, :) !温位
  real(8), allocatable  :: ss_ExnerDisturb(:, :)   !エクスナー関数
  real(8), allocatable  :: fs_VelXDisturb(:, :)    !速度 u
  real(8), allocatable  :: sf_VelZDisturb(:, :)    !速度 w
                                                                 !=end

contains
                                                                 !=begin
  !== Procedure Interface 
  !
  !=== Initialize module and acquire NAMELIST
  !
  !NAMELIST から擾乱場の初期値生成プログラム disturb の設定を読み取る. 
  !
  subroutine disturbset_calc_init(cfgfile)
    !=== Dependency
    use dc_trace,   only: BeginSub, EndSub 
    use dc_message, only: MessageNotify
    use gridset,    only: s_X, s_Z, f_X, f_Z, Xmin, Xmax, Zmin, Zmax, &
         &                DimXMin, DimXMax, DimZMin, DimZMax, RegZMin
    use physset,    only: MolWtDry, CpDry, Pi
    use fileset,    only: RandomFile
    
    !=== Input
    character(*), intent(in) :: cfgfile
                                                                 !=end
    !=== Work
    character(80) :: ExnerDisturbType !無次元圧力の擾乱 (DisturbType.f90 参照)
    real(8)       :: ExnerSigmaX      !擾乱の幅 (X 方向)
    real(8)       :: ExnerSigmaZ      !擾乱の幅 (Z 方向)
    real(8)       :: ExnerMax         !擾乱の最大値
    real(8)       :: ExnerCenterX     !擾乱の中心 (X 方向)
    real(8)       :: ExnerCenterZ     !擾乱の中心 (X 方向)
    
    character(80) :: PotTempDisturbType  !温位の擾乱 (DisturbType.f90 参照)
    real(8)       :: PotTempSigmaX       !擾乱の幅 (X 方向)
    real(8)       :: PotTempSigmaZ       !擾乱の幅 (X 方向)
    real(8)       :: PotTempMax          !擾乱の最大値
    real(8)       :: PotTempCenterX      !擾乱の中心 (X 方向)
    real(8)       :: PotTempCenterZ      !擾乱の中心 (X 方向)
    
    character(80) :: VelDisturbType   !速度の擾乱 (DisturbType.f90 参照)
    real(8)       :: VelX             !速度 u 
    real(8)       :: VelZ             !速度 w

    integer :: i, k
    real(8) :: mu
    real(8) :: sigma
    real(8) :: random
                                                                 !=end
                                                                 !=begin
    !=== NAMELIST
    NAMELIST /disturbset/ &
         & ExnerDisturbType, ExnerMax, ExnerSigmaX, ExnerSigmaZ, &
         &                  ExnerCenterX, ExnerCenterZ, &
         & PotTempDisturbType, PotTempMax, PotTempSigmaX, PotTempSigmaZ, &
         &                  PotTempCenterX, PotTempCenterZ, &
         & VelDisturbType, VelX, VelZ
                                                                 !=end    
    call BeginSub("disturbset_calc_init", &
&                 fmt="%c",       &
&                 c1="Initialize values of disturbance. ")


    open (10, FILE=cfgfile)
    read(10, NML=disturbset)
    close(10)

    !==== 変数の値のチェック

    if ( ExnerSigmaX > 1 )  then 
       call MessageNotify("Message", "disturbset_calc_init", &
&                         "SigmaX is > 1")
       stop
    end if

    if ( ExnerSigmaZ > 1 )  then
       call MessageNotify("Message", "disturbset_calc_init", &
&                         "SigmaZ is > 1")
       stop
    end if

    if ( PotTempSigmaX > 1 )  then 
       call MessageNotify("Message", "disturbset_calc_init", &
&                         "SigmaX is > 1")
       stop
    end if

    if ( PotTempSigmaZ > 1 )  then
       call MessageNotify("Message", "disturbset_calc_init", &
&                         "SigmaZ is > 1")
       stop
    end if

    if ( ExnerCenterX > 1 ) then 
       call MessageNotify("Message", "disturbset_calc_init", &
&                         "CenterX is > 1")
       stop
    end if

    if ( ExnerCenterZ > 1 ) then
       call MessageNotify("Message", "disturbset_calc_init", &
&                         "CenterZ is > 1")
       stop
    end if

    if ( PotTempCenterX > 1 ) then 
       call MessageNotify("Message", "disturbset_calc_init", &
&                         "CenterX is > 1")
       stop
    end if

    if ( PotTempCenterZ > 1 ) then
       call MessageNotify("Message", "disturbset_calc_init", &
&                         "CenterZ is > 1")
       stop
    end if


    !==== 配列の割り当て, 初期化
    allocate( ss_PotTempDisturb(DimXMin:DimXMax, DimZMin:DimZMax), &
         &    ss_ExnerDisturb(DimXMin:DimXMax, DimZMin:DimZMax),   &
         &    fs_VelXDisturb(DimXMin:DimXMax, DimZMin:DimZMax),    &
         &    sf_VelZDisturb(DimXMin:DimXMax, DimZMin:DimZMax)  )

    ss_ExnerDisturb = 0.0d0
    ss_PotTempDisturb = 0.0d0
    fs_VelXDisturb = 0.0d0
    sf_VelZDisturb = 0.0d0
    

    !-----------------------------------------------------------
    !--- エクスナー関数の擾乱
    !-----------------------------------------------------------
    select case (ExnerDisturbType)
       
    case ("Exner_Gauss_X")
       !--- X 方向のガウシアン分布
       mu = ( Xmax - Xmin ) * ExnerCenterX
       sigma = ( Xmax - Xmin ) * ExnerSigmaX
       
       do i = DimXMin, DimXMax
          ss_ExnerDisturb(i,:) = &
               & ExnerMax * dexp(-((s_X(i) - mu) / sigma )**2.0d0 / 2.0d0) 
       end do
       
    case ("Exner_Gauss_Z")    
       !--- Z 方向のガウシアン分布
       mu = ( Zmax - Zmin ) * ExnerCenterZ
       sigma = ( Zmax - Zmin ) * ExnerSigmaZ
       
       do k = DimZMin, DimZMax
          ss_ExnerDisturb(:,k) = &
               & ExnerMax * dexp(-((s_Z(k) - mu) / sigma)**2.0d0 / 2.0d0) 
       end do
       
    case ("Exner_Gauss_XZ")
       !--- 同心円状のガウシアン分布
       mu = ( Xmax - Xmin ) * ExnerCenterX
       sigma = ( Xmax - Xmin ) * ExnerSigmaX
       
       do i = DimXMin, DimXMax
          ss_ExnerDisturb(i,:) = &
               & ExnerMax * dexp(-((s_X(i) - mu) / sigma)**2.0d0 / 2.0d0) 
       end do
       
       mu = ( Zmax - Zmin ) * ExnerCenterZ
       sigma = ( Zmax - Zmin ) * ExnerSigmaZ
       
       do k = DimZMin, DimZMax
          ss_ExnerDisturb(:,k) = ss_ExnerDisturb(:,k) &
               & * dexp(-((s_Z(k) - mu) / sigma)**2.0d0 / 2.0d0) 
       end do
    end select

    !-----------------------------------------------------------
    !--- 仮温位の擾乱
    !-----------------------------------------------------------
    select case (PotTempDisturbType)
       
    case ("PotTemp_Gauss_X")
       !--- X 方向のガウシアン分布
       mu = ( Xmax - Xmin ) * PotTempCenterX
       sigma = ( Xmax - Xmin ) * PotTempSigmaX
       
       do i = DimXMin, DimXMax
          ss_PotTempDisturb(i,:) = &
               & PotTempMax * dexp(-((s_X(i) - mu) / sigma)**2.0d0 / 2.0d0) 
       end do
       
    case ("PotTemp_Gauss_Z")
       !--- Z 方向のガウシアン分布
       mu = ( Zmax - Zmin ) * PotTempCenterX
       sigma = ( Zmax - Zmin ) * PotTempSigmaZ
       
       do k = DimZMin, DimZMax
          ss_PotTempDisturb(:,k) = &
               & PotTempMax * dexp(-((s_Z(k) - mu) / sigma)**2.0d0 / 2.0d0) 
       end do
       
    case ("PotTemp_Gauss_XZ")     
       !--- 同心円状のガウシアン分布
       mu = ( Xmax - Xmin ) * PotTempCenterX
       sigma = ( Xmax - Xmin ) * PotTempSigmaX
       
       do i = DimXMin, DimXMax
          ss_PotTempDisturb(i,:) = &
               & PotTempMax * dexp(-((s_X(i) - mu) / sigma)**2.0d0 / 2.0d0) 
       end do
       
       mu = ( Zmax - Zmin ) * PotTempCenterZ
       sigma = ( Zmax - Zmin ) * PotTempSigmaZ
       
       do k = DimZMin, DimZMax
          ss_PotTempDisturb(:,k) =  &
               & ss_PotTempDisturb(:,k) &
               &   * dexp(-((s_Z(k) - mu) / sigma)**2.0d0 / 2.0d0) 
       end do
       
!       where (ss_PotTempDisturb < PotTempMax * 1.0d-4)
!          ss_PotTempDisturb = 0.0d0 
!       end where
       
    case ("RandomLowLayer")
        !--- 最下層にランダムな擾乱を与える

       write(*,*) RandomFile
       open(20,file=RandomFile)

       do i = DimXMin, DimXMax

          read(20,*) random
          ss_PotTempDisturb(i,RegZMin + 1) = &
               & PotTempMax * random
       end do

       close(20)

    case default
       
       ss_PotTempDisturb = 0.0d0
       
    end select

     
    !----------------------------------------------------------
    !--- 速度
    !----------------------------------------------------------
    select case (VelDisturbType)
       
    case ("Zero")
       fs_VelXDisturb = 0.0d0
       sf_VelZDisturb = 0.0d0
       
    case ("VelX_test")
       fs_VelXDisturb = VelX
       fs_VelXDisturb(-5:10, :) = 0.0d0
       fs_VelXDisturb(90:105, :) = 0.0d0
       sf_VelZDisturb = 0.0d0
       
    case ("VelX_Uniform")
       !--- 速度 u 一定
       fs_VelXDisturb = VelX
       sf_VelZDisturb = 0.0d0
       
    case ("VelZ_Uniform")
       !--- 速度 w 一定
       fs_VelXDisturb = 0.0d0
       sf_VelZDisturb = VelZ
       
    case ("Solid_Rotation")
       !--- 領域中心の剛体回転
       
       do i = DimXMin, DimXMax   
          fs_VelXDisturb(i,:) = &
               & 2.0d0*pi / 1000.0d0 * ( f_Z(:) - (Zmax + Zmin) / 2.0D0 )
       end do
       
       do k = DimZMin, DimZMax
          sf_VelZDisturb(:,k) = - 2.0d0 * pi / 1000.0d0 &
               &                 * ( s_X(:) - (Xmax + Xmin) / 2.0D0 )
       end do
       
    end select
    
    call EndSub("disturbset_calc_init")
    
  end subroutine disturbset_calc_init
  
end module disturbset
