!= Module BasicEnvInit
!
! Authors::   SUGIYAMA Koichiro, ODAKA Masatsugu
! Version::   $Id: initialdata_takemi2007.f90,v 1.1 2011-04-26 04:08:10 sugiyama Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview 
!
!デフォルトの基本場を設定するための変数参照型モジュール
!   * BasicEnvFile_init: 基本場の値を netCDF ファイルから取得
!   * BasicEnvCalc_Init: 基本場の情報を Namelist から取得して値を計算
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!== Future Plans
!
module initialdata_takemi2007
  !
  !デフォルトの基本場を設定するためのサブルーチン. 
  !基本場を計算し, BasicSet モジュールの値を初期化する. 
  !
  !コンパイルの順序の問題から, 基本場の値(hogeBasicZ な変数)を
  !計算する部分をBasicSet モジュールから切り離している. 
  !ECCM 始め, BasicSet 自体に依存するが hogeBasicZ は use しない
  !外部サブルーチンを利用するためである. 
  !

  !モジュール読み込み
  use dc_types,   only: STRING, DP
  use dc_message, only: MessageNotify
  use gridset,  only: DimXMin,       &!配列の X 方向の下限
    &                 DimXMax,       &!配列の X 方向の上限
    &                 DimZMin,       &!配列の Z 方向の下限
    &                 DimZMax,       &!配列の Z 方向の上限
    &                 RegZMin,       &!配列の Z 方向の下限
    &                 RegZMax,       &!配列の Z 方向の上限
    &                 SpcNum,        &!凝縮成分の数
    &                 z_Z,           &!スカラー格子点での高度
    &                 DelZ            !Z 方向の格子点間隔
  use basicset, only: BasicSetArray_Init, &!
    &                 PressBasis,    &!温位の基準圧力
    &                 GasRDry,       &!乾燥成分の定圧比熱
    &                 CpDry,         &!乾燥成分の定圧比熱
    &                 CvDry,         &!乾燥成分の定積比熱
    &                 MolWtDry,      &!乾燥成分の分子量
    &                 SpcWetMolFr,   &!凝縮成分の初期モル比
    &                 MolWtWet,      &!凝縮成分の分子量
    &                 Grav,          &!重力加速度
    &                 Tropopause,    &!対流圏界面高度
    &                 TempSfc,       &!地表面温度
    &                 PressSfc,      &!地表面圧力
    &                 GasRUniv        !気体定数
  use chemcalc, only: SvapPress       !
 
  !暗黙の型宣言禁止
  implicit none

  !デフォルトは private
  private

  !初期化だけ公開
  public  initialdata_takemi2007_basic
  public  initialdata_takemi2007_wind

contains

