!---------------------------------------------------------------------
!     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.17 2014/07/08 00:59:08 sugiyama Exp $ 
!   * Tag Name: $Name:  $
!   * Change History: 


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, xyz_X, XMax,    &! X 座標軸(スカラー格子点)
    &                  y_Y, xyz_Y, YMax,    &! X 座標軸(スカラー格子点)
    &                  z_Z, xyz_Z, ZMax      ! 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_cosXZ
  public initialdata_disturb_cosXY
  public initialdata_disturb_cosYZ
  public initialdata_disturb_cosXYZ
  public initialdata_disturb_coneXZ
  public initialdata_disturb_coneXY
  public initialdata_disturb_coneYZ
  public initialdata_disturb_dryreg
  public initialdata_disturb_square
  public initialdata_disturb_moist
  public initialdata_disturb_tanh
  public initialdata_disturb_tanh_sin
  public initialdata_disturb_circleXZ

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
    kpos    = 1

    ! 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

!!!
!!! Klemp and Wilhelmson (1978) の初期値
!!! 

!  subroutine initialdata_disturb_kw1978(DelMax, Xc, Xr, Yc, Yr, Zc, Zr, xyz_Var)
!    !
!    ! Klemp and Wilhelmson (1978) の初期値
!
!    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)
!    real(DP)               :: xyz_beta(imin:imax, jmin:jmax, kmin:kmax)
!    integer                :: i, j, k
!
!    do k = kmin, kmax
!      do j = jmin, jmax
!        do i = imin, imax
!          beta(i,j,k) =                               &
!            & (                                       &
!            &      ( ( x_X(i) - Xc ) / Xr ) ** 2.0d0  &
!            &    + ( ( y_Y(j) - Yc ) / Yr ) ** 2.0d0  &
!            &    + ( ( z_Z(k) - Zc ) / Zr ) ** 2.0d0  &
!            &  ) ** 5.0d-1
!        end do
!      end do
!    end do
!    
!    where ( beta < 1.0d0 )
!      xyz_Var = DelMax * ( dcos( Pi * 5.0d-1 * beta ) ** 2.0d0 )
!    end where
!
!  end subroutine initialdata_disturb_kw1978



!!!
!!! tanh
!!! 

  subroutine initialdata_disturb_tanh(VarMean, VarDel, Zc, Zr, aaz_Var, aaz_VarBZ)
    !
    ! tanh 型のシア
    !   A(z) = A0 + A1 \tanh( (z - Zc) / h )

    implicit none

    real(DP), intent(in)            :: VarMean, VarDel, Zc, Zr
    real(DP), intent(in), optional  :: aaz_VarBZ(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), intent(out)           :: aaz_Var(imin:imax, jmin:jmax, kmin:kmax)
    integer                         :: i, j, k
    
    do k = kmin, kmax
      do j = jmin, jmax
        do i = imin, imax
          aaz_Var(i,j,k) = VarMean + VarDel * tanh( (z_Z(k) - Zc) / Zr ) 
        end do
      end do
    end do

    if ( present( aaz_VarBZ) ) then 
      aaz_Var = aaz_Var - aaz_VarBZ
    end if

  end subroutine initialdata_disturb_tanh


  subroutine initialdata_disturb_tanh_sin(VarMean, VarDel, Zc, Zr, aaz_Var, aaz_VarBZ)
    !
    ! tanh 型のシア
    !   A(z) = A0 + A1 \tanh( (z - Zc) / h )

    implicit none

    real(DP), intent(in)            :: VarMean, VarDel, Zc, Zr
    real(DP), intent(in), optional  :: aaz_VarBZ(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), intent(out)           :: aaz_Var(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), parameter             :: Pi = 3.1415926535897932385d0 
    integer                         :: i, j, k
    
    do k = kmin, kmax
      do j = jmin, jmax
        do i = imin, imax
          aaz_Var(i,j,k) =                                  &
            &   VarMean                                     &
            & + VarDel * tanh( (z_Z(k) - Zc) / Zr )         &
            & + VarDel * sin( x_X(i) * 2.0d0 * Pi / XMax )  &
            &          * sin( z_Z(k) * 2.0d0 * Pi / ZMax )
        end do
      end do
    end do

    if ( present( aaz_VarBZ) ) then 
      aaz_Var = aaz_Var - aaz_VarBZ
    end if

  end subroutine initialdata_disturb_tanh_sin


