!= 簡単雲モデル
!
!= Simple cloud model
!
! Authors::   Yoshiyuki O. Takahashi
! Version::   $Id: cloud_simple.f90,v 1.9 2015/01/29 12:06:43 yot Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module cloud_simple
  !
  != 簡単雲モデル
  !
  != Simple cloud model
  !
  ! <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_simple_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 :: CloudSimple
  public :: CloudSimpleCalcPRCPKeyLLTemp
  public :: CloudSimpleCalcPRCPKeyLLTemp3D
  public :: CloudSimpleCalcPRCPStepPC
  public :: CloudSimpleWithIce
  public :: CloudSimpleCalcCloudCover
  public :: CloudSimpleInit

  public :: CloudSimpleDivideWatAndIce


  ! 公開変数
  ! Public variables
  !


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

  integer , save        :: IDCloudCoverMethod
  integer , parameter   :: IDCloudCoverMethodConst    = 1
  integer , parameter   :: IDCloudCoverMethodRH       = 2
  integer , parameter   :: IDCloudCoverMethodRHLin    = 3

  integer , save        :: IDSnowMethod
  integer , parameter   :: IDSnowMethodKeyLLTemp      = 11
  integer , parameter   :: IDSnowMethodStepPC         = 12

  real(DP), save        :: CloudLifeTime
  real(DP), save        :: CloudWatLifeTime
  real(DP), save        :: CloudIceLifeTime

  real(DP), save        :: CloudCover
  real(DP), save        :: CloudCoverRHCrtl
  real(DP), save        :: CloudCoverMin


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


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

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

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

