!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004, 2005. All rights reserved.
!---------------------------------------------------------------------
!= Module DisturbEnv
!
!   * Developer: SUGIYAMA Ko-ichiro, ODAKA Masatsugu
!   * Version: $Id: initialdata_disturb.f90,v 1.1 2011-04-26 04:08:10 sugiyama Exp $ 
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview 
!
! 擾乱のデフォルト値を与えるための基本関数群. 
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!== Future Plans
!
!

module initialdata_disturb
  !
  !擾乱のデフォルト値を与えるためのルーチン. 
  !
  
  !モジュール読み込み
  use dc_types,   only: STRING, DP
  use dc_message, only: MessageNotify
  use gridset_3d, only: XMin, XMax, YMin, YMax, ZMin, ZMax, &
    &                  imin,         &! 配列の X 方向の下限
    &                  imax,         &! 配列の X 方向の上限
    &                  jmin,         &! 配列の Y 方向の下限
    &                  jmax,         &! 配列の Y 方向の上限
    &                  kmin,         &! 配列の Z 方向の下限
    &                  kmax,         &! 配列の Z 方向の上限
    &                  nx,           &! 配列の Z 方向の下限
    &                  ny,           &! 配列の Z 方向の上限
    &                  nz,           &! 計算領域のマージン
    &                  na,           &! 計算領域のマージン
    &                  x_X,             &! X 座標軸(スカラー格子点)
    &                  y_Y,             &! X 座標軸(スカラー格子点)
    &                  z_Z               ! Z 座標軸(スカラー格子点)
  use basicset_3d,   only:                 &
    &                  SpcWetMolFr,     &!凝縮成分の初期モル比
    &                  MolWtWet,        &!凝縮成分の分子量
    &                  MolWtDry,        &!乾燥成分の分子量
    &                  xyz_TempBasicZ,   &! 基本場の温度
    &                  xyz_PressBasicZ,  &! 基本場の圧力
    &                  xyza_MixRtBasicZ   ! 基本場の混合比
  use eccm_3d,  only:  eccm_molfr

  !暗黙の型宣言禁止
  implicit none

  !属性
  private

  !変数定義
  ! Moist 用
  real(DP), save   :: Humidity = 0.0d0 !相対湿度
  ! Thermal-Gauss 用
  real(DP), save   :: DelMax = 0.0d0   !擾乱の最大値
  real(DP), save   :: XrRate = 0.0d0   !擾乱の中心位置(水平方向)の領域に対する割合
  real(DP), save   :: YrRate = 0.0d0   !擾乱の半径(水平方向)の領域に対する割合
  real(DP), save   :: ZrRate = 0.0d0   !擾乱の半径(鉛直方向)の領域に対する割合
  real(DP), save   :: XcRate = 0.0d0   !擾乱の中心位置(水平方向)の領域に対する割合
  real(DP), save   :: YcRate = 0.0d0   !擾乱の半径(水平方向)の領域に対する割合
  real(DP), save   :: ZcRate= 0.0d0    !擾乱の中心位置(鉛直方向)の領域に対する割合
  real(DP), save   :: Xc = 0.0d0       !擾乱の中心位置(X方向)
  real(DP), save   :: Yc = 0.0d0       !擾乱の中心位置(Y方向)
  real(DP), save   :: Zc = 0.0d0       !擾乱の中心位置(鉛直方向)
  real(DP), save   :: Xr = 0.0d0       !擾乱の半径(X方向)
  real(DP), save   :: Yr = 0.0d0       !擾乱の半径(Y方向)
  real(DP), save   :: Zr = 0.0d0       !擾乱の半径(鉛直方向)
  ! Therma-Random 用
  real(DP), save   :: Zpos = 0.0d0     !擾乱の Z 座標 [m] 
  ! DryRegion 用
  real(DP), save   :: XposMin = 0.0d0    !乾燥域の X 座標 [m] 
  real(DP), save   :: YposMin = 0.0d0    !乾燥域の Y 座標 [m] 
  real(DP), save   :: ZposMin = 0.0d0    !乾燥域の Z 座標 [m] 
  real(DP), save   :: XposMax = 0.0d0    !乾燥域の X 座標 [m] 
  real(DP), save   :: YposMax = 0.0d0    !乾燥域の Y 座標 [m] 
  real(DP), save   :: ZposMax = 0.0d0    !乾燥域の Z 座標 [m] 
  ! FLAGS
  integer, save, public :: ID_FlagInitialdataDisturbPotTemp = 0
  integer, save, public :: ID_FlagInitialdataDisturbMixRt = 0
  integer, parameter, public :: ID_FlagInitialDataDisturbPotTemp_GaussXY = 1
  integer, parameter, public :: ID_FlagInitialDataDisturbPotTemp_GaussXZ = 2
  integer, parameter, public :: ID_FlagInitialDataDisturbPotTemp_GaussXYZ = 3
  integer, parameter, public :: ID_FlagInitialDataDisturbPotTemp_random = 4
  integer, parameter, public :: ID_FlagInitialDataDisturbMixRt_dryregion = 1
  integer, parameter, public :: ID_FlagInitialDataDisturbMixRt_moist = 2

  public initialdata_disturb_nml
  public initialdata_disturb_random
  public initialdata_disturb_gaussXZ
  public initialdata_disturb_gaussXY
  public initialdata_disturb_gaussXYZ
  public initialdata_disturb_dryregion
  public initialdata_disturb_moist

contains

  subroutine initialdata_disturb_nml(cfgfile)
    ! nml の読み込み
    !
    use dc_iounit,     only : FileOpen    

    implicit none

    character(STRING), intent(in) :: cfgfile    
    integer                       :: unit     !装置番号
    character(STRING)  :: FlagInitialdataDisturbPotTemp = ""
    character(STRING)  :: FlagInitialdataDisturbMixRt = ""

    NAMELIST /initialdata_disturb/ &
      & FlagInitialdataDisturbPotTemp, FlagInitialdataDisturbMixRt, &
      & DelMax, XrRate, XcRate, YrRate, YcRate, ZrRate, ZcRate, Zpos, &
      & XposMin, YposMin, ZposMin, XposMax, YposMax, ZposMax, Humidity
    
    call FileOpen(unit, file=cfgfile, mode='r')
    read(unit, NML=initialdata_disturb)
    close(unit)

    if (FlagInitialdataDisturbPotTemp == "GaussXY") then 
      ID_FlagInitialdataDisturbPotTemp = ID_FlagInitialdataDisturbPotTemp_GaussXY
    elseif (FlagInitialdataDisturbPotTemp == "GaussXZ") then 
      ID_FlagInitialdataDisturbPotTemp = ID_FlagInitialdataDisturbPotTemp_GaussXZ
    elseif (FlagInitialdataDisturbPotTemp == "GaussXYZ") then 
      ID_FlagInitialdataDisturbPotTemp = ID_FlagInitialdataDisturbPotTemp_GaussXYZ
    elseif (FlagInitialdataDisturbPotTemp == "random") then 
      ID_FlagInitialdataDisturbPotTemp = ID_FlagInitialdataDisturbPotTemp_random
    end if

    if (FlagInitialdataDisturbMixRt == "dryregion") then 
      ID_FlagInitialdataDisturbMixRt = ID_FlagInitialdataDisturbMixRt_dryregion
    elseif (FlagInitialdataDisturbMixRt == "moist") then 
      ID_FlagInitialdataDisturbMixRt = ID_FlagInitialdataDisturbMixRt_moist
    end if
    
    ! 位置の決定
    !  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 )

  end subroutine initialdata_disturb_nml

  
  subroutine initialdata_disturb_random( xyz_Var )
    
    implicit none
    
    real(DP), intent(out) :: xyz_Var(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)              :: Random           !ファイルから取得した乱数
    real(DP)              :: Random1(imin:imax, jmin:jmax)
    real(DP)              :: Random2(imin:imax, jmin:jmax)
    integer :: i, j, k, kpos

    !  0.0--1.0 の擬似乱数発生
    do j = jmin, jmax
      do i = imin, imax
        call random_number(random)
        Random1(i,j) = random
      end do
    end do

    ! 擾乱が全体としてはゼロとなるように調整. 平均からの差にする. 
    do j = jmin, jmax
      do i = imin, imax
        Random2(i,j) = Random1(i,j) - sum( Random1(1:nx,j) ) / real(nx, 8) 
      end do
    end do

    ! 指定された高度の配列添字を用意
    do k = kmin, kmax
      if ( z_Z(k) >= Zpos ) then 
        kpos = k
        exit
      end if
    end do

    ! 指定された高度に擾乱を置く. 
    do j = 1, ny
      do i = 1, nx
        xyz_Var(i, j, kpos) = DelMax * Random2(i,j) / sum( Random2(1:nx,j) )
      end do
    end do

  end subroutine initialdata_disturb_random
  
  
  subroutine initialdata_disturb_gaussXZ(xyz_Var)

    implicit none
    
    real(DP), intent(out) :: xyz_Var(imin:imax, jmin:jmax, kmin:kmax)
    integer         :: i, j, k, n

    do k = kmin, kmax
      do j = jmin, jmax
        do i = imin, imax
          xyz_Var(i,j,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 do
    
  end subroutine initialdata_disturb_gaussXZ
  

  subroutine initialdata_disturb_gaussXY(xyz_Var)
    
    implicit none
    
    real(DP), intent(out) :: xyz_Var(imin:imax, jmin:jmax, kmin:kmax)
    integer         :: i, j, k, n

    do k = kmin, kmax
      do j = jmin, jmax
        do i = imin, imax
          xyz_Var(i,j,k) = &
            & DelMax * dexp( - ( (x_X(i) - Xc) / Xr )**2.0d0 * 5.0d-1   &
            &                - ( (y_Y(k) - Yc) / Yr )**2.0d0 * 5.0d-1 )
        end do
      end do
    end do
    
  end subroutine initialdata_disturb_gaussXY


  subroutine initialdata_disturb_gaussXYZ(xyz_Var)
    
    implicit none
    
    real(DP), intent(out) :: xyz_Var(imin:imax, jmin:jmax, kmin:kmax)
    integer         :: i, j, k, n

    do k = kmin, kmax
      do j = jmin, jmax
        do i = imin, imax
          xyz_Var(i,j,k) = &
            & DelMax * dexp( - ( (x_X(i) - Xc) / Xr )**2.0d0 * 5.0d-1   &
            &                - ( (y_Y(k) - Yc) / Yr )**2.0d0 * 5.0d-1   &
            &                - ( (z_Z(k) - Zc) / Zr )**2.0d0 * 5.0d-1 ) 
        end do
      end do
    end do
    
  end subroutine initialdata_disturb_gaussXYZ


  subroutine initialdata_disturb_dryregion(xyza_MixRt)
    
    implicit none
    
    real(DP), intent(out) :: xyza_MixRt(imin:imax, jmin:jmax, kmin:kmax, 1:na)
    integer         :: i, j, k, s
    
    ! XposMin:XposMax,ZposMin:ZposMax で囲まれた領域の初期の湿度をゼロにするために
    ! 基本場と逆符号の水蒸気擾乱を与える
    do s = 1, na
      do k = kmin,kmax  
        do j = jmin, jmax
          do i = imin,imax
            if (z_Z(k) >= ZposMin .AND. z_Z(k) < ZposMax &
              & .AND. y_Y(j) >= YposMin .AND. y_Y(k) < YposMax &
              & .AND. x_X(i) >= XposMin .AND. x_X(i) < XposMax) then
              xyza_MixRt(i,j,k,s) = - xyza_MixRtBasicZ(i,j,k,s)
            end if
          end do
        end do
      end do
    end do
    
  end subroutine initialdata_disturb_dryregion
  
  
  subroutine initialdata_disturb_moist(xyza_MixRt)
    
    implicit none
    
    real(DP), intent(out) :: xyza_MixRt(imin:imax, jmin:jmax, kmin:kmax, 1:na)
    real(DP)        :: za_MolFr(kmin:kmax, 1:na)
    integer         :: i, j, k, s
  
    ! 湿度ゼロなら何もしない
    if ( Humidity == 0.0d0 ) return

    ! 水平一様なので, i=0 だけ計算. 
    i = 1
    j = 1
    call eccm_molfr( SpcWetMolFr(1:na), Humidity, xyz_TempBasicZ(i,j,:), &
      &              xyz_PressBasicZ(i,j,:), za_MolFr )
    
    !気相のモル比を混合比に変換
    do s = 1, na
      do k = 1, nz
        do j = 1, ny
          do i = 1, nx
            xyza_MixRt(i,j,k,s) = za_MolFr(k,s) * MolWtWet(s) / MolWtDry - xyza_MixRtBasicZ(i,j,k,s)
          end do
        end do
      end do
    end do
    
  end subroutine initialdata_disturb_moist
  
end module initialdata_disturb
