!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2005. All rights reserved.
!---------------------------------------------------------------------
!= Subroutine 
!
!   * Developer: SUGIYAMA Ko-ichiro
!   * Version: $Id: warmrainprm.f90,v 1.1.2.1 2006/01/13 06:08:20 sugiyama Exp $
!   * 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 )

  !== 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 differentiate_center4, 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) 
                                            !凝結による水蒸気から雲水への変換. 

  !== 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) :: DMixRtSatDPotTemp_tmp
  real(8) :: DMixRtSatDSvapPress_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

  
!!!
!!!相変化を表現する部分を除いて温位, 混合比の式を解く
!!!
  !=== 時間積分. 時間方向にリープフロッグ, 空間方向に中心差分
  !==== ここでは雲・降水過程は除く
  
  !=== 温位
  ss_PotTemp_al =                                                     &
    &   ss_PotTemp_bl                                                 &
    & + DelTime                                                       &
    &   * (                                                           &
    &     - ss_AdvX_ss_fs(ss_PotTemp_nl, fs_VelX_nl)                  &
    &     - ss_AdvZ_ss_sf(ss_PotTemp_nl, sf_VelZ_nl)                  &
    &     - ss_AdvZ_ss_sf(ss_PotTempBasicZ, sf_VelZ_nl)               &
    &     + ss_TurbXZ_ss_ss(ss_PotTemp_bl, ss_Kh_bl)                  &
    &     + ss_TurbZ_ss_ss(ss_PotTempBasicZ, ss_Kh_bl)                &
    &     + ss_NumDiffXZ_ss(ss_PotTemp_bl)                            &
    &    )
  
  !=== 水蒸気混合比
  ss_MixRtV_al =                                                      &
    &   ss_MixRtV_bl                                                  &  
    & + DelTime                                                       &
    &   * (                                                           &
    &     - ss_AdvX_ss_fs(ss_MixRtV_nl, fs_VelX_nl)                   &
    &     - ss_AdvZ_ss_sf(ss_MixRtV_nl, sf_VelZ_nl)                   &
    &     - ss_AdvZ_ss_sf(ss_MixRtVBasicZ, sf_VelZ_nl)                &
    &     + ss_TurbXZ_ss_ss(ss_MixRtV_bl, ss_Kh_bl)                   & 
    &     + ss_TurbZ_ss_ss(ss_MixRtVBasicZ, ss_Kh_bl)                 & 
    &     + ss_NumDiffXZ_ss(ss_MixRtV_bl)                             &
    &    )

  !=== 雲粒混合比
  ss_MixRtC_al =                                                      &     
    &   ss_MixRtC_bl                                                  &
    & + DelTime                                                       &
    &   * (                                                           &
    &     - ss_AdvX_ss_fs(ss_MixRtC_nl, fs_VelX_nl)                   &
    &     - ss_AdvZ_ss_sf(ss_MixRtC_nl, sf_VelZ_nl)                   &
    &     - ss_AdvXZ_ss_fs_sf(ss_MixRtC_nl, fs_VelX_nl, sf_VelZ_nl)   & 
    &     + ss_TurbXZ_ss_ss  (ss_MixRtC_bl, ss_Kh_bl)                 &
    &     + ss_NumDiffXZ_ss  (ss_MixRtC_bl)                           &
    &    )
  
  !=== 雨粒混合比
  ss_MixRtR_al  =                                                     &
    &   ss_MixRtR_bl                                                  &
    & + DelTime                                                       &
    &   * (                                                           &
    &     - ss_AdvX_ss_fs(ss_MixRtR_nl, fs_VelX_nl)                   &
    &     - ss_AdvZ_ss_sf(ss_MixRtR_nl, sf_VelZ_nl)                   &
    &     + ss_TurbXZ_ss_ss(ss_MixRtR_bl, ss_Kh_bl)                   &
    &     + ss_NumDiffXZ_ss( ss_MixRtR_bl )                           &
    &     )

  !確認
  write(*,*) "MixRtC: ", minval(ss_MixRtC_al), maxval(ss_MixRtC_al), &
    & sum(ss_DensBasicZ(1:RegXMax,1:RegZMax) * ss_MixRtC_al(1:RegXMax,1:RegZMax))
  

  !=== 値が負になったときの処理
  !==== 全量を求める