!!!------------------------------------------------------------------------------!!!
  subroutine HUM_Takemi2007(z_RH,Ztr,TR,LL)
    !
    !== 概要
    ! * deepconv の地球用のテスト計算としてTakemi(2007)の再現計算を
    !   するための湿度の基本場を作成する
    !   * Takemi(2007) の熱帯場と中緯度場の相対湿度の"BASE CASE"の
    !     高度 1.5km 以上の湿度を計算している(詳細はTakemi(2007)の図2を参照のこと)
    !     * Takemi(2007)では高度1.5km以上はある関数で与えているが,
    !       高度1.5km以下は一定の混合比で与えているため,
    !       高度1.5km以下の湿度は温度や気圧が求まらないと得られないので後で求める
    !
    implicit none
    
    real(DP),intent(out) :: z_RH(DimZMin:DimZMax) !Relative Humidity(相対湿度)
    real(DP),intent(out) :: Ztr                   !対流圏界面の高度
    integer,intent(out) :: TR                    !対流圏界面の高度の格子番号
    integer,intent(out) :: LL                    !高度 1.5 km 直下の格子番号 (Low Level)
    integer :: k
    
    z_RH = 0.0d0
    Ztr  = 0.0d0

    ! Ztrに対流圏界面高度の値を与える
    ! Takemi(2007)では12kmだが, 必ずしも格子点がちょうど12kmには乗らないため,
    ! 12kmか,12kmの直上の位置を対流圏界面にする
    do k = RegZMin+1,DimZMax

      ! z_Z が ちょうど 12 km なら 12 km が対流圏界面
      if (z_Z(k).eq.12000) then
        Ztr = z_Z(k)
        
      ! z_Z がちょうど12kmではなかった場合,12kmの直上のz_Zが対流圏界面Ztrになる
      else if (z_Z(k).gt.12000) then
        Ztr = minval( z_Z, 1, z_Z > 12000 )
      end if
    end do

    ! 対流圏界面高度の格子点を求める (Ztr = z_Z(TR))
    do k = RegZMin+1,DimZMax
      if (z_Z(k).ge.12000.and.z_Z(k).le.(12000 + DelZ)) then
        TR = k
      end if
    end do
    
    ! 高度1.5kmもしくはその直下の格子番号を求める
    do k = RegZMin+1,DimZMax
      if (z_Z(k).le.1500.and.z_Z(k).ge.(1500 - DelZ)) then
        LL = k
      end if
    end do
    
    ! Takemi(2007)の式を用いて相対湿度を求める
    ! とりあえず高度1.5km以下の湿度はゼロにしておく
    do k = RegZMin+1, DimZMax
      if (k.le.LL) then
        z_RH(k) = 0.0d0

      else if (k.le.TR.and.k.gt.LL) then
        z_RH(k) = 1.0d0 - 0.750d0 * (z_Z(k) / Ztr)**1.25

      else if (k.gt.TR) then
        z_RH(k) = 0.25
      end if
    end do
  end subroutine HUM_Takemi2007

!!!------------------------------------------------------------------------------!!!
  subroutine HUM_Takemi2007_TDRY(z_RH,Ztr,TR,LL)
    !
    !== 概要
    ! * deepconv の地球用のテスト計算としてTakemi(2007)の再現計算を
    !   するための湿度の基本場を作成する
    !   * Takemi(2007) の熱帯場と中緯度場の相対湿度の"DRY CASE"の
    !     高度 1.5km 以上の湿度を計算している(詳細はTakemi(2007)の図2を参照のこと)
    !     * Takemi(2007)では高度1.5km以上はある関数で与えているが,
    !       高度1.5km以下は一定の混合比で与えているため,
    !       高度1.5km以下の湿度は温度や気圧が求まらないと得られないので後で求める
    !     * 中緯度場とは乾燥のさせ方が違う
    !

    implicit none

    real(DP),intent(out) :: z_RH(DimZMin:DimZMax) !Relative Humidity(相対湿度)
    real(DP),intent(out) :: Ztr                   !対流圏界面の高度
    integer,intent(out) :: TR                    !対流圏界面の高度の格子番号
    integer,intent(out) :: LL                    !高度 1.5 km 直下の格子番号 (Low Level)
    integer :: k

    z_RH = 0.0d0
    Ztr  = 0.0d0

    ! Ztrに対流圏界面高度の値を与える
    ! Takemi(2007)では12kmだが, 必ずしも格子点がちょうど12kmには乗らないため,
    ! 12kmか,12kmの直上の位置を対流圏界面にする
    do k = RegZMin+1,DimZMax
      
      ! z_Z が ちょうど 12 km なら 12 km が対流圏界面
      if (z_Z(k).eq.12000) then
        Ztr = z_Z(k)

      ! z_Z がちょうど12kmではなかった場合,12kmの直上のz_Zが対流圏界面Ztrになる
      else if (z_Z(k).gt.12000) then
        Ztr = minval( z_Z, 1, z_Z > 12000 )
      end if
    end do

    ! 対流圏界面高度の格子点を求める (Ztr = z_Z(TR))
    do k = RegZMin+1,DimZMax
      if (z_Z(k).ge.12000.and.z_Z(k).le.(12000 + DelZ)) then
        TR = k
      end if
    end do

    ! 高度1.5kmもしくはその直下の格子番号を求める
    do k = RegZMin+1,DimZMax
      if (z_Z(k).le.1500.and.z_Z(k).ge.(1500 - DelZ)) then
        LL = k
      end if
    end do

    ! Takemi(2007)の式を用いて相対湿度を求める
    ! とりあえず高度1.5km以下の湿度はゼロにしておく
    ! さらに, 乾燥大気のエントレインメントを考慮した熱帯の"DRY CASE"では
    ! 高度 2.5 km, 5.0 km, 7.5 km のどこかで相対湿度を 20 % 減らすDRY1,2,3ケースを用意する
    do k = RegZMin+1, DimZMax
      if (k.le.LL) then
        z_RH(k) = 0.0d0

      ! 1.5 km から (2.5 km, 5.0 km, 7.5 km) まではそのまま変化

