!---------------------------------------------------------------------
!     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.12 2012-01-12 11:53:08 sugiyama Exp $ 
!   * Tag Name: $Name: arare5-20120229 $
!   * 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 mpi_wrapper,only: myrank, nprocs
  use axesset,   only: &
    &                  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             ! 計算領域のマージン

  !暗黙の型宣言禁止
  implicit none

  !属性
  private

  public initialdata_disturb_random
  public initialdata_disturb_gaussXZ
  public initialdata_disturb_gaussXY
  public initialdata_disturb_gaussYZ
  public initialdata_disturb_gaussXYZ
  public initialdata_disturb_dryreg
  public initialdata_disturb_rectangle
  public initialdata_disturb_moist

contains
    
  subroutine initialdata_disturb_random( DelMax, Zpos, xyz_Var )
    
    implicit none
    
    real(DP), intent(in)  :: DelMax, Zpos
    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, ix, jy

    ! 初期化
    xyz_Var = 0.0d0

    ! 0.0--1.0 の擬似乱数発生
    !  mpi の場合に, 各 CPU の持つ乱数が異なるよう調整している.  
    !
    do j = jmin, jmax + ( ny * nprocs )
      do i = imin, imax * ( nx * nprocs )
        call random_number(random)
        if (imin + nx * myrank <= i .AND. i <= imax + nx * myrank) then 
          if (jmin + ny * myrank <= j .AND. j <= jmax + ny * myrank) then 
            ix = i - nx * myrank
            jy = j - ny * myrank
            Random1(ix,jy) = random
          end if
        end if
      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(DelMax, Xc, Xr, Zc, Zr, xyz_Var)

    implicit none

    real(DP), intent(in)  :: DelMax, Xc, Xr, Zc, Zr
    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(DelMax, Xc, Xr, Yc, Yr, xyz_Var)
    
    implicit none

    real(DP), intent(in)  :: DelMax, Xc, Xr, Yc, Yr
    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_gaussYZ(DelMax, Yc, Yr, Zc, Zr, xyz_Var)
    
    implicit none

    real(DP), intent(in)  :: DelMax, Yc, Yr, Zc, Zr
    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( - ( (z_Z(k) - Zc) / Zr )**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_gaussYZ


  subroutine initialdata_disturb_gaussXYZ(DelMax, Xc, Xr, Yc, Yr, Zc, Zr, xyz_Var)
    
    implicit none

    real(DP), intent(in)  :: DelMax, Xc, Xr, Yc, Yr, Zc, Zr
    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_dryreg( &
    & XposMin, XposMax, YposMin, YposMax, ZposMin, ZposMax, &
    & xyzf_QMix)

    use basicset, only: xyzf_QMixBZ
    
    implicit none

    real(DP), intent(in)  ::XposMin, XposMax, YposMin, YposMax, ZposMin, ZposMax
    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_dryreg
  
  
  subroutine initialdata_disturb_rectangle( &
    & DelMax, XposMin, XposMax, YposMin, YposMax, ZposMin, ZposMax, &
    & xyz_Var)

    implicit none

    real(DP), intent(in)  :: DelMax, XposMin, XposMax, YposMin, YposMax, ZposMin, ZposMax
    real(DP), intent(out) :: xyz_Var(imin:imax, jmin:jmax, kmin:kmax)
    integer         :: i, j, k

    write(*,*) DelMax, XposMin, XposMax, YposMin, YposMax, ZposMin, ZposMax
    
    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
            xyz_Var(i,j,k) = DelMax
            write(*,*) x_X(i), y_Y(j), z_Z(k), xyz_Var(i,j,k)
          end if
        end do
      end do
    end do
    
  end subroutine initialdata_disturb_rectangle
  
  
  subroutine initialdata_disturb_moist(Hum, xyzf_QMix)
    
    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

    real(DP), intent(in)  :: Hum
    real(DP), intent(out) :: xyzf_QMix(imin:imax, jmin:jmax, kmin:kmax, 1:ncmax)
    real(DP)              :: zf_MolFr(kmin:kmax, 1:ncmax)
    integer               :: i, j, k, s
  
    ! 湿度ゼロなら何もしない
    if ( Hum == 0.0d0 ) return

    ! 水平一様なので, i=0 だけ計算. 
    i = 1
    j = 1
    call eccm_molfr( SpcWetMolFr(1:ncmax), Hum, xyz_TempBZ(i,j,:), &
      &              xyz_PressBZ(i,j,:), zf_MolFr )
    
    !気相のモル比を混合比に変換
    do s = 1, ncmax
      do k = 1, nz
        do j = 1, ny
          do i = 1, nx
            xyzf_QMix(i,j,k,s) = zf_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