!  ss_PotTempAll_al = ss_PotTemp_al + ss_PotTempBasicZ
!  ss_MixRtVAll_al  = ss_MixRtV_al  + ss_MixRtVBasicZ

  !==== 負になった部分の穴埋めを行う
!  call FillNegative( ss_PotTempAll_al )
!  call FillNegative( ss_MixRtVAll_al )
  call FillNegative( ss_MixRtC_al )
!  call FillNegative( ss_MixRtR_al )

!  where ( ss_PotTempAll_al <= Min )
!    ss_PotTempAll_al = Min
!  end where
!  where ( ss_MixRtVAll_al <= Min )
!    ss_MixRtVAll_al = Min
!  end where
  
  !==== 擾乱の値に戻す
!  ss_PotTemp_al = ss_PotTempAll_al - ss_PotTempBasicZ
!  ss_MixRtV_al  = ss_MixRtVAll_al  - ss_MixRtVBasicZ
  

  !確認
  write(*,*) "MixRtC: ", minval(ss_MixRtC_al), maxval(ss_MixRtC_al), &
    & sum(ss_DensBasicZ(1:RegXMax,1:RegZMax) * ss_MixRtC_al(1:RegXMax,1:RegZMax))



!!!
!!! 雲・降水過程の考慮. 値が負にならないように気をつける
!!! 
  !=== 通常処理
!  where (ss_MixRtC_al - DelTime * (ss_CNcr_nl + ss_CLcr_nl) >= Min) 
    !=== 雨粒混合比
!    ss_MixRtR_al = ss_MixRtR_al + DelTime * (ss_CNcr_nl + ss_CLcr_nl)

    !=== 雲粒混合比
!    ss_MixRtC_al = ss_MixRtC_al - DelTime * (ss_CNcr_nl + ss_CLcr_nl)  
!  end where

  !=== 例外処理: 雲粒混合比以上に凝縮が生じてしまう場合
!  where (ss_MixRtC_al - DelTime * (ss_CNcr_nl + ss_CLcr_nl) < Min) 
    !=== 雨粒混合比
!    ss_MixRtR_al  = ss_MixRtR_al + ss_MixRtC_al  

    !=== 雲粒混合比
!    ss_MixRtC_al = 0.0d0 
!  end where


  !=== 通常処理
!  where (ss_MixRtR_al - DelTime * (ss_EVrv_nl - ss_PRr_nl) >= Min) 
    !=== 温位
!    ss_PotTemp_al = ss_PotTemp_al - DelTime * ss_Gamma_nl * ss_EVrv_nl   
  
    !=== 水蒸気混合比
!    ss_MixRtV_al = ss_MixRtV_al + DelTime * ss_EVrv_nl   
    
    !=== 雨粒混合比
!    ss_MixRtR_al  = ss_MixRtR_al - DelTime * (ss_EVrv_nl - ss_PRr_nl)
!  end where
  
  !=== 例外処理: 雨粒混合比以上に蒸発が生じる場合
!  where (ss_MixRtR_al - DelTime * (ss_EVrv_nl - ss_PRr_nl) < Min)
    !=== 温位
!    ss_PotTemp_al = ss_PotTemp_al - ss_Gamma_nl * ss_MixRtR_al 
  
    !=== 水蒸気混合比
!    ss_MixRtV_al = ss_MixRtV_al + ss_MixRtR_al 
    
    !=== 雨粒混合比
!    ss_MixRtR_al  = 0.0d0
!  end where