!     else if (z_Z(k).lt.2500.and.k.gt.LL) then   ! DRY1
!     else if (z_Z(k).lt.5000.and.k.gt.LL) then   ! DRY2
      else if (z_Z(k).lt.7500.and.k.gt.LL) then   ! DRY3
        z_RH(k) = 1.0d0 - 0.750d0 * (z_Z(k) / Ztr)**1.25

      ! (2.5 km, 5.0 km, 7.5 km)で20 % 湿度を減少させる 
!     else if (z_Z(k).ge.2500) then   ! DRY1
!     else if (z_Z(k).ge.5000) then   ! DRY2
      else if (z_Z(k).ge.7500) then   ! DRY3
        z_RH(k) = 1.0d0 - 0.750d0 * (z_Z(k) / Ztr)**1.25 - 0.20d0

      end if
    end do

    ! 上空で相対湿度が 25% 以下になる場合は 25% で固定する
    do k = RegZMin+1, DimZMax
      if (z_RH(k).le.(0.25).and.k.gt.LL) then
        z_RH(k) = 0.250d0
      end if
    end do

  end subroutine HUM_Takemi2007_TDRY

!!!------------------------------------------------------------------------------!!!
  subroutine HUM_Takemi2007_MDRY(z_RH,Ztr,TR,LL)
    !
    !== 概要
    ! * deepconv の地球用のテスト計算としてTakemi(2007)の再現計算を
    !   するための湿度の基本場を作成する
    !   * Takemi(2007) の熱帯場と中緯度場の相対湿度の"DRY CASE"の
    !     高度 1.5km 以上の湿度を計算している(詳細はTakemi(2007)の図2を参照のこと)
    !     * Takemi(2007)では高度1.5km以上はある関数で与えているが,
    !       高度1.5km以下は一定の混合比で与えているため,
    !       高度1.5km以下の湿度は温度や気圧が求まらないと得られないので後で求める
    !     * 熱帯場とは乾燥のさせ方が違う

    implicit none

    real(DP),intent(out) :: z_RH(DimZMin:DimZMax) !Relative Humidity(相対湿度)
    real(DP),intent(out) :: Ztr                   !対流圏界面の高度
    integer,intent(out) :: TR                    !対流圏界面の高度の格子番号
    integer,intent(out) :: LL                    !高度 1.5 km 直下の格子番号 (Low Level)
    integer :: k

    z_RH = 0.0d0
    Ztr  = 0.0d0

    ! Ztrに対流圏界面高度の値を与える
    ! Takemi(2007)では12kmだが, 必ずしも格子点がちょうど12kmには乗らないため,
    ! 12kmか,12kmの直上の位置を対流圏界面にする
    do k = RegZMin+1,DimZMax

      ! z_Z が ちょうど 12 km なら 12 km が対流圏界面
      if (z_Z(k).eq.12000) then
        Ztr = z_Z(k)

      ! z_Z がちょうど12kmではなかった場合,12kmの直上のz_Zが対流圏海面Ztrになる
      else if (z_Z(k).gt.12000) then
        Ztr = minval( z_Z, 1, z_Z > 12000 )
      end if
    end do

    ! 対流圏界面高度の格子点を求める (Ztr = z_Z(TR))
    do k = RegZMin+1,DimZMax
      if (z_Z(k).ge.12000.and.z_Z(k).le.(12000 + DelZ)) then
        TR = k
      end if
    end do

    ! 高度1.5kmもしくはその直下の格子番号を求める
    do k = RegZMin+1,DimZMax
      if (z_Z(k).le.1500.and.z_Z(k).ge.(1500 - DelZ)) then
        LL = k
      end if
    end do


    ! Takemi(2007)の式を用いて相対湿度を求める
    ! とりあえず高度1.5kmまでの湿度はゼロにしておく
    ! 乾燥大気のエントレインメントを考慮した中緯度場の"DRY CASE"では
    ! 高度 2.5 km で相対湿度を 13 % もしくは 30 % 減らす DRY1,2 ケースを用意する.
    do k = RegZMin+1, DimZMax
      if (k.le.LL) then
        z_RH(k) = 0.0d0

      ! 1.5 km から2.5 km まではそのまま変化
      else if (z_Z(k).lt.2500.and.k.gt.LL) then
        z_RH(k) = 1.0d0 - 0.750d0 * (z_Z(k) / Ztr)**1.25

      ! 2.5 km 以降は 13 % or 30 % 減らす
      else if (z_Z(k).ge.2500) then
