!= Subroutine WarmRainPrm
!
! Authors::   SUGIYAMA Ko-ichiro, ODAKA Masatsugu
! Version::   $Id: warmrainprm.f90,v 1.16 2006/09/30 05:53:54 odakker Exp $
! Tag Name::  $Name: arare4-20061019 $
! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview
!
!暖かい雨のバルク法を用いた, 水蒸気と雨, 雲と雨の混合比の変換係数を求める.
!   * 中島健介 (1994) で利用した定式をそのまま利用. 
! 
!== Error Handling
!
!== Bugs
!
!== Note
!
!== Future Plans
!
!

module WarmRainPrm
  !
  !暖かい雨のバルク法を用いた, 水蒸気と雨, 雲と雨の混合比の変換係数を求める.
  !   * 中島健介 (1994) で利用した定式をそのまま利用. 
  ! 
  
  !モジュール読み込み
  use dc_message, only: MessageNotify

  use gridset, only:  DimXMin,           &! x 方向の配列の下限
    &                 DimXMax,           &! x 方向の配列の上限
    &                 DimZMin,           &! z 方向の配列の下限
    &                 DimZMax,           &! z 方向の配列の上限
    &                 SpcNum              !化学種の数
  use basicset, only: PressBasis,        &!温位の基準圧力 
    &                 CpDry,             &!乾燥成分の比熱
    &                 MolWtWet,          &!
    &                 MolWtDry,          &!
    &                 SpcWetID,          &!
    &                 SpcWetSymbol,      &!
    &                 xz_DensBasicZ,     &!基本場の密度
    &                 xz_PotTempBasicZ,  &!基本場の温位
    &                 xz_ExnerBasicZ,    &!基本場の無次元圧力
    &                 xza_MixRtBasicZ,   &!基本場の混合比
    &                 GasRDry             !乾燥成分の気体定数 
  use average,  only: xz_avr_xr        
  use differentiate_center2, only: xr_dz_xz
  use ChemData, only: ChemData_OneSpcID
  use ChemCalc, only: xz_SvapPress, xz_LatentHeat, ReactHeatNH4SH
  use MoistFunc,only: xz_DelMixRtNH4SH
  use StoreSet, only: StoreCond
  
  !暗黙の型宣言禁止
  implicit none

  !属性の指定
  private

  !関数を public にする
  public WarmRainPrm_Init
  public WarmRainPrm_prm
  public xz_Rain2GasHeat
  public xz_Rain2GasHeatNH4SH
  public xza_Rain2Gas
  public xza_Rain2GasNH4SH
  public xza_Cloud2Rain
  public xz_FallRain

  real(8)      :: FactorJ      = 1.0d0 !雲物理過程のパラメータ
                                       !木星では 3.0d0
                                       !地球では 1.0d0 とする
  real(8)      :: AutoConvTime = 1.0d3 !併合成長の時定数 [sec]
  real(8)      :: MixRt_AutoConvCr = 1.0d-3 
                                       !併合成長を生じる臨界混合比 [kg/kg]
  integer      :: LoopNum      = 0
  integer      :: LoopNum2     = 0
  integer      :: GasNum(10)   = 0
  integer      :: CloudNum(10) = 0
  integer      :: RainNum(10)  = 0
  integer      :: NH3Num   = 0
  integer      :: H2SNum   = 0
  integer      :: NH4SHCloudNum = 0
  integer      :: NH4SHRainNum = 0
  real(8), allocatable :: RainSW(:)

  save FactorJ, AutoConvTime, MixRt_AutoConvCr
  save LoopNum, LoopNum2, RainNum, CloudNum, GasNum
  save NH3Num, H2SNum, NH4SHCloudNum, NH4SHRainNum
  save RainSW

contains  

!!!=================================================================================!!!
  subroutine WarmRainPrm_Init( cfgfile )

    !暗黙の型宣言禁止
    implicit none

    !変数定義
    integer                  :: s
    integer                  :: n1, n2    

    !入力変数
    character(*), intent(in) :: cfgfile

    !-----------------------------------------------------------
    ! NAMELIST から情報を取得
    !-----------------------------------------------------------
    !
    NAMELIST /warmrainprm/                              &
      & FactorJ, AutoConvTime, MixRt_AutoConvCr

    open (10, FILE=cfgfile)
    read(10, NML=warmrainprm)
    close(10)

    !-----------------------------------------------------------
    ! 雨粒と雲粒と気体の ID の組を作る
    !-----------------------------------------------------------
    !初期化
    allocate( RainSW(SpcNum) )
    LoopNum  = 0
    LoopNum2 = 0
    RainSW   = 0.0d0

    !化学種の中から雨粒を作るものを選び, その配列添え字と分子量を保管.
    SelectCloud: do s = 1, SpcNum

      !'Cloud' という文字列が含まれるものの個数を数える
      n1 = index(SpcWetSymbol(s), '-Cloud' )
      if (n1 /= 0) then
        LoopNum           = LoopNum + 1
        GasNum(LoopNum)   = minloc(SpcWetID, 1, SpcWetID == ChemData_OneSpcID(SpcWetSymbol(s)(1:n1-3) // '-g'))
        CloudNum(LoopNum) = s
      end if

      !'Rain' という文字列が含まれるものの個数を数える
      n2 = index(SpcWetSymbol(s), '-Rain' )
      if (n2 /= 0) then
        LoopNum2          = LoopNum2 + 1
        RainNum(LoopNum2) = s
        RainSW(s)         = 1.0d0
      end if

      ! NH4SH が存在する場合は LoopNum を 1 つ減らす
      if ( trim(SpcWetSymbol(s)) == 'NH4SH-s-Cloud' ) then 
        LoopNum = LoopNum - 1
      end if

    end do SelectCloud
    
    !-----------------------------------------------------------
    ! 硫化アンモニウム, およびアンモニアと硫化水素の ID を取得
    !   'Cloud' の方が 'Rain' よりも先に存在すると仮定している
    !-----------------------------------------------------------
    NH3Num   = minloc(SpcWetID, 1, SpcWetID == ChemData_OneSpcID('NH3-g'))
    H2SNum   = minloc(SpcWetID, 1, SpcWetID == ChemData_OneSpcID('H2S-g'))
    NH4SHCloudNum = minloc(SpcWetID, 1, SpcWetID == ChemData_OneSpcID('NH4SH-s'))
    NH4SHRainNum  = NH4SHCloudNum + 1

    !-----------------------------------------------------------
    ! 確認
    !-----------------------------------------------------------
    if ( LoopNum == 0 ) then 
      call MessageNotify( "W", &
        & "WarmRainPrm", "CloudNum = 0, please comment out of WarmRainPrm" )
!      stop
    end if

    call MessageNotify( "M", &
      &  "WarmRainPrm_Init", "LoopNum = %d",  i=(/LoopNum/) )
    call MessageNotify( "M", &
      & "WarmRainPrm_Init", "LoopNum2 = %d", i=(/LoopNum2/) )    
    call MessageNotify( "M", &
      & "WarmRainPrm_Init", "GasNum = %d",   i=(/GasNum/)   ) 
    call MessageNotify( "M", &
      & "WarmRainPrm_Init", "CloudNum = %d", i=(/CloudNum/) )
    call MessageNotify( "M", &
      & "WarmRainPrm_Init", "RainNum = %d",  i=(/RainNum/)  ) 
    call MessageNotify( "M", &
      & "WarmRainPrm_Init", "RainSW = %f",   d=(/RainSW/)   )
    call MessageNotify( "M", &
      & "WarmRainPrm_Init", "NH3Num = %d",   i=(/NH3Num/)   )
    call MessageNotify( "M", &
      & "WarmRainPrm_Init", "H2SNum = %d",   i=(/H2SNum/)   )
    call MessageNotify( "M", &
      & "WarmRainPrm_Init", "NH4SHCloudNum = %d", i=(/NH4SHCloudNum/) )
    call MessageNotify( "M", &
      & "WarmRainPrm_Init", "NH4SHRainNum = %d", i=(/NH4SHRainNum/) )

  end subroutine WarmRainPrm_Init


  subroutine WarmRainPrm_prm( M1, MA, MB, M2, M3, M4, M5, M6 )
    
    implicit none
    integer, intent(out) ::  M1, MA, MB(10), M2(10), M3(10), M4, M5, M6

    M1 = LoopNum
    MA = LoopNum2
    MB = RainNum
    M2 = CloudNum
    M3 = GasNum    
    M4 = NH3Num
    M5 = H2SNum
    M6 = NH4SHCloudNum

!    write(*,*) M1, MA, MB, M2, M3, M4, M5, M6 
  end subroutine WarmRainPrm_prm

!!!=================================================================================!!!  
  function xza_Rain2Gas(xz_Exner, xz_PotTemp, xza_MixRt, DelTime)
    !
    ! 雨粒から蒸気への変換量を計算するためのルーチン
    !
    ! 変換量および, 蒸気と雨粒の混合比は正の量なので, 計算の途中途中で
    ! 値が正になることを保証している. また, 元々存在する以上の雨粒が
    ! 蒸気に変換されないように, 元々の雨粒混合比を変換量の上限としている.
    !
    ! 木星の場合は, FactorJ で変換量を加速する. 
    !
        
    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(8), intent(in) :: xz_PotTemp(DimXMin:DimXMax, DimZMin:DimZMax)
                                          !温位の擾乱成分
    real(8), intent(in) :: xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax)
                                          !温度の擾乱成分 
    real(8), intent(in) :: xza_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                          !混合比の擾乱成分
    real(8)             :: DelTime        !時間刻み
    real(8)             :: xza_Rain2Gas(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                          !
    real(8)             :: xz_TempAll(DimXMin:DimXMax, DimZMin:DimZMax)
                                          !温度の擾乱成分 + 平均成分
    real(8)             :: xz_PressAll(DimXMin:DimXMax, DimZMin:DimZMax)
                                          !全圧
    real(8)             :: xza_MixRtAll(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                          !混合比の擾乱成分 + 平均成分
    real(8)             :: xz_NonSaturate(DimXMin:DimXMax, DimZMin:DimZMax)
                                          !未飽和度(飽和混合比と蒸気の混合比の差)
    integer             :: s

    !温度, 圧力, 混合比の全量を求める
    !擾乱成分と平均成分の足し算
    xz_TempAll   = ( xz_PotTemp + xz_PotTempBasicZ ) * ( xz_Exner + xz_ExnerBasicZ )
    xz_PressAll  = PressBasis * ((xz_Exner + xz_ExnerBasicZ) ** (CpDry / GasRDry))
    xza_Rain2Gas = 0.0d0
 
    !混合比は正の値であることを保証
    !移流拡散計算で負になることがあり得るので. 
    xza_MixRtAll = max( 0.0d0, xza_MixRt + xza_MixRtBasicZ )
    
    do s = 1, LoopNum
          
      !飽和蒸気圧と混合比の差(飽和度)を計算. 
      !  雨から蒸気への変換量は飽和度に比例する.
      xz_NonSaturate =                                       &
        & max(                                               &
        &   0.0d0,                                           &
        &   xz_SvapPress(SpcWetID(CloudNum(s)), xz_TempAll)  &
        &     * MolWtWet(CloudNum(s))                        &
        &     / ( MolWtDry * xz_PressAll)                    &
        &     - xza_MixRtAll(:,:,GasNum(s))                  &
        &    )
      
      !雨の変換量
      !  元々の雨粒の混合比以上に蒸発が生じないように上限値を設定
      xza_Rain2Gas(:,:,RainNum(s)) =                                       &
        & - min(                                                           &
        &    DelTime * 4.85d-2 * FactorJ * xz_NonSaturate                  &
        &     * ( xza_MixRtAll(:,:,RainNum(s)) * xz_DensBasicZ )** 0.65d0, &
        &    xza_MixRtAll(:,:,RainNum(s))                                  &
        &   ) 
      
      !蒸気の変換量
      !  雨粒の変換量とは符号が逆となる
      xza_Rain2Gas(:,:,GasNum(s)) = - xza_Rain2Gas(:,:,RainNum(s)) 
    end do
    
  end function xza_Rain2Gas
  

!!!=================================================================================!!!  
  function xza_Rain2GasNH4SH(xz_Exner, xz_PotTemp, xza_MixRt, DelTime)
    !
    ! 雨粒から蒸気への変換量を計算するためのルーチン
    !
    ! 変換量および, 蒸気と雨粒の混合比は正の量なので, 計算の途中途中で
    ! 値が正になることを保証している. また, 元々存在する以上の雨粒が
    ! 蒸気に変換されないように, 元々の雨粒混合比を変換量の上限としている.
    !
        
    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(8), intent(in) :: xz_PotTemp(DimXMin:DimXMax, DimZMin:DimZMax)
                                          !温位の擾乱成分
    real(8), intent(in) :: xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax)
                                          !温度の擾乱成分 
    real(8), intent(in) :: xza_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                          !混合比の擾乱成分
    real(8), intent(in) :: DelTime        !時間刻み
    real(8)             :: xza_Rain2GasNH4SH(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                          !
    real(8)             :: xz_TempAll(DimXMin:DimXMax, DimZMin:DimZMax)
                                          !温度の擾乱成分 + 平均成分
    real(8)             :: xz_PressAll(DimXMin:DimXMax, DimZMin:DimZMax)
                                          !圧力の擾乱成分 + 平均成分
    real(8)             :: xza_MixRtAll(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                          !混合比の擾乱成分 + 平均成分
    real(8)             :: xz_NonSaturate(DimXMin:DimXMax, DimZMin:DimZMax)
                                          !未飽和度(飽和混合比と蒸気の混合比の差)

    !温度, 圧力, 混合比の全量を求める
    !擾乱成分と平均成分の足し算
    xz_TempAll   = ( xz_PotTemp + xz_PotTempBasicZ ) *  ( xz_Exner + xz_ExnerBasicZ )
    xz_PressAll  = PressBasis * ((xz_Exner + xz_ExnerBasicZ) ** (CpDry / GasRDry))    
    xza_Rain2GasNH4SH = 0.0d0

    !混合比は正の値であることを保証
    !移流拡散計算で負になることがあり得るので. 
    xza_MixRtAll = max( 0.0d0, xza_MixRt + xza_MixRtBasicZ )
        
    !飽和蒸気圧と混合比の差(飽和度)を計算. 
    !  雨から蒸気への変換量は飽和度に比例する.
    !  未飽和度を求めたいので, マイナスをかけ算している
    !  (DelMixRtNH4SH は, NH4SH が増加する方向, すなわち飽和度を正としている)
    xz_NonSaturate =                                                 &
      & max(                                                         &
      &  0.0d0,                                                      &
      &   - xz_DelMixRtNH4SH(                                        &
      &       xz_TempAll, xz_PressAll,                               &
      &       xza_MixRtAll(:,:,NH3Num), xza_MixRtAll(:,:,H2SNum),    &
      &       MolWtWet(NH3Num), MolWtWet(H2SNum)                     &
      &     )                                                        &
      &  )
        
    !雨の変換量
    !  元々の雨粒の混合比以上に蒸発が生じないように上限値を設定
    xza_Rain2GasNH4SH(:,:,NH4SHRainNum) =                         &
      & - min(                                                    &
      &     DelTime * 4.85d-2 * FactorJ * xz_NonSaturate          &
      &      * (                                                  &
      &          xza_MixRtAll(:,:,NH4SHRainNum) * xz_DensBasicZ   &
      &         ) ** 0.65d0,                                      &
      &     xza_MixRtAll(:,:,NH4SHRainNum)                        &
      &    ) 
        
    !蒸気の変換量
    !  雨粒の変換量とは符号が逆となる
    xza_Rain2GasNH4SH(:,:,NH3Num) =                              &
      & - xza_Rain2GasNH4SH(:,:,NH4SHRainNum) * MolWtWet(NH3Num) &
      &   / MolWtWet(NH4SHRainNum)
    xza_Rain2GasNH4SH(:,:,H2SNum) =                              &
      & - xza_Rain2GasNH4SH(:,:,NH4SHRainNum) * MolWtWet(H2SNum) &
      &   / MolWtWet(NH4SHRainNum)
    
  end function xza_Rain2GasNH4SH
    

!!!=================================================================================!!!  
  function xz_Rain2GasHeatNH4SH(xz_Exner, xza_DelMixRt)
    !
    ! 雨粒から蒸気への変換量を計算するためのルーチン
    !
    ! 変換量および, 蒸気と雨粒の混合比は正の量なので, 計算の途中途中で
    ! 値が正になることを保証している. また, 元々存在する以上の雨粒が
    ! 蒸気に変換されないように, 元々の雨粒混合比を変換量の上限としている.
    !
        
    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(8), intent(in) :: xza_DelMixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                          !混合比の変化量
    real(8), intent(in) :: xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax)
                                          !温度の擾乱成分 
    real(8)             :: xz_Rain2GasHeatNH4SH(DimXMin:DimXMax, DimZMin:DimZMax)
                                          !
    real(8)             :: xz_ExnerAll(DimXMin:DimXMax, DimZMin:DimZMax)

    !圧力の全量を求める
    !擾乱成分と平均成分の足し算
    xz_ExnerAll    = xz_Exner + xz_ExnerBasicZ

    !雨から蒸気への相変化に伴う発熱
    xz_Rain2GasHeatNH4SH =                                &
      & ReactHeatNH4SH * xza_DelMixRt(:,:,NH4SHRainNum)   &
      &  / (xz_ExnerAll * CpDry)

    !値の保管
    call StoreCond( xz_Rain2GasHeatNH4SH ) 

  end function xz_Rain2GasHeatNH4SH
  

!!!=================================================================================!!!  
  function xz_Rain2GasHeat(xz_PotTemp, xz_Exner, xza_DelMixRt)
    !
    ! 雨粒から蒸気への変換量を計算するためのルーチン
    !
    ! 変換量および, 蒸気と雨粒の混合比は正の量なので, 計算の途中途中で
    ! 値が正になることを保証している. また, 元々存在する以上の雨粒が
    ! 蒸気に変換されないように, 元々の雨粒混合比を変換量の上限としている.
    !
        
    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(8), intent(in) :: xz_PotTemp(DimXMin:DimXMax, DimZMin:DimZMax)
                                          !温位の擾乱成分
    real(8), intent(in) :: xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax)
                                          !温度の擾乱成分 
    real(8), intent(in) :: xza_DelMixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                          !混合比の変化
    real(8)             :: xz_Rain2GasHeat(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)             :: xza_LatentHeat(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
    real(8)             :: xz_TempAll(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)             :: xz_ExnerAll(DimXMin:DimXMax, DimZMin:DimZMax)
    integer             :: s

    !温度, 圧力, 混合比の全量を求める
    !擾乱成分と平均成分の足し算
    xz_ExnerAll    = xz_Exner + xz_ExnerBasicZ
    xz_TempAll     = ( xz_PotTemp + xz_PotTempBasicZ ) *  ( xz_Exner + xz_ExnerBasicZ )
    xza_LatentHeat = 0.0d0

    !雨から蒸気への相変化に伴う発熱    
    do s = 1, LoopNum
      xza_LatentHeat(:,:,s) =                                &
        & xz_LatentHeat( SpcWetID(RainNum(s)), xz_TempAll )  &
        &  * xza_DelMixRt(:,:,RainNum(s))                    &
        &  / (xz_ExnerAll * CpDry) 
    end do
    xz_Rain2GasHeat = sum( xza_LatentHeat, 3 )

    !値の保管
    call StoreCond( xz_Rain2GasHeat ) 
    
  end function xz_Rain2GasHeat
  

!!!=================================================================================!!!  
  function xza_Cloud2Rain( xza_MixRt, DelTime )
    !
    ! 雲粒から雨粒への変換量を計算するためのルーチン
    ! 併合成長は Berry (1968) のパラメタリゼーションを利用し, 
    ! 衝突合体成長は Kessler (1969) のパラメタリゼーションを利用する. 
    !
    ! 変換量および, 雲粒と雨粒の混合比は正の量なので, 計算の途中途中で
    ! 値が正になることを保証している. また, 元々存在する以上の雲粒が
    ! 雨粒に変換されないように, 元々の雲粒混合比を変換量の上限としている.
    ! 正の値を保証するために, 引数として時間刻みが必要となる. 
    ! (AutoConv, Collect は時間刻み幅での積分値を計算)
    !
    ! このルーチンでは, 凝縮物質と反応生成物とを区別する必要が全くないので, 
    ! ループを回す回数を LoopNum2 回としている. 
    !
    
    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(8), intent(in) :: xza_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                          !混合比の擾乱成分
    real(8)             :: DelTime        !時間刻み
    real(8)             :: xza_Cloud2Rain(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                          !雲から雨への変換量
    real(8)             :: xza_MixRtAll(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                          !混合比の擾乱成分 + 平均成分
    real(8)             :: xz_AutoConv(DimXMin:DimXMax, DimZMin:DimZMax)
                                          !飽和混合比
    real(8)             :: xz_Collect(DimXMin:DimXMax, DimZMin:DimZMax)
                                          !規格化された潜熱
    real(8), parameter  :: N0 = 5.0d7 
    real(8), parameter  :: D0 = 3.66d-1
    integer             :: s


    xza_Cloud2Rain  = 0.0d0

    !混合比は正の値を保証
    !移流拡散で負になることもあり得るので. 
    xza_MixRtAll = max( 0.0d0, xza_MixRt + xza_MixRtBasicZ )


    do s = 1, LoopNum2
      xz_AutoConv = 0.0d0
      xz_Collect  = 0.0d0
      
      !併合成長
!     !  Kessler (1969) のパラメタリゼーション        

!      xz_AutoConv =                                                 &
!        & DelTime / AutoConvTime                                    &
!        & * max( 0.0d0, ( xza_MixRtAll(:,:,CloudNum(s)) - MixRt_AutoConvCr) )

      !  Berry (1968) のパラメタリゼーション      
      xz_AutoConv =                                                   &
        & DelTime                                                     &
        & * xz_DensBasicZ                                             &
        & * ( xza_MixRtAll(:,:,CloudNum(s)) ** 3.0d0  ) * 1.0d6       &
        & / ( 60.0d0                                                  &
        &     * (                                                     &
        &         2.0d0 * xza_MixRtAll(:,:,CloudNum(s))               &
        &       + 2.66d-8 * N0 / ( xz_DensBasicZ * D0 )               &
        &      )                                                      &
        &   )

      !衝突合体成長
      !  Kessler (1969) のパラメタリゼーション    
      xz_Collect =                                               &
        &  DelTime                                               &
        &  * 2.2d0 * FactorJ * xza_MixRtAll(:,:,CloudNum(s))     &
        &  * (                                                   &
        &       xza_MixRtAll(:,:,RainNum(s)) * xz_DensBasicZ     &
        &     ) ** 0.875d0  
      
      !雲の変換量: 併合成長と合体衝突の和
      !  元々の変化量を上限値として設定する. 負の値となる.
      xza_Cloud2Rain(:,:,CloudNum(s)) =                         &
        & - min(                                                &
        &         xza_MixRtAll(:,:,CloudNum(s)),                &
        &         ( xz_AutoConv + xz_Collect )                  &
        &       )
      
      !雨の変換量. 符号は雲の変換量とは反対. 
      xza_Cloud2Rain(:,:,RainNum(s)) = - xza_Cloud2Rain(:,:,CloudNum(s)) 
          
    end do

!    write(*,*) 'C2R: ', minval(xza_Cloud2Rain(:,:,1)), maxval(xza_Cloud2Rain(:,:,1))
!    write(*,*) 'C2R: ', minval(xza_Cloud2Rain(:,:,2)), maxval(xza_Cloud2Rain(:,:,2))
!    write(*,*) 'C2R: ', minval(xza_Cloud2Rain(:,:,3)), maxval(xza_Cloud2Rain(:,:,3))
    
  end function xza_Cloud2Rain
  

!!!=================================================================================!!!
  function xz_FallRain( xz_MixRt, num )
    !
    ! 雨粒の落下による移流を求める. 
    ! 
    ! このルーチンの引数 xz_MixRt は 2 次元配列である. 
    ! 引数に与えられた混合比に対し, 移流を計算する. 
    ! もしも, xza_MixRt が雨粒の混合比を表すなら(RainSW = 1)なら
    ! その値を出力する. 
    ! もしも, xza_MixRt が蒸気もしくは雲粒の混合比を表すなら(RainSW = 0)
    ! ならば, 落下による移流はゼロとする. 
    !

    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    integer, intent(in) :: num
    real(8), intent(in) :: xz_MixRt(DimXMin:DimXMax, DimZMin:DimZMax)  
                                                 !蒸気混合比(擾乱)
    real(8)             :: xz_MixRtAll(DimXMin:DimXMax, DimZMin:DimZMax)  
                                                 !蒸気混合比(擾乱 + 平均場)
    real(8)             :: xz_FallRain(DimXMin:DimXMax, DimZMin:DimZMax)
                                                 !雨粒の落下効果
    real(8)             :: xz_VelZRain(DimXMin:DimXMax, DimZMin:DimZMax)
                                                 !雨粒落下速度

    xz_MixRtAll = max( 0.0d0, xz_MixRt + xza_MixRtBasicZ(:,:,num) )
    xz_FallRain = 0.0d0
    xz_VelZRain = 0.0d0

    
    !雨粒終端速度
    xz_VelZRain = 12.2d0 * FactorJ * ( xz_MixRtAll ** 0.125d0 )
      
    !落下による移流
    !  Dens の avr を取ってから割ると, ゼロ割が生じるので注意
    xz_FallRain =                                                     &
      &  xz_avr_xr(                                                   &
      &         xr_dz_xz(xz_DensBasicZ * xz_VelZRain * xz_MixRtAll)  &
      &      ) / xz_DensBasicZ                                        &
      &        * RainSW(num)

!    write(*,*) 'MixRt: ', minval( xz_MixRt    ),  maxval( xz_MixRt    )
!    write(*,*) 'Fall:  ', minval( xz_FallRain ),  maxval( xz_FallRain )

  end function xz_FallRain
  
end module WarmRainPrm