!!!
!!! 湿潤飽和調節法を用いた水蒸気--雲粒間の変換
!!!
  MoistAjastX: do i = DimXMin, DimXMax

    !!!TEST!!!
    !NO CLOUD
    cycle MoistAjastX

    MoistAjastZ: do k = DimZMin, DimZMax
      
      !=== 今までに得られた擾乱成分の値を暫定量とみなす. 添え字を変更
      MixRtV_tmp   = ss_MixRtV_al(i,k)
      MixRtC_tmp   = ss_MixRtC_al(i,k)
      PotTemp_tmp  = ss_PotTemp_al(i,k)
      
      !=== 湿潤飽和法では圧力は変化させない. 
      !==== 湿潤飽和調節によって水蒸気モル比が変化すれば, 
      !==== 比熱と気体定数の値が変化し, エクスナー関数の値は変化するはず. 
      ExnerAll_adj = ss_ExnerAll_nl(i,k)
      PressAll_adj = ss_PressAll_nl(i,k)
      
      !=== ループ回数の初期化
      num = 0
      
      Adjasting: do 
        !=== 擾乱成分と平均成分の和を求める        
        PotTempAll_tmp = PotTemp_tmp    + ss_PotTempBasicZ(i,k) 
        MixRtVAll_tmp  = MixRtV_tmp     + ss_MixRtVBasicZ(i,k) 
        TempAll_tmp    = PotTempAll_tmp * ExnerAll_adj
        
        !==== 飽和蒸気圧
        SvapPress_tmp = s_svap_s(TempAll_tmp)
  
        !==== 飽和蒸気圧から飽和混合比を計算
        MixRtVSat_tmp = &
          & ((MolWtWet/MolWtDry) * SvapPress_tmp) &
          & / (PressSfc * (ExnerAll_adj**(CpDry / GasRDry)) - SvapPress_tmp)
        
        !=== 変数の初期化
        Gamma_tmp = s_LatentHeat_s(TempAll_tmp) / (ExnerAll_adj * CpDry)
        
        !=== 水蒸気分圧が飽和蒸気圧より多い or 雲粒混合比が正
        if (MixRtVAll_tmp - MixRtVSat_tmp > Min  .OR.        &
          & MixRtC_tmp > Min) then 

          DMixRtSatDSvapPress_tmp =  &
            & PressSfc                              &
            & * (MolWtWet / MolWtDry)               &
            & * (ExnerAll_adj ** (CpDry / GasRDry)) &
            & / ((PressSfc * (ExnerAll_adj ** (CpDry / GasRDry)) &
            &     - SvapPress_tmp) ** 2.0d0)
          
          DMixRtSatDPotTemp_tmp = &
            & DMixRtSatDSvapPress_tmp * s_DSvapDT_s(TempAll_tmp) * ExnerAll_adj
          
          DelMixRtV_adj =  &
            & (MixRtVAll_tmp - MixRtVSat_tmp) &
            &  / ( 1.0d0 + Gamma_tmp * DMixRtSatDPotTemp_tmp )          
          
          PotTemp_adj = PotTemp_tmp + Gamma_tmp * DelMixRtV_adj 
          MixRtV_adj  = MixRtV_tmp - DelMixRtV_adj             
          MixRtC_adj  = MixRtC_tmp + DelMixRtV_adj            

          !=== 雲粒混合比が収束していたら繰り返しを中止する
          if (abs((MixRtC_adj - MixRtC_tmp) / MixRtC_adj) <= Min) then 
            
            ss_PotTemp_al(i,k) = PotTemp_adj                 
            ss_MixRtV_al(i,k)  = MixRtV_adj
            ss_MixRtC_al(i,k)  = MixRtC_adj
            
            !=== 雲粒混合比が負になったときは, 雲粒混合比を正にする. 