!       z_RH(k) = 1.0d0 - 0.750d0 * (z_Z(k) / Ztr)**1.25 - 0.13d0   ! DRY1
        z_RH(k) = 1.0d0 - 0.750d0 * (z_Z(k) / Ztr)**1.25 - 0.30d0   ! DRY2

      end if
    end do

    ! 上空で相対湿度が 25% 以下になる場合は 25% で固定する
    do k = RegZMin+1, DimZMax
      if (z_RH(k).le.(0.25).and.k.gt.LL) then
        z_RH(k) = 0.250d0
      end if
    end do

  end subroutine HUM_Takemi2007_MDRY


!!!------------------------------------------------------------------------------!!!
  subroutine  initialdata_takemi2007_basic( z_Temp, z_Press, za_MolFr )
    !
    !== 概要
    ! * deepconv の地球用のテスト計算としてTakemi(2007)の再現計算を
    !   するための湿度の基本場を作成する
    !   * 基本場の温度の式が温位で与えられているため, 温度に変換する必要がある
    !

    implicit none

    real(DP), intent(out):: z_Press(DimZMin:DimZMax)!圧力
    real(DP), intent(out):: z_Temp(DimZMin:DimZMax) !温度
    real(DP), intent(out):: za_MolFr(DimZMin:DimZMax, 1:SpcNum) !モル比
    !
    ! ECCM_Takemiでのみ必要な変数
    !
    real(DP) :: z_TakemiRH(DimZMin:DimZMax)         !高度 1.5km 以上の相対湿度
    real(DP) :: z_TakemiTheta(DimZMin:DimZMax)      !基本場の温位
    real(DP) :: Theta_0                             ! 地表面の温位
    real(DP) :: Theta_tr                            ! 対流圏界面の温位
    real(DP) :: TakemiZtr                           ! 対流圏界面の高度
    real(DP) :: HumSfc                              ! 地表面の相対湿度
    real(DP) :: a_MolFrSfc                          ! 地表面のモル比
    real(DP) :: z_LLhumidity(DimZMin:DimZMax)       ! 1.5 km までの相対湿度
