!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004, 2005. All rights reserved.
!---------------------------------------------------------------------
!= Module DisturbEnv
!
!   * Developer: SUGIYAMA Ko-ichiro, ODAKA Masatsugu
!   * Version: $Id: distset.f90,v 1.1 2011-03-30 13:53:06 sugiyama Exp $ 
!   * Tag Name: $Name: arare4-20120511 $
!   * Change History: 
!
!== Overview 
!
!擾乱のデフォルト値を与えるためのルーチン. 
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!== Future Plans
!
!  * 設定可能な擾乱のタイプを増やす
!  * エラー処理に gf4f90io を利用
!  * randomfile はこのサブルーチンで定義
!

module DistSet
  !
  !擾乱のデフォルト値を与えるためのルーチン. 
  !
  
  !モジュール読み込み
  use dc_message, only: MessageNotify
  use gridset,    only:DimXMin,         &! 配列の X 方向の下限
    &                  DimXMax,         &! 配列の X 方向の上限
    &                  RegXMin,         &! 配列の X 方向の下限
    &                  RegXMax,         &! 配列の X 方向の上限
    &                  DimZMin,         &! 配列の Z 方向の下限
    &                  DimZMax,         &! 配列の Z 方向の上限
    &                  RegZMin,         &! 配列の Z 方向の下限
    &                  RegZMax,         &! 配列の Z 方向の上限
    &                  MarginX,         &! 計算領域のマージン
    &                  MarginZ,         &! 計算領域のマージン
    &                  DelZ,            &! Z 方向の刻み幅
    &                  SpcNum,          &! 凝縮性化学種の数
    &                  XMin, XMax,      &! 
    &                  ZMin, ZMax,      &! 
    &                  x_X,             &! X 座標軸(スカラー格子点)
    &                  z_Z               ! Z 座標軸(スカラー格子点)
  use fileset,    only:RandomFile        ! 乱数ファイル
  use basicset,   only:                 &
    &                  SpcWetMolFr,     &!凝縮成分の初期モル比
    &                  MolWtWet,        &!凝縮成分の分子量
    &                  MolWtDry,        &!乾燥成分の分子量
    &                  xz_TempBasicZ,   &! 基本場の温度
    &                  xz_PressBasicZ,  &! 基本場の圧力
    &                  xz_ExnerBasicZ,  &! 無次元圧力
    &                  xza_MixRtBasicZ   ! 基本場の混合比
  use Boundary, only:  BoundaryXCyc_xza , &!
    &                  BoundaryZSym_xza,  & 
    &                  BoundaryXCyc_xz , &!
    &                  BoundaryZSym_xz , &
    &                  BoundaryXCyc_pz , &!
    &                  BoundaryZSym_pz 
  use ECCM,     only:  ECCM_Wet


  !暗黙の型宣言禁止
  implicit none

  !属性
  private

  !変数定義
  real(8), parameter, save :: Pi = 3.1415926535897932385d0 !円周率

  !public
  public distset_pottemp
  public distset_exner
  public distset_mixrt
  public distset_velx