contains

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

  subroutine CloudSimple(                             &
    & xyr_Press, xyz_Press,                           & ! (in)
    & xyz_Temp,                                       & ! (inout)
!!$    & xyz_DQH2OLiqDtCum, xyz_DQH2OLiqDtLSC,           & ! (in)
    & xyz_QH2OVap, xyz_QH2OLiq,                       & ! (inout)
    & xy_Rain, xy_Snow                                & ! (out)
    & )

    ! USE statements
    !

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


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


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

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

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

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

    ! 実行文 ; Executable statement
    !

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


    ! tentative treatment
    xyz_DQH2OLiqDtCum = 0.0_DP
    xyz_DQH2OLiqDtLSC = 0.0_DP


    ! Numerical solution

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


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

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


    ! Analytical solution

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

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

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



    select case ( IDSnowMethod )
    case ( IDSnowMethodKeyLLTemp )

      call CloudSimpleCalcPRCPKeyLLTemp3D(      &
        & xyr_Press, xyz_Press, xyz_DQRainDt,   &  ! (in )
        & xyz_Temp, xyz_QH2OVap,                &  ! (in )
        & xy_Rain, xy_Snow                      &  ! (out)
        & )

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

    case ( IDSnowMethodStepPC )

      xyz_DQSnowDt = 0.0_DP

      call CloudSimpleCalcPRCPStepPC(      &
        & xyr_Press, xyz_Press,                  & ! (in   )
        & xyz_DQRainDt, xyz_DQSnowDt,            & ! (in   )
        & xyz_Temp, xyz_QH2OVap,                 & ! (inout)
        & xy_Rain, xy_Snow                       & ! (out  )
        & )

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

    end select


  end subroutine CloudSimple

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

  subroutine CloudSimpleCalcPRCPKeyLLTemp(  &
    & xyz_Temp, xy_PRCP,                   &  ! (in )
    & xy_SurfRainFlux, xy_SurfSnowFlux     &  ! (out)
    & )


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


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


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


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


    if ( FlagSnow ) then

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

    else

      xy_SurfRainFlux = xy_PRCP
      xy_SurfSnowFlux = 0.0_DP

    end if


  end subroutine CloudSimpleCalcPRCPKeyLLTemp

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

  subroutine CloudSimpleCalcPRCPKeyLLTemp3D(  &
    & xyr_Press, xyz_Press, xyz_DQH2OLiqDt,   &  ! (in )
    & xyz_Temp, xyz_QH2OVap,                  &  ! (in )
    & xy_SurfRainFlux, xy_SurfSnowFlux        &  ! (out)
    & )

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only:  &
      & CpDry,            &
                              ! $ C_p $ [J kg-1 K-1].
                              ! 乾燥大気の定圧比熱.
                              ! Specific heat of air at constant pressure
      & Grav,             &
                              ! $ g $ [m s-2].
                              ! 重力加速度.
                              ! Gravitational acceleration
      & LatentHeat,       &
                              ! $ L $ [J kg-1] .
                              ! 蒸発の潜熱.
                              ! Latent heat
      & 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_CalcQVapSatOnLiq

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


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

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

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

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


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


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


    if ( FlagPRCPEvap ) then

      xyz_DQRainDt = xyz_DQH2OLiqDt

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

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

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

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

          end do
        end do
      end do

      xy_PRCP = xy_SurfRainFlux

    else

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

    end if


    call CloudSimpleCalcPRCPKeyLLTemp(     &
      & xyz_Temp, xy_PRCP,                 &  ! (in )
      & xy_SurfRainFlux, xy_SurfSnowFlux   &  ! (out)
      & )


  end subroutine CloudSimpleCalcPRCPKeyLLTemp3D

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

  subroutine CloudSimpleWithIce(                      &
    & xyr_Press, xyz_Press,                           & ! (in)
    & xyz_Temp,                                       & ! (inout)
!!$    & xyz_DQH2OLiqDtCum, xyz_DQH2OLiqDtLSC,           & ! (in)
    & xyz_QH2OVap, xyz_QH2OLiq, xyz_QH2OSol,          & ! (inout)
    & xy_Rain, xy_Snow                                & ! (out)
    & )

    ! USE statements
    !

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

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

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


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


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

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

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


    ! 実行文 ; Executable statement
    !

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

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


    xyz_DQH2OLiqDt = 0.0_DP
    xyz_DQH2OSolDt = 0.0_DP

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


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

    k_loop : do k = kmax, 1, -1

      ! Freezing/melting and evaporation of precipitation
      !
      do j = 1, jmax
        do i = 0, imax-1
          call CloudUtilsPRCPStepPC1Grid(  &
            & xyr_Press(i,j,k-1), xyr_Press(i,j,k),       & ! (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_QH2OLiq(:,:,k) =                                                    &
        &   xyz_QH2OLiq(:,:,k)                                                &
        &     * exp( - 2.0_DP * DelTime / ( CloudWatLifeTime + 1.0e-100_DP ) )&
        & + xyz_DQH2OLiqDt(:,:,k) * CloudWatLifeTime                          &
        &   * ( 1.0_DP - exp( - 2.0_DP * DelTime / ( CloudWatLifeTime + 1.0e-100_DP ) ) )

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


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

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

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


    end do k_loop


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


  end subroutine CloudSimpleWithIce

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

  subroutine CloudSimpleWithIceOld(                   &
    & xyr_Press, xyz_Press,                           & ! (in)
    & xyz_Temp,                                       & ! (inout)
!!$    & xyz_DQH2OLiqDtCum, xyz_DQH2OLiqDtLSC,           & ! (in)
    & xyz_QH2OVap, xyz_QH2OLiq, xyz_QH2OSol,          & ! (inout)
    & xy_Rain, xy_Snow                                & ! (out)
    & )

    ! USE statements
    !

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


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


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

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

    ! 実行文 ; Executable statement
    !

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


    xyz_DQH2OLiqDt = 0.0_DP
    xyz_DQH2OSolDt = 0.0_DP


    ! Numerical solution

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


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

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


    ! Analytical solution

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


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

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


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

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



    select case ( IDSnowMethod )
    case ( IDSnowMethodKeyLLTemp )

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

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

    case ( IDSnowMethodStepPC )

      call CloudSimpleCalcPRCPStepPC(      &
        & xyr_Press, xyz_Press,                  & ! (in   )
        & xyz_DQRainDt, xyz_DQSnowDt,            & ! (in   )
        & xyz_Temp, xyz_QH2OVap,                 & ! (inout)
        & xy_Rain, xy_Snow                       & ! (out  )
        & )

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

    end select


  end subroutine CloudSimpleWithIceOld

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

  subroutine CloudSimpleCalcPRCPStepPC(      &
    & xyr_Press, xyz_Press,                  & ! (in   )
    & xyz_DQRainDt, xyz_DQSnowDt,            & ! (in   )
    & xyz_Temp, xyz_QH2OVap,                 & ! (inout)
    & xy_SurfRainFlux, xy_SurfSnowFlux       & ! (out  )
    & )

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only:  &
      & CpDry,            &
                              ! $ C_p $ [J kg-1 K-1].
                              ! 乾燥大気の定圧比熱.
                              ! Specific heat of air at constant pressure
      & Grav,             &
                              ! $ g $ [m s-2].
                              ! 重力加速度.
                              ! Gravitational acceleration
      & LatentHeat,       &
                              ! $ L $ [J kg-1] .
                              ! 蒸発の潜熱.
                              ! Latent heat
      & LatentHeatFusion, &
                              ! $ L $ [J kg-1] .
                              ! 融解の潜熱.
                              ! Latent heat of fusion
      & EpsV
                              ! $ \epsilon_v $ .
                              ! 水蒸気分子量比.
                              ! Molecular weight of water vapor

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

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


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


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

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

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


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


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

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

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

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


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

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

        end do
      end do
    end do


  end subroutine CloudSimpleCalcPRCPStepPC

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

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

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: &
      & Grav, &
                              ! $ g $ [m s-2].
                              ! 重力加速度.
                              ! Gravitational acceleration
      & GasRDry
                              ! $ R $ [J kg-1 K-1].
                              ! 乾燥大気の気体定数.
                              ! Gas constant of air


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


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

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

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

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

    real(DP) :: DelZ


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

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

    DelZ = DelMass / Dens


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

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

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

    DelPRCPFlux = min( DelPRCPFlux, PRCP )


  end subroutine CloudSimpleEvap1Grid

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

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


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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: &
      & Grav, & 
                              ! $ g $ [m s-2]. 
                              ! 重力加速度. 
                              ! Gravitational acceleration
      & CpDry, &
                              ! $ C_p $ [J kg-1 K-1]. 
                              ! 乾燥大気の定圧比熱. 
                              ! Specific heat of air at constant pressure
      & LatentHeat, &
                              ! $ L $ [J kg-1] . 
                              ! 凝結の潜熱. 
                              ! Latent heat of condensation
      & LatentHeatFusion
                              ! $ L $ [J kg-1] .
                              ! 融解の潜熱.
                              ! Latent heat of fusion

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

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


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

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

    xy_SumB = xy_Sum

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

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



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

    xy_SumB = xy_Sum

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

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


  end subroutine CloudSimpleConsChk

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

  subroutine CloudSimpleCalcCloudCover(  &
    & xyz_Press, xyz_Temp, xyz_QH2OTot,  & ! (in )
    & xyz_CloudCover                     & ! (out)
    & )

    ! USE statements
    !

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

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


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

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


    ! 実行文 ; Executable statement
    !

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


    select case ( IDCloudCoverMethod )
    case ( IDCloudCoverMethodConst )

      xyz_CloudCover = CloudCover

    case ( IDCloudCoverMethodRH )

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

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

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

    case ( IDCloudCoverMethodRHLin )

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

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

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

    end select


  end subroutine CloudSimpleCalcCloudCover

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

  subroutine CloudSimpleDivideWatAndIce(  &
    & xyz_Temp,                           & ! (in )
    & xyz_QH2OWatAndIce,                  & ! (in )
    & xyz_QH2OWat, xyz_QH2OIce            & ! (out)
    & )

    ! USE statements
    !

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

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


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

    ! 実行文 ; Executable statement
    !

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


    call SaturateWatFraction(      &
      & xyz_Temp,                  & ! (in )
      & xyz_WatFrac                & ! (out)
      & )

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


  end subroutine CloudSimpleDivideWatAndIce

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

  subroutine CloudSimpleInit( &
    & ArgFlagSnow             &
    & )

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

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

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

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

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

    ! 宣言文 ; Declaration statements
    !

    logical, intent(in) :: ArgFlagSnow


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

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

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

    ! 実行文 ; Executable statement
    !

    if ( cloud_simple_inited ) return


    FlagSnow = ArgFlagSnow


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

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

    SnowMethod          = 'KeyLLTemp'

    CloudCover          = 1.0_DP

    CloudCoverRHCrtl    = 0.0_DP

    CloudCoverMin       = 0.0_DP

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


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

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

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


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


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



    ! Initialization of modules used in this module
    !

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

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


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



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


    cloud_simple_inited = .true.

  end subroutine CloudSimpleInit

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

end module cloud_simple
