!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2005. All rights reserved.
!---------------------------------------------------------------------
!= Subroutine
!
!   * Developer: SUGIYAMA Ko-ichiro
!   * Version: $Id: $
!   * Tag Name: $Name:  $
!   * Change History:
!
!== Overview
!
!暖かい雨のバルク法を用いた, 水蒸気と雨, 雲と雨の混合比の変換係数を求める.
!
!== Error Handling
!
!== Bugs
!
!== Note
!
!== Future Plans
!
!
subroutine WarmRainPrm( &
  & DelTime, &
  & fs_VelX_nl, sf_VelZ_nl, ss_Exner_nl, &
  & ss_PotTemp_nl, ss_PotTemp_bl, &
  & ss_MixRtV_nl, ss_MixRtV_bl, &
  & ss_MixRtC_nl, ss_MixRtC_bl, &
  & ss_MixRtR_nl,  ss_MixRtR_bl, &
  & ss_Kh_bl, &
  & ss_PotTemp_al, &
  & ss_MixRtV_al, &
  & ss_MixRtC_al, &
  & ss_MixRtR_al, &
  & ss_EVrv_nl, &
  & ss_CNcr_nl,&
  & ss_CLcr_nl,&
  & ss_VelZRain_nl,&
  & ss_PRr_nl, &
  & ss_MixRtVSat_al    )

  !== Dependency
  use dc_trace, only: BeginSub, EndSub
  use gridset,  only: DimXMin, DimXMax, DimZMin, DimZMax, &
    &                 RegXMin, RegXMax, RegZMin, RegZMax
  use bcset,    only: ss_BC
  use average,  only: ss_avr_fs, ss_avr_sf, fs_avr_ss, sf_avr_ss
  use basicset, only: ss_ExnerBasicZ, &
       &              ss_DensBasicZ,  &
       &              ss_PotTempBasicZ, &
       &              ss_MixRtVBasicZ,  &
       &              ss_PressBasicZ
  use physset,  only: PressSfc, GasRUniv, GasRDry, MolWtDry, MolWtWet, CpDry, epsilon
  use arareset, only: Min
  use satvappress, only: ss_Svap_ss, Svap, ss_LatentHeat_ss, LatentHeat, DSvapDTemp
  use differentiate_center2, only: sf_dz_ss

  !--- 暗黙の型宣言禁止
  implicit none

  !== Input
  real(8), intent(in) :: DelTime !積分時間
  real(8), intent(in) :: ss_Exner_nl(DimXMin:DimXMax, DimZMin:DimZMax)
                                                          !エクスナー関数
  real(8), intent(in) :: ss_PotTemp_nl(DimXMin:DimXMax, DimZMin:DimZMax)
                                                          !温位
  real(8), intent(in) :: ss_PotTemp_bl(DimXMin:DimXMax, DimZMin:DimZMax)
                                                          !温位
  real(8), intent(in) :: ss_MixRtV_nl(DimXMin:DimXMax, DimZMin:DimZMax)
                                                          !水蒸気量
  real(8), intent(in) :: ss_MixRtV_bl(DimXMin:DimXMax, DimZMin:DimZMax)
                                                          !水蒸気量
  real(8), intent(in) :: ss_MixRtC_nl(DimXMin:DimXMax, DimZMin:DimZMax)
                                                          !雲粒混合比
  real(8), intent(in) :: ss_MixRtC_bl(DimXMin:DimXMax, DimZMin:DimZMax)
                                                          !雲粒混合比
  real(8), intent(in) :: ss_MixRtR_nl(DimXMin:DimXMax, DimZMin:DimZMax)
                                                          !雨粒混合比
  real(8), intent(in) :: ss_MixRtR_bl(DimXMin:DimXMax, DimZMin:DimZMax)
                                                          !雨粒混合比
  real(8), intent(in) :: fs_VelX_nl(DimXMin:DimXMax, DimZMin:DimZMax)
                                                          !速度 u
  real(8), intent(in) :: sf_VelZ_nl(DimXMin:DimXMax, DimZMin:DimZMax)
                                                          !速度 w
  real(8), intent(in) :: ss_Kh_bl(DimXMin:DimXMax, DimZMin:DimZMax)
                                                          !渦粘性係数

  !== Out
  real(8), intent(Out)  :: ss_PotTemp_al(DimXMin:DimXMax, DimZMin:DimZMax)
                                            !温位
  real(8), intent(Out)  :: ss_MixRtV_al(DimXMin:DimXMax, DimZMin:DimZMax)
                                            !水蒸気混合比
  real(8), intent(Out)  :: ss_MixRtC_al(DimXMin:DimXMax, DimZMin:DimZMax)
                                            !雲粒混合比
  real(8), intent(Out)  :: ss_MixRtR_al(DimXMin:DimXMax, DimZMin:DimZMax)
                                            !雨粒混合比
  real(8), intent(Out)  :: ss_EVrv_nl(DimXMin:DimXMax, DimZMin:DimZMax)
                                            !蒸発による雨水から水蒸気への変換
  real(8), intent(Out)  :: ss_CNcr_nl(DimXMin:DimXMax, DimZMin:DimZMax)
                                            !併合成長による雲水から雨水への変換
  real(8), intent(Out)  :: ss_CLcr_nl(DimXMin:DimXMax, DimZMin:DimZMax)
                                            !凝結による水蒸気から雲水への変換
  real(8), intent(Out)  :: ss_VelZRain_nl(DimXMin:DimXMax, DimZMin:DimZMax)
                                            !雨粒終端速度
  real(8), intent(Out)  :: ss_PRr_nl(DimXMin:DimXMax, DimZMin:DimZMax)
                                            !凝結による水蒸気から雲水への変換.
  real(8), intent(Out)  :: ss_MixRtVSat_al(DimXMin:DimXMax, DimZMin:DimZMax)
                                            !飽和蒸気圧

  !== Work
  real(8) :: ss_Gamma_nl(DimXMin:DimXMax, DimZMin:DimZMax)
                                             ! 潜熱 / ( 比熱 * エクスナー関数 )
  real(8) :: ss_PressAll_nl(DimXMin:DimXMax, DimZMin:DimZMax)
                                             ! 圧力(擾乱+平均)
  real(8) :: ss_ExnerAll_nl(DimXMin:DimXMax, DimZMin:DimZMax)
                                             ! 無次元圧力(擾乱+平均)
  real(8) :: ss_TempAll_nl(DimXMin:DimXMax, DimZMin:DimZMax)
                                             ! 温度(擾乱+平均)
  real(8) :: ss_PotTempAll_nl(DimXMin:DimXMax, DimZMin:DimZMax)
                                             ! 温位(擾乱+平均)
  real(8) :: ss_MixRtVAll_nl(DimXMin:DimXMax, DimZMin:DimZMax)
                                             ! 水蒸気混合比(擾乱+平均)
  real(8) :: ss_MixRtVSat_nl(DimXMin:DimXMax, DimZMin:DimZMax)
                                             ! 飽和蒸気圧混合比
  real(8) :: ss_SvapPress_nl(DimXMin:DimXMax, DimZMin:DimZMax)
                                             ! 飽和蒸気圧モル比
  real(8) :: ss_PotTempAll_al(DimXMin:DimXMax, DimZMin:DimZMax)
                                             ! 全温位
  real(8) :: ss_MixRtVAll_al(DimXMin:DimXMax, DimZMin:DimZMax)
                                             ! 全水蒸気

  real(8) :: MixRtV_tmp
  real(8) :: MixRtC_tmp
  real(8) :: PotTemp_tmp

  real(8) :: PressAll_adj
  real(8) :: ExnerAll_adj
  real(8) :: DelMixRtV_adj

  real(8) :: TempAll_tmp
  real(8) :: PotTempAll_tmp
  real(8) :: MixRtVAll_tmp

  real(8) :: SvapPress_tmp
  real(8) :: MixRtVSat_tmp
  real(8) :: Gamma_tmp

  real(8) :: PotTemp_adj
  real(8) :: MixRtV_adj
  real(8) :: MixRtC_adj

  real(8) :: DMixRtVSatDPotTemp_tmp
  real(8) :: DMixRtVSatDSvapPress_tmp

  real(8) :: DensGround         ! 地表面密度
  integer :: i, k, num          ! 添え字
  real(8), parameter :: K1 = 1.0d-3  !ケスラーの係数
  real(8), parameter :: K2 = 2.2d0   !ケスラーの係数
  real(8), parameter :: A  = 1.0d-3  !雲水量の閾値
  real(8), parameter :: N0 = 5.0d7
  real(8), parameter :: D0 = 3.66d-1

  !=== サブルーチンの初期化
  call BeginSub("Scalar", &
    &           fmt="%c",              &
    &           c1="Time integration of scalar component.")

  !=== 変数に零を代入
  call warmrainprm_init

  !=== 温度, 圧力, 混合比の全量を求める
  !==== 擾乱成分と平均成分の足し算
  ss_PotTempAll_nl = ss_PotTemp_nl  + ss_PotTempBasicZ
  ss_ExnerAll_nl   = ss_Exner_nl    + ss_ExnerBasicZ
  ss_MixRtVAll_nl  = ss_MixRtV_nl   + ss_MixRtVBasicZ
  ss_TempAll_nl    = ss_ExnerAll_nl * ss_PotTempAll_nl

  !=== 飽和蒸気圧に伴う物理量の計算
  !==== 飽和蒸気圧のモル比
  ss_SvapPress_nl = ss_Svap_ss(ss_TempAll_nl)

  !==== 飽和蒸気圧から飽和混合比を計算
  ss_MixRtVSat_nl = &
    & ((MolWtWet/MolWtDry) * ss_SvapPress_nl) &
    & / (PressSfc * (ss_ExnerAll_nl **(CpDry / GasRDry)) - ss_SvapPress_nl)

  !=== 変数の初期化
  ss_Gamma_nl   = &
    & ss_LatentHeat_ss(ss_TempAll_nl) / (ss_ExnerAll_nl * CpDry)
  DensGround = ss_DensBasicZ(RegXMin, 0)   ! X 座標はどこでも同じ

