subroutine CloudSimple( xyr_Press, xyz_Press, xyz_Temp, xyz_QH2OVap, xyz_QH2OLiq, xy_Rain, xy_Snow )

    ! USE statements
    !

    ! 時刻管理
    ! Time control
    !
    use timeset, only: DelTime            ! $ \Delta t $ [s]


    real(DP), intent(in   ) :: xyr_Press        ( 0:imax-1, 1:jmax, 0:kmax )
    real(DP), intent(in   ) :: xyz_Press        ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(inout) :: xyz_Temp         ( 0:imax-1, 1:jmax, 1:kmax )
!!$    real(DP), intent(in   ) :: xyz_DQH2OLiqDtCum( 0:imax-1, 1:jmax, 1:kmax )
!!$    real(DP), intent(in   ) :: xyz_DQH2OLiqDtLSC( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(inout) :: xyz_QH2OVap      ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(inout) :: xyz_QH2OLiq      ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(out  ) :: xy_Rain          ( 0:imax-1, 1:jmax )
    real(DP), intent(out  ) :: xy_Snow          ( 0:imax-1, 1:jmax )


    ! Tentative
    real(DP) xyz_DQH2OLiqDtCum( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) xyz_DQH2OLiqDtLSC( 0:imax-1, 1:jmax, 1:kmax )

    real(DP) :: xyz_TempB   ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_QH2OVapB( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_QH2OLiqB( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_DQRainDt( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_DQSnowDt( 0:imax-1, 1:jmax, 1:kmax )

    real(DP) :: xyz_QH2OSolB( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_QH2OSol ( 0:imax-1, 1:jmax, 1:kmax )

!!$    real(DP) :: xyz_DTempDtPrcpPCCum( 0:imax-1, 1:jmax, 1:kmax )
!!$    real(DP) :: xyz_DTempDtPrcpPCLsc( 0:imax-1, 1:jmax, 1:kmax )

    ! 実行文 ; Executable statement
    !

    ! 初期化確認
    ! Initialization check
    !
    if ( .not. cloud_simple_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    ! tentative treatment
    xyz_DQH2OLiqDtCum = 0.0_DP
    xyz_DQH2OLiqDtLSC = 0.0_DP


    ! Numerical solution

!!$      xyz_DQCloudWaterDt = xyz_DQCloudWaterDt &
!!$        & - xyz_QCloudWater / ( CloudLifeTime + 1.0d-100 )


!      ( X_{t+1} - X_{t-1} ) / ( 2 \Delta t ) = Q - X_{t+1} / \tau
!
!      X_{t+1} / ( 2 \Delta t )  + X_{t+1} / \tau = X_{t-1} / ( 2 \Delta t ) + Q
!      ( 1 / ( 2 \Delta t )  + 1 / \tau ) X_{t+1} = X_{t-1} / ( 2 \Delta t ) + Q
!      X_{t+1} = ( X_{t-1} / ( 2 \Delta t ) + Q ) / ( 1 / ( 2 \Delta t )  + 1 / \tau ) 

!!$    xyz_QH2OLiq =                                                           &
!!$      &   (                                                                 &
!!$      &       xyz_QH2OLiq / ( 2.0_DP * DelTime )                            &
!!$      &     + xyz_DQH2OLiqDtCum + xyz_DQH2OLiqDtLSC                         &
!!$      &   )                                                                 &
!!$      & / ( 1.0_DP / ( 2.0_DP * DelTime ) + 1.0_DP / ( CloudLifeTime + 1.0d-100 ) )
!!$
!!$    call CloudUtilsCalcPRCPKeyLLTemp3D(          &
!!$      & xyr_Press, xyz_Temp, xyz_DQH2OLiqDtCum,  &  ! (in )
!!$      & xy_RainCum, xy_SnowCum                   &  ! (out)
!!$      & )
!!$    call CloudUtilsCalcPRCPKeyLLTemp3D(          &
!!$      & xyr_Press, xyz_Temp, xyz_DQH2OLiqDtLsc,  &  ! (in )
!!$      & xy_RainLsc, xy_SnowLsc                   &  ! (out)
!!$      & )
    !-----


    ! Analytical solution

    ! save values before adjustment
    xyz_TempB    = xyz_Temp
    xyz_QH2OVapB = xyz_QH2OVap
    xyz_QH2OLiqB = xyz_QH2OLiq

    xyz_QH2OLiq = xyz_QH2OLiq * exp( - 2.0_DP * DelTime / ( CloudLifeTime + 1.0e-100_DP ) ) + ( xyz_DQH2OLiqDtCum + xyz_DQH2OLiqDtLSC ) * CloudLifeTime * ( 1.0_DP - exp( - 2.0_DP * DelTime / ( CloudLifeTime + 1.0e-100_DP ) ) )

    xyz_DQRainDt = xyz_QH2OLiqB + ( xyz_DQH2OLiqDtCum + xyz_DQH2OLiqDtLSC ) * 2.0_DP * DelTime - xyz_QH2OLiq
    xyz_DQRainDt = xyz_DQRainDt / ( 2.0_DP * DelTime )



    select case ( IDSnowMethod )
    case ( IDSnowMethodKeyLLTemp )

      call CloudSimpleCalcPRCPKeyLLTemp3D( xyr_Press, xyz_Press, xyz_DQRainDt, xyz_Temp, xyz_QH2OVap, xy_Rain, xy_Snow )

      xyz_QH2OSolB = 0.0_DP
      xyz_QH2OSol  = 0.0_DP
      call CloudSimpleConsChk( .false., xyr_Press, xyz_TempB, xyz_QH2OVapB, xyz_QH2OLiqB, xyz_QH2OSolB, xyz_Temp , xyz_QH2OVap , xyz_QH2OLiq , xyz_QH2OSol , xy_Rain, xy_Snow )

    case ( IDSnowMethodStepPC )

      xyz_DQSnowDt = 0.0_DP

      call CloudSimpleCalcPRCPStepPC( xyr_Press, xyz_Press, xyz_DQRainDt, xyz_DQSnowDt, xyz_Temp, xyz_QH2OVap, xy_Rain, xy_Snow )

      xyz_QH2OSolB = 0.0_DP
      xyz_QH2OSol  = 0.0_DP
      call CloudSimpleConsChk( .true., xyr_Press, xyz_TempB, xyz_QH2OVapB, xyz_QH2OLiqB, xyz_QH2OSolB, xyz_Temp , xyz_QH2OVap , xyz_QH2OLiq , xyz_QH2OSol , xy_Rain, xy_Snow )

    end select


  end subroutine CloudSimple

  subroutine CloudSimpleCalcCloudCover( xyz_Press, xyz_Temp, xyz_QH2OTot, xyz_CloudCover )

    ! USE statements
    !

    ! 飽和比湿の算出
    ! Evaluate saturation specific humidity
    !
    use saturate, only: xyz_CalcQVapSat

    real(DP), intent(in ) :: xyz_Press     ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(in ) :: xyz_Temp      ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(in ) :: xyz_QH2OTot   ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(out) :: xyz_CloudCover( 0:imax-1, 1:jmax, 1:kmax )


    real(DP) :: xyz_RH(0:imax-1, 1:jmax, 1:kmax)

    integer :: i
    integer :: j
    integer :: k


    ! 実行文 ; Executable statement
    !

    ! 初期化確認
    ! Initialization check
    !
    if ( .not. cloud_simple_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    select case ( IDCloudCoverMethod )
    case ( IDCloudCoverMethodConst )

      xyz_CloudCover = CloudCover

    case ( IDCloudCoverMethodRH )

      ! see Sundqvist et al. (1989), Del Genio et al. (1996)
      xyz_RH = xyz_QH2OTot / xyz_CalcQVapSat( xyz_Temp, xyz_Press )
      xyz_RH = min( xyz_RH, 1.0_DP )

      xyz_CloudCover = 1.0_DP - sqrt( ( 1.0_DP - xyz_RH ) / ( 1.0_DP - CloudCoverRHCrtl ) )

      xyz_CloudCover = max( xyz_CloudCover, CloudCoverMin )
      xyz_CloudCover = min( xyz_CloudCover, 1.0_DP )

    case ( IDCloudCoverMethodRHLin )

      xyz_RH = xyz_QH2OTot / xyz_CalcQVapSat( xyz_Temp, xyz_Press )
      xyz_RH = min( xyz_RH, 1.0_DP )

      if ( CloudCoverRHCrtl < 1.0_DP ) then
!!$      xyz_CloudCover = 2.0_DP * xyz_RH - 1.0_DP
        xyz_CloudCover = xyz_RH           / ( 1.0_DP - CloudCoverRHCrtl ) - CloudCoverRHCrtl / ( 1.0_DP - CloudCoverRHCrtl )
      else
        do k = 1, kmax
          do j = 1, jmax
            do i = 0, imax-1
              if ( xyz_RH(i,j,k) >= 1.0_DP ) then
                xyz_CloudCover(i,j,k) = 1.0_DP
              else
                xyz_CloudCover(i,j,k) = 0.0_DP
              end if
            end do
          end do
        end do
      end if

      xyz_CloudCover = max( xyz_CloudCover, CloudCoverMin )
      xyz_CloudCover = min( xyz_CloudCover, 1.0_DP )

    end select


  end subroutine CloudSimpleCalcCloudCover

  subroutine CloudSimpleCalcPRCPKeyLLTemp( xyz_Temp, xy_PRCP, xy_SurfRainFlux, xy_SurfSnowFlux )


    ! 雪と海氷の定数の設定
    ! Setting constants of snow and sea ice
    !
    use constants_snowseaice, only: TempCondWater


    real(DP), intent(in ) :: xyz_Temp       ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(in ) :: xy_PRCP        ( 0:imax-1, 1:jmax )
    real(DP), intent(out) :: xy_SurfRainFlux( 0:imax-1, 1:jmax )
    real(DP), intent(out) :: xy_SurfSnowFlux( 0:imax-1, 1:jmax )


    ! 作業変数
    ! Work variables
    !
    integer:: i               ! 経度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in longitude
    integer:: j               ! 緯度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in latitude


    ! 初期化確認
    ! Initialization check
    !
    if ( .not. cloud_simple_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    if ( FlagSnow ) then

      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_Temp(i,j,1) > TempCondWater ) then
            xy_SurfRainFlux(i,j) = xy_PRCP(i,j)
            xy_SurfSnowFlux(i,j) = 0.0_DP
          else
            xy_SurfRainFlux(i,j) = 0.0_DP
            xy_SurfSnowFlux(i,j) = xy_PRCP(i,j)
          end if
        end do
      end do

    else

      xy_SurfRainFlux = xy_PRCP
      xy_SurfSnowFlux = 0.0_DP

    end if


  end subroutine CloudSimpleCalcPRCPKeyLLTemp

  subroutine CloudSimpleCalcPRCPKeyLLTemp3D( xyr_Press, xyz_Press, xyz_DQH2OLiqDt, xyz_Temp, xyz_QH2OVap, xy_SurfRainFlux, xy_SurfSnowFlux )

    ! 時刻管理
    ! Time control
    !
    use timeset, only: DelTime            ! $ \Delta t $ [s]

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: CpDry, Grav, LatentHeat, LatentHeatFusion, EpsV
                              ! $ \epsilon_v $ .
                              ! 水蒸気分子量比.
                              ! Molecular weight of water vapor

    ! 飽和比湿の算出
    ! Evaluate saturation specific humidity
    !
    use saturate, only: xyz_CalcQVapSat, xyz_CalcQVapSatOnLiq

    real(DP), intent(in   ) :: xyr_Press       ( 0:imax-1, 1:jmax, 0:kmax )
    real(DP), intent(in   ) :: xyz_Press       ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(in   ) :: xyz_DQH2OLiqDt  ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(inout) :: xyz_Temp        ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(inout) :: xyz_QH2OVap     ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(out  ) :: xy_SurfRainFlux ( 0:imax-1, 1:jmax )
    real(DP), intent(out  ) :: xy_SurfSnowFlux ( 0:imax-1, 1:jmax )


    ! 作業変数
    ! Work variables
    !
    real(DP) :: xyz_DelMass ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_DQRainDt( 0:imax-1, 1:jmax, 1:kmax )

    real(DP) :: VirTemp
    real(DP) :: aaa_QH2OVapSat(1,1,1)
    real(DP) :: QH2OVapSat
    real(DP) :: PRCPFlux
    real(DP) :: DelPRCPFlux
    real(DP) :: DelQH2OVap
    real(DP) :: LatentHeatLocal
    character(STRING) :: CharPhase

    real(DP) :: xy_PRCP( 0:imax-1, 1:jmax )

    integer  :: i
    integer  :: j
    integer  :: k


    ! 初期化確認
    ! Initialization check
    !
    if ( .not. cloud_simple_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    do k = 1, kmax
      xyz_DelMass(:,:,k) = ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav
    end do


    if ( FlagPRCPEvap ) then

      xyz_DQRainDt = xyz_DQH2OLiqDt

      xy_SurfRainFlux = 0.0_DP
      do j = 1, jmax
        do i = 0, imax-1
          do k = kmax, 1, -1

            ! This is moved below.
!!$            xy_SurfRainFlux(i,j) = xy_SurfRainFlux(i,j) &
!!$              & + xyz_DQRainDt(i,j,k) * xyz_DelMass(i,j,k)

            CharPhase = 'liquid'
            aaa_QH2OVapSat(1:1,1:1,1:1) = xyz_CalcQVapSatOnLiq( xyz_Temp(i:i,j:j,k:k), xyz_Press(i:i,j:j,k:k) )
            PRCPFlux = xy_SurfRainFlux(i,j)
            QH2OVapSat = aaa_QH2OVapSat(1,1,1)
            VirTemp = xyz_Temp(i,j,k) * ( 1.0_DP + ((( 1.0_DP / EpsV ) - 1.0_DP ) * xyz_QH2OVap(i,j,k)) )
            call CloudSimpleEvap1Grid( CharPhase, xyz_DelMass(i,j,k), xyz_Press(i,j,k), xyz_QH2OVap(i,j,k), QH2OVapSat, VirTemp, PRCPFlux, DelPRCPFlux )
            xy_SurfRainFlux(i,j) = PRCPFlux - DelPRCPFlux
            LatentHeatLocal      = LatentHeat
            DelQH2OVap = DelPRCPFlux * ( 2.0_DP * DelTime ) / xyz_DelMass(i,j,k)
            xyz_QH2OVap(i,j,k) = xyz_QH2OVap(i,j,k) + DelQH2OVap
            xyz_Temp(i,j,k) = xyz_Temp(i,j,k) - LatentHeatLocal * DelQH2OVap / CpDry

            xy_SurfRainFlux(i,j) = xy_SurfRainFlux(i,j) + xyz_DQRainDt(i,j,k) * xyz_DelMass(i,j,k)

          end do
        end do
      end do

      xy_PRCP = xy_SurfRainFlux

    else

      xy_PRCP = 0.0d0
      do k = kmax, 1, -1
        xy_PRCP = xy_PRCP + xyz_DQH2OLiqDt(:,:,k) * ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav
      end do

    end if


    call CloudSimpleCalcPRCPKeyLLTemp( xyz_Temp, xy_PRCP, xy_SurfRainFlux, xy_SurfSnowFlux )


  end subroutine CloudSimpleCalcPRCPKeyLLTemp3D

  subroutine CloudSimpleCalcPRCPStepPC( xyr_Press, xyz_Press, xyz_DQRainDt, xyz_DQSnowDt, xyz_Temp, xyz_QH2OVap, xy_SurfRainFlux, xy_SurfSnowFlux )

    ! 時刻管理
    ! Time control
    !
    use timeset, only: DelTime            ! $ \Delta t $ [s]

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: CpDry, Grav, LatentHeat, LatentHeatFusion, EpsV
                              ! $ \epsilon_v $ .
                              ! 水蒸気分子量比.
                              ! Molecular weight of water vapor

    ! 雪と海氷の定数の設定
    ! Setting constants of snow and sea ice
    !
    use constants_snowseaice, only: TempCondWater

    ! 飽和比湿の算出
    ! Evaluate saturation specific humidity
    !
    use saturate, only: xyz_CalcQVapSat, xyz_CalcQVapSatOnLiq, xyz_CalcQVapSatOnSol


    real(DP), intent(in   ) :: xyr_Press       ( 0:imax-1, 1:jmax, 0:kmax )
    real(DP), intent(in   ) :: xyz_Press       ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(in   ) :: xyz_DQRainDt    ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(in   ) :: xyz_DQSnowDt    ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(inout) :: xyz_Temp        ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(inout) :: xyz_QH2OVap     ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(out  ) :: xy_SurfRainFlux ( 0:imax-1, 1:jmax )
    real(DP), intent(out  ) :: xy_SurfSnowFlux ( 0:imax-1, 1:jmax )


    ! 作業変数
    ! Work variables
    !
    real(DP) :: xyz_DelMass( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: MassMaxFreezeRate
    real(DP) :: MassFreezeRate
    real(DP) :: MassMaxMeltRate
    real(DP) :: MassMeltRate

    real(DP) :: VirTemp
    real(DP) :: aaa_QH2OVapSat(1,1,1)
    real(DP) :: QH2OVapSat
    real(DP) :: PRCPFlux
    real(DP) :: DelPRCPFlux
    real(DP) :: DelQH2OVap
    real(DP) :: LatentHeatLocal
    character(STRING) :: CharPhase

    integer  :: i
    integer  :: j
    integer  :: k
    integer  :: l


    ! 初期化確認
    ! Initialization check
    !
    if ( .not. cloud_simple_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    do k = 1, kmax
      xyz_DelMass(:,:,k) = ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav
    end do

    ! Freezing and melting switching at temperature of TempCondWater
    xy_SurfRainFlux = 0.0_DP
    xy_SurfSnowFlux = 0.0_DP
    do j = 1, jmax
      do i = 0, imax-1
        do k = kmax, 1, -1

          ! These are moved below.
!!$          xy_SurfRainFlux(i,j) = xy_SurfRainFlux(i,j) &
!!$            & + xyz_DQRainDt(i,j,k) * xyz_DelMass(i,j,k)
!!$          xy_SurfSnowFlux(i,j) = xy_SurfSnowFlux(i,j) &
!!$            & + xyz_DQSnowDt(i,j,k) * xyz_DelMass(i,j,k)

          MassMaxFreezeRate = CpDry * ( TempCondWater - xyz_Temp(i,j,k) ) * xyz_DelMass(i,j,k) / LatentHeatFusion / ( 2.0_DP * DelTime )
          if ( MassMaxFreezeRate >= 0.0_DP ) then
            ! freezing
            if ( xy_SurfRainFlux(i,j) >= MassMaxFreezeRate ) then
              MassFreezeRate = MassMaxFreezeRate
            else
              MassFreezeRate = xy_SurfRainFlux(i,j)
            end if
            xy_SurfRainFlux(i,j) = xy_SurfRainFlux(i,j) - MassFreezeRate
            xy_SurfSnowFlux(i,j) = xy_SurfSnowFlux(i,j) + MassFreezeRate
            xyz_Temp(i,j,k) = xyz_Temp(i,j,k) + LatentHeatFusion * MassFreezeRate * 2.0_DP * DelTime / ( CpDry * xyz_DelMass(i,j,k) )
          else
            ! melting
            MassMaxMeltRate = - MassMaxFreezeRate
            if ( xy_SurfSnowFlux(i,j) >= MassMaxMeltRate ) then
              MassMeltRate = MassMaxMeltRate
            else
              MassMeltRate = xy_SurfSnowFlux(i,j)
            end if
            xy_SurfRainFlux(i,j) = xy_SurfRainFlux(i,j) + MassMeltRate
            xy_SurfSnowFlux(i,j) = xy_SurfSnowFlux(i,j) - MassMeltRate
            xyz_Temp(i,j,k) = xyz_Temp(i,j,k) - LatentHeatFusion * MassMeltRate * 2.0_DP * DelTime / ( CpDry * xyz_DelMass(i,j,k) )
          end if


          if ( FlagPRCPEvap ) then
!!$            do l = 0, 0   ! for test
            do l = 1, 2
              select case ( l )
              case ( 0 ) ! mixture
                CharPhase = 'mixture'
                aaa_QH2OVapSat(1:1,1:1,1:1) = xyz_CalcQVapSat     ( xyz_Temp(i:i,j:j,k:k), xyz_Press(i:i,j:j,k:k) )
                PRCPFlux = xy_SurfRainFlux(i,j)
              case ( 1 ) ! liquid
                CharPhase = 'liquid'
                aaa_QH2OVapSat(1:1,1:1,1:1) = xyz_CalcQVapSatOnLiq( xyz_Temp(i:i,j:j,k:k), xyz_Press(i:i,j:j,k:k) )
                PRCPFlux = xy_SurfRainFlux(i,j)
              case ( 2 ) ! solid
                CharPhase = 'solid'
                aaa_QH2OVapSat(1:1,1:1,1:1) = xyz_CalcQVapSatOnSol( xyz_Temp(i:i,j:j,k:k), xyz_Press(i:i,j:j,k:k) )
                PRCPFlux = xy_SurfSnowFlux(i,j)
              case default
                call MessageNotify( 'E', module_name, 'Unexpected number for select case.' )
              end select
              QH2OVapSat = aaa_QH2OVapSat(1,1,1)
              VirTemp = xyz_Temp(i,j,k) * ( 1.0_DP + ((( 1.0_DP / EpsV ) - 1.0_DP ) * xyz_QH2OVap(i,j,k)) )
              call CloudSimpleEvap1Grid( CharPhase, xyz_DelMass(i,j,k), xyz_Press(i,j,k), xyz_QH2OVap(i,j,k), QH2OVapSat, VirTemp, PRCPFlux, DelPRCPFlux )
              select case ( l )
              case ( 0 ) ! mixture
                xy_SurfRainFlux(i,j) = PRCPFlux - DelPRCPFlux
                LatentHeatLocal      = LatentHeat
              case ( 1 ) ! liquid
                xy_SurfRainFlux(i,j) = PRCPFlux - DelPRCPFlux
                LatentHeatLocal      = LatentHeat
              case ( 2 ) ! solid
                xy_SurfSnowFlux(i,j) = PRCPFlux - DelPRCPFlux
                LatentHeatLocal      = LatentHeat + LatentHeatFusion
              end select
              DelQH2OVap = DelPRCPFlux * ( 2.0_DP * DelTime ) / xyz_DelMass(i,j,k)
              xyz_QH2OVap(i,j,k) = xyz_QH2OVap(i,j,k) + DelQH2OVap
              xyz_Temp(i,j,k) = xyz_Temp(i,j,k) - LatentHeatLocal * DelQH2OVap / CpDry
            end do
          end if

          xy_SurfRainFlux(i,j) = xy_SurfRainFlux(i,j) + xyz_DQRainDt(i,j,k) * xyz_DelMass(i,j,k)
          xy_SurfSnowFlux(i,j) = xy_SurfSnowFlux(i,j) + xyz_DQSnowDt(i,j,k) * xyz_DelMass(i,j,k)

        end do
      end do
    end do


  end subroutine CloudSimpleCalcPRCPStepPC

  subroutine CloudSimpleDivideWatAndIce( xyz_Temp, xyz_QH2OWatAndIce, xyz_QH2OWat, xyz_QH2OIce )

    ! USE statements
    !

    ! 飽和比湿の算出
    ! Evaluate saturation specific humidity
    !
    use saturate, only : SaturateWatFraction

    real(DP), intent(in ) :: xyz_Temp         (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ) :: xyz_QH2OWatAndIce(0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(out) :: xyz_QH2OWat      (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(out) :: xyz_QH2OIce      (0:imax-1, 1:jmax, 1:kmax)


    real(DP) :: xyz_WatFrac(0:imax-1, 1:jmax, 1:kmax)

    ! 実行文 ; Executable statement
    !

    ! 初期化確認
    ! Initialization check
    !
    if ( .not. cloud_simple_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    call SaturateWatFraction( xyz_Temp, xyz_WatFrac )

    xyz_QH2OWat = xyz_QH2OWatAndIce * xyz_WatFrac
    xyz_QH2OIce = xyz_QH2OWatAndIce * ( 1.0_DP - xyz_WatFrac )


  end subroutine CloudSimpleDivideWatAndIce

  subroutine CloudSimpleInit( ArgFlagSnow )

    ! ファイル入出力補助
    ! File I/O support
    !
    use dc_iounit, only: FileOpen

    ! NAMELIST ファイル入力に関するユーティリティ
    ! Utilities for NAMELIST file input
    !
    use namelist_util, only: namelist_filename, NmlutilMsg, NmlutilAryValid

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoAddVariable

    ! 飽和比湿の算出
    ! Evaluate saturation specific humidity
    !
    use saturate, only: SaturateInit

    ! 大規模凝結 (非対流性凝結)
    ! Large scale condensation (non-convective condensation)
    !
    use lscond, only : LScaleCondInit

    ! 宣言文 ; Declaration statements
    !

    logical, intent(in) :: ArgFlagSnow


    character(STRING) :: CloudCoverMethod
    character(STRING) :: SnowMethod

    integer:: unit_nml        ! NAMELIST ファイルオープン用装置番号.
                              ! Unit number for NAMELIST file open
    integer:: iostat_nml      ! NAMELIST 読み込み時の IOSTAT.
                              ! IOSTAT of NAMELIST read

    ! NAMELIST 変数群
    ! NAMELIST group name
    !
    namelist /cloud_simple_nml/ CloudLifeTime, CloudWatLifeTime, CloudIceLifeTime, CloudCoverMethod, SnowMethod, CloudCover, CloudCoverRHCrtl, CloudCoverMin, FlagPRCPEvap, PRCPArea, PRCPEvapArea
          !
          ! デフォルト値については初期化手続 "cloud_simple#CloudSimpleInit"
          ! のソースコードを参照のこと.
          !
          ! Refer to source codes in the initialization procedure
          ! "cloud_simple#CloudSimpleInit" for the default values.
          !

    ! 実行文 ; Executable statement
    !

    if ( cloud_simple_inited ) return


    FlagSnow = ArgFlagSnow


    ! デフォルト値の設定
    ! Default values settings
    !
    CloudLifeTime       = 3600.0_DP
    CloudWatLifeTime    = 3600.0_DP
    CloudIceLifeTime    = 3600.0_DP

    CloudCoverMethod    = 'Const'
!!$    CloudCoverMethod    = 'RH'

    SnowMethod          = 'KeyLLTemp'

    CloudCover          = 1.0_DP

    CloudCoverRHCrtl    = 0.0_DP

    CloudCoverMin       = 0.0_DP

    FlagPRCPEvap        = .false.
!!$    PRCPEvapArea        = 0.5_DP
    PRCPArea            = 1.0_DP
!!$    PRCPArea            = 0.5_DP
    PRCPEvapArea        = 1.0_DP
!!$    PRCPEvapArea        = 0.5_DP


    ! NAMELIST の読み込み
    ! NAMELIST is input
    !
    if ( trim(namelist_filename) /= '' ) then
      call FileOpen( unit_nml, namelist_filename, mode = 'r' ) ! (in)

      rewind( unit_nml )
      read( unit_nml, nml = cloud_simple_nml, iostat = iostat_nml )             ! (out)
      close( unit_nml )

      call NmlutilMsg( iostat_nml, module_name ) ! (in)
    end if


    select case ( CloudCoverMethod )
    case ( 'Const' )
      IDCloudCoverMethod = IDCloudCoverMethodConst
    case ( 'RH' )
      IDCloudCoverMethod = IDCloudCoverMethodRH
    case ( 'RHLin' )
      IDCloudCoverMethod = IDCloudCoverMethodRHLin
    case default
      call MessageNotify( 'E', module_name, 'CloudCoverMethod=<%c> is not supported.', c1 = trim(CloudCoverMethod) )
    end select


    select case ( SnowMethod )
    case ( 'KeyLLTemp' )
      IDSnowMethod = IDSnowMethodKeyLLTemp
    case ( 'StepPC' )
      IDSnowMethod = IDSnowMethodStepPC
    case default
      call MessageNotify( 'E', module_name, 'SnowMethod=<%c> is not supported.', c1 = trim(SnowMethod) )
    end select



    ! Initialization of modules used in this module
    !

    ! 飽和比湿の算出
    ! Evaluate saturation specific humidity
    !
    call SaturateInit

    ! 大規模凝結 (非対流性凝結) (Manabe, 1965)
    ! Large scale condensation (non-convective condensation) (Le Treut and Li, 1991)
    !
    call LScaleCondInit( FlagSnow )


    ! ヒストリデータ出力のためのへの変数登録
    ! Register of variables for history data output
    !
!!$    call HistoryAutoAddVariable( 'EffCloudCover', &
!!$      & (/ 'lon ', 'lat ', 'time' /), &
!!$      & 'effective cloud cover', '1' )



    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, 'CloudLifeTime       = %f', d = (/ CloudLifeTime /) )
    call MessageNotify( 'M', module_name, 'CloudWatLifeTime    = %f', d = (/ CloudWatLifeTime /) )
    call MessageNotify( 'M', module_name, 'CloudIceLifeTime    = %f', d = (/ CloudIceLifeTime /) )
    call MessageNotify( 'M', module_name, 'CloudCoverMethod    = %c', c1 = trim(CloudCoverMethod) )
    call MessageNotify( 'M', module_name, 'SnowMethod          = %c', c1 = trim(SnowMethod) )
    call MessageNotify( 'M', module_name, 'CloudCover          = %f', d = (/ CloudCover /) )
    call MessageNotify( 'M', module_name, 'CloudCoverRHCrtl    = %f', d = (/ CloudCoverRHCrtl /) )
    call MessageNotify( 'M', module_name, 'CloudCoverMin       = %f', d = (/ CloudCoverMin /) )
    call MessageNotify( 'M', module_name, 'FlagPRCPEvap        = %b', l = (/ FlagPRCPEvap /) )
    call MessageNotify( 'M', module_name, 'PRCPArea            = %f', d = (/ PRCPArea /) )
    call MessageNotify( 'M', module_name, 'PRCPEvapArea        = %f', d = (/ PRCPEvapArea /) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )


    cloud_simple_inited = .true.

  end subroutine CloudSimpleInit

  subroutine CloudSimpleWithIce( xyr_Press, xyz_Press, xyz_Temp, xyz_QH2OVap, xyz_QH2OLiq, xyz_QH2OSol, xy_Rain, xy_Snow )

    ! USE statements
    !

    ! 時刻管理
    ! Time control
    !
    use timeset, only: DelTime            ! $ \Delta t $ [s]

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: Grav
                              ! $ g $ [m s-2].
                              ! 重力加速度.
                              ! Gravitational acceleration

    ! 雲関系ルーチン
    ! Cloud-related routines
    !
    use cloud_utils, only : CloudUtilsPRCPStepPC1Grid, CloudUtilsPRCPEvap1Grid, CloudUtilConsChk


    real(DP), intent(in   ) :: xyr_Press        ( 0:imax-1, 1:jmax, 0:kmax )
    real(DP), intent(in   ) :: xyz_Press        ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(inout) :: xyz_Temp         ( 0:imax-1, 1:jmax, 1:kmax )
!!$    real(DP), intent(in   ) :: xyz_DQH2OLiqDtCum( 0:imax-1, 1:jmax, 1:kmax )
!!$    real(DP), intent(in   ) :: xyz_DQH2OLiqDtLSC( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(inout) :: xyz_QH2OVap     ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(inout) :: xyz_QH2OLiq      ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(inout) :: xyz_QH2OSol      ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(out  ) :: xy_Rain          ( 0:imax-1, 1:jmax )
    real(DP), intent(out  ) :: xy_Snow          ( 0:imax-1, 1:jmax )


    real(DP) :: xyz_TempB   ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_QH2OVapB( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_QH2OLiqB( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_QH2OSolB( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_DQH2OLiqDt( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_DQH2OSolDt( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xy_DQRainDt( 0:imax-1, 1:jmax )
    real(DP) :: xy_DQSnowDt( 0:imax-1, 1:jmax )

    real(DP) :: xyz_DelMass( 0:imax-1, 1:jmax, 1:kmax )

    integer :: i
    integer :: j
    integer :: k


    ! 実行文 ; Executable statement
    !

    ! 初期化確認
    ! Initialization check
    !
    if ( .not. cloud_simple_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if

    ! save values before adjustment
    xyz_TempB    = xyz_Temp
    xyz_QH2OVapB = xyz_QH2OVap
    xyz_QH2OLiqB = xyz_QH2OLiq
    xyz_QH2OSolB = xyz_QH2OSol


    xyz_DQH2OLiqDt = 0.0_DP
    xyz_DQH2OSolDt = 0.0_DP

    do k = 1, kmax
      xyz_DelMass(:,:,k) = ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav
    end do


    ! Rain and snow at the surface
    xy_Rain = 0.0_DP
    xy_Snow = 0.0_DP

    k_loop : do k = kmax, 1, -1

      ! Freezing/melting and evaporation of precipitation
      !
      do j = 1, jmax
        do i = 0, imax-1
          call CloudUtilsPRCPStepPC1Grid( xyr_Press(i,j,k-1), xyr_Press(i,j,k), xyz_Temp(i,j,k), xy_Rain(i,j), xy_Snow(i,j) )
        end do
      end do
      do j = 1, jmax
        do i = 0, imax-1
          call CloudUtilsPRCPEvap1Grid( xyz_Press(i,j,k), xyr_Press(i,j,k-1), xyr_Press(i,j,k), PRCPArea, PRCPEvapArea, xyz_Temp(i,j,k), xyz_QH2OVap(i,j,k), xy_Rain(i,j), xy_Snow(i,j) )
        end do
      end do


      xyz_QH2OLiq(:,:,k) = xyz_QH2OLiq(:,:,k) * exp( - 2.0_DP * DelTime / ( CloudWatLifeTime + 1.0e-100_DP ) ) + xyz_DQH2OLiqDt(:,:,k) * CloudWatLifeTime * ( 1.0_DP - exp( - 2.0_DP * DelTime / ( CloudWatLifeTime + 1.0e-100_DP ) ) )

      xy_DQRainDt = xyz_QH2OLiqB(:,:,k) + xyz_DQH2OLiqDt(:,:,k) * 2.0_DP * DelTime - xyz_QH2OLiq(:,:,k)
      xy_DQRainDt = xy_DQRainDt / ( 2.0_DP * DelTime )


      xyz_QH2OSol(:,:,k) = xyz_QH2OSol(:,:,k) * exp( - 2.0_DP * DelTime / ( CloudIceLifeTime + 1.0e-100_DP ) ) + xyz_DQH2OSolDt(:,:,k) * CloudIceLifeTime * ( 1.0_DP - exp( - 2.0_DP * DelTime / ( CloudIceLifeTime + 1.0e-100_DP ) ) )

      xy_DQSnowDt = xyz_QH2OSolB(:,:,k) + xyz_DQH2OSolDt(:,:,k) * 2.0_DP * DelTime - xyz_QH2OSol(:,:,k)
      xy_DQSnowDt = xy_DQSnowDt / ( 2.0_DP * DelTime )

      xy_Rain = xy_Rain + xy_DQRainDt * xyz_DelMass(:,:,k)
      xy_Snow = xy_Snow + xy_DQSnowDt * xyz_DelMass(:,:,k)


    end do k_loop


    call CloudUtilConsChk( "CloudSimpleWithIce", xyr_Press, xyz_TempB, xyz_QH2OVapB, xyz_QH2OLiqB, xyz_QH2OSolB, xyz_Temp , xyz_QH2OVap , xyz_QH2OLiq , xyz_QH2OSol , xy_Rain, xy_Snow )


  end subroutine CloudSimpleWithIce

  subroutine CloudSimpleConsChk( FlagIncludeIcePhaseChange, xyr_Press, xyz_TempB, xyz_QH2OVapB, xyz_QH2OLiqB, xyz_QH2OSolB, xyz_Temp , xyz_QH2OVap , xyz_QH2OLiq , xyz_QH2OSol , xy_Rain, xy_Snow )


    ! 時刻管理
    ! Time control
    !
    use timeset, only: DelTime            ! $ \Delta t $ [s]

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: Grav, CpDry, LatentHeat, LatentHeatFusion
                              ! $ L $ [J kg-1] .
                              ! 融解の潜熱.
                              ! Latent heat of fusion

    logical , intent(in) :: FlagIncludeIcePhaseChange
    real(DP), intent(in) :: xyr_Press   (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(in) :: xyz_TempB   (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in) :: xyz_QH2OVapB(0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in) :: xyz_QH2OLiqB(0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in) :: xyz_QH2OSolB(0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in) :: xyz_Temp    (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in) :: xyz_QH2OVap (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in) :: xyz_QH2OLiq (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in) :: xyz_QH2OSol (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in) :: xy_Rain     (0:imax-1, 1:jmax)
    real(DP), intent(in) :: xy_Snow     (0:imax-1, 1:jmax)

    ! Local variables
    !
    real(DP) :: xyz_DelMass(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xy_Val(0:imax-1, 1:jmax)
    real(DP) :: xy_SumB(0:imax-1, 1:jmax)
    real(DP) :: xy_Sum(0:imax-1, 1:jmax)
    real(DP) :: xy_Ratio(0:imax-1, 1:jmax)
    integer  :: i
    integer  :: j
    integer  :: k


    do k = 1, kmax
      xyz_DelMass(:,:,k) = ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav
    end do

    xy_Sum = 0.0_DP
    do k = kmax, 1, -1
      xy_Val =   CpDry * xyz_TempB(:,:,k) + LatentHeat * xyz_QH2OVapB(:,:,k) - LatentHeatFusion * xyz_QH2OSolB(:,:,k)
      xy_Sum = xy_Sum + xy_Val * xyz_DelMass(:,:,k)
    end do

    xy_SumB = xy_Sum

    xy_Sum = 0.0_DP
    do k = kmax, 1, -1
      xy_Val =   CpDry * xyz_Temp (:,:,k) + LatentHeat * xyz_QH2OVap (:,:,k) - LatentHeatFusion * xyz_QH2OSol (:,:,k)
      xy_Sum = xy_Sum + xy_Val * xyz_DelMass(:,:,k)
    end do
    if ( FlagIncludeIcePhaseChange ) then
      xy_Sum = xy_Sum - LatentHeatFusion * xy_Snow * 2.0_DP * DelTime
    end if

    xy_Ratio = ( xy_Sum - xy_SumB ) / ( xy_Sum + 1.0d-100 )
    do j = 1, jmax
      do i = 0, imax-1
        if ( abs( xy_Ratio(i,j) ) > 1.0d-10 ) then
          call MessageNotify( 'M', module_name, 'Modified condensate static energy is not conserved, %f.', d = (/ xy_Ratio(i,j) /) )
        end if
      end do
    end do



    xy_Sum = 0.0_DP
    do k = kmax, 1, -1
      xy_Val = xyz_QH2OVapB(:,:,k) + xyz_QH2OLiqB(:,:,k) + xyz_QH2OSolB(:,:,k)
      xy_Sum = xy_Sum + xy_Val * xyz_DelMass(:,:,k)
    end do

    xy_SumB = xy_Sum

    xy_Sum = 0.0_DP
    do k = kmax, 1, -1
      xy_Val = xyz_QH2OVap (:,:,k) + xyz_QH2OLiq (:,:,k) + xyz_QH2OSol (:,:,k)
      xy_Sum = xy_Sum + xy_Val * xyz_DelMass(:,:,k)
    end do
    xy_Sum = xy_Sum + ( xy_Rain + xy_Snow ) * 2.0_DP * DelTime

    xy_Ratio = ( xy_Sum - xy_SumB ) / ( xy_Sum + 1.0d-100 )
    do j = 1, jmax
      do i = 0, imax-1
        if ( abs( xy_Ratio(i,j) ) > 1.0d-10 ) then
          call MessageNotify( 'M', module_name, 'H2O mass is not conserved, %f.', d = (/ xy_Ratio(i,j) /) )
        end if
      end do
    end do


  end subroutine CloudSimpleConsChk

  subroutine CloudSimpleEvap1Grid( CharPhase, DelMass, Press, QH2OVap, QH2OVapSat, VirTemp, PRCP, DelPRCPFlux )

    ! 物理・数学定数設定
    ! Physical and mathematical constants settings
    !
    use constants0, only: PI                    ! $ \pi $ .
                              ! 円周率.  Circular constant

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: Grav, GasRDry
                              ! $ R $ [J kg-1 K-1].
                              ! 乾燥大気の気体定数.
                              ! Gas constant of air


    character(*), intent(in ) :: CharPhase
    real(DP)    , intent(in ) :: DelMass
    real(DP)    , intent(in ) :: Press
    real(DP)    , intent(in ) :: QH2OVap
    real(DP)    , intent(in ) :: QH2OVapSat
    real(DP)    , intent(in ) :: VirTemp
    real(DP)    , intent(in ) :: PRCP
    real(DP)    , intent(out) :: DelPRCPFlux


    ! Parameters for evaporation of rain
    real(DP), parameter :: DensWater            = 1.0d3
    !                            rho_w
    !   Values below are from Kessler (1969)
    real(DP), parameter :: PRCPFallVelFactor0        = 130.0d0
    !                            K
    real(DP), parameter :: MedianDiameterFactor      = 3.67d0
    !                            C'
    real(DP), parameter :: PRCPDistFactor            = 1.0d7
    !                            N0
    real(DP), parameter :: PRCPEvapRatUnitDiamFactor = 2.24d3
    !                            C
    real(DP), parameter :: H2OVapDiffCoef            = 1.0d-5
    !                            Kd

    real(DP) :: PRCPFallVelRatio
    real(DP) :: PRCPFallVelFactor

    real(DP) :: Dens0
    !                            rho_0
    real(DP) :: V00
    !                            V_{00}
    real(DP) :: PRCPEvapFactor

    real(DP) :: Dens
    !                           rho
    real(DP) :: DensPRCP
    !                           (rho q_r)
!!$    real(DP) :: RainArea
!!$    !                           a_p
!!$    real(DP) :: RainEvapArea
!!$    !                           A = max( a_p - a, 0 )
!!$    real(DP) :: xy_CloudCover   (0:imax-1, 1:jmax)
!!$    !                           a
    real(DP) :: PRCPEvapRate

    real(DP) :: DelZ


    select case ( CharPhase )
    case ( 'liquid' )
      ! for liquid water
      PRCPFallVelRatio = 1.0_DP
    case ( 'solid' )
      ! for solid water (ice)
      PRCPFallVelRatio = 0.5_DP
    case ( 'mixture' )
      ! for mixture, this is only for test
      PRCPFallVelRatio = 1.0_DP
    case default
      call MessageNotify( 'E', module_name, 'Unexpected character for select case.' )
    end select
    !
    PRCPFallVelFactor = PRCPFallVelFactor0 * PRCPFallVelRatio

    ! Parameters for evaporation of rain
    Dens0 = 1013.0d2 / ( GasRDry * 300.0_DP )
    V00 = PRCPFallVelFactor * sqrt( MedianDiameterFactor ) / ( PI * DensWater * PRCPDistFactor )**(1.0d0/8.0d0)
    PRCPEvapFactor = PRCPEvapRatUnitDiamFactor * 1.429624558860304d0 * H2OVapDiffCoef * PRCPDistFactor**(7.0d0/20.0d0) / ( PI * DensWater )**(13.0d0/20.0d0)
    ! Values for evaporation of rain
    Dens = Press / ( GasRDry * VirTemp )

    DelZ = DelMass / Dens


!!$    RainArea   = RainArea
!!$    xy_CloudCover = CloudCover
!!$    xy_RainEvapArea = max( xy_RainArea - xy_CloudCover, 0.0_DP )
!!$    RainEvapArea = RainEvapArea

    DensPRCP = ( PRCP / ( PRCPArea + 1.0d-10 ) / ( V00 * sqrt( Dens0 / Dens ) ) )**(8.0d0/9.0d0)
    PRCPEvapRate = Dens * PRCPEvapArea * PRCPEvapFactor * max( QH2OVapSat - QH2OVap, 0.0_DP ) * DensPRCP**(13.0d0/20.0d0)

    ! PRCPEvapRate (kg m-3 s-1)
    ! DelZ         (m)
    ! DelPRCPFlux  (kg m-2 s-1)
    DelPRCPFlux = PRCPEvapRate * DelZ

    DelPRCPFlux = min( DelPRCPFlux, PRCP )


  end subroutine CloudSimpleEvap1Grid

  subroutine CloudSimpleWithIceOld( xyr_Press, xyz_Press, xyz_Temp, xyz_QH2OVap, xyz_QH2OLiq, xyz_QH2OSol, xy_Rain, xy_Snow )

    ! USE statements
    !

    ! 時刻管理
    ! Time control
    !
    use timeset, only: DelTime            ! $ \Delta t $ [s]


    real(DP), intent(in   ) :: xyr_Press        ( 0:imax-1, 1:jmax, 0:kmax )
    real(DP), intent(in   ) :: xyz_Press        ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(inout) :: xyz_Temp         ( 0:imax-1, 1:jmax, 1:kmax )
!!$    real(DP), intent(in   ) :: xyz_DQH2OLiqDtCum( 0:imax-1, 1:jmax, 1:kmax )
!!$    real(DP), intent(in   ) :: xyz_DQH2OLiqDtLSC( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(inout) :: xyz_QH2OVap     ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(inout) :: xyz_QH2OLiq      ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(inout) :: xyz_QH2OSol      ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(out  ) :: xy_Rain          ( 0:imax-1, 1:jmax )
    real(DP), intent(out  ) :: xy_Snow          ( 0:imax-1, 1:jmax )


    real(DP) :: xyz_TempB   ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_QH2OVapB( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_QH2OLiqB( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_QH2OSolB( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_DQH2OLiqDt( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_DQH2OSolDt( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_DQRainDt( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_DQSnowDt( 0:imax-1, 1:jmax, 1:kmax )

!!$    real(DP) :: xyz_DTempDtPrcpPCCum( 0:imax-1, 1:jmax, 1:kmax )
!!$    real(DP) :: xyz_DTempDtPrcpPCLsc( 0:imax-1, 1:jmax, 1:kmax )

    ! 実行文 ; Executable statement
    !

    ! 初期化確認
    ! Initialization check
    !
    if ( .not. cloud_simple_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    xyz_DQH2OLiqDt = 0.0_DP
    xyz_DQH2OSolDt = 0.0_DP


    ! Numerical solution

!!$      xyz_DQCloudWaterDt = xyz_DQCloudWaterDt &
!!$        & - xyz_QCloudWater / ( CloudLifeTime + 1.0e-100_DP )


!      ( X_{t+1} - X_{t-1} ) / ( 2 \Delta t ) = Q - X_{t+1} / \tau
!
!      X_{t+1} / ( 2 \Delta t )  + X_{t+1} / \tau = X_{t-1} / ( 2 \Delta t ) + Q
!      ( 1 / ( 2 \Delta t )  + 1 / \tau ) X_{t+1} = X_{t-1} / ( 2 \Delta t ) + Q
!      X_{t+1} = ( X_{t-1} / ( 2 \Delta t ) + Q ) / ( 1 / ( 2 \Delta t )  + 1 / \tau ) 

!!$    xyz_QH2OLiq =                                                           &
!!$      &   (                                                                 &
!!$      &       xyz_QH2OLiq / ( 2.0_DP * DelTime )                            &
!!$      &     + xyz_DQH2OLiqDtCum + xyz_DQH2OLiqDtLSC                         &
!!$      &   )                                                                 &
!!$      & / ( 1.0_DP / ( 2.0_DP * DelTime ) + 1.0_DP / ( CloudLifeTime + 1.0d-100 ) )
!!$
!!$    call CloudUtilsCalcPRCPKeyLLTemp3D(          &
!!$      & xyr_Press, xyz_Temp, xyz_DQH2OLiqDtCum,  &  ! (in )
!!$      & xy_RainCum, xy_SnowCum                   &  ! (out)
!!$      & )
!!$    call CloudUtilsCalcPRCPKeyLLTemp3D(          &
!!$      & xyr_Press, xyz_Temp, xyz_DQH2OLiqDtLsc,  &  ! (in )
!!$      & xy_RainLsc, xy_SnowLsc                   &  ! (out)
!!$      & )
    !-----


    ! Analytical solution

    ! save values before adjustment
    xyz_TempB    = xyz_Temp
    xyz_QH2OVapB = xyz_QH2OVap
    xyz_QH2OLiqB = xyz_QH2OLiq
    xyz_QH2OSolB = xyz_QH2OSol


    xyz_QH2OLiq = xyz_QH2OLiq * exp( - 2.0_DP * DelTime / ( CloudWatLifeTime + 1.0e-100_DP ) ) + xyz_DQH2OLiqDt * CloudWatLifeTime * ( 1.0_DP - exp( - 2.0_DP * DelTime / ( CloudWatLifeTime + 1.0e-100_DP ) ) )

    xyz_DQRainDt = xyz_QH2OLiqB + xyz_DQH2OLiqDt * 2.0_DP * DelTime - xyz_QH2OLiq
    xyz_DQRainDt = xyz_DQRainDt / ( 2.0_DP * DelTime )


    xyz_QH2OSol = xyz_QH2OSol * exp( - 2.0_DP * DelTime / ( CloudIceLifeTime + 1.0e-100_DP ) ) + xyz_DQH2OSolDt * CloudIceLifeTime * ( 1.0_DP - exp( - 2.0_DP * DelTime / ( CloudIceLifeTime + 1.0e-100_DP ) ) )

    xyz_DQSnowDt = xyz_QH2OSolB + xyz_DQH2OSolDt * 2.0_DP * DelTime - xyz_QH2OSol
    xyz_DQSnowDt = xyz_DQSnowDt / ( 2.0_DP * DelTime )



    select case ( IDSnowMethod )
    case ( IDSnowMethodKeyLLTemp )

      xyz_DQRainDt = xyz_DQRainDt + xyz_DQSnowDt
      !
      call CloudSimpleCalcPRCPKeyLLTemp3D( xyr_Press, xyz_Press, xyz_DQRainDt, xyz_Temp, xyz_QH2OVap, xy_Rain, xy_Snow )

      call CloudSimpleConsChk( .false., xyr_Press, xyz_TempB, xyz_QH2OVapB, xyz_QH2OLiqB, xyz_QH2OSolB, xyz_Temp , xyz_QH2OVap , xyz_QH2OLiq , xyz_QH2OSol , xy_Rain, xy_Snow )

    case ( IDSnowMethodStepPC )

      call CloudSimpleCalcPRCPStepPC( xyr_Press, xyz_Press, xyz_DQRainDt, xyz_DQSnowDt, xyz_Temp, xyz_QH2OVap, xy_Rain, xy_Snow )

      call CloudSimpleConsChk( .true., xyr_Press, xyz_TempB, xyz_QH2OVapB, xyz_QH2OLiqB, xyz_QH2OSolB, xyz_Temp , xyz_QH2OVap , xyz_QH2OLiq , xyz_QH2OSol , xy_Rain, xy_Snow )

    end select


  end subroutine CloudSimpleWithIceOld