!!--------------------------------------------------
    integer             :: k, s, i, m, b, T, L

    z_TakemiRH = 0.0d0
    z_TakemiTheta = 0.0d0

    ! 湿度場の選択
    call HUM_Takemi2007(z_TakemiRH,TakemiZtr,T,L)
!   call HUM_Takemi2007_TDRY(z_TakemiRH,TakemiZtr,T,L)
!   call HUM_Takemi2007_MDRY(z_TakemiRH,TakemiZtr,T,L)

    ! 温度・湿度・気圧を求める

    !地表面の温位は300K
    Theta_0 = 300.0d0

    ! 中緯度場と熱帯場の選択
    ! 対流圏界面の温位で区別している
!   Theta_tr = 343.0d0     ! 中緯度
    Theta_tr = 358.0d0     ! 熱帯

    ! 大気最下層の温位
    z_TakemiTheta(RegZMin+1) = Theta_0 + &
                           & (Theta_tr - Theta_0)*(z_Z(RegZMin+1) / TakemiZtr)**1.25

    z_Temp  = 1.0d-60
    z_Press = 1.0d-60
    z_LLhumidity = 0.0d0

    ! 大気最下層の気温・気圧・水蒸気混合比・分子量の計算

    ! 地表面のモル比
    ! 10g/kg〜18g/kgの混合比で与えられているためモル比に直す
    ! なお熱帯場はQ18のみ使用
    a_MolFrSfc = 0.0180d0 * MolWtDry / MolWtWet(1)  ! Q18
!   a_MolFrSfc = 0.0160d0 * MolWtDry / MolWtWet(1)  ! Q16
!   a_MolFrSfc = 0.0140d0 * MolWtDry / MolWtWet(1)  ! Q14
!   a_MolFrSfc = 0.0120d0 * MolWtDry / MolWtWet(1)  ! Q12
!   a_MolFrSfc = 0.0100d0 * MolWtDry / MolWtWet(1)  ! Q10 

    ! 地表面の相対湿度
    ! モル比から逆算. 気圧を求めるのに必要
    HumSfc = a_MolFrSfc * PressSfc / SvapPress(6,TempSfc)

    ! 大気最下層の気圧
    z_Press(RegZMin+1) = &
           & PressSfc - (Grav * PressSfc * DelZ) &
           & /((GasRUniv/ &
           & (MolWtDry + a_MolFrSfc * (MolWtWet(1) - MolWtDry))) &
           & *TempSfc)

    ! 大気最下層の温度.　温位を元に計算
    ! 地表面気圧は 1000hPa と仮定する
    z_Temp(RegZMin+1) = z_TakemiTheta(RegZMin+1) &
                        &* (z_Press(RegZMin+1) / 100000.0d0)**(GasRDry / CpDry)

    ! 大気最下層のモル比
    ! 熱帯場はQ18のみ使用
    za_MolFr(RegZMin+1,1) = 0.0180d0 * MolWtDry / MolWtWet(1)  ! Q18
!   za_MolFr(RegZMin+1,1) = 0.0160d0 * MolWtDry / MolWtWet(1)  ! Q16
!   za_MolFr(RegZMin+1,1) = 0.0140d0 * MolWtDry / MolWtWet(1)  ! Q14
!   za_MolFr(RegZMin+1,1) = 0.0120d0 * MolWtDry / MolWtWet(1)  ! Q12
!   za_MolFr(RegZMin+1,1) = 0.0100d0 * MolWtDry / MolWtWet(1)  ! Q10
 
    ! 大気最下層の相対湿度
    ! モル比から逆算. 気圧を求めるのに必要
    z_LLhumidity(RegZMin+1) = za_MolFr(RegZMin+1,1) &
                              &* z_Press(RegZMin+1) / SvapPress(6,z_Temp(RegZMin+1))

    ! 大気最下層の種々の物理量をもとに温度・湿度気圧の基本場を作成する
    ! 高度 0 - 1.5km と 1.5km - で湿度の与え方が異なり,
    ! さらに, 対流圏界面以上では温度の式が異なるので, 以下の 3パターンに分けて
    ! 基本場を求める
    ! 高度 0 - 1.5km, 1.5km - 対流圏界面(Ztr), Ztr以上 
    do k = RegZMin+1, DimZMax

    !
    ! 基本場(I)
    ! 大気最下層から高度 1.5 km 以下(k が L-1 以下)の基本場
    !
     if (k.le.L-1) then

      ! モル比
      ! 高度 1.5km まではモル比が任意の値で固定されている
      ! ただし, 熱帯場は 18 g/kg のみ
      za_MolFr(k+1,1) = 0.0180d0 * MolWtDry / MolWtWet(1)