!!!
!!! 円錐形
!!! 

  subroutine initialdata_disturb_ConeXY(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 * (1.0 - sqrt( ((x_X(i) - Xc) / Xr )**2.0d0   &
            &                    +  ((y_Y(j) - Yc) / Yr )**2.0d0) )
          xyz_Var(i,j,k) = MAX(xyz_Var(i,j,k), 0.0d0)
        end do
      end do
    end do

  end subroutine initialdata_disturb_ConeXY


  subroutine initialdata_disturb_ConeXZ(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 * (1.0 - sqrt( ((x_X(i) - Xc) / Xr )**2.0d0   &
            &                    +  ((z_Z(k) - Zc) / Zr )**2.0d0) )

          xyz_Var(i,j,k) = MAX(xyz_Var(i,j,k), 0.0d0)
        end do
      end do
    end do

  end subroutine initialdata_disturb_ConeXZ 


  subroutine initialdata_disturb_ConeYZ(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 * (1.0 - sqrt( ((y_Y(i) - Yc) / Yr )**2.0d0   &
            &                    +  ((z_Z(k) - Zc) / Zr )**2.0d0) )

          xyz_Var(i,j,k) = MAX(xyz_Var(i,j,k), 0.0d0)
        end do
      end do
    end do

  end subroutine initialdata_disturb_ConeYZ

  
!!!
!!! 一定
!!!  
  
  subroutine initialdata_disturb_circleXZ(DelMax, Xc, Xr, Zc, xyz_Var)

    implicit none

    real(DP), intent(in)  :: DelMax, Xc, Zc, Xr
    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
          if ( ( (x_X(i) - Xc) ** 2.0d0 + (z_Z(k) - Zc) ** 2.0d0 ) < Xr*Xr ) then 
            xyz_Var(i,j,k) = DelMax
          else
            xyz_Var(i,j,k) = 0.0d0
          end if
        end do
      end do
    end do
    
  end subroutine initialdata_disturb_circleXZ


!!!
!!! ガウシアン
!!!  
  
  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_cosXZ(DelMax, Xc, Xr, Zc, Zr, xyz_Var)
    !
    ! 重力流の計算に利用する擾乱
    ! A [cos(πL) + 1]*0.5 ( L < 1.0 ) or 0.0
    
    implicit none
    
    real(DP), intent(in)  :: DelMax, Xc, Xr, Zc, Zr
    real(DP), intent(out) :: xyz_Var(imin:imax, jmin:jmax, kmin:kmax)
    real(DP)              :: xyz_Var2(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), parameter   :: Pi = 3.1415926535897932385d0 

    xyz_Var = 0.0d0   
    xyz_Var2 = sqrt(                       &
         &   + ((xyz_X - Xc) / Xr )**2.0d0 &
         &   + ((xyz_Z - Zc) / Zr )**2.0d0 &
         &  )
    
    where ( xyz_Var2 <= 1.0d0 )
       xyz_Var = DelMax * ( cos( Pi * xyz_Var2 ) + 1.0d0 ) * 0.5d0
    end where
    
  end subroutine initialdata_disturb_cosXZ
  
  
  subroutine initialdata_disturb_cosXY(DelMax, Xc, Xr, Yc, Yr, xyz_Var)
    !
    ! 重力流の計算に利用する擾乱
    ! A [cos(πL) + 1]*0.5 ( L < 1.0 ) or 0.0
    
    implicit none

    real(DP), intent(in)  :: DelMax, Xc, Xr, Yc, Yr
    real(DP), intent(out) :: xyz_Var(imin:imax, jmin:jmax, kmin:kmax)
    real(DP)              :: xyz_Var2(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), parameter   :: Pi = 3.1415926535897932385d0 
    
    xyz_Var = 0.0d0   
    xyz_Var2 = sqrt(                       &
         &   + ((xyz_X - Xc) / Xr )**2.0d0 &
         &   + ((xyz_Y - Yc) / Yr )**2.0d0 &
         &  )
    
    where ( xyz_Var2 <= 1.0d0 )
       xyz_Var = DelMax * ( cos( Pi * xyz_Var2 ) + 1.0d0 ) * 0.5d0
    end where
    
  end subroutine initialdata_disturb_cosXY
  

  subroutine initialdata_disturb_cosYZ(DelMax, Yc, Yr, Zc, Zr, xyz_Var)
    !
    ! 重力流の計算に利用する擾乱
    ! A [cos(πL) + 1]*0.5 ( L < 1.0 ) or 0.0
    
    implicit none

    real(DP), intent(in)  :: DelMax, Yc, Yr, Zc, Zr
    real(DP), intent(out) :: xyz_Var(imin:imax, jmin:jmax, kmin:kmax)
    real(DP)              :: xyz_Var2(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), parameter   :: Pi = 3.1415926535897932385d0 

    xyz_Var = 0.0d0   
    xyz_Var2 = sqrt(                       &
         &   + ((xyz_Y - Yc) / Yr )**2.0d0 &
         &   + ((xyz_Z - Zc) / Zr )**2.0d0 &
         &  )

    where ( xyz_Var2 <= 1.0d0 )
       xyz_Var = DelMax * ( cos( Pi * xyz_Var2 ) + 1.0d0 ) * 0.5d0
    end where
    
  end subroutine initialdata_disturb_cosYZ
  
  
  subroutine initialdata_disturb_cosXYZ(DelMax, Xc, Xr, Yc, Yr, Zc, Zr, xyz_Var)
    !
    ! 重力流の計算に利用する擾乱
    ! A [cos(πL) + 1]*0.5 ( L < 1.0 ) or 0.0

    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)
    real(DP)              :: xyz_Var2(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), parameter   :: Pi = 3.1415926535897932385d0 
    
    xyz_Var = 0.0d0   
    xyz_Var2 = sqrt(                       &
         &   + ((xyz_X - Xc) / Xr )**2.0d0 &
         &   + ((xyz_Y - Yc) / Yr )**2.0d0 &
         &   + ((xyz_Z - Zc) / Zr )**2.0d0 &
         &  )
    
    where ( xyz_Var2 <= 1.0d0 )
       xyz_Var = DelMax * ( cos( Pi * xyz_Var2 ) + 1.0d0 ) * 0.5d0
    end where
    
  end subroutine initialdata_disturb_cosXYZ

!!!
!!!---------------------------------------------------------------------
!!!

  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_square( &
    & 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_square
  
  
  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