!= Tiedtke (1993) に基づく雲モデル
!
!= Cloud model based on Tiedtke (1993)
!
! Authors::   Yoshiyuki O. Takahashi
! Version::   $Id: cloud_T1993base.f90,v 1.6 2015/01/29 12:20:40 yot Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module cloud_T1993base
  !
  != Tiedtke (1993) に基づく雲モデル
  !
  != Cloud model based on Tiedtke (1993)
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! 簡単雲モデルによる雲の計算.
  !
  ! In this module, the amount of cloud is calculated by use of a simple
  ! cloud model.
  !
  !== Procedures List
  !
!!$  ! RadiationFluxDennouAGCM :: 放射フラックスの計算
!!$  ! ------------            :: ------------
!!$  ! RadiationFluxDennouAGCM :: Calculate radiation flux
  !
  !== NAMELIST
  !
  ! NAMELIST#cloud_T1993base_nml
  !

  ! モジュール引用 ; USE statements

  !
  ! Kind type parameter
  !
  use dc_types, only: DP, &      ! Double precision.
    &                 STRING, &  ! Strings.
    &                 TOKEN      ! Keywords.

  ! メッセージ出力
  ! Message output
  !
  use dc_message, only: MessageNotify

  ! 格子点設定
  ! Grid points settings
  !
  use gridset, only: imax, & ! 経度格子点数.
                             ! Number of grid points in longitude
    &                jmax, & ! 緯度格子点数.
                             ! Number of grid points in latitude
    &                kmax    ! 鉛直層数.
                             ! Number of vertical level

  implicit none

  private


  ! 公開手続き
  ! Public procedure
  !
  public :: CloudT1993base
  public :: CloudT1993baseWithIce
  public :: CloudT1993baseInit


  ! 公開変数
  ! Public variables
  !


  ! 非公開変数
  ! Private variables
  !
  logical , save        :: FlagSnow
                           ! A flag for snow

  real(DP), save        :: RHThresholdCrtl
  real(DP), save        :: RHThresholdSigmaMin
  real(DP), save        :: RHThresholdOrd

  real(DP), save        :: CloudLifeTime0

  real(DP), save        :: CloudWatLifeTime0
  real(DP), save        :: CloudIceLifeTime0

  real(DP), save        :: QCloudWatEffConv0
  real(DP), save        :: QCloudIceEffConv0

  real(DP), save        :: TempBFEffectSat
                           ! Temperature at which Bergeron-Findeisen effect 
                           ! saturates.


  real(DP), save :: FactBlsThreshold = 1.0e-10_DP

  real(DP), save :: PRCPArea
    !                           a_p
  real(DP), save :: PRCPEvapArea
    !                           A = max( a_p - a, 0 )

  real(DP), save :: PRCPColFactor
  !                             This is C1 in Sundqvist et al. (1989)


  logical, save :: cloud_T1993base_inited = .false.
                              ! 初期設定フラグ.
                              ! Initialization flag

  character(*), parameter:: module_name = 'cloud_T1993base'
                              ! モジュールの名称.
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: cloud_T1993base.f90,v 1.6 2015/01/29 12:20:40 yot Exp $'
                              ! モジュールのバージョン
                              ! Module version

  !--------------------------------------------------------------------------------------

contains

  !----------------------------------------------------------------------------
  ! This subroutine, SetCloudCloudWaterAdjustWithCloudCoverV4, is 
  ! created at 2012/08/13.

