!= Module BasicEnvInit
!
! Authors::   SUGIYAMA Koichiro, ODAKA Masatsugu
! Version::   $Id: bascset.f90,v 1.1 2011-03-30 13:53:06 sugiyama Exp $
! Tag Name::  $Name: arare4-20120911 $
! 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 Bscset
  !
  !デフォルトの基本場を設定するためのサブルーチン. 
  !基本場を計算し, BasicSet モジュールの値を初期化する. 
  !
  !コンパイルの順序の問題から, 基本場の値(hogeBasicZ な変数)を
  !計算する部分をBasicSet モジュールから切り離している. 
  !ECCM 始め, BasicSet 自体に依存するが hogeBasicZ は use しない
  !外部サブルーチンを利用するためである. 
  !

  !モジュール読み込み
  use dc_message, only: MessageNotify
  use gtool_history   
  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 Boundary,only:  BoundaryXCyc_xz,  &!
    &                 BoundaryZSym_xz,  &!  
    &                 BoundaryXCyc_xza, &!
    &                 BoundaryZSym_xza   !  
  use ECCM,     only: ECCM_Dry,              &!
    &                 ECCM_Wet              
  use chemcalc, only: SvapPress,        &!
    &                 LatentHeatPerMol, &!
    &                 ReactHeatNH4SHPerMol
    
  !暗黙の型宣言禁止
  implicit none

  !デフォルトは private
  private

  real(8), save  :: DryHeight = 0.0d0

  !初期化だけ公開
  public bscset_init