!     za_MolFr(k+1,1) = 0.0160d0 * MolWtDry / MolWtWet(1)
!     za_MolFr(k+1,1) = 0.0140d0 * MolWtDry / MolWtWet(1)
!     za_MolFr(k+1,1) = 0.0120d0 * MolWtDry / MolWtWet(1)
!     za_MolFr(k+1,1) = 0.0100d0 * MolWtDry / MolWtWet(1)
 
      ! 気圧
      z_Press(k+1) = &
           & z_Press(k)-(Grav*z_Press(k)*DelZ) &
           & /((GasRUniv/ &
           & (MolWtDry + za_MolFr(k,1) * (MolWtWet(1) - MolWtDry))) &
           & *z_Temp(k))

      ! 温位から温度を計算
      z_TakemiTheta(k+1) = Theta_0 + (Theta_tr - Theta_0)*(z_Z(k+1) / TakemiZtr)**1.25
      z_Temp(k+1) = z_TakemiTheta(k+1) * (z_Press(k+1) / 100000.0d0)**(GasRDry / CpDry)

      ! モル比から相対湿度を計算
      z_LLhumidity(k+1) = za_MolFr(k+1,1) * z_Press(k+1) / SvapPress(6,z_Temp(k+1))

      !
      !== Takemi(2007)の基本場の使用上の注意
      !  * 低層で混合比を固定し続けると高度 1.5km までに相対湿度が100%を超える高度がある
      !  * 100% を超えてしまった高度の湿度は, その直下の高度の湿度にすると
      !    Takemi(2007)と同じ湿度の鉛直プロファイルが描ける
      !    * なぞ
      !  * MQ16 は 湿度 99.・・・% という値が出て, 100% を超えなくてもほぼ100%
      !    な値が出てしまう
      !    * Takemi(2007)ではこのような値があるようには見えない
      !      しょうがないのでMQ18 とMQ14 の 100% を超える直前の値の間を取った値を入れておく
      !      * 0.94511507230d0という怪しい数字が上記の値
      if (z_LLhumidity(k+1).gt.1) then
        z_LLhumidity(k+1) = z_LLhumidity(k)
