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

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


    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_QCloudWaterConvThreshold(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

    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
      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) >= CloudProdRHThreshold * 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_QCloudWaterConvThreshold = QCloudWaterEffConv / ( 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_QCloudWaterConvThreshold ) )**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_DTempDtLSC, xyz_DQVapDtLSC, 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 CloudUtilsCalcPRCPKeyLLTemp( 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

    ! 大規模凝結 (非対流性凝結)
    ! 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/ CloudProdRHThreshold, CloudLifeTime0, QCloudWaterEffConv, TempBFEffectSat
          !
          ! デフォルト値については初期化手続 "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
    !
    CloudProdRHThreshold = 0.5_DP
    CloudLifeTime0       = 1000.0_DP
    QCloudWaterEffConv   = 4.0d-4
    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.

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

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


    ! ヒストリデータ出力のためのへの変数登録
    ! 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, 'CloudProdRHThreshold = %f', d = (/ CloudProdRHThreshold /) )
    call MessageNotify( 'M', module_name, 'CloudLifeTime0       = %f', d = (/ CloudLifeTime0 /) )
    call MessageNotify( 'M', module_name, 'QCloudWaterEffConv   = %f', d = (/ QCloudWaterEffConv /) )
    call MessageNotify( 'M', module_name, 'TempBFEffectSat      = %f', d = (/ TempBFEffectSat /) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )


    cloud_T1993base_inited = .true.

  end subroutine CloudT1993baseInit