subroutine RelaxedArakawaSchubert( xy_SurfTemp, xyz_Press, xyr_Press, xyz_Exner, xyr_Exner, xyz_Temp, xyz_QH2OVap, xyz_DTempDt, xyz_DQVapDt, xyz_DQH2OLiqDt, xyz_MoistConvDetTend, xyz_MoistConvSubsidMassFlux )
    !
    ! relaxed Arakawa-Schubert スキームにより, 温度と比湿を変化. 
    !
    ! Change temperature and specific humidity by relaxed Arakawa-Schubert scheme
    !

    ! モジュール引用 ; USE statements
    !

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

    ! 時刻管理
    ! Time control
    !
    use timeset, only: DelTime, TimeN, TimesetClockStart, TimesetClockStop

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

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

    ! Arakawa-Schubert scheme by Lord et al. (1982)
    ! Arakawa-Schubert scheme by Lord et al. (1982)
    !
    use arakawa_schubert_L1982, only : ArakawaSchubertL1982CalcCWFCrtl


    ! 宣言文 ; Declaration statements
    !

    real(DP), intent(in   ) :: xy_SurfTemp (0:imax-1, 1:jmax)
                              ! Pressure
    real(DP), intent(in   ) :: xyz_Press   (0:imax-1, 1:jmax, 1:kmax)
                              ! Pressure
    real(DP), intent(in   ) :: xyr_Press   (0:imax-1, 1:jmax, 0:kmax)
                              ! Pressure
    real(DP), intent(in   ) :: xyz_Exner   (0:imax-1, 1:jmax, 1:kmax)
                              ! Exner function
    real(DP), intent(in   ) :: xyr_Exner   (0:imax-1, 1:jmax, 0:kmax)
                              ! Exner function
    real(DP), intent(inout) :: xyz_Temp    (0:imax-1, 1:jmax, 1:kmax)
                              ! Temperature
    real(DP), intent(inout) :: xyz_QH2OVap (0:imax-1, 1:jmax, 1:kmax)
                              ! $ q $ .     比湿. Specific humidity