contains
  
  subroutine distset_pottemp(cfgfile, xz_PotTemp)
    ! 暗黙の型宣言禁止
    implicit none

    !変数定義
    character(*), intent(in) :: cfgfile
    real(8), intent(out)     :: xz_PotTemp(DimXMin:DimXMax, DimZMin:DimZMax)
    character(20)    :: Type =""
    real(8)         :: DelMax = 0.0d0   !擾乱の最大値
    real(8)         :: XcRate = 0.0d0   !擾乱の中心位置(水平方向)の領域に対する割合
    real(8)         :: XrRate = 0.0d0   !擾乱の半径(水平方向)の領域に対する割合
    real(8)         :: YcRate = 0.0d0   !擾乱の中心位置(水平方向)の領域に対する割合
    real(8)         :: YrRate = 0.0d0   !擾乱の半径(水平方向)の領域に対する割合
    real(8)         :: ZcRate= 0.0d0    !擾乱の中心位置(鉛直方向)の領域に対する割合
    real(8)         :: ZrRate = 0.0d0   !擾乱の半径(鉛直方向)の領域に対する割合
    real(8)         :: Zpos = 0.0d0     !擾乱の Z 座標 [m] (Therma-Random 用)
    real(8)         :: Xc               !擾乱の中心位置(X方向)
    real(8)         :: Yc               !擾乱の中心位置(Y方向)
    real(8)         :: Zc               !擾乱の中心位置(鉛直方向)
    real(8)         :: Xr               !擾乱の半径(X方向)
    real(8)         :: Yr               !擾乱の半径(Y方向)
    real(8)         :: Zr               !擾乱の半径(鉛直方向)
    real(8)         :: Random           !ファイルから取得した乱数
    real(8)         :: beta(DimXMin:DimXMax, DimZMin:DimZMax)
                                        !擾乱の最大値に対する割合
    real(8)         :: RandomNum(DimXMin:DimXMax)
    real(8)         :: RandomNum2(DimXMin:DimXMax)
    integer         :: i, j, k, n

    ! NAMELIST ファイルの定義
    NAMELIST /disturbenv_pottemp/ Type, &
      & DelMax, XrRate, XcRate, YrRate, YcRate, ZrRate, ZcRate, Zpos

    ! namelist の読み込み
    open (10, FILE=cfgfile)
    read(10, NML=disturbenv_pottemp)
    close(10)

    ! 位置の決定
    !  XMax には MPI を含めた全領域が入っている
    Xr = minval( x_X, 1, x_X > (XMax - XMin) * XrRate )
    Xc = minval( x_X, 1, x_X > (XMax - XMin) * XcRate )
