!= Module WarmRainPrm_3d
!
! Authors::   SUGIYAMA Ko-ichiro, ODAKA Masatsugu
! Version::   $Id: warmrainprm_3d.f90,v 1.2 2011-03-28 01:48:54 sugiyama Exp $
! Tag Name::  $Name: arare4-20111215 $
! 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_3d
  !
  !暖かい雨のバルク法を用いた, 水蒸気と雨, 雲と雨の混合比の変換係数を求める.
  !   * 中島健介 (1994) で利用した定式をそのまま利用. 
  ! 
  
  !モジュール読み込み
  use dc_types,   only : DP
  use dc_iounit,  only : FileOpen
  use dc_message, only : MessageNotify

  use gridset_3d, only : DimXMin,           &!x 方向の配列の下限
    &                    DimXMax,           &!x 方向の配列の上限
    &                    DimYMin,           &!y 方向の配列の上限
    &                    DimYMax,           &!y 方向の配列の上限
    &                    DimZMin,           &!z 方向の配列の下限
    &                    DimZMax,           &!z 方向の配列の上限
    &                    SpcNum              !化学種の数
  use basicset_3d,only : PressBasis,        &!温位の基準圧力 
    &                    CpDry,             &!乾燥成分の比熱
    &                    MolWtWet,          &!
    &                    MolWtDry,          &!
    &                    SpcWetID,          &!
    &                    xyz_DensBasicZ,    &!基本場の密度
    &                    xyz_PotTempBasicZ, &!基本場の温位
    &                    xyz_ExnerBasicZ,   &!基本場の無次元圧力
    &                    xyza_MixRtBasicZ,  &!基本場の混合比
    &                    GasRDry             !乾燥成分の気体定数 
  use moistset_3d, only : CondNum,          &!凝結過程の数
    &                  IdxCG,            &!凝結過程(蒸気)の配列添え字
    &                  IdxCC,            &!凝結過程(雲)の配列添え字
    &                  IdxCR,            &!凝結過程(雲)の配列添え字
    &                  CloudNum,         &!雲の数
    &                  RainNum,          &!雨の数
    &                  IdxC,             &!雲の配列添え字
    &                  IdxR,             &!雨の配列添え字
    &                  IdxNH3,           &!NH3(蒸気)の配列添え字
    &                  IdxH2S,           &!H2S(蒸気)の配列添え字
    &                  IdxNH4SHr          !NH4SH(雨)の配列添え字

  use xyz_base_module, only : xyz_avr_xyr
  use xyz_deriv_module,only : xyr_dz_xyz

  use ChemCalc_3d,  only : xyz_SvapPress, xyz_LatentHeat, ReactHeatNH4SH
  use MoistFunc_3d, only : xyz_DelMixRtNH4SH
  use StoreMixRt_3d,only : StoreMixRtRain
  
  !暗黙の型宣言禁止
  implicit none

  !属性の指定
  private

  !関数を public にする
  public WarmRainPrm_Init
  public xyz_Rain2GasHeat
  public xyz_Rain2GasHeatNH4SH
  public xyza_Rain2Gas
  public xyza_Rain2GasNH4SH
  public xyza_Cloud2Rain
  public xyza_FallRain

  real(DP)      :: FactorJ      = 1.0d0 !雲物理過程のパラメータ
                                       !木星では 3.0d0
                                       !地球では 1.0d0 とする
  real(DP)      :: AutoConvTime = 1.0d3 !併合成長の時定数 [sec]
  real(DP)      :: MixRt_AutoConvCr = 1.0d-3 
                                       !併合成長を生じる臨界混合比 [kg/kg]

  save FactorJ, AutoConvTime, MixRt_AutoConvCr

contains  

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

    !暗黙の型宣言禁止
    implicit none

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

    !内部変数
    integer                  :: unit    !装置番号

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

    call FileOpen(unit, file=cfgfile, mode='r')
    read(unit, NML=warmrainprm)
    close(unit)

    call MessageNotify( "M", &
      &  "WarmRainPrm_Init", "FactorJ = %f",  d=(/FactorJ/) )
    call MessageNotify( "M", &
      &  "WarmRainPrm_Init", "AutoConvTime = %f",  d=(/AutoConvTime/) )
    call MessageNotify( "M", &
      &  "WarmRainPrm_Init", "MixRt_AutoConvCr = %f",  d=(/MixRt_AutoConvCr/) )
    
  end subroutine WarmRainPrm_Init


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

    !変数定義
    real(DP), intent(in) :: xyz_PotTemp(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                          !温位の擾乱成分
    real(DP), intent(in) :: xyz_Exner(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                          !温度の擾乱成分 
    real(DP), intent(in) :: xyza_MixRt(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum)
                                          !混合比の擾乱成分
    real(DP)             :: DelTime        !時間刻み
    real(DP)             :: xyza_Rain2Gas(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum)
                                          !
    real(DP)             :: xyz_TempAll(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                          !温度の擾乱成分 + 平均成分
    real(DP)             :: xyz_PressAll(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                          !全圧
    real(DP)             :: xyza_MixRtAll(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum)
                                          !混合比の擾乱成分 + 平均成分
    real(DP)             :: xyz_NonSaturate(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                          !未飽和度(飽和混合比と蒸気の混合比の差)
    integer             :: s

    !温度, 圧力, 混合比の全量を求める
    !擾乱成分と平均成分の足し算
    xyz_TempAll   = ( xyz_PotTemp + xyz_PotTempBasicZ ) * ( xyz_Exner + xyz_ExnerBasicZ )
    xyz_PressAll  = PressBasis * ((xyz_Exner + xyz_ExnerBasicZ) ** (CpDry / GasRDry))
    xyza_Rain2Gas = 0.0d0
 
    !混合比は正の値であることを保証
    !移流拡散計算で負になることがあり得るので. 
    xyza_MixRtAll = max( 0.0d0, xyza_MixRt + xyza_MixRtBasicZ )
    
    do s = 1, CondNum
          
      !飽和蒸気圧と混合比の差(飽和度)を計算. 
      !  雨から蒸気への変換量は飽和度に比例する.
      xyz_NonSaturate =                                      &
        & max(                                               &
        &   0.0d0,                                           &
        &   xyz_SvapPress(SpcWetID(IdxCC(s)), xyz_TempAll)   &
        &     * MolWtWet(IdxCG(s))                           &
        &     / ( MolWtDry * xyz_PressAll)                   &
        &     - xyza_MixRtAll(:,:,:,IdxCG(s))                &
        &    )
      
      !雨の変換量
      !  元々の雨粒の混合比以上に蒸発が生じないように上限値を設定
      xyza_Rain2Gas(:,:,:,IdxCR(s)) =                                       &
        & - min(                                                            &
        &    DelTime * 4.85d-2 * FactorJ * xyz_NonSaturate                  &
        &     * ( xyza_MixRtAll(:,:,:,IdxCR(s)) * xyz_DensBasicZ )** 0.65d0,&
        &    xyza_MixRtAll(:,:,:,IdxCR(s))                                  &
        &   ) 
      
      !蒸気の変換量
      !  雨粒の変換量とは符号が逆となる
      xyza_Rain2Gas(:,:,:,IdxCG(s)) = - xyza_Rain2Gas(:,:,:,IdxCR(s)) 
    end do
    
  end function xyza_Rain2Gas
  

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

    !変数定義
    real(DP), intent(in) :: xyz_PotTemp(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                          !温位の擾乱成分
    real(DP), intent(in) :: xyz_Exner(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                          !温度の擾乱成分 
    real(DP), intent(in) :: xyza_MixRt(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum)
                                          !混合比の擾乱成分
    real(DP), intent(in) :: DelTime        !時間刻み
    real(DP)             :: xyza_Rain2GasNH4SH(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum)
                                          !
    real(DP)             :: xyz_TempAll(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                          !温度の擾乱成分 + 平均成分
    real(DP)             :: xyz_PressAll(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                          !圧力の擾乱成分 + 平均成分
    real(DP)             :: xyza_MixRtAll(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum)
                                          !混合比の擾乱成分 + 平均成分
    real(DP)             :: xyz_NonSaturate(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                          !未飽和度(飽和混合比と蒸気の混合比の差)

    !温度, 圧力, 混合比の全量を求める
    !擾乱成分と平均成分の足し算
    xyz_TempAll = ( xyz_PotTemp + xyz_PotTempBasicZ ) *  ( xyz_Exner + xyz_ExnerBasicZ )
    xyz_PressAll= PressBasis * ((xyz_Exner + xyz_ExnerBasicZ) ** (CpDry / GasRDry))    
    xyza_Rain2GasNH4SH = 0.0d0

    !混合比は正の値であることを保証
    !移流拡散計算で負になることがあり得るので. 
    xyza_MixRtAll = max( 0.0d0, xyza_MixRt + xyza_MixRtBasicZ )
        
    !飽和蒸気圧と混合比の差(飽和度)を計算. 
    !  雨から蒸気への変換量は飽和度に比例する.
    !  未飽和度を求めたいので, マイナスをかけ算している
    !  (DelMixRtNH4SH は, NH4SH が増加する方向, すなわち飽和度を正としている)
    xyz_NonSaturate =                      &
      & max(                               &
      &  0.0d0,                            &
      &   - xyz_DelMixRtNH4SH(             &  
      &       xyz_TempAll, xyz_PressAll,   &
      &       xyza_MixRtAll(:,:,:,IdxNH3), &
      &       xyza_MixRtAll(:,:,:,IdxH2S), &
      &       MolWtWet(IdxNH3),            &
      &       MolWtWet(IdxH2S)             &
      &     )                                                         &
      &  )
        
    !雨の変換量
    !  元々の雨粒の混合比以上に蒸発が生じないように上限値を設定
    xyza_Rain2GasNH4SH(:,:,:,IdxNH4SHr) =                        &
      & - min(                                                   &
      &     DelTime * 4.85d-2 * FactorJ * xyz_NonSaturate        &
      &      * (                                                 &
      &          xyza_MixRtAll(:,:,:,IdxNH4SHr) * xyz_DensBasicZ &
      &         ) ** 0.65d0,                                     &
      &     xyza_MixRtAll(:,:,:,IdxNH4SHr)                       &
      &    ) 
        
    !蒸気の変換量
    !  雨粒の変換量とは符号が逆となる
    xyza_Rain2GasNH4SH(:,:,:,IdxNH3) =                           &
      & - xyza_Rain2GasNH4SH(:,:,:,IdxNH4SHr) * MolWtWet(IdxNH3) &
      &   / MolWtWet(IdxNH4SHr)

    xyza_Rain2GasNH4SH(:,:,:,IdxH2S) =                           &
      & - xyza_Rain2GasNH4SH(:,:,:,IdxNH4SHr) * MolWtWet(IdxH2S) &
      &   / MolWtWet(IdxNH4SHr)
    
  end function xyza_Rain2GasNH4SH
    

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

    !変数定義
    real(DP), intent(in) :: xyza_DelMixRt(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum)
                                          !混合比の変化量
    real(DP), intent(in) :: xyz_Exner(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                          !温度の擾乱成分 
    real(DP)             :: xyz_Rain2GasHeatNH4SH(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                          !
    real(DP)             :: xyz_ExnerAll(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)

    !圧力の全量を求める
    !擾乱成分と平均成分の足し算
    xyz_ExnerAll    = xyz_Exner + xyz_ExnerBasicZ

    !雨から蒸気への相変化に伴う発熱
    xyz_Rain2GasHeatNH4SH =                                             &
      & ReactHeatNH4SH * xyza_DelMixRt(:,:,:,IdxNH4SHr) / (xyz_ExnerAll * CpDry)

  end function xyz_Rain2GasHeatNH4SH
  

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

    !変数定義
    real(DP), intent(in) :: xyz_PotTemp(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                          !温位の擾乱成分
    real(DP), intent(in) :: xyz_Exner(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                          !温度の擾乱成分 
    real(DP), intent(in) :: xyza_DelMixRt(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum)
                                          !混合比の変化
    real(DP)             :: xyz_Rain2GasHeat(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
    real(DP)             :: xyza_LatentHeat(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum)
    real(DP)             :: xyz_TempAll(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
    real(DP)             :: xyz_ExnerAll(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
    integer             :: s

    !温度, 圧力, 混合比の全量を求める
    !擾乱成分と平均成分の足し算
    xyz_ExnerAll = xyz_Exner + xyz_ExnerBasicZ
    xyz_TempAll  = ( xyz_PotTemp + xyz_PotTempBasicZ ) * ( xyz_Exner + xyz_ExnerBasicZ )
    xyza_LatentHeat = 0.0d0
    
    !雨から蒸気への相変化に伴う発熱    
    do s = 1, CondNum
      xyza_LatentHeat(:,:,:,s) =                              &
        & xyz_LatentHeat( SpcWetID(IdxCR(s)), xyz_TempAll )  &
        &  * xyza_DelMixRt(:,:,:,IdxCR(s))                    &
        &  / (xyz_ExnerAll * CpDry) 
    end do
    xyz_Rain2GasHeat = sum( xyza_LatentHeat, 4 )

  end function xyz_Rain2GasHeat
  

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

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


    xyza_Cloud2Rain  = 0.0d0

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


    do s = 1, CloudNum
      xyz_AutoConv = 0.0d0
      xyz_Collect  = 0.0d0
      
      !併合成長
      !  Kessler (1969) のパラメタリゼーション        
      xyz_AutoConv =                                                 &
        & DelTime / AutoConvTime                                    &
        & * max( 0.0d0, ( xyza_MixRtAll(:,:,:,IdxC(s)) - MixRt_AutoConvCr) )

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

      !衝突合体成長
      !  Kessler (1969) のパラメタリゼーション    
      xyz_Collect =                                           &
        &  DelTime                                            &
        &  * 2.2d0 * FactorJ * xyza_MixRtAll(:,:,:,IdxC(s))   &
        &  * (                                                &
        &       xyza_MixRtAll(:,:,:,IdxR(s)) * xyz_DensBasicZ &
        &     ) ** 0.875d0  
      
      !雲の変換量: 併合成長と合体衝突の和
      !  元々の変化量を上限値として設定する. 負の値となる.
      xyza_Cloud2Rain(:,:,:,IdxC(s)) =                       &
        & - min(                                             &
        &         xyza_MixRtAll(:,:,:,IdxC(s)),              &
        &         ( xyz_AutoConv + xyz_Collect )             &
        &       )
      
      !雨の変換量. 符号は雲の変換量とは反対. 
      xyza_Cloud2Rain(:,:,:,IdxR(s)) = - xyza_Cloud2Rain(:,:,:,IdxC(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 xyza_Cloud2Rain
  

!!!=================================================================================!!!
  function xyza_FallRain( xyza_MixRt )
    !
    ! 雨粒の落下による移流を求める. 
    ! 

    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP), intent(in) :: xyza_MixRt(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum) 
                                                 !蒸気混合比(擾乱)
    real(DP)  :: xyza_MixRtAll(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum)
                                                 !蒸気混合比(擾乱 + 平均場)
    real(DP)  :: xyza_FallRain(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum)
                                                 !雨粒の落下効果
    real(DP)  :: xyz_VelZRain(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                                 !雨粒落下速度
    integer  :: s

    xyza_MixRtAll = max( 0.0d0, xyza_MixRt + xyza_MixRtBasicZ )
    xyza_FallRain = 0.0d0
    xyz_VelZRain = 0.0d0
    
    !落下による移流
    !  Dens の avr を取ってから割ると, ゼロ割が生じるので注意
    do s = 1, RainNum
      !雨粒終端速度
      xyz_VelZRain = 12.2d0 * FactorJ * ( xyza_MixRtAll(:,:,:,IdxR(s)) ** 0.125d0 )
      
      xyza_FallRain(:,:,:,IdxR(s)) =                   &
        &  xyz_avr_xyr(                                &
        &     xyr_dz_xyz(xyz_DensBasicZ                &
        &               * xyz_VelZRain                 &
        &               * xyza_MixRtAll(:,:,:,IdxR(s)) &
        &              )                               &
        &      ) / xyz_DensBasicZ                      
    end do
    
    call StoreMixRtRain( xyza_FallRain )
    
  end function xyza_FallRain
  
  
end module WarmRainPrm_3d
