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_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_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_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_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_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_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_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_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_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_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

  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_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_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