contains


  subroutine ECCM_init(cfgfile)
    !
    ! ECCM モジュールの初期化ルーチン
    !

    !暗黙の型宣言禁止
    implicit none

    !
    !NAMELIST から放射強制の設定を取得
    !

    character(*), intent(in) :: cfgfile
    
    ! NAMELIST から情報を取得
    NAMELIST /eccm/ DryHeight
    
    open (10, FILE=cfgfile)
    read(10, NML=eccm)
    close(10)
    
    !確認
    call MessageNotify( "M", "ECCM_init", "DryHeight = %f", d=(/DryHeight/) )
    
  end subroutine ECCM_init


  subroutine bscset_init(cfgfile)

    !暗黙の型宣言禁止
    implicit none
    
    !変数の定義
    character(*), intent(in) :: cfgfile
    real(8)                  :: xz_DensBasicZ(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)                  :: xz_PressBasicZ(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)                  :: xz_ExnerBasicZ(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)                  :: xz_TempBasicZ(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)                  :: xz_PotTempBasicZ(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)                  :: xz_VelSoundBasicZ(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)                  :: xza_MixRtBasicZ(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
    real(8)                  :: xz_EffMolWtBasicZ(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)                  :: z_TempBasicZ(DimZMin:DimZMax)
    real(8)                  :: z_PressBasicZ(DimZMin:DimZMax)
    real(8)                  :: MolFrIni(SpcNum)
    real(8)                  :: xza_MolFr(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
    real(8)                  :: za_MolFr(DimZMin:DimZMax, SpcNum)
    real(8)                  :: xza_MixRtDivMolWt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
    real(8)                  :: Humidity = 0.0d0          !相対湿度 
    character(20)            :: Type = ""                 !基本場の温度設定, 'Dry' or 'Moist'   
    real(8)                  :: TempStrat = 200.0d0       !成層圏の温度 [k]
    real(8)                  :: Dhight = 5.0d3            !重み関数のパラメータ [m]  
    real(8)                  :: z_llhum(DimZMin:DimZMax)
    integer                  :: i, s
     
    !---------------------------------------------------------------
    ! NAMELIST の定義
    !---------------------------------------------------------------
    NAMELIST /basicenv/ Type, Humidity, TempStrat, Dhight
    
    !---------------------------------------------------------------
    ! 配列の初期化
    !---------------------------------------------------------------
    xz_PressBasicZ    = 0.0d0
    xz_ExnerBasicZ    = 0.0d0
    xz_TempBasicZ     = 0.0d0
    xz_PotTempBasicZ  = 0.0d0
    xz_VelSoundBasicZ = 0.0d0
    xza_MixRtBasicZ   = 0.0d0
    xz_EffMolWtBasicZ = 0.0d0
    z_TempBasicZ      = 0.0d0
    z_PressBasicZ     = 0.0d0
    za_MolFr          = 0.0d0
    
    !---------------------------------------------------------------
    ! Type を元に, 温度, 圧力, 組成を決める
    !---------------------------------------------------------------
    MolFrIni = SpcWetMolFr(1:SpcNum) 
    open (10, FILE=cfgfile)
    read(10, NML=basicenv)
    close(10)
    
    select case(Type)
    case("Dry") ! 乾燥断熱的な初期場
      call ECCM_Dry( MolFrIni, Humidity, z_TempBasicZ, z_PressBasicZ, za_MolFr )
      call bscset_upper(Humidity, TempStrat, Dhight, z_TempBasicZ, z_PressBasicZ)
      
    case("Wet") ! 湿潤断熱的な初期場
      call ECCM_Wet( MolFrIni, Humidity, z_TempBasicZ, z_PressBasicZ, za_MolFr )
      call bscset_upper(Humidity, TempStrat, Dhight, z_TempBasicZ, z_PressBasicZ) 
      
    case("Yamasaki1983") ! Yamasaki(1983)の温度と相対湿度の観測値を使用する場合 
      call ECCM_Ymazaki1983( z_TempBasicZ, z_PressBasicZ, za_MolFr )
      
    case("Takemi2007") ! Takemi(2007)の基本場を使用する場合
      call ECCM_Takemi2007( z_TempBasicZ, z_PressBasicZ, za_MolFr, z_llhum )
    end select
    
    ! 2 次元配列に格納
    do i = DimXMin, DimXMax
      xz_TempBasicZ(i,:)  = z_TempBasicZ  
      xz_PressBasicZ(i,:) = z_PressBasicZ  
    end do
    
    !境界条件
    call BoundaryXCyc_xz( xz_TempBasicZ )
    call BoundaryZSym_xz( xz_TempBasicZ )
    call BoundaryXCyc_xz( xz_PressBasicZ )
    call BoundaryZSym_xz( xz_PressBasicZ )
    
    !---------------------------------------------------------------
    ! 混合比
    !---------------------------------------------------------------
    !水平方向には一様
    do i = DimXMin, DimXMax      
      xza_MolFr(i,:,:) = za_MolFr
    end do
    
    !気相のモル比を混合比に変換
    do s = 1, SpcNum
      xza_MixRtBasicZ(:,:,s) = xza_MolFr(:,:,s) * MolWtWet(s) / MolWtDry
    end do
    
    !境界条件
    call BoundaryXCyc_xza( xza_MixRtBasicZ )
    call BoundaryZSym_xza( xza_MixRtBasicZ )
    
    !---------------------------------------------------------------
    ! 分子量の効果
    !---------------------------------------------------------------
    do s = 1, SpcNum
      xza_MixRtDivMolWt(:,:,s) = xza_MixRtBasicZ(:,:,s) / MolWtWet(s)
    end do
    
    xz_EffMolWtBasicZ = &
      & (1.0d0 + sum(xza_MixRtBasicZ,3) ) &
      & / ( MolWtDry * ((1.0d0 / MolWtDry) + sum(xza_MixRtDivMolWt,3)) )
    
    !境界条件  
    call BoundaryXCyc_xz( xz_EffMolWtBasicZ )
    call BoundaryZSym_xz( xz_EffMolWtBasicZ )  
    
    !---------------------------------------------------------------    
    ! 温位
    !---------------------------------------------------------------
    xz_PotTempBasicZ = &
      & xz_TempBasicZ * (PressBasis / xz_PressBasicZ) ** (GasRDry / CpDry) 
    
    !境界条件  
    call BoundaryXCyc_xz( xz_PotTempBasicZ )
    call BoundaryZSym_xz( xz_PotTempBasicZ )
    
    !---------------------------------------------------------------    
    ! エクスナー関数
    !---------------------------------------------------------------
    xz_ExnerBasicZ = xz_TempBasicZ / xz_PotTempBasicZ    
    
    !境界条件
    call BoundaryXCyc_xz( xz_ExnerBasicZ )
    call BoundaryZSym_xz( xz_ExnerBasicZ )
    
    !---------------------------------------------------------------    
    ! 密度
    !---------------------------------------------------------------
    xz_DensBasicZ = &
      & PressBasis * (xz_ExnerBasicZ ** (CvDry / GasRDry)) &
      &  / (GasRDry * xz_PotTempBasicZ / xz_EffMolWtBasicZ)
    
    !境界条件
    call BoundaryXCyc_xz( xz_DensBasicZ )
    call BoundaryZSym_xz( xz_DensBasicZ )
    
    !---------------------------------------------------------------    
    ! 音速
    !---------------------------------------------------------------
    xz_VelSoundBasicZ = &
      & sqrt ( &
      &   CpDry * GasRDry * xz_ExnerBasicZ * xz_PotTempBasicZ &
      &   / (CvDry * xz_EffMolWtBasicZ) &
      & )
    
    !境界条件
    call BoundaryXCyc_xz( xz_VelSoundBasicZ )
    call BoundaryZSym_xz( xz_VelSoundBasicZ )
    
    !----------------------------------------------------------
    ! BasicSet モジュールに値を設定
    !----------------------------------------------------------
    call BasicSetArray_Init(                                  &
      & xz_PressBasicZ,    xz_ExnerBasicZ, xz_TempBasicZ,     &
      & xz_PotTempBasicZ,  xz_DensBasicZ,  xz_VelSoundBasicZ, &
      & xza_MixRtBasicZ, xz_EffMolWtBasicZ )

  end subroutine bscset_init


  subroutine bscset_upper(Humidity, TempStrat, Dhight, z_TempBasicZ, z_PressBasicZ)
    !
    !
    !デフォルトの基本場を設定するためのサブルーチン. 
    !基本場を計算し, BasicSet モジュールの値を初期化する. 
    !
    !コンパイルの順序の問題から, 基本場の値(hogeBasicZ な変数)を
    !計算する部分をBasicSet モジュールから切り離している. 
    !ECCM 始め, BasicSet 自体に依存するが hogeBasicZ は use しない
    !外部サブルーチンを利用するためである. 
    !
    
    implicit none
    
    integer :: k
    real(8), intent(in)    :: Humidity      !基本場の相対湿度
    real(8), intent(in)    :: TempStrat     !成層圏の温度 [k]
    real(8), intent(in)    :: Dhight        !重み関数のパラメータ [m]
    real(8), intent(inout) :: z_TempBasicZ(DimZMin:DimZMax)
    real(8), intent(inout) :: z_PressBasicZ(DimZMin:DimZMax)
    real(8) :: z_DTempDZ(DimZMin:DimZMax)
    real(8) :: DTempDZ
    real(8) :: Weight
    
    
    do k = RegZMin+1, DimZMax-1     
      z_DTempDZ(k) = (z_TempBasicZ(k) - z_TempBasicZ(k-1)) / DelZ
    end do
    
    ! 対流圏界面より上の扱い
    !   圏界面より上は等温大気とする. 等温位大気から等温大気への遷移は 
    !   tanh を用いてなめらかにつなぐ
    do k = RegZMin+1, DimZMax
      
      !重みつけの関数を用意. tanh を用いる
      Weight = ( tanh( (z_Z(k) - tropopause ) / Dhight ) + 1.0d0 ) * 5.0d-1
      
      !仮値として温度を計算する. 圏界面より上では TempStrat の等温大気に近づける
      z_TempBasicZ(k) = z_TempBasicZ(k) * ( 1.0d0 - Weight ) + TempStrat * Weight
      
      !温度減率が断熱温度減率より小さくならないように
      DTempDZ = &
        & max( z_DTempDZ(k), (z_TempBasicZ(k) - z_TempBasicZ(k-1)) / DelZ )
      
      !基本場の温度を決める
      z_TempBasicZ(k) = z_TempBasicZ(k-1) + DTempDZ * DelZ 
      
      !圧力を静水圧平衡から計算
      z_PressBasicZ(k) =                                                    &
        &  z_PressBasicZ(k-1) * ( ( z_TempBasicZ(k-1) / z_TempBasicZ(k) )   &
        &    ** (Grav * MolWtDry / ( DTempDZ * GasRUniv ) ) )
    end do
    
  end subroutine Bscset_upper

!!------------------------------------------------------------------------------!!!
  subroutine ECCM_Yamasaki1983( z_Temp, z_Press, za_MolFr )
    !
    !== 概要
    !  * Yamasaki, 1983 の温度と相対湿度の観測値で基本場を作成する
    !    * 温度の基本場
    !      * 観測データをNetCDFファイル化したものから値を読み込む
    !        * 読み込んだ値を線形補間して温度の基本場を作成する
    !    * 湿度の基本場
    !      * 相対湿度は subroutine HUM で作成済み
    !        * ここでは相対湿度をモル比に変換して湿度の基本場を作成する 
    !    * 気圧の基本場
    !      * 湿潤成分と乾燥成分の分子量差を考慮した静水圧平衡から計算する
    !

    implicit none

    real(8), intent(out):: z_Press(DimZMin:DimZMax)!圧力
    real(8), intent(out):: z_Temp(DimZMin:DimZMax) !温度
    real(8), intent(out):: za_MolFr(DimZMin:DimZMax, 1:SpcNum) !モル比
    real(8)             :: Ob_altitude(1:89)       !観測データに対応する高度
    real(8)             :: Ob_HumZ(1:88)           !相対湿度の観測値(1次元)
    real(8)             :: z_HumZ(DimZMin:DimZMax) !相対湿度の基本場の値(%)
    real(8)        :: Observed_Temp(1:100,1:89)!NetCDF から読み込んだ気温の観測値(2次元)
    real(8)        :: Ob_TempZ(1:89)           !NetCDF から読み込んだ気温の観測値(1次元)
    real(8)             :: z_PressSaturated(DimZMin:DimZMax)!Tetens(1930) の飽和蒸気圧
    integer             :: k, s, i, m, b

    z_HumZ = 0.0d0

    ! subroutine HUMから相対湿度の基本場を取得
    call HUM(z_HumZ, Ob_altitude, Ob_HumZ)

    ! 温度データをNetCDFから取得し線形補間する

    ! NetCDF ファイルからの気温の観測値の読み込み
    call HistoryGet('Ob_Temp.nc', 'Ob_Temp', array=Observed_Temp)
    
    ! 温度の観測値のファイル"Ob_Temp.nc"は2次元データ(x,z)なので
    ! 1次元データ(z)にする
    do m=1,89
      Ob_TempZ(m) = Observed_Temp(1,m)
    end do

    ! Ob_TempZ を線形補間して温度の基本場 z_Temp を作成する
    ! kansoku_altitude(m) < z_Z(k) < kansoku_altitude(m+1),
    ! kansoku_altitude = z_Z(k),
    ! kansoku_altitude < z_Z(k)
    ! の3つで場合分け

    ! 初期化
    z_Temp = 1.0d-60      

    do k = RegZMin+1,DimZMax
     do m = 1,88

      ! z_Z が Ob_altitude のある区間に挟まれるとき,
      ! その区間の Obaltitude(m), Obaltitude(m+1) を結ぶ
      ! 直線で Ob_TempZ を線形補間する
      if (Ob_altitude(m).ne.z_Z(k).and.z_Z(k).gt.Ob_altitude(m) &
          & .and.z_Z(k).lt.Ob_altitude(m+1)) then
        z_Temp(k) = Ob_TempZ(m)+((Ob_TempZ(m+1)-Ob_TempZ(m)) &
        & /(Ob_altitude(m+1)-Ob_altitude(m)))*(z_Z(k)-Ob_altitude(m))

      ! z_Z(k) = Obaltitude(m) の時は, Ob_TempZ(m) = z_Temp(k) である
      else if (Ob_altitude(m).eq.z_Z(k)) then
        z_Temp(k) = Ob_TempZ(m)

      ! z_Z(k) > Ob_altitude では観測データが無いので等温大気にする
      else if (Ob_altitude(m).lt.z_Z(k)) then
        z_Temp(k) = z_Temp(k-1)

      end if
     end do
    end do

    ! 気圧とモル比を計算する

    ! 初期化
    z_Press = 1.0d-60

    ! 地表面気圧と温度から大気最下層の気圧を求める
    ! 静水圧の式 dP/dz = - \rho * g を使用
    ! 簡単のため, \rho = 1 kg/m^3 の乾燥大気を仮定している
    z_Press(RegZMin+1) = PressSfc - (Grav * DelZ)

    ! 大気最下層のモル比を計算する
    za_MolFr(RegZMin+1,1) = SvapPress( 6,z_Temp(RegZMin+1)) * &
                           &(z_HumZ(RegZMin+1)/100)/z_Press(RegZMin+1)

    ! Tetens(1930) の飽和蒸気圧も使ってみたいので以下で計算
    do k=RegZMin+1,DimZMax
      z_PressSaturated(k) = 6.11*(10**(7.5*(z_Temp(k)-273)/(z_Temp(k)-35.7)))*100
    end do

    ! 大気最下層の値を元に, 気圧と湿度の基本場を求める
    ! 気圧は湿潤成分と乾燥成分の分子量差を考慮した静水圧平衡から計算する
    ! 相対湿度はモル比に変換する
    ! 以下のif文のDryHeightは基本場に乾燥した層を入れる場合に必要な値
    ! NAMELISTで鉛直領域以上の値を指定しておけば乾燥した層は入らない

    ! 湿度の基本場の計算
    do k = RegZMin+1,DimZMax-1
     if (z_Z(k+1).lt.DryHeight) then
        za_MolFr(k+1,1) = SvapPress( 6,z_Temp(k)) * (z_HumZ(k)/100)/z_Press(k)

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

      !DryHeightより上に領域があるときは上で求めたモル比をゼロにして
      !乾燥大気における静水圧の式から気圧を求める
      else if (z_Z(k+1).ge.DryHeight) then
        za_MolFr(k+1,1) = 0.0d0

        z_Press(k+1) = &
           & z_Press(k)-(Grav*z_Press(k)*DelZ) &
           & /((GasRUniv/MolWtDry)*z_Temp(k))
      end if
    end do


  end subroutine ECCM_YAMASAKI1983

!!!------------------------------------------------------------------------------!!!
  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(8),intent(out) :: z_RH(DimZMin:DimZMax) !Relative Humidity(相対湿度)
    real(8),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(8),intent(out) :: z_RH(DimZMin:DimZMax) !Relative Humidity(相対湿度)
    real(8),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(8),intent(out) :: z_RH(DimZMin:DimZMax) !Relative Humidity(相対湿度)
    real(8),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 ECCM_Takemi2007( z_Temp, z_Press, za_MolFr, z_LLhumidity )
    !
    !== 概要
    ! * deepconv の地球用のテスト計算としてTakemi(2007)の再現計算を
    !   するための湿度の基本場を作成する
    !   * 基本場の温度の式が温位で与えられているため, 温度に変換する必要がある
    !

    implicit none

    real(8), intent(out):: z_Press(DimZMin:DimZMax)!圧力
    real(8), intent(out):: z_Temp(DimZMin:DimZMax) !温度
    real(8), intent(out):: za_MolFr(DimZMin:DimZMax, 1:SpcNum) !モル比
    real(8), intent(out) :: z_LLhumidity(DimZMin:DimZMax)  ! 1.5 km までの相対湿度
    !
    ! ECCM_Takemiでのみ必要な変数
    !
    real(8) :: z_TakemiRH(DimZMin:DimZMax)         !高度 1.5km 以上の相対湿度
    real(8) :: z_TakemiTheta(DimZMin:DimZMax)      !基本場の温位
    real(8) :: Theta_0                             ! 地表面の温位
    real(8) :: Theta_tr                            ! 対流圏界面の温位
    real(8) :: TakemiZtr                           ! 対流圏界面の高度
    real(8) :: HumSfc                              ! 地表面の相対湿度
    real(8) :: a_MolFrSfc                          ! 地表面のモル比
!!--------------------------------------------------
    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 ECCM_Takemi2007

!!!------------------------------------------------------------------------------!!!
  subroutine HUM_Takemi2007ALL(z_RHall)
    !
    !== 概要
    ! Takemi(2007) での湿度場の与え方の都合上, 高度1.5km以下と高度1.5km以上で
    ! 相対湿度を別々に求めた
    ! 前者は HUM_Takemi2007 or HUM_Takemi2007_TDRY or HUM_Takemi2007_MDRY で求め
    ! 後者は ECCM_Takemi2007 で求めている.
    ! そのため, 全高度での相対湿度 z_RHall を確認のためにも求めるために,
    ! このサブルーチンで両者を結合している
    ! 値の確認が主な目的なため, 基本的に必要ないが, 仮にECCM_MolFrが必要なときに
    ! 相対湿度が必要なこともあり, 一応求めている

    implicit none

    real(8),intent(out) :: z_RHall(DimZMin:DimZMax)
    real(8) :: a_MolFrIn(1:SpcNum)
    real(8) :: Humidit
    real(8) :: z_Tem(DimZMin:DimZMax)
    real(8) :: z_Pre(DimZMin:DimZMax)
    real(8) :: z_MolMea(DimZMin:DimZMax)
    real(8) :: za_MolF(DimZMin:DimZMax)
    real(8) :: ZTR
    real(8) :: z_humidity1(DimZMin:DimZMax) !高度0 - 1.5 km の相対湿度
    real(8) :: z_humidity2(DimZMin:DimZMax) !高度1.5 km 以上の相対湿度
    integer :: k, T, L

    !高度1.5km以下の湿度を取得
    call ECCM_Takemi2007(z_Tem, z_Pre, za_MolF, z_humidity1)

    !高度1.5km以上の湿度を取得
    !三つから選択
    call HUM_Takemi2007(z_humidity2, ZTR, T, L)
!   call HUM_Takemi2007_TDRY(z_humidity2, ZTR, T, L)
!   call HUM_Takemi2007_MDRY(z_humidity2, ZTR, T, L)

    do k = RegZMin+1, DimZMax
      if (k.le.L) then
        z_RHall(k) = z_humidity1(k)

      else if (k.gt.L) then
        z_RHall(k) = z_humidity2(k)

      end if
    end do

  end subroutine HUM_Takemi2007ALL

!!------------------------------------------------------------------------------!!!
  subroutine HUM(z_humid, kansoku_altitude, kansoku_HumZ)
    !
    !== 概要
    !
    !  * Yamasaki, 1983 の湿度の観測データから基本場を作成する 
    !    * 相対湿度の観測データをNetCDFファイル化したものから値を読み込む
    !    * 読み込んだ値を線形補間して相対湿度の基本場を作成する
    !

    implicit none

    real(8) :: Observed_Humidity(1:100,1:88)        !NetCDF から読み込んだ相対湿度の観測値
                                                    !水平一様な2次元データ
    real(8),intent(out) :: kansoku_HumZ(1:88)       !鉛直1次元にした相対湿度の観測値の代入先
    real(8),intent(out) :: z_humid(DimZMin:DimZMax) !線形補間した相対湿度
    real(8),intent(out) :: kansoku_altitude(1:89)   !観測データに対応する高度
    integer :: k,m

    z_humid = 0.0d0
    kansoku_altitude = 0.0d0
    kansoku_HumZ = 0.0d0

    ! 観測データに対応する高度を与える
    ! これは線形補間の際に使用する
    ! なお, Yamsaki(1983)にある元データはデータ間隔がばらばらであったので
    ! 250m間隔に直したものを使用している
    do m = 1,89
     kansoku_altitude(m) = 25.0d1*(m-1)
    end do

    ! NetCDF ファイルから相対湿度を読み込む
    call HistoryGet('Ob_Humidity.nc', 'Ob_Humidity', array=Observed_Humidity)

    ! 相対湿度の観測値のファイル"Ob_Humidity.nc"は2次元データ(x,z)であるため
    ! 1次元データにする
    do m=1,88
      kansoku_HumZ(m)  = Observed_Humidity(1,m)
    end do
    
    ! kansoku_HumZ を線形補間してz_humidを作成する
    ! kansoku_altitude(m) < z_Z(k) < kansoku_altitude(m+1),
    ! kansoku_altitude = z_Z(k),
    ! kansoku_altitude < z_Z(k)
    ! の3つで場合分け
    do k = RegZMin+1,DimZMax
     do m = 1,87
      if (kansoku_altitude(m).ne.z_Z(k).and.z_Z(k).gt.kansoku_altitude(m) &
         & .and.z_Z(k).lt.kansoku_altitude(m+1)) then
        z_humid(k) = kansoku_HumZ(m)+((kansoku_HumZ(m+1) - kansoku_HumZ(m)) &
        & /(kansoku_altitude(m+1)-kansoku_altitude(m)))*(z_Z(k)-kansoku_altitude(m))

      else if (kansoku_altitude(m).eq.z_Z(k)) then
        z_humid(k) = kansoku_HumZ(m)

      else if (kansoku_altitude(m).lt.z_Z(k)) then
        z_humid(k) = z_humid(k-1)
      end if
     end do
    end do

  end subroutine HUM


  subroutine Yamazaki1983
    real(8) :: Temp(1:89) ! 温度の観測データ
    real(8) :: Hum(1:89)  ! 相対湿度の観測データ
    
    Temp = (/ &
      & 300.00, 298.00, 296.00, 294.22, 293.15, 291.08, 289.80, 288.55, 287.30, 286.50, &
      & 284.78, 283.49, 282.16, 280.81, 279.45, 278.07, 276.70, 275.33, 273.95, 272.58, &
      & 271.20, 269.83, 268.45, 267.06, 265.63, 264.19, 262.70, 261.20, 259.68, 258.11, &
      & 256.55, 254.86, 253.18, 251.41, 249.60, 247.76, 245.83, 243.89, 241.90, 240.30, &
      & 237.90, 235.90, 233.90, 231.90, 229.90, 227.90, 225.90, 223.91, 221.97, 220.03, &
      & 218.09, 216.16, 214.24, 212.42, 210.64, 208.86, 207.37, 205.96, 204.55, 203.35, &
      & 202.48, 201.60, 200.73, 200.36, 200.12, 199.89, 199.65, 199.91, 200.30, 200.70, &
      & 201.09, 201.55, 202.30, 203.05, 203.80, 204.55, 205.30, 206.05, 206.80, 207.55, &
      & 208.30, 209.05, 209.69, 210.31, 210.93, 211.54, 212.16, 212.77, 212.99 /)
    
    Hum = (/ &
      & 83.00, 88.85, 89.70, 91.70, 93.70, 94.35, 95.00, 94.35, 93.70, 91.85, 90.00, &
      & 87.50, 85.00, 82.90, 80.80, 78.80, 76.70, 74.35, 72.00, 71.00, 70.00, 68.85, &
      & 67.70, 66.85, 66.00, 65.60, 65.20, 64.75, 64.30, 63.90, 63.50, 63.10, 62.70, &
      & 62.25, 61.80, 61.40, 61.00, 61.00, 61.00, 61.35, 61.70, 61.85, 62.00, 62.15, &
      & 62.30, 62.65, 63.00, 63.00, 63.00, 63.35, 63.70, 64.50, 65.30, 66.60, 67.90, &
      & 70.10, 72.30, 73.30, 74.30, 74.60, 74.90, 68.50, 62.10, 54.05, 46.00, 41.45, &
      & 36.90, 33.25, 29.60, 24.65, 19.70, 15.75, 11.80, 9.00, 6.20, 5.25, 4.30, 3.35, &
      & 2.40, 3.35, 4.30, 2.95, 1.60, 1.05, 0.50, 0.45, 0.40, 0.35 /)


 end subroutine Yamazaki1983


end module Bscset


