subroutine CloudT1993base( xyz_Press, xyr_Press, xyz_VirTemp, xyz_DQCloudWaterDtCum, xyz_MoistConvDetTend, xyz_OMG, xyz_MoistConvSubsidMassFlux, xyz_DTempDtPhy, xyz_Temp, xyz_QH2OVap, xyz_QCloudWater, xyz_CloudCover, xy_SurfRainFlux, xy_SurfSnowFlux )

    ! USE statements
    !

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

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: Grav, CpDry, GasRDry, GasRWet, LatentHeat
                              ! $ L $ [J kg-1] .
                              ! 凝結の潜熱.
                              ! Latent heat of condensation

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

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

    ! 簡単雲モデル
    ! Simple cloud model
    !
    use cloud_simple, only : CloudSimpleCalcPRCPKeyLLTemp


    real(DP), intent(in   ) :: xyz_Press                  (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in   ) :: xyr_Press                  (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(in   ) :: xyz_VirTemp                (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in   ) :: xyz_DQCloudWaterDtCum      (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in   ) :: xyz_MoistConvDetTend       (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in   ) :: xyz_OMG                    (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in   ) :: xyz_MoistConvSubsidMassFlux(0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in   ) :: xyz_DTempDtPhy             (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(inout) :: xyz_QCloudWater            (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(inout) :: xyz_CloudCover             (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)


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

    real(DP) :: xyz_QH2OVapSat       (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_DQH2OVapSatDPress(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_DQH2OVapSatDTemp (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_DQH2OVapSatDt    (0:imax-1, 1:jmax, 1:kmax)

    real(DP) :: xyz_ZeroOneCloudProd(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_ZeroOneCloudLoss(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_DelCloudCoverStr(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_FactA           (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_FactA1          (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_FactA2          (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_FactB           (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_FactC           (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_FactC1          (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_FactC2          (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_FactD           (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_FactE           (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_FactE1          (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_FactE2          (0:imax-1, 1:jmax, 1:kmax)

    real(DP) :: xy_Rain                  (0:imax-1, 1:jmax)
    real(DP) :: xyz_Rain                 (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xy_RainConvFactor        (0:imax-1, 1:jmax)
    real(DP) :: xy_FactCo                (0:imax-1, 1:jmax)
    real(DP) :: xy_FactBF                (0:imax-1, 1:jmax)
    real(DP) :: xy_QCloudWatConvThreshold(0:imax-1, 1:jmax)

    real(DP) :: xyz_DTempDtLSC             (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_DQVapDtLSC             (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_RainLSC                (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xy_RainLSC                 (0:imax-1, 1:jmax)


    real(DP), parameter :: MixCoef                  = 1.0d-6
    real(DP), parameter :: QCloudWaterEvapThreshold = 1.0d-10
!!$    real(DP), parameter :: RHThreshold              = 0.999_DP
    real(DP), parameter :: RHThreshold              = 1.0_DP - 1.0d-10
    ! Values below are obtained from Sundqvist et al. (1989), but some of 
    ! those are arbitrarily selected.
    real(DP)            :: RainConvFactor0
    real(DP), parameter :: C1                       = 300.0_DP
    real(DP), parameter :: C2                       = 0.5_DP


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

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

    real(DP) :: xyz_Dens        (0:imax-1, 1:jmax, 1:kmax)
    !                           rho
    real(DP) :: xy_DensRain     (0:imax-1, 1:jmax)
    !                           (rho q_r)
    real(DP) :: xy_RainArea     (0:imax-1, 1:jmax)
    !                           a_p
    real(DP) :: xy_RainEvapArea (0:imax-1, 1:jmax)
    !                           A = max( a_p - a, 0 )
    real(DP) :: xyz_RainEvapRate(0:imax-1, 1:jmax, 1:kmax)

    real(DP) :: xy_DelRain      (0:imax-1, 1:jmax)
    real(DP) :: DelQH2OVap

    real(DP) :: RainFallVel

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

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


    ! 実行文 ; Executable statement
    !

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


    ! Parameters for evaporation of rain
    Dens0 = 1013.0d2 / ( GasRDry * 300.0_DP )
    V00 = RainFallVelFactor * sqrt( MedianDiameterFactor ) / ( PI * DensWater * RainDistFactor )**(1.0d0/8.0d0)
    RainEvapFactor = RainEvapRatUnitDiamFactor * 1.429624558860304d0 * H2OVapDiffCoef * RainDistFactor**(7.0d0/20.0d0) / ( PI * DensWater )**(13.0d0/20.0d0)
    ! Values for evaporation of rain
    xyz_Dens = xyz_Press / ( GasRDry * xyz_VirTemp )


    ! Save cloud water amount
    !
    xyz_QCloudWaterB = xyz_QCloudWater


    xyz_QH2OVapSat       = xyz_CalcQVapSat( xyz_Temp, xyz_Press )
    xyz_DQH2OVapSatDTemp = xyz_CalcDQVapSatDTemp( xyz_Temp, xyz_QH2OVapSat )

    xyz_DQH2OVapSatDPress = xyz_QH2OVapSat / xyz_Press * ( LatentHeat * GasRDry * xyz_VirTemp - CpDry * GasRWet * xyz_Temp**2 ) / ( xyz_CloudCover * LatentHeat**2 * xyz_QH2OVapSat + CpDry * GasRWet * xyz_Temp**2 )

    xyz_DQH2OVapSatDt = xyz_DQH2OVapSatDPress * ( xyz_OMG + Grav * xyz_MoistConvSubsidMassFlux ) + xyz_DQH2OVapSatDTemp / ( 1.0_DP + xyz_CloudCover * LatentHeat / CpDry * xyz_DQH2OVapSatDTemp ) * xyz_DTempDtPhy


    ! set zero-one flag
    do k = 1, kmax
      xyz_Sigma(:,:,k) = xyz_Press(:,:,k) / xyr_Press(:,:,0)
    end do
    xyz_CloudProdRHThreshold = RHThresholdCrtl + ( 1.0_DP - xyz_QH2OVap / xyz_QH2OVapSat ) * (   max( ( xyz_Sigma - RHThresholdSigmaMin ), 0.0_DP ) / ( xyz_Sigma - RHThresholdSigmaMin )                )**RHThresholdOrd
    do k = 1, kmax
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_DQH2OVapSatDt(i,j,k) < 0.0_DP ) then
            if ( xyz_QH2OVap(i,j,k) >= xyz_CloudProdRHThreshold(i,j,k) * xyz_QH2OVapSat(i,j,k) ) then
              xyz_ZeroOneCloudProd(i,j,k) = 1.0_DP
            else
              xyz_ZeroOneCloudProd(i,j,k) = 0.0_DP
            end if
            xyz_ZeroOneCloudLoss(i,j,k) = 1.0_DP
          else
            xyz_ZeroOneCloudProd(i,j,k) = 0.0_DP
            xyz_ZeroOneCloudLoss(i,j,k) = 0.0_DP
          end if
        end do
      end do
    end do


    ! Rain at the surface
    xy_Rain     = 0.0_DP
    ! Rain area
    xy_RainArea = 0.0_DP

    k_loop : do k = kmax, 1, -1

      ! Evaporation of rain
      !
      if ( k == kmax ) then
        xy_RainArea = 0.0_DP
      else
        do j = 1, jmax
          do i = 0, imax-1
            if ( xyz_Rain(i,j,k+1) > 0.0_DP ) then
              if ( xyz_CloudCover(i,j,k+1) > xy_RainArea(i,j) ) then
                xy_RainArea(i,j) = xyz_CloudCover(i,j,k+1)
              end if
            end if
          end do
        end do
      end if
      xy_DensRain = ( xy_Rain / ( xy_RainArea + 1.0d-10 ) / ( V00 * sqrt( Dens0 / xyz_Dens(:,:,k) ) ) )**(8.0d0/9.0d0)
      xy_RainEvapArea = max( xy_RainArea - xyz_CloudCover(:,:,k), 0.0_DP )
      xyz_RainEvapRate(:,:,k) = xyz_Dens(:,:,k) * xy_RainEvapArea * RainEvapFactor * max( xyz_QH2OVapSat(:,:,k) - xyz_QH2OVap(:,:,k), 0.0_DP ) * xy_DensRain**(13.0d0/20.0d0)
      do j = 1, jmax
        do i = 0, imax-1
          RainFallVel = V00 * sqrt( Dens0 / xyz_Dens(i,j,k) ) * xy_DensRain(i,j)**(1.0d0/8.0d0)
          if ( xy_RainArea(i,j) * RainFallVel * xyz_RainEvapRate(i,j,k) * 2.0_DP * DelTime > xy_Rain(i,j) ) then
            xyz_RainEvapRate(i,j,k) = xy_Rain(i,j) / ( ( xy_RainArea(i,j) + 1.0d-10 ) * RainFallVel * 2.0_DP * DelTime )
            xy_DelRain(i,j) = - xy_Rain(i,j)
          else
            xy_DelRain(i,j) = - xy_RainArea(i,j) * RainFallVel * xyz_RainEvapRate(i,j,k) * ( 2.0_DP * DelTime )
          end if
          xy_Rain(i,j) = xy_Rain(i,j) + xy_DelRain(i,j)
!!$            xyz_QH2OVap(i,j,k) = xyz_QH2OVap(i,j,k)              &
!!$              & + xyz_RainEvapRate(i,j,k) * ( 2.0_DP * DelTime ) &
!!$              &     / xyz_Dens(i,j,k)
!!$            xyz_Temp(i,j,k) = xyz_Temp(i,j,k)                    &
!!$              & - xyz_RainEvapRate(i,j,k) * ( 2.0_DP * DelTime ) &
!!$              &     / xyz_Dens(i,j,k)                            &
!!$              &     * LatentHeat / CpDry
          ! DelRain = DelQH2OVap * DelPress / Grav / ( 2.0_DP * DelTime )
          DelQH2OVap = - xy_DelRain(i,j) * Grav / ( xyr_Press(i,j,k-1) - xyr_Press(i,j,k) )
          xyz_QH2OVap(i,j,k) = xyz_QH2OVap(i,j,k) + DelQH2OVap
          xyz_Temp(i,j,k) = xyz_Temp(i,j,k) - DelQH2OVap * LatentHeat / CpDry
        end do
      end do


      xyz_FactC1(:,:,k) = xyz_MoistConvDetTend(:,:,k)
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_QH2OVap(i,j,k) < RHThreshold * xyz_QH2OVapSat(i,j,k) ) then
            xyz_FactC2(i,j,k) = - ( 1.0_DP - xyz_CloudCover(i,j,k) ) / ( xyz_QH2OVapSat(i,j,k) - xyz_QH2OVap(i,j,k) ) * xyz_DQH2OVapSatDt(i,j,k) * xyz_ZeroOneCloudProd(i,j,k)
          else
!!$              xyz_FactC2(i,j,k) = 0.0_DP
            xyz_FactC2(i,j,k) = 1.0_DP / ( 2.0_DP * DelTime )
          end if
        end do
      end do
      xyz_FactC(:,:,k) = xyz_FactC1(:,:,k) + xyz_FactC2(:,:,k)
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_QH2OVapSat(i,j,k) > xyz_QH2OVap(i,j,k) ) then
            xyz_FactD(i,j,k) = 2.0_DP * xyz_CloudCover(i,j,k) * MixCoef * ( xyz_QH2OVapSat(i,j,k) - xyz_QH2OVap(i,j,k) ) / ( xyz_QCloudWater(i,j,k) + 1.0d-100 )
          else
            xyz_FactD(i,j,k) = 0.0_DP
          end if
        end do
      end do
      !
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_QH2OVap(i,j,k) < RHThreshold * xyz_QH2OVapSat(i,j,k) ) then
            xyz_FactE1(i,j,k) = ( 1.0_DP - xyz_CloudCover(i,j,k) )**2 / ( 2.0_DP * ( xyz_QH2OVapSat(i,j,k) - xyz_QH2OVap(i,j,k) ) ) * xyz_DQH2OVapSatDt(i,j,k) * xyz_ZeroOneCloudProd(i,j,k)
          else
            xyz_FactE1(i,j,k) = 0.0_DP
          end if
        end do
      end do
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_QH2OVapSat(i,j,k) > xyz_QH2OVap(i,j,k) ) then
            xyz_FactE2(i,j,k) = + xyz_CloudCover(i,j,k)**2 * MixCoef * ( xyz_QH2OVapSat(i,j,k) - xyz_QH2OVap(i,j,k) ) / ( xyz_QCloudWater(i,j,k) + 1.0d-100 )
          else
            xyz_FactE2(i,j,k) = 0.0_DP
          end if
        end do
      end do
      xyz_FactE(:,:,k) = xyz_FactE1(:,:,k) + xyz_FactE2(:,:,k)
      !
      xyz_DelCloudCoverStr(:,:,k) = xyz_FactC2(:,:,k) * 2.0_DP * DelTime - xyz_FactC2(:,:,k) / ( xyz_FactC(:,:,k) + xyz_FactD(:,:,k) + 1.0d-100 ) * ( xyz_FactC(:,:,k) * 2.0_DP * DelTime + ( xyz_CloudCover(:,:,k) - xyz_FactC(:,:,k) / ( xyz_FactC(:,:,k) + xyz_FactD(:,:,k) + 1.0d-100 ) ) * ( 1.0_DP - exp( - ( xyz_FactC(:,:,k) + xyz_FactD(:,:,k) ) * 2.0_DP * DelTime ) ) ) + xyz_FactE1(:,:,k) * 2.0_DP * DelTime
      !
      xyz_FactA1(:,:,k) = xyz_DQCloudWaterDtCum(:,:,k)
      xyz_FactA2(:,:,k) = - xyz_CloudCover(:,:,k) * xyz_DQH2OVapSatDt(:,:,k) * xyz_ZeroOneCloudProd(:,:,k) - xyz_DelCloudCoverStr(:,:,k) * xyz_DQH2OVapSatDt(:,:,k) * xyz_ZeroOneCloudProd(:,:,k) / 2.0_DP - xyz_CloudCover(:,:,k) * xyz_DQH2OVapSatDt(:,:,k) * ( 1.0_DP - xyz_ZeroOneCloudLoss(:,:,k) ) - xyz_CloudCover(:,:,k) * MixCoef * ( xyz_QH2OVapSat(:,:,k) - xyz_QH2OVap(:,:,k) )
      !    The value of xyz_FactA2 is checked, and is updated.
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_FactA2(i,j,k) * 2.0_DP * DelTime > xyz_QH2OVap(i,j,k) ) then
            xyz_FactA2(i,j,k) = xyz_QH2OVap(i,j,k) / ( 2.0_DP * DelTime ) * ( 1.0_DP - 1.0d-14 )
          end if
        end do
      end do
      xyz_FactA(:,:,k) = xyz_FactA1(:,:,k) + xyz_FactA2(:,:,k)

!!$        xy_RainConvFactor = 1.0_DP / ( CloudLifeTime0 + 1.0d-100 )
      !
      xy_FactCo = 1.0_DP + C1 * sqrt( xy_Rain )
      ! Factor for Bergeron-Findeisen effect
      !   Sundqvist et al. (1989)
!!$      xy_FactBF = 1.0_DP &
!!$        & + C2 * sqrt( max( 268.0_DP - xyz_Temp(:,:,k), 0.0_DP ) )
      !   An original equation following that by IFS CY38r1
      !   (p.82 of http://www.ecmwf.int/research/ifsdocs/CY38r1/IFSPart4.pdf)
!!$      xy_FactBF = 1.0_DP                                               &
!!$        & + C2 * sqrt( min(                                            &
!!$        &                   max( 268.0_DP - xyz_Temp(:,:,k), 0.0_DP ), &
!!$        &                   max( 268.0_DP - TempBFEffectSat, 0.0_DP )  &
!!$        &                  )                                           &
!!$        &            )
      !   Constant (no effect)
      xy_FactBF = 1.0_DP
      !
      RainConvFactor0 = 1.0_DP / ( CloudLifeTime0 + 1.0d-100 )
!!$      xy_QCloudWatConvThreshold = &
!!$        & QCloudWatEffConv0 / ( xy_FactCo * xy_FactBF )
      xy_QCloudWatConvThreshold = ( QCloudWatEffConv0 + 1.0e-100_DP ) / ( xy_FactCo * xy_FactBF )
      xy_RainConvFactor = RainConvFactor0 * xy_FactCo * xy_FactBF * ( 1.0_DP - exp( - ( xyz_QCloudWater(:,:,k) / ( ( xyz_CloudCover(:,:,k) + 1.0d-100 ) * xy_QCloudWatConvThreshold )   )**2 ) )
      !
      xyz_FactB(:,:,k) = xy_RainConvFactor
      !
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_FactB(i,j,k) < 1.0_DP / 1.0d10 ) then
            xyz_FactB(i,j,k) = 1.0_DP / 1.0d10
          end if
        end do
      end do


      ! Values at next time step are calculated.
      !
      xyz_QCloudWater(:,:,k) = xyz_QCloudWater(:,:,k) * exp( - xyz_FactB(:,:,k) * 2.0_DP * DelTime ) + xyz_FactA(:,:,k) / xyz_FactB(:,:,k) * ( 1.0_DP - exp( - xyz_FactB(:,:,k) * 2.0_DP * DelTime ) )
      !   The value of cloud water amount is checked, and xyz_FactA 
      !   is updated.
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_QCloudWater(i,j,k) < 0.0_DP ) then
            xyz_FactA2(i,j,k) = - xyz_FactA1(i,j,k) - xyz_FactB(i,j,k) * xyz_QCloudWaterB(i,j,k) * exp( - xyz_FactB(i,j,k) * 2.0_DP * DelTime ) / ( 1.0_DP - exp( - xyz_FactB(i,j,k) * 2.0_DP * DelTime ) )
            xyz_QCloudWater(i,j,k) = 0.0_DP
          end if
        end do
      end do
      xyz_FactA(:,:,k) = xyz_FactA1(:,:,k) + xyz_FactA2(:,:,k)
      !
      xyz_CloudCover(:,:,k) = xyz_CloudCover(:,:,k) * exp( - ( xyz_FactC(:,:,k) + xyz_FactD(:,:,k) ) * 2.0_DP * DelTime ) + ( xyz_FactC(:,:,k) + xyz_FactE(:,:,k) ) / ( xyz_FactC(:,:,k) + xyz_FactD(:,:,k) + 1.0d-100 ) * ( 1.0_DP - exp( - ( xyz_FactC(:,:,k) + xyz_FactD(:,:,k) ) * 2.0_DP * DelTime ) )
      !
      xyz_QH2OVap(:,:,k) = xyz_QH2OVap(:,:,k) - xyz_FactA2(:,:,k) * 2.0_DP * DelTime
      !
      xyz_Temp(:,:,k) = xyz_Temp(:,:,k) + xyz_FactA2(:,:,k) * 2.0_DP * DelTime * LatentHeat / CpDry


      ! Rain
      !
      xyz_Rain(:,:,k) = xyz_FactA(:,:,k) * 2.0_DP * DelTime + xyz_QCloudWaterB(:,:,k) * ( 1.0_DP - exp( - xyz_FactB(:,:,k) * 2.0_DP * DelTime ) ) - xyz_FactA(:,:,k) / xyz_FactB(:,:,k) * ( 1.0_DP - exp( - xyz_FactB(:,:,k) * 2.0_DP * DelTime ) )
      xyz_Rain(:,:,k) = xyz_Rain(:,:,k) / ( 2.0_DP * DelTime )

      ! Rain at the surface
      xy_Rain = xy_Rain + xyz_Rain(:,:,k) * ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav

      ! Evaporation
      !
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_QCloudWater(i,j,k) < QCloudWaterEvapThreshold ) then
            xyz_CloudCover(i,j,k) = 0.0_DP
          end if
        end do
      end do

      ! Cloud cover is restricted.
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_CloudCover(i,j,k) > 1.0_DP ) then
            xyz_CloudCover(i,j,k) = 1.0_DP
          else if ( xyz_CloudCover(i,j,k) < 0.0_DP ) then
            xyz_CloudCover(i,j,k) = 0.0_DP
          end if
        end do
      end do


      ! Check values
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_QH2OVap(i,j,k) < 0.0_DP ) then
            write( 6, * ) 'QH2OVap is negative', i, j, k, xyz_QH2OVap(i,j,k)
          end if
          if ( xyz_QCloudWater(i,j,k) < 0.0_DP ) then
            write( 6, * ) 'QCloudWater is negative', i, j, k, xyz_QCloudWater(i,j,k)
          end if
        end do
      end do

    end do k_loop


    xyz_DTempDtLSC = 0.0_DP
    xyz_DQVapDtLSC = 0.0_DP

    ! 大規模凝結 (非対流性凝結) (Manabe, 1965)
    ! Large scale condensation (non-convective condensation) (Manabe, 1965)
    !
    call LScaleCond( xyz_Temp, xyz_QH2OVap, xyz_Press, xyr_Press, xyz_RainLSC )

    xy_RainLSC = 0.0_DP
    do k = kmax, 1, -1
      xy_RainLSC = xy_RainLSC + xyz_RainLSC(:,:,k) * ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav
    end do


    xy_Rain = xy_Rain + xy_RainLSC


    call CloudSimpleCalcPRCPKeyLLTemp( xyz_Temp, xy_Rain, xy_SurfRainFlux, xy_SurfSnowFlux )


  end subroutine CloudT1993base

  subroutine CloudT1993baseInit( 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


    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_T1993base_nml/ RHThresholdCrtl, RHThresholdSigmaMin, RHThresholdOrd, CloudLifeTime0, CloudWatLifeTime0, CloudIceLifeTime0, QCloudWatEffConv0, QCloudIceEffConv0, TempBFEffectSat, PRCPArea, PRCPEvapArea, PRCPColFactor
          !
          ! デフォルト値については初期化手続 "cloud_T1993base#CloudT1993baseInit"
          ! のソースコードを参照のこと.
          !
          ! Refer to source codes in the initialization procedure
          ! "cloud_T1993base#CloudT1993baseInit" for the default values.
          !

    ! 実行文 ; Executable statement
    !

    if ( cloud_T1993base_inited ) return


    FlagSnow = ArgFlagSnow


    ! デフォルト値の設定
    ! Default values settings
    !
    RHThresholdCrtl      = 0.5_DP
!!$    RHThresholdCrtl      = 0.8_DP ! ECMWF IFS
    RHThresholdSigmaMin  = 1.0_DP
!!$    RHThresholdSigmaMin  = 0.8_DP ! ECMWF IFS
    RHThresholdOrd       = 2.0_DP ! ECMWF IFS
    CloudLifeTime0       = 1000.0_DP
    CloudWatLifeTime0    = 1000.0_DP
    CloudIceLifeTime0    = 1000.0_DP
    QCloudWatEffConv0    = 1.0e-4_DP
    QCloudIceEffConv0    = 1.0e-3_DP
    TempBFEffectSat      = 250.0_DP
      !   This value follows that by IFS CY38r1
      !   (p.82 of http://www.ecmwf.int/research/ifsdocs/CY38r1/IFSPart4.pdf)
      !   Actually, IFS CY38r1 uses 250.16 K.
    PRCPArea            = 1.0_DP
!!$    PRCPArea            = 0.5_DP
    PRCPEvapArea        = 1.0_DP

    PRCPColFactor       = 300.0_DP
                          ! This value comes from Sundqvist et al. (1989)
                          ! IFS CY38r1 uses 100, but in addition , 
                          ! precipitation area has to be considered.

    ! 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_T1993base_nml, iostat = iostat_nml )             ! (out)
      close( unit_nml )

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


    ! 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, 'RHThresholdCrtl      = %f', d = (/ RHThresholdCrtl /) )
    call MessageNotify( 'M', module_name, 'RHThresholdSigmaMin  = %f', d = (/ RHThresholdSigmaMin /) )
    call MessageNotify( 'M', module_name, 'RHThresholdOrd       = %f', d = (/ RHThresholdOrd /) )
    call MessageNotify( 'M', module_name, 'CloudLifeTime0       = %f', d = (/ CloudLifeTime0 /) )
    call MessageNotify( 'M', module_name, 'CloudIceLifeTime0    = %f', d = (/ CloudIceLifeTime0 /) )
    call MessageNotify( 'M', module_name, 'CloudWatLifeTime0    = %f', d = (/ CloudWatLifeTime0 /) )
    call MessageNotify( 'M', module_name, 'QCloudWatEffConv0    = %f', d = (/ QCloudWatEffConv0 /) )
    call MessageNotify( 'M', module_name, 'QCloudIceEffConv0    = %f', d = (/ QCloudIceEffConv0 /) )
    call MessageNotify( 'M', module_name, 'TempBFEffectSat      = %f', d = (/ TempBFEffectSat /) )
    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_T1993base_inited = .true.

  end subroutine CloudT1993baseInit

  subroutine CloudT1993baseWithIce( xyz_Press, xyr_Press, xyz_VirTemp, xyz_DQCloudWatDtCum, xyz_DQCloudIceDtCum, xyz_MoistConvDetTend, xyz_OMG, xyz_MoistConvSubsidMassFlux, xyz_DTempDtPhy, xyz_Temp, xyz_QH2OVap, xyz_QCloudWat, xyz_QCloudIce, xyz_CloudCover, xy_SurfRainFlux, xy_SurfSnowFlux )

    ! USE statements
    !

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

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

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

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

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

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

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


    real(DP), intent(in   ) :: xyz_Press                  (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in   ) :: xyr_Press                  (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(in   ) :: xyz_VirTemp                (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in   ) :: xyz_DQCloudWatDtCum        (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in   ) :: xyz_DQCloudIceDtCum        (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in   ) :: xyz_MoistConvDetTend       (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in   ) :: xyz_OMG                    (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in   ) :: xyz_MoistConvSubsidMassFlux(0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in   ) :: xyz_DTempDtPhy             (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(inout) :: xyz_QCloudWat              (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(inout) :: xyz_QCloudIce              (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(inout) :: xyz_CloudCover             (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)


    ! Local variables
    !
    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_QCloudWatB       (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_QCloudIceB       (0:imax-1, 1:jmax, 1:kmax)

    real(DP) :: LatentHeatSubl

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

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

!!$    real(DP) :: xyz_QH2OVapSatOnLiq  (0:imax-1, 1:jmax, 1:kmax)
!!$    real(DP) :: xyz_QH2OVapSatOnIce  (0:imax-1, 1:jmax, 1:kmax)

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

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

    real(DP) :: xyz_ZOCloudProd   (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_ZOCloudLoss   (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_DelCloudCoverStr(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_FactAl          (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_FactAs          (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_FactAl1         (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_FactAs1         (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_FactA20         (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_FactAl2         (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_FactAs2         (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_FactBl          (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_FactBs          (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_FactC           (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_FactC1          (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_FactC2          (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_FactD           (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_FactE           (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_FactE1          (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_FactE2          (0:imax-1, 1:jmax, 1:kmax)

    real(DP) :: xy_Rain                    (0:imax-1, 1:jmax)
    real(DP) :: xyz_Rain                   (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xy_Snow                    (0:imax-1, 1:jmax)
    real(DP) :: xyz_Snow                   (0:imax-1, 1:jmax, 1:kmax)

    real(DP) :: xy_FactCo                (0:imax-1, 1:jmax)
    real(DP) :: xy_FactBF                (0:imax-1, 1:jmax)
    real(DP) :: xy_QCloudWatConvThreshold(0:imax-1, 1:jmax)
    real(DP) :: xy_QCloudIceConvThreshold(0:imax-1, 1:jmax)
    real(DP) :: xy_FactIceTempDep        (0:imax-1, 1:jmax)
    real(DP) :: xy_FactConvThreshold     (0:imax-1, 1:jmax)

    real(DP) :: xy_RainConvFactor          (0:imax-1, 1:jmax)
    real(DP) :: xy_SnowConvFactor          (0:imax-1, 1:jmax)

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


    real(DP), parameter :: MixCoef                  = 1.0d-6
    real(DP), parameter :: QCloudWaterEvapThreshold = 1.0d-10
!!$    real(DP), parameter :: RHThreshold              = 0.999_DP
    real(DP), parameter :: RHThreshold              = 1.0_DP - 1.0d-10
    ! Values below are obtained from Sundqvist et al. (1989), but some of 
    ! those are arbitrarily selected.
    real(DP)            :: RainConvFactor0
    real(DP)            :: SnowConvFactor0
    real(DP), parameter :: C2                       = 0.5_DP


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

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

    real(DP) :: xyz_Dens        (0:imax-1, 1:jmax, 1:kmax)
    !                           rho
    real(DP) :: xy_DensRain     (0:imax-1, 1:jmax)
    !                           (rho q_r)
    real(DP) :: xy_RainArea     (0:imax-1, 1:jmax)
    !                           a_p
    real(DP) :: xy_RainEvapArea (0:imax-1, 1:jmax)
    !                           A = max( a_p - a, 0 )
    real(DP) :: xyz_RainEvapRate(0:imax-1, 1:jmax, 1:kmax)

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

    real(DP) :: QCloudWatTentative
    real(DP) :: QCloudIceTentative

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

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


    ! 実行文 ; Executable statement
    !

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


    ! Save current values
    !
    xyz_TempB      = xyz_Temp
    xyz_QH2OVapB   = xyz_QH2OVap
    xyz_QCloudWatB = xyz_QCloudWat
    xyz_QCloudIceB = xyz_QCloudIce


    ! Latent heat for sublimation (sum of evaporation and fusion)
    LatentHeatSubl = LatentHeat + LatentHeatFusion


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

    ! Parameters for evaporation of rain
    Dens0 = 1013.0d2 / ( GasRDry * 300.0_DP )
    V00 = RainFallVelFactor * sqrt( MedianDiameterFactor ) / ( PI * DensWater * RainDistFactor )**(1.0d0/8.0d0)
    RainEvapFactor = RainEvapRatUnitDiamFactor * 1.429624558860304d0 * H2OVapDiffCoef * RainDistFactor**(7.0d0/20.0d0) / ( PI * DensWater )**(13.0d0/20.0d0)
    ! Values for evaporation of rain
    xyz_Dens = xyz_Press / ( GasRDry * xyz_VirTemp )


    ! Calculate liquid water fraction and ice fraction
    !
    call SaturateWatFraction( xyz_Temp, xyz_WatFrac )
    xyz_IceFrac = 1.0_DP - xyz_WatFrac


    xyz_QH2OVapSat       = xyz_CalcQVapSat( xyz_Temp, xyz_Press )
    xyz_DQH2OVapSatDTemp = xyz_CalcDQVapSatDTemp( xyz_Temp, xyz_QH2OVapSat )

!!$    xyz_QH2OVapSatOnLiq       = xyz_CalcQVapSatOnLiq( xyz_Temp, xyz_Press )
!!$    xyz_DQH2OVapSatOnLiqDTemp = xyz_CalcDQVapSatOnLiqDTemp( xyz_Temp, xyz_QH2OVapSatOnLiq )
!!$    xyz_QH2OVapSatOnIce       = xyz_CalcQVapSatOnIce( xyz_Temp, xyz_Press )
!!$    xyz_DQH2OVapSatOnIceDTemp = xyz_CalcDQVapSatOnIceDTemp( xyz_Temp, xyz_QH2OVapSatOnIce )

    xyz_DQH2OVapSatDPressDenom = 1.0_DP + xyz_CloudCover / CpDry * ( xyz_WatFrac * LatentHeat + xyz_IceFrac * LatentHeatSubl ) * xyz_DQH2OVapSatDTemp

    xyz_DQH2OVapSatDPress = ( - xyz_QH2OVapSat + GasRDry * xyz_VirTemp / CpDry * xyz_DQH2OVapSatDTemp ) / ( xyz_Press * xyz_DQH2OVapSatDPressDenom )

    xyz_DQH2OVapSatDt = xyz_DQH2OVapSatDPress * ( xyz_OMG + Grav * xyz_MoistConvSubsidMassFlux ) + xyz_DQH2OVapSatDTemp / xyz_DQH2OVapSatDPressDenom * xyz_DTempDtPhy


    ! set zero-one flag
    do k = 1, kmax
      xyz_Sigma(:,:,k) = xyz_Press(:,:,k) / xyr_Press(:,:,0)
    end do
    xyz_CloudProdRHThreshold = RHThresholdCrtl + ( 1.0_DP - xyz_QH2OVap / xyz_QH2OVapSat ) * (   max( ( xyz_Sigma - RHThresholdSigmaMin ), 0.0_DP ) / ( xyz_Sigma - RHThresholdSigmaMin )                )**RHThresholdOrd
    do k = 1, kmax
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_DQH2OVapSatDt(i,j,k) < 0.0_DP ) then
            if ( xyz_QH2OVap(i,j,k) >= xyz_CloudProdRHThreshold(i,j,k) * xyz_QH2OVapSat(i,j,k) ) then
              xyz_ZOCloudProd(i,j,k) = 1.0_DP
            else
              xyz_ZOCloudProd(i,j,k) = 0.0_DP
            end if
            xyz_ZOCloudLoss(i,j,k) = 1.0_DP
          else
            xyz_ZOCloudProd(i,j,k) = 0.0_DP
            xyz_ZOCloudLoss(i,j,k) = 0.0_DP
          end if
        end do
      end do
    end do


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

    ! Rain area
    xy_RainArea = 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_FactC1(:,:,k) = xyz_MoistConvDetTend(:,:,k)
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_QH2OVap(i,j,k) < RHThreshold * xyz_QH2OVapSat(i,j,k) ) then
            xyz_FactC2(i,j,k) = - ( 1.0_DP - xyz_CloudCover(i,j,k) ) / ( xyz_QH2OVapSat(i,j,k) - xyz_QH2OVap(i,j,k) ) * xyz_DQH2OVapSatDt(i,j,k) * xyz_ZOCloudProd(i,j,k)
          else
!!$              xyz_FactC2(i,j,k) = 0.0_DP
            xyz_FactC2(i,j,k) = 1.0_DP / ( 2.0_DP * DelTime )
          end if
        end do
      end do
      xyz_FactC(:,:,k) = xyz_FactC1(:,:,k) + xyz_FactC2(:,:,k)
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_QH2OVapSat(i,j,k) > xyz_QH2OVap(i,j,k) ) then
            xyz_FactD(i,j,k) = 2.0_DP * xyz_CloudCover(i,j,k) * MixCoef * ( xyz_QH2OVapSat(i,j,k) - xyz_QH2OVap(i,j,k) ) / ( xyz_QCloudWat(i,j,k) + 1.0d-100 )
          else
            xyz_FactD(i,j,k) = 0.0_DP
          end if
        end do
      end do
      !
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_QH2OVap(i,j,k) < RHThreshold * xyz_QH2OVapSat(i,j,k) ) then
            xyz_FactE1(i,j,k) = ( 1.0_DP - xyz_CloudCover(i,j,k) )**2 / ( 2.0_DP * ( xyz_QH2OVapSat(i,j,k) - xyz_QH2OVap(i,j,k) ) ) * xyz_DQH2OVapSatDt(i,j,k) * xyz_ZOCloudProd(i,j,k)
          else
            xyz_FactE1(i,j,k) = 0.0_DP
          end if
        end do
      end do
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_QH2OVapSat(i,j,k) > xyz_QH2OVap(i,j,k) ) then
            xyz_FactE2(i,j,k) = + xyz_CloudCover(i,j,k)**2 * MixCoef * ( xyz_QH2OVapSat(i,j,k) - xyz_QH2OVap(i,j,k) ) / ( xyz_QCloudWat(i,j,k) + 1.0d-100 )
          else
            xyz_FactE2(i,j,k) = 0.0_DP
          end if
        end do
      end do
      xyz_FactE(:,:,k) = xyz_FactE1(:,:,k) + xyz_FactE2(:,:,k)
      !
      xyz_DelCloudCoverStr(:,:,k) = xyz_FactC2(:,:,k) * 2.0_DP * DelTime - xyz_FactC2(:,:,k) / ( xyz_FactC(:,:,k) + xyz_FactD(:,:,k) + 1.0d-100 ) * ( xyz_FactC(:,:,k) * 2.0_DP * DelTime + ( xyz_CloudCover(:,:,k) - xyz_FactC(:,:,k) / ( xyz_FactC(:,:,k) + xyz_FactD(:,:,k) + 1.0d-100 ) ) * ( 1.0_DP - exp( - ( xyz_FactC(:,:,k) + xyz_FactD(:,:,k) ) * 2.0_DP * DelTime ) ) ) + xyz_FactE1(:,:,k) * 2.0_DP * DelTime
      !
      !
      !
      !
      xyz_FactAl1(:,:,k) = xyz_DQCloudWatDtCum(:,:,k)
      xyz_FactAs1(:,:,k) = xyz_DQCloudIceDtCum(:,:,k)
      !
      !
      !
      xyz_FactA20(:,:,k) = - xyz_CloudCover(:,:,k) * xyz_DQH2OVapSatDt(:,:,k) * xyz_ZOCloudProd(:,:,k) - xyz_DelCloudCoverStr(:,:,k) * xyz_DQH2OVapSatDt(:,:,k) * xyz_ZOCloudProd(:,:,k) / 2.0_DP - xyz_CloudCover(:,:,k) * xyz_DQH2OVapSatDt(:,:,k) * ( 1.0_DP - xyz_ZOCloudLoss(:,:,k) ) - xyz_CloudCover(:,:,k) * MixCoef * ( xyz_QH2OVapSat(:,:,k) - xyz_QH2OVap(:,:,k) )
      !
      !    The value of xyz_FactA2 is checked, and is updated.
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_FactA20(i,j,k) * 2.0_DP * DelTime > xyz_QH2OVap(i,j,k) ) then
            xyz_FactA20(i,j,k) = xyz_QH2OVap(i,j,k) / ( 2.0_DP * DelTime ) * ( 1.0_DP - 1.0e-14_DP )
          end if
        end do
      end do
      xyz_FactAl2(:,:,k) = xyz_WatFrac(:,:,k) * xyz_FactA20(:,:,k)
      xyz_FactAs2(:,:,k) = xyz_IceFrac(:,:,k) * xyz_FactA20(:,:,k)


      xyz_FactAl(:,:,k) = xyz_FactAl1(:,:,k) + xyz_FactAl2(:,:,k)
      xyz_FactAs(:,:,k) = xyz_FactAs1(:,:,k) + xyz_FactAs2(:,:,k)


      RainConvFactor0 = 1.0_DP / ( CloudWatLifeTime0 + 1.0d-100 )
      ! Factor for collection
      xy_FactCo = 1.0_DP + PRCPColFactor * sqrt( xy_Rain + xy_Snow )
      ! Factor for Bergeron-Findeisen effect
      !   Sundqvist et al. (1989)
!!$      xy_FactBF = 1.0_DP &
!!$        & + C2 * sqrt( max( 268.0_DP - xyz_Temp(:,:,k), 0.0_DP ) )
      !   An original equation following that by IFS CY38r1
      !   (p.82 of http://www.ecmwf.int/research/ifsdocs/CY38r1/IFSPart4.pdf)
!!$      xy_FactBF = 1.0_DP                                               &
!!$        & + C2 * sqrt( min(                                            &
!!$        &                   max( 268.0_DP - xyz_Temp(:,:,k), 0.0_DP ), &
!!$        &                   max( 268.0_DP - TempBFEffectSat, 0.0_DP )  &
!!$        &                  )                                           &
!!$        &            )
      !   Constant (no effect)
      xy_FactBF = 1.0_DP
      !
!!$      xy_QCloudWatConvThreshold = &
!!$        & QCloudWatEffConv0 / ( xy_FactCo * xy_FactBF )
      xy_QCloudWatConvThreshold = ( QCloudWatEffConv0 + 1.0e-100_DP ) / ( xy_FactCo * xy_FactBF )
!!$      xy_QCloudWatConvThreshold = &
!!$        & QCloudWatEffConv0 / xy_FactCo
      xy_FactConvThreshold = 1.0_DP - exp( - ( xyz_QCloudWat(:,:,k) / ( xyz_CloudCover(:,:,k) + 1.0e-10_DP ) / xy_QCloudWatConvThreshold              )**2 )
      !
!!$      xy_RainConvFactor = RainConvFactor0 * xy_FactCo * xy_FactBF
      xy_RainConvFactor = RainConvFactor0 * xy_FactCo * xy_FactBF * xy_FactConvThreshold


      SnowConvFactor0 = 1.0_DP / ( CloudIceLifeTime0 + 1.0d-100 )
      xy_FactIceTempDep = exp( 0.025_DP * ( xyz_Temp(:,:,k) - 273.15_DP ) )
!!$      xy_QCloudIceConvThreshold = QCloudIceEffConv0
      xy_QCloudIceConvThreshold = QCloudIceEffConv0 + 1.0e-100_DP
      xy_FactConvThreshold = 1.0_DP - exp( - ( xyz_QCloudIce(:,:,k) / ( xyz_CloudCover(:,:,k) + 1.0e-10_DP ) / xy_QCloudIceConvThreshold              )**2 )
!!$      xy_SnowConvFactor = SnowConvFactor0 * xy_FactIceTempDep
      xy_SnowConvFactor = SnowConvFactor0 * xy_FactIceTempDep * xy_FactConvThreshold

      !
      !
      !
!!$      xyz_FactBl(:,:,k) = 1.0_DP / ( CloudWatLifeTime0 + 1.0d-100 )
      xyz_FactBl(:,:,k) = xy_RainConvFactor
      !
!!$      xyz_FactBs(:,:,k) = 1.0_DP / ( CloudIceLifeTime0 + 1.0d-100 )
      xyz_FactBs(:,:,k) = xy_SnowConvFactor
      !
      !
      !


      ! Values at next time step are calculated.
      !
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_FactBl(i,j,k) >= FactBlsThreshold ) then
            xyz_QCloudWat(i,j,k) = xyz_QCloudWat(i,j,k) * exp( - xyz_FactBl(i,j,k) * 2.0_DP * DelTime ) + xyz_FactAl(i,j,k) / xyz_FactBl(i,j,k) * ( 1.0_DP - exp( - xyz_FactBl(i,j,k) * 2.0_DP * DelTime ) )
          else
            xyz_QCloudWat(i,j,k) = xyz_QCloudWat(i,j,k) + xyz_FactAl(i,j,k) * ( 2.0_DP * DelTime )
          end if
        end do
      end do
      !   The value of cloud water amount is checked, and xyz_FactA is updated.
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_QCloudWat(i,j,k) < 0.0_DP ) then
            if ( xyz_FactBl(i,j,k) >= FactBlsThreshold ) then
              xyz_FactAl2(i,j,k) = - xyz_FactAl1(i,j,k) - xyz_FactBl(i,j,k) * xyz_QCloudWatB(i,j,k) * exp( - xyz_FactBl(i,j,k) * 2.0_DP * DelTime ) / ( 1.0_DP - exp( - xyz_FactBl(i,j,k) * 2.0_DP * DelTime ) )
            else
!!$              xyz_FactAl2(i,j,k) = - xyz_FactAl1(i,j,k)
              xyz_FactAl2(i,j,k) = - (   xyz_QCloudWatB(i,j,k) + xyz_FactAl1(i,j,k) * ( 2.0_DP * DelTime ) ) / ( 2.0_DP * DelTime )

!!$              call MessageNotify( 'E', module_name, 'Unexpected negative QCloudWat.' )
            end if
            xyz_QCloudWat(i,j,k) = 0.0_DP
          end if
        end do
      end do
      !
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_FactBs(i,j,k) >= FactBlsThreshold ) then
            xyz_QCloudIce(i,j,k) = xyz_QCloudIce(i,j,k) * exp( - xyz_FactBs(i,j,k) * 2.0_DP * DelTime ) + xyz_FactAs(i,j,k) / xyz_FactBs(i,j,k) * ( 1.0_DP - exp( - xyz_FactBs(i,j,k) * 2.0_DP * DelTime ) )
          else
            xyz_QCloudIce(i,j,k) = xyz_QCloudIce(i,j,k) + xyz_FactAs(i,j,k) * ( 2.0_DP * DelTime )
          end if
        end do
      end do
      !   The value of cloud water amount is checked, and xyz_FactA is updated.
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_QCloudIce(i,j,k) < 0.0_DP ) then
            if ( xyz_FactBs(i,j,k) >= FactBlsThreshold ) then
              xyz_FactAs2(i,j,k) = - xyz_FactAs1(i,j,k) - xyz_FactBs(i,j,k) * xyz_QCloudIceB(i,j,k) * exp( - xyz_FactBs(i,j,k) * 2.0_DP * DelTime ) / ( 1.0_DP - exp( - xyz_FactBs(i,j,k) * 2.0_DP * DelTime ) )
            else
!!$              xyz_FactAs2(i,j,k) = - xyz_FactAs1(i,j,k)
              xyz_FactAs2(i,j,k) = - (   xyz_QCloudIceB(i,j,k) + xyz_FactAs1(i,j,k) * ( 2.0_DP * DelTime ) ) / ( 2.0_DP * DelTime )

!!$              call MessageNotify( 'E', module_name, 'Unexpected negative QCloudIce.' )
            end if
            xyz_QCloudIce(i,j,k) = 0.0_DP
          end if
        end do
      end do


      ! Al and As are updated because Al2 and As2 were updated above
      xyz_FactAl(:,:,k) = xyz_FactAl1(:,:,k) + xyz_FactAl2(:,:,k)
      xyz_FactAs(:,:,k) = xyz_FactAs1(:,:,k) + xyz_FactAs2(:,:,k)



      xyz_CloudCover(:,:,k) = xyz_CloudCover(:,:,k) * exp( - ( xyz_FactC(:,:,k) + xyz_FactD(:,:,k) ) * 2.0_DP * DelTime ) + ( xyz_FactC(:,:,k) + xyz_FactE(:,:,k) ) / ( xyz_FactC(:,:,k) + xyz_FactD(:,:,k) + 1.0d-100 ) * ( 1.0_DP - exp( - ( xyz_FactC(:,:,k) + xyz_FactD(:,:,k) ) * 2.0_DP * DelTime ) )


!!$      xyz_QH2OVap(:,:,k) = xyz_QH2OVap(:,:,k) &
!!$        & - xyz_FactA2(:,:,k) * 2.0_DP * DelTime
      xyz_QH2OVap(:,:,k) = xyz_QH2OVap(:,:,k) - ( xyz_FactAl2(:,:,k) + xyz_FactAs2(:,:,k) ) * 2.0_DP * DelTime
      !
!!$      xyz_Temp(:,:,k) = xyz_Temp(:,:,k) &
!!$        & + xyz_FactA2(:,:,k) * 2.0_DP * DelTime * LatentHeat / CpDry
      xyz_Temp(:,:,k) = xyz_Temp(:,:,k) + (   xyz_FactAl2(:,:,k) * LatentHeat + xyz_FactAs2(:,:,k) * LatentHeatSubl ) * 2.0_DP * DelTime / CpDry


      ! Rain
      !
      do j = 1, jmax
        do i = 0, imax-1
!!$          if ( xyz_FactBl(i,j,k) >= FactBlsThreshold ) then
!!$            xyz_Rain(i,j,k) =                                                 &
!!$              &   xyz_FactAl(i,j,k) * 2.0_DP * DelTime                        &
!!$              & + xyz_QCloudWatB(i,j,k)                                       &
!!$              &    * ( 1.0_DP - exp( - xyz_FactBl(i,j,k) * 2.0_DP * DelTime ) )  &
!!$              & - xyz_FactAl(i,j,k) / xyz_FactBl(i,j,k)                       &
!!$              &    * ( 1.0_DP - exp( - xyz_FactBl(i,j,k) * 2.0_DP * DelTime ) )
!!$            xyz_Rain(i,j,k) = xyz_Rain(i,j,k) / ( 2.0_DP * DelTime )
!!$          else
!!$            xyz_Rain(i,j,k) = 0.0_DP
!!$          end if
          xyz_Rain(i,j,k) = xyz_QCloudWatB(i,j,k) + xyz_FactAl(i,j,k) * 2.0_DP * DelTime - xyz_QCloudWat (i,j,k)
          xyz_Rain(i,j,k) = xyz_Rain(i,j,k) / ( 2.0_DP * DelTime )
        end do
      end do

      ! Snow
      !
      do j = 1, jmax
        do i = 0, imax-1
!!$          if ( xyz_FactBs(i,j,k) >= FactBlsThreshold ) then
!!$            xyz_Snow(i,j,k) =                                               &
!!$              &   xyz_FactAs(i,j,k) * 2.0_DP * DelTime                      &
!!$              & + xyz_QCloudIceB  (i,j,k)                                   &
!!$              &    * ( 1.0_DP - exp( - xyz_FactBs(i,j,k) * 2.0_DP * DelTime ) )  &
!!$              & - xyz_FactAs(i,j,k) / xyz_FactBs(i,j,k)                     &
!!$              &    * ( 1.0_DP - exp( - xyz_FactBs(i,j,k) * 2.0_DP * DelTime ) )
!!$            xyz_Snow(i,j,k) = xyz_Snow(i,j,k) / ( 2.0_DP * DelTime )
!!$          else
!!$            xyz_Snow(i,j,k) = 0.0_DP
!!$          end if
          xyz_Snow(i,j,k) = xyz_QCloudIceB  (i,j,k) + xyz_FactAs(i,j,k) * 2.0_DP * DelTime - xyz_QCloudIce   (i,j,k)
          xyz_Snow(i,j,k) = xyz_Snow(i,j,k) / ( 2.0_DP * DelTime )
        end do
      end do

      ! Rain at the surface
      xy_Rain = xy_Rain + xyz_Rain(:,:,k) * xyz_DelMass(:,:,k)
      ! Snow at the surface
      xy_Snow = xy_Snow + xyz_Snow(:,:,k) * xyz_DelMass(:,:,k)


      ! Treatment of supersaturation
      do j = 1, jmax
        do i = 0, imax-1
          QCloudWatTentative = xyz_QCloudWat(i,j,k)
          QCloudIceTentative = xyz_QCloudIce(i,j,k)
          call LScaleCond1Grid( xyz_Temp(i,j,k), xyz_QH2OVap(i,j,k), QCloudWatTentative, QCloudIceTentative, xyz_Press(i,j,k), xyr_Press(i,j,k-1), xyr_Press(i,j,k), xyz_DQH2OLiqDtLSC(i,j,k), xyz_DQH2OSolDtLSC(i,j,k) )

!!$          xy_Rain(i,j) = xy_Rain(i,j) &
!!$            & + ( xyz_QCloudWat(i,j,k) - QCloudWatTentative ) &
!!$            &     / ( 2.0_DP * DelTime ) &
!!$            &     * ( 2.0_DP * DelTime ) &
!!$            &     * xyz_DelMass(i,j,k)
!!$          xy_Snow(i,j) = xy_Snow(i,j) &
!!$            & + ( xyz_QCloudIce(i,j,k) - QCloudIceTentative ) &
!!$            &     / ( 2.0_DP * DelTime ) &
!!$            &     * ( 2.0_DP * DelTime ) &
!!$            &     * xyz_DelMass(i,j,k)
!!$          xyz_QCloudWat(i,j,k) = xyz_QCloudWat(i,j,k) &
!!$            & - ( xyz_QCloudWat(i,j,k) - QCloudWatTentative ) &
!!$            &     / ( 2.0_DP * DelTime ) &
!!$            &     * ( 2.0_DP * DelTime )
!!$          xyz_QCloudIce(i,j,k) = xyz_QCloudIce(i,j,k) &
!!$            & - ( xyz_QCloudIce(i,j,k) - QCloudIceTentative ) &
!!$            &     / ( 2.0_DP * DelTime ) &
!!$            &     * ( 2.0_DP * DelTime )
        end do
      end do

      xy_Rain = xy_Rain + xyz_DQH2OLiqDtLSC(:,:,k) * xyz_DelMass(:,:,k)
      xy_Snow = xy_Snow + xyz_DQH2OSolDtLSC(:,:,k) * xyz_DelMass(:,:,k)


      ! Evaporation
      !
      do j = 1, jmax
        do i = 0, imax-1
          if ( ( xyz_QCloudWat(i,j,k) < QCloudWaterEvapThreshold ) .and. ( xyz_QCloudIce(i,j,k) < QCloudWaterEvapThreshold ) ) then
            xyz_CloudCover(i,j,k) = 0.0_DP
          end if
        end do
      end do

      ! Cloud cover is restricted.
      xyz_CloudCover(:,:,k) = min( max( xyz_CloudCover(:,:,k), 0.0_DP ), 1.0_DP )

      ! Check values
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_QH2OVap(i,j,k) < 0.0_DP ) then
            write( 6, * ) 'QH2OVap is negative', i, j, k, xyz_QH2OVap(i,j,k)
          end if
          if ( xyz_QCloudWat(i,j,k) < 0.0_DP ) then
            write( 6, * ) 'QCloudWat is negative', i, j, k, xyz_QCloudWat(i,j,k)
          end if
          if ( xyz_QCloudIce  (i,j,k) < 0.0_DP ) then
            write( 6, * ) 'QCloudIce is negative', i, j, k, xyz_QCloudIce  (i,j,k)
          end if
        end do
      end do

    end do k_loop


    ! 大規模凝結 (非対流性凝結) (Manabe, 1965)
    ! Large scale condensation (non-convective condensation)
    ! (Manabe, 1965)
    !
    ! It should be noted that H2OLiq and H2OSol have updated in above
    ! subroutine.
    !
    ! This routine is commented out, since this is done in an above k_loop.

!!$    call LScaleCond1D3DWrapper(                            &
!!$      & xyz_Temp, xyz_QH2OVap,          & ! (inout)
!!$      & xyz_QCloudWat, & ! (inout)
!!$      & xyz_QCloudIce, & ! (inout)
!!$      & xyz_Press, xyr_Press,                              & ! (in)
!!$      & xyz_DQH2OLiqDtLSC, xyz_DQH2OSolDtLSC               & ! (out)
!!$      & )


    xy_SurfRainFlux = xy_Rain
    xy_SurfSnowFlux = xy_Snow


    call CloudUtilConsChk( "CloudT1993baseWithIce", xyr_Press, xyz_TempB, xyz_QH2OVapB, xyz_QCloudWatB, xyz_QCloudIceB, xyz_Temp , xyz_QH2OVap , xyz_QCloudWat , xyz_QCloudIce , xy_SurfRainFlux, xy_SurfSnowFlux )


  end subroutine CloudT1993baseWithIce