!     if (z_LLhumidity(k+1).gt.0.99) then      !MQ16 用
!       z_LLhumidity(k+1) = 0.94511507230d0    !MQ16 用
     end if

     !
     ! 基本場(II)
     ! 高度 1.5km より上で, 対流圏界面以下(L < k < T)の基本場
     !
     else if (k.le.T.and.k.gt.L) then

      ! k = L での相対湿度
      ! この if 文の初めの k は L+1 から始まり z_TakemiRH(L) が区間外なので値を与える
      z_TakemiRH(L) = z_LLhumidity(L)

      ! 気圧の基本場
      z_Press(k) = &
           & z_Press(k-1)-(Grav*z_Press(k-1)*DelZ) &
           & /((GasRUniv/ &
           & (MolWtDry + za_MolFr(k+1,1) * (MolWtWet(1) - MolWtDry))) &
           & *z_Temp(k-1))

       ! 温位から温度を計算
       z_TakemiTheta(k) = Theta_0 + (Theta_tr - Theta_0)*(z_Z(k) / TakemiZtr)**1.25
       z_Temp(k) = z_TakemiTheta(k) * (z_Press(k) / 100000.0d0)**(GasRDry / CpDry)

       ! モル比と大気の平均分子量
       za_MolFr(k,1) = SvapPress( 6,z_Temp(k))* z_TakemiRH(k)/z_Press(k)
     
       !
       ! 基本場(III)
       ! 対流圏界面より上 (k > T) の場
       !
       else if (k.gt.T) then

       ! 気圧の基本場
       z_Press(k) = &
           & z_Press(k-1)-(Grav*z_Press(k-1)*DelZ) &
           & /((GasRUniv/ &
           & (MolWtDry + za_MolFr(k,1) * (MolWtWet(1) - MolWtDry))) &
           & *z_Temp(k-1))

       ! 温位から温度を計算
       z_TakemiTheta(k) = Theta_tr * exp(Grav * (z_Z(k) - TakemiZtr) / (CpDry * z_Temp(T)))
       z_Temp(k) = z_TakemiTheta(k) * (z_Press(k) / 100000.0d0)**(GasRDry / CpDry)

       ! モル比と大気の平均分子量
       za_MolFr(k,1) = SvapPress( 6,z_Temp(k))* z_TakemiRH(k)/z_Press(k)

     end if
    end do

  end subroutine Initialdata_takemi2007_basic


  subroutine initialdata_takemi2007_wind

    !-------------------------------------------------------------!
    ! シアーの設定 (Takemi,2007)                                  !
    !-------------------------------------------------------------!
    !
    != 概要
    !* case "Takemi2007" での計算時に鉛直シアーのある風を与える時に使用する
    !* 風の与え方には, 以下のようなバリエーションがある
    !  (1) シアーを与える高度を変える
    !  (2) シアーのある風の最大風速 (U_s) を変える
    !
    !  (1) については, (a) 0 - 2.5 km, (b) 2.5 - 5.0 km, (c) 5.0 - 7.5 km の
    !  三パターンがある
    !  (2) については, Takemi (2007) では熱帯場と中緯度場の温度場毎に
    !  異なる値を設定している
    !
    !  その強度(Us)は, 以下の通り
    !  <熱帯場>   (1) 5 m/s, (2) 10 m/s, (3) 15 m/s
    !  <中緯度場> (1) 10 m/s, (2) 15 m/s, (3) 20 m/s
    !
    !* シアーの形の模式図 (Takemi, 2007)   |
    !                                     /| 7.5 km
    !                                    / |
    !                                   /  |
    !                                  / ←|
    !                                 ｜  /| 5.0 km
    !                                 ｜ / |
    !                                 ｜/  |
    !                                  / ←|
    !                                 ｜  /| 2.5 km
    !                                 ｜ / |
    !                                 ｜/  |
    !                                  / ←|
    !---------------------------------+------------- 0.0 km
    !                                Us (m/s)
    !
    
    ! シアーを与える高度を決める
    ! 傾き始めの高度の格子点 (B) と 傾きの終わりの格子点 (U) で決定
    
    do k = RegZMin+1, DimZMax
      if (s_Z(k).ge.Hb.and.s_Z(k).le.(Hb + DelZ)) then
        B = k
      end if
      
      if (s_Z(k).ge.Hu.and.s_Z(k).le.(Hu + DelZ)) then
        U = k
      end if
    end do
    
    ! 上記で決めた格子点BからUの間に鉛直シアーのある風を与える
    
    do k = RegZMin+1, U
      do i = DimXMin, DimXMax
        if (k.lt.B) then
          pz_VelX(i,k) = - Us
        else
          pz_VelX(i,k) = - Us + (Us / (s_Z(U) - s_Z(B))) * (s_Z(k) - s_Z(B))
        end if
      end do
    end do
    
  end subroutine initialdata_takemi2007_wind
  
end module initialdata_takemi2007
