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)  :: xyrf_FallRainFlux(imin:imax,jmin:jmax,kmin:kmax, ncmax) 
                                                 !雨粒落下フラックス
    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

    xyrf_FallRainFlux = 0.0d0
    do s = 1, RainNum
      ! 雨粒終端速度
      xyz_VelZRain = - 12.2d0 * FactorJ * ( xyzf_QMixAll(:,:,:,IdxR(s)) ** 0.125d0 )

      xyrf_FallRainFlux(:,:,:,IdxR(s)) = xyr_avr_xyz( xyz_DensBZ * xyzf_QMixAll(:,:,:,IdxR(s)) * xyz_VelZRain )
    end do
    ! 上端のフラックスはゼロ
    do s = 1, RainNum
      xyrf_FallRainFlux(:,:,nz,IdxR(s)) = 0.0d0
    end do

    ! 雨粒落下による時間変化率
    xyzf_FallRain = 0.0d0
    do s = 1, RainNum
      xyzf_FallRain(:,:,:,IdxR(s)) = - xyz_dz_xyr( xyrf_FallRainFlux(:,:,:,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))
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(l))//'_FallFluxAtLB', xyrf_FallRainFlux(1:nx, 1:ny, 0, l))
    end do

    ! SetMargin
    !
!    call SetMargin_xyzf(xyzf_DQMixDt)

  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", "Cloudphys_K1969_Init", "Planet = %c",  c1=trim(Planet))
      call MessageNotify( "M", "Cloudphys_K1969_Init", "FactorJ = %f",  d=(/FactorJ/) )
      call MessageNotify( "M", "Cloudphys_K1969_Init", "AutoConvTime = %f",  d=(/AutoConvTime/) )
      call MessageNotify( "M", "Cloudphys_K1969_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='float')

    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='float')

      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='float')

      call HistoryAutoAddVariable( varname=trim(SpcWetSymbol(l))//'_FallFluxAtLB', dims=(/'x','y','t'/), longname='Falling Rain Flux ' //trim(SpcWetSymbol(l)), units='kg.m-2.s-1', xtype='float')
    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
    integer                        :: iG, iC, iR

    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)
    real(DP)             :: xyz_PressDry(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    

    ! 全エクスナー関数・全圧を計算. サブルーチン内では変化しない.
    !
    xyz_PressAll = PressBasis * ((xyz_ExnerNl + xyz_ExnerBZ) ** (CpDry / GasRDry))
    xyz_ExnerAll = xyz_ExnerNl + xyz_ExnerBZ

    !------------------------------------------    
    ! 暖かい雨のパラメタリゼーション.
    ! * 雲<-->雨 の変換を行う.
    !
    ! 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

      ! 値を保管
      !
      iC = IdxC(s)
      iR = IdxR(s)

      !併合成長
      !
      xyz_AutoConv = DelTime / AutoConvTime * max( 1.0d-60, ( xyzf_QMixAll(:,:,:,iC) - QMixCr) )

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

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

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

    ! Set Margin
    !
!    call SetMargin_xyzf(xyzf_QMixAl)

    !-------------------------------------------
    ! 湿潤飽和調節
    ! * 蒸気<-->雲の変換を行う.
    !
    ! Moist adjustment.
    ! * Conversion from vapor to cloud.
    !
    xyzf_QMixAll = max( 1.0d-60, xyzf_QMixAl + xyzf_QMixBZ )
    xyz_PressDry = xyz_PressDry_xyzf_xyz( xyzf_QMixAll, xyz_PressAll )

    call MoistAdjustSvapPress( xyz_PressDry, xyz_ExnerNl, xyz_PTempAl, xyzf_QMixAl )
    if (IdxNH4SHr /= 0) then 
      call MoistAdjustNH4SH( xyz_PressDry, 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_TempAll   = ( xyz_PTempAl + xyz_PTempBZ ) * ( xyz_ExnerNl + xyz_ExnerBZ )
    xyzf_QMixAll  = max( 1.0d-60, xyzf_QMixAl + xyzf_QMixBZ )
    xyz_PressDry  = xyz_PressDry_xyzf_xyz( xyzf_QMixAll, xyz_PressAll )

    xyzf_Rain2Gas = 0.0d0
    xyzf_DelPTemp = 0.0d0
    xyzf_Rain2GasNH4SH = 0.0d0
    xyz_DelPTempNH4SH  = 0.0d0

    do s = 1, CondNum

       ! 値を保管
       !
       iG = IdxCG(s)
       iC = IdxCC(s)
       iR = IdxCR(s)

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

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

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

    end do

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

    if (IdxNH4SHr /= 0) then 
      xyz_NonSaturate = max( 1.0d-60, - xyz_DelQMixNH4SH( xyz_TempAll, xyz_PressDry, 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

    !------------------------------------------
    ! 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

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

  end subroutine Cloudphys_K1969_forcing