!            if ( MixRtC_adj <= Min) then 
!              ss_PotTemp_al(i,k) = PotTemp_tmp - Gamma_tmp * MixRtC_tmp
!              ss_MixRtV_al(i,k)  = MixRtV_tmp + MixRtC_tmp 
!              ss_MixRtC_al(i,k)  = 0.0d0
!            end if

            exit Adjasting            
          end if
          
          !=== ループ回数が回りすぎたらエラーを吐いて続ける
         if ( num > 20 ) then 
           write(*,*) "Looping is too much, num > 20"
            write(*,*) "ss_PotTemp_al(i,k) =     ", PotTemp_adj
            write(*,*) "ss_MixRtV_al(i,k) = ", MixRtV_adj
            write(*,*) "ss_MixRtC_al(i,k) = ", MixRtC_adj
            ss_PotTemp_al(i,k) = PotTemp_adj     
            ss_MixRtV_al(i,k)  = MixRtV_adj
            ss_MixRtC_al(i,k)  = MixRtC_adj
            exit Adjasting
          end if

        else 

          !=== 収束したとして終了
          ss_PotTemp_al(i,k) = PotTemp_tmp - Gamma_tmp * MixRtC_tmp
          ss_MixRtV_al(i,k) = MixRtV_tmp + MixRtC_tmp
          ss_MixRtC_al(i,k) = 0.0d0          
          exit Adjasting
        end if

         !---- 繰り返しのための変数定義
        PotTemp_tmp = PotTemp_adj
        MixRtV_tmp  = MixRtV_adj
        MixRtC_tmp  = MixRtC_adj
        num         = num + 1
       
      end do Adjasting
      
    end do MoistAjastZ
  end do MoistAjastX
  

  !==== 飽和蒸気圧のモル比
  ss_SvapPress_nl = &
    & ss_svap_ss(ss_ExnerAll_nl  * (ss_PotTemp_al(i,k) + ss_PotTempBasicZ(i,k))) 

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


  !=== 境界条件
  call boundary(ss_BC, ss_PotTemp_al)
  call boundary(ss_BC, ss_MixRtV_al)
  call boundary(ss_BC, ss_MixRtC_al)
  call boundary(ss_BC, ss_MixRtR_al)  


  
  call EndSub("Scalar")

  
contains

  function ss_AdvXZ_ss_fs_sf(ss_Scalar, fs_VelX, sf_VelZ)
    use gridset, only: DimXMin, DimXMax, DimZMin, DimZMax
    use average,  only: ss_avr_fs, ss_avr_sf
    use differentiate_center4, only: fs_dx_ss, sf_dz_ss
    implicit none
    
    real(8), intent(in) :: fs_VelX(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8), intent(in) :: sf_VelZ(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8), intent(in) :: ss_Scalar(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)             :: ss_AdvXZ_ss_fs_sf(DimXMin:DimXMax, DimZMin:DimZMax)
    
    ss_AdvXZ_ss_fs_sf = 0.0d0  !初期化
    ss_AdvXZ_ss_fs_sf = &
      &   ss_avr_fs(fs_VelX * fs_dx_ss(ss_Scalar)) &
      & + ss_avr_sf(sf_VelZ * sf_dz_ss(ss_Scalar))
    
  end function ss_AdvXZ_ss_fs_sf


  function ss_AdvZ_ss_sf(ss_Scalar, sf_VelZ)
    use gridset, only: DimXMin, DimXMax, DimZMin, DimZMax
    use average,  only: ss_avr_fs, ss_avr_sf
    use differentiate_center4, only: sf_dz_ss
    implicit none
    
    real(8), intent(in) :: sf_VelZ(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8), intent(in) :: ss_Scalar(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)             :: ss_AdvZ_ss_sf(DimXMin:DimXMax, DimZMin:DimZMax)
    
    ss_AdvZ_ss_sf = 0.0d0  !初期化
    ss_AdvZ_ss_sf = ss_avr_sf(sf_VelZ * sf_dz_ss(ss_Scalar))