!!!
!!! 計算に必要となる, 水蒸気--雲粒--雨粒間の変換係数を求める
!!!
  !=== 併合成長による雲水から雨水への変換
  !==== Berry (1968) のパラメタリゼーション
!  where (ss_MixRtC_nl > A)
!    ss_CNcr_nl = K1 * (ss_MixRtC_nl - A)
!  end where
  where (ss_MixRtC_nl > Min)
    ss_CNcr_nl =  &
      & ss_DensBasicZ * (ss_MixRtC_nl ** 3.0d0) * 1.0d6 &
      & / (60.0d0 * 2.0d0 * ss_MixRtC_nl                &
      &    + 60.0d0 * 2.66d-8 * N0 / (ss_DensBasicZ * D0))
  end where

  !=== 凝結による水蒸気から雲水への変換
  !==== Kessler (1969) のパラメタリゼーション
  where (ss_MixRtC_nl > Min .AND. ss_MixRtR_nl > Min)
    ss_CLcr_nl = K2 * ss_MixRtC_nl * (ss_MixRtR_nl ** 0.875d0)
  end where

  !=== 蒸発による雨水から水蒸気への変換
  !====== 中島健介 (1994) で利用した定式.
  where (ss_MixRtVSat_nl > ss_MixRtVAll_nl  .AND. ss_MixRtR_nl > Min )
    ss_EVrv_nl = &
      &  4.85d-2                                        &
      &  * (ss_MixRtVSat_nl - ss_MixRtVAll_nl )         &
      &  * ((ss_DensBasicZ * ss_MixRtR_nl) ** 0.65d0 )
  end where

  !=== 凝結による水蒸気から雲水への変換. 終端速度を計算してから評価.
  !==== 雨粒終端速度
  where (ss_MixRtR_nl > Min)
    ss_VelZRain_nl = 12.2d0 * (ss_MixRtR_nl ** 0.125d0)
  end where

  !==== 凝結による水蒸気から雲水への変換.
  !===== Dens の avr を取ってから割ると, ゼロ割が生じるので注意
  where (ss_MixRtR_nl > Min)
    ss_PRr_nl &
      & = ss_avr_sf(  &
      &     sf_dz_ss(ss_DensBasicZ * ss_VelZRain_nl * ss_MixRtR_nl) &
      &   ) / ss_DensBasicZ
  end where

end subroutine WarmRainPrm