!    Yr = minval( y_Y, 1, y_Y > (YMax - YMin) * YrRate )
!    Yc = minval( y_Y, 1, y_Y > (YMax - YMin) * YcRate )
    Zr = minval( z_Z, 1, z_Z > (ZMax - ZMin) * ZrRate )
    Zc = minval( z_Z, 1, z_Z > (ZMax - ZMin) * ZcRate )
    
    ! 擾乱のタイプの設定
    select case(Type)
      
      ! 指定された高度にランダムな擾乱を与える
    case ("Random")      

      ! 乱数発生
      do i = RegXMin, RegXMax
        call random_number(random)
        RandomNum(i) = random
      end do
      
      ! 擾乱が全体としてはゼロとなるように調整. 平均からの差にする. 
      do i = RegXMin, RegXmax
        RandomNum2(i) = RandomNum(i) - sum( RandomNum(RegXMin+1:RegXMax) ) / real(RegXMax - RegXMin, 8) 
      end do
    
      ! 指定された高度に擾乱を置く. DelMax は温度なので, 温位に変換する. 
      do i = RegXMin, RegXMax
        xz_PotTemp(i, maxloc(z_Z, z_Z <= Zpos) - MarginZ - 1) = &
          & DelMax * RandomNum2(i) / xz_ExnerBasicZ(i, maxloc(z_Z, z_Z <= Zpos) - MarginZ - 1)
      end do
      
      ! 指定された座標を中心としたガウス分布の温位擾乱を与える. 
    case ("Gauss")
      do k = DimZMin, DimZMax
        do i = DimXMin, DimXMax
          xz_PotTemp(i,k) = &
            & DelMax * dexp( - ( (x_X(i) - Xc) / Xr )**2.0d0 * 5.0d-1   &
            &                - ( (z_Z(k) - Zc) / Zr )**2.0d0 * 5.0d-1 ) &
            & / xz_ExnerBasicZ(i,k)
        end do
      end do
    
      ! 指定された領域内に温位擾乱を与える (Klemp and Wilhelmson, 1978) 
    case ("KW1978")
      do k = DimZMin, DimZMax
        do i = DimXMin, DimXMax
          beta(i,k) =                                 &
            & (                                       &
            &      ( ( x_X(i) - Xc ) / Xr ) ** 2.0d0  &
            &    + ( ( z_Z(k) - Zc ) / Zr ) ** 2.0d0  &
            &  ) ** 5.0d-1
        end do
      end do
      
      where ( beta < 1.0d0 )
        xz_PotTemp = DelMax * ( dcos( Pi * 5.0d-1 * beta ) ** 2.0d0 ) &
          &          / xz_ExnerBasicZ
      end where
      
      ! Hueso and Sanchez-Lavega を模した設定
    case ("HS2001")
      i = ( DimXMax - DimXMin - 10) / 2 
      k = minloc( z_Z, 1, z_Z > 2.5d4 ) - MarginZ
      n = int( 5.0d3 / DelZ )
      
      xz_PotTemp(i-n:i,k-n:k) = DelMax
      
      ! Skamarock and Klemp (1989) の Cold-bubble 実験
    case ("SK1989")
      do k = DimZMin, DimZMax
        do i = DimXMin, DimXMax
          beta(i,k) =                                 &
            & sqrt(                                   &
            &      ( ( x_X(i) - Xc ) / Xr ) ** 2.0d0  &
            &    + ( ( z_Z(k) - Zc ) / Zr ) ** 2.0d0  &
            &  )
        end do
      end do
      
      where ( beta < 1.0d0 )
        xz_PotTemp = 0.5d0*DelMax*(cos(pi*beta) + 1.0d0)
      end where
    end select
    
    ! 境界条件
    call BoundaryXCyc_xz( xz_PotTemp )
    call BoundaryZSym_xz( xz_PotTemp )
    
  end subroutine distset_pottemp

  subroutine distset_exner(cfgfile, xz_Exner)

    ! 暗黙の型宣言禁止
    implicit none

    !変数定義
    character(*), intent(in) :: cfgfile
    real(8), intent(out)     :: xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax)
    character(20)   :: Type =""
    real(8)         :: DelMax = 0.0d0   !擾乱の最大値
    real(8)         :: XcRate = 0.0d0   !擾乱の中心位置(水平方向)の領域に対する割合
    real(8)         :: XrRate = 0.0d0   !擾乱の半径(水平方向)の領域に対する割合
    real(8)         :: YcRate = 0.0d0   !擾乱の中心位置(水平方向)の領域に対する割合
    real(8)         :: YrRate = 0.0d0   !擾乱の半径(水平方向)の領域に対する割合
    real(8)         :: ZcRate= 0.0d0    !擾乱の中心位置(鉛直方向)の領域に対する割合
    real(8)         :: ZrRate = 0.0d0   !擾乱の半径(鉛直方向)の領域に対する割合
    real(8)         :: Xc               !擾乱の中心位置(X方向)
    real(8)         :: Yc               !擾乱の中心位置(Y方向)
    real(8)         :: Zc               !擾乱の中心位置(鉛直方向)
    real(8)         :: Xr               !擾乱の半径(X方向)
    real(8)         :: Yr               !擾乱の半径(Y方向)
    real(8)         :: Zr               !擾乱の半径(鉛直方向)
    integer         :: i, j, k

    ! NAMELIST ファイルの定義
    NAMELIST /disturbenv_exner/ Type, &
      & DelMax, XrRate, XcRate, YrRate, YcRate, ZrRate, ZcRate

    ! namelist の読み込み    
    open (10, FILE=cfgfile)
    read(10, NML=disturbenv_exner)
    close(10)

    ! 位置の決定
    !  XMax には MPI を含めた全領域が入っている
    Xr = minval( x_X, 1, x_X > (XMax - XMin) * XrRate )
    Xc = minval( x_X, 1, x_X > (XMax - XMin) * XcRate )
!    Yr = minval( y_Y, 1, y_Y > (YMax - YMin) * YrRate )
!    Yc = minval( y_Y, 1, y_Y > (YMax - YMin) * YcRate )
    Zr = minval( z_Z, 1, z_Z > (ZMax - ZMin) * ZrRate )
    Zc = minval( z_Z, 1, z_Z > (ZMax - ZMin) * ZcRate )

    ! エクスナー関数の設定. 
    select case (Type)    

      ! 指定された場所に, ガウシアンなエクスナー関数の擾乱を与える. 
    case ("Gauss")
      do k = DimZMin, DimZMax
        do i = DimXMin, DimXMax
          xz_Exner(i,k) =                                             &
            & DelMax * dexp( - ( (x_X(i) - Xc) / Xr )**2.0d0 * 5.0d-1 &
            &                - ( (z_Z(k) - Zc) / Zr )**2.0d0 * 5.0d-1 ) 
        end do
      end do
    end select
    
    ! 境界条件
    call BoundaryXCyc_xz( xz_Exner )
    call BoundaryZSym_xz( xz_Exner )    
    
  end subroutine distset_exner


  subroutine distset_mixrt(cfgfile, xza_MixRt)
    ! 暗黙の型宣言禁止
    implicit none

    !変数定義
    character(*), intent(in) :: cfgfile
    real(8), intent(out)     :: xza_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, 1:SpcNum)
    character(20)   :: Type =""
    real(8)         :: DelMax = 0.0d0   !擾乱の最大値
    real(8)         :: XcRate = 0.0d0   !擾乱の中心位置(水平方向)の領域に対する割合
    real(8)         :: XrRate = 0.0d0   !擾乱の半径(水平方向)の領域に対する割合
    real(8)         :: YcRate = 0.0d0   !擾乱の中心位置(水平方向)の領域に対する割合
    real(8)         :: YrRate = 0.0d0   !擾乱の半径(水平方向)の領域に対する割合
    real(8)         :: ZcRate= 0.0d0    !擾乱の中心位置(鉛直方向)の領域に対する割合
    real(8)         :: ZrRate = 0.0d0   !擾乱の半径(鉛直方向)の領域に対する割合
    real(8)         :: XposMin = 0.0d0  !
    real(8)         :: XposMax = 0.0d0  !
    real(8)         :: ZposMin = 0.0d0  !
    real(8)         :: ZposMax = 0.0d0  !
    real(8)         :: Humidity = 0.0d0 !相対湿度
    real(8)         :: Xc               !擾乱の中心位置(X方向)
    real(8)         :: Yc               !擾乱の中心位置(Y方向)
    real(8)         :: Zc               !擾乱の中心位置(鉛直方向)
    real(8)         :: Xr               !擾乱の半径(X方向)
    real(8)         :: Yr               !擾乱の半径(Y方向)
    real(8)         :: Zr               !擾乱の半径(鉛直方向)
    integer         :: i, j, k, s
    real(8)         :: MolFrIni(SpcNum)
    real(8)         :: za_MolFr(DimZMin:DimZMax, SpcNum)

    ! NAMELIST ファイルの定義
    NAMELIST /disturbenv_mixrt/ Type, &
      & DelMax, XrRate, XcRate, YrRate, YcRate, ZrRate, ZcRate, &
      & XposMin, XposMax, ZposMin, ZposMax, Humidity

    ! namelist の読み込み    
    open (10, FILE=cfgfile)
    read(10, NML=disturbenv_mixrt)
    close(10)
    
    ! 位置の決定
    !  XMax には MPI を含めた全領域が入っている
    Xr = minval( x_X, 1, x_X > (XMax - XMin) * XrRate )
    Xc = minval( x_X, 1, x_X > (XMax - XMin) * XcRate )
!    Yr = minval( y_Y, 1, y_Y > (YMax - YMin) * YrRate )
!    Yc = minval( y_Y, 1, y_Y > (YMax - YMin) * YcRate )
    Zr = minval( z_Z, 1, z_Z > (ZMax - ZMin) * ZrRate )
    Zc = minval( z_Z, 1, z_Z > (ZMax - ZMin) * ZcRate )

    ! 擾乱の設定
    select case (Type)
    
      ! 指定された座標を中心としたガウス分布の混合比擾乱を与える. 
    case ("Gauss")
      do k = DimZMin, DimZMax
        do i = DimXMin, DimXMax
          do s = 1, SpcNum
            xza_MixRt(i,:,s) =                                     &
              & DelMax                                             &
              & * dexp( - ( (x_X(i) - Xc) / Xr )**2.0d0 * 5.0d-1   &
              &         - ( (z_Z(k) - Zc) / Zr )**2.0d0 * 5.0d-1 )
          end do
        end do
      end do
    
      ! XposMin:XposMax,ZposMin:ZposMax で囲まれた領域の初期の湿度をゼロにするために
      ! 基本場と逆符号の水蒸気擾乱を与える
    case ("Dry")
      do s = 1, SpcNum
        do k = DimZMin,DimZMax  
          do i = DimXMin,DimXMax
            if (z_Z(k) >= ZposMin .AND. z_Z(k) < ZposMax &
              & .AND. x_X(i) >= XposMin .AND. x_X(i) < XposMax) then
              xza_MixRt(i,k,s) = - xza_MixRtBasicZ(i,k,s)
            end if
          end do
        end do
      end do
      
      ! 与えられた相対湿度に設定する.
    case ("Humidity")
      MolFrIni = SpcWetMolFr(1:SpcNum) 

      ! 水平一様なので, i=0 だけ計算. 
      call ECCM_Wet( MolFrIni, Humidity, xz_TempBasicZ(1,:), xz_PressBasicZ(1,:),  za_MolFr )
      
      !気相のモル比を混合比に変換
      do s = 1, SpcNum
        do k = RegZMin, RegZMax
          do i = RegXMin, RegXMax
            xza_MixRt(i,k,s) = za_MolFr(k,s) * MolWtWet(s) / MolWtDry - xza_MixRtBasicZ(i,k,s)
          end do
        end do
      end do

    end select

    ! 境界条件
    call BoundaryXCyc_xza( xza_MixRt )
    call BoundaryZSym_xza( xza_MixRt )

  end subroutine distset_mixrt
  

  subroutine distset_VelX(cfgfile, pz_VelX)
    ! 暗黙の型宣言禁止
    implicit none

    !変数定義
    character(*), intent(in) :: cfgfile
    real(8), intent(out)     :: pz_VelX(DimXMin:DimXMax, DimZMin:DimZMax)
    character(20)   :: Type =""
    real(8)         :: Umax = 0.0d0     !擾乱の中心位置(水平方向)の領域に対する割合
    real(8)         :: ZposMin = 0.0d0  !擾乱の Z 座標 [m] (下端)
    real(8)         :: ZposMax = 0.0d0  !擾乱の Z 座標 [m] (上端)
    integer         :: i, j, k, b, u

    ! namelist の定義
    NAMELIST /disturbenv_VelX/ Type, Umax, ZposMin, ZposMax

    ! 値を読み出す
    open (10, FILE=cfgfile)
    read(10, NML=disturbenv_velx)
    close(10)
    
    select case(Type)
    
    case ("Takemi2007")
      !
      != 概要
      !* case "Takemi2007" での計算時に鉛直シアーのある風を与える時に使用する
      !* 風の与え方には, 以下のようなバリエーションがある
      !  (1) シアーを与える高度を変える
      !  (2) シアーのある風の最大風速 (U_s) を変える
      !
      !  (1) については, (a) 0 - 2.5 km, (b) 2.5 - 5.0 km, (c) 5.0 - 7.5 km の
      !  三パターンがある
      !  (2) については, Takemi (2007) では熱帯場と中緯度場の温度場毎に
      !  異なる値を設定している
      !
      !  その強度(Us)は, 以下の通り
      !  <熱帯場>   (1) 5 m/s, (2) 10 m/s, (3) 15 m/s
      !  <中緯度場> (1) 10 m/s, (2) 15 m/s, (3) 20 m/s
      !
      !* シアーの形の模式図 (Takemi, 2007)   |
      !                                     /| 7.5 km
      !                                    / |
      !                                   /  |
      !                                  / ←|
      !                                 ｜  /| 5.0 km
      !                                 ｜ / |
      !                                 ｜/  |
      !                                  / ←|
      !                                 ｜  /| 2.5 km
      !                                 ｜ / |
      !                                 ｜/  |
      !                                  / ←|
      !---------------------------------+------------- 0.0 km
      !                                Us (m/s)
      !
      
      ! シアーを与える高度を決める
      ! 傾き始めの高度の格子点 (B) と 傾きの終わりの格子点 (U) で決定
      do k = RegZMin+1, DimZMax
        if (z_Z(k) <= ZposMin .AND. ZposMin < z_Z(k+1)) then
          B = k
        end if
        
        if (z_Z(k) <= ZposMax .AND. ZposMax < z_Z(k+1)) then      
          U = k
        end if
      end do
      
      ! 上記で決めた格子点BからUの間に鉛直シアーのある風を与える
      do k = RegZMin+1, U
        do i = DimXMin, DimXMax
          if (k < B) then
            pz_VelX(i,k) =  Umax
          else
            pz_VelX(i,k) =  Umax - (Umax / (z_Z(U) - z_Z(B)) * (z_Z(k) - z_Z(B)))
          end if
        end do
      end do
    end select

    !境界条件
    call BoundaryXCyc_pz( pz_VelX )
    call BoundaryZSym_pz( pz_VelX )

  end subroutine Distset_VelX

end module DistSet