!!$    real(DP), intent(inout) :: xy_Rain (0:imax-1, 1:jmax)
!!$                              ! 降水量. 
!!$                              ! Precipitation
    real(DP), intent(inout) :: xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax)
                              ! 温度変化率. 
                              ! Temperature tendency
    real(DP), intent(inout) :: xyz_DQVapDt (0:imax-1, 1:jmax, 1:kmax)
                              ! 比湿変化率. 
                              ! Specific humidity tendency

    real(DP), intent(out  ) :: xyz_DQH2OLiqDt(0:imax-1, 1:jmax, 1:kmax)

    real(DP), intent(out  ), optional :: xyz_MoistConvDetTend       (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(out  ), optional :: xyz_MoistConvSubsidMassFlux(0:imax-1, 1:jmax, 1:kmax)

    ! 作業変数
    ! Work variables
    !
    real(DP) :: xyz_Height  (0:imax-1, 1:jmax, 1:kmax)
                              !
                              ! Height
    real(DP) :: xyr_Height  (0:imax-1, 1:jmax, 0:kmax)
                              !
                              ! Height
    real(DP) :: xy_RainCumulus (0:imax-1, 1:jmax)
                              ! 降水量. 
                              ! Precipitation
    real(DP) :: xyz_DTempDtCumulus (0:imax-1, 1:jmax, 1:kmax)
                              ! 温度変化率. 
                              ! Temperature tendency
    real(DP) :: xyz_DQVapDtCumulus (0:imax-1, 1:jmax, 1:kmax)
                              ! 比湿変化率. 
                              ! Specific humidity tendency

    real(DP) :: xyz_DelPress(0:imax-1, 1:jmax, 1:kmax)
                              ! $ \Delta p $
                              !
    real(DP) :: xyz_PotTemp (0:imax-1, 1:jmax, 1:kmax)
                              ! Potential temperature
                              !
    real(DP) :: xyz_QH2OVapSat       (0:imax-1, 1:jmax, 1:kmax)
                              ! 飽和比湿. 
                              ! Saturation specific humidity. 

    ! Dry and moist static energy in environment (Env) and cloud (Cld)
    !
    real(DP) :: xyz_EnvDryStaticEne     (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyr_EnvDryStaticEne     (0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xyz_EnvMoistStaticEne   (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyr_EnvMoistStaticEne   (0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xyz_EnvMoistStaticEneSat(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyr_EnvMoistStaticEneSat(0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xyr_CldMoistStaticEne   (0:imax-1, 1:jmax, 0:kmax)

    real(DP) :: xy_Kernel               (0:imax-1, 1:jmax)
                   ! Tendency of cloud work function by cumulus convection, kernel
    real(DP) :: xy_CWF                  (0:imax-1, 1:jmax)
                   ! Cloud work function
    real(DP) :: xyz_CWF                 (0:imax-1, 1:jmax, 1:kmax)
                   ! Cloud work function
                   ! (variable for output)
    real(DP) :: xy_DCWFDtLS             (0:imax-1, 1:jmax)
                   ! Tendency of cloud work function by large scale motion
    real(DP) :: xyz_DCWFDtLS            (0:imax-1, 1:jmax, 1:kmax)
                   ! Tendency of cloud work function by large scale motion
                   ! (variable for output)
    real(DP) :: xy_CldMassFluxBottom    (0:imax-1, 1:jmax)
                   ! Cloud mass flux at cloud bottom

    real(DP) :: xyz_Beta                (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_BetaCldTop          (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_Gamma               (0:imax-1, 1:jmax, 1:kmax)

    real(DP) :: xyz_GammaDSE            (0:imax-1, 1:jmax, 1:kmax)
                                  ! Tendency of dry static energy per unit mass flux
    real(DP) :: xyz_GammaMSE            (0:imax-1, 1:jmax, 1:kmax)
                                  ! Tendency of moist static energy per unit mass flux

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

    real(DP) :: xy_PressCldBase         (0:imax-1, 1:jmax)
                                  ! Pressure of cloud base
    real(DP) :: xyz_CWFCrtl             (0:imax-1, 1:jmax, 1:kmax)
                                  ! "Critical value" of cloud work function
    real(DP) :: xyz_RainFactor          (0:imax-1, 1:jmax, 1:kmax)

    real(DP) :: xy_EntParam             (0:imax-1, 1:jmax)
                                  ! Entrainment factor
    real(DP) :: xyz_EntParam            (0:imax-1, 1:jmax, 1:kmax)
                                  ! Entrainment factor (variable for output)
    real(DP) :: xy_EntParamLL           (0:imax-1, 1:jmax)
                                  ! Entrainment factor for a cloud with top at one layer
                                  ! higher level
    real(DP) :: xy_EntParamUL           (0:imax-1, 1:jmax)
                                  ! Entrainment factor for a cloud with top at one layer
                                  ! lower level

    ! Difference of normalized mass flux between layer interface
    real(DP) :: xyz_DelNormMassFlux     (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xy_DelNormMassFluxCldTop(0:imax-1, 1:jmax)
    ! Normalized mass flux at layer interface and cloud top
    real(DP) :: xyr_NormMassFlux        (0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xy_NormMassFluxCldTop   (0:imax-1, 1:jmax)

    ! Liquid water at cloud top
    real(DP) :: xy_CldQH2OLiqCldTop     (0:imax-1, 1:jmax)

    ! Mass flux distribution function
    real(DP) :: xyz_MassFluxDistFunc    (0:imax-1, 1:jmax, 1:kmax)


    real(DP) :: xyz_DelH2OMass  (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xy_H2OMassB     (0:imax-1, 1:jmax)
    real(DP) :: xy_H2OMassA     (0:imax-1, 1:jmax)

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

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

    logical  :: xy_FlagCrossSatEquivPotTemp(0:imax-1, 1:jmax)
                              ! 
                              ! Flag showing whether a parcel in cloud has moist static 
                              ! energy larger than environment's

    real(DP) :: xyr_QH2OVapSat       (0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xyr_TempAdiabAscent  (0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xy_SurfPotTemp       (0:imax-1, 1:jmax)

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


    ! Variables for looking for top of mixed layer
    !
    logical  :: xy_FlagMixLayTopFound (0:imax-1, 1:jmax)
    integer  :: xy_IndexMixLayTop     (0:imax-1, 1:jmax)


    ! Variables for modification of cloud mass flux
    !
    real(DP) :: xyz_QH2OVapTentative   (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: CldMassFluxCorFactor
    real(DP) :: xy_CldMassFluxCorFactor(0:imax-1, 1:jmax)

    real(DP) :: xyz_TempB   (0:imax-1, 1:jmax, 1:kmax)
                              ! 調節前の温度.
                              ! Temperature before adjustment
    real(DP) :: xyz_QH2OVapB(0:imax-1, 1:jmax, 1:kmax)
                              ! 調節前の比湿.
                              ! Specific humidity before adjustment

    ! Flags for modification of
    !
    logical  :: xy_FlagKernelNegative (0:imax-1, 1:jmax)
    logical  :: xy_FlagNegH2OLiqCldTop(0:imax-1, 1:jmax)


    ! Variables for subsidence mass flux between updrafts
    !
    real(DP) :: DelNormMassFluxHalfLayer
    real(DP) :: NormMassFlux


    ! Variables for debug
    !
!!$    real(DP) :: xyz_DelVal(0:imax-1, 1:jmax, 1:kmax)
!!$    real(DP) :: xy_SumValB(0:imax-1, 1:jmax)
!!$    real(DP) :: xy_SumValA(0:imax-1, 1:jmax)
!!$    real(DP) :: Ratio


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


    integer :: i               ! 経度方向に回る DO ループ用作業変数
                               ! Work variables for DO loop in longitude
    integer :: j               ! 緯度方向に回る DO ループ用作業変数
                               ! Work variables for DO loop in latitude
    integer :: k               ! 鉛直方向に回る DO ループ用作業変数
                               ! Work variables for DO loop in vertical direction
    integer :: l
    integer :: m
    integer :: n



    ! 実行文 ; Executable statement
    !

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

    ! 計算時間計測開始
    ! Start measurement of computation time
    !
    call TimesetClockStart( module_name )


    ! 調節前 "Temp", "QH2OVap" の保存
    ! Store "Temp", "QH2OVap" before adjustment
    !
    xyz_TempB    = xyz_Temp
    xyz_QH2OVapB = xyz_QH2OVap


    ! Preparation of variables
    !
    !
    !   Auxiliary variables
    !     Pressure difference between upper and lower interface of the layer
    do k = 1, kmax
      xyz_DelPress(:,:,k) = xyr_Press(:,:,k-1) - xyr_Press(:,:,k)
    end do
    !     beta
    do k = 1, kmax
      xyz_Beta(:,:,k)       = CpDry / Grav * ( xyr_Exner(:,:,k-1) - xyr_Exner(:,:,k) )
    end do
    do k = 1, kmax
      xyz_BetaCldTop(:,:,k) = CpDry / Grav * ( xyr_Exner(:,:,k-1) - xyz_Exner(:,:,k) )
    end do
    !
    ! Search for top of mixed layer (lifting condensation level) based on 
    !   a description in p.684 of Arakawa and Shubert (1974). 
    !
    call RelaxedArakawaSchubertHeight( xyz_Temp, xyz_Exner, xyz_Beta, xyz_BetaCldTop, xyz_Height, xyr_Height )
    !
    !====================================
    !
!!$    xyz_TempAdiabAscent(:,:,1) = xyz_Temp(:,:,1)
!!$    do k = 2, kmax
!!$      xyz_TempAdiabAscent(:,:,k) = &
!!$        & xyz_Temp(:,:,1) - Grav / CpDry * ( xyz_Height(:,:,k) - xyz_Height(:,:,1) )
!!$    end do
!!$    xyz_TempAdiabAscent = max( xyz_TempAdiabAscent, 1.0_DP )
!!$    xyz_QH2OVapSat = xyz_CalcQVapSat( xyz_TempAdiabAscent, xyz_Press )
!!$    xy_IndexMixLayTop     = 1
!!$    xy_FlagMixLayTopFound = .false.
!!$    do k = 2, kmax
!!$      do j = 1, jmax
!!$        do i = 0, imax-1
!!$          if ( ( xyz_QH2OVap(i,j,1) >= xyz_QH2OVapSat(i,j,k) ) .and. &
!!$            &  ( .not. xy_FlagMixLayTopFound(i,j) )                 ) then
!!$            xy_IndexMixLayTop    (i,j) = k - 1
!!$            xy_FlagMixLayTopFound(i,j) = .true.
!!$          end if
!!$        end do
!!$      end do
!!$    end do
    !
    !------------------------------------
    !
!!$    xyr_TempAdiabAscent(:,:,0) = xy_SurfTemp
!!$    do k = 1, kmax
!!$      xyr_TempAdiabAscent(:,:,k) = &
!!$        & xy_SurfTemp - Grav / CpDry * ( xyr_Height(:,:,k) - 0.0_DP )
!!$    end do
!!$    xyr_TempAdiabAscent = max( xyr_TempAdiabAscent, 1.0_DP )
!!$
    xyr_TempAdiabAscent(:,:,0) = xy_SurfTemp
    xy_SurfPotTemp = xy_SurfTemp / xyr_Exner(:,:,0)
    do k = 1, kmax
      xyr_TempAdiabAscent(:,:,k) = xy_SurfPotTemp * xyr_Exner(:,:,k)
    end do
    !
    xyr_QH2OVapSat(:,:,0       ) = 1.0d100
    xyr_QH2OVapSat(:,:,1:kmax-1) = xyz_CalcQVapSat( xyr_TempAdiabAscent(:,:,1:kmax-1), xyr_Press(:,:,1:kmax-1) )
    xyr_QH2OVapSat(:,:,kmax    ) = xyr_QH2OVapSat(:,:,kmax-1)
    !
    xy_IndexMixLayTop     = 1
    xy_FlagMixLayTopFound = .false.
    do k = 2, kmax
      do j = 1, jmax
        do i = 0, imax-1
          if ( ( xyz_QH2OVap(i,j,1) >= xyr_QH2OVapSat(i,j,k) ) .and. ( .not. xy_FlagMixLayTopFound(i,j) )                 ) then
            xy_IndexMixLayTop    (i,j) = k - 1
            xy_FlagMixLayTopFound(i,j) = .true.
          end if
        end do
      end do
    end do
    !
    !====================================
    !
    do j = 1, jmax
      do i = 0, imax-1
        if ( .not. xy_FlagMixLayTopFound(i,j) ) then
          xy_IndexMixLayTop(i,j) = kmax - 1
        end if
      end do
    end do
    !
    !   Critical cloud work function
    !
    if ( FlagZeroCrtlCWF ) then
      xyz_CWFCrtl = 0.0_DP
    else
      do j = 1, jmax
        do i = 0, imax-1
          xy_PressCldBase(i,j) = xyr_Press(i,j,xy_IndexMixLayTop(i,j))
        end do
      end do
      call ArakawaSchubertL1982CalcCWFCrtl( xy_PressCldBase, xyz_Press, xyz_CWFCrtl )
    end if
    !
    !   Rain conversion factor
    !
    if ( RainConversionFactor < 0.0_DP ) then
      do k = 1, kmax
        do j = 1, jmax
          do i = 0, imax-1
            if ( xyz_Press(i,j,k) < 500.0d2 ) then
              xyz_RainFactor(i,j,k) = 1.0_DP
            else if ( xyz_Press(i,j,k) < 800.0d2 ) then
              xyz_RainFactor(i,j,k) = 0.8_DP + ( 800.0d2 - xyz_Press(i,j,k) ) / 1500.0d2
            else
              xyz_RainFactor(i,j,k) = 0.8_DP
            end if
          end do
        end do
      end do
    else
      xyz_RainFactor = RainConversionFactor
    end if



    xyz_RainCumulus     (:,:,1) = 0.0_DP
    xyz_EntParam        (:,:,1) = 0.0_DP
    xyz_CWF             (:,:,1) = 0.0_DP
    xyz_DCWFDtLS        (:,:,1) = 0.0_DP
    xyz_MassFluxDistFunc(:,:,1) = 0.0_DP


    if ( present( xyz_MoistConvDetTend ) ) then
      xyz_MoistConvDetTend(:,:,1) = 0.0_DP
    end if
    if ( present( xyz_MoistConvSubsidMassFlux ) ) then
      ! Subsidence mass flux between the updrafts
      ! Initialization
      !
      xyz_MoistConvSubsidMassFlux = 0.0_DP
    end if


    loop_cloud_top : do l = 2, kmax


      call RelaxedArakawaSchubertHeight( xyz_Temp, xyz_Exner, xyz_Beta, xyz_BetaCldTop, xyz_Height, xyr_Height )


      !   Potential temperature
      !
      xyz_PotTemp = xyz_Temp / xyz_Exner

      !   Saturation mixing ratio
      !
      xyz_QH2OVapSat = xyz_CalcQVapSat( xyz_Temp, xyz_Press )

      !   Calculation of dry and moist static energies
      !
      xyz_EnvDryStaticEne      = CpDry * xyz_Temp + Grav * xyz_Height
      xyz_EnvMoistStaticEne    = xyz_EnvDryStaticEne + LatentHeat * xyz_QH2OVap
      !
      k = 0
      xyr_EnvDryStaticEne  (:,:,k) = 1.0d100
      xyr_EnvMoistStaticEne(:,:,k) = 1.0d100
      do k = 1, kmax-1
        xyr_EnvDryStaticEne  (:,:,k) = ( xyz_EnvDryStaticEne  (:,:,k) + xyz_EnvDryStaticEne  (:,:,k+1) ) / 2.0_DP
        xyr_EnvMoistStaticEne(:,:,k) = ( xyz_EnvMoistStaticEne(:,:,k) + xyz_EnvMoistStaticEne(:,:,k+1) ) / 2.0_DP
      end do
      k = kmax
      xyr_EnvDryStaticEne  (:,:,k) = xyz_EnvDryStaticEne  (:,:,k)
      xyr_EnvMoistStaticEne(:,:,k) = xyz_EnvMoistStaticEne(:,:,k)

      !   Calculation of saturated moist static energy
      !
      xyz_EnvMoistStaticEneSat = xyz_EnvDryStaticEne + LatentHeat * xyz_QH2OVapSat
      !
      k = 0
      xyr_EnvMoistStaticEneSat(:,:,k) = 1.0d100
      do k = 1, kmax-1
        xyr_EnvMoistStaticEneSat(:,:,k) = ( xyz_EnvMoistStaticEneSat(:,:,k) + xyz_EnvMoistStaticEneSat(:,:,k+1) ) / 2.0_DP
      end do
      k = kmax
      xyr_EnvMoistStaticEneSat(:,:,k) = xyz_EnvMoistStaticEneSat(:,:,k)


      !   Auxiliary variables
      !
      xyz_Gamma = LatentHeat / CpDry * xyz_CalcDQVapSatDTemp( xyz_Temp, xyz_QH2OVapSat )
      !
      k = 1
      xyz_Mu (:,:,k) = 1.0d100
      xyz_Eps(:,:,k) = 1.0d100
      do k = 2, kmax
        xyz_Mu (:,:,k) = ( xyz_Exner(:,:,k  ) - xyr_Exner(:,:,k) ) / ( xyz_Exner(:,:,k) * ( 1.0_DP + xyz_Gamma(:,:,k) ) )
        xyz_Eps(:,:,k) = ( xyr_Exner(:,:,k-1) - xyz_Exner(:,:,k) ) / ( xyz_Exner(:,:,k) * ( 1.0_DP + xyz_Gamma(:,:,k) ) )
      end do


      ! Entrainment parameter
      !
      call RelaxedArakawaSchubertEntParam( l, xyz_Beta, xyz_BetaCldTop, xyz_PotTemp, xyz_EnvMoistStaticEne, xyz_EnvMoistStaticEneSat, xy_IndexMixLayTop, xy_EntParam )
      if ( l >= 3 ) then
        call RelaxedArakawaSchubertEntParam( l-1, xyz_Beta, xyz_BetaCldTop, xyz_PotTemp, xyz_EnvMoistStaticEne, xyz_EnvMoistStaticEneSat, xy_IndexMixLayTop, xy_EntParamLL )
      else
        xy_EntParamLL = 1.0d100
      end if
      if ( l <= kmax-1 ) then
        call RelaxedArakawaSchubertEntParam( l+1, xyz_Beta, xyz_BetaCldTop, xyz_PotTemp, xyz_EnvMoistStaticEne, xyz_EnvMoistStaticEneSat, xy_IndexMixLayTop, xy_EntParamUL )
      else
        xy_EntParamUL = 1.0d100
      end if
      !   for output
      xyz_EntParam(:,:,l) = xy_EntParam


      ! Difference of normalized mass flux
      !
      !   difference of normalized mass flux between layer bottom and top
      !
      xyz_DelNormMassFlux(:,:,1) = 1.0d100
      do k = 2, l-1
        xyz_DelNormMassFlux(:,:,k) = - xy_EntParam * xyz_Beta(:,:,k) * xyz_PotTemp(:,:,k)
      end do
      do k = l, kmax
        xyz_DelNormMassFlux(:,:,k) = 1.0d100
      end do
      !
      !   difference of normalized mass flux between layer bottom and mid-point
      !
      xy_DelNormMassFluxCldTop = - xy_EntParam * xyz_BetaCldTop(:,:,l) * xyz_PotTemp(:,:,l)


      ! Normalized mass flux
      !
      !   normalized mass flux at layer interface
      !
      xyr_NormMassFlux(:,:,0) = 0.0_DP
      do k = 1, l-1
        do j = 1, jmax
          do i = 0, imax-1
            if ( k < xy_IndexMixLayTop(i,j) ) then
              xyr_NormMassFlux(i,j,k) = 0.0_DP
            else if ( k == xy_IndexMixLayTop(i,j) ) then
              xyr_NormMassFlux(i,j,k) = 1.0_DP
            else
              xyr_NormMassFlux(i,j,k) = xyr_NormMassFlux(i,j,k-1) - xyz_DelNormMassFlux(i,j,k)
            end if
          end do
        end do
      end do
      do k = l, kmax
        xyr_NormMassFlux(:,:,k) = 0.0_DP
      end do
      !
      !   normalized mass flux at cloud top (at layer mid-point)
      !
      xy_NormMassFluxCldTop = xyr_NormMassFlux(:,:,l-1) - xy_DelNormMassFluxCldTop


      ! Liquid water content at top of clouds
      !   If l is less than xy_IndexMixLayTop(i,j), i.e. the cloud top is below top of 
      !   mixed layer, xy_SumTmp is zero, then, xy_CldQH2OLiqCldTop is also zero.
      !
      do j = 1, jmax
        do i = 0, imax-1

          if ( l > xy_IndexMixLayTop(i,j) ) then
            xy_SumTmp(i,j) = xyz_QH2OVap(i,j,xy_IndexMixLayTop(i,j))
            do k = xy_IndexMixLayTop(i,j)+1, l-1
              xy_SumTmp(i,j) = xy_SumTmp(i,j) - xyz_DelNormMassFlux(i,j,k) * xyz_QH2OVap(i,j,k)
            end do
            xy_SumTmp(i,j) = xy_SumTmp(i,j) - xy_DelNormMassFluxCldTop(i,j) * xyz_QH2OVap(i,j,l)
          else
            xy_SumTmp(i,j) = 0.0_DP
          end if

        end do
      end do
      xy_CldQH2OLiqCldTop = xy_SumTmp / ( xy_NormMassFluxCldTop + 1.0d-100 ) - xyz_QH2OVapSat(:,:,l)

      ! Check whether kernel is positive or negative.
      !
      do j = 1, jmax
        do i = 0, imax-1
          if ( xy_CldQH2OLiqCldTop(i,j) < 0.0_DP ) then
            xy_FlagNegH2OLiqCldTop(i,j) = .true.
          else
            xy_FlagNegH2OLiqCldTop(i,j) = .false.
          end if
        end do
      end do

      !   avoid negative value
      xy_CldQH2OLiqCldTop = max( xy_CldQH2OLiqCldTop, 0.0_DP )


      ! Moist static energy in clouds
      !
      xyr_CldMoistStaticEne(:,:,0) = 1.0d100
      do k = 1, l-1
        do j = 1, jmax
          do i = 0, imax-1
            if ( k < xy_IndexMixLayTop(i,j) ) then
              xyr_CldMoistStaticEne(i,j,k) = 1.0d100
            else if ( k == xy_IndexMixLayTop(i,j) ) then
              xyr_CldMoistStaticEne(i,j,k) = xyz_EnvMoistStaticEne(i,j,xy_IndexMixLayTop(i,j))
            else
              xyr_CldMoistStaticEne(i,j,k) = ( xyr_NormMassFlux(i,j,k-1) * xyr_CldMoistStaticEne(i,j,k-1) - xyz_DelNormMassFlux(i,j,k) * xyz_EnvMoistStaticEne(i,j,k)  ) / xyr_NormMassFlux(i,j,k)
            end if
          end do
        end do
      end do
      do k = l, kmax
        xyr_CldMoistStaticEne(:,:,k) = 1.0d100
      end do


      !###############################################
      ! Check whether a parcel in cloud has moist static energy larger than environment's
      !
!!$      xy_FlagCrossSatEquivPotTemp = .false.
!!$      do j = 1, jmax
!!$        do i = 0, imax-1
!!$          do k = xy_IndexMixLayTop(i,j), l-1
!!$            if ( xyr_EnvMoistStaticEneSat(i,j,k) < xyr_CldMoistStaticEne(i,j,k) ) then
!!$              xy_FlagCrossSatEquivPotTemp(i,j) = .true.
!!$            end if
!!$          end do
!!$        end do
!!$      end do
      !###############################################


      ! Cloud work function
      !
      xy_CWF = 0.0_DP
      do k = 2, l-1
        do j = 1, jmax
          do i = 0, imax-1
            if ( k > xy_IndexMixLayTop(i,j) ) then
              xy_CWF(i,j) = xy_CWF(i,j) + xyz_Mu (i,j,k) * xyr_NormMassFlux(i,j,k  ) * ( xyr_CldMoistStaticEne(i,j,k  ) - xyz_EnvMoistStaticEneSat(i,j,k) ) + xyz_Eps(i,j,k) * xyr_NormMassFlux(i,j,k-1) * ( xyr_CldMoistStaticEne(i,j,k-1) - xyz_EnvMoistStaticEneSat(i,j,k) )
            end if
          end do
        end do
      end do
      k = l
      do j = 1, jmax
        do i = 0, imax-1
          if ( k > xy_IndexMixLayTop(i,j) ) then
            xy_CWF(i,j) = xy_CWF(i,j) + xyz_Eps(i,j,k) * xyr_NormMassFlux(i,j,k-1) * ( xyr_CldMoistStaticEne(i,j,k-1) - xyz_EnvMoistStaticEneSat(i,j,k) )
          end if
        end do
      end do

      !   for save
      xyz_CWF(:,:,l) = xy_CWF

      ! Time derivative of cloud work function by large scale motion
      !
      xy_DCWFDtLS = ( xy_CWF - xyz_CWFCrtl(:,:,l) ) / ( 2.0_DP * DelTime )
      !   for save
      xyz_DCWFDtLS(:,:,l) = xy_DCWFDtLS

      ! Tendency of dry static energy per unit mass flux
      !
      do k = 1, l
        xyz_GammaDSE(:,:,k) = - Grav / xyz_DelPress(:,:,k) * (   xyr_NormMassFlux(:,:,k-1) * ( xyr_EnvDryStaticEne(:,:,k-1) - xyz_EnvDryStaticEne(:,:,k) ) + xyr_NormMassFlux(:,:,k  ) * ( xyz_EnvDryStaticEne(:,:,k  ) - xyr_EnvDryStaticEne(:,:,k) ) )
      end do
      k = l
      xyz_GammaDSE(:,:,k) = xyz_GammaDSE(:,:,k) - Grav / xyz_DelPress(:,:,k) * LatentHeat * xy_CldQH2OLiqCldTop * xy_NormMassFluxCldTop * ( 1.0_DP - xyz_RainFactor(:,:,k) )
      do k = l+1, kmax
        xyz_GammaDSE(:,:,k) = 0.0_DP
      end do


      ! Tendency of moist static energy per unit mass flux
      !
      do k = 1, l
        xyz_GammaMSE(:,:,k) = - Grav / xyz_DelPress(:,:,k) * (   xyr_NormMassFlux(:,:,k-1) * ( xyr_EnvMoistStaticEne(:,:,k-1) - xyz_EnvMoistStaticEne(:,:,k) ) + xyr_NormMassFlux(:,:,k  ) * ( xyz_EnvMoistStaticEne(:,:,k  ) - xyr_EnvMoistStaticEne(:,:,k) ) )
      end do
      k = l
      xyz_GammaMSE(:,:,k) = xyz_GammaMSE(:,:,k) + Grav / xyz_DelPress(:,:,k) * xy_NormMassFluxCldTop * ( xyz_EnvMoistStaticEneSat(:,:,k) - xyz_EnvMoistStaticEne(:,:,k) )
      do k = l+1, kmax
        xyz_GammaMSE(:,:,k) = 0.0_DP
      end do


      ! Kernel, time derivative of cloud work function by cumulus convection per unit 
      ! mass flux
      !
      do j = 1, jmax
        do i = 0, imax-1

          xy_Kernel(i,j) = xyz_Eps(i,j,xy_IndexMixLayTop(i,j)+1) * xyz_GammaMSE(i,j,xy_IndexMixLayTop(i,j)) - xyz_Eps(i,j,l) * xyr_NormMassFlux(i,j,l-1) * ( 1.0_DP + xyz_Gamma(i,j,l) ) * xyz_GammaDSE(i,j,l)
          do n = xy_IndexMixLayTop(i,j)+1, l-1
            xy_SumTmp(i,j) = 0.0_DP
            do m = xy_IndexMixLayTop(i,j)+1, n
              xy_SumTmp(i,j) = xy_SumTmp(i,j) + xyz_DelNormMassFlux(i,j,m) * xyz_GammaMSE(i,j,m)
            end do
            xy_Kernel(i,j) = xy_Kernel(i,j) + ( xyz_Eps(i,j,n+1) + xyz_Mu(i,j,n) ) * ( xyz_GammaMSE(i,j,xy_IndexMixLayTop(i,j)) - xy_SumTmp(i,j) ) - (   xyz_Eps(i,j,n) * xyr_NormMassFlux(i,j,n-1) + xyz_Mu (i,j,n) * xyr_NormMassFlux(i,j,n  ) ) * ( 1.0_DP + xyz_Gamma(i,j,n) ) * xyz_GammaDSE(i,j,n)
          end do

        end do
      end do

      ! Check whether kernel is positive or negative.
      !
      do j = 1, jmax
        do i = 0, imax-1
          if ( xy_Kernel(i,j) < 0.0_DP ) then
            xy_FlagKernelNegative(i,j) = .true.
          else
            xy_FlagKernelNegative(i,j) = .false.
          end if
        end do
      end do


      ! Load et al. (1982), p.108
      xy_Kernel = min( xy_Kernel, -5.0d-3 )


      ! Cloud mass flux at cloud bottom
      !
      xy_CldMassFluxBottom = - xy_DCWFDtLS / xy_Kernel
      !
      !   mass flux has to be zero or positive
      xy_CldMassFluxBottom = max( xy_CldMassFluxBottom, 0.0_DP )
      !   mass flux is zero if entrainment parameter is zero or negative
      do j = 1, jmax
        do i = 0, imax-1
          if ( xy_EntParam(i,j) <= 0.0_DP ) then
            xy_CldMassFluxBottom(i,j) = 0.0_DP
          end if
        end do
      end do
!!$      !   mass flux is zero if it is below lifting condensation level
!!$      do j = 1, jmax
!!$        do i = 0, imax-1
!!$          if ( .not. xy_FlagCrossSatEquivPotTemp(i,j) ) then
!!$            xy_CloudMassFluxBottom(i,j) = 0.0_DP
!!$          end if
!!$        end do
!!$      end do
      !   mass flux is zero in an 'inversion layer'
      !    Lord et al. (1982), p.112,
      !    Strictly speaking, the process below is different from that 
      !    proposed by Lord et al. (1982). Lord et al. (1982) compare 
      !    entrainment parameters at 3 levels. But, entrainment 
      !    parameters at 2 levels are compared below, because comparison 
      !    of values between 2 levels seems to be sufficient.
!!$      if ( ( 3 <= l ) .and. ( l <= kmax-1 ) ) then
!!$        do j = 1, jmax
!!$          do i = 0, imax-1
!!$            if ( ( xy_EntParamLL(i,j) < xy_EntParam  (i,j) ) .and. &
!!$              &  ( xy_EntParam  (i,j) < xy_EntParamUL(i,j) ) ) then
!!$              if ( ( xy_EntParamLL(i,j) > 0.0_DP ) .and. &
!!$                &  ( xy_EntParam  (i,j) > 0.0_DP ) .and. &
!!$                &  ( xy_EntParamUL(i,j) > 0.0_DP ) ) then
!!$                xy_CloudMassFluxBottom(i,j) = 0.0_DP
!!$              end if
!!$            end if
!!$          end do
!!$        end do
!!$      end if
      do j = 1, jmax
        do i = 0, imax-1
!!$          if ( xy_IndexMixLayTop(i,j) == l ) then
!!$            if ( ( xy_EntParam  (i,j) > 0.0_DP ) .and. &
!!$              &  ( xy_EntParamUL(i,j) > 0.0_DP ) ) then
!!$              if ( xy_EntParam  (i,j) < xy_EntParamUL(i,j) ) then
!!$                xy_CloudMassFluxBottom(i,j) = 0.0_DP
!!$              end if
!!$            end if
!!$          else if ( ( xy_IndexMixLayTop(i,j) < l ) .and. ( l <= kmax-1 ) ) then
!!$            if ( ( xy_EntParamLL(i,j) > 0.0_DP ) .and. &
!!$              &  ( xy_EntParam  (i,j) > 0.0_DP ) .and. &
!!$              &  ( xy_EntParamUL(i,j) > 0.0_DP ) ) then
!!$              if ( ( xy_EntParamLL(i,j) < xy_EntParam  (i,j) ) .and. &
!!$                &  ( xy_EntParam  (i,j) < xy_EntParamUL(i,j) ) ) then
          if ( ( xy_IndexMixLayTop(i,j) <= l ) .and. ( l <= kmax-1 ) ) then
            if ( ( xy_EntParam  (i,j) > 0.0_DP ) .and. ( xy_EntParamUL(i,j) > 0.0_DP ) ) then
              if ( xy_EntParam  (i,j) < xy_EntParamUL(i,j) ) then
                xy_CldMassFluxBottom(i,j) = 0.0_DP
              end if
            end if
          end if
        end do
      end do
      !
      !   mass flux is zero unless kernel is negative
      !
      do j = 1, jmax
        do i = 0, imax-1
          if ( .not. xy_FlagKernelNegative(i,j) ) then
            xy_CldMassFluxBottom(i,j) = 0.0_DP
          end if
        end do
      end do
      !
      !   mass flux is zero if liquid water at a cloud top is negative
      !
      do j = 1, jmax
        do i = 0, imax-1
          if ( xy_FlagNegH2OLiqCldTop(i,j) ) then
            xy_CldMassFluxBottom(i,j) = 0.0_DP
          end if
        end do
      end do
      !
      !   multiply factor
      !
      xy_CldMassFluxBottom = xy_CldMassFluxBottom * min( 2.0_DP * DelTime / AdjTimeConst, 1.0_DP )
      !
      !   for output
      xyz_MassFluxDistFunc(:,:,l) = xy_CldMassFluxBottom



      ! Check values of cloud mass flux
      !   If water vapor amount transported by convection is larger than that in a 
      !   column, cloud mass flux is reduced.
      !
      !   tendency of specific humidity is calculated tentatively
      do k = 1, kmax
        xyz_DQVapDtCumulus(:,:,k) = + xy_CldMassFluxBottom * ( xyz_GammaMSE(:,:,k) - xyz_GammaDSE(:,:,k) ) / LatentHeat
      end do
      !   total H2O mass in a vertical column after RAS
      xyz_QH2OVapTentative = xyz_QH2OVap + xyz_DQVapDtCumulus * 2.0_DP * DelTime
      xy_CldMassFluxCorFactor = 1.0_DP
      do k = 1, kmax
        do j = 1, jmax
          do i = 0, imax-1
            if ( xyz_QH2OVapTentative(i,j,k) < 0.0_DP ) then
              CldMassFluxCorFactor = xyz_QH2OVap(i,j,k) / ( xyz_QH2OVap(i,j,k) - xyz_QH2OVapTentative(i,j,k) )
            else
              CldMassFluxCorFactor = 1.0_DP
            end if
            if ( CldMassFluxCorFactor < xy_CldMassFluxCorFactor(i,j) ) then
              xy_CldMassFluxCorFactor(i,j) = CldMassFluxCorFactor
            end if
          end do
        end do
      end do
      !   modify cloud mass flux
      xy_CldMassFluxBottom = xy_CldMassFluxCorFactor * xy_CldMassFluxBottom



!!$      do k = 1, kmax
!!$        xyz_DQVapDtCumulus(:,:,k) =                                                  &
!!$          & + xy_CloudMassFluxBottom * ( xyz_GammaMSE(:,:,k) - xyz_GammaDSE(:,:,k) ) &
!!$          &     / LatentHeat
!!$      end do
!!$      !   total H2O mass in a vertical column before RAS
!!$      xyz_DelH2OMass = xyz_QH2OVap * xyz_DelPress / Grav
!!$      xy_H2OMassB = 0.0_DP
!!$      do k = kmax, 1, -1
!!$        xy_H2OMassB = xy_H2OMassB + xyz_DelH2OMass(:,:,k)
!!$      end do
!!$      !   total H2O mass in a vertical column after RAS
!!$      xyz_QH2OVapTentative = xyz_QH2OVap + xyz_DQVapDtCumulus * 2.0_DP * DelTime
!!$      xyz_DelH2OMass = xyz_QH2OVapTentative * xyz_DelPress / Grav
!!$      xy_H2OMassA = 0.0_DP
!!$      do k = kmax, 1, -1
!!$        xy_H2OMassA = xy_H2OMassA + xyz_DelH2OMass(:,:,k)
!!$      end do
!!$      !   modify cloud mass flux
!!$      do j = 1, jmax
!!$        do i = 0, imax-1
!!$          if ( xy_H2OMassA(i,j) < 0.0_DP ) then
!!$            ! A safety factor ( 1.0_DP + 1.0d-5 ) is arbitrary.
!!$            xy_CloudMassFluxBottom(i,j) = xy_CloudMassFluxBottom(i,j)                 &
!!$              & * xy_H2OMassB(i,j)                                                    &
!!$              &     / ( ( xy_H2OMassB(i,j) - xy_H2OMassA(i,j) ) * ( 1.0_DP + 1.0d-5 ) )
!!$          end if
!!$        end do
!!$      end do



      ! Tendencies of specific temperature and humidity
      !
      do k = 1, kmax
        xyz_DTempDtCumulus(:,:,k) = + xy_CldMassFluxBottom * xyz_GammaDSE(:,:,k) / CpDry
        xyz_DQVapDtCumulus(:,:,k) = + xy_CldMassFluxBottom * ( xyz_GammaMSE(:,:,k) - xyz_GammaDSE(:,:,k) ) / LatentHeat
      end do
!!$      !
!!$      !   modification of tendency of temperature and water vapor in the mixed layer
!!$      !
!!$      if ( FlagUniformMixedLayer ) then
!!$        xy_SumTmp = 0.0_DP
!!$        do k = 1, kmax
!!$          do j = 1, jmax
!!$            do i = 0, imax-1
!!$              if ( k <= xy_IndexMixLayTop(i,j) ) then
!!$                xy_SumTmp(i,j) = xy_SumTmp(i,j) &
!!$                  & + xyz_DTempDtCumulus(i,j,k) &
!!$                  &     * ( xyr_Press(i,j,k-1) - xyr_Press(i,j,k) )
!!$              end if
!!$            end do
!!$          end do
!!$        end do
!!$        do j = 1, jmax
!!$          do i = 0, imax-1
!!$            xy_SumTmp(i,j) = xy_SumTmp(i,j) &
!!$              & / ( xyr_Press(i,j,0) - xyr_Press(i,j,xy_IndexMixLayTop(i,j)) )
!!$          end do
!!$        end do
!!$        do k = 1, kmax
!!$          do j = 1, jmax
!!$            do i = 0, imax-1
!!$              if ( k <= xy_IndexMixLayTop(i,j) ) then
!!$                xyz_DTempDtCumulus(i,j,k) = xy_SumTmp(i,j)
!!$              end if
!!$            end do
!!$          end do
!!$        end do
!!$        !
!!$        xy_SumTmp = 0.0_DP
!!$        do k = 1, kmax
!!$          do j = 1, jmax
!!$            do i = 0, imax-1
!!$              if ( k <= xy_IndexMixLayTop(i,j) ) then
!!$                xy_SumTmp(i,j) = xy_SumTmp(i,j) &
!!$                  & + xyz_DQVapDtCumulus(i,j,k) &
!!$                  &     * ( xyr_Press(i,j,k-1) - xyr_Press(i,j,k) )
!!$              end if
!!$            end do
!!$          end do
!!$        end do
!!$        do j = 1, jmax
!!$          do i = 0, imax-1
!!$            xy_SumTmp(i,j) = xy_SumTmp(i,j) &
!!$              & / ( xyr_Press(i,j,0) - xyr_Press(i,j,xy_IndexMixLayTop(i,j)) )
!!$          end do
!!$        end do
!!$        do k = 1, kmax
!!$          do j = 1, jmax
!!$            do i = 0, imax-1
!!$              if ( k <= xy_IndexMixLayTop(i,j) ) then
!!$                xyz_DQVapDtCumulus(i,j,k) = xy_SumTmp(i,j)
!!$              end if
!!$            end do
!!$          end do
!!$        end do
!!$      end if



      ! add tendencies to temperature and specific humidity
      !
      xyz_Temp    = xyz_Temp    + xyz_DTempDtCumulus * 2.0_DP * DelTime
      xyz_QH2OVap = xyz_QH2OVap + xyz_DQVapDtCumulus * 2.0_DP * DelTime


      ! Precipitation rate at cloud top level
      !   unit is kg m-2 s-1
      !
      xyz_RainCumulus(:,:,l) = xy_CldMassFluxBottom * xyz_RainFactor(:,:,l) * xy_NormMassFluxCldTop * xy_CldQH2OLiqCldTop



      ! mass fix
      !
      xyz_DelH2OMass = xyz_QH2OVap * xyz_DelPress / Grav
      !   total H2O mass in a vertical column
      xy_H2OMassB = 0.0_DP
      do k = kmax, 1, -1
        xy_H2OMassB = xy_H2OMassB + xyz_DelH2OMass(:,:,k)
      end do
      do j = 1, jmax
        do i = 0, imax-1
          if ( xy_H2OMassB(i,j) < 0.0_DP ) then
            call MessageNotify( 'E', module_name, 'Mass of water vapor in a column is negative (%d,%d), %f.', i = (/i,j/), d = (/xy_H2OMassB(i,j)/) )
          end if
        end do
      end do
      !   negative mass is borrowed from above
      do k = 1, kmax-1
        do j = 1, jmax
          do i = 0, imax-1
            if ( xyz_DelH2OMass(i,j,k) < 0.0_DP ) then
              xyz_DelH2OMass(i,j,k+1) = xyz_DelH2OMass(i,j,k+1) + xyz_DelH2OMass(i,j,k)
              xyz_DelH2OMass(i,j,k  ) = 0.0_DP
            end if
          end do
        end do
      end do
      k = kmax
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_DelH2OMass(i,j,k) < 0.0_DP ) then

!!$            call MessageNotify( 'E', module_name,                                   &
!!$              & 'Mass of water vapor in the top layer is negative (%d,%d,%d), %f.', &
!!$              & i = (/i,j,k/), d = (/xyz_DelH2OMass(i,j,k)/) )
!!$
!!$            xyz_RainCumulus(i,j,l) = xyz_RainCumulus(i,j,l) &
!!$              & - xyz_DelH2OMass(i,j,k) / ( 2.0_DP * DelTime )
!!$            xyz_Temp       (i,j,k) = xyz_Temp(i,j,k)                                 &
!!$              & - LatentHeat * xyz_DelH2OMass(i,j,k) / ( xyz_DelPress(i,j,k) / Grav )&
!!$              &    / CpDry

            xyz_DelH2OMass (i,j,k) = 0.0_DP
          end if
        end do
      end do
!!$      do j = 1, jmax
!!$        do i = 0, imax-1
!!$          if ( xyz_RainCumulus(i,j,l) < 0.0_DP ) then
!!$            call MessageNotify( 'E', module_name, &
!!$              & 'Mass of water vapor is insufficient at (%d,%d,%d), %f.', &
!!$              & i = (/i,j,k/), d = (/xyz_RainCumulus(i,j,l)/) )
!!$          end if
!!$        end do
!!$      end do


      !   total H2O mass in a vertical column, again
      xy_H2OMassA = 0.0_DP
      do k = kmax, 1, -1
        xy_H2OMassA = xy_H2OMassA + xyz_DelH2OMass(:,:,k)
      end do
      !   total mass in a vertical column is adjusted
      do j = 1, jmax
        do i = 0, imax-1
          if ( xy_H2OMassA(i,j) > 0.0_DP ) then
!!$            write( 6, * ) i, j, xy_H2OMassB(i,j), xy_H2OMassB(i,j) / xy_H2OMassA(i,j)
            do k = 1, kmax
              xyz_DelH2OMass(i,j,k) = xyz_DelH2OMass(i,j,k) * xy_H2OMassB(i,j) / xy_H2OMassA(i,j)
            end do
          else
            do k = 1, kmax
              xyz_DelH2OMass(i,j,k) = 0.0_DP
            end do
          end if
        end do
      end do
      xyz_QH2OVap = xyz_DelH2OMass / ( xyz_DelPress / Grav )


      ! Detrainment mass tendency per unit mass (kg m-3 s-1 / ( kg m-3 ) = s-1).
      ! This corresponds to condensation rate (kg m-2 s-1) divided by layer thickness (m)
      ! and density (kg m-3), in other words.

      ! kg m-2 s-1 / ( Pa / ( m s-2 ) )
      ! = kg m-2 s-1 Pa-1 m s-1 = kg m-2 (kg m s-2 m-2)-1 m s-2
      ! = kg m-2 s-1 kg-1 m-1 s2 m2 m s-2 = s-1

      if ( present( xyz_MoistConvDetTend ) ) then
        xyz_MoistConvDetTend(:,:,l) = xy_CldMassFluxBottom * xy_NormMassFluxCldTop / ( xyz_DelPress(:,:,l) / Grav )
      end if

      if ( present( xyz_MoistConvSubsidMassFlux ) ) then
        ! Subsidence mass flux between the updrafts
        do k = 1, l-1
          do j = 1, jmax
            do i = 0, imax-1
              if ( k > xy_IndexMixLayTop(i,j) ) then
                DelNormMassFluxHalfLayer = - xy_EntParam(i,j) * xyz_BetaCldTop(i,j,k) * xyz_PotTemp(i,j,k)
                NormMassFlux = xyr_NormMassFlux(i,j,k-1) - DelNormMassFluxHalfLayer
                xyz_MoistConvSubsidMassFlux(i,j,k) = xyz_MoistConvSubsidMassFlux(i,j,k) + xy_CldMassFluxBottom(i,j) * NormMassFlux
              end if
            end do
          end do
        end do
      end if


    end do loop_cloud_top




    ! 温度変化率, 比湿変化率
    ! Calculate specific humidity tendency and temperature tendency
    !   (In fact, temperature tendency does not need to calculate, here.)
    !
    xyz_DTempDtCumulus = ( xyz_Temp    - xyz_TempB    ) / ( 2.0_DP * DelTime )
    xyz_DQVapDtCumulus = ( xyz_QH2OVap - xyz_QH2OVapB ) / ( 2.0_DP * DelTime )


    xyz_DTempDt = xyz_DTempDt + xyz_DTempDtCumulus
    xyz_DQVapDt = xyz_DQVapDt + xyz_DQVapDtCumulus



    ! Precipitation rate at the surface
    !   unit is kg m-2 s-1
    !
!!$    xy_RainCumulus = 0.0d0
!!$    do k = kmax, 1, -1
!!$      xy_RainCumulus = xy_RainCumulus + xyz_RainCumulus(:,:,k)
!!$    end do


    xyz_DQH2OLiqDt = xyz_RainCumulus / ( xyz_DelPress / Grav )

!!$    xyz_RainCumulus = xyz_DQH2OLiqDt * ( xyz_DelPress / Grav )
!!$    xy_RainCumulus = 0.0d0
!!$    do k = kmax, 1, -1
!!$      xy_RainCumulus = xy_RainCumulus + xyz_RainCumulus(:,:,k)
!!$    end do
!!$
!!$    xy_Rain     = xy_Rain     + xy_RainCumulus



    ! Calculation for debug
    !   check of conservation of water amount and internal energy
    !
!!$    xyz_DelVal = xyz_QH2OVapB * xyz_DelPress / Grav
!!$    xy_SumValB = 0.0_DP
!!$    do k = kmax, 1, -1
!!$      xy_SumValB = xy_SumValB + xyz_DelVal(:,:,k)
!!$    end do
!!$    !
!!$    xyz_DelVal = xyz_QH2OVap * xyz_DelPress / Grav
!!$    xy_SumValA = 0.0_DP
!!$    do k = kmax, 1, -1
!!$      xy_SumValA = xy_SumValA + xyz_DelVal(:,:,k)
!!$    end do
!!$    !
!!$    xy_SumValA = xy_SumValA + xy_RainCumulus * 2.0_DP * DelTime
!!$    !
!!$    do j = 1, jmax
!!$      do i = 0, imax-1
!!$        Ratio = ( xy_SumValA(i,j) - xy_SumValB(i,j) ) &
!!$          & / max( xy_SumValA(i,j), 1.0d-100 )
!!$        if ( abs( Ratio ) > 1.0d-14 ) then
!!$          write( 6, * ) 'H2O: ', i, j, &
!!$            & xy_SumValB(i,j), xy_SumValA(i,j), &
!!$            & xy_RainCumulus(i,j) * 2.0_DP * DelTime, &
!!$            & Ratio
!!$        end if
!!$      end do
!!$    end do
!!$    !
!!$    !
!!$    xyz_DelVal = CpDry * xyz_TempB * xyz_DelPress / Grav
!!$    xy_SumValB = 0.0_DP
!!$    do k = kmax, 1, -1
!!$      xy_SumValB = xy_SumValB + xyz_DelVal(:,:,k)
!!$    end do
!!$    !
!!$    xyz_DelVal = CpDry * xyz_Temp * xyz_DelPress / Grav
!!$    xy_SumValA = 0.0_DP
!!$    do k = kmax, 1, -1
!!$      xy_SumValA = xy_SumValA + xyz_DelVal(:,:,k)
!!$    end do
!!$    !
!!$    xy_SumValA = xy_SumValA - LatentHeat * xy_RainCumulus * 2.0_DP * DelTime
!!$    !
!!$    do j = 1, jmax
!!$      do i = 0, imax-1
!!$        Ratio = ( xy_SumValA(i,j) - xy_SumValB(i,j) ) &
!!$          & / max( xy_SumValA(i,j), 1.0d-100 )
!!$        if ( abs( Ratio ) > 1.0d-14 ) then
!!$          write( 6, * ) 'CpT: ', i, j, &
!!$            & xy_SumValB(i,j), xy_SumValA(i,j), &
!!$            & - LatentHeat * xy_RainCumulus(i,j) * 2.0_DP * DelTime, &
!!$            & Ratio
!!$        end if
!!$      end do
!!$    end do


    ! calculation for output
    xyz_RainCumulus = xyz_DQH2OLiqDt * ( xyz_DelPress / Grav )
    xy_RainCumulus = 0.0d0
    do k = kmax, 1, -1
      xy_RainCumulus = xy_RainCumulus + xyz_RainCumulus(:,:,k)
    end do


    ! ヒストリデータ出力
    ! History data output
    !
    call HistoryAutoPut( TimeN, 'RainCumulus'        , xy_RainCumulus * LatentHeat )
    call HistoryAutoPut( TimeN, 'DTempDtCumulus'     , xyz_DTempDtCumulus          )
    call HistoryAutoPut( TimeN, 'DQVapDtCumulus'     , xyz_DQVapDtCumulus          )
    call HistoryAutoPut( TimeN, 'RASMassFluxDistFunc', xyz_MassFluxDistFunc        )
    call HistoryAutoPut( TimeN, 'RASEntParam'        , xyz_EntParam                )
    call HistoryAutoPut( TimeN, 'RASCWF'             , xyz_CWF                     )
    call HistoryAutoPut( TimeN, 'RASCWFCrtl'         , xyz_CWFCrtl                 )
    call HistoryAutoPut( TimeN, 'RASDCWFDtLS'        , xyz_DCWFDtLS                )
    call HistoryAutoPut( TimeN, 'RASMixLayTopIndex'  , real( xy_IndexMixLayTop )   )



!!$    if ( present( xyz_DQH2OLiqDt ) ) then
!!$
!!$      !   unit is kg m-2 s-1
!!$      xyz_DDelLWDtCCPLV = xyz_RainCumulus
!!$
!!$      ! Negative cloud production rate is filled with values in lower layers.
!!$      !
!!$      xy_NegDDelLWDt = 0.0d0
!!$      do k = kmax, 1, -1
!!$        do j = 1, jmax
!!$          do i = 0, imax-1
!!$            xyz_DDelLWDtCCPLV(i,j,k) = xyz_DDelLWDtCCPLV(i,j,k) + xy_NegDDelLWDt(i,j)
!!$            if ( xyz_DDelLWDtCCPLV(i,j,k) < 0.0d0 ) then
!!$              xy_NegDDelLWDt(i,j) = xyz_DDelLWDtCCPLV(i,j,k)
!!$              xyz_DDelLWDtCCPLV(i,j,k) = 0.0d0
!!$            end if
!!$          end do
!!$        end do
!!$      end do
!!$
!!$      !   unit is s-1
!!$      xyz_DQH2OLiqDt = xyz_DDelLWDtCCPLV / ( xyz_DelPress / Grav )
!!$
!!$    end if


    ! 計算時間計測一時停止
    ! Pause measurement of computation time
    !
    call TimesetClockStop( module_name )

  end subroutine RelaxedArakawaSchubert

  subroutine RelaxedArakawaSchubertInit
    !
    ! moist_conv_adjust モジュールの初期化を行います. 
    ! NAMELIST#moist_conv_adjust_nml の読み込みはこの手続きで行われます. 
    !
    ! "moist_conv_adjust" module is initialized. 
    ! "NAMELIST#moist_conv_adjust_nml" is loaded in this procedure. 
    !

    ! モジュール引用 ; USE statements
    !

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

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

    ! 文字列操作
    ! Character handling
    !
    use dc_string, only: StoA

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

    ! Arakawa-Schubert scheme by Lord et al. (1982)
    ! Arakawa-Schubert scheme by Lord et al. (1982)
    !
    use arakawa_schubert_L1982, only : ArakawaSchubertL1982Init

    ! 宣言文 ; Declaration statements
    !
    implicit none

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

    ! NAMELIST 変数群
    ! NAMELIST group name
    !
    namelist /relaxed_arakawa_schubert/ AdjTimeConst, RainConversionFactor, FlagZeroCrtlCWF
          ! デフォルト値については初期化手続 "moist_conv_adjust#CumAdjInit" 
          ! のソースコードを参照のこと. 
          !
          ! Refer to source codes in the initialization procedure
          ! "moist_conv_adjust#MoistConvAdjustInit" for the default values. 
          !

    ! 実行文 ; Executable statement
    !

    if ( relaxed_arakawa_schubert_inited ) return


    ! デフォルト値の設定
    ! Default values settings
    !
!!$    FlagUse               = .true.
!!$    FlagUniformMixedLayer = .false.
    AdjTimeConst          = 7200.0_DP
    RainConversionFactor  = -1.0_DP
    FlagZeroCrtlCWF       = .false.

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

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


    ! Check values
    !
    if ( RainConversionFactor > 1.0_DP ) then
      call MessageNotify( 'E', module_name, 'RainConversionFactor is %f, but it must be less than or equal to 1', d = (/ RainConversionFactor /) )
    end if


    ! ヒストリデータ出力のためのへの変数登録
    ! Register of variables for history data output
    !
    call HistoryAutoAddVariable( 'RainCumulus', (/ 'lon ', 'lat ', 'time' /), 'precipitation by cumulus scheme, RAS', 'W m-2' )
    call HistoryAutoAddVariable( 'DTempDtCumulus', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'cumulus condensation heating, RAS', 'K s-1' )
    call HistoryAutoAddVariable( 'DQVapDtCumulus', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'cumulus condensation moistening, RAS', 'kg kg-1 s-1' )
    call HistoryAutoAddVariable( 'RASMassFluxDistFunc', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'mass flux distribution function, RAS', 'kg m-2 s-1' )
    call HistoryAutoAddVariable( 'RASEntParam', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'entrainment parameter', 'm-1' )
    call HistoryAutoAddVariable( 'RASCWF', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'cloud work function', 'J kg-1' )
    call HistoryAutoAddVariable( 'RASCWFCrtl', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'critical cloud work function', 'J kg-1' )
    call HistoryAutoAddVariable( 'RASDCWFDtLS', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'time derivative of cloud work function by large scale', 'J kg-1 s-1' )
    call HistoryAutoAddVariable( 'RASMixLayTopIndex', (/ 'lon ', 'lat ', 'time' /), 'index of top of mixed layer', '1' )


    ! Initialization of modules used in this module
    !

    ! Arakawa-Schubert scheme by Lord et al. (1982)
    ! Arakawa-Schubert scheme by Lord et al. (1982)
    !
    call ArakawaSchubertL1982Init


    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
!!$    call MessageNotify( 'M', module_name, '  FlagUse               = %b', l = (/ FlagUse /) )
!!$    call MessageNotify( 'M', module_name, '  FlagUniformMixedLayer = %b', l = (/ FlagUniformMixedLayer /) )
    call MessageNotify( 'M', module_name, '  AdjTimeConst          = %f', d = (/ AdjTimeConst /) )
    call MessageNotify( 'M', module_name, '  RainConversionFactor  = %f', d = (/ RainConversionFactor /) )
    call MessageNotify( 'M', module_name, '  FlagZeroCrtlCWF       = %b', l = (/ FlagZeroCrtlCWF /) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )


    relaxed_arakawa_schubert_inited = .true.

  end subroutine RelaxedArakawaSchubertInit

  subroutine RelaxedArakawaSchubertEntParam( l, xyz_Beta, xyz_BetaCldTop, xyz_PotTemp, xyz_EnvMoistStaticEne, xyz_EnvMoistStaticEneSat, xy_IndexMixLayTop, xy_EntParam )
    !
    ! エントレインメントパラメータの計算
    !
    ! Calculation of entrainment parameter
    !

    ! モジュール引用 ; USE statements
    !

    ! 宣言文 ; Declaration statements
    !

    integer , intent(in ) :: l
    real(DP), intent(in ) :: xyz_Beta                (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ) :: xyz_BetaCldTop          (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ) :: xyz_PotTemp             (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ) :: xyz_EnvMoistStaticEne   (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ) :: xyz_EnvMoistStaticEneSat(0:imax-1, 1:jmax, 1:kmax)
    integer , intent(in ) :: xy_IndexMixLayTop       (0:imax-1, 1:jmax)
    real(DP), intent(out) :: xy_EntParam             (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
    integer  :: k               ! 鉛直方向に回る DO ループ用作業変数
                                ! Work variables for DO loop in vertical direction


    ! 実行文 ; Executable statement
    !


    ! Entrainment parameter
    !
    xy_EntParam = 0.0_DP
    do k = 2, l-1

      do j = 1, jmax
        do i = 0, imax-1

          if ( k > xy_IndexMixLayTop(i,j) ) then
            xy_EntParam(i,j) = xy_EntParam(i,j) + xyz_Beta(i,j,k) * xyz_PotTemp(i,j,k) * ( xyz_EnvMoistStaticEneSat(i,j,l) - xyz_EnvMoistStaticEne(i,j,k) )
          end if

        end do
      end do

    end do

    do j = 1, jmax
      do i = 0, imax-1

        if ( l > xy_IndexMixLayTop(i,j) ) then
          xy_EntParam(i,j) = xy_EntParam(i,j) + xyz_BetaCldTop(i,j,l) * xyz_PotTemp(i,j,l) * ( xyz_EnvMoistStaticEneSat(i,j,l) - xyz_EnvMoistStaticEne(i,j,l) )
          xy_EntParam(i,j) = ( xyz_EnvMoistStaticEne(i,j,xy_IndexMixLayTop(i,j)) - xyz_EnvMoistStaticEneSat(i,j,l) ) / ( xy_EntParam(i,j) + 1.0d-100 )
        end if

      end do
    end do


  end subroutine RelaxedArakawaSchubertEntParam

  subroutine RelaxedArakawaSchubertHeight( xyz_Temp, xyz_Exner, xyz_Beta, xyz_BetaCldTop, xyz_Height, xyr_Height )
    !
    ! 高度の計算
    !
    ! Calculation of height
    !

    ! モジュール引用 ; USE statements
    !

    ! 宣言文 ; Declaration statements
    !

    real(DP), intent(in ) :: xyz_Temp      (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ) :: xyz_Exner     (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ) :: xyz_Beta      (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ) :: xyz_BetaCldTop(0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(out) :: xyz_Height    (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(out) :: xyr_Height    (0:imax-1, 1:jmax, 0:kmax)

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

    integer  :: k               ! 鉛直方向に回る DO ループ用作業変数
                               ! Work variables for DO loop in vertical direction


    ! 実行文 ; Executable statement
    !
    xyz_PotTemp = xyz_Temp / xyz_Exner

    xyr_Height(:,:,0) = 0.0_DP
    do k = 1, kmax
      xyz_Height(:,:,k) = xyr_Height(:,:,k-1) + xyz_BetaCldTop(:,:,k) * xyz_PotTemp(:,:,k)
      xyr_Height(:,:,k) = xyr_Height(:,:,k-1) + xyz_Beta      (:,:,k) * xyz_PotTemp(:,:,k)
    end do


  end subroutine RelaxedArakawaSchubertHeight