!!$  call CloudT1993base(                                        &
!!$    & xyz_Press, xyr_Press, xyz_VirTemp,                      & ! (in)
!!$    & xyz_DQH2OLiqDtCum,                                      & ! (in)
!!$    & xyz_MoistConvDetTend,                                   & ! (in)
!!$    & xyz_OMG, xyz_MoistConvSubsidMassFlux, xyz_DTempDt,      & ! (in)
!!$    & xyz_TempA, xyzf_QMixA(:,:,:,IndexH2OVap),               & ! (inout)
!!$    & xyzf_QMixA(:,:,:,CompositionInqIndex('H2OLiq')),        & ! (inout)
!!$    & xyzf_QMixA(:,:,:,CompositionInqIndex('CloudCoverTMP')), & ! (inout)
!!$    & xy_Rain, xy_Snow                                        & ! (out)
!!$    & )


  subroutine CloudT1993base(                                    &
    & xyz_Press, xyr_Press, xyz_VirTemp,                        & ! (in)
    & xyz_DQCloudWaterDtCum,                                    & ! (in)
    & xyz_MoistConvDetTend,                                     & ! (in)
    & xyz_OMG, xyz_MoistConvSubsidMassFlux, xyz_DTempDtPhy,     & ! (in)
    & xyz_Temp, xyz_QH2OVap,                                    & ! (inout)
    & xyz_QCloudWater,                                          & ! (inout)
    & xyz_CloudCover,                                           & ! (inout)
    & xy_SurfRainFlux, xy_SurfSnowFlux                          & ! (out)
    & )

    ! 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, &
                              ! $ g $ [m s-2].
                              ! 重力加速度.
                              ! Gravitational acceleration
      & CpDry, &
                              ! $ C_p $ [J kg-1 K-1].
                              ! 乾燥大気の定圧比熱.
                              ! Specific heat of air at constant pressure
      & GasRDry, &
                              ! $ R $ [J kg-1 K-1].
                              ! 乾燥大気の気体定数.
                              ! Gas constant of air
      & GasRWet, &
                              ! $ R_v $ [J kg-1 K-1].
                              ! 凝結成分の気体定数.
                              ! Gas constant of condensible elements
      & 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 * gamma( 13.0d0/5.0d0 ) &
      & 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,              & ! (inout)
      & xyz_Press, xyr_Press,               & ! (in)
      & xyz_RainLSC                         & ! (out)
      & )

    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,                 &  ! (in )
      & xy_SurfRainFlux, xy_SurfSnowFlux   &  ! (out)
      & )


  end subroutine CloudT1993base

  !----------------------------------------------------------------------------

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

    ! 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, &
                              ! $ g $ [m s-2].
                              ! 重力加速度.
                              ! Gravitational acceleration
      & CpDry, &
                              ! $ C_p $ [J kg-1 K-1].
                              ! 乾燥大気の定圧比熱.
                              ! Specific heat of air at constant pressure
      & GasRDry, &
                              ! $ R $ [J kg-1 K-1].
                              ! 乾燥大気の気体定数.
                              ! Gas constant of air
      & GasRWet, &
                              ! $ R_v $ [J kg-1 K-1].
                              ! 凝結成分の気体定数.
                              ! Gas constant of condensible elements
      & LatentHeat, &
                              ! $ L $ [J kg-1] .
                              ! 凝結の潜熱.
                              ! Latent heat of condensation
      & LatentHeatFusion, &
                              ! $ L $ [J kg-1] .
                              ! 融解の潜熱.
                              ! Latent heat of fusion
      & 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: &
!!$      & LScaleCond1D3DWrapper, &
      & 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 * gamma( 13.0d0/5.0d0 ) &
      & 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,             & ! (in )
      & xyz_WatFrac           & ! (out)
      & )
    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),       & ! (in   )
            & xyz_Temp(i,j,k),                            & ! (inout)
            & xy_Rain(i,j), xy_Snow(i,j)                  & ! (inout)
            & )
        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), & ! (in)
            & PRCPArea, PRCPEvapArea,                                 & ! (in)
            & xyz_Temp(i,j,k), xyz_QH2OVap(i,j,k),                    & ! (inout)
            & xy_Rain(i,j), xy_Snow(i,j)                              & ! (inout)
            & )
        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), &  ! (inout)
!!$            & xyz_QCloudWat(i,j,k), xyz_QCloudIce(i,j,k),           & ! (inout)
            & QCloudWatTentative, QCloudIceTentative,               & ! (inout)
            & xyz_Press(i,j,k), xyr_Press(i,j,k-1), xyr_Press(i,j,k), & ! (in)
            & xyz_DQH2OLiqDtLSC(i,j,k), xyz_DQH2OSolDtLSC(i,j,k)      & ! (out)
            & )

!!$          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",                                 &
!!$      & FlagIncludeIcePhaseChange,                               &
!!$      & .true.,                                                  &
      & 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

  !--------------------------------------------------------------------------------------

  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, &          ! (out)
        & namelist_filename, mode = 'r' ) ! (in)

      rewind( unit_nml )
      read( unit_nml,                     & ! (in)
        & nml = cloud_T1993base_nml,      & ! (out)
        & 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

  !--------------------------------------------------------------------------------------

end module cloud_T1993base