!---------------------------------------------------------------------
!     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.7 2011-06-03 16:31:45 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 axesset,   only: XMin, XMax, YMin, YMax, ZMin, ZMax, &
    &                  x_X,             &! X 座標軸(スカラー格子点)
    &                  y_Y,             &! X 座標軸(スカラー格子点)
    &                  z_Z               ! Z 座標軸(スカラー格子点)
  use gridset,   only: &
    &                  imin,         &! 配列の X 方向の下限
    &                  imax,         &! 配列の X 方向の上限
    &                  jmin,         &! 配列の Y 方向の下限
    &                  jmax,         &! 配列の Y 方向の上限
    &                  kmin,         &! 配列の Z 方向の下限
    &                  kmax,         &! 配列の Z 方向の上限
    &                  nx,           &! 配列の Z 方向の下限
    &                  ny,           &! 配列の Z 方向の上限
    &                  nz,           &! 計算領域のマージン
    &                  ncmax             ! 計算領域のマージン
  use basicset,   only:                 &
    &                  xyz_TempBZ,   &! 基本場の温度
    &                  xyz_PressBZ,  &! 基本場の圧力
    &                  xyzf_QMixBZ   ! 基本場の混合比
  use composition,   only:                 &
    &                  MolWtWet,        &!凝縮成分の分子量
    &                  SpcWetMolFr       !凝縮成分の初期モル比
  use constants, only: MolWtDry          !乾燥成分の分子量
  use eccm,       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_FlaginitdataPTemp = 0
  integer, save, public :: ID_FlaginitdataExner = 0
  integer, save, public :: ID_FlaginitdataQMix = 0
  integer, parameter, public :: ID_Flaginitdata_GaussXY = 1
  integer, parameter, public :: ID_Flaginitdata_GaussXZ = 2
  integer, parameter, public :: ID_Flaginitdata_GaussXYZ = 3
  integer, parameter, public :: ID_Flaginitdata_random = 4
  integer, parameter, public :: ID_FlaginitdataQMix_dryregion = 5
  integer, parameter, public :: ID_FlaginitdataQMix_moist = 6

  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)  :: FlaginitdataPTemp = ""
    character(STRING)  :: FlaginitdataExner = ""
    character(STRING)  :: FlaginitdataQMix = ""

    NAMELIST /initialdata_disturb/ &
      & FlaginitdataPTemp, FlaginitdataExner, FlaginitdataQMix, &
      & 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 (FlaginitdataPTemp == "GaussXY") then 
      ID_FlaginitdataPTemp = ID_Flaginitdata_GaussXY
    elseif (FlaginitdataPTemp == "GaussXZ") then 
      ID_FlaginitdataPTemp = ID_Flaginitdata_GaussXZ
    elseif (FlaginitdataPTemp == "GaussXYZ") then 
      ID_FlaginitdataPTemp = ID_Flaginitdata_GaussXYZ
    elseif (FlaginitdataPTemp == "random") then 
      ID_FlaginitdataPTemp = ID_Flaginitdata_random
    end if

    if (FlaginitdataExner == "GaussXY") then 
      ID_FlaginitdataExner = ID_Flaginitdata_GaussXY
    elseif (FlaginitdataExner == "GaussXZ") then 
      ID_FlaginitdataExner = ID_Flaginitdata_GaussXZ
    elseif (FlaginitdataExner == "GaussXYZ") then 
      ID_FlaginitdataExner = ID_Flaginitdata_GaussXYZ
    elseif (FlaginitdataExner == "random") then 
      ID_FlaginitdataExner = ID_Flaginitdata_random
    end if

    if (FlaginitdataQMix == "GaussXY") then 
      ID_FlaginitdataQMix = ID_Flaginitdata_GaussXY
    elseif (FlaginitdataQMix == "GaussXZ") then 
      ID_FlaginitdataQMix = ID_Flaginitdata_GaussXZ
    elseif (FlaginitdataQMix == "GaussXYZ") then 
      ID_FlaginitdataQMix = ID_Flaginitdata_GaussXZ
    elseif (FlaginitdataQMix == "dryregion") then 
      ID_FlaginitdataQMix = ID_FlaginitdataQMix_dryregion
    elseif (FlaginitdataQMix == "moist") then 
      ID_FlaginitdataQMix = ID_FlaginitdataQMix_moist
    end if
    
    ! 位置の決定
    !  XMax には MPI を含めた全領域が入っている
    Xc = (XMax - XMin) * XcRate 
    Yc = (YMax - YMin) * YcRate 
    Zc = (ZMax - ZMin) * ZcRate 
    Xr = (XMax - XMin) * XrRate 
    Yr = (YMax - YMin) * YrRate 
    Zr = (ZMax - ZMin) * ZrRate 
!    Xr = minval( x_X, 1, x_X > (XMax - XMin) * XrRate )
!    Yr = minval( y_Y, 1, y_Y > (YMax - YMin) * YrRate )
!    Zr = minval( z_Z, 1, z_Z > (ZMax - ZMin) * ZrRate )
!    Xc = minval( x_X, 1, x_X > (XMax - XMin) * XcRate )
!    Yc = minval( y_Y, 1, y_Y > (YMax - YMin) * YcRate )
!    Zc = minval( z_Z, 1, z_Z > (ZMax - ZMin) * ZcRate )

!    call MessageNotify( "M", "initialdata_disturb", "Xr = %f", d=(/Xr/) )
!    call MessageNotify( "M", "initialdata_disturb", "Xc = %f", d=(/Xc/) )
!    call MessageNotify( "M", "initialdata_disturb", "Yr = %f", d=(/Yr/) )
!    call MessageNotify( "M", "initialdata_disturb", "Yc = %f", d=(/Yc/) )
!    call MessageNotify( "M", "initialdata_disturb", "Zr = %f", d=(/Zr/) )
!    call MessageNotify( "M", "initialdata_disturb", "Zc = %f", d=(/Zc/) )
    
  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)
    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 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 * (Random1(i,j) - sum( Random1(1:nx,1:ny) ) / real((nx * ny),8))
      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

    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

!    where ( xyz_Var < DelMax * 1.0d-2) 
!      xyz_Var = 0.0d0
!    end where
    
  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
    
    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(j) - Yc) / Yr )**2.0d0 * 5.0d-1 )
        end do
      end do
    end do

!    where ( xyz_Var < DelMax * 1.0d-2) 
!      xyz_Var = 0.0d0
!    end where
    
  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
    
    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(j) - Yc) / Yr )**2.0d0 * 5.0d-1   &
            &                - ( (z_Z(k) - Zc) / Zr )**2.0d0 * 5.0d-1 ) 
        end do
      end do
    end do
    
!    where ( xyz_Var < DelMax * 1.0d-2) 
!      xyz_Var = 0.0d0
!    end where

  end subroutine initialdata_disturb_gaussXYZ


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

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