!    !!!!test!!!!
!    ss_AdvZ_ss_sf(:,DimZMin:1) = 0.0d0
    
  end function ss_AdvZ_ss_sf
  
  
  function ss_AdvX_ss_fs(ss_Scalar, fs_VelX)
    use gridset, only: DimXMin, DimXMax, DimZMin, DimZMax
    use average,  only: ss_avr_fs, ss_avr_fs
    use differentiate_center4, only: fs_dx_ss
    implicit none
    
    real(8), intent(in) :: fs_VelX(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8), intent(in) :: ss_Scalar(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)             :: ss_AdvX_ss_fs(DimXMin:DimXMax, DimZMin:DimZMax)
    
    ss_AdvX_ss_fs = 0.0d0  !初期化
    ss_AdvX_ss_fs = ss_avr_fs(fs_VelX * fs_dx_ss(ss_Scalar))
    
  end function ss_AdvX_ss_fs
  
  
  function ss_TurbXZ_ss_ss(ss_Scalar, ss_Kh)
    use gridset, only: DimXMin, DimXMax, DimZMin, DimZMax
    use average,  only: fs_avr_ss, sf_avr_ss
    use differentiate_center4, only: ss_dx_fs, fs_dx_ss, ss_dz_sf, sf_dz_ss
    implicit none
    
    real(8), intent(in) :: ss_Kh(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8), intent(in) :: ss_Scalar(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)             :: ss_TurbXZ_ss_ss(DimXMin:DimXMax, DimZMin:DimZMax)

    ss_TurbXZ_ss_ss = 0.0d0  !初期化
    ss_TurbXZ_ss_ss = &
      &   ss_dx_fs(fs_avr_ss(ss_Kh) * fs_dx_ss(ss_Scalar)) &
      & + ss_dz_sf(sf_avr_ss(ss_Kh) * sf_dz_ss(ss_Scalar))
    
  end function ss_TurbXZ_ss_ss


  function ss_TurbZ_ss_ss(ss_Scalar, ss_Kh)
    use gridset, only: DimXMin, DimXMax, DimZMin, DimZMax
    use average,  only: fs_avr_ss, sf_avr_ss
    use differentiate_center4, only: ss_dx_fs, fs_dx_ss, ss_dz_sf, sf_dz_ss
    implicit none

    real(8), intent(in) :: ss_Kh(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8), intent(in) :: ss_Scalar(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)             :: ss_TurbZ_ss_ss(DimXMin:DimXMax, DimZMin:DimZMax)

    ss_TurbZ_ss_ss = 0.0d0  !初期化
    ss_TurbZ_ss_ss =  ss_dz_sf(sf_avr_ss(ss_Kh) * sf_dz_ss(ss_Scalar))
    
  end function ss_TurbZ_ss_ss


  function ss_NumDiffXZ_ss(ss_Scalar)
    use gridset, only: DimXMin, DimXMax, DimZMin, DimZMax
    use differentiate_center4, only: ss_dx_fs, fs_dx_ss, ss_dz_sf, sf_dz_ss
    use arareset, only: NuH, NuV
    implicit none
    
    real(8), intent(in) :: ss_Scalar(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)             :: ss_NumDiffXZ_ss(DimXMin:DimXMax, DimZMin:DimZMax)
    
    ss_NumDiffXZ_ss = 0.0d0  !初期化
    ss_NumDiffXZ_ss = &
      &   NuH * (ss_dx_fs(fs_dx_ss( ss_Scalar ))) &
      & + NuV * (ss_dz_sf(sf_dz_ss( ss_Scalar ))) 
    
  end function ss_NumDiffXZ_ss


  function s_DSvapDT_s(s_temp)
    ! 暗黙の型宣言禁止
    implicit none

    ! 入出力変数
    real(8)             :: s_DSvapDT_s          !飽和蒸気圧 [Pa]
    real(8), intent(in) :: s_temp               !温度 [K]
    
    ! Antoine の式を利用. H2O の係数は化学工学便覧参照. 
    real(8), parameter  :: A = 23.1964d0
    real(8), parameter  :: B = 3816.44d0
    real(8), parameter  :: C = -46.13d0
    
    s_DSvapDT_s = dexp( A - B / (C + s_temp) ) * (B / ((C + s_temp) ** 2.0d0))
    
  end function s_DSvapDT_s

  
  function s_svap_s(s_temp)
    ! 暗黙の型宣言禁止
    implicit none

    ! 入出力変数
    real(8)             :: s_svap_s             !飽和蒸気圧 [Pa]
    real(8), intent(in) :: s_temp               !温度 [K]
    
    ! Antoine の式を利用. H2O の係数は化学工学便覧参照. 
    real(8), parameter  :: A = 23.1964d0
    real(8), parameter  :: B = 3816.44d0
    real(8), parameter  :: C = -46.13d0
    
    s_svap_s = dexp( A - B / (C + s_temp) ) !Antoine の式
    
  end function s_svap_s

  
  function ss_svap_ss(ss_temp)
    ! 暗黙の型宣言禁止
    implicit none

    ! 入出力変数
    real(8)             :: ss_svap_ss(DimXMin:DimXMax, DimZMin:DimZMax)
                                                        !飽和蒸気圧 [Pa]
    real(8), intent(in) :: ss_temp(DimXMin:DimXMax, DimZMin:DimZMax) 
                                                               !温度 [K]
    
    ! Antoine の式を利用. H2O の係数は化学工学便覧参照. 
    real(8), parameter  :: A = 23.1964d0
    real(8), parameter  :: B = 3816.44d0
    real(8), parameter  :: C = -46.13d0
    
    ss_svap_ss = dexp( A - B / (C + ss_temp) ) !Antoine の式

  end function ss_svap_ss

  
  function s_LatentHeat_s(temp)
    use physset, only: GasRWet
    implicit none

    real(8), intent(in)    :: temp         !温度
    real(8)                :: s_LatentHeat_s   !潜熱
    real(8)                :: svap_log_dt

    ! Antoine の式を利用. H2O の係数は化学工学便覧参照. 
    real(8), parameter  :: A = 23.1964d0
    real(8), parameter  :: B = 3816.44d0
    real(8), parameter  :: C = -46.13d0

    svap_log_dt = B / ((C + temp) ** 2.0d0)
    s_LatentHeat_s  = svap_log_dt * GasRWet * ( temp ** 2.0d0 )
    
  end function s_LatentHeat_s
  
  function ss_LatentHeat_ss(ss_temp)
    use gridset, only: DimXMin, DimXMax, DimZMin, DimZMax
    use physset, only: GasRWet
    implicit none
    
    real(8), intent(in) :: ss_temp(DimXMin:DimXMax, DimZMin:DimZMax) !温度
    real(8)             :: ss_LatentHeat_ss(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)             :: svap_log_dt(DimXMin:DimXMax, DimZMin:DimZMax)
    
    ! Antoine の式を利用. H2O の係数は化学工学便覧参照. 
!    real(8), parameter  :: A = 23.1964d0
    real(8), parameter  :: B = 3816.44d0
    real(8), parameter  :: C = -46.13d0

    svap_log_dt = B / ((C + ss_temp)** 2.0d0)
    ss_LatentHeat_ss  = svap_log_dt * GasRWet * ( ss_Temp ** 2.0d0 ) 
    
  end function ss_LatentHeat_ss


  
    
  subroutine warmrainprm_init
    !=== 出力変数
    ss_PotTemp_al = 0.0d0
    ss_MixRtV_al = 0.0d0
    ss_MixRtC_al = 0.0d0
    ss_MixRtR_al = 0.0d0
    ss_EVrv_nl = 0.0d0
    ss_CNcr_nl = 0.0d0
    ss_CLcr_nl = 0.0d0
    ss_VelZRain_nl = 0.0d0
    ss_PRr_nl = 0.0d0

    !=== 内部変数
    ss_Gamma_nl = 0.0d0
    ss_PressAll_nl = 0.0d0
    ss_ExnerAll_nl = 0.0d0
    ss_TempAll_nl = 0.0d0
    ss_PotTempAll_nl = 0.0d0
    ss_MixRtVAll_nl = 0.0d0
    ss_MixRtVSat_nl = 0.0d0
    ss_PotTempAll_al = 0.0d0
    ss_SvapPress_nl = 0.0d0

    MixRtV_tmp = 0.0d0
    MixRtC_tmp = 0.0d0
    PotTemp_tmp = 0.0d0
    
    PressAll_adj = 0.0d0
    ExnerAll_adj = 0.0d0
    
    TempAll_tmp = 0.0d0
    PotTempAll_tmp = 0.0d0
    MixRtVAll_tmp = 0.0d0
    
    SvapPress_tmp = 0.0d0
    MixRtVSat_tmp = 0.0d0
    Gamma_tmp = 0.0d0
    
    PotTemp_adj = 0.0d0
    MixRtV_adj = 0.0d0
    MixRtC_adj = 0.0d0
    DensGround = 0.0d0

  end subroutine warmrainprm_init
  

end subroutine WarmRainPrm
