subroutine CloudPhys_K1969_FallRain(xyzf_QMix, xyzf_DQMixDt )
    !
    ! 雨粒の落下による移流を求める. 
    ! 

    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP), intent(in) :: xyzf_QMix(imin:imax,jmin:jmax,kmin:kmax, ncmax) 
                                                 !蒸気混合比(擾乱)
    real(DP), intent(inout) :: xyzf_DQMixDt(imin:imax,jmin:jmax,kmin:kmax, ncmax) 
                                                 !蒸気混合比の変化量
    real(DP)  :: xyzf_QMixAll(imin:imax,jmin:jmax,kmin:kmax, ncmax)
                                                 !蒸気混合比(擾乱 + 平均場)
    real(DP)  :: xyzf_FallRain(imin:imax,jmin:jmax,kmin:kmax, ncmax)
                                                 !雨粒の落下効果
    real(DP)  :: xyz_VelZRain(imin:imax,jmin:jmax,kmin:kmax)
                                                 !雨粒落下速度
    real(DP)  :: xyzf_DQMixDtOrig(imin:imax,jmin:jmax,kmin:kmax, ncmax) 
                                                 !蒸気混合比の変化量
    integer  :: s, l

    xyzf_QMixAll = max( 0.0d0, xyzf_QMix + xyzf_QMixBZ )
    xyzf_DQMixDtOrig = xyzf_DQMixDt
    xyzf_FallRain = 0.0d0
    xyz_VelZRain = 0.0d0
    
    !落下による移流
    !  Dens の avr を取ってから割ると, ゼロ割が生じるので注意
    do s = 1, RainNum
      !雨粒終端速度
      xyz_VelZRain = 12.2d0 * FactorJ * ( xyzf_QMixAll(:,:,:,IdxR(s)) ** 0.125d0 )
      
      xyzf_FallRain(:,:,:,IdxR(s)) = xyz_avr_xyr( xyr_dz_xyz(xyz_DensBZ * xyz_VelZRain * xyzf_QMixAll(:,:,:,IdxR(s)) ) ) / xyz_DensBZ                      
    end do

    xyzf_DQMixDt = xyzf_DQMixDtOrig + xyzf_FallRain
    
    do l = 1, ncmax
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(l))//'_Fall', xyzf_FallRain(1:nx, 1:ny, 1:nz, l))
    end do
    
  end subroutine CloudPhys_K1969_FallRain

  subroutine Cloudphys_K1969_Init

    !暗黙の型宣言禁止
    implicit none

    !内部変数
    integer  :: unit    !装置番号
    integer  :: l
    character(STRING) :: Planet = ""

    !-----------------------------------------------------------
    ! NAMELIST から情報を取得
    !-----------------------------------------------------------
    ! NAMELIST から情報を取得
    NAMELIST /cloudphys_k1969_nml/ Planet, FactorJ, AutoConvTime, QMixCr

    call FileOpen(unit, file=namelist_filename, mode='r')
    read(unit, NML=cloudphys_k1969_nml)
    close(unit)

    if (trim(Planet) == "Earth") then 
      FactorJ = 1.0d0
    elseif (trim(Planet) == "Jupiter") then 
      FactorJ = 3.0d0
    end if

    if (myrank == 0) then 
      call MessageNotify( "M", "WarmRainPrm_Init", "Planet = %c",  c1=trim(Planet))
      call MessageNotify( "M", "WarmRainPrm_Init", "FactorJ = %f",  d=(/FactorJ/) )
      call MessageNotify( "M", "WarmRainPrm_Init", "AutoConvTime = %f",  d=(/AutoConvTime/) )
      call MessageNotify( "M", "WarmRainPrm_Init", "QMixCr = %f",  d=(/QMixCr/) )
    end if

    call HistoryAutoAddVariable( varname='PTempCond', dims=(/'x','y','z','t'/), longname='Latent heat term of potential temperature', units='K.s-1', xtype='double')

    do l = 1, ncmax
      call HistoryAutoAddVariable( varname=trim(SpcWetSymbol(l))//'_Cond', dims=(/'x','y','z','t'/), longname='Condensation term of ' //trim(SpcWetSymbol(l))//' mixing ratio', units='kg.kg-1.s-1', xtype='double')
      
      call HistoryAutoAddVariable( varname=trim(SpcWetSymbol(l))//'_Fall', dims=(/'x','y','z','t'/), longname='Fall Rain term of ' //trim(SpcWetSymbol(l))//' mixing ratio', units='kg.kg-1.s-1', xtype='double')
    end do
    
  end subroutine Cloudphys_K1969_Init

  subroutine Cloudphys_K1969_forcing(xyz_ExnerNl, xyz_PTempAl, xyzf_QMixAl)

    implicit none

    real(DP), intent(in)           :: xyz_ExnerNl(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), intent(inout)        :: xyz_PTempAl(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), intent(inout)        :: xyzf_QMixAl(imin:imax, jmin:jmax, kmin:kmax, ncmax)
    real(DP)                       :: xyz_PTempOrig(imin:imax, jmin:jmax, kmin:kmax)
    real(DP)                       :: xyz_PTempWork(imin:imax, jmin:jmax, kmin:kmax)
    real(DP)                       :: xyz_DelPTemp(imin:imax, jmin:jmax, kmin:kmax)
    real(DP)                       :: xyz_Del(imin:imax, jmin:jmax, kmin:kmax)
    real(DP)                       :: xyzf_QMixOrig(imin:imax, jmin:jmax, kmin:kmax, ncmax)
    real(DP)                       :: xyzf_QMixWork(imin:imax, jmin:jmax, kmin:kmax, ncmax)
    real(DP)                       :: xyzf_DelQMix(imin:imax, jmin:jmax, kmin:kmax, ncmax)
    real(DP)                       :: xyzf_Del(imin:imax, jmin:jmax, kmin:kmax, ncmax)
    real(DP)                       :: DelTime
    integer                        :: l, s

    real(DP)             :: xyzf_Cloud2Rain(imin:imax,jmin:jmax,kmin:kmax, ncmax)
                                          !雲から雨への変換量
    real(DP)             :: xyz_AutoConv(imin:imax,jmin:jmax,kmin:kmax)
                                          !飽和混合比
    real(DP)             :: xyz_Collect(imin:imax,jmin:jmax,kmin:kmax)
                                          !規格化された潜熱

    real(DP)             :: xyzf_QMixAll(imin:imax,jmin:jmax,kmin:kmax, ncmax)
                                          !混合比の擾乱成分 + 平均成分
    real(DP)             :: xyz_TempAll(imin:imax,jmin:jmax,kmin:kmax)
                                          !温度の擾乱成分 + 平均成分
    real(DP)             :: xyz_PressAll(imin:imax,jmin:jmax,kmin:kmax)
                                          !全圧
    real(DP)             :: xyz_ExnerAll(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyz_NonSaturate(imin:imax,jmin:jmax,kmin:kmax)
                                          !未飽和度(飽和混合比と蒸気の混合比の差)
    real(DP)             :: xyzf_Rain2Gas(imin:imax,jmin:jmax,kmin:kmax, ncmax)
    real(DP)             :: xyzf_Rain2GasNH4SH(imin:imax,jmin:jmax,kmin:kmax, ncmax)
    real(DP)             :: xyzf_DelPTemp(imin:imax,jmin:jmax,kmin:kmax, ncmax)
    real(DP)             :: xyz_DelPTempNH4SH(imin:imax,jmin:jmax,kmin:kmax)

    !-----------------------------------------
    ! 時間刻み幅. Leap-frog なので, 2 \del t
    !
    DelTime = 2.0d0 * DelTimeLong

    !------------------------------------------
    ! 初期値を保管 Store Initial Value
    !
    xyz_PTempOrig = xyz_PTempAl
    xyzf_QMixOrig  = xyzf_QMixAl    

    !------------------------------------------    
    ! 暖かい雨のパラメタリゼーション.
    ! * 雲<-->雨 の変換を行う.
    !
    ! Warm rain parameterization.
    ! * Conversion from cloud to rain.
    
    !これまでの値を作業配列に保管
    ! Previous values are stored to work area.
    !
    xyzf_QMixWork = xyzf_QMixAl
    
    !雨への変化量を計算
    ! Conversion values are calculated.
    !    
    xyzf_QMixAll = max( 1.0d-60, xyzf_QMixAl + xyzf_QMixBZ )

    xyzf_Cloud2Rain = 0.0d0

    do s = 1, CloudNum
      xyz_AutoConv = 0.0d0
      xyz_Collect  = 0.0d0
      
      !併合成長
      !
      xyz_AutoConv = DelTime / AutoConvTime * max( 1.0d-60, ( xyzf_QMixAll(:,:,:,IdxC(s)) - QMixCr) )

      !衝突合体成長
      !
      xyz_Collect = DelTime * 2.2d0 * FactorJ * xyzf_QMixAll(:,:,:,IdxC(s)) * (xyzf_QMixAll(:,:,:,IdxR(s)) * xyz_DensBZ) ** 0.875d0  

      !雲の変換量: 併合成長と合体衝突の和
      !  元々の変化量を上限値として設定する. 負の値となる.
      !
      xyzf_Cloud2Rain(:,:,:,IdxC(s)) = - min( xyzf_QMixAll(:,:,:,IdxC(s)), ( xyz_AutoConv + xyz_Collect ) )
      
      !雨の変換量. 符号は雲の変換量とは反対. 
      xyzf_Cloud2Rain(:,:,:,IdxR(s)) = - xyzf_Cloud2Rain(:,:,:,IdxC(s)) 
    end do

    ! 変化量を足し込む
    !
    xyzf_QMixAl = xyzf_QMixWork + xyzf_Cloud2Rain

    ! Set Margin
    !
    call SetMargin_xyzf(xyzf_QMixAl)

    !-------------------------------------------
    ! 湿潤飽和調節
    ! * 蒸気<-->雲の変換を行う.
    !
    ! Moist adjustment.
    ! * Conversion from vapor to cloud.
    !
    call MoistAdjustSvapPress( xyz_ExnerNl, xyz_PTempAl, xyzf_QMixAl )
    if (IdxNH4SHr /= 0) then 
      call MoistAdjustNH4SH( xyz_ExnerNl, xyz_PTempAl, xyzf_QMixAl )
    end if

    ! Set Margin
    !
    call SetMargin_xyz(xyz_PTempAl)
    call SetMargin_xyzf(xyzf_QMixAl)

    !-------------------------------------------    
    ! 暖かい雨のパラメタリゼーション.
    ! * 蒸気<-->雨 の変換を行う
    !
    ! Warm rain parameterization.
    ! * Conversion from rain to vapor.
    
    !これまでの値を作業配列に保管
    ! Previous values are stored to work area.
    !
    xyz_PTempWork = xyz_PTempAl
    xyzf_QMixWork  = xyzf_QMixAl
    
    ! 雨から蒸気への混合比変化を求める
    ! * 温位の計算において, 混合比変化が必要となるため, 
    !   混合比変化を 1 つの配列として用意する.
    !
    ! Conversion values are calculated.
    !

    !温度, 圧力, 混合比の全量を求める
    !擾乱成分と平均成分の足し算
    !
    xyz_ExnerAll  = xyz_ExnerNl + xyz_ExnerBZ
    xyz_TempAll   = ( xyz_PTempAl + xyz_PTempBZ ) * ( xyz_ExnerNl + xyz_ExnerBZ )
    xyz_PressAll  = PressBasis * ((xyz_ExnerNl + xyz_ExnerBZ) ** (CpDry / GasRDry))
    xyzf_QMixAll = max( 1.0d-60, xyzf_QMixAl + xyzf_QMixBZ )
    xyzf_Rain2Gas = 0.0d0
    xyzf_Rain2GasNH4SH = 0.0d0
    xyz_NonSaturate = 0.0d0
    xyzf_DelPTemp = 0.0d0

    do s = 1, CondNum
      !飽和蒸気圧と混合比の差(飽和度)を計算. 
      !  雨から蒸気への変換量は飽和度に比例する.
      !
      xyz_NonSaturate = max( 1.0d-60, xyz_SvapPress(SpcWetID(IdxCC(s)), xyz_TempAll) * MolWtWet(IdxCG(s)) / ( MolWtDry * xyz_PressAll) - xyzf_QMixAll(:,:,:,IdxCG(s)) )

      !雨の変換量
      !  元々の雨粒の混合比以上に蒸発が生じないように上限値を設定
      !
      xyzf_Rain2Gas(:,:,:,IdxCR(s)) = - min( DelTime * 4.85d-2 * FactorJ * xyz_NonSaturate * ( xyzf_QMixAll(:,:,:,IdxCR(s)) * xyz_DensBZ )** 0.65d0, xyzf_QMixAll(:,:,:,IdxCR(s)) ) 

      !蒸気の変換量
      !  雨粒の変換量とは符号が逆となる
      !
      xyzf_Rain2Gas(:,:,:,IdxCG(s)) = - xyzf_Rain2Gas(:,:,:,IdxCR(s)) 
    
      ! xyzf_DelQMix を元に潜熱を計算
      !
      xyzf_DelPTemp(:,:,:,s) = xyz_LatentHeat( SpcWetID(IdxCR(s)), xyz_TempAll ) * xyzf_Rain2Gas(:,:,:,IdxCR(s)) / (xyz_ExnerAll * CpDry) 

    end do

    !飽和蒸気圧と混合比の差(飽和度)を計算. 
    !  雨から蒸気への変換量は飽和度に比例する.
    !  未飽和度を求めたいので, マイナスをかけ算している
    !  (DelQMixNH4SH は, NH4SH が増加する方向, すなわち飽和度を正としている)
    !
    xyz_NonSaturate    = 0.0d0
    xyzf_Rain2GasNH4SH = 0.0d0
    xyz_DelPTempNH4SH  = 0.0d0

    if (IdxNH4SHr /= 0) then 
      xyz_NonSaturate = max( 1.0d-60, - xyz_DelQMixNH4SH( xyz_TempAll, xyz_PressAll, xyzf_QMixAll(:,:,:,IdxNH3), xyzf_QMixAll(:,:,:,IdxH2S), MolWtWet(IdxNH3), MolWtWet(IdxH2S) ) )

      !雨の変換量
      !  元々の雨粒の混合比以上に蒸発が生じないように上限値を設定
      !
      xyzf_Rain2GasNH4SH(:,:,:,IdxNH4SHr) = - min( DelTime * 4.85d-2 * FactorJ * xyz_NonSaturate * (xyzf_QMixAll(:,:,:,IdxNH4SHr) * xyz_DensBZ) ** 0.65d0, xyzf_QMixAll(:,:,:,IdxNH4SHr) ) 
     
      !蒸気の変換量
      !  雨粒の変換量とは符号が逆となる
      !
      xyzf_Rain2GasNH4SH(:,:,:,IdxNH3) = - xyzf_Rain2GasNH4SH(:,:,:,IdxNH4SHr) * MolWtWet(IdxNH3) / MolWtWet(IdxNH4SHr)
      xyzf_Rain2GasNH4SH(:,:,:,IdxH2S) = - xyzf_Rain2GasNH4SH(:,:,:,IdxNH4SHr) * MolWtWet(IdxH2S) / MolWtWet(IdxNH4SHr)

      xyz_DelPTempNH4SH = ReactHeatNH4SH * xyzf_Rain2GasNH4SH(:,:,:,IdxNH4SHr) / (xyz_ExnerAll * CpDry)

    end if

    !変化量を足し算
    !
    xyzf_DelQMix = xyzf_Rain2Gas + xyzf_Rain2GasNH4SH
    xyz_DelPTemp = sum(xyzf_DelPTemp, 4) + xyz_DelPTempNH4SH

    ! 温位と混合比の計算. 雨から蒸気への変換分を追加
    !
    xyz_PTempAl = xyz_PTempWork + xyz_DelPTemp
    xyzf_QMixAl = xyzf_QMixWork + xyzf_DelQMix

    ! Set Margin
    !
    call SetMargin_xyz(xyz_PTempAl)
    call SetMargin_xyzf(xyzf_QMixAl)

    !------------------------------------------
    ! Output
    !
    xyz_Del  = (xyz_PTempAl - xyz_PTempOrig) / DelTime
    xyzf_Del = (xyzf_QMixAl - xyzf_QMixOrig) / DelTime

    call HistoryAutoPut(TimeN, 'PTempCond', xyz_Del(1:nx, 1:ny, 1:nz) / DelTime)
    do l = 1, ncmax
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(l))//'_Cond', xyzf_Del(1:nx, 1:ny, 1:nz, l) / DelTime)
    end do

  end subroutine Cloudphys_